Plotting the P-V curve and calculating thermal efficiency of an Otto cycle (Using Python)

Objective: To plot the P-V graph of an Otto cycle and calculate its thermal efficiency using Python.

Assumptions

1. The temperature at state 3, i.e., after heat is added to the cycle at constant volume, is given as an input as it is assumed to be known initially.
2. All air standard cycle assumptions are valid throughout the cycle.

Procedure

Step 1: Two major kinds of parameters are given as inputs, namely the state variables (pressure and temperature at state point 1, and the temperature at state point 3) along with the ratio of specific heats and also the engine geometric parameters (bore diameter, stroke length, and length of the connecting rod). Additional parameters given as inputs are the compression ratio as well as the initial and final angular positions of the crank.

Step 2: The various state points 2,3 and 4 can be calculated by using the isentropic relation and constant volume relation by the following formulae:

`pV^gamma = constant`

Which is valid for the processes 1-2 and 3-4, which are isentropic compression and expansion processes respectively.

`Volume = constant`

Which is valid for the processes 2-3 and 4-1, which are constant volume heat addition and heat rejection processes respectively.

In order to correctly capture the isentropic processes on the final P-V plot, a function called Piston_Kinematics() is used, in which the relation between the volume covered by the piston during the compression and expansion processes and the clearance volume is expressed as a function of the crank angle `theta`. The relation is given by:

`V/V_c = 1 + 1/2 (r_c -1)[R + 1 - costheta - (R^2 - sin^2theta)^(1/2)]`

Where `V` is the volume covered by the piston, `V_c` is the clearance volume, `r_c` is the compression ratio, `theta` is the angular position of the crank, and `R` is defined as

`R = (2* Conrod)/ (Stroke)`

Where `Conrod` is the length of the connecting rod and `Stroke` is the stroke length.

After obtaining the corresponding values of `V` for both the compression and expansion processes, the corresponding values of `P` are calculated by using the isentropic relation mentioned earlier. Various values of `theta` taken into consideration by using the linspace function. Thus, the corresponding calculation for `V_compression` and `P_compression` are taken as arrays.

Step 3: All the state points are plotted on a `P-V` diagram, with `(P_compression, V_compression)` and `(P_expansion, V_expansion)` plotted instead of `(p_3, V_3)` and `(p_4, V_4)` respectively.

The thermal efficiency is also calculated using the formula:

`eta = 1 - 1/(r_c)^(gamma-1)`

Where `eta` is the thermal efficiency, `r_c` is the compression ratio (denoted throughout in the code as `cr`) and `gamma` is the ratio of specific heats.

Common errors committed and how they are overcome

1. Not giving the starting and ending crank angles in radians. This can cause issues in calculating the total volume, which in turn, would affect the state points 2,3 and 4. This can be avoided by using the math.radians().

2. The incorrect order of inputs given to the engine_kinematics(). This can give incorrect outputs, especially when the knowledge of the compression and expansion strokes is not known to the user.

3. Returning the value of V within the for-loop. This causes Python to return nothing, as there is no output for V. Hence, V must be returned outside the for-loop but within the engine_kinematics() function definition.

4. Typing errors can be avoided by careful viewing of the code.

Conclusion

Thus, it can be concluded that the Otto cycle follows the plot as shown below, and the thermal efficiency of the same is also calculated using Python by the equation provided to it. It can also be seen that the thermal efficiency of the Otto cycle is only dependent on the compression ratio and the ratio of specific heats, which is an indicator that it does not depend on the initial and final thermodynamic states of the cycle.

#Program to simulate an Otto cycle

import math
import matplotlib.pyplot as plt 

def engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank):
	#Engine kinematics
	import math

	#Geometric parameters

	bore = 0.1
	stroke = 0.1
	con_rod = 0.15
	cr = 12
	a = stroke/2
	R = con_rod/a

	# Volume parameters
	V_s = math.pi/4*pow(bore,2)*stroke
	V_c = V_s/(cr-1)

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

	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)
		term4 = pow(term3,0.5)
		V.append((1 + term1*(term2 - term4))*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

#p2v2^gamma = p1v1^gamma
p2 = p1*pow(v1,gamma)/pow(v2,gamma)

#p2v2/t2 = p1v1/t1 | rhs = p1v1/t1 | p2v2/t2 = rhs | t2 = p2v2/rhs
rhs = p1*v1/t1
t2 = p2*v2/rhs
V_compression = engine_kinematics(bore,stroke,con_rod,cr,180,0)

constant_comp = p1*pow(v1,gamma)
#constant = P_compression*(V_compression^gamma)

P_compression = []
print(V_compression)
for V in V_compression:
	P_compression.append(constant_comp/pow(V,gamma))

#state point 3
v3 = v2

#p3v3/t3 = p2v2/t2 | rhs = p2v2/t2 | p3 = rhs*t3/v3

rhs = p2*v2/t2
p3 = rhs*t3/v3

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

constant_exp = p3*pow(v3,gamma)
#constant = P_compression*(V_compression^gamma)

P_expansion = []

for V in V_expansion:
	P_expansion.append(constant_exp/pow(V,gamma))

#state point 4
v4 = v1

#p4v4^gamma = p3v3^gamma

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

#p4v4/t4 = p3v3/t3
t4 = p4*v4/rhs

#Thermal efficiency of the engine
frac = pow(1/cr,gamma-1)
therm_eff = 1 - frac
print('The thermal efficiency of the engine is {} percentage'.format(therm_eff*100))

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

#Function definition of the engine_kinematics()
import math

#Geometric parameters

bore = 0.1
stroke = 0.1
con_rod = 0.15
cr = 12
a = stroke/2
R = con_rod/a

#Volume parameters
V_s = math.pi/4*pow(bore,2)*stroke
V_c = V_s/(cr-1)

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 = (1 + term1*(term2 - term3))*V_c
print(V_s)
print(V)