Air Standard Cycle (Otto cycle) in PYTHON

AIM : Program to write a code for the PV diagram of an Air Standard Cycle (Otto Cycle) and calculate the thermal efficiency of the engine.

THEORY : An Otto cycle is an idealized thermodynamic cycle that describes the functioning of a typical spark ignition piston engine. It is the thermodynamic cycle most commonly found in automobile engines.The Otto cycle is a description of what happens to a mass of gas as it is subjected to changes of pressure, temperature, volume, addition of heat, and removal of heat.

The processes are described by :

• Process 0–1 a mass of air is drawn into piston/cylinder arrangement at constant pressure.
• Process 1–2 is an adiabatic (isentropic) compression of the charge as the piston moves from bottom dead centre (BDC) to top dead centre (TDC).
• Process 2–3 is a constant-volume heat transfer to the working gas from an external source while the piston is at top dead centre.
• Process 3–4 is an adiabatic (isentropic) expansion (power stroke).
• Process 4–1 completes the cycle by a constant-volume process in which heat is rejected
• Process 1–0 the mass of air is released to the atmosphere in a constant pressure process.

To caluclate all the state variables we use the following formulae :

p2v2^gamma = p1v1^gamma which gives us p2=p1*(v1/v2)^gamma = p1*cr^gamma

we know p1v1/t1=p2v2/t2 so we get t2 = p2*v2*t1/(p1*v1)

we know p3v3/t3=p2v2/tv so we get p3=p2*t3/t2

we know p3v3^gamma=p4v4^gamma so we get p4 = p3*(v3/v4)^gamma

we know p3v3/t3 = p4v4/t4 so we get t4 = p4*v4*t3/(p3*v3)

thermal_efficiency = 1-(1/(cr-1)^(gamma-1))

PYTHON CODE :

""" PROGRAM TO WRITE A CODE FOR PV DIAGRAM OF OTTO CYCLE & TO CALCULATE THE THERMAL EFFICIENCY"""

import matplotlib.pyplot as plt 
import math

def engine_k(stroke,bore,cr,con_rod,start_crank,end_crank):       # function to calculate the required parameters
	a = stroke/2
	R = con_rod/a
	swept_vol = math.pi/4*(pow(bore,2))*stroke                    # to calculate the swept volume  
	clearance_vol = swept_vol/(cr-1)                              # to calculate the clearance volume
	sc = math.radians(start_crank)
	ec = math.radians(end_crank)
	np = 50
	dtheta = (ec-sc)/(np-1)
	v=[]

	for i in range(0,np):
		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((1+term1*(term2-term3))*clearance_vol)
	return(v)

# INPUTS 

gamma = 1.4
p1 = 101325                                                       # state variables
t1 = 500                                                          # state variables
t3 = 2300


# ENGINE PARAMETERS 

bore = 0.001
stroke = 0.04
con_rod = 0.1
cr = 10
swept_vol = (math.pi/4)*pow(bore,2)*stroke
clearance_vol = swept_vol/(cr-1)
v1 = clearance_vol + swept_vol

# AT STATE POINT 2

v2 = clearance_vol
p2 = (p1*pow(v1,gamma))/(pow(v2,gamma))                           # p2v2^gamma = p1v1^gamma which gives us p2=p1*(v1/v2)^gamma = p1*cr^gamma
t2 = (p2*v2*t1)/(p1*v1)                                           # we know p1v1/t1=p2v2/t2 so we get t2 = p2*v2*t1/(p1*v1)
t2 = (p2*v2*t1)/(p1*v1)                 
const = p1*pow(v1,gamma)
v_comp = engine_k(stroke,bore,cr,con_rod,0,180)
p_comp = []
for i in v_comp:
	p_comp.append(const/(pow(i,gamma)))

# AT STATE POINT 3 

v3 = v2
p3 = (p2*t3)/t2                                                   # we know p3v3/t3=p2v2/tv so we get p3=p2*t3/t2

# AT STATE POINT 4 

v4 = v1
p4 = (p3*pow(v3,gamma))/(pow(v4,gamma))                           # we know p3v3^gamma=p4v4^gamma so we get p4 = p3*(v3/v4)^gamma
t4 = (p4*v4*t3)/(p3*v3)                                           # we know p3v3/t3 = p4v4/t4 so we get t4 = p4*v4*t3/(p3*v3)
const1 = p3*pow(v3,gamma)
v_exp = engine_k(stroke,bore,cr,con_rod,180,0)
p_exp = []
for v in v_exp:
	p_exp.append(const1/(pow(v,gamma)))

thermal_efficiency = 1-(1/(pow((cr-1),(gamma-1))))                # calculating thermal efficiency 
print('THERMAL EFFICIENCY = %s'%thermal_efficiency)

# GRAPH PLOTS 

plt.plot([v2,v3],[p2,p3])
plt.plot([v4,v1],[p4,p1])
plt.plot(v_comp,p_comp)
plt.plot(v_exp,p_exp)
plt.xlabel('volume')
plt.ylabel('pressure')
plt.show()

GRAPH :

OUTPUT :

ERRORS :

CORRECTION :thermal_efficiency = 1-(1/(pow((cr-1),(gamma-1))))