Air standard Cycle Using Python

Aim : To plot PV diagram for standard otto cycle and calculate the thermal efficiency for the same

Theoretical Explanation of the otto cycle:

As we can see, otto cycle in P-V diagram has four basic steps:

1-->2 is isoentropic comprssion. No heat exchange takes place during this process. The cylinder piston moves from BDC to TDC and the volume decreases from V1 to V2. Also, the pressure increases from P1 to P2.

2-->3 is constant volume heat addition. So V2 = V3 but pressure increases from P2 to P3. In this step, the piston stays at TDC for a while and heat addition takes place.

3-->4 is isoentropic expansion. In this step, the piston moves from TDC to BDC. VOlume increase from V3 to V4 whereas the pressure increases from P3 to P4.

4-->1 is constant volume heat rejection.V4 = v1 but pressre decreases from P4 to P1.In this step , the piston stays at BDC for a short amount of time.

Steps involved in the calcuation and plotting along with the code:

#declaring input parameters

p1 = 10000 pa
t1 = 500 K

#declaring engine parameters

bore = 0.1
stroke = 0.1
c_rod = 0.15 #connecting rod length
c_ratio = 12 #compression ratio
a = stroke/2 #crank pin radius
r = c_rod/a
# calculating parameters at point 1
v_swept = (math.pi/4)*pow(bore,2)*stroke

v_clearance = v_swept/(c_ratio-1)

v1 = v_swept + v_clearance
v2 = v_clearance

# calculation parameters at point 2
gamma = 1.4;
p2 = p1*pow(c_ratio,gamma)
t2 = (p2*v2*t1)/(p1*v1)

#calculating parameters at point 3
v3 = v2
t3 = 2200
p3 = p2*(t3/t2)

#caculating parameters at point 4
v4 = v1
p4 = p3*pow((v3/v4),gamma)

Hence, we will have all the points required for plotting the PV graph.

(p1,v1) (p2,v2) (p3,v3) (p4,v4)

but, we cant simply plot these points and join them because he resulting curve will be along a straight line .We will have to use a different formula for getting curve for (volume, pressure) expansion and compression .

For this, we will create a function named engine_kinematics , to make our code look neat and to save the hassle of writing same code again and again

Volume = 1 +1/2 (c_ratio -1) {(r + 1-cosd(theta) - (r^2 - sind(theta)^2)^0.5}

where a = stroke/2

r = c_rod/a

for simplification, we will use substitution method in the formula:

term1 = 0.5*(c_ratio-1);
term2 = r+1-math.cos(theta[i]);
term3 = (r**2-math.sin(theta[i])**2)**0.5

Therefore, Volume = (1 + term1*(term2-term3))*v_clearance

and for getting multiple values for Volume expansion and volume compression, we will use different values for theta, which shall create an array of values.

for V_expansion, theta = linspace(0,180,120) [linspace function will create 120 values of theta]

hence,P_expansion = constant_c./V_expansion**gamma

Similarly for V_compression , theta = linspace(180,0,120)

and we can get P_compression also.

At the end, we shall be using plot function for plotting all the values

Now, we can calculate the thermal efficiency of the otto cycle , by using the formula:

efficiency = 1-1/c_ratio**(gamma-1), and the efficiency came out to be 62.99 %

import math
import matplotlib.pyplot as plt

def engine_kinematics(bore,stroke,con_rod,cr,start_angle,end_angle):
	"""
	Engine Kinematics
	"""
	# Engine Parameters
	a = stroke/2
	R = con_rod/a

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

	sc = math.radians(start_angle)
	ec = math.radians(end_angle)

	num_values = 50
	dtheta = (sc-ec)/(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
p1 = 101325
t1 = 500
gamma = 1.4
t3 = 2300

# engine parameters
bore = 0.1
stroke = 0.1
con_rod = 0.15
cr = 8

# Volume Computation
v_s = (math.pi/4)*pow(bore,2)*stroke
v_c = v_s/(cr-1)
v1 = v_s + v_c

# Stage2
v2 = v_c

p2 = p1*pow(v1,gamma)/pow(v2,gamma)
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))

# Stage3
v3 = v2
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))

# Stage4
v4 = v1

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

t4 = p4*v4/rhs

plt.plot([v2,v3],[p2,p3])
plt.plot(v_compression,p_compression)
plt.plot(v_expansion,p_expansion)
plt.plot([v4,v1],[p4,p1])
plt.show()

# Thermal efficiency
n = 1 - (pow(cr,(1-gamma)))
print(n)

Output: