1. Objective:-

To write code in Matlab and that can solve an otto cycle and Plot PV Diagram.
i. The code should create a PV diagram.
ii. Finding the thermal efficiency of the engine.

2. Basic:-

2.1 Ideal Otto Cycle:-
Air standard Otto cycle is used for the SI Engines (Spark Ignition Engine). The cycle Contains Mainly four Processes. In the overall cycle, suction and Exhaust are neglected, also power consumption is not considered.
The Ideal Otto cycle consists of the four processes which are as below.

1–2 Reversible Adiabatic compression (Or) Isentropic Compression.
2–3 Constant Volume Heat Addition.
3–4 Reversible Adiabatic Expansion (Or) Isentropic Expansion.
4–1 Constant Volume Heat Rejection.

When the piston in the cylinder of the Internal Combustion Engine (SI Engine) moves from BDC (i.e Bottom Dead Centre) to TDC (i.e Top Dead Centre) the Pressure and volume Relatively change, which is represented on the Diagram known as PV-Diagram. The ideal PV-Diagram can be seen below.

The terms related to the otto cycle is as discussed below.

2.2 Terms in the Otto cycle:
As we know the basic otto cycle from above, there are two terms in induced TDC and BDC. Where TDC is called Top dead Centre which is Extreme position of the piston at the top end of cylinder similarly BDC is Extreme position at the bottom end of the cylinder. other important terms related to the Otto cycle is Explain further.

i. Compression ratio (cr):- It is the ratio of volume before compression to the volume after compression.

Compression ratio (cr) = Volume before compression/volume after compression

also, Compression ratio (cr) = (Vs+Vc)/(Vc)

Where:-

V1= volume before compression = V_swept + V_clearance

V2 = volume after compression = V_clearance

Vs = V_swept , Vc = V_clearance

Therefore, Compression ratio (cr) = (V1)/(V2)

ii. Expansion ratio (er):- It is the ratio of volume after Expansion to the volume before compression.

Similarly, as the compression ratio, we can get an Expansion ratio.

Therefore, Expansion ratio (er) = (V4)/(V3)

Note:- for Otto cycle compression Ratio = Expansion ratio i.e (cr=er).

iii. Swept volume and Clearance volume:- Swept Volume (V_s) is a volume Before compression stroke. Also known as stoke volume and Displacement volume. Similarly, Clearance volume (V_c) is a volume after compression stroke.

3. Defining State Variable:-

3.1 State variable at Point-1:

i. Initial Pressure & Temperature:

P1= 100 Kpa, T1= 450 K

ii. Engine Geometry parameters:

Bore (D) = 0.1 m

Stroke (L) = 0.1 m

Connceting rod Length (con_rod) = 0.15 m

And Compression Ratio (cr) = 9 and 8

iii. Swept volume and Clearance volume:

Clearance volume V_c = V_s/(cr-1)

Swept Volume V_s = (pi/4)*D²*L

Where:-
cr = compression ratio, D = Bore and L = Stroke

3.2 State Variable at point-2:

i. Calculating Pressure at point 2 (P2):

Due to isentropic compression, we can apply the PV Relation here.

P1V1^gamma = P2V2^gamma

Where, P2 is Unkown to us which can be calculated with the help of the relation between the volume, pressure and compression ratio

Also, we can calculate the P2 as follow

Cr = V1/V2 = (P2/P1)^gamma

Therefore, P2 = P1 * ( cr )^gamma

ii. Calculating temp at point 2 (T2):

To get the value of temp T2 we are using the ideal gas equation as follows.

P1V1 = mRT1 , mR = P1V1/T1

P2V2 =mRT2 , mR = P2V2/T2

Therefore, P1V1/T1 = P2V2/T2

Hence, T2 = P2*V2*(T1)/(P1*V1)

3.3 State Variable at point-3:

i. Calculating Volume at point 3 (V3):

V3=V2

ii. Calculating Pressure at point 3 (P3):

By using the ideal gas equation

P3 = P2 * (T3)/(T2)

3.4 State Variable at point-4:

i. Calculating Volume at point 4 (V4):

V1= V4

ii. Calculating Pressure at point 3 (P3):

Also from P3V3^gamma = P4V4^gamma

P4 = P3 *((V3)/(V4))^gamma

4. Detailed Procedure and Explanation:-

4.1 Procedure and Explanation for Matlab:

Step-1: Initial Program

After Defining all the state variable we are writing the first program as follows

close all
clear all
clc% Inputs
 gamma= 1.4
 T3=2778
 
 % State variable P1 in kpa and T1 in K
 P1 = 100
 T1= 450
 
 %engine geometry parameter
 bore = 0.08
 stroke =0.12
 con_rod = 0.15
 cr = 9
 
 % swept volume and clearance volume 
 V_swept = (pi/4) * bore^2 * stroke % bore = D and stroke = L
 V_clearance = V_swept/(cr-1)
 V1=V_swept+V_clearance
 V2=V_clearance
 
 %state variable at point 2
 %p1v1^gamma = p2v2^gamma 
 P2=P1*cr^gamma
  T2 = P2*V2*T1/(P1*V1)
 
 %state variable at point 3
 V3=V2
 P3=P2*T3/T2 
 %state variable at point 4
 V4=V1
 P4= P3*(V3/V4)^gamma
 
 % thermal efficiency
 Percentage_Efficiency = (1-(cr^(gamma-1))^-1)*100 % Plotting 
 hold on
 
