Air standard Cycle using Python

Aim: To plot PV diagram for air standard otto cycle and calculate thermal efficiency using Python

Theory:

Thermodynamic cycles are the most fundamantal part of any powerplant or an engine. Different thermodynamic cycles are otto cylce, diesel cycle, dual cycle. rankine cycle etc. Among all these PV diagram of otto cycle is being visualized with Matlab here.

Otto cycle consists of four process: two constant pressure process and two adiabatic process. Here ideal cycle is considered for PV diagram visualization. The Four process are as follows:

1. Process(1-2) -- Adiabatic Compression

2. Process(2-3) -- Constant pressure heat addition

3. Process(3-4) -- Adiabatic Expansion

4. Process(4-1) -- Constant pressure heat rejection

Basic Equations:

For Adiabatic process:

`P*V^gamma` = Constant

Ideal gas law:

`P*V/T` = Constant

Variation of Volume in cylinder with crank angle is given by:

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

where, R = connecting Rod lenght / crank pin radius

crank pin radius = stroke/2

Python Program:

"""
    Air standard otto cycle simulator
"""

import math
import matplotlib.pyplot as plt
import numpy as np

def engine_kinematics(bore, stroke, connecting_rod, cr, start_crank, end_crank):
    """
        Engine Kinematics
    """

    a = stroke/2
    R = connecting_rod/a

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

    num_values = 100

    theta = np.linspace(sc, ec, num_values)

    v = []

    for i in theta:
        term1 = 0.5*(cr-1)
        term2 = R + 1 - math.cos(i)
        term3 = pow(R,2) - pow(math.sin(i),2)
        term3 = pow(term3,0.5)
        v.append((1 + term1*(term2 - term3))*v_c)

    return v
    

# inputs
p1 = 101325
t1 = 298
gamma = 1.4
t3 = 2800

# geometric parameters
bore = 0.1
stroke = 0.1
connecting_rod = 0.15
cr = 12

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

v_compression = engine_kinematics(bore, stroke, connecting_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 2
v2 = v_c
# p1*v1^gamma = p2*v2^gamma
rhs = p1*pow(v1,gamma)
p2 = rhs/pow(v2,gamma)
# p1*v1/t1 = p2*v2/t2
rhs = (p1*v1)/t1
t2 = rhs/(p2*v2)

# state point 3
v3 = v2
# p2*v2/t2 = p3*v3/t3
p3 = rhs*t3/v3

v_expansion = engine_kinematics(bore, stroke, connecting_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
# p3*v3^gamma = p4*v4^gamma
rhs = p3*pow(v3,gamma)
p4 = rhs/pow(v4,gamma)
# p3*v3/t3 = p4*v4/t4
rhs = (p3*v3)/t3
t4 = rhs/(p4*v4)

# plot results
plt.plot(v_compression, p_compression)
plt.plot([v2,v3],[p2,p3])
plt.plot(v_expansion,p_expansion)
plt.plot([v1,v4],[p1,p4])
plt.xlabel('Volume(m^3)')
plt.ylabel('Pressure(Pa)')
plt.title('Air standard otto cycle')
plt.savefig('Air standard otto cycle')

# Thermal Efficiency
therm_eff = 1 - 1/pow(cr,gamma-1)
print('Thermal Efficiency = ' + str(therm_eff))

A function engine_kinematics is defined which returns volume for different crank angle during compression and expansion process. It takes bore, stroke, connecting rod length, starting carnk angle and end crank angle as input and return volume values as an array.

math.radians converts degree angle into radians. np.linspace is used to create equally spaced crank angle values between start_crank and end_crank. A for loop is defined to calculate volume for diiferent crank angle using formulae discussed earlier.

Inputs and geometric parameters are defiend and volume is computated at point 1 using formulae. Volume variation during compression process is determined by calling enigne_kinematics function with crank angle from 180 to 0 deg. Pressure values are determined for compression process using adiabatic relationship using for loop. Temprature at point 2 is determined using ideal gas law.

Pressure and point 3 is defined using ideal gas law. Volume variation during expansion process is determined by calling enigne_kinematics function with crank angle from 0 to 180 deg. Pressure values are determined for expansion process using adiabatic relationship using for loop. Temprature at point 4 is determined using ideal gas law.

Cycle is plotted using pressure and volume values by using plt.plot . Axis labels and chart title are defined and plot is saved as png using plt.savefig . Thermal efficiency of the cylce is calculated and printed using print command.

Output:

Thermal Efficiency = 0.6298928275128466