Visualisation of OTTO Cycle using MATLAB and to calculate its Thermal Efficiency

AIR STANDARD CYCLE : OTTO CYCLE

Otto Cycle is an Constant Volume thermodynamic Cycle typically used in Spark Ignition Engines commonly found in Automobiles.

A typical Otto cycle consists of

• Pair of Isentropic Processes which are reversible Adiabatic.
• Pair of Isochoric Processes which are Constant Volume Cycles.

The typical PV diagram of the otto cycle is as follows,

• 0-1 :- Intake Stroke
• 1-2 :- Compression Stroke
• 2-3 :- Ignition Phase
• 3-4 :- Expansion Stroke
• 4-1 :- Heat Rejection Phase
• 1-0 :- Exhaust Stroke

Using OCTAVE, To Solve the Otto Cycle and to plot PV Diagram.

Let's Consider an example,

• In a SI Engine of compression ratio=9, pressure and temp at beginning of cycle 1 bar and 30C, also peak pressure being 40 bar
  (1bar = 100000 Pa , 1atm= 101325 Pa); gamma=1.4 since air. Calculate the Thermal Efficiency & plot PV Diagram.

OBJECTIVE :

To Find the State Variables of the Otto Cycle and to plot the PV Diagram using MATLAB program and calculate the Thermal Efficiency.

PROGRAM :

clear all
close all
clc

gamma=1.4;

% Engine Specs
bore=0.1;
stroke=0.1;
con_rod=0.15;
cr=9;

%state 1
p1=100000;
t1=303;

v_swept=(pi/4)*bore^2*stroke;
v_clr=v_swept/(cr-1);

v1=v_clr+v_swept;
v2=v_clr;

%state 2
%isentropic Compresssion p2v2^gamma=p1v1^gamma
%p2=p1*(v1/v2)^gamma=p1*cr^gamma
p2=p1*cr^gamma;

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

%state 3
v3=v2;
p3=4000000;

%(p3v3/t3)=(p2v2/t2)
t3=p3*v3*t2/(p2*v2);


%state 4
v4=v1;

% isentropic Expansion p3v3^gamma=p4v4^gamma
%p4=p3(v3/v4)^gamma
p4=p3*(v3/v4)^gamma;
t4=p4*v4*t3/(p3*v3);

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

The MATLAB program begin with 'clear all', 'close all' & 'clc' to clear and close all the previous data to avoid manipulation of our values.

All the input variables are defined,

since Air Cycle.

The input values contains the Engine Specification,

• Bore Dia (d) = 0.1m
• Stroke Length (l) =0.1m
• Connecting rod Length=0.15m
• Compression ratio (cr)= 9

State 1 variables are given in data,

• p1=100000 pa
• t1=303 K
• Swept Volume (v_swept) :- It is the volume swept by the piston in cylinder while moving between TDC & BDC.
• Clearance Volume (v_clr) :- It is the volume left in the cylinder above the piston, when the piston is at TDC.

• v1=v_clr+v_swept ie; the total volume.
• v2=v_clr

State 2 :- Since process 1-2 is Isentropic Compression ,

also,

Since,

State 3 :- Constant Volume Process

• v3=v2
• p3=4000000 ;given in data

State 4 :- Isentropic Expansion

• v4=v1

Equations for all the State Variables are defined, and when the program is run it gives output of the Variables

• To plot the State Variables,

All the 4 state Variables are marked with 'asterik' using plot comand & hold on comand is used to to add all the 4 plots together in a graph and to avoid the rewrite of plot comand over the other.

Now we have got the plot of state variables, joining them will give the PV diagram.

...continued
plot([v1 v2 v3 v4 v1],[p1 p2 p3 p4 p1])

using the above comand will result in a straight line joining the State Variables, which is not a perfect solution as the piston never moves linearly within the cylinder.

To overcome this, Piston Kinematics has to be considered,

The combustion or the volume change inside the Engine cylinder follows the Piston kinematics and is defined by equation

`V/V_C = 1+1/2(cr-1)[R+1-cos theta-(R^2-sin^2theta)^(1/2)]`

CODE FOR PISTON KINEMATICS

