Otto cycle simulation using Engine Kinematics in python

AIM:
To generate the PV diagram and to determine the thermal efficiency of the engine.

THEORY:

Otto Cycle is an idealized thermodynamic cycle that describes the working of a spark-ignition engine. The cycle consists of 4 processes as illustrated in the figure below:

it consists of 4 strokes namely:

1.Suction Stroke(1-2): The charge is compressed isentropically to a smaller volume and hence the pressure and temperature increase to p2,t2 from p1,t1.

2.Compression Stroke(2-3): In this stage, the heat is added at constant volume hence the pressure and temperature increases to p3,t3

3.Expansion Stroke(3-4): The charge is expanded isentropically hence the pressure and temperature reduces to p4,t4

4. Exhaust stroke(4-1): The heat is removed keeping the volume constant, hence reduces the pressure and temperature to p1,t1.

Engine kinematics is used to plot changing values of volume as the piston is moving from Top Dead Centre(TDC) to Bottom Dead Centre(BDC). Hence it is a function of bore, stroke, length of connecting rod and crank angle

Various geometric parameters of the engine are used to calculate Clearance volume(V_c), Swept volume(V_s) and other coordinates of properties they are:

1.Bore(bore):Diameter of cylinder

2.Stroke(stroke): Length that the piston travels from TDC to BDC

3.Length of connecting rod(con_rod)

4.Compression ratio(cr): Degree to which the fuel mixture is compressed

5.Crankpin radius(a)=stroke/2

6.A constant, R=con_rod/a

The isentropic curve is plotted by finding out the volume trace which is given by :

`V/V_c=1+1/2(r_c-1)[R+1-costheta-(R^2-sin^2 theta)^(1/2)]`

PROGRAM:

"""
otto cycle simulator
"""
import math
import matplotlib.pyplot as plt

def engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank):
	"""
		Engine Kinematics
	"""
	
	#Geometric parameters

	
	a=stroke/2
	R=con_rod/a

	#volume parameter
	V_s=math.pi*(1/4)*pow(bore,2)*stroke
	V_c=V_s/(cr-1)

	sc=math.radians(start_crank)
	ec=math.radians(end_crank)

	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((1+term1*(term2-term3))*V_c)
	return V
		




#inputs
p1=101325
t1=500
gamma=1.4
t3=2300

#geometric parameters
bore=0.1
stroke=0.1
con_rod=0.15
cr=8

#Volume computation
v_s=(math.pi/4)*pow(bore,2)*stroke
v_c=v_s/(cr-1)
v1=v_c+v_s


#state point 2
v2=v_c

#p2v2^gamma=p1v1^gamma
p2=p1*pow(v1,gamma)/pow(v2,gamma)

#p2v2/t2=p1v1/t1| rhs= p1v1/t1|p2v2/t2=rhs|t2=p2v2/rhs
rhs=p1*v1/t1
t2=p2*v2/rhs
V_compression=engine_kinematics(bore,stroke,con_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
v3=v2

#p3v3/t3=p2v2/t2|rhs=p2v2/t2|p3=rhs*t3/v3
rhs=p2*v2/t2
p3=rhs*t3/v3

V_expansion=engine_kinematics(bore,stroke,con_rod,cr,0,180)

constant=p3*pow(v3,gamma)

P_expansion=[]

for v in V_expansion: 
	P_expansion.append(constant/pow(v,gamma))


#state point 4
v4=v1

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

#p4v4/t4=rhs
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(cubic metre)')
plt.ylabel('Pressure(Pa)')
plt.show()

#thermal efficiency
Thermal_efficiency=(1-(1/(pow(cr,gamma-1))))*100
print(f'{Thermal_efficiency}%')

PROCEDURE:

1.The math,matplotlib.pyplot modules were imported.

2. The engine_kinematics function returns the volume(V) as an array of volumes which gradually changes during the isentropic processes otherwise the isentropic curve would be linear which is not possible.

3. Initial conditions were defined as

p1=101325 Pa

t1=500K

gamma=1.4(Air)

t3=2300K

4. Geometric parameters required for calculations are defined and volumes are computed.

5. Values of parameters at all states are computed using formulas:

For constant volume process:`(PV)/(T)=constant`

For isentropic processes:`(PV^gamma)=constant`

6. Using plot function the pressures vs volume graph is plotted

7.Thermal efficiency =[1-(1/cr^(gamma-1))]*100`[1-(1/(cr)^(gamma-1))]*100`

OUTPUT:

P V PLOT

The PV plot is generated as shown below

Thermal Efficiency

Thermal efficiency computed is 56.47%

CONCLUSION:

PV diagram was generated and the thermal efficiency of the engine was determined using python.