Python: Calculate state variable and plot P-V diagram of Air Standard Cycle (Otto Cycle)

AIM: To write code in Python to plot P-V diagram for air standard cycle (Otto cycle).

OBJECTIVES:

• Calculate state variable and plot a P-V diagram of air standard cycle (Otto Cycle).
• Calculate output thermal efficiency of the engine.

GOVERNING EQUATION:

Ideal Gas Law: `"PV/T=constant"`

Adiabatic Process: `"(PV)^γ=constant"`

Piston Volume = `V/V_c =1+1/2 (r_c-1)[R+1-cos⁡θ-[(R^2-sin^2 θ)^(1/2) ]`

Thermal Efficiency (`η`) = `1-1/(compressionratio)^(γ-1)`

`"Swept Volume= π/4*(Bore)^2*Stroke"`

`"Clearance Volume=(Swept Volume)/(Compression Ratio-1)"`

THEORY:

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 cycle consists of following four processes -

• Process 1–2 is an adiabatic (isentropic) compression of the charge as the piston moves from bottom dead center (BDC) to top dead center (TDC). Fig.(a)
• Process 2–3 is a constant-volume heat transfer to the working gas from an external source while the piston is at top dead center. This process is intended to represent the ignition of the fuel-air mixture and the subsequent rapid burning. Fig.(b)
• Process 3–4 is an adiabatic (isentropic) expansion (power stroke). Fig.(c)
• 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 center. Fig.(d)

Figure 3.11. Processes in a reciprocating engine (i.e. a diesel or petrol engine). (a) compression (b) combustion (c) expansion (d) exhaust and induction.

STEPS:

First, we import important libraries require for our program. In this program we import math and matplotlib.pyplot (as plt) library.

We need to define a function from which we calculate piston volume which will help us generate PV plot. The full equation is big and complex so we write it in parts (terms) and append them. We calculate the value of piston volume for 100 different angles of theta and store in array named ‘V’ for further calculations.

After successfully defining the function for calculating piston volume we define values for parameters for 4 state conditions. For 1 state condition we assume some values like value of gamma is taken as 1.4 (for air), ambient temperature = 350K and atmospheric pressure = 101325Pa. Also, the values for engine geometric properties are inputted. From this swept volume, clearance volume, v1 and v2.

The maximum temperature inside the cylinder is known. Other state variables for remaining states are calculated using ideal gas equation and adiabatic process equation.

Now that we have calculated all the values, we plot the PV diagram but the plot obtained is incorrect. The compression and expansion stroke are represented by straight lines which is incorrect. To solve this problem a different function is created for storing the different values of crank angle rotation as compression and expansion occurring in different strokes or cycles. This function is then called to the main program and accordingly the plot is drawn from various points obtained. For this number of values for which the program will iterate is 100, and accordingly, the V array is set as null and after every iteration using for loop, a new value will be stored as an array using 'append' command and will be returned back after receiving desired values. Now the plot is generated and proper x and y axis limit, labels and plot title is provided.

Finally, we calculate the thermal efficiency and print it as an output.

CODE:

''' A program for Air Standard Cycle (Otto Cycle)
 By Utkarsh Palwekar'''

import math
import matplotlib.pyplot as plt

#Defining Function
def engine_kinematics(bore, stroke, connecting_rod, c_r, start_crank_angle, end_crank_angle):
    a = stroke / 2
    R = connecting_rod / a
    v_swept = (math.pi / 4) * pow(bore, 2) * stroke
    v_clearance = v_swept / (c_r - 1)

    theta_startcrank = math.radians(start_crank_angle)
    theta_endcrank = math.radians(end_crank_angle)

    num_val = 100
    dtheta = (theta_endcrank - theta_startcrank) / (num_val - 1)
    V = [] #Piston Volume
    for i in range(0, num_val):
        theta = theta_startcrank + i * dtheta
        term_1 = 0.5 * (c_r - 1)
        term_2 = R + 1 - math.cos(theta)
        term_3 = (pow(R, 2) - pow(math.sin(theta), 2))
        term_3 = pow(term_3, 0.5)
        V.append((1 + term_1 * (term_2 - term_3)) * v_clearance)

    return V


# Inputs varibale

gamma = 1.4  # for air gamma= 1.4
T3 = 2000  # max temp inside engine cyc

# State variable
# ambient atmospheric conditions assumed
p1 = 101325
T1 = 350

# Engine geometry properties

bore = 0.15
stroke = 0.2
connecting_rod = 0.8
c_r = 10

# Calculating swept volume and clearance volume

v_swept = (math.pi / 4) * pow(bore, 2) * stroke
v_clearance = v_swept / (c_r - 1)
v1 = v_swept + v_clearance
v2 = v_clearance

# state variables at state point 2
# p2*v2^gamma = p1*v1^gamma --- Isentropic compression

p2 = p1 * pow(c_r, gamma)

# (p1*v1)/T1 = (p2*v2)/T2 | T2 = p2*v2*T1/(p1*v1)

T2 = p2 * v2 * T1 / (p1 * v1)

V_comp = engine_kinematics(bore, stroke, connecting_rod, c_r, 180, 0)
C_const = p1 * pow(v1, gamma)
P_comp = []
for v in V_comp:
    P_comp.append(C_const / pow(v, gamma))

# State variables at state point 3

v3 = v2  # Isochoric process

''' (p3*v3)/T3 = (p2*v2)/T2 | p3 = p2*T3/(T2) '''

p3 = p2 * T3 / (T2)

V_expp = engine_kinematics(bore, stroke, connecting_rod, c_r, 0, 180)
C_const = p3 * pow(v3, gamma)
P_expp = []
for v in V_expp:
    P_expp.append(C_const / pow(v, gamma))

# state variables at state point 4
v4 = v1  # Isochoric process

''' p3*v3^gamma = p4*v4^gamma '''  # Isentropic expansion process
p4 = p3 * pow(v3, gamma) / pow(v4, gamma)

''' (p3*v3)/T3 = (p4*v4)/T4 | T4 = p4*v4*T3/(p3*v3) '''
T4 = p4 * v4 * T3 / (p3 * v3)

#Plotting PV Diagram
plt.plot(V_comp, P_comp)
plt.plot([v2, v3], [p2, p3])
plt.plot(V_expp, P_expp)
plt.plot([v4, v1], [p4, p1])
plt.ylim([-1000000, 8000000])
plt.xlim([0, 0.005])
plt.xlabel('Volume (m^3)')
plt.ylabel('Pressure(Pa)')
plt.title('PV diagram for Air Standard-Otto Cycle')
plt.text(0.00385, -250000, r'$  1  $')
plt.text(0.0002, 2000000, r'$  2  $')
plt.text(0.0002, 5850000, r'$  3  $')
plt.text(0.00385, 400000, r'$  4  $')
plt.grid()
plt.show()

# Calculating thermal efficiency

thermal_efficincy = 1 - (1 / pow((c_r), (gamma - 1)))
print('Thermal efficiency =', thermal_efficincy)

RESULTS/OUTPUT:

PV diagram for Air Standard Cycle (Otto Cycle)

Thermal Efficiency of the Air Standard Cycle (Otto Cycle)

CONCLUSION:

The PV diagram for air standard cycle (Otto cycle) is plotted successfully and efficiency is calculated.

ERRORS & SOLUTIONS:

The compression and expansion stroke were represented by straight lines which was incorrect. To solve this problem a different function was created for storing the different values of crank angle rotation as compression and expansion occurring in different strokes or cycles. This function was then called to the main program and accordingly the plot was drawn from various points obtained.