Air standard cycle using Python

Aim:

Write a program to plot the PV diagram of air standard Otto cycle for an engine and calculate the thermal efficiency of the engine.

Theory:

Otto cycle is the gas power cycle used in SI engine. This cycle consists of four processes. Two isentropic processes and two isochoric processes.

Process 1: Isentropic compression

The piston moves from BDC to TDC, while the air in the engine cylinder undergoes reversible adiabatic compression. The properties of air while the piston is at BDC gives the state point 1, properties of air at TDC gives the state point 2. From point 1-2, the volume decreases and pressure increases.

`p1.v1^gamma =p2.v2^gamma`

`frac{PV}{T} = Constant (Ideal gas equation)`

Where,

p1, v1 & t1 are properties at state point 1

p2, v2 & t2 are properties at state point 2

Process 2: Constant volume of heat addition

The piston remains at TDC, while spark is initiated while external heat is added by combustion of the air-fuel mixture. The temperature and pressure values of the gas increase while volume is constant. Since volume is constant, pressure and temperature increase linearly. State point 3 gives the final position in this process. Volume is equal to clearance volume.

`Heat (IN) =m.Cv.(t3-t2)`

Where,

m = Mass of the gas

Cv = Specific heat at constant volume

Process 3: Isentropic expansion

The piston moves from TDC to BDC, while the air in the engine cylinder undergoes reversible adiabatic expansion. This process is shown between state points 3 to 4. From point 3-4, the volume increases while pressure decreases.

`p3.v3^gamma =p4.v4^gamma`

Where,

P3, v3 & t3 are properties at state point 3

P4, v4 & t4 are properties at state point 4

Process 4: Constant volume of heat rejection

The piston remains at BDC, the hot gasses are flown out from the cylinder while rejecting the heat. The volume is equal to total volume (Clearance volume + Swept volume), and this is kept constant. Since volume is constant, temperature and pressure vary linearly.

`Heat (OUT) =m.Cv.(t4-t1)`

Where,

m = Mass of the gas

Cv = Specific heat at constant volume

To calculate the volume of the cylinder, while piston moves can be given by the following equation.

`V = Vc.(1+frac{1}{2}.(cr-1)[R+1-costheta-(R^2-sin^2theta)^frac{1}{2}])`

Where,

V = Volume of cylinder

Vc = clearance volume

cr = Compression ratio

R = Length of connecting rod/Crankpin radius

θ = Angle between connecting rod and crank pin.

From the volume values and isentropic equation, pressure values can be calculated.

For process 1, the piston moves from BDC to TDC having the value of the angle between crankpin and connecting rod varies from 180 deg to 0 deg. Similarly, for process 3, the angle varies from 0deg to 180 deg.

The thermal efficiency of the engine:

Thermal efficiency is the ratio of work done to the heat supplied.

`eta = 1-frac{(t4-t1)}{(t3-t2)} =1-frac{1}{cr^(gamma-1)`

Program:

Main program:

# Engine kinematics
# by Manikanteshwara reddy 

# Importing modules
import math
import matplotlib.pyplot as plt

# Defining the function for calculating volume and pressure values for isentropic process
def cylinder_volume(bore,stoke,c_rod,cr,start_crank,end_crank):
	
	# crank pin radius
	a = stroke/2
	# ratio of connecting rod and crank pin
	R = c_rod/a

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

	sc = math.radians(start_crank)
	ec = math.radians(end_crank)
	# Vaiable for smoothen the curve
	num_values = 50
	dtheta =(ec-sc)/(num_values-1)

	v = []
	for i in range(0,num_values):
		theta = sc + 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)

		v.append(v_c*(1+term1*(term2-term3)))
	return v

# Engine design parameters
# Bore of engine cylinder (m)
bore = 0.1
# Stroke length of piston (m)
stroke = 0.1
# Length of connecting rod (m)
c_rod = 0.15
# Compression ratio (total volume/clearance volume)
cr = 12
# Swept volume (m3)
v_s = (math.pi/4)*pow(bore,2)*stroke
# Clearance volume (m3)
v_c = v_s/(cr-1)
# gamma = Cp/Cv | for air
gamma = 1.4 

# State point 1 inputs
# Pressure (Pascals)
p1 = 101325
# Initial and combustion temperature (K)
t1 = 500
t3 = 2300

# Volume (m3) piston @ BDC
v1 = v_s+v_c

# By ideal gas law (p*v/t = constant), say c1
c1 = (p1*v1)/t1

# State point 2 (1@2> isentropic compression process) (p*v^gamma = constant)
# Piston @ TDC
v2 = v_c
p2 = p1*pow(v1,gamma)/(pow(v2,gamma))
t2 =(p2*v2)/c1

# Volume values over the isentropic curve | Governing equation, p*v^gamma = constant |
v_compression = cylinder_volume(bore,stroke,c_rod,cr,180,0)
constant = p1*pow(v1,gamma)

p_compression = []
for v in v_compression:
	p_compression.append(constant/pow(v,gamma))

# State point 3 (2@3> constant volume process, heat addition)
v3 = v2
p3 = (c1*t3)/v3

# State point 4 (3@4> isentropic expansion process)
# Piston @ BDC
v4 = v1
p4 = p3*pow(v3,gamma)/pow(v4,gamma)
t4 = (p4*v4)/c1

v_expansion = cylinder_volume(bore,stroke,c_rod,cr,0,180)
constant = p3*pow(v3,gamma)

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

# (4@1> constant volume process, heat rejection)

# Thermal efficiency | n_th = 1-(t4-t1)/(t3-t2) |
n_th = 1-(1/pow(cr,gamma-1))
print('Thermal efficiency =',n_th)

# plots
plt.figure()
# Plotting the state points 1, 2, 3 & 4
plt.plot([v1,v2,v3,v4],[p1,p2,p3,p4],'*')
# plot of isentropic compression
plt.plot(v_compression,p_compression)
# plot of constant volume heat addition
plt.plot([v2,v3],[p2,p3])
# plot of isentropic expansion
plt.plot(v_expansion,p_expansion)
# plot of constant volume heat rejection
plt.plot([v4,v1],[p4,p1])
plt.savefig('Otto cycle PV diagram')
plt.show()

Supplementary program to define the function for the volume of the cylinder with respect to the angle theta:

import math
import matplotlib.pyplot as pyplot
bore = 0.1
stroke = 0.1
c_rod = 0.15
cr = 12
# crank pin radius
a = stroke/2
R = c_rod/a

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

theta = math.radians(180)
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 = v_c*(1+term1*(term2-term3))

After defining the function the returned value V i.e output of the function is an array. So if the function is called in the program, the output will be an array of volume with respect to the angle theta.

Appending the value V in function inside the for loop is to be done to store all the values of volumes. If we append V outside the for loop and return the value V to the function, THe output array of volume will have the last calculated value of volume only.

To convert the angle input from degrees to radians, math.radians has to be used.

All the steps and information are commented in the program for a better explanation.

The above program is successfully run in the python.

PV diagram for the above program with a compression ration of 12 is shown below.

Results:

All four process in an Otto cycle is plotted successfully. Two Isochoric processes are linear while the other two processes are Isentropic following a pattern of curve.

Since the compression ratio is high, the volume is observed to be more compressed from the PV diagram. For example, if we lower the compression ratio to 5, compression of volume from the PV diagram can be seen below.

The thermal efficiency of air standard cycle for an engine having the compression ratio of 12 is found to be 0.629