Week 2 Air standard Cycle

Aim:
To write a python program for an Otto cycle 

Objectives:
1)To solve and find out different state variables in the Otto cycle and plot its p-V diagram.
2)To calculate the thermal efficiency of the given Otto cycle.

Requirements:
The equations of Otto cycle involved in this project are as follows
Compression ratio (cr) , cr=Vs+Vc/Vc,

Swept volume (Vs) = π/4⋅d2⋅h, where d is bore diameter and h is the stroke length,
Piston Kinematic equation(V/Vc) = 1 + 1/2(cr−1)[R+1−cosθ−(R2−sin^2θ)^0.5], 
          where R = length of connecting rod/crank radius ,
Adiabatic process equation = p⋅v^γ = constant,
Ideal Gas equation = p⋅v=n⋅R⋅T,
Thermal efficiency of Otto cycle(η)= 1−1/(r^γ−1),
       r=compression ration which is the ratio of maximum volume to minimum volume 

Procedure:
>Intially import the values from the python dictionary such as import numpy as np,
import math and import matplotlib.pyplot as plt to perform mathematical operations and to plot the graphs.
>Take the input values of an air standard cycles such as Temperature(T), Pressure(Pa) and specific heat(gamma=1.4)
>As the geometrical parameters such as bore,stroke,connecting rod,compression ratio,are known then the swept volume and clearance volume can be calculated.
>start crank angle and end angle are taken are taken in the range of 0 to 50.
>As the state values before compression and the maximum temperature in the chamber are known then with the help of gas equation and adiabatic equation we can calcualte the other state variables.
>With the help of plot command (plt.plot)if we plot the P-V diagram we can observe a straight line for compression as well as expansion which is unusual.
>So to overcome this error we have to create an empty volume array for storing the different values of crank angle rotation for both compression and expansion. 
>This fuction is then called to the main function to plot various points for both compression and expansion.
>The number of values for which the program to be iterated is considered as 50 and an empty array for V is created and using for loop a new value will be generated and stored in the array using the append command.After every iteration a new V value will be returned back with the return command until the anticipated value is reached.
>The theramal efficency of the engine is calculated from the stated formula above.


Python code:
#Engine_kinematics
import numpy as np
import math
import matplotlib.pyplot as plt


def Engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank):

	#geometric_parameters

	a=stroke/2  #crank pin radius
	R=con_rod/a  # constant 

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

	sc=math.radians(start_crank)
	ec=math.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-math.cos(Theta)
		term3=pow(R,2)-pow(math.sin(Theta),2)
		term3=pow(term3,0.5)

		V.append((1+term1*(term2-term3))*V_c)

	return V

#inputs

P_1=101325
T_1=500
gamma=1.4
T_3=2300

#geometric_parameters
bore=0.1
stroke=0.1
con_rod=0.15
cr=12
a=stroke/2  #crank pin radius
R=con_rod/a  # constant 

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

#at state point 1

V_1=V_c+V_s #volume at state point1 which is sum of clearance volume and swept volume

# state point 2

V_2=V_c # after compression the volume at state point at 2 is equal to clearance volume 

P_2=P_1*pow(V_1,gamma)/pow(V_2,gamma)

#by ideal gas law p2v2/t2=p1v1/t1
rhs=P_1*V_1/T_1
T_2=P_2*V_2/rhs

#to plot the compression curve 

V_compression=Engine_kinematics(bore, stroke, con_rod, cr, 180 ,0)
constant=P_1*pow(V_1,gamma)
P_compression=[]
for v in V_compression:
	P_compression.append(constant/pow(v,gamma))

print(P_compression)
	
# state point 3
V_3=V_2           # constant volume heat addition after combustion so V3=V2
#by ideal gas law p3v3/t3=p2v2/t2 | p2v2/t2=rhs
P_3=(rhs*T_3)/V_3

#to plot the expansion curve

V_expansion=Engine_kinematics(bore, stroke, con_rod, cr, 0,180)
constant=P_3*pow(V_3,gamma)
P_expansion=[]
for v in V_expansion:
	P_expansion.append(constant/pow(v,gamma))

print(P_expansion)

# state point 4

V_4=V_1           # after expansion the volume at state point at 4 is equal to clearance volume +swept volume (V1)

P_4=P_3*pow(V_3,gamma)/pow(V_4,gamma)

#by ideal gas law p4v4/t4=p3v3/t3
rhs=P_3*V_3/T_3
T_4=P_4*V_4/rhs

# thermal efficiency calcualtion
thermal_efficiency=1-(1/pow((cr),(gamma-1)))
print('Thermal efficiency of otto cycle is :',thermal_efficiency)

plt.plot([V_2,V_3],[P_2,P_3])
plt.plot(V_compression,P_compression)
plt.plot(V_expansion,P_expansion)
plt.plot([V_4,V_1],[P_4,P_1])
plt.ylim([-500000,6000000])
plt.xlim([0,0.001])
plt.xlabel('Volume (m^3)')
plt.ylabel('Pressure(Pa)')
plt.title('P-V diagram of an Otto Cycle')
plt.grid()
plt.show()

Output:
1)P-V diagram for an Otto cycle can be observed.
2)Thermal efficiency of otto cycle is : 0.6298928275128466