Air standard Cycle using Python

AIM

To write a program to solve an Otto Cycle as per the given inputs and create plot of the PV diagram as

well as display the thermal efficiency of the cycle.

Theory:

Air Standard Cycles work on assumptions and are different from the Actual Cycles.

Assumptions:

Gas mixture is air.

Air properties are used.

Closed cycle is assumed instead of the actual open cycle.

Almost constant pressure intake and exhaust processes are replaced by perfectly constant pressure intake and exhaust.

Compression and Expansion strokes are isentropic.

All processes are reversible.

Otto Cycle

Thermodynamic processes

0-1: Intake Stroke.

1-2: Compression Stroke.

2-3: Combustion Process (Const. Vol. Heat Addition).

3-4: Expansion Stroke (Power Stroke).

4-1: Heat Rejection.

1-0: Exhaust at Constant Pressure.

Thermodynamic Relations for all the processes

0-1 Intake Stroke:

P1 = P0.

1-2 Isentropic Compression Stroke:

Q1-2 = 0.

P2/P1 = (V1/V2)^k.

T2/T1 = (V1/V2)^(k-1).

2-3 Constant Volume Heat Addition:

V2=V3.

T3= Tmax.

P3 = Pmax.

3-4 Isentropic Expansion Stroke:

Q3-4 = 0.

P4/P3 = (V3/V4)^k.

T4/T3 = (V3/V4)^(k-1).

4-1 Constant Volume Heat Rejection:

V4 = V1.

1-0 Exhaust at Constant Pressure:

P0 =P1.

Inputs Given for the Otto Cycle are as follows:

State Variables:P1 = 1 atm.

T1 = 500 K.

γ = 1.4

Engine Parameters:

Bore = 0.1 m.

Connecting Rod Length = 0.15 m.

Stroke = 0.1 m.

Compression Ratio = 12.

Engine Kinematics:

Swept Volume = (π/4)*(bore)^2*stroke

Clearance Volume = V_Swept/(CR-1)

Governing Equation for piston kinematics:

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

Calculations of State Variables

State Point 2:

1-2 : Isentropic Compression

P1(V1)γ = P2(V2)γ

P1V1/T1 = P2V2/T2

State Point 3:

2-3 : Constant Volume Heat Addition

V2 = V3

T3 = 2300 K

P2V2/T2 = P3V3/T3

P3 > P2

State Point 4:

3-4 : Isentropic Expansion

P4(V4)γ = P3(V3)γ

Code of the Project:

import math
import matplotlib.pyplot as plt

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


  
	  a=stroke/2
	  n=con_rod/a
      # Volume swept by the piston that is swept volume
	  v_swept=p1*pow(bore,2)*stroke/4
	  #Clearance volume
	  v_c=v_swept/(cr-1)
	  #Defining the starrting crank angle
	  starting_crank=math.radians(start_crank)

	  #Defining the ending crank angle
	  ending_crank=math.radians(end_crank)
	  #Number of values of theta
	  num_values=100
	  #For different values of theta we need to make an empty array of v and then append it each tiome the loop iterates
	  v=[]

	  for i in range(0,num_values):
	    theta=starting_crank+i*(ending_crank-starting_crank)/(num_values-1)
	    
	 
	# Using v/v_c =1+0.5*(cr-1)*[n+1-cos(theta)-(n^2-sin^2(theta))^0.5]
	# whre v is the volume inside the cylinder at any crank angle theta
	#v_c is the clearance volume

	#breaking the equation into smaller parts
	    t1=0.5*(cr-1)
	    t2=n+1-math.cos(theta)
	    t3=pow(pow(n,2)-pow(math.sin(theta),2),0.5)
	    v.append((1+t1*(t2-t3))*v_c)
	    #In afunction we need to return the value of the entity for further use
	  return v
  
 


#Otto cycle
#Inputs
#Pressure at state 1
p1=101325
#Temperature at state 1
t1=500
#Maximum temperature of the cycle
t3=2300
gamma=1.4
bore=0.1
stroke=0.2
#Length of connecting rod
con_rod=0.15
#compression ratio
cr=12

v_swept=p1*pow(bore,2)*stroke/4
volume_clearance=v_swept/(cr-1)
#Initially the piston is at Bottom Dead centre
v1=v_swept+volume_clearance
#When the piston is at Top dead centre
v2=volume_clearance
#Since 1-2 is an adiabatic process PV^gamma=constant
p2=p1*pow(v1/v2,gamma)
#Using the ideal gas equation
t2=p2*v2*t1/(p1*v1)

#At state point 3
v3=v2
p3=p2*v2*t3/(v3*t2)
#At state point 4

v4=v1
#Since 3-4 is an adiabatic process
p4=p3*pow(v3/v4,gamma)
#Using the ideas gas equation
t4=p4*v4*t3/(p3*v3)
#Calling the function to calculate the value of v_compression for different crank angles
v_compression=engine_kinematics(0.1,0.2,0.15,12,180,0)
#Calculating the value of constant using constant=P1*V1^gamma
c1=p1*pow(v1,gamma)
#For each vale of volume compression calculating the value of p_compression
p_compression=[]
for volume in v_compression:
     p_compression.append(c1/pow(volume,gamma))
     
#Similarly calling the function to calculate the value of v_expansion for different crank angles
v_expansion=engine_kinematics(0.1,0.2,0.15,12,0,180)
c3=p3*pow(v3,gamma)
p_expansion=[]
#For each vale of volume expansion calculating the value of p_expansion
for volume in v_expansion:
  p_expansion.append(c3/pow(volume,gamma))
#Plotting the graph to visualise the cycle
plt.plot(v_compression,p_compression)
plt.plot([v2,v3],[p2,p3])
plt.plot(v_expansion,p_expansion)
plt.plot([v4,v1],[p4,p1])
plt.xlabel('Volume')
plt.ylabel('Pressure')
#Shows the plot
plt.show()
Efficiency_of_cycle=1-pow(cr,-(gamma-1))
print('Efficiency_of_cycle=',Efficiency_of_cycle)

PLOT and Result