Solving Otto Cycle using MatLab

OBJECTIVE:

To write a MatLab program that can solve the Otto Cycle, generate a PV diagram for the air standard cycle and output the thermal efficiency.

OVERVIEW:

Air standard cycles are reference cycles that give an appropriation to the performance of internal combustion engines. Otto cycle is one of the commonly used air standard cycles.

Otto Cycle corresponds to the four-stroke petrol or gasoline-powered engines also called as Spark Ignition (SI) engines. It has a piston-cylinder arrangement. The piston is connected to the crankshaft by a connecting rod. The rotational motion of the crankshaft makes the piston slide up and down in the cylinder. Initially, the piston is at the top dead center of the cylinder.

First, the crankshaft rotates from 0° to 180°, and as it rotates, the piston slides down from the top dead center(TDC) to the bottom dead center (BDC), while the air-fuel mixture enters the cylinder.

Second, as the crankshaft rotates from 180° to 360°, the piston compresses the air-fuel mixture in the cylinder. This process is called Isentropic compression. When the piston is at the TDC, the air-fuel mixture is compressed and the pressure inside the cylinder is increased, at this point, the mixture is ignited by a spark plug. The mixture explodes (combustion of the mixture). This step is called the Constant volume heat addition.

Due to the combustion of the mixture, the piston is pushed down making a 360° to 540° rotation of the crankshaft. This step is the Isentropic expansion. After the combustion, the temperature and the pressure inside the cylinder lowers, making this Constant volume heat rejection.

PROGRAM CODE:

Otto cycle program.

%OTTO CYCLE
clear all
close all
clc

%Inputs
gamma = 1.4;
t3 = 2300;

%State variables at state point 1 (in Pa,K,m^3)
p1 = 101325; % 1 atm = 101325 Pa
t1 = 500;

%Engine geometric parameters (in meters)
bore = 0.073;
stroke = 0.0715;
con_rod = 0.1;
r = 11;

%Calculating Swept Volume and Clearance Volume (in m^3)
v_s = (pi/4)*bore^2*stroke;
v_c = v_s/(r-1);

v1 = v_s + v_c; %Volume when piston is at BDC
v2 = v_c; %Volume when piston is at TDC

%State variables at state point 2
%p2*v2^gamma=p1*v1^gamma | p2 = p1*(v1/v2)^gamma
p2 = p1*r^gamma;

%Ideal gas law
%p1*v1/t1 = p2*v2/t2
t2 = p2*v2*t1/(p1*v1);

%Calculating compression volume and pressure
const_c = p1*v1^gamma;
V_comp = engine_kinematics(bore,stroke,con_rod,r,180,360);
P_comp = const_c./V_comp.^gamma;

%State variables at state point 3
v3 = v2;

%p3*v3/t3 = p2*v2/t2
p3 = p2*t3/t2;

%Calculating expansion volume and pressure
const_c = p3*v3^gamma;
V_exp = engine_kinematics(bore,stroke,con_rod,r,360,540);
P_exp = const_c./V_exp.^gamma;

%State variables at state point 4
v4 = v1;

%p3*v3^gamma = p4*v4^gamma
p4 = p3*(v3/v4)^gamma;

%p1*v1/t1 = p4*v4/t4
t4 = p4*t1/p1;

Efficiency_Otto_cycle = 1 - ((t4-t1)/(t3-t2))

%P-V plot
figure(1)
hold on
plot(v1,p1,'*','color','k');
plot(V_comp,P_comp,'linewidth',2,'color',[ 0.91 0.41 0.17]);
plot(v2,p2,'*','color','k');
plot([v2 v3],[p2 p3],'linewidth',2,'color','r');
plot(v3,p3,'*','color','k');
plot(V_exp,P_exp,'linewidth',2,'color','b');
plot(v4,p4,'*','color','k');
plot([v4 v1],[p4 p1],'linewidth',2,'color',[0 0.5 0.8]);
xlabel('Volume (in m^3)');
ylabel('Pressure (in Pa)');
title('P-V diagram of Otto Cycle');

Volume trace function program.