clear all
close all
clc

bore=0.1;
stroke=0.1;
con_rod=0.15;
cr=9;

v_swept=(pi/4)*bore^2*stroke;
v_clr=v_swept/(cr-1);

a=stroke/2;
R=con_rod/a;

theta=linspace(0,720,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_clr;

plot(theta,v,'linewidth',2)

This new program contains the same engine spec inputs

• R is the ratio of Connecting rod length to Cranking radius
• a is cranking radius;a=stroke/2
• theta is defined from an interval of (0 to 720) degrees as it is 4-stroke cycle and each 180 degree represents a stroke in an ideal otto cycle as shown in fig.

The equation is divided into 3 different terms for the coding feasiblity

Since theta is an array, Element-wise operator or Dot Product is used in the Equation and its terms.

plot comand is used to plot 'theta' values against 'v' , 'linewidth' defines thickness of plot.

When this program is executed, we get the plot of piston movement from TDC to BDC and hence back as shown below which proves that piston movement is not linear.

Now, calling this 2nd program (engine_kinematics.m) as a function in our previous program (otto_cycle.m) with the constant values defined, will plot the actual PV Diagram.

engine_kinematics.m (modified)

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

v_swept=(pi/4)*bore^2*stroke;
v_clr=v_swept/(cr-1);

a=stroke/2;
R=con_rod/a;

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_clr;

end

otto_cycle.m (final)

clear all
close all
clc
gamma=1.4; %data

% Engine Specs
bore=0.1;
stroke=0.1;
con_rod=0.15;
cr=9;

%state 1
p1=100000; %data
t1=303; %data

v_swept=(pi/4)*bore^2*stroke;
v_clr=v_swept/(cr-1);

v1=v_clr+v_swept;
v2=v_clr;

%state 2
%isentropic Compresssion p2v2^gamma=p1v1^gamma
%p2=p1*(v1/v2)^gamma=p1*cr^gamma
p2=p1*cr^gamma;

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

constant_c=p1*v1^gamma;
v_compression=engine_kinematics(bore,stroke,con_rod,cr,180,0);
p_compression=constant_c./v_compression.^gamma


%state 3
v3=v2;
p3=4000000; %data

%(p3v3/t3)=(p2v2/t2)
t3=p3*v3*t2/(p2*v2);

constant_c=p3*v3^gamma;
v_expansion=engine_kinematics(bore,stroke,con_rod,cr,0,180);
p_expansion=constant_c./v_expansion.^gamma

%state 4
v4=v1;

% isentropic Expansion p3v3^gamma=p4v4^gamma
%p4=p3(v3/v4)^gamma
p4=p3*(v3/v4)^gamma;
t4=p4*v4*t3/(p3*v3);

%to find eficiency of otto cycle engine
efficiency=1-(v2/v1)^(gamma-1)

%plotting
hold on
plot(v1,p1,'*','color','r')
plot(v_compression,p_compression,'color','g','linewidth',2)
plot(v2,p2,'*','color','r')
plot([v2 v3], [p2 p3],'color','r','linewidth',2)
plot(v3,p3,'*','color','r')
plot(v_expansion,p_expansion,'color','b','linewidth',2)
plot(v4,p4,'*','color','r')
plot([v4 v1], [p4 p1],'color','r','linewidth',2)
xlabel('Volume')
ylabel('Pressure')

• v_compression, p_compression are used to plot the compression curve by using the function defined in engine_kinematics.m program
• the function is called out to calculate v_compression where parameters remains same for engine but the start_crank is 180 and end_crank is 360 or 0.
• p_compression is find out by using the constant_c & v_compression.
• Similarly v_expansion & p_expansion are calculated by calling out the function again for another set of values from 0 to 180 degrees.
• Efficiency of the cycle is also calculated within the program by using formula

• Finally series of plot comands with hold on comand are used in order to get the set of 4 State variables, compression curve, expansion curve and 2 parallel Constant Volume curves with different 'color' and thickness of plot is defined by 'linewidth'.

Result

The PV Diagram for the Otto cycle is plotted and the Efficiency is found to be 58.47%