Python code to plot P-V diagram of Otto cycle and calculate it's thermal efficiency

Otto cycle processes:-

Process 1-2 (Isentropic compression)

Process 2-3 (Contant volume heat addition)

Process 3-4 (Isentropic expansion)

Process 4-1 (Contant volume heat rejection)

Themal efficiency, `eta = 1 - (1/(r^(gamma-1)))`

This code will initialize the inputs. It will also initilize the array containing crank angle at step size of `1^o` Then, it will calculate all the unknown state variables. Further, it will calculate in process(1-2) and (3-4) states at various crank angles in arrays 'theta12' and 'theta34' respectively. Finally, t will calculate thermal efficiency and plot the P-V diagram.

"""
This code plots the P-V diagram for Otto cycle
and displays it's thermal efficiency.


Processes:-

1-2  Isentropic Compression
2-3  Constant volume heat addition
3-4  Isentropic Expansion
4-1  Constant volume heat rejection
"""

import math
import numpy as np
import matplotlib.pyplot as plt

# Inputs
 
p1 = 101325        	 # Pressure at state 1 (Pa)
t1 = 500      	 # Temperature at state 1 (K)
t3 = 2300     	 # Temperature at state 3 (K)
cr = 12       	 # Compression ratio
d = 0.1       	 # Bore diameter (m)
L = 0.1       	 # Stroke length (m)
gamma = 1.4   	 # Ratio of specific heats
con_rod = 0.15    # Connecting rod length (m)
crank = L/2 # Crank length
R = con_rod/crank
theta12 = np.radians(np.arange(180,361,1)) # Crank angles for process 1-2 in steps of 1 degree
theta34 = np.radians(np.arange(0,181,1)) # Crank angles for process 3-4 in steps of 1 degree

# State Calculations

Vs = (math.pi/4)*d*d*L                  # Swept volume (m^3)
Vc = Vs/(cr-1)        				    # Compressed volume (m^3)
V1 = Vs+Vc            					# Volume at state 1 (m^3)
V4 = V1									# Volume at state 4 (m^3)
V2 = Vc                					# Volume at state 2 (m^3)
p2 = p1*math.pow((V1/V2),gamma) 		# Pressure at state 2 (Pa)
t2 = t1*(p2/p1)*(V2/V1)                 # Temperature at state 2 (K)
p3 = p2*(t3/t2)                         # Pressure at state 3 (Pa)
V3 = V2 								# Volume at state 3 (m^3)
p4 = p3*math.pow((V3/V4),gamma)			# Pressure at state 4 (Pa)
t4 = t3*(p4/p3)*math.pow((V4/V3),gamma) # Temperature at state 4 (K)

# Process Calculations

# Process 1-2
term1 = R+1-np.cos(theta12)              # R + 1 - cos(theta)
term2 = R*R - np.square(np.sin(theta12)) # R^2 - sin(theta)^2
term3 = term1 - np.sqrt(term2)
V12 = Vc * (1 + 0.5 * (cr - 1) * term3)
p12 = p1*np.power((V1/V12),gamma) 

# Process 3-4
term1 = R+1-np.cos(theta34)              # R + 1 - cos(theta)
term2 = R*R - np.square(np.sin(theta34)) # R^2 - sin(theta)^2
term3 = term1 - np.sqrt(term2)
V34 = Vc * (1 + 0.5 * (cr - 1) * term3)
p34 = p3*np.power((V3/V34),gamma)


# Efficiency Calculation
eff = (1 - (1/(math.pow(cr,gamma-1)))) * 100
print('Efficiency = {} %'.format(eff))

# Plot
plt.plot(V12,p12)
plt.plot([V2,V3],[p2,p3])
plt.plot(V34,p34)
plt.plot([V4,V1],[p4,p1])
plt.plot([V1,V2,V3,V4],[p1,p2,p3,p4],'bo')
plt.annotate('1',xy=(V1,p1),horizontalalignment='left', verticalalignment='top', fontsize=20)         
plt.annotate('2',xy=(V2,p2),horizontalalignment='right', verticalalignment='top', fontsize=20)     
plt.annotate('3',xy=(V3,p3),horizontalalignment='right', verticalalignment='bottom', fontsize=20)    
plt.annotate('4',xy=(V4,p4),horizontalalignment='left', verticalalignment='bottom', fontsize=20)   
plt.xlabel('Volume (m^3)')
plt.ylabel('Pressure (Pa)')
plt.title('Compression Ratio={}, Bore={} m, Stroke={} m'.format(cr,d,L))
plt.show()

Output:-