Python program - Otto cycle with piston kinematics

OBJECTIVE:-

Write a program in python to solve an otto cycle, create a pressure-volume plot and calculating the thermal efficiency of the engine.

THEORY:-

Otto cycle is an air standard cycle which forms as primary fundamental for the operation of an SI engine.

Fig(1) Operational sequence of the otto cycle.

SOURCE:- http://www.learneasy.info/MDME/iTester/tests/10907_Engine_Performance/images/4-stroke.html

The operations in an Otto cycle areas follows:-
Process 0-1 Intake stroke:- The piston moves from TDC (top dead center) to BDC (bottom dead center) while sucking the air-fuel mixture inside the chamber. An increase in volume is observed.

2. Process 1-2 Isentropic compression stroke:- The piston moves from BDC to TDC compressing the fuel-air mixture. An increase in pressure with a proportional decrease in volume is observed.

3. Process 2-3 Constant pressure Heat addition:- The piston is at TDC and the fuel-air mixture is ignited by means of a spark via a spark plug. A rise in pressure is observed which is at its peak at point 3 during the time of ignition.

4. Process 3-4 Isentropic expansion stroke:- The piston moves from top to bottom during the stroke. The pressure decreases along with a proportional increase in the volume.

5. Process 4-1 Constant pressure heat rejection:- In this process the remants of the fuel air mixture are rejected out of the system as the piston rests at BDC for the moment.

ASSUMPTIONS:-

1. Working fluid is approximated as Air (Ideal Gas).

2. The compression and expansion processes are reversible and adiabatic.

3. Heat added during combustion is at constant volume (piston at TDC).

4. Heat is rejected at constant volume when piston at BDC.

5. The piston displacement for intake and exhaust strokes is excluded. (Ref 1)

6. Piston at TDC; crank angle = 0°.

7. Piston at BDC; crank angle = 180°.

8. P0 = 101325 PASCAL (ATMOSPHERIC PRESSURE), V0 = Clearance volume.

MATHEMATICAL FORMULAE USED:-

1. For reversible and adiabatic (isentropic processes) `P * V^gamma = constant`

2. For an ideal gas; `(P*V)/T` = constant

3. For thermal efficiency:

`eta_(th) =`(Work done/ Heat Supplied)

work done = heat supplied - heat rejected. (Heat is supplied from process 2 to 3 and Heat is rejected from process 3 to 4.)

Heat Supplied:- `Q_s = m * C_v * (t3 - t2)`.

Heat Rejected:- `Q_r = m * C_v * (t4 - t1)`.

`eta_(th) = (m * C_v * (t3 - t2) - m * C_v * (t4 - t1))/ (m * C_v *(t3 - t2)`.

On simplifying we get;

`eta_(th) = 1 - ((t4 - t1)/(t3 - t2))`.

SOLUTION:-

A program was implemented in python to generate an ideal Pressure Volume plot for Otto cycle.

Code:-

''' 
Program:- Otto cycle
Author: Harsh Shah
'''

#Import libraries
import math
import matplotlib.pyplot as plt

# Function defintion 
def engine_kinematics(bore, stroke, con_rod, cr, start_crank, end_crank):
    '''
    ENGINE KINEMATICS
    '''
	
    # Geomteric parameters of the engine
    a = stroke/2
    R = con_rod/a

    # Volume parameters
    V_s = math.pi/4*pow(bore,2)*stroke
    V_c = V_s/(cr-1)
	
    # Specifying angles
    sc = math.radians(start_crank)
    ec = math.radians(end_crank)

    # No. of crank angles considered 
    num_values = 30
	
    # Spacing between the crank angles
    dtheta = (ec-sc)/(num_values-1)

    # Storing volume values for corresponding crank angles
    V = []

    for i in range(0,num_values):
	theta = sc + i* dtheta
		
	# Solving the volume trace equation
	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)
	
        # Storing calculated volume values
	V.append((1 + term1*(term2-term3))* V_c)
		
    return V  # returning value calculated by the function 

# Input variables 
p1 = 101325    # PRESSURE FOR STATE 1
t1 = 300       # TEMPERATURE FOR STATE 1
gamma = 1.4    # CONSTANT
t3 =2500       # TEMPERATURE FOR STATE 3

# Geometric parameters
bore = 0.0825
stroke = 0.0928
con_rod = 0.154
cr = 11.5

# Swept Volume and Clearance volume
v_s = (math.pi/4)*pow(bore,2)*stroke
v_c = v_s/(cr -1)
v1 = v_c + v_s # Volume @ 1

# State point 2
v2 = v_c # Volume @ 2

# p2.v2^gamma = p1.v1^gamma
p2 = p1 * pow(v1,gamma)/pow(v2,gamma) #Pressure @ 2

# (p2.v2)/t2 = (p1.v1)/t1 | ideal gas law
rhs = p1*v1/t1
t2 = p2*v2/rhs # Temperature @ 2

