P-V Diagram and Thermal Efficiency for an OTTO CYCLE in Python

Aim:

To write and execute the program for creating the PV diagram for Otto Cycle and to determine the thermal efficiency of the cycle

Description:

Otto Cycle is the thermodynamic cycle which describes the fuctioning of Spark Ignition Engines,commonly found in automobile engines.This cycle consists of four processes which depicts four state points for the construction of the PV diagram.

State Point 1 to 2 --- Isentropic Compression Process

State Point 2 to 3 --- Constant Volume Heat Addition Process

State Point 3 to 4 --- Isentropic Expansion Process

State Point 4 to 1 --- Constant Volume Heat Rejection Process

Once the state points are calculated, the Thermal Efficiency of the Otto Cycle can be calculated from,

Thermal Efficiency = `1-((t4-t1)/(t3-t2))`

Coding:

#otto cycl e simulator

import math
import matplotlib.pyplot as plt


def engine_kinematics(bore,stroke,con_rod,rc,initial_crankangle,final_crankangle):
	
	r = stroke/2
	R = con_rod/r

	V_s = math.pi/4*pow(bore,2)*stroke
	V_c = V_s/(rc-1)

	sc = math.radians(initial_crankangle)
	ec = math.radians(final_crankangle)

	num = 60

	dtheta = (ec-sc)/(num-1)

	V=[]

	for i in range(0,num):
		theta = sc + i*dtheta
		term1 = 0.5*(rc-1)
		term2 = R+1-math.cos(theta)
		term3 = pow(R,2)-pow(math.sin(theta),2)
		term4 = pow(term3,0.5)
			   
		V.append((1+term1*(term2-term4))*V_c)
	return V

#Inputs
p1 =100000
t1 =400
t3 = 2300
gamma = 1.4

#Engine Geometric Parameter
bore = 0.8
stroke = 0.8
con_rod = 1.0
rc = 15

#Calculation Of Swept Volume and Clearance Volume
V_s = math.pi/4*pow(bore,2)*stroke
V_c = V_s/(rc-1)

#State points in Stage 1
v1 = V_c + V_s
const = p1*pow(v1,gamma)

#State Points in Stage 2
v2 = V_c
p2 = const/pow(v2,gamma)
t2 = ((p2*v2)/(p1*v1))*t1 


V_compression=engine_kinematics(bore,stroke,con_rod,rc,180,0)

P_compression =[]

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


#State Points in Stage 3
v3 = v2
p3 = (t3/t2)*p2

V_expansion=engine_kinematics(bore,stroke,con_rod,rc,0,180)
const = p3*pow(v3,gamma)

P_expansion =[]

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

#State Points in Stage 4
v4 = v1
p4 = const/(pow(v4,gamma))
t4 = ((p4*v4)/(p3*v3))*t3


#Plotting
plt.plot([v2,v3],[p2,p3])
plt.plot(V_compression,P_compression)
plt.plot(V_expansion,P_expansion)
plt.plot([v4,v1],[p4,p1])
plt.xlabel(\'VOLUME in m^3\')
plt.ylabel(\'PRESSURE in N/m^2\')
plt.show()

#Calculating Thermal Efficiency
thermal_eff = (1-((t4-t1)/(t3-t2)))*100
print(\'The Thermal Efficiency of the Otto Cycle is\', thermal_eff,\'%\')

Coding Explaination:

import math
import matplotlib.pyplot as plt


def engine_kinematics(bore,stroke,con_rod,rc,initial_crankangle,final_crankangle):
	
	r = stroke/2
	R = con_rod/r

	V_s = math.pi/4*pow(bore,2)*stroke
	V_c = V_s/(rc-1)

	sc = math.radians(initial_crankangle)
	ec = math.radians(final_crankangle)

	num = 60

	dtheta = (ec-sc)/(num-1)

	V=[]

	for i in range(0,num):
		theta = sc + i*dtheta
		term1 = 0.5*(rc-1)
		term2 = R+1-math.cos(theta)
		term3 = pow(R,2)-pow(math.sin(theta),2)
		term4 = pow(term3,0.5)
			   
		V.append((1+term1*(term2-term4))*V_c)
	return V

In MATLAB, we should have to write the program code for Function Definition in the seperate tab.But in Python, we can define the function within the main program code itself.Here the function for variation Volume of piston with respect to the crank angle during expansion and expression is defined.The Formula is,

In Compression, the initial crank angle and the final crank angle is setted between 180 to 0 degrees and in the case of Expansion, it is setted up 0 to 180 degrees.Then the crank angle is divided into so many values and the volume of piston at those values is calculated and the empty Volume array is created in order to dump that claculated values in that array using \'append\' operator.And \'return\' command is used to return that array ,where the function is called.

#Inputs
p1 =100000
t1 =400
t3 = 2300
gamma = 1.4

#Engine Geometric Parameter
bore = 0.8
stroke = 0.8
con_rod = 1.0
rc = 15

The inputs of the Air Standard Cycle such as initial pressure(p1),initial temperature(t1),temperature after combustion(t3) are defined.Then the dimension of the engine,that is,Engine Geometric Parameters sucha as bore,stroke,connecting rod length,compression ration are defined.

#Calculation Of Swept Volume and Clearance Volume
V_s = math.pi/4*pow(bore,2)*stroke
V_c = V_s/(rc-1)

The Swept Volume and the Clearance Volume are calculated from the inputs.Swept Volume is the volume through which piston moves and the Clearance Volume is the volume which is not swept by the piston.

#State points in Stage 1
v1 = V_c + V_s
const = p1*pow(v1,gamma)

#State Points in Stage 2
v2 = V_c
p2 = const/pow(v2,gamma)
t2 = ((p2*v2)/(p1*v1))*t1 


V_compression=engine_kinematics(bore,stroke,con_rod,rc,180,0)

P_compression =[]

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

The Pressure and Volume for the state point 1 is calculated and for the state point 2, and inbetween this the variation of Volume in compression is calculated by calling the Function named \'engine_kinematics\' and returned Volume value from the function is stored in the array \'V_compression\' and consequently Pressure for that compression is calculated from the adiabatic relation,P*(V)^gamma

#State Points in Stage 3
v3 = v2
p3 = (t3/t2)*p2

V_expansion=engine_kinematics(bore,stroke,con_rod,rc,0,180)
const = p3*pow(v3,gamma)

P_expansion =[]

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

#State Points in Stage 4
v4 = v1
p4 = const/(pow(v4,gamma))
t4 = ((p4*v4)/(p3*v3))*t3

Next, the pressure and the Volume for state point 3 and for state point 4.In the similar way, that we did it in the compression, the variation of volume and pressure is calculated for expansion.

#Plotting
plt.plot([v2,v3],[p2,p3])
plt.plot(V_compression,P_compression)
plt.plot(V_expansion,P_expansion)
plt.plot([v4,v1],[p4,p1])
plt.xlabel(\'VOLUME in m^3\')
plt.ylabel(\'PRESSURE in N/m^2\')
plt.show()

In this section, PV diagram is plotted from the determined values of Pressure,Volume for all the state points using \'plt.plot()\' command,labels are specifien and \'plt.show()\' command is used to show the figure.

#Calculating Thermal Efficiency
thermal_eff = (1-((t4-t1)/(t3-t2)))*100
print(\'The Thermal Efficiency of the Otto Cycle is\', thermal_eff,\'%\')

Then, Finally the Thermal Efficiency of the Otto Cycle is determined by the temperature limit values and represented in percentage(%).

Output: