“AIR STANDARD CYCLE PROGRAM USING MATLAB AND PYTHON.”

13 min readJul 22, 2020

1. Objective:-

To write code in Matlab and Python 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.

The Behavior of the piston and crank cause changes in the process of the Otto cycle PV diagram is as shown below.

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.

Thermodynamic relation during compression
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% 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
 constant_c = P1*V1^gamma;
 V_compression = piston_kinematics(bore,stroke,con_rod,cr,180,0)
 P_compression = constant_c./V_compression.^gamma; T2 = P2*V2*T1/(P1*V1)
 
 %state variable 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 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(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)')

Where:-

Plot Command is used to get the final plot but we need to plot the two graphs in the same figure. Hence we are using subplot command to do it. subplot(2,1,1) indicates the first plot and subplot(2,1,2) indicates the second plot.

xlable() and ylabel() command used to define the labels in the plots. title() command is used to define the plot title.

4.2 Procedure and Explanation for Python:

Step-1: Importing Libraries.
Before Going to Do program we are using import command as follows.

import math
import matplotlib.pyplot as plt

Where:-
The “import math” and “import matplotlib.pyplot as plt” are used to use the mathematical constants and plotting the graph respectively.

Step-2: Basic Program.
First, we are going to do the basic program with the help of defined state variable points.

# Basic Program For Otto Cycleimport math
import matplotlib.pyplot as plt# Inputs 
p1 = 101325
t1 = 500
gamma = 1.4
t3 = 2300# Geometric parameters
bore = 0.1
stroke = 0.1
con_rod = 0.15
cr = 8# Volume Computation
v_s = (math.pi/4) * pow (bore,2)*stroke
v_c = v_s / (cr-1)
v1 = v_c + v_s# at state point 2
v2= v_c#p2v2^gamma = p1 v1^gamma
p2 = p1*pow(v1,gamma)/pow(v2,gamma)# p2v2/t2 = p1v1/t1 ; rhs = p2v2/t2 ; t2 = p2v2/rhsrhs = p1*v1/t1
t2 = p2*v2/rhs# at point 3
v3 =v2# p3v3/t3 = p2v2/t2 ; rhs = p2v2/t2 ; p3 = rhs*t3/v3rhs = p2*v2/t2
p3 = rhs* t3/v3# at point 4
v4 = v1
p4 = p3*pow(v3,gamma)/pow(v4,gamma)# p4v4/t4 = rhs
t4 = p4*v4/rhsplt.plot([v1,v2,v3,v4,v1],[p1,p2,p3,p4,p1])
plt.title('Basic PV Diagram for Otto Cycle')
plt.xlabel('Volume')
plt.ylabel('Pressure')
plt.savefig('Basic PV Diagram for Otto Cycle')
plt.show()

The Output of the Basic Program:

From the Above Plot, we can see that the compression and Expansion processes are shown linear, which is practically not possible. Hence to make the proper plot we need the Engine kinematics equation as discussed below.

Step-3: Defining Engine Kinematics function.

To plot correctly, we need the following two things.

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

Here, we are taking all the parameters the same as discussed in step 1 and 2, 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. The Formula for the Engine kinematics is as follows.

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

To use this formula in the program we simplified it. Hence, To simplify the formula we divided formula into three sub formulas (terms) as follows.

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

At last, After combining the terms to get the final formula

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

After that our Engine kinematics formula is used as a function in the program is as shown below.

def engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank):	# Geometric Parameters
	a = stroke/2
	R = con_rod/a	# Volume Parameters
	V_s = math.pi*(1/4)*pow(bore,2)*stroke
	V_c = V_s/(cr-1)	# Position of Starting Crank and Ending Crank
	sc = math.radians(start_crank)
	ec = math.radians(end_crank)
	
	num_values = 50	dtheta = (ec-sc)/(num_values-1)	V = []	for i in range(0,num_values):
		theta = sc + i*dtheta
		term1 = 0.5*(cr-1)
		term2 = R + 1 - math.cos(theta)
		term3 = pow(R,2) - pow(math.sin(theta),2)
		term3 = pow(term3,0.5)		V.append((1+ term1 * (term2 - term3))*V_c)
	return V

Where:-
i. We need the formula of engine kinematics to create the volume Arrays concerning the compression and Expansion Processes. After that, we can use it in our Otto cycle Program. Also, we need to calculate the pressure that corresponds to the volume. Therefore, We defined the Engine kinematics as the function.

ii. Since we have a function so we need not pre-define the value of the parameters related to it such as values of bore, stroke, connecting rod length, etc.

iii. Now defining the sc and ec (i.e starting crank angle and ending crank angle). Because we want to get several values between the starting crank and ending crank position and each of those values we need to get the volume. hence using ec and sc and assigning it to the dtheta which is known as a crank angle spacing. Also defined the number of values is equal to 50 as follows.

“# Position of Starting Crank and Ending Crank
sc = math.radians(start_crank)
ec = math.radians(end_crank)

num_values = 50
dtheta = (ec-sc)/(num_values-1)”

iv. Now we can form a crank angle array. to do that we are using for-loop. For loop use to repeat the set of actions multiple times. The empty volume array indicates that when for loop executed each time the value calculated in the for-loop gets store in the empty array.