function [V] = engine_kinematics(bore,stroke,con_rod,r,crank_start,crank_end)

a = stroke/2; %Crank pin radius
R = con_rod/a; 

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

theta = linspace(crank_start,crank_end,180);

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

V = (1 +term1*(term2-term3)).*v_c; %Volume Trace

end

EXPLANATION:

Otto cycle coding:

Below are the generic variables used in the MatLab program.

bore - bore diameter of the cylinder (in meters)

stroke - height of the cylinder (in meters)

con_rod - length of the connecting rod (in meters)

r - compression ratio

p1, v1, t1 - state variables at state point 1 (in Pa, m^3, K)

p2, v2, t2 - state variables at state point 2 (in Pa, m^3, K)

p3, v3, t3 - state variables at state point 3 (in Pa, m^3, K)

p4, v4, t4 - state variables at state point 4 (in Pa, m^3, K)

v_s - Swept volume is the volume of the cylinder from the TDC to BDC (in m^3)

v_c - Clearance volume is the volume of the cylinder when the piston is at the TDC (in m^3)

Initially, the input variables namely the State variables at state point 1 are defined and values are assigned.

The Swept volume is calculated by defining the variable v_s as the area of the cylinder in terms of bore and stroke. v_s = pi/4*bore^2*stroke

The Clearance volume is derived from the compression ratio formula and is assigned to the variable v_c = v_s/(r-1).

Volume at state point 1 (v1) is the sum of Swept volume and Clearance volume. v1 = v_s + v_c

Volume at state point 2 (v2) is the volume of the cylinder at TDC. Hence, v2 = v_c

%OTTO CYCLE
clear all
close all
clc

%Inputs
gamma = 1.4;
t3 = 2300;

%State variables at state point 1 (in Pa,K,m^3)
p1 = 101325; % 1 atm = 101325 Pa
t1 = 500;

%Engine geometric parameters (in meters)
bore = 0.073;
stroke = 0.0715;
con_rod = 0.1;
r = 11;

%Calculating Swept Volume and Clearance Volume (in m^3)
v_s = (pi/4)*bore^2*stroke;
v_c = v_s/(r-1);

v1 = v_s + v_c; %Volume when piston is at BDC
v2 = v_c; %Volume when piston is at TDC

State point 2 calculation is shown below.

p2 is calculated using the relation between the pressure, volume and specific heat ratio (gamma 'γ') at state points 1 & 2.

t2 is calculated using the Ideal gas law with respect to the state points 1 & 2.

%State variables at state point 2
%p2*v2^gamma=p1*v1^gamma | p2 = p1*(v1/v2)^gamma
p2 = p1*r^gamma;

%Ideal gas law
%p1*v1/t1 = p2*v2/t2
t2 = p2*v2*t1/(p1*v1);

State point 3 calculation is shown below.

v2 = v3 (∵ volume is constant at State points 2 &3).

p3 is calculated using the Ideal gas law with respect to the state points 2 &3.

%State variables at state point 3
v3 = v2;

%p3*v3/t3 = p2*v2/t2
p3 = p2*t3/t2;

State point 4 calculation is shown below.

Similarly, volume at state points 4 & 1 is constant. ∴ v4 = v1.

p4 calculated using the relation between the pressure, volume and specific heat ratio (gamma 'γ') at state points 3 & 4.

t4 is calculated using the Ideal gas law with respect to the state points 3 & 4.

%State variables at state point 4
v4 = v1;

%p3*v3^gamma = p4*v4^gamma
p4 = p3*(v3/v4)^gamma;

%p1*v1/t1 = p4*v4/t4
t4 = p4*t1/p1;

Otto Cycle efficiency calculation is shown below.

Now that the temperatures of all the state points are calculated, the efficiency of the Otto cycle can be computed using the formula η = 1−(t4−t1)/(t3−t2). This will be printed as output as the semicolon is not used in this line.

Efficiency_Otto_cycle = 1 - ((t4-t1)/(t3-t2))

Function coding:

