Plotting the PV Diagram for an Otto Cycle and Calculating its Thermal Efficiency

Objective: To write a program in python which can calculate the state properties at different states of the otto cycle and hence, plot the PV diagram. The program should also be able to display the thermal efficiency of the cycle.

Procedure:-

• We write a program in which based on the input and geometric parameters, the state properties like pressure, volume and temperature are calculated at different states of the otto cycle.
• By using these state properties, the PV diagram is plotted.
• Finally, the thermal efficiency of the engine is evaluated using the temperature values. We know, the thermal efficiency,
  • `eta = 1 - frac{Q_c}{Q_h}` where `Q_c`is the heat dissipated to the environment and `Q_h`is the heat supplied through the spark plug.
  • But we know that `Q_c = mc_v(T_4 - T_1)` and `Q_h = mc_v(T_3 - T_2)`.
  • Therefore, `eta = 1 - frac{T_4 - T_1}{T_3 - T_2}`
• The code along with the explanation of different steps has been described below:

# Otto cycle simulator

import math
import matplotlib.pyplot as plt

# Defining a function which returns the values of volume and pressure at different points on the adiabatic process line
def engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank):
	"""
	Engine kinematics
	"""

	# Geometric parameters
	a = stroke/2
	R = con_rod/a

	# Volume parameters

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

	# Converting the start_crank and end_crank angles to radians
	sc = math.radians(start_crank)
	ec = math.radians(end_crank)

	# Calculating volume and pressure values at 25 different points on the adiabatic process line

	num_values = 25

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

	V = []

	for i in range(0,num_values):
		theta = sc + i*dtheta
		# Using the formula to evaluate volume using the crank angle and the geometric parameters
		term1 = 0.5*(cr - 1)
		term2 = R + 1 - math.cos(theta)
		term3 = pow(R,2) - pow(math.sin(theta),2)
		Volume = V_c*(1 + term1*(term2 - pow(term3,0.5)))
		V.append(Volume)

	return V


# Inputs

# Pressure: unit in Pascal
# Lets take this pressure as atmospheric pressure which is 101325 Pa.
p1 = 101325
# Temperature at state 1: unit in Kelvin
# Lets take this temperature as 300 K
t1 = 300
gamma = 1.4
# Temperature at state 3: unit in Kelvin
# Lets take this temperature as 2500
t3 = 2500

# Geometrical parameters
# Lets take bore, stroke and connecting rod to be 0.14 m, 0.18 m and 0.25 m respectively.
bore = 0.14
stroke = 0.18
# connecting rod
con_rod = 0.25
# compression ratio: Lets take this as 9
cr = 9

# Volume computation

# stroke volume
v_s = math.pi*(1/4)*pow(bore,2)*stroke
# clearance volume
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 | p1v1/t1 = rhs | t2 = p2v2/rhs
rhs = p1*v1/t1
t2 = p2*v2/rhs

# Volume and pressure at different points on the compression process line (State 1 -> State 2)
V_compression = engine_kinematics(bore,stroke,con_rod,cr,180,360)

P_compression = []

constant = p1*pow(v1,gamma)

for V in V_compression:
	Pressure = constant/pow(V,gamma)
	P_compression.append(Pressure)

# state point 3
v3 = v2

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

# state point 4
v4 = v1

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

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

# Volume and pressure at different points on the expansion process line (State 3 -> State 4)

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

P_expansion = []

constant = p3*pow(v3,gamma)

for V in V_expansion:
	Pressure = constant/pow(V,gamma)
	P_expansion.append(Pressure)

# Plotting each process line one by one

# Adiabatic compression from state 1 to state 2
plt.plot(V_compression,P_compression)
# Isochoric process from state 2 to 3
plt.plot([v2,v3],[p2,p3])
# Adiabatic expansion process from state 3 to 4
plt.plot(V_expansion,P_expansion)
# Isochoric process from state 4 to 1
plt.plot([v4,v1],[p4,p1])
plt.xlabel('Volume (in m^3)')
plt.ylabel('Pressure (in Pa)')
plt.show()

thermal_efficiency = 1 - ((t4 - t1)/(t3 - t2))
print('The thermal efficiency of the otto cycle is:')
print(thermal_efficiency)

Results:-

• We get the PV diagram for the otto cycle:-

• We also get the thermal efficiency of the otto cycle as can been seen in this screenshot:-