Simulation of Air Standard Otto Cycle using MATLAB

1. AIM: To write a MATLAB programme for simulation of Air standard Otto cycle.

2. INTRODUCTION:

• Otto cycle consist four processes, in which two are constant volume (isochoric) and remaining two are isentropic processes. In this cycle heat addition is done at constant volume process so this cycle is also known as constant volume cycle.
• We know that actual analysis is complicated so we use some assumption to analysis the petrol engine. The main assumption is that the working fluid is ideal air and heat is transfer by heat reservoir so this cycle known as air standard cycle.
• This cycle was given by Nicolas Otto so known as Otto cycle.

3. GOVERNING EQUATIONS:

3.1 Diagram & Details of Processes in Otto cycle

Fig. 1: p-V diagram of Otto cycle

Let,

At the beginning of compression i.e. at state point 1

P₁, V₁, T₁ = Pressure (N/m²), Volume (m³) & Temperature (K)

At the end of compression stroke i.e. at state point 2

P₂, V₂, T₂ = Pressure (N/m²), Volume (m³) & Temperature (K)

At the beginning of expansion stroke i.e. at state point 3

P₃, V₃, T₃ = Pressure (N/m²), Volume (m³) & Temperature (K)

At the end of expansion stroke i.e. at state point 4

P₄, V₄, T₄ = Pressure (N/m²), Volume (m³) & Temperature (K)

3.2 Initial engine geometric Parameters

1. Swept volume: It is the volume between top dead center (TDC) and bottom dead center (BDC) the volume between top dead center (TDC) and bottom dead center (BDC).

2. Clearance volume: Clearance volume is a volume between the cylinder head and the piston top when the piston is at top dead center (TDC)

3. v₁= Swept volume + Clearance volume

4. v₂ = Clearance volume

5. Equation for change in volume at changing crank angle:

3.3 Otto cycle consists four processes which are as follows:

1. Isentropic Compression Process (1–2): At pt. 1 cylinder is full of air with volume V₁, pressure P₁ and temp. T₁. Piston moves from BDC to TDC and an ideal gas (air) is compressed isentropically to state point 2 through compression ratio (cr).

2. Constant Volume Heat Addition Process (2 – 3): Heat is added at constant volume from an external heat source. The pressure rises and the ratio is called expansion ratio or pressure ratio.

3. Isentropic Expansion Process (3 – 4): The increased high pressure exerts a greater amount of force on the piston and pushes it towards the BDC.

4. Constant Volume Heat Rejection Process (4 – 1): Heat is rejected to the external sink at constant volume. This process is so controlled that ultimately the working fluid comes to its initial state 1 and the cycle is repeated.

3.4 Efficiency of Otto Cycle.

4. OBJECTIVES OF THE PROJECT:

1. To determine thermodynamic properties such as pressure, temperature and volume at state 1, 2, 3 & 4 of air standard Otto cycle.
2. To plot PV diagram using state variable as well as using equation for change in volume at different crank angle.
3. To estimate thermal efficiency of Otto cycle.

5. MAIN PROGRAM:

Software used: MATLAB R2020a trial version

Following is the program used to simulate Otto cycle and estimate thermal efficency of the same.

Program includes comments which explains commands used in program and various other parameters.

% To write a MATLAB programme for simulation of Air standard Otto cycle

clear all 
close all
clc

%Inputs
gamma = 1.4                                 % Ratio of specific heat (Cp/Cv)

% Given state varaibles
p1 = 101325                                 % Pressure at the begining of compression stroke (N/m^2)
T1 = 500                                    % Temperature at the begining of compression stroke (K)
T3 = 2300                                   % Temperature at the end of compression stroke (K)

% Engine geometric parameters
bore = 0.1;                                 % Bore of cylinder (m)
stroke = 0.1;                               % Length of stroke of cylinder (m) 
con_rod = 0.15;                             % Length of connecting rod (m)
cr = 12;                                    % Compression ratio (unit less)

