Week 2 Air standard Cycle
import math
import matplotlib.pyplot as plt
import numpy as np
#function to perform path tracing during compression and expansion
def engine_kinematics(c_r,v_clearance,R,cs,ce):
theta= np.linspace(ce,cs,180)
volume=[]
#the formula to calculate the volume during compression
#and expansion for path tracing is split into 3 parts
for t in theta:
term1= 0.5 * (c_r-1)
term2= R + 1 - math.cos(t)
term3= pow((pow(R,2)-pow(math.sin(t),2)),0.5)
volume.append((1 + (term1 * (term2 - term3))) * v_clearance)
return volume
#engine geometric parameters
#all dimensions in metre
bore=0.08
stroke= 0.0795
#compression ratio
c_r= 12
#ratio of connecting rod length to crank pin radius
R= 3.5
#as the crank angle changes during piston movement
#defining starting and ending crank angle
cs= 0
ce= math.pi
#specific heat ratio
gamma= 1.4
#calculation of swept volume and clearance volume
v_swept= math.pi/4 * pow(bore,2) * stroke
v_clearance= v_swept/(c_r-1)
#variables at state point 1
#pressure(in pascals) and temperature(in kelvin) are assumed
p1= 101325
t1= 300
v1= v_clearance + v_swept
#variables at state point 2
v2= v_clearance
p2= (p1 * pow(v1,gamma))/(pow(v2,gamma))
t2= (p2 * v2 * t1)/(p1 * v1)
#calling the function to trace the path during compression
v_compression= engine_kinematics(c_r,v_clearance,R,math.pi,0)
constant= p1 * pow(v1,gamma)
p_compression=[]
for v in v_compression:
p_compression.append(constant/(pow(v,gamma)))
#variables at state point 3
#assuming temperature(in kelvin) after compression
t3= 1700
v3= v_clearance
p3= (p2 * t3)/t2
#calling the function to trace the path during expansion
v_expansion= engine_kinematics(c_r,v_clearance,R,0,math.pi)
constant= p3 * pow(v3,gamma)
p_expansion=[]
for v in v_expansion:
p_expansion.append(constant/(pow(v,gamma)))
#variables at state point 4
v4= v1
p4= (p3 * pow(v3,gamma))/(pow(v4,gamma))
t4= (p4 * v4 * t3)/(p3 * v3)
#plotting the otto cycle
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.xlabel('volume(m^3)')
plt.ylabel('pressure(Pa)')
plt.show()
#calculating the thermal efficiency
therm_eff= (1- (1/(pow(c_r,0.4))))* 100
print('thermal efficiency of the cycle is ',therm_eff,' %')