Week 2 Air standard Cycle

Objective: To find the thermal efficiency of the Otto cycle using Python code simulation

Problem statement:

A square engine cylinder of SI engine has bore diameter 10 cm whose connecting rod length is 17.5 cm started its compression at a pressure of 1 bar at a temperature of 300 K. Peak temperature attained in the cycle seemingly around 1350 K at a compression ratio of 7.

Solve PV diagram and thermal efficiency of the Otto cycle

Solution procedure:

Fig: Schematic of Otto cycle PV diagram

Otto cycle mainly has four processes:

(1-2) Isentropic compression [ P₁ V₁^Y= P₂V₂^Y]:

piston moves from BDC to TDC i.e. crank angle shifts from 180⁰ to 0⁰

(2-3) Constant volume heat addition (V₂ = V₃)

(3-4) Isentropic expansion [ P₃ V₃^Y= P₄V₄^Y]:

Piston moves from TDC to BDC i.e. crank angle shifts from 0⁰ to 180⁰

(4-1) Constant volume heat expansion (V₄ =V₁)

V2 = clearance volume (Vc), as piston is at BDC

V1 = swept volume (Vs) + clearance volume (Vc), piston is at TDC

Piston kinematics is key factor to consider in the design of both compression and expansion processes and is explained by following equation

Non-dimensional form of cylinder volume to clearance volume ratio:

Here, r_c = compression ratio

In general,

Crank pin radius (a) = 0.5 x stroke length

R = connecting rod length

`theta` is the crank angle which defines the position of the piston and the process variables such as pressure, volume, and Temperature.

When crank angle theta = 0⁰ `=>` V = Vc, piston moves from BDC to TDC (expansion process)

At theta = 180⁰ `=>`V= Vs + Vc, piston moves from TDC to BDC (compression process)

Solution procedure:

Python code is used to simulate the Otto cycle at the given conditions

Piston kinematics is linked to compression process using function definition to get the exact PV curve.

import math
import matplotlib.pyplot as plt 

'''

Unit system:

Dimensions are in 'cm'

pressure in pascal

Temperature is in Kelvin
'''

def engine_kinematics(bore,stroke,con_rod,cr,crank_start,crank_end):

	a = stroke/2                                      # a is crank radius in general stroke length (stroke) = 2 X crank radius(a)
	R = con_rod/a                                     # R is the ratio of connecting rod length(con_rod) to crank radius(a)

	v_s = math.pi/4 * pow(bore,2)*stroke              # v_s = swept volume
	v_c = v_s/(cr-1)                                  # v_c = bore volume

	EC = math.radians(crank_end)                      # EC denotes crank ending angle  
	SC = math.radians(crank_start)                    # SC denotes crank starting angle

	num_values = 50                                   # assumed 50 divisions for crank movement between 0 to 180 and vice-versa 
	dtheta = (EC-SC)/(num_values-1)                   # dtheta is the angle of each divison

 	## solving engine kinematics equation

	v = []

	for i in range(0,num_values):

 		theta = SC + i*dtheta

 		term1 = 0.5 * (cr-1)
 		term2 = R + 1- math.cos(theta) 

 		K = pow(R,2) - pow(math.sin(theta),2) 

 		term3 = pow(K,0.5)

 		v.append((1 + term1*(term2-term3))*v_c)         # V is the volume of the cylinder occupied with respect to crank angle

	return v

#Inputs

p1 = 1                                                  # p1 is pressure at starting of the cycle
t1 = 300                                               # t1 is starting point temperature 
t3 = 1350

gamma = 1.4

stroke = 10
bore = 10  

cr = 7                                        #compression ratio                             
con_rod = 17.5                                    # length of the connecting rod 


# volume computation

v_s = math.pi/4 * pow(bore,2) *stroke

v_c = v_s/(cr-1)                                   # clearance volume



# (1-2) Isentropic compression

v1= v_c + v_s                                      # v1 indicates cylinder total volume as compression starts from BDC

# solving state point 2 variables
v2= v_c                                            # v2 indicates cylinder volume at the end of compression process->pistion reaches TDC
p2 = p1* pow(v1,gamma)/pow(v2,gamma)               


# p1v1/t1 = p2v2/t2 | rhs = p2v2/t2 = p1v1/t1 = p3v3/t3 = p4v4/t4    

rhs = p1*v1/t1

t2=p2*v2/rhs

v_compression = engine_kinematics(bore, stroke, con_rod,cr, 180,0)      # pistion moving from BDC to TDC as V_compression = f(theta)

constant = p1* pow(v1,gamma) # p_compression* v_compression^ gamma = constant

p_compression = []

for V in v_compression:
	p_compression.append(constant/pow(V,gamma))            # pressure with respect to volume compressed



# state point 3

v3 = v2

# p3v3/t3 = p2v2/t2 = rhs
rhs = p2*v2/t2

p3 = rhs *t3/v3

v_expansion = engine_kinematics(bore, stroke, con_rod,cr, 0,180) # piston moving from TDC to BDC, crank angle tends from 180 to 0 degrees

constant = p3* pow(v3,gamma) # p_compression* v_compression^ gamma = constant

p_expansion = []

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

# state point 4

v4 = v1

# rhs = p3v3/t3 = p4v4/t4

p4 = p3*pow(v3,gamma)/pow(v4,gamma)

t4 = p4*v4/rhs

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 ------> cm^3')                             # volume in cubic centimeters

plt.ylabel('Pressure ------> bar')

plt.title('P-V diagram of otto cycle')

plt.show()

# solving thermal efficiency (n= 1-1/[r^n-1])

y = gamma-1

therm_eff = 1- 1/pow(cr,y)

print('Thermal efficiency of otto cycle :',therm_eff*100,"%")

Results: