Simulation of OTTO CYCLE and calculation of its Thermal Efficiency Using Python

Code to Simulate the PV diagram of an OTTO CYCLE, and to calculate the THERMAL EFFICIENCY with the results obtained.

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 consists 4 Processes, which can be seen in the diagram as given below:
Please Note: This figure is taken from wikipedia, however we will write a code to generate the same.

Figure: Otto Cycle PV Diagram

Processes:

• 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 Center (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. This process is intended to represent the ignition of the fuel-air mixture and the subsequent rapid burning.
• 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 from the air while the piston is at bottom dead centre.
• Process 1–0 the mass of air is released to the atmosphere in a constant pressure process.

Processes 1-2, and 3-4 which are Adiabatic Compression and Adiabatic Expansion if must be simulated we must define a formula that relates the piston motion and the volume of the fuel in the Cylinder, this formula is as given below.

Where:
rc = Compression Ratio
theta is the angles between the start and end crank
Vc = Clearance Volume
V = Volume of Fuel in the Cylinder

CODING:

Engine Parameters Assumed:
p1 = 101325 #Pressure1 in PA#
t1 = 500 #Temperature in Kelvin#
gamma = 1.4 #Ratio of Specific Heats
t3 =2300 #Temperature in Kelvin#
bore = 0.1 #Bore Diameter in meters#
stroke = 0.1 #Length of stroke in meters#
con_rod = 0.15 #Length of Connecting Rod in meters#
CR = 3 #Compression Ratio#

The Code for the function which simulates the adiabatic process, and returns the values of Volumes is as given below:

# -*- coding: utf-8 -*-
\"\"\"
Created on Sun Jul  1 00:07:13 2018

@author: shabarish
\"\"\"
#Defining Engine Kinematics for Isentropic Processes#
import math
def Engine_Kinematics(bore, stroke, con_rod, CR, start_crank, end_crank):
    #Geometric Parameters#    
    a = stroke/2 #crank pin radius#
    R = con_rod/a
    #Volume Parameters#
    V_S = math.pi*(1/4)*pow(bore,2)*stroke
    V_C = V_S/(CR-1)
    #Defining starting and ending Crank#
    sc = math.radians(start_crank)
    ec = math.radians(end_crank)
    num_values = 50
    dtheta = (ec - sc)/(num_values - 1)
    V = []  #Valume trace which will be defined later#
    #We need to solve an Equation given below to define the isentropic process#    
    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

The Main Program which plots the PV diagram for the OTTO Cycle is as given below:

# -*- coding: utf-8 -*-
\"\"\"Created on Fri Jun 22 13:33:28 2018
@author: shabarish
\"\"\"
#Defining Engine Kinematics for Isentropic Processes#
import math
import Engine_Kinematics as EK
#Simulation of Otto Cycle#
import matplotlib.pyplot as plt
#Inputs#
p1 = 101325             #Pressure1 in PA#
t1 = 500                #Temperature in Kelvin#
gamma = 1.4             #Ratio of Specific Heats
t3 =2300                #Temperature in Kelvin#
#Geametric Parameters#
bore = 0.1              #Bore Diameter in meters#
stroke = 0.1            #Length of stroke in meters#
con_rod = 0.15          #Length of Connecting Rod in meters#
CR = 3                 #Compression Ratio#
#Volume Coputation#
V_S = ((math.pi)/4)*pow(bore,2)*stroke #Swept Volume in m^3#
V_C = V_S/(CR-1)
V1 = V_S + V_C          #State point 1 variables V1,p1,t1 have been calculated#
#State Point 2 Variables to calculate#
V2 = V_C
# p1*V1^gamma = p2*V2^gamma #
p2 = p1*pow(V1,gamma)/pow(V2, gamma)
# Use Ideal Gas Law to find t2 #
# (p2*V2)/t2 = (p1*V1)/t1 #
rhs = (p1*V1)/t1
t2 = (p2*V2)/rhs
constant = p1 * pow(V1, gamma)
V_compression = EK.Engine_Kinematics(bore, stroke, con_rod, CR, 180, 0)
p_compression = []
for v1 in V_compression:
    p_compression.append(constant/ pow(v1, gamma))
# To calculate state point 3 Variables #
# Constant Volume heat addition #
V3 = V2
# Ideal Gas Law to find P3 #
# (p2*V2)/t2 = (p3*V3)/t3 #
rhs = (p2*V2)/t2
p3 = (rhs*t3)/V3
constant2 = p3 * pow(V3, gamma)
V_expansion = EK.Engine_Kinematics(bore, stroke, con_rod, CR, 0, 180)
p_expansion = []
for v2 in V_expansion:
    p_expansion.append(constant2/ pow(v2, gamma))
# To calculate state point 4 Variables #
V4 = V1
# p4*V4^gamma = p3*V3^gamma #
p4 = p3*pow(V3,gamma)/pow(V4,gamma)
# Ideal Gas Law to find t4 #
# (p3*V3)/t3 = (p4*V4)/t4 #
rhs = (p3*V3)/t3
t4 = (p4*V4)/rhs
#Plotting#
plt.plot([V2,V3],[p2,p3])
plt.plot(V_compression, p_compression)
plt.xlabel(\'Volume [m^3]\')
plt.ylabel(\'Pressure [Kg/m^2]\')
plt.plot(V_expansion, p_expansion)
plt.plot([V4,V1],[p4,p1])
plt.show()
a = (t4 - t1)/(t3 - t2)
TEff = 1 - a
print(\'Thermal Efficiency = \', TEff)

Using this Code, when we simulate the OTTO Cycle we will get the PV Diagram as shown Below:

Figure: PV Diagram of OTTO Cycle.

The thermal efficiency of the engine, for the given parameters as Outputted by the Program is: 0.35560598502274576