Week 2 Air standard Cycle

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

#User-defined function to calculate the volume variation during compression and expansion process using piston kinematics equation
def eng_kinematics(bore_dia,stroke,con_rod,cr,start_crank,end_crank):
    a=stroke/2
    R=con_rod/a
    vs=M.pi*(1/4) *pow(bore_dia,2)*stroke
    vc=vs/(cr-1)

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

    theta=np.linspace(sc,ec,50)
    
    v=[]
    
    for i in range(0,50):
        Theta=theta[i]
        term1=0.5*(cr-1)
        term2=R+1-M.cos(Theta)
        term3=pow(R,2)-pow(M.sin(Theta),2)
        term3=pow(term3,0.5)
        v.append((1+term1*(term2-term3))*vc)
    return v

#Engine geometric parameters
bore_dia=0.1
stroke=0.1
con_rod=0.15
cr=12
gamma=1.4

#Pressure and temperature values
p1=101325
t1=500
t3=2300

#Calculation of swept volume and clearance volume
vs=(M.pi/4)*(bore_dia*bore_dia)*stroke
vc=vs/(cr-1)
v1=vs+vc
v2=vc

p2=p1*pow(cr,gamma)
t2=(p2*v2*t1)/(p1*v1)

#Calling the function for compression process
V_comp=eng_kinematics(bore_dia,stroke,con_rod,cr,180,0)
C=p1*pow(v1,gamma)
P_comp=[]

for j in V_comp:
    P_comp.append(C/pow(j,gamma))
    

v3=v2
p3=(p2*t3)/t2


#Calculation of p4,v4,t4
v4=v1
p4=p3*pow((v3/v1),gamma)
t4=(p4*v4*t3)/(p3*v3)

Th_eff=1-(t4-t1)/(t3-t2)
print(Th_eff)
#Calling the function for expansion process
V_power=eng_kinematics(bore_dia,stroke,con_rod,cr,0,180)
C=p3*pow(v3,gamma)
P_power=[]

for k in V_power:
    P_power.append(C/pow(k,gamma))
    

#Plotting the theoretical pv- diagram
plt.figure()
plt.plot([v1,v2,v3,v4,v1],[p1,p2,p3,p4,p1])
plt.xlabel('Volume')
plt.ylabel('Pressure')
plt.title('Theoretical p-v diagram of an otto cycle')
plt.grid()
plt.show()

#Plotting the Actual pv-diagram 
plt.figure()
plt.plot([v2, v3],[p2, p3])
plt.plot([v4, v1],[p4, p1])
plt.plot(V_power, P_power)
plt.plot(V_comp, P_comp)
plt.xlabel('Volume')
plt.ylabel('Pressure')
plt.title('Actual p-v diagram of an otto cycle')
plt.legend(['Constant volume heat addition','Constant volume heat rejection','reversible adiabatic expansion','reversible adiabatic compression'])
plt.grid()
plt.show()