Solve an otto cycle and make plots for it

Objective : write code that can solve an otto cycle and make plots for it.

Problem Statement : Here are the requirements

1. Your code should create a PV diagram
2. You should output the thermal efficiency of the engine.

Solution Procedure :

It is thermodynamic cycle which is used in internal combustion engine to produce work.

There are different type of processes in Otto cycle :-

1. Adiabatic compression (1-2)- air and fuel mixture is compressed isentropically inside the cylinder.
2. Constant volume heat addition (2-3)- combustion take place inside the cylinder due to external source.
3. Adiabatic expansion (3-4) - high pressure and high temperature gases expand isentropically inside the cylinder.
4. Constant volume heat rejection (4-1):- gases come out from the cylinder and cycle is completed.

we have to give required inputs and by using basic thermodynamic formulas type in the code and save the required output. and to make a adiabatic curve by using the calling the function .

""" 
otto cycle simulator
"""
import math
import matplotlib.pyplot as plt
import numpy as np

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

	#geometric parameters
	#crank pin radius
	a=stroke/2
	R=con_rod/a

	#volume comutation
	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_value=50

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

	V=[]

	for i in range(0,num_value):
		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=500
gamma=1.4
t3=2300

bore=0.1
stroke=0.1
con_rod=0.15
cr=13


#volume comutation
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

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

#p2v2/t2=p1v1/t1 |  rhs=p1v1/t1 | p2v2/t2=rhs |

rhs=p1*v1/t1

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))

# state point 3

v3=v2

#p3v3/t3=p2v2/t2 | rhs=p2v2/t2 | p3v3/t3=rhs |

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=p3v3/t3 | rhs=p3v3/t3  | p4v4/t4=rhs |

rhs=p3*v3/t3 

t4=p4*v4/rhs
plt.figure(1)
plt.plot(V_compression,P_compression,label='Adiabatic compression')
plt.plot([v2,v3],[p2,p3],label='Constant volume heat addition')
plt.plot(V_expansion,P_expansion,label='Adiabatic expansion')
plt.plot([v4,v1],[p4,p1],label='Constant volume heat rejection')
plt.xlabel('volume')
plt.ylabel('pressure')
plt.legend()
plt.savefig('otto cycle pv diagram')
plt.show()


#thermal efficiency

eff=1-(1/pow(cr,(gamma-1)))
eff=str(eff*100)
print('thermal efficiency =' + eff + '%')

results :