MATLAB code to plot an Air standard cycle (OTTO)

Aim: To write a program in MATLAB to plot p-v diagram of Otto cycle.

Governing equations:

Otto cycle is a gas power cycle that is used in spark-ignition internal combustion engines . This cycle was named after a German Engineer - Dr. Nikolaus Otto.

An Otto cycle consists of four processes:

1. Two isentropic (reversible adiabatic) processes

Process 1 - 2 : isentropic compression

Process 3 - 4 : isentropic expansion

2. Two isochoric (constant volume) processes

Process 2 - 3 : constant volume heat addition

Process 3 - 4 : constant volume heat rejection

These processes can be easily understood if we understand p-V (Pressure-Volume) diagram of Otto cycle as shown in thev above figure.

Consider a piston-cylinder with bore size (B) of 0.1m , stroke (L) of length 0.1m, connecting rod of length 0.15m, compression ratio (cr) of 7.

Swept volume is given by Vs = `π/4 B^2 L` -------------(1)

Clearance volume Vc is calculated as-

`cr = (V1 )/ (V2)`

V1 = Vclearance + Vswept ,

V2 = Vclearance ,

`cr = (Vc + Vs)/(Vc)`

`cr = 1 + ( Vs)/(Vc)`

`cr -1 = ( Vs)/(Vc)`

`Vc = ( Vs)/(cr-1)` -------------(2)

Now initially we know the temperature, pressure and volume at point 1.

At point 1:

T1 = 500 K

P1 = 101325 Pa

V1 = Vclearance + Vswept

At point 2:

V2 = Vclearance

also for process 1 - 2 we know that it follows the relation `P1〖V1〗^γ=P2〖V2〗^γ`

`P2=P1* 〖cr〗^γ` -------------(3)

Also `(T2)/(T1)=〖(P2)/(P1)〗^((γ-1)/γ)= 〖(V1)/(V2)〗^(γ-1)`

`T2= 〖cr〗^(γ-1) . T1` -------------(4)

At point 3:

V3 = V2

`P3=P2* (T2)/(T3)` -------------(5)

At point 4:

V4 = V1

`P4 = (P3)/(〖cr〗)^γ` -------------(6)

Along with this the volume of the piston cylinder can be determined as a function of crank angle from the compression ratio, the stroke, bore and connecting rod length. The geometric parameters of the piston cylinder are represented in figure below

The relation is given by

`V/(Vc)=1+0.5(rc-1)( (R+1-cosθ+〖(R^2-sinθ)〗^0.5)` -------------(7)

Objective of the project:

1. To plot p-v diagram for otto cycle using matlab code.
2. To calculate the thermnal efficiency of the cycle

Code:

Code of Function to be called later:

function [V] = engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank)
a = stroke/2;
R = con_rod/a;

Vs = (pi/4)*bore^2*stroke;
Vc = Vs/(cr-1);

% crank angle is divide into 100 equal divisions between start crank angle and end crank angle
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)).* Vc;

end

Main program code:

%a program for plotting otto cycle
clear all
close all
clc

%inputs pressure in pascal, temperature in kelvin 
P1 = 101325;
T1 = 308;
T3 = 1100;
gamma = 1.4;

%Geometry of piston and cylinder ( length in meters )
bore = 0.1;
stroke = 0.1;
cr = 7;
con_rod = 0.15;

V_swept = (pi/4)*bore^2*stroke;
V_clearance = V_swept/(cr-1);

%thermal efficiency of otto cycle is given by this formula
thermal_efficiency = (1 - 1/cr^(gamma-1))*100; 

%state variables at point 1
V1 = V_swept + V_clearance;

constant_c1 = P1*V1^gamma; % as PV^gamma = constant
% calling function to get the variation of volume at each and every instant
V_compression = engine_kinematics(bore,stroke,con_rod,cr,180,0);
P_compression = constant_c1./V_compression.^gamma; 

%state variables at point 2
V2 = V_clearance;
P2 = P1 * (cr)^gamma;
T2 = cr^(gamma-1) * T1;

%state variables at point 3
V3 = V2;
P3 = P2 * (T2/T3);

constant_c2 = P3*V3^gamma;  % as PV^gamma = constant
% calling function to get the variation of volume at each and every instant
V_expansion = engine_kinematics(bore,stroke,con_rod,cr,0,180);
P_expansion = constant_c2./V_expansion.^gamma;

%state variables at point 4
V4 = V1;
P4 = P3/cr^gamma;

%plotting
figure(1);
hold on;
plot(V1,P1,'*','color','r');          % to get point 1 on plot
plot(V_compression,P_compression);    % this gives the exact variation during process 1-2
plot([V2 V3],[P2 P3],'color','r'); 
plot(V2,P2,'*','color','r');          % to get point 2 on plot
plot(V3,P3,'*','color','r');          % to get point 3 on plot
plot(V_expansion,P_expansion);        % this gives the exact variation during process 3-4
plot([V4 V1],[P4 P1],'color','r');
plot(V4,P4,'*','color','r');          % to get point 3 on plot
title('P-V diagram of Otto cycle');
xlabel('Volume (m^3)');
ylabel('Pressure (pascal)');

Explanation of code:

1. To begin with the parameters known to us like bore diameter, stroke length, connecting rod length, compression ratio, pressure and temperature at beginning of compression, temperature at the end of heat addition and the value of isentropic index (gamma) was assigned appropriate value using variables as used in the codes.
2. Then the swept and clearance volume was calculated using formula derived in the governing equation section.
3. The thermal efficiency of otto cycle was calculated using the given parameters.
4. The only unknown state variable at point 1 was volume which is now calculated.
5. Now for the compression process 1-2 and expansion process 3-4, to track the exact path of otto cycle variation of volume with respect to crank angle was calculated by creating function engine kinematics which was called here to calculate the value of volume at different crank angle from to (during compression) and to (during expansion).
6. In the function [V] = engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank), the formula was divided into 3 terms for simplicity. The equation (7) was divided

As term1 = 0.5*(cr-1);

term2 = R + 1 - cosd(theta);

term3 = (R^2 - sind(theta).^2).^0.5

Also the theta between to was divided into 100 equal divisions using linspace command.

6. Temperature , pressure and volume are calcukated using equations (3), (4), (5) and (6) discussed in governing equations section
7. Finally plotting command was used to plot p-v diagram

Errors Faced:

The code was executing and showing results in workspace but I was getting this kind of plot. So after going through the code again I found out one simple mistake as shown below:

In line number 25, I mistakenly used V1 = V_swept but as we know it should be V1 = V_swept + V_clearance. Once i made this correction then the problem was solved. Apart from this error was in the use of elementwise operator in engine kinematics function which was solved easily.

Result:

1. The p-v diagram for otto cycle was plotted using MATLAB code.
2. The use of function in MATLAB was studied and implemented.
3. The thermal efficiency of otto cycle is about 54 % as shown in the workspace.