A function program is created separately and saved in the same folder as the main program. This is done to optimize the program and also the same function can be used in some other similar program as well. This function is used to calculate several values of the volume during the compression and expansion process of the Otto cycle.

The function is named 'engine_kinematics', it is stored in an array called '[V]'.

Different input variables that are needed to compute [V] are the bore, stroke, con_rod, r, crank starting angle and crank ending angle. And varied values of [V] are calculated by using the linspace() command to the angle (theta) as shown in the coding below.

The formula to calculate the volume 'V' is derived from the compression ratio formula:

V = [1+(0.5*(r-1))*(R+1-cosθ-(R^2-sin^2θ)^0.5)]*v_c

R = con_rod/a

where,

θ - crank angle made by the crankshaft with the line normal to the piston.

a - crank radius

r - compression ratio

R - the ratio of the length of the connecting rod to the crank radius

function [V] = engine_kinematics(bore,stroke,con_rod,r,crank_start,crank_end)

a = stroke/2; %Crank pin radius
R = con_rod/a; 

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

theta = linspace(crank_start,crank_end,180);

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

V = (1 +term1*(term2-term3)).*v_c; %Volume Trace

end

Compression volume values are calculated by calling the function 'engine_kinematics' and storing them in the array 'V_comp'. Using each value from the V_comp array, the compression pressure 'P_comp' is computed. This compression volume & pressure function calling program is typed after the 'State variables at state point 2' computation part.

%Calculating compression volume and pressure 
const_c = p1*v1^gamma;
V_comp = engine_kinematics(bore,stroke,con_rod,r,180,360);
P_comp = const_c./V_comp.^gamma;

Expansion volume values are calculated by calling the function 'engine_kinematics' and storing them in the array 'V_exp'. Using each value from the V_exp array, the expansion pressure values 'P_exp' are computed. This expansion volume & pressure function calling program is typed after the 'State variables at state point 3' computation part.

%Calculating expansion volume and pressure
const_c = p3*v3^gamma;
V_exp = engine_kinematics(bore,stroke,con_rod,r,360,540);
P_exp = const_c./V_exp.^gamma;

Plotting:

The Pressure-Volume diagram is generated using the 'plot() command'.

• The coordinates of the State point 1,2,3 & 4 are plotted in black color using the values (v1,p1), (v2,p2), (v3,p3) & (v4,p4) in four separate 'plot() commands' as shown in the program snippet below.
• The compression and expansion curves are generated to connect the points 1 to 2 and 3 to 4 using the arrays V_comp, P_comp, V_exp & P_exp in two separate 'plot() commands' as shown in the program snippet below.
• The constant volume heat addition &constant volume heat rejection lines are generated to connect the points 2 to 3 & 4 to 1 using the points (v2,v3), (p2,p3), (v4,v1) & (p4,p1) in two separate 'plot() commands' as shown in the program snippet below.
• The points, curves and lines generation plot() commands are arranged in the order given below.
• The X & Y axes are labeled as 'Volume (in m^3)' and 'Pressure (in Pa)' using xlabel() and ylabel() commands as shown in the program snippet below.
• A title is given to the P-V Diagram using the title() command.
• The 'hold on' command is used so that all the plot commands are executed in the same figure.

%P-V plot
figure(1)
hold on
plot(v1,p1,'*','color','k');
plot(V_comp,P_comp,'linewidth',2,'color',[ 0.91 0.41 0.17]);
plot(v2,p2,'*','color','k');
plot([v2 v3],[p2 p3],'linewidth',2,'color','r');
plot(v3,p3,'*','color','k');
plot(V_exp,P_exp,'linewidth',2,'color','b');
plot(v4,p4,'*','color','k');
plot([v4 v1],[p4 p1],'linewidth',2,'color',[0 0.5 0.8]);
xlabel('Volume (in m^3)');
ylabel('Pressure (in Pa)');
title('P-V diagram of Otto Cycle');

Below is the workspace generated upon the execution of the Otto Cycle MatLab program.

OUTPUT:

1. The efficiency output of the Otto Cycle output.

2. P-V Diagram of the Otto Cycle.