Solving the Otto Cycle to obtain the thermal efficiency of an engine and generate its PV Diagram at different crank angles using Python

Given below is a program code that can solve the otto cycle and make plots for the thermal efficiency of an engine and generate a PV diagram of it taking in input values of :

1. Bore Diameter
2. Stroke Length
3. Connecting Rod length
4. The Compression Ratio
5. Start Crank Angle
6. End Crank Angle
7. The specific heat ratio, gamma
8. The temperatures at intake and at combustion (Power Stroke)

Thermal efficiency for an otto cycle engine is given by the formula:

where r is the Compression ratio and gamma is the specific heat ratio.

# Otto Cycle Simulator

import math
import matplotlib.pyplot as plt

def engine_kinematics(bore, stroke, con_rod, cr, start_crank, end_crank):
	\"\"\"
	Engine Kinematics

	\"\"\"

	#crank pin radius
	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)

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

	num_values = 100
	V=[]

	#crank angle spacing
	dtheta=(ec-sc)/(num_values-1)

	for i in range(0,num_values-1):	
		theta=sc + i*dtheta
		

		#V/Vc = 1 + 0.5(rc-1)[R+1-cos(theta)-(R^2 - sin^2(theta))^0.5]
		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 parameter (meters)
bore=0.1
stroke=0.1
con_rod=0.15
cr=2

#Thermal efficiency, n_th = 1- 1/(cr^(gamma-1))
term7=pow(cr,gamma-1)
n_th = 1 -1/term7 
print(n_th)

#Volume Computation, swept vol, clearence vol
v_s= (math.pi/4)*pow(bore,2)*stroke
v_c=v_s/(cr-1)

#State point 2
v1=v_c + v_s
 		
#State point 2
v2=v_c

# p2v2^gamma = p1v1^gamma
p2=p1*pow(v1,gamma)/pow(v2,gamma)

#p2v2/t2 = p1v1/t1 | rhs = p1v1/t1 | p2v2=rhs | t2 = p2v2/rhs
rhs=p1*v1/t1
t2=p2*v2/rhs

V_compresion=engine_kinematics(bore, stroke, con_rod, cr, 180,0)
constant=p1*pow(v1,gamma)

P_compression = []

for v in V_compresion:
	P_compression.append(constant/pow(v,gamma))


#state point 3
v3=v2

#p3v3/t3 = p2v2/t2 | rhs = p2v2/t2 | p3=rhs*t3/v3
rhs=p2*v2/t2
p3=rhs*t3/v3

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

P_expansion = []

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

#state point 4
v4=v1

#p4v4^gamma = p3v3^gamma
p4=p3*pow(v3,gamma)/pow(v4,gamma)

#p4v4/t4 = rhs
t4= p4*v4/rhs

f1 = plt.figure(\'PV Diagram\')
plt.plot([v2,v3],[p2, p3])
plt.plot(V_compresion,P_compression)
plt.plot(V_expansion,P_expansion)
plt.plot([v4,v1],[p4, p1])
plt.title(\'PV Diagram\')
plt.xlabel(\'Volume (m^3)\')
plt.ylabel(\'Pressure (Pa)\')
plt.show()

For, (in meters)

1. Bore Diameter = 0.1
2. Stroke Length = 0.1
3. Connecting Rod length = 0.15
4. The Compression Ratio = 2
5. Start Crank Angle = 0 Degrees
6. End Crank Angle = 180 Degrees
7. The specific heat ratio, gamma = 1.4
8. The temperatures at intake (t1=500) and at combustion (Power Stroke) (t3=2300)

The resulting PV diagram:

The thermal efficiency of the engine = 0.24214