Otto cycle simulation using python

Based on the concepts of engine kinematics. The code for simulation of otto cycle is given below. With the help of this we can simulate other cycles too.

First we will need to have some inputs.

Then defining the function we can create a pv plot for otto cylce and get the thermal efficiency.

import matplotlib.pyplot as plt
import math


def engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank):                  #defining function

	#engine kinematics

	a = stroke/2                #Area
	R = con_rod/a

	#volume calculations
	v_s = (math.pi/4)*pow(bore,2)*stroke
	v_c = v_s/(cr-1)
	
	sc = math.radians(start_crank)    #starting crank
	ec = math.radians(end_crank)      #ending crank

	num_values = 200                   #numer of values

	dtheta = (ec - sc)/(num_values-2)

	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 = 270
gamma = 1.4
t3 = 1800

#geometric parameters
bore = 0.1
stroke = 0.1
con_rod = 0.15
cr = 12

#volume computation
v_s = (math.pi/4)*pow(bore,2)*stroke
v_c = v_s/(cr-1)
v1 = v_c + v_s

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


rhs = p1*v1/t1                           #p2v2/t2 = p1v1/t1| rhs = p1v1/t1 | p2v2/t2 = rhs | t2 = p2v2/rhs
t2 = p2*v2/rhs
V_compression = engine_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))            #P_compression * V_comression ^ gamma = constant

#state point 3
v3 = v2


rhs = p2*v2/t2                      #p3v3/t3 = p2v2/t2 | rhs = p2v2/t2 | p3 = rhs*t3/v3
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
p4 = p3*pow(v3,gamma)/pow(v4,gamma)       #p4v4^gamma = p3v3^gamma

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


#thermal efficiency of otto cycle
efficiency=1-(1/pow(cr,gamma-1))          #formula for thermal efficiency
print('THe thermal efficiency of otto cycle is',efficiency)

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')                         #x axis label
plt.ylabel('pressure')                       #y axis label
plt.title('OTTO CYCLE SIMULATION')
plt.show()