% Calculating the swept volume and the clearance volume
v_swept = (pi/4)*bore^2*stroke;             % Swept volume: Volume travel by piston in single stroke (m^3)
v_clearance = v_swept/(cr-1);               % Clearance volume: Volume b/w cylinder & piston head at TDC (m^3)

v1 = v_swept + v_clearance                  % Total volume of cylinder as shwon in figure 1.
v2 = v_clearance                            % Volume of cylinder when piston is at TDC as shwon in figure 1

% State varaiable at state point 2
% p2v2^gamma = p1v1^gamma
% p2 = p1(v1/v2)^gamma = p1(cr)^gamma
p2 = p1*cr^gamma;                           % Pressure at state 2 i.e. at the end of compression stroke (N/m^2)

% p1v1/T1 = p2v2/T2 | T2 = p2*v2*T1/(p1*v1)
T2 = p2*v2*T1/(p1*v1);                      % Temperature at the begining of const. vol. heat addition (K)

constant_c = p1*v1^gamma;
V_compression = engine_kinematics(bore, stroke, con_rod, cr,0 ,180); % compression 180 to 0 degree

%P_compression*V_compression^gamma = constant_c
P_compression = constant_c./V_compression.^gamma;
% State varaibale at state point 3
v3 = v2                                     % Since process 2 to 3 is a constant volume process

% p3v3/T3 = p2v2/T2 | p3 = p2*T3/T2
p3 = p2*T3/T2;                              % Pressure at state 3 i.e. at the end of constant volume process (N/m^2)

%V_heat_addition = engine_kinematics(bore,stroke,con_rod,cr,0,0);
V_heat_addition = linspace(v2,v3,180);
P_heat_addition = linspace(p2,p3,180);

constant_c = p3*v3^gamma;
V_expansion = engine_kinematics(bore, stroke, con_rod, cr,0 ,180);      % expansion 0 to 180 degree

%P_expansion*V_expansion^gamma = constant_c
P_expansion = constant_c./V_expansion.^gamma;
% State varaibale at state point 4
v4 = v1                                     % Since process 4 to 1 is a constant volume process

% p3v3^gamma = p4v4^gamma | p4 = p3* (v3/v4)^gamma
p4 = p3* (v3/v4)^gamma                      % Pressure at state 4 i.e. at the end of expansion process (N/m^2)

% State variable for heat rejection process (4-1)
V_heat_rejection = linspace(v4,v1,180);
P_heat_rejection = linspace(p4,p1,180);

% plotting state varaiables
figure(1)                                       % to create figure 1
hold on                                         % to avoid over write of previous plot by recent plot command
plot(v1,p1,'.','color','r')                     % to plot state 1 on figure (1)
plot(V_compression,P_compression,'color','r')
plot(v2,p2,'.','color','r')                     % to plot state 2 on figure (1) 
plot(V_heat_addition,P_heat_addition)
plot(v3,p3,'.','color','r')                     % to plot state 3 on figure (1)
plot(V_expansion,P_expansion,'color','g')
plot(v4,p4,'.','color','r')                     % to plot state 4 on figure (1)
plot(V_heat_rejection,P_heat_rejection)
xlabel('Volume (m^3)'); 
ylabel('Pressure (N/m^2)'); 
title('PV Diagram for Otto Cycle');

% labeling the figure
text(v1,p1, '1','fontsize',6)                   % to represent state number on PV diagram
text(v2,p2, '2','fontsize',12)
text(v3,p3, '3','fontsize',12)
text(v4,p4, '4','fontsize',6)
text(V_compression(100),P_compression(100),'Isentropic Compression','fontsize',12,'color','r')
text(V_expansion(120),P_expansion(120),'Isentropic Expansion','fontsize',12, 'color','g')
grid on

% Thermal Efficiency of otto Cycle
eta = (1 - 1/((cr)^(gamma-1)))*100

Program to plot state of PV digram and connect those states: To connect state points of PV diagram a program is written from the graph given below we can observe that state points 1, 2, 3 & 4 are connected by straight line which is not a feasible solution.

