Air standard Cycle Using Python

OBJECTIVE

• To create the PV diagram of Otto Cycle for the given inputs such as engine dimensions, pressure, and temperature, using Python programming.
• To obtain the thermal efficiency of the engine.

PROCEDURE

# Otto Cycle Simulator
import math
import matplotlib.pyplot as pyplot	

def engine_kinematics(bore, stroke, con_rod, cr, start_crank, end_crank):
	"""
	Engine_Kinematics
	V/Vc = 1+0.5*(cr-1)*[R + 1 - math.cos(theta) - pow(pow(r,2)-pow(math.sin(theta),2),0.5)]
	"""
	#crank_radius,a 
	a = stroke/2
	#Ratio of connecting rod to crank, R
	R = con_rod/a

	#Volume Parameters
	Vs = math.pi*(1/4)*pow(bore,2)*stroke
	Vc = Vs/(cr - 1)

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

	num = 50
	dtheta = (ec - sc)/(num - 1)

	V  = []

	for i in range(0,num):
		theta = sc + i*dtheta

		term_1 =  0.5*(cr - 1)
		term_2 = R + 1 - math.cos(theta)	
		term_3 = pow((pow(R,2)-pow(math.sin(theta),2)),0.5)

		V.append((1 + term_1*(term_2 - term_3))*Vc)
	return V

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

# Geometric Properties
bore = 0.1
stroke =  0.1
con_rod = 0.15
cr = 12

# Volume Computation
v_s = (math.pi/4)*pow(bore,2)*stroke

# Compression Ratio = (Vs + Vc)/Vc
v_c = v_s/(cr-1)

# Point 1
v1 = v_s + v_c

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

# Point 3
v3 = v2

# p3*v3/t3 = p2*v2/t2
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))

# Point 4
v4 = v1

# p4*v4^gamma = p3*v3^gamma
p4 = p3*pow(v3,gamma)/pow(v4,gamma)

# p4*v4/t4 = p3*v3/t3
t4 = p4*v4/rhs

#Thermal Efficiency = (1 - (1/pow(cr,(gamma-1))))
thermal_eff = 1 - (1/pow(cr,(gamma-1))) 
print(thermal_eff)

pyplot.plot([v2,v3],[p2,p3])
pyplot.plot(V_compression, P_compression)
pyplot.plot(V_expansion, P_expansion)
pyplot.plot([v4,v1],[p4,p1])
pyplot.suptitle('PV Diagram')
pyplot.xlabel('Volume')
pyplot.ylabel('Pressure')
pyplot.show()

OUTPUT

Thermal_eff = 0.6298928275128466