PYTHON PROGRAM TO BUILD AN OTTO CYCLE AND PLOT THE P-V DIAGRAM

PYTHON PROGRAM TO BUILD AN OTTO CYCLE AND PLOT THE P-V DIAGRAM

Objective:

To write a Python program to build an Otto cycle and plot its P-V diagram.Also to calculate its thermal efficiency.

Physics involved:

Otto cycle:

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 of isentropic compression, heat addition at constant volume, isentropic expansion, and rejection of heat at constant volume. In the case of a four-stroke Otto cycle, technically there are two additional processes: one for the exhaust of waste heat and combustion products at constant pressure (isobaric), and one for the intake of cool oxygen-rich air also at constant pressure.

• 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 centre(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.

Procedure:

1. Libraries used:

The program requires a plotting function called pyplot that is imported using the library, matplotlib.The alias name used for this function is plt.

Here, a graph of pressure versus volume across all the processes is plotted.

Also math function is used in order to apply trigonometric functions used in the program.

2. Variable declaration and initialization:

The following inputs have been declared and initialized respectively.

P1 = 100 ( supplied pressure in kPa)

T1 = 308 ( supplied temperature in Kelvin)

gamma = 1.4 ( specific heat ratio of air)

T3 = 3632 ( peak temperature of the cycle)

Geometric parameters are

bore = 0.1 (piston diameter in meters)

stroke = 0.1 (distance moved by the cylinder in meters)

con_rod = 0.15 (length of connecting rod in meters)

Cr = 4 (compression ratio)

a = stroke / 2 (half stroke length)

R = con_rod / a (crank pin radius)

Formula used in computing volume parameters are

Vs = (22/7) * (bore)^2 * stroke (Swept volume -Volume moved by the cylinder inside the piston)

Vc = Vs/(Cr - 1) (Clearance volume -Volume not used by the cylinder)

V1 = Vc + Vs (Volume at state point 1)

V2 = Vc (Volume at state point 2)

3. Calculation:

A new function called engine kinematics is defined in order to call the function values to be used in case of isentropic compression and expansion processes.

Variables like SC and EC are initialized in order to declare the starting and ending crank angles respectively .
The cycle is set to run with number of initial values equal to 50.

The crank angle spacing is calculated using the formula

dtheta = (EC - SC) / (no. of values -1)

A for loop is used to find the values of theta and volume with respect to number of initial values.

Relationship between crank angle and volume present inside the combustion chamber is given by the following equation,

Volume = Vc * [1 + 0.5 * (Cr - 1) * (R + 1 - cos(theta) - [(R²) - .sin²(theta)])]

An empty array V stores all the values computed inside the for loop.

All the above variables are defined inside the function, engine kinematics.

Process 1-2 (Isentropic compression process)

Pressure and temperature at state point 2 are calculated using the ideal gas law PV = RT, where R is the universal gas constant.

P2 = (P1 * V1^gamma) / (V2^gamma)

T2 = ((P2 * V2) * T1) / (P1 * V1)

Process 2-3 (Constant volume heat addition process)

Pressure and volume at state point 3 are calculated using the following formula,

P3 = ((P2 * V2) * T3) / (T2 * V3)

V3 = V2

Process 3-4 (Isentropic expansion process)

Pressure at state point 4 is calculated using the following formula,

P4 = P3 *(V3^gamma) / V4^gamma

Process 4-1 (Constant volume heat addition process)

Volume at state point 4 is calculated using the following formula,

V4 = V1

During isentropic compression process, the engine kinematics function is called with starting crank angle input as 180 deg and ending crank angle input as 0 deg and it stored in an array called V_comp.

The corresponding pressure values are calculated using the isentropic relation and stored in an array called P_comp.

Similarly, during isentropic expansion process, the engine kinematics function is called with starting crank angle input as 0 deg and ending crank angle input as 180 deg and it stored in an array called V_exp.

The corresponding pressure values are calculated using the isentropic relation and stored in an array called P_exp.

4.Output:

The efficiency of Otto cycle is calculated using the following formula,

Efficiency = 1 - [1/(Cr^gamma-1)] * 100

The thermal efficiency of Otto cycle is found to be 42.56 %.

A graph is plotted between the Volume and Pressure across all the processes of Otto cycle.

The x-axis is labeled with volume in kg/m³ and y-axis is labeled with pressure in kPa.

Program:

# To build an Otto cycle

import matplotlib.pyplot as plt
import math

def engine_kinematics(bore, stroke, con_rod, Cr, start_crank, end_crank):
    a = stroke / 2
    R = con_rod / a  #crank pin radius

    Vs = (math.pi / 4) * pow(bore, 2) * stroke  # swept volume
    Vc = Vs / (Cr - 1)  # clearance volume

    # crank angle
    sc = math.radians(start_crank)
    ec = math.radians(end_crank)

    n_values = 50

    # crank angle spacing
    dtheta = (ec - sc) / (n_values - 1)
    V = []

    for i in range(0, n_values):
        theta = start_crank + (i * dtheta)

        # Relationship between crank angle and volume present inside the combustion chamber
        Volume = Vc * (1 + 0.5 * (Cr - 1) * (R + 1 - math.cos(theta) - pow((pow(R, 2) - pow(math.sin(theta), 2)), 0.5)))
        V.append(Volume)

    return V

# inputs
P1 = 100      # supplied pressure in kPa
T1 = 308      # supplied temperature in Kelvin
gamma = 1.4   # specfic heat ratio of air
T3 = 3632     # peak temperature

# Geometric parameters
bore = 0.1
stroke = 0.1
con_rod = 0.15
Cr = 4

# volume computation
Vs = (math.pi / 4) * pow(bore, 2) * stroke  # swept volume
Vc = Vs / (Cr - 1)  # clearance volume
V1 = Vc + Vs  # state point 1
V2 = Vc  # state point 2

# isentropic compression process between state points 1 and 2
P2 = (P1 * pow(V1, gamma)) / pow(V2, gamma)
T2 = ((P2 * V2) * T1) / (P1 * V1)

V_comp = engine_kinematics(bore, stroke, con_rod, Cr, 0, 180)
P_comp = []

for i in V_comp:
    P_comp.append((P1 * pow(V1, gamma))/pow(i, gamma))
 
# constant volume process between state points 2 and 3
V3 = V2  # state point 3

P3 = ((P2 * V2) * T3) / (T2 * V3)

V_exp = engine_kinematics(bore, stroke, con_rod, Cr, 0, 180)
P_exp = []

for i in V_exp:
    P_exp.append((P3 * pow(V3, gamma))/pow(i, gamma))

# constant volume process between state points 4 and 1
V4 = V1  # state point 4
P4 = P3 * pow(V3, gamma) / pow(V4, gamma)
T4 = ((P4 * V4) * T2) / (P2 * V2)

Efficiency = 1 - (1/(pow(Cr,gamma-1)))
print('Efficiency=', Efficiency*100)

#plotting
plt.plot(V_comp, P_comp)
plt.plot([V2, V3], [P2, P3])
plt.plot(V_exp, P_exp)
plt.plot([V4, V1], [P4, P1])
plt.xlabel('Volume (kg/m^3)')
plt.ylabel('Pressure (kPa)')
plt.title('P-V Diagram of Otto cycle')
plt.savefig('P-V Diagram of Otto cycle')
plt.show()