iv. But to generate the multiple values we need one more command called as append.”.append” is used to continuously append (increase) the value and store it in the empty array.

v. using “return V” to get multiple output values of volume.

Step-4: Defining the Compression and Expansion processes.

To define the compression and Expansion processes in the program we took help to the engine kinematic function as shown below.

# for Compression process at State point 2
V_comp  = engine_kinematics(bore,stroke,con_rod,cr,180,0)const = p1*pow(v1,gamma)P_comp =[]for v in V_comp:
	P_comp.append(const/pow(v,gamma))
# for Expansion process at State point 4
V_exp  = engine_kinematics(bore,stroke,con_rod,cr,0,180)const = p3*pow(v3,gamma)P_exp =[]for v in V_exp:
	P_exp.append(const/pow(v,gamma))

Where:-
For loop is used to get multiple values of the volume and pressure and store it in the empty array using Append command.

For compression, the taken starting crank angle is at 180 degrees and the ending crank angle is 0 degree (i.e piston moves from BDC to TDC). Similarly, For Expansion, the taken starting crank angle is at 0 degrees and the ending crank angle is 180 degrees (i.e piston moves from TDC to BDC).

Note:- Here We took the crank angle for reference to the compression and expansion process not a movement of the piston. Also, Suction and exhaust are neglected.

Step-5: Calculating Thermal Efficiency.

we are using the efficiency formula for the ideal otto cycle in percentage we can determine the efficiency of the otto cycle. The Thermal efficiency formula in percentage is as following.,

Percentage_Efficiency = `1–1/(cr^(gamma-1))*100`

# Calculating Thermal Efficiency
n = pow(pow(cr,(gamma-1)),-1)
Percentage_Efficiency = (1-n)*100
print(Percentage_Efficiency)

Step-6: Creating Plot.

After defining all the terms we need the output so to do that we are using the plot command and plotted all the processes.

After properly arranging the program we get out the desired program and output which is shown below.

Final Program in Python:

# Final Program For Otto Cycleimport math
import matplotlib.pyplot as plt def engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank):	# Geometric Parameters
	a = stroke/2
	R = con_rod/a	# Volume Parameters
	V_s = math.pi*(1/4)*pow(bore,2)*stroke
	V_c = V_s/(cr-1)	# Position of Starting Crank and Ending Crank
	sc = math.radians(start_crank)
	ec = math.radians(end_crank)
	
	num_values = 50
	dtheta = (ec-sc)/(num_values-1)	V = []	for i in range(0,num_values):
		theta = sc + i*dtheta
		term1 = 0.5*(cr-1)
		term2 = R + 1 - math.cos(theta)
		term3 = pow(R,2) - pow(math.sin(theta),2)
		term3 = pow(term3,0.5)		V.append((1+ term1 * (term2 - term3))*V_c)
	return V#----------------------------------------------------------
# Inputs 
p1 = 101325
t1 = 500
gamma = 1.4
t3 = 2300# Geometric parameters
bore = 0.1
stroke = 0.1
con_rod = 0.15
cr = 8# Volume Computation
v_s = (math.pi/4) * pow (bore,2)*stroke
v_c = v_s / (cr-1)
v1 = v_c + v_s#-----------------------------------------------------------
# at state point 2 
v2= v_c#p2v2^gamma = p1 v1^gamma
p2 = p1*pow(v1,gamma)/pow(v2,gamma)# p2v2/t2 = p1v1/t1 ; rhs = p2v2/t2 ; t2 = p2v2/rhs
rhs = p1*v1/t1
t2 = p2*v2/rhsV_comp  = engine_kinematics(bore,stroke,con_rod,cr,180,0)
const = p1*pow(v1,gamma)P_comp =[]for v in V_comp:
	P_comp.append(const/pow(v,gamma))#-----------------------------------------------------------
# at state point 3
v3 =v2# p3v3/t3 = p2v2/t2 ; rhs = p2v2/t2 ; p3 = rhs*t3/v3
rhs = p2*v2/t2
p3 = rhs* t3/v3#-----------------------------------------------------------
# at state point 4
v4 = v1
p4 = p3*pow(v3,gamma)/pow(v4,gamma)# p4v4/t4 = rhs
t4 = p4*v4/rhsV_exp  = engine_kinematics(bore,stroke,con_rod,cr,0,180)
const = p3*pow(v3,gamma)
P_exp =[]for v in V_exp:
	P_exp.append(const/pow(v,gamma))#-----------------------------------------------------------
# Calculating Thermal Efficiency
n = pow(pow(cr,(gamma-1)),-1)
Percentage_Efficiency = (1-n)*100
print(Percentage_Efficiency)#-----------------------------------------------------------
# Potiing
plt.plot([v2,v3],[p2,p3])
plt.plot(V_comp,P_comp)
plt.plot(V_exp,P_exp)
plt.plot([v4,v1],[p4,p1])
plt.title('Final PV Diagram for Otto Cycle')
plt.xlabel('Volume')
plt.ylabel('Pressure')
plt.savefig('Final PV Diagram for Otto Cycle')
plt.show()