Visualization of P-V diagram of Air standard Otto cycle using Python & MATLAB

OBJECTIVE:

To obtain a p-v diagram of air standard otto cycle and find it's Thermal efficiency using Python & MATLAB program.

THEORY:

Otto cycle

Otto cycle is a constant volume thermodynamic cycle which is used in Spark Ignition engines to produce power. It consists of 4 strokes. They are

1. Sucktion stroke (1-2:Isentropic compression)
2. Compression stroke (2-3:Constant volume heat addition)
3. Expansion stroke (3-4:Isentropic Expansion)
4. Exhaust stroke (4-1:Constant volume heat rejection)

Engine Nomenclature of otto cycle is given below

1. Bore
2. Stroke
3. Connecting rod
4. Compression ratio

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~MATLAB~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

MAIN PROGRAM for Otto Cycle

clear all
close all
clc

%PV DIAGRAM FOR OTTO CYCLE

%INPUT PARAMTERS ASSUMED
gamma = 1.4;
p1 = 101325;
t1 = 500;
t3 = 2300;

%ENGINE PARAMETERS TO FIND v1 & v2
bore = 0.1;
stroke = 0.1;
connecting_rod = 0.15;
compression_ratio = 12;

%Formula for SWEPT VOLUME & CLEARANCE VOLUME
swept_volume = pi/4*bore^2*stroke;
clearance_volume = swept_volume/(compression_ratio-1);

%Determination of v1 & v2
v1 = swept_volume + clearance_volume;
v2 = clearance_volume;


%STATE VAIRABLE AT POINT 1
    %v1 ,p1 ,t1 are already assumed & calculated
    

%STATE VARIABLE AT POINT 2
    %to find p2 and t2
    %wkt p2v2^gamma = p1v1^gamma (i.e) p2 = p1(v1/v2)^gamma
    %cr = v1/v2
p2 = p1*compression_ratio^gamma;
    
    %wkt p1*v1/t1 = p2*v2/t2
t2 = p2*v2*t1/(p1*v1);
    
    %Using engine kinematics equation gives value for different crank angles at point 2-3
constant_c = p1*v1^gamma;
V_compression = engine_kinematics(bore,stroke,connecting_rod,compression_ratio,180,0);
P_compression = constant_c./V_compression.^gamma;


%STATE VARIABLE AT POINT 3
    %to find p3
    %t3 is assigned
    %due to constant volume heat addition, v3 = v2
v3 = v2;
    
    %wkt p3*v3/t3 = p2*v2/t2
    %v3 = v2
p3 = (p2*t3)/t2;
    
    %Using engine kinematics equation gives value for different crank angles at point 3-4
constant_c = p3*v3^gamma;
V_expansion = engine_kinematics(bore,stroke,connecting_rod,compression_ratio,0,180);
P_expansion = constant_c./V_expansion.^gamma;
    

%STATE VARIABLE AT POINT 4
    %to find p4
    %due to constant volume heat rejection, v4 = v1
v4 = v1;

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


%PLOTTING
hold on
title('Otto cycle PV diagram')
xlabel('Volume(v)')
ylabel('Pressure(p)')

plot(v1,p1,'marker','*','color','k')
plot(V_compression,P_compression,'color','g','linewidth',1)
plot(v2,p2,'marker','*','color','k')
plot(v3,p3,'marker','*','color','k')
plot([v2 v3],[p2 p3],'color','b','linewidth',1)
plot(V_expansion,P_expansion,'color','r','linewidth',1)
plot(v4,p4,'marker','*','color','k')
plot([v4 v1],[p4 p1],'color','k','linewidth',1)


%TO CALCULATE THERMAL EFFICIENCY
Thermal_efficiency = (1 - (1/(compression_ratio)^(gamma - 1)))*100;

Explanation:

