Simulation of Otto Cycle with P-V curve plotting using Python

OBJECTIVE

The key objective of this program is to make OTTO cycle simulator.
Calculate the unknown pressure & temperature and plot the PV diagram using thermal equations and find the thermal efficiency with the help of Python scripting.

Otto cycle: The air-standard-Otto cycle is the idealized cycle for the spark-ignition internal combustion engines.

An Otto cycle consists of four processes:

1. Two isentropic (reversible adiabatic) processes
2. Two isochoric (constant volume) processes

PROCEDURE

1. Consider some given data:

p1 = 101325 pa
T1 = 500 k
T3 = 2300 k
gamma = 1.4
Stroke length = 0.1 m
Bore dia. = 0.1 m
connecting rod = 0.15 m
compression ratio (cr = 12)
Crankpin radius = stroke length/2
R = constant = con_rod/crankpin radius

2. Some Formulas used:

swept vol. (v_s) = π/4 d²L

Clearance vol. (v_c) = swept vol./(cr -1)

V = v_c + v_s

p1v1/t1 = p2v2/t2 | Ideal gas equation

V/Vc= 1 + 1 (rc - 1)[R + 1 - cos θ - (R² – sin² θ)1/2]

This eq. used to describe a relation between crank angle & vol. inside the combustion chamber. By this eq. we get volume trace and we need this to finish our program.

3. Import the modules whichever required.

• import matplotlib.pyplot as plt
  This module helps us to plot the graph of the given result.
• import math
  This module helps us to use mathematical functions like pie, sine, cos, etc.
• import numpy
  This module helps us to do mathematical operations on an element of an array.
  Whenever we want to define more than one value in a single variable we have to use an array and numpy.linspace is also the same thing it stores value in a range like this; numpy.linspace (lower limit, upper limit, no. of division)
  This works acc. to a formula:
  Size of division = (upper limit-lower limit)/(no. of division-1)
  

Use an empty array for storing a value that we get from FOR loop.

Append is a command that will continuously append the value of a variable. It keeps on adding the values into an existing array.

4. Define the function: We can call the function whenever and wherever we required.

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

sc = start crank angle: when piston at B.D.C.

ec = End crank angle: when piston at T.D.C.

Here return v is used to return the value by ending the execution of function call to the function.

Volume doesn\'t vary linearly in Isentropic processes. so, that should be curved, not a straight line. That\'s why we use V_compression, P_compression, V_expansion & P_expansion with function to give more accurate results to plot the graph.

5. Use plt.plot() command to plot the graph of the given result and use plt.figure() command. For proper animation we use .plt.xlim & plt.ylim command.

6. For calculating Thermal efficiency, eth = (1-((t4-t1)/(t3-t2)))*100

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
	\"\"\"

	#geometic parameters
	a = stroke/2  # a= crank pin radius = stroke length/2
	R = con_rod/a # R= constant = con_rod/a

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

	num_values = 25
	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	

#input
p1 = 101325    #in pascal = 1 atm
t1 = 500       #in kelvin
t3 = 2300
gamma= 1.4

#geometic parameters
bore = 0.1 # dia in m
stroke = 0.1 # length in m
con_rod = 0.15 # length of connecting rod in m
cr = 12  # compression ratio

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

#state pt. 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 pt. 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 pt. 4
v4 = v1

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

#thermal efficiency
eth = (1-((t4-t1)/(t3-t2)))*100
print(eth)

plt.plot([v3,v2],[p3,p2])
plt.plot([v1,v4], [p1,p4])
plt.plot(V_compression, P_compression)
plt.plot(V_expansion,P_expansion)
plt.show()

This graph we get from this program

CONCLUSION

• Mainly we learn about the Functions. With the help of this, we can get the result at any crank angle.
• P-V curve successfully plotted with the help of python scripting.
• The thermal efficiency of the engine we find by scripting.
• Learned the proper use of the FOR loop while using Functions.