plot(V1,P1,'*','color','r')
plot(V2,P2,'*','color','r')
plot(V3,P3,'*','color','r')
plot(V4,P4,'*','color','r')
xlabel('Volume (m^3)')
ylabel('Pressure (KPa)')

Where:-

Clear all command used to clear text from the Command Window before running a command.

Close all command used to close all the plots window.

clc command is used to Clear Command Window

The input and state points are taken as discussed.

For plotting, we are using plot command and we get a PV diagram having four-point in it which Represents State points. after Executing the Program we get our required output as,

We cannot join them form line because isentropic processes do not have constant slope. If we try, we get the result as follows by using the following command, which is an Error.

plot([V1 V2 V3 V4 V1],[P1 P2 P3 P4 P1], ‘color’,’r’)

Step-2: Error Elimination.

To plot correctly, we need the following two things.

1. Thermodynamic relation during compression
2. Piston kinematics to determine How the volume change from TDC to BDC and how combustion volume changes as a function of the crank angle.

Program to get Piston kinematics Relation is as follows :

function [V] = piston_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)
 
angle = linspace(start_crank,end_crank,100)
 
%formula
term1 = 0.5*(cr-1);
term2= R+1-cosd(angle);
term3=(R^2-sind(angle).^2).^0.5;
 
V=(1+term1*(term2-term3)).*V_c;
end

Here, we are taking all the parameters same as above but we need a different formula for the calculation of the volume change inside the combustion chamber. The required relationship we can get by the following formula.

V/Vc = 1+ 0.5(cr-10)* [ R+1 -cos(angle) — (R² -sin²(angle)⁰.5]

To simplify the formula we divided formula into three sub formulas as follows

`term1 = 0.5*(cr-1)`;
`term2= R+1-cosd(theta)`;
`term3=(R²-sind²(theta))⁰.5`;

At last, we get our final formula as

`V=(1+term1*(term2-term3))*V_c`

Where:-

R= Length of connecting rod/Crankpin radius

Now, we are going to use this Program in our first program for getting PV Diagram correctly, But using total program to our first program is non-convenient therefore we are introducing the predefine command in Matlab known as function, which is a piece of code that we can read by calling its name.

Function [V] =piston_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank)

Note:- Name of the function should be the same as the name of file otherwise we get an error as follows.

So in our first program, we are going to use this function by calling its name.

Now, in our main program, we are using this function for the Pressure and volume changes in the compression and Expansion stroke, for the volume we can define the function as follows.

V_expansion = piston_kinematics(bore,stroke,con_rod,cr,180,360);

But for the pressure, we know that P*V^gamma = constant

Hence, we are defining this relation and using it as, constant_c = P1*V1^gamma; and P_compression = constant_c./V_compression.^gamma;

similarly, for Expansion, we are using the above relation as,

constant_c = P3*V3^gamma;
V_expansion = piston_kinematics(bore,stroke,con_rod,cr,180,360);
P_expansion = constant_c./V_expansion.^gamma;

Step-3: Calculating Thermal Efficiency:-

But in our Program, we need thermal efficiency of the cycle so to that

we are using the efficiency formula for the ideal otto cycle in percentage.

Thermal efficiency formula,

Percentage_Efficiency = (1-(cr^(gamma-1))^-1)*100

Final Matlab Program:

close all
clear all
clc
% Input
gamma= 1.4
 T3=2778

%engine geometry parameter
 bore = 0.1;
 stroke =0.1;
 con_rod = 0.15;
 cr = 12; 

 % State variables at point 1
 V_swept = (pi/4) * bore^2 * stroke % bore = D and stroke = L
 V_clearance = V_swept/(cr-1)
 V1=V_swept+V_clearance
 P1 = 100
 T1= 450

 %State point at point 2
 %p1v1^gamma = p2v2^gamma 
 V2=V_clearance;
 P2=P1*cr^gamma;
 T2 = P2*V2*T1/(P1*V1);
 constant_c = P1*V1^gamma;
 V_compression = piston_kinematics(bore,stroke,con_rod,cr,0,180)
 P_compression = constant_c./V_compression.^gamma;
 
 % State point at point 3
 V3=V2;
 P3=P2*T3/T2 ;
 constant_c = P3*V3^gamma;
 V_expansion = piston_kinematics(bore,stroke,con_rod,cr,180,360);
 P_expansion = constant_c./V_expansion.^gamma;

 % State point at point 4 
 V4=V1;
 P4= P3*(V3/V4)^gamma;

 % thermal efficiency
 Percentage_Efficiency = (1-(cr^(gamma-1))^-1)*100 % Plotting 

 % plotting
 hold on
 plot(V_compression,P_compression,'linewidth',1.5,'color', 'c')
 plot(V_expansion,P_expansion,'linewidth',1.5,'color', 'c' )
 plot([V2 V3],[P2 P3],'linewidth',1.5,'color', 'g')
 plot([V4 V1],[P4 P1],'linewidth',1.5,'color', 'g')
 plot(V1,P1,'*','color','r')
 plot(V2,P2,'*','color','r')
 plot(V3,P3,'*','color','r')
 plot(V4,P4,'*','color','r')
 xlabel('Volume (m^3)')
 ylabel('Pressure (KPa)')​