figure (2)                                                % to create figure 2

plot([v1,v2,v3,v4,v1],[p1,p2,p3,p4,p1],'color','g')       % to plot pv-diagram all four states are connected by line.

xlabel('Volume (m^3)'); 
ylabel('Pressure (N/m^2)'); 

title('PV-Diagram for Otto Cycle');

grid on

% labeling the figure

text(v1,p1, '1','fontsize',6)                               % to represent state number on PV diagram
text(v2,p2, '2','fontsize',12)
text(v3,p3, '3','fontsize',12)
text(v4,p4, '4','fontsize',6)

grid on

Figure No. 02 : PV diagram for Otto cycle (Not a correct representation)

from figure no. 02 we can observe that from state 1 to state 2 the varaition of volume with respect to pressure is linear also from state 3 to 4 the varaition of volume with respect to pressure is linear. In actual PV digram from state 1 to 2 i.e. Isentropic compression & from state 3 to 4 i.e. Isentropic expansion p*v^gamma= constant equation follows so can clearly conclude that above PV digram is not a feasible representation of PV diagram of Air standard Otto cycle.

Function used in Main Program: In above program variation of volume with respect to geometric parameters of engine not considered so we were getting straight line connecting the state points 1,2,3 & 4.To track variation volume during compression and expansion stroke with respect to crank angle, function named "engine_kinematics" is used.

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

a = stroke/2;                                   % crank pin radius
R = con_rod/a;                                  % Ratio of connecting rod length to crank pin radius

% Calculating the swept volume and the clearance volume
v_swept = (pi/4)*bore^2*stroke;                 % Swept volume: Volume travel by piston in single stroke (m^3)
v_clearance = v_swept/(cr-1);                   % Clearance volume: Volume b/w cylinder & piston head at TDC (m^3)

theta = linspace(start_crank,end_crank,720);    % Creating an array of theta from 0 to 720 for 720 values
                                                % When theta is zero V corresponds TDC volume i.e. Clearance volume

term1 = 0.5*(cr-1);
term2 = R +1-cosd(theta);               
term3 = (R^2 - sind(theta).^2).^0.5;            % elemnetwise operator as theta is an array

V = (1 + term1*(term2-term3)).*v_clearance;
end

from above program we get the variation of volume with engine geometric parameter

For crank angle 0 to 180 degree

Figure No. 03: Varaiation of volume with cank angle (0 to 180 degree)

For crank angle 0 to 720 degree

Figure No. 04 : Varaiation of volume with cank angle (0 to 720 degree)

From figure no. 03 & 04 it is evident that with the variation of crank angle volume varaition is not linear. Therefore, while plotting Isentropic compression processes & Isentropic expansion process on PV digram variation of volume with crank angle must be considered.

Output: After considering variation of volume with respect to crank angle and to track variation during compression and expansion stroke function named "engine_kinematics" is called in main program which gives the following outputs.

Output State Vaiable: Values of thermodynamic properties of states i.e. pressure, temperarure & volume at state 1, 2, 3 & 4.

Figure No.05: Sate parameters & Efficiency of Otto cycle

Output-- PV Diagram:

Figure No. 06: PV Digram for Otto Cycle

Error occured during formulation of Program:

Error 1: Pressure at state 1 is represented by p1 but in program it was typed as P1.

Error 2: Compression ratio in the program is represented by cr but in program it was typed as rc.

CONCLUSIONs: Following conclusions can be made from the above program

1. In order to construct the pv-diagram completely we must know how the volume changes from state 1 to state 2 i.e. during Isentropic compression process and from state 3 to state 4 i.e. during Isentropic expansion process. This variation actually depends on engine kinematics.
2. With the help of engine parameters kinematic equation of a piston can be solved and volume trace can be computed. From this volume trace, swept volume and clearance volume and ultimetly v₁ & v₂ can be determined.
3. If engine parameters not considered in calculations thanvaraiation of pressure and volume will be linear which is not in line with the actual PV diagram.