• The program begins with usual clear all, close all, clc to avoid manipulation of the previous values.
• The input variables are defined such as gamma= 1.4 and for state point 1 we assign the data for p1 = 101325, t1 = 500 and t3 = 2300 for state point 3.
• The engine parameters for finding v1 & v2 are given as bore = 0.1, stroke = 0.1, connecting rod = 0.15, compression ratio = 12.
• Formula for Swept volume & Clearance volume is now defined from which we can calculate v1 & v2.
• Now we have to calculate all the state vairables at state points 1 to 4 which is explained from line 28 to line 70 with the respective equations and conditions known.
• v_compression & p_compression is used to plot the curve (1-2) by the function defined in the first program ((i.e) engine_kinematics.m) given in the line 43,44 & 45 along with constant_c known standard value`
• Similarly to plot the curve (3-4), v_expansion & p_expansion is defined from the same program with constant_c given in line 59,60 & 61.
• Here for the compression the start & end angle is set at 180 & 0 degree respectively, since the compression should end with 0 degree whereas for the expansion it is vice versa ((i.e)start - 0,end-180).
• Plot command is now used to plot the 4 points which is executed four time for the 4 coordinates. Color, maker, linewidth are used within plot command to make the graph more appealable.
• To connect the points 1 & 2 and 3 & 4, v_compression & p_compression and v_expansion & p_expansion is used.
• Since it is a constant volume process, the points 2 & 3 and 4 & 1 are connected by including appropriate v and p values in pairs.

Function program of engine kinematics

%Determination of the function "engine kinematics" with the equation V

function[V] = engine_kinematics(bore,stroke,connecting_rod,compression_ratio,start_crank,end_crank)

%Formula for a(crank pin radius) & R
a = stroke/2;
R = connecting_rod/a;

swept_volume = [pi/4]*bore^2*stroke;
clearance_volume = swept_volume/(compression_ratio-1);
thetha = linspace(start_crank,end_crank,100);

%Breaking down the equation into 3 terms
term1 = 0.5*(compression_ratio-1);
term2 = R + 1 - cosd(thetha);
term3 = (R^2 - sind(thetha).^2).^0.5;

%Whole Equation
V = (1 + term1*(term2 - term3)).*clearance_volume

end

Explanation:

• This function is determined for representing the isentropic curve
• The function engine kinematics is assigned for input values such as bore,stroke,connecting rod,compression ratio,start crank and end crank.
• The inputs for this program will be a(crank pin radius) and R(ratio of connecting rod & a)
• Swept volume and Clearance volume is given by the formula.
• thetha is defined from the interval of 100 with start and end crank angle being inputs.
• The above equation is broken down into three terms and defined.
• The whole equation V is now determined.

Output:

As we can see the input for the function "engine kinematics" gives the above output of 100 values.

RESULT:

Output of P-V diagram

Calculated values

The P-V diagram for the Air standard Otto cycle has been plotted and the Thermal efficiency has been found to be 62.98%.

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~PYTHON~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

MAIN PROGRAM

import math
import matplotlib.pyplot as plt

#PV DIAGRAM FOR OTTO CYCLE
#Determination of the function \"engine_kinematics\"

def engine_kinematics(bore,stroke,connecting_rod,compression_ratio,start_crank,end_crank):

	#Formula for a(crank pin radius) & R
	a = stroke/2
	R = connecting_rod/a

	#Formula for SWEPT VOLUME & CLEARANCE VOLUME
	swept_volume = math.pi*(1/4)*pow(bore,2)*stroke
	clearance_volume = swept_volume/(compression_ratio-1)

	sc = math.radians(start_crank)
	ec = math.radians(end_crank)

	number_of_values = 50
	dt = (ec - sc)/(number_of_values - 1)
	V = []

	for i in range(0,number_of_values):
		theta = sc + i*dt
		
		#Breaking down the equation into 3 terms
		term1 = 0.5*(compression_ratio - 1)
		term2 = R + 1 - math.cos(theta)
		term3 = pow((pow(R,2) - pow(math.sin(theta),2)),0.5)

		#Whole Equation
		V.append((1 + term1*(term2 - term3))*clearance_volume)
	return V


#INPUT PARAMTERS ASSUMED
gamma = 1.4
p1 = 101325
t1 = 500
t3 = 2300

#ENGINE PARAMETERS TO FIND v1 & v2
bore = 0.1
stroke = 0.1
connecting_rod = 0.15
compression_ratio = 12

#Formula for SWEPT VOLUME & CLEARANCE VOLUME
swept_volume = (math.pi/4)*pow(bore,2)*stroke
clearance_volume = swept_volume/(compression_ratio - 1)

#Determination of v1 & v2
v1 = swept_volume + clearance_volume
v2 = clearance_volume


#STATE VAIRABLE AT POINT 1
#v1 ,p1 ,t1 are already assumed & calculated
    

#STATE VARIABLE AT POINT 2
#to find p2 and t2
#wkt p2v2^gamma = p1v1^gamma (i.e) p2 = p1(v1/v2)^gamma
p2 = p1*pow(v1,gamma)/pow(v2,gamma)
    
#wkt p1*v1/t1 = p2*v2/t2 = rhs
#rhs = right hand side
rhs = p1*v1/t1
t2 = (p2*v2)/rhs

#For start of expansion stroke
#Using engine kinematics equation gives value for different crank
#angles at point 2-3
constant1 = p1*pow(v1,gamma)
V_compression = engine_kinematics(bore,stroke,connecting_rod,compression_ratio,180,0)
P_compression = []

for v in V_compression:
	P_compression.append(constant1/pow(v,gamma))


#STATE VARIABLE AT POINT 3
#to find p3
#t3 is assigned
#due to constant volume heat addition, v3 = v2
v3 = v2
rhs = p2*v2/t2
    
#wkt p3*v3/t3 = p2*v2/t2
#v3 = v2
p3 = rhs*t3/v3
    
#Using engine kinematics equation gives value for different crank angles at point 3-4
constant2 = p3*pow(v3,gamma)
V_expansion = engine_kinematics(bore,stroke,connecting_rod,compression_ratio,0,180)
P_expansion = []

for v in V_expansion:
	P_expansion.append(constant2/pow(v,gamma))


#STATE VARIABLE AT POINT 4
#to find p4
#due to constant volume heat rejection, v4 = v1
v4 = v1

#wkt p4*v4^gamma = p3*v3^gamma
p4 = p3*pow(v3,gamma)/pow(v4,gamma)
t4 = p4*v4/rhs


#PLOTTING
plt.figure()
plt.plot([v2,v3],[p2,p3])
plt.title(\'Otto cycle PV diagram\')
plt.xlabel(\'Volume(v)\')
plt.ylabel(\'Pressure(p)\')
plt.plot(V_compression,P_compression)
plt.plot(V_expansion,P_expansion)
plt.plot([v4,v1],[p4,p1])
plt.show()


#TO CALCULATE THERMAL EFFICIENCY
Thermal_efficiency = (1 - pow(1/(compression_ratio),(gamma - 1)))*100
print(Thermal_efficiency)

Explanation:

• The program begins with importing math & matplotlib.pyplot to calculate math values & plot graphs respectively.
• The following formula is determined for representing the isentropic curve for equation V.
  
• The function engine kinematics is assigned for input values such as bore,stroke,connecting rod,compression ratio,start crank and end crank.
• The inputs for this program will be a(crank pin radius) and R(ratio of connecting rod & a)
• Swept volume and Clearance volume is given by the formula.
• theta is defined from the interval of 50 with start and end crank angle being inputs.
• The input variables are defined such as gamma= 1.4 and for state point 1 we assign the data for p1 = 101325, t1 = 500 and t3 = 2300 for state point 3.
• The engine parameters for finding v1 & v2 are given as bore = 0.1, stroke = 0.1, connecting rod = 0.15, compression ratio = 12.
• Formula for Swept volume & Clearance volume is now defined from which we can calculate v1 & v2.
• Now we have to calculate all the state vairables at state points 1 to 4 which is explained from line 28 to line 70 with the respective equations and conditions known.
• v_compression & p_compression is used to plot the curve (1-2) by the function defined ((i.e) engine_kinematics.m) given along with constant1 known standard value.
• Similarly to plot the curve (3-4), v_expansion & p_expansion is defined from the same program with constant2 given in line 59,60 & 61.
• Here for the compression the start & end angle is set at 180 & 0 degree respectively, since the compression should end with 0 degree whereas for the expansion it is vice versa ((i.e)start - 0,end-180).
• Plot command plt is now used to plot the 4 points which is executed four time for the 4 coordinates.
• To connect the points 1 & 2 and 3 & 4, v_compression & p_compression and v_expansion & p_expansion is used.
• Since it is a constant volume process, the points 2 & 3 and 4 & 1 are connected by including appropriate v and p values in pairs.

OUTPUT:

Plot of PV diagram

Thermal efficiency of engine

RESULT:

The P-V diagram for the Air standard Otto cycle has been plotted and the Thermal efficiency has been found to be 62.98%.

CONCLUSION:

We can conclude that using both Python & MATLAB yielded same results with Thermal efficiency being identical.