# Calculating V_compression, P_compression for certain crank angles
# p1* v1^gamma = constant
k = p1 * pow(v1,gamma)

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

P_compression=[]

for volume in V_compression:
    P_compression.append(k/pow(volume,gamma))

# State Point 3
v3 = v2 # volume @ 3

# (p3.v3)/t3 = (p2.v2)/t2
rhs = p2*v2/t2
p3 = rhs * t3/v3 # pressure @ 3

# Calculating V_expansion, P_expansion for certain crank angles
k = p3 * pow(v3,gamma)

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

P_expansion=[]

for volume in V_expansion:
    P_expansion.append(k/pow(volume,gamma))

# State point 4
v4 = v1 # volume @ point 4

# (p4.v4^gamma)=(p3.v3^gamma)
p4 = p3 * pow(v3,gamma)/pow(v4,gamma) # pressure @ 4

# (p4.v4)/t4 = (p3.v3)/t3 = (p2.v2)/t2 = rhs | Ideal Gas Law
t4 = (p4*v4)/rhs # temperature @ 4

# Generating plots

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.grid()
plt.show()

# Calculating thermal efficiency
term4 = (t4-t1)/(t3-t2)
therm_eff = (1 - term4)
print("Thermal efficiency of the cycle is: {}" .format(therm_eff))

Code analysis:-

1. Importing libraries for plotting and using mathematical commands.

2. Defining a function. A function is set of code which you can access to perform certain tasks which need to be done repetitively within a program.

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

• def is specifies to suggest that a function is being defined.
• engine_kinematics is the name of the function which shall be passed to a variable to calculate and store certain values in the program further.
• (x,y,z) are varying parameters which need to be specified mandatory. The function calculates as per the paramters passed.
• : is to be specified mandatory just like for loop. The staemnt preceding the : needs ot be provided with indentation.

Then parameters for volume and geometric parameters are provided as inputs.

Line 24,25:- While specifying crank angles; radians() command is used to convert the angle values pecified to radians. As the trignometric commands from math library take default input in radians.

Then we calculate trace volumes for varying crank angles and they are stored in an empty array using append command.

Line 47:- return V returns the value of calculated trace volume when the funcftion is called.

Now we specify some Geometric parameters of the engine like bore, stroke, etc and then we specify variables like pressure, volume and temperature at point 1.

Then we proceed to calculate swept and clearance volumes.

We calculate parameters like pressure, volume and temperature using the mathematical formulae specified earlier.

Now after calculating varaibles for points 1 and 2 (i.e. isentropic compression process) on

Line 78:-

k = p1 * pow(v1,gamma)

• As `P* V^gamma =` constant

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

• V_compression is the variable which shall store the values calculated by the function.
• engine_kinematics calls the function to calculate the trace volume in our case as per the arguements passed on by the user.
• () the values are arguments to be passed which shall be used for calculating the trace volume. Parameters like bore,stroke, con_rod, cr are specified at the start.
• The start and end crank angle are passed as 180 and 0 repsectively as the piston moves from BDC to TDC.

P_compression=[]

• It is an empty array to calculate pressure values for the varying volume values during compression stroke.

for volume in V_compression:
P_compression.append(k/pow(volume,gamma))

• The specified for loop stores a value from V_compression array value and uses it to calculate the the P_compression.
• While calculating and appending the value to the array we need to be careful to pass volume as an arguement instead of V_compression to calculate the pressure; as V_compression is an array while volume is a variable storing the value.

Similarly we calculate values for P_expansion with crank angles varying. During expansion piston moves from TDC to BDC and start crank angle = 0 and end crank angle =180.

Then post calculating the state variables and varying pressure values for varying volume values as per the crank angle variation we plot the values.

Line 115, 116:-

plt.plot(V_compression,P_compression)
plt.plot(V_expansion,P_expansion)

This are mentioned seperately as the variation is not linear.

as per the above specified formula.

print("Thermal efficiency of the cycle is: {}" .format(therm_eff))

"" specifies a string to be printed in the output window.

.format(therm_eff) returns the value of calculated thermal efficiency

OUTPUT RESULT:-

The output for the PV plot is as shown below:-

Fig(2) - PV plot

The result for thermal efficiency is as follows:-

Fig(3)- Thermal Efficency output

CONCLUSIONS:-

The PV diagram is plotted and the thermal efficiency of the engine is output. We can conclude that the PV diagram plot is dependent on state variables as well as the geometric dimensions of the engine which may differ and hence make it possible to calculate theoretical limits for the engine.

REFERENCES:-

1. Hershey, R. L., Eberhardt, J. E., and Hottel, H. C. (1936). Thermodynamic properties of the working fluid in internal combustion engines, SAE Journal, 39, 409–419, October.

2. https://en.wikipedia.org/wiki/Lamborghini_V10 (FOR GEOMETRIC DIMENSIONS OF ENGINE).