Air standard Cycle(Otto cycle) using PYTHON

Objective: -  

Develop a code that can solve an otto cycle and make plots for it.

Principle: -

An Otto cycle is an idealized thermodynamic cycle that describes the functioning of a typical spark-ignition engine. It is the thermodynamic cycle most commonly found in automobile engines. The Otto Cycle is an air-standard cycle that approximates the processes in petrol or diesel engines. It is based on constant volume heat addition (combustion) and heat rejection processes, and Isentropic compression and expansion

There are three air-standard cycles:

• constant volume ‘combustion’ (Otto)
• constant pressure ‘combustion’ (Diesel)
• dual ‘combustion’ – this is a combination of constant volume and constant pressure combustion and results in a slightly more realistic cycle.

Solution Procedure: -

The basic Otto cycle has four processes:

•Isentropic compression

•Constant volume heat addition

•Isentropic expansion 

•Constant volume heat rejection

Processes: -

• Process 0 - 1 Intake Stroke
• Process 1 - 2 Compression Stroke
• Process 2 - 3 Ignition Phase
• Process 3 - 4 Expansion Stroke
• Process 4 - 1 Idealized Heat Rejection
• Process 1 - 0 Exhaust Stroke

Governing Equations: -

• Change of Volume Inside the combustion chamber of an engine -

        `(V)/(Vc)=1+0.5(CR-1)[R+1-cosθ-(R^2-sin^2 θ)^0.5]`

• Isentropic Process

        `PV^r=constant`

• Isenchoric Process

        `P/T = constant`

• Compression Ratio

       `CR = (Vs + Vc )/(Vc)`

• Thermal Efficiency 

         `η=1-1/(RC)^(γ-1)`

 CODE for the Function used: - 

def engine_kinematics(bore , stroke, con_rod , cr ,start_crank , end_crank):
	
	# a is crank pin radius
	a=stroke/2

	#R is the ratio for the connecting rod to the crank pin radius
	R = con_rod/a

    #swepth volume and clearance volume
	v_swept =math.pi*(1/4)*pow(bore,2)*stroke
	v_c= v_swept/(cr-1)


	st_cr=math.radians(start_crank)
	en_cr = math.radians(end_crank)

	num_val= 50

	dtheta = (en_cr - st_cr)/(num_val - 1)
	V=[]

	for i in range (0,num_val):
		theta = st_cr + i*dtheta

		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)
        #Equation
		V.append((1+term1*(term2 -term3))*v_c)
		
	return V

Explanation of the above code: -

• def engine_kinematics(bore, stroke, con_rod, cr,start_crank, end_crank) -  function is used to make the code easy to implement, it  is a group of statements that together  perform a task

• # - comment makes a program user friendly

• a=stroke/2 - a is crankpin radius 

• R = con_rod/a - R is the ratio for the connecting rod to the crank pin radius

• #swepth volume and clearance volume - comments

• v_swept =math.pi*(1/4)*pow(bore,2)*stroke - calculation of the swept volume(m^3)

• v_c= v_swept/(cr-1)  - calculation of the clearance volume (m^3)

• st_cr=math.radians(start_crank) - initializing value in radians to crank start angle

• en_cr = math.radians(end_crank) -  initializing value in radians to crank end angle

• num_val= 50 - initializing the num val

• dtheta = (en_cr - st_cr)/(num_val - 1) 

• V=[] - assigning a volume array

• for i in range (0,num_val): - execution of the statements in the for loop till the value

• theta = st_cr + i*dtheta - initializing the value of the theta

• term1 = 0.5*(cr-1) - term 1 according to the equation 

• term2 = R + 1-math.cos(theta)  - term 2 according to the equation 

• term3 = pow(R,2) - pow(math.sin(theta),2)  - term 3 according to the equation 

• term3 = pow(term3,0.5)  - term 3 according to the equation 
• #Equation - comments
• V.append((1+term1*(term2 -term3))*v_c) - append the calculated value in the V array 
• return V - return the value of the V array as output when the function is called

Code of the Otto cycle Program including the above function: -

#Otto cycle simulator
import math
import matplotlib.pyplot as plt

def engine_kinematics(bore , stroke, con_rod , cr ,start_crank , end_crank):
	
	# a is crank pin radius
	a=stroke/2

	#R is the ratio for the connecting rod to the crank pin radius
	R = con_rod/a

    #swepth volume and clearance volume
	v_swept =math.pi*(1/4)*pow(bore,2)*stroke
	v_c= v_swept/(cr-1)


	st_cr=math.radians(start_crank)
	en_cr = math.radians(end_crank)

	num_val= 50

	dtheta = (en_cr - st_cr)/(num_val - 1)
	V=[]

	for i in range (0,num_val):
		theta = st_cr + i*dtheta

		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)
        #Equation
		V.append((1+term1*(term2 -term3))*v_c)
		
	return V

#inputs

gamma =1.4
t3 =2300
#State variable 1
p1 = 101325
t1 = 500

#Engine geometric parameters
bore= 0.1
stroke = 0.1
con_rod = 0.15 #lenght of the connecting rod
cr= 12

#calculating the swept volume and clearence volume

v_swept =(math.pi/4)*pow(bore,2)*stroke
v_clearance= v_swept/(cr-1)
v1= v_clearance+v_swept

#state point 2
v2 = v_clearance

#v1/v2 is the compression ratio

'''p2v2^gamma = p1v1^gamma 
p2 = p1*(v1/v2)^gamma = p1*cr^gamma'''

p2= p1*pow(cr,gamma)


#p1v1/t1 = p2v2/t2 | t2 =p2*v2*t1/(p1*v1) 
#temperature at state point 2 is t2
t2 =p2*v2*t1/(p1*v1)
v_compression = engine_kinematics(bore,stroke,con_rod,cr,180,0)
constant_c = p1*pow(v1,gamma)

p_compression = []

for v in v_compression:
	p_compression.append(constant_c/pow(v,gamma))


#state point 3 determing pressure
v3=v2

#p3v3/t3=p2v2/t2 | p3 =p2*t3/t2
p3 = p2*t3/t2

v_expansion = engine_kinematics(bore,stroke,con_rod,cr,0,180)
constant_exp = p3*pow(v3,gamma)

p_expansion = []

for v in v_expansion:
	p_expansion.append(constant_exp/pow(v,gamma))

# state variables at state point 4
v4=v1

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

#p4v4/t4 = p3v3/t3

t4 = (p4*v4*t3)/(p3*v3)

#Thermal Efficiency
x=1-gamma
thermal_effi = 1-pow(cr,x)
print('Thermal efficiency = ',thermal_effi)

#plotting

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')
plt.ylabel('Pressure')
plt.title('Air standard cycle(otto cycle)')
plt.show()

Explanation of the above code: -

• #Otto cycle simulator - comments
• import math - importing the math 
• import matplotlib.pyplot as plt - importing matplot library
• #inputs -comments
• gamma =1.4 - initializing gamma value as 1.4
• t3 =2300 - initializing to as 2300 kelvin

#State variable 1 - comments

• p1 = 101325 - 101325  initializing p1 as 101325 pascal
• t1 = 500 -  initializing t1 as 500 kelvin

#Engine geometric parameters - comments

• bore= 0.1 - initializing with 0.1 meter
• stroke = 0.1 - initializing with 0.1 meter
• con_rod = 0.15 #lenght of the connecting rod - initializing connecting rod length as 0.15 meter
• cr= 12 - initializing compression ration as 12

#calculating the swept volume and clearence volume - comments 

• v_swept =(math.pi/4)*pow(bore,2)*stroke -  calculating the swepth volume
• v_clearance= v_swept/(cr-1) - calculating the clearance volume
• v1= v_clearance+v_swept - calculating v1 volume

#state point 2

• v2 = v_clearance -calculating v2 volume

#v1/v2 is the compression ratio - comments

'''p2v2^gamma = p1v1^gamma
p2 = p1*(v1/v2)^gamma = p1*cr^gamma''' - comments 

• p2= p1*pow(cr,gamma) - calculating p2 pressure

#p1v1/t1 = p2v2/t2 | t2 =p2*v2*t1/(p1*v1) -  comments for making the program user friendly
#temperature at state point 2 is t2 - comments

• t2 =p2*v2*t1/(p1*v1) - calculating temperature t2
• v_compression = engine_kinematics(bore,stroke,con_rod,cr,180,0) -  calling of the function to calculate compression volume
• constant_c = p1*pow(v1,gamma) - calculating constant c
• p_compression = [] - assigning the pressure compression array
• for v in v_compression: - execution of the statements in the for loop till the value
• p_compression.append(constant_c/pow(v,gamma)) - append the calculated value in the compression pressure array 

#state point 3 determing pressure - comments

• v3=v2 initializing volume v3  to volume v2 

#p3v3/t3=p2v2/t2 | p3 =p2*t3/t2 - comments

• p3 = p2*t3/t2 - calculating p3 pressure
• v_expansion = engine_kinematics(bore,stroke,con_rod,cr,0,180) - calling of the function to calculate expansion volume
• constant_exp = p3*pow(v3,gamma) -calculating constant
• p_expansion = [] - assigning the pressure expansion array
• for v in v_expansion: - for executing a given statements
• p_expansion.append(constant_exp/pow(v,gamma)) - append the calculated value in the expansion pressure array 

# state variables at state point 4

• v4=v1 initializing volume v4  to volume v1 

#p3v3^gamma = p4v4^gamma | p4 = p3(v3/v4)^gamma

• k=v3/v4 - ratio of v3 to v4 used to calculate p4 pressure
• p4 = p3*pow(k,gamma) -  calculating p4 pressure

#p4v4/t4 = p3v3/t3

• t4 = (p4*v4*t3)/(p3*v3) - t4 temperature calculation

#Thermal Efficiency

• x=1-gamma - initializing x as 1-gamma to use in calculating thermal efficiency
• thermal_effi = 1-pow(cr,x) - calculating thermal efficiency
• print('Thermal efficiency = ',thermal_effi) - printing thermal efficiency

#plotting

• plt.plot([v2,v3],[p2,p3]) - plotting the volume and pressure v2 and v3 and p2 and p3
• plt.plot(v_compression, p_compression) -  plotting the compression volume and compression pressure 
• plt.plot(v_expansion, p_expansion) - plotting the expansion volume and expansion pressure 
• plt.plot([v4,v1],[p4,p1]) - plotting the volume and pressure v1 and v4 and p1 and p4
• plt.xlabel('Volume') - labeling the x axis
• plt.ylabel('Pressure') -  labeling the y axis
• plt.title('Air standard cycle(otto cycle)') -  labeling the title of the plot
• plt.show() - showing the plot in the output screen

Output: -

 Thermal efficiency in the command window 

Thermal efficiency = 0.6298 or 62.98%

The output of the otto cycle plot: -

Conclusion and Discussion: -

• Thermal Efficiency was calculated to be 0.6299 or 62.99%
• The graph was plotted between Pressure and volume for an otto cycle