Week 2 Air standard Cycle

import math
import matplotlib.pyplot as plt

def engine_kinematics(bore, strock, con_rod, cr, start_crank, end_crank):
	a = strock/2
	R = con_rod/a 

	V_s = math.pi*(1/4)*pow(bore,2)*strock
	V_c = V_s/(cr-1)

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

	num_values = 500

	dtheta = (ec-sc)/(num_values-1)
	V = []
	for i in range (0,num_values):
		theta = sc + i*dtheta
		term_1 = 0.5* (cr -1)
		term_2 = R + 1 - math.cos(theta)
		term_3 = pow(R,2) - pow(math.sin (theta), 2)
		term_3 = pow(term_3,0.5)
		V. append ((1 + term_1* (term_2 -term_3))*V_c)
	return V	


p1 = 101325
t1 = 500
gamma = 1.4 
t3 = 2300

bore = 0.1
strock = 0.1
con_rod = 0.15
cr = 8		

# volume computation 
v_s = (math.pi/4)*pow(bore,2)*strock
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)

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

V_commpression = engine_kinematics(bore,strock,con_rod,cr,180,0)
constant = p1*pow(v1,gamma)
P_commpression = []
for v in V_commpression:
	P_commpression.append(constant/pow(v,gamma))

#stat point 3
v3 = v2
rhs = p2*v2/t2
p3 = rhs*t3/v3
V_expansion = engine_kinematics(bore,strock,con_rod,cr,0,180)
constant = p3*pow(v3,gamma)
P_expansion = []
for v in V_expansion:
	P_expansion.append(constant/pow(v,gamma))

#stat point 4
v4 = v1
p4 = p3*pow(v3,gamma)/pow(v4,gamma)
t4 = p4*v4/rhs

plt.plot((v2,v3),(p2,p3))
plt.plot(V_commpression,P_commpression)
plt.plot(V_expansion,P_expansion)
plt.plot((v4,v1),(p4,p1))
plt.show()


# Thermal efficincy 
thermal_efficincy = 1- (1/pow( (cr), (gamma-1)))
print(thermal_efficincy)​