A Program in MATLAB to solve Air Standard Cycle

A program in Matlab to solve Air Standard Cycle

AIM:- To write a program in MATLAB to solve Air Standard Cycle (Otto Cycle).

OBJECTIVE:-

• To plot P-V diagram of the air standard cycle.
• To calculate the thermal efficiency of the engine (Gasoline engine).

INTRODUCTION:-

• Internal Combustion Engine: An internal combustion engine (ICE) is a heat engine in which the combustion of a fuel occurs with an oxidizer (usually air) in a combustion chamber that is an integral part of the working fluid flow circuit. In an internal combustion engine, the expansion of the high-temperature and high-pressure gases produced by combustion applies direct force to some component of the engine. The force is applied typically to pistons, turbine blades, rotor or a nozzle. This force moves the component over a distance, transforming chemical energy into useful work.

• Air Standard Cycle: Air standard cycles are reference cycles which give an approximation to the performance of internal combustion engines.

• Otto cycle: The Otto cycle is an air-standard cycle which approximates the processes in petrol or diesel engines. It is based on constant volume heat addition (combustion) and heat rejection processes, and isentropic compression and expansion.

fig.1: Otto Cycle P-V diagram

BODY OF CONTENT:-

• MATLAB programing language is used.

clear all
close all
clc

% Inputs
gamma = 1.4
t3 = 2300

% State variables
p1 = 101325
t1 = 500

%Engine geoetric parameters
bore = 0.1;
stroke = 0.1;
con_rod = 0.15;
cr = 12;

% Calculating the swept volume and the clearance volume
v_swept = (pi/4)*bore^2*stroke;
v_clearance = v_swept/(cr-1);

v1 = v_swept + v_clearance
v2 = v_clearance

% State variables at state point 2
% p2v2^gamma = p1v1^gamma
% p2 = p1*(v1/v2)^gama = p1*cr^gamma
p2 = p1*cr^gamma;

%p1v1/t1 = p2v2/t2 | t2 = p2*v2*t1/(p1*v1)
t2 = p2*v2*t1/(p1*v1);

% State variables at state point 3
v3 = v2

% p3v3/t3 = p2v2/t2 | p3 = p2*t3/t2
p3 = p2*t3/t2

% State variables at state point 4
v4 = v1

% p3v3^gamma = p4v4^gamma | p4 = p3(v3/v4)^gamma
p4 = p3*(v3/v4)^gamma

figure(1)
hold on

plot(v1,p1,'*','color','r')
plot(v2,p2,'*','color','r')
plot(v3,p3,'*','color','r')
plot(v4,p4,'*','color','r')

xlabel('volume,V')
ylabel('pressure,P')

From this code we get the below P-V plot :

fig.2: Point P-V diagram

But in this P-V plot we get only location of points.

Now, we have to join these points in a proper manner by following the below code applying further in the above written code for finding the points location-

figure(2)
plot([v1 v2 v3 v4 v1], [p1 p2 p3 p4 p1], 'color', 'b')

xlabel('volume,V')
ylabel('pressure,P')

After by this code, we get the P-V plot as below-

fig.3: Straight line P-V diagram

Here, we get straight lines in the P-V plot diagram. Which is right for process 2-3 & 4-1 only because of constant volume, but wrong for the processes 1-2 & 4-1 because volume is variable.

Now we are going to use Piston Kinematics for correction of our plot.

Piston Kinematics equation is,

• Piston Kinematics formula determines how the volume inside combustion chamber is going to change.

Now next step is to code a program to calculate volume information during compression and expansion stroke.

So we have to use “function” to call the input data for calculation.

For that, the program will be

function [v] = engine_kinematics(bore, stroke, con_rod, cr, start_crank, end_crank)


a = stroke/2
R = con_rod/a

v_s = pi/4*bore^2*stroke;
v_c = v_s/(cr-1);

theta = linspace(start_crank, end_crank, 100);

term1 = 0.5*(cr-1);
term2 = R + 1 - cosd(theta);
term3 = (R^2 - sind(theta).^2).^0.5;

v = (1 + term1*(term2 - term3)).*v_c;

end

Final code-

clear all
close all
clc

% Inputs
gamma = 1.4
t3 = 2300

% State variables
p1 = 101325
t1 = 500

%Engine geoetric parameters
bore = 0.1;
stroke = 0.1;
con_rod = 0.15;
cr = 12;

% Calculating the swept volume and the clearance volume
v_swept = (pi/4)*bore^2*stroke;
v_clearance = v_swept/(cr-1);

v1 = v_swept + v_clearance
v2 = v_clearance

% State variables at state point 2
% p2v2^gamma = p1v1^gamma
% p2 = p1*(v1/v2)^gama = p1*cr^gamma
p2 = p1*cr^gamma;

%p1v1/t1 = p2v2/t2 | t2 = p2*v2*t1/(p1*v1)
t2 = p2*v2*t1/(p1*v1);

constant_c = p1*v1^gamma;
v_compression = engine_kinematics(0.1,0.1,0.15,12,0,180);
p_compression = constant_c./v_compression.^gamma;

% State variables at state point 3
v3 = v2

% p3v3/t3 = p2v2/t2 | p3 = p2*t3/t2
p3 = p2*t3/t2

constant_c = p3*v3^gamma;
v_expansion = engine_kinematics(0.1,0.1,0.15,12,180,0);
p_expansion = constant_c./v_expansion.^gamma;

% State variables at state point 4
v4 = v1

% p3v3^gamma = p4v4^gamma | p4 = p3(v3/v4)^gamma
p4 = p3*(v3/v4)^gamma


figure(1)
hold on
plot(v1,p1,'*','color','r')
plot(v_compression, p_compression)
plot(v2,p2,'*','color','r')
plot([v2 v3],[p2 p3])
plot(v3,p3,'*','color','r')
plot(v_expansion, p_expansion)
plot(v4,p4,'*','color','r')
plot([v4 v1],[p4 p1])

%THERMAL EFFICIENCY
Thermal_eff = (1-(1/cr^(gamma-1)))*100

title('Air Standard Cycle P-V diagram')
xlabel('volume,V')
ylabel('pressure,P')
grid on


figure(2)
plot([v1 v2 v3 v4 v1], [p1 p2 p3 p4 p1], 'color', 'b')
xlabel('volume,V')
ylabel('pressure,P')

RESULT:-