Air cycle

# otto cycle
import math
import matplotlib.pyplot as plt

def air_cycle(bore,stroke,con_rod,cr,start_crank,end_crank):

	


	#geometric parameters 

	
	a=stroke/2
	R=con_rod/a

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

	detheta = (ec-sc)/(num-1)
	V=[]
	for i in range(0,num):
	  theta = sc + i*detheta
	  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 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

#state point 2
v2 = v_c
p2=p1*pow(v1,gamma)/pow(v2,gamma)
rhs=p1*v1/t1
t2 = p2*v2/rhs

V_com = air_cycle(bore,stroke,con_rod,cr,180,0)

constant = p1*pow(v1,gamma)

P_com = []
for v in V_com:
	P_com.append(constant/pow(v,gamma))


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

V_exp = air_cycle(bore,stroke,con_rod,cr,0,180)

constant = p3*pow(v3,gamma)

P_exp = []
for v in V_exp:
	P_exp.append(constant/pow(v,gamma))
#state point 4
v4= v1
p4 = p3*pow(v3,gamma)/pow(v4,gamma)
t4 = p4*v4/rhs

#thermal efficiency
TE = 1 - (1/pow(cr,(gamma-1)))
TE = TE*100
print(\"Thermal_efficiency =\",TE)

plt.plot([v2,v3],[p2,p3])
plt.plot(V_com ,P_com)
plt.plot(V_exp ,P_exp)
plt.plot([v4,v1],[p4,p1])
plt.xlabel(\'VOLUME(Kg/m^3)\')
plt.ylabel(\'PRESSURE(N/m^2)\')
plt.show()