Air Standard cycle Simulator

# Otto cycle pv diagram
import math
import matplotlib.pyplot as plt

# defining engine kinematics function

def engine_kinematics(bore, stroke, con_rod, cr, start_crank, end_crank):
	a = stroke/2
	r = con_rod/a

	# swept volume and clearence
	v_s = math.pi*(1/4)*pow(bore,2)*stroke
	v_c = v_s/(cr-1)

	# starting and ending crank angle

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

	num_values = 65

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

	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
t3 = 2300
gamma = 1.4
 
#state 1 variables||p1,v1,t1
p1 = 101325
t1 = 500

# Engine Geometry
bore = 0.1
stroke = 0.1
con_rod = 0.15
cr = 12
 
# to calculate v1
# swept volume = (pi/4)*bore^2*stroke
v_s = (math.pi/4)*pow(bore,2)*stroke
v_c = v_s/(cr-1)
v1 = v_s + v_c

# state 2 variable || p2,v2,t2
v2 = v_c
# to calculate p2 | p2*v2^gamma = p1*v1^gamma,  p2 = p1*(cr^gamma)
p2 = p1*pow(cr,gamma)
# to calculate t2 | p2*v2/t2 = p1*v1/t1,  t2 = p2*v2*t1/p1*v1
t2 = p2*v2*t1/(p1*v1)

# actual pressure & volume while compression
constant = p1*pow(v1,gamma)
v_comp = engine_kinematics(bore, stroke, con_rod, cr, 180,0)
p_comp = []

for v in v_comp:
	p_comp.append(constant/pow(v,gamma))


# state 3 Vraiable || p3,v3,t3
v3  = v2
# to calculate p3 || p2*t3 = p3*t2,   p3 = p2*t3/t2
p3 = p2*t3/t2

# actual pressure & volume while expansion
constant = p3*pow(v3,gamma)

v_expn = engine_kinematics(bore, stroke, con_rod, cr, 0,180)

p_expn = []

for v in v_expn:
	p_expn.append(constant/pow(v,gamma))


# state 4 variables  || p4,v4,t4
v4 = v1
# to calculate p4 | p4*v4^gamma = p3*v3^gamma,  p4 = p3*((v3/v4)^gamma)
p4 = p3*pow((v3/v4),gamma)
# to calculate t4 | p4*v4/t4 = p3*v3/t3,  t4 = p4*v4*t3/p3*v3
t4 = p4*v4*t3/p3*v3


#ploting   figure 1 

plt.plot([v2,v3],[p2,p3])
plt.plot([v4,v1],[p4,p1])
plt.plot(v_comp,p_comp)
plt.plot(v_expn,p_expn)
plt.show()

coding for otto cycle pv diagram,

step 1)

importing math module and matplotlib.pyplot

step 2)

asuming predefined inputs

t3 = ignition tempeture= 2300 in kelvin

gamma = 1.4 = constant

here engine geometry are

bore = diameter of cylinder = 0.1 in metre

stroke = lenght of cylinder = 0.1 in metre

con_rod = connecting rod length = .15 in metre

cr = compression ratio = v1/clearence volume

a = crank pin radius = stroke/2

r = con_rod/2

step 3)

swept volume = pi/4*bore^2*stroke

clearence volume = swept volume/(cr-1)

v1 = volume at state 1 = swept volume + clearence volume.

step 4)

As otto cycle has four process

1) isentropic compression

state 1 variables are

p1 = pressure in pascal = 101325 which is atmospheric pressure

t1 = temperature in kelvin = 500

v1= volume in m^3 = swept volume + clearence volume

2) constant volume heat addition

state 2 variables are

p2 = pressure

t2 = temperature

v2= volume = clearence volume

to calculate p2 we can use p1*v1^gamma = p2*v2^gamma.

p2 = p1*(cr^gamma)

to calcuate t2 we can use (p1*v1)/t1 = (p2*v2)/t2 as rest of all varaibles aren known

t2 = t2 = p2*v2*t1/(p1*v1)

(3) isenttropics expansion

state 3 variables are

p3 = pressure

t3 = temperature in kelvin = 2300

v3= volume

4) constant volume heat rejection.

state 4 variables are

p4 = pressure

t4 = temperature

v4= volume

now state 3 & state 4 varaibles can be calculated same as p2 and t2.

step 5)

now to obatin actual curve we have to define

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

i.e engine_kinematics is a function of bore,stroke,con_rod,cr,start_crank,end_crank

here starting and ending crank angle is a array , for each process crank rotate 180deg

so

starting and ending crank angle

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

interval between 2 consecative values can be calculate as follow

num_values = 65

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

now the below term can be break into 3 parts to make it simple

geometry v/vc = 1 + 0.5(rc-1)[r+1-cosd(theta)-(r^2-sind(theta)^2)^0.5]

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)

volume
V = (1 + term1*(term2-term3))*v_c

now we can write volume array as

v.append(((1 + term1*(term2 - term3))*v_c))

step 6)

actual compression curve

we need v_compression and p_compression

v_comp = engine_kinematics(bore,stroke,con_rod,cr,180,0)

as in compression stroke crank rotate from 180 to 0 deg

here we need constant to calculate

constant = p1*v1.^gamma

now we can calculate p_compression using

p_comp = (constant_c./v_comp.^gamma)

can be calculated as follow

p_comp = []

for v in v_comp:
p_comp.append(constant/pow(v,gamma))

same we can do for actual expansion curve

actual expansion curve

we need v_expansion and p_expansion

v_expn = engine_kinematics(bore,stroke,con_rod,cr,0,180)

as in expansion stroke crank rotate from 0 to 180 deg

here we need constant to calculate

constant = p3*v3^gamma

now we can calculate p_expansion using

p_expn = (constant_/v_expn^gamma)

can be calculated as follow

p_expn = []

for v in v_expn:
p_expn.append(constant/pow(v,gamma))

step 7) plotting

we can plot

v_compression and p_compression and v_expansion and p_expansion

and we can connect straight line between them.

plt.plot([v2,v3],[p2,p3])
plt.plot([v4,v1],[p4,p1])
plt.plot(v_comp,p_comp)
plt.plot(v_expn,p_expn)
plt.show()

step 8)

Thermal efficiency can be calculated using

|| Thermal_eff = 1-1/cr^(gamma-1) ||

Thermal_eff = (1-1/cr^(gamma-1))*100

= 62.98 %