Otto cycle using python

The program written is for simulating an Air standard cycle (Otto cycle). The Otto cycle consists of following four processes.

• 1-2 Isentropic compression
• 2-3 Constant volume heat addition
• 3-4 Isentropic expansion
• 4-1 Constant volume heat rejection

This requires inputs like V1, T1, Compression ratio(cr),T3 from which P1,V2,P2,T2,P3 and V3 can be calculated using the formula (P1*V1/T1)=(P2*V2/T2) and (P*V^gamma)=constant. This also requires parameters like Bore, Stroke, Connecting rod length, Compression ratio.

The computer program for plotting PV diagram and thermal efficiency of the Engine is as follows.

"""
  otto cycle simulator
"""

import math
import matplotlib.pyplot as plt

def eng_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank):
    
    #geometric parameters
	a=stroke/2
	R=con_rod/a

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

	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=.5*(cr-1)
		term2=R+1-math.cos(theta)
		term3=pow(pow(R,2)-pow(math.sin(theta),2),.5)

		v.append((1+term1*(term2-term3))*v_c)
	
	return v


#inputs
gamma=1.4
p1=101.325
t1=500
t3=2300

# eng geometric parameters
bore=0.1
stroke=0.1
con_rod=.15
cr=12

#calculating the swept vol and clearence vol
v_s=(math.pi/4)*pow(bore,2)*stroke
v_c=v_s/(cr-1);

v1=v_s+v_c

#state point 2
v2=v_c

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

#p2*v2/t2=p1*v1/t1 | rhs=p1*v1/t1 | lhs=p2*v2/t2
t2=t1*p2*v2/(p1*v1)
v_compression=eng_kinematics(bore,stroke,con_rod,cr,180,0)
constant=p1*pow(v1,gamma)

p_compression=[]

for v in v_compression:
	p_compression.append(constant/pow(v,gamma))

#state point 3
v3=v2
p3=p2*t3/t2

v_expansion=eng_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
p4=p3*pow((v3/v4),gamma)
t4=t3/pow(cr,gamma-1);


thermal_efficiency=1-(1/pow(cr,gamma-1))
print('Thermal efficiency is=',thermal_efficiency*100)


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('v')
plt.ylabel('p')
plt.show()

plot is as shown below

Output is