Air Standard Cycle Simulation using PYTHON
Simulation of Air Standard Cycle (Otto cycle) using PYTHON
Objective - To simulate Air Standard cycle by creating PV diagram and to calculate output thermal efficiency
Theory -
Processes of Air Standard Cycle-
1) Process 1-2 : Isentropic Compression
2) Process 2-3 : Constant Volume Heat Addition
3) Process 3-4 : Isentropic Expansion
4) Process 4-1 : Constant Volume Heat Rejection
Ideal Gas Equation -
PV`/`T = constant
Equation for Isentropic Compression and Expansion Process -
PV`^``gamma` = constant
Engine Geometric Parameters -
Note : All dimensions are in meter
bore of the cylinder, b = 0.1
stroke of the cylinder, s = 0.1
connecting rod length, l = 0.2
compression ratio, CR = 12
swept volume, v_swept = `pi``/`4`*`b`^`2`*`s
clearance volume, v_clearance = v_swept`/`(CR - 1)
crank radius, r = s`/`2
ratio of connecting rod length to crank radius, R = l`/`r
Equation for calculation of volume during Isentropic Compression and Expansion Process -
v = [1 + 1`/`2`*`(CR - 1)[R + 1 - cos`theta` - `sqrt`(R`^`2 - (sin`theta`)`^`2)]]`*`v_clearance
PYTHON Program -
Code for plotting of P-V diagram and calculating Output Thermal Efficiency:-
# Program to plot P-V diagram and calculate the thermal efficiency of Air Standard Cycle
# Importing the basic libraries/module
import math
import matplotlib.pyplot as plt
import numpy as np
# Inputs for Simulation of Otto Cycle
print(\'Inputs for Simulation of Otto Cycle :-\')
print(\'\\n\')
# A) Engine geometric parameters
print(\'A) Engine Geometric Parameters :-\')
print(\'\\t\')
# Bore of the cylinder (meter)
b = 0.1
print(\'Bore of the cylinder (meter), b = \', b)
print(\'\\t\')
# Stroke of the cylinder (meter)
s = 0.1
print(\'Stroke of the cylinder (meter), s = \', s)
print(\'\\t\')
# Connecting rod length (meter)
l = 0.2
print(\'Connecting rod length (meter), l = \', l)
print(\'\\t\')
# Compression ratio
CR = 12
print(\'Compression ratio, CR = \', CR)
print(\'\\t\')
# Swept volume (m^3)
v_swept = (math.pi/4)*(pow(b, 2))*s
print(\'Swept Volume (m^3), v_swept = \', v_swept)
print(\'\\t\')
# Clearance volume (m^3)
v_clearance = v_swept/(CR - 1)
print(\'Clearance Volume (m^3), v_clearance = \', v_clearance)
print(\'\\n\')
# B)Properties of Air
# Polytropic Index of Air
gamma = 1.4
print(\'B) Properties of Air :-\')
print(\'\\t\')
print(\'Polytropic Index of Air, gamma = \', gamma)
print(\'\\n\')
# Engine Kinematics Function for finding all values of volume during the isentropic compression and expansion process
def engine_kinematics(s, l, crank_start_angle, crank_end_angle, CR, v_clearance):
# Crank radius (meter)
r = s/2
# Ratio of connecting rod length to crank radius
R = l/r
# Creation of an array \'theta\' for calculating various values of cylinder volume at various values of theta
theta = np.linspace(math.radians(crank_start_angle), math.radians(crank_end_angle), 100)
# Number of elements in array \'theta\'
n = len(theta)
# Creating an empty array for storing different values of cylinder volume for corresponding values of theta
volume = []
# For loop for calculating values of cylinder volume using the theta values
for i in range(0, n):
THETA = theta[i]
term1 = 0.5*(CR - 1)
term2 = R + 1 - (math.cos(THETA))
term3 = pow(((pow(R, 2) - pow((math.sin(THETA)), 2))), 0.5)
v = (1 + (term1*(term2 - term3)))*v_clearance
volume.append(v)
return(volume)
# State and Processes
print(\'State and Processes :-\')
print(\'\\t\')
# State 1 Variables
print(\'State 1 Variables :-\')
print(\'\\t\')
v1 = v_swept + v_clearance
t1 = 300
p1 = 101325
print(\'Volume at State 1 (m^3), v1 =\', v1)
print(\'\\t\')
print(\'Temperature at State 1 (K), t1 =\', t1)
print(\'\\t\')
print(\'Pressure at State 1 (N/m^2), p1 =\', p1)
print(\'\\n\')
# Isentropic Compression Process (1 -> 2)
print(\'Isentropic Compression Process (1 -> 2) Variables :-\')
print(\'\\t\')
v_compression = engine_kinematics(s, l, 180, 0, CR, v_clearance)
p_compression = p1*(pow((v1/(np.array(v_compression))), gamma))
print(\'Volume during Compression Process (m^3), v_compression =\', v_compression)
print(\'\\t\')
print(\'Pressure during Compression Process (N/m^2), p_compression =\', p_compression)
print(\'\\n\')
# State 2 Variables
print(\'State 2 Variables :-\')
print(\'\\t\')
v2 = v_clearance
t2 = t1*(pow(CR, (gamma - 1)))
p2 = p1*(pow(CR, gamma))
print(\'Volume at State 2 (m^3), v2 =\', v2)
print(\'\\t\')
print(\'Temperature at State 2 (K), t2 =\', t2)
print(\'\\t\')
print(\'Pressure at State 2 (N/m^2), p2 =\', p2)
print(\'\\n\')
# State 3 Variables
print(\'State 3 Variables :-\')
print(\'\\t\')
v3 = v2
t3 = 1800
p3 = p2*(t3/t2)
print(\'Volume at State 3 (m^3), v3 =\', v3)
print(\'\\t\')
print(\'Temperature at State 3 (K), t3 =\', t3)
print(\'\\t\')
print(\'Pressure at State 3 (N/m^2), p3 =\', p3)
print(\'\\n\')
# Isentropic Expansion Process (3 -> 4)
print(\'Isentropic Expansion Process (3 -> 4) Variables :-\')
print(\'\\t\')
v_expansion = engine_kinematics(s, l, 0, 180, CR, v_clearance)
p_expansion = p3*(pow((v3/(np.array(v_expansion))), gamma))
print(\'Volume during Expansion Process (m^3), v_expansion =\', v_expansion)
print(\'\\t\')
print(\'Pressure during Expansion Process (N/m^2), p_expansion =\', p_expansion)
print(\'\\n\')
# State 4 Variables
print(\'State 4 Variables :-\')
print(\'\\t\')
v4 = v1
t4 = t3*(pow((1/CR), (gamma - 1)))
p4 = p3*(pow((1/CR), gamma))
print(\'Volume at State 4 (m^3), v4 =\', v4)
print(\'\\t\')
print(\'Temperature at State 4 (K), t4 =\', t4)
print(\'\\t\')
print(\'Pressure at State 4 (N/m^2), p4 =\', p4)
print(\'\\n\')
# Output Thermal Efficiency
thermal_efficiency = (1 - (1/(pow(CR, (gamma - 1)))))*100
print(\'Thermal Efficiency = \' + str(thermal_efficiency) + \'%\')
# Plotting
# Creation of empty figure window
plt.figure()
# Plotting of State Variables
plt.plot(v1, p1, marker = \'o\', markerfacecolor = \'w\', markeredgecolor = \'k\')
plt.plot(v2, p2, marker = \'o\', markerfacecolor = \'w\', markeredgecolor = \'k\')
plt.plot(v3, p3, marker = \'o\', markerfacecolor = \'w\', markeredgecolor = \'k\')
plt.plot(v4, p4, marker = \'o\', markerfacecolor = \'w\', markeredgecolor = \'k\')
# Plotting of Processes
plt.plot(v_compression, p_compression, color = \'r\', linewidth = 2)
plt.plot([v2, v3], [p2, p3], color = \'y\', linewidth = 2)
plt.plot(v_expansion, p_expansion, color = \'b\', linewidth = 2)
plt.plot([v4, v1], [p4, p1], color = \'m\', linewidth = 2)
# Giving title to the P-V diagram
plt.title(\'P-V Diagram for Otto Cycle\', fontsize = 20, fontweight = \'bold\')
# Labelling the x and y axes
plt.xlabel(\'Volume (m^3)\', fontsize = 12, fontweight = \'bold\')
plt.ylabel(\'Pressure (N/m^2)\', fontsize = 12, fontweight = \'bold\')
# Inserting the required text 1, 2, 3, 4 in the figure
plt.text(0.000857313, -127545, \'1\', fontweight = \'bold\')
plt.text(0.0000637116, 3073080, \'2\', fontweight = \'bold\')
plt.text(0.0000637116, 7350000, \'3\', fontweight = \'bold\')
plt.text(0.000857313, 279000, \'4\', fontweight = \'bold\')
# Inserting the text Thermal Efficiency = 62.98928275128466% in the figure window
plt.text(0.000331141, 4748000, \'Thermal Efficiency = \' + str(thermal_efficiency) + \'%\', color = \'g\', fontweight = \'bold\', fontsize = 15)
# Saving the figure window
plt.savefig(\'P-V Diagram for Otto Cycle.png\')
# Making the plot visible
plt.show()
Explanation of Code :-
1) First Step - Importing the required basic modules
Description:-
import math % used to import the \'math\' module if in case any equation in the code contains any trigonometric function or/and any mathematical constant (like pi)
import matplotlib.pyplot as plt % used to import the \'plotting\' module which is used to obtain any plots as per the requirement
import numpy as np % used to import the \'numpy\' module which is used to perform any mathematical operations on each element of an array or between the array and also to use the \'linspace\' command for generating an array
Output of the step :- All modules get imported.
2)Sceond Step - Defining inputs and and getting them print
Description :-
b= 0.1 % Bore (diameter) of the cylinder (meters)
s= 0.1 % Stroke of the cylinder (meters)
l = 0.2 % Length of connecting rod (meters)
CR = 12 % Compression ratio
v_swept % Swept volume of the cylinder defined using geometric relations (m^3)
v_clearance % Clerance volume of the cylinder defined using geometric relations (m^3)
gamma = 1.4 % Polytropic index of the air
Output of the step :- All input parameters i.e engine geometric parameters(bore, stroke, connecting rod length, compression ratio) and properties of air(polytropic index) gets defined. Engine geometric parameters like swept volume and clearance volume aslo gets calculated and defined. Also, all the inputs get printed in the command window using the \'print\' command.
3)Third Step - Defining function named as engine_kinematics further used for calculating the volume trace during isentropic compression and isentropic expansion process
Description :-
def engine_kinematics(s, l, crank_start_angle, crank_end_angle, CR, v_clearance): % Syntax for defining the function
engine_kinematics %Name of the function
s %Stroke of the cylinder in meter passed as an input argument in the function
l %Connecting rod length in meter passed as an input argument in the function
crank_start_angle %Starting value of the crank angle for the process like 180 degrees and 0 degrees for isentropic compression and an isentropic expansion process respectively passed as an input argument in the function (in degrees)
crank_end_angle %End value of the crank angle for the process like 0 degrees and 180 degrees for isentropic compression and an isentropic expansion process respectively passed as an input argument in the function (in degrees)
CR %Compression ratio for the engine passed as an input argument in the function
v_clearance %Clearance volume passed as an input argument in the function in m^3
r = s/2 %Crank radius in meter
R = l/r %Ratio of connecting rod length to crank radius
theta = np.linspace(math.radians(crank_start_angle), math.radians(crank_end_angle), 100) % An array containing different values of crank angle during the compression and expansion process created using the the command \'linspace\' and the \'numpy\' module
math.radians(crank_start_angle) %Command used to convert the value of crank_start_angle from degrees into radians
math.radians(crank_end_angle) %Command used to convert the value of crank_start_angle from degrees into radians
100 %Required number of elements
n = len(theta) %Number of elements in an array \'theta\' calculated using the command \'len\' and passing \'theta\' as an input argument
volume = [] %An empty array created for storing different values of cylinder volume for all values of theta during both isentropic compression and expansion processes
for i in range(0, n): %Syntax for defining for loop for calculating values of cylinder volume using the theta values having counter variable \'i\' varying from 0 to n
THETA = theta[i] %Variable used for storing the value of theta fetched during every iteration of for loop
v %Volume trace for every value of theta during each iteration
volume.append(v) %Command used to append an empty array \'volume\' (defined earlier) with value \'v\'
return(volume) %Command used to return the value of an array \'volume\'
4)Fourth Step - State 1 variables are defined
Description:-
v1 = v_swept + v_clearance %Volume (meter^3) at state 1 which is calculated as sum of swept and clearance volume i.e total volume as at state 1 the angle of crank is 180 deg which represents that piston is at the bottom dead centre
t1 = 300 %Temperature (degree kelvin) at state 1 which is defined by considering air to be at normal atmospheric condition
p1 = 101325 %Pressure (N/m^2) at state 1 which is defined by considering air intake to be at normal atmospheric condition
Output of the step:- Volume at state 1 gets calculated. Temperature and pressure at state 1 gets defined. Also, all the variables get printed using the \'print\' command.
5)Fifth Step - Isentropic compression process (1 -> 2) variables are calculated by calling the function named engine_kinematics
Description:-
v_compression %An array containing volume in m^3 during compression process (crank angle varies from 180 to 0 degree) calculated by calling the function \'engine_kinematics\'
p_compression %An array containing pressure values in N/m^2 for the corresponding values of volume during compression is calculated using values like v1, p1, gamma and array \'v_compression\' and commands like \'np.array()\' and \'pow\'
Output of the step :- Volume at hundred different values of angle of crank and corresponding pressure during compression process both gets calculated and presented in a row vector form. While calculating the volume, the crank angle varies from 180 degree to 0 degree. Also, all the variables get printed using the \'print\' command.
6)Sixth Step - State 2 variables are defined
Description:-
v2 = v_clearance %Volume (m^3) at state 2 which is calculated as equal to clearance volume as at state 2 the angle of crank is 0 deg which represents that piston is at the top dead centre
t2 = t1*(pow(CR, (gamma - 1))) %Temperature (degree kelvin) at state 2 is calculated from the thermodynamic relations using values like t1, compression ratio and polytropic index and command \'pow\'
p2 = p1*(pow(CR, gamma)) %Pressure (N/m^2) at state 2 is calculated from the thermodynamic relations using values like p1, compression ratio and polytropic index and command \'pow\'
Output of the step:- Volume at state 2 gets defined. Temperature and pressure at state 2 gets calculated. Also, all the variables get printed using the \'print\' command.
7)Seventh Step - State 3 variables are defined
Description:-
v3 = v2 %Volume (m^3) at state 3 which is defined as equal to volume at state 2 as process 2 -> 3 is isochoric (i.e constant volume) process
t3 = 1800 %Temperature (degree kelvin) at state 3 which is defined as 1800 degree kelvin
p3 = p3 = p2*(t3/t2) %Pressure (N/m^2) at state 3 which is calculated from the thermodynamic relations using values like p2, t2 and t3
Output of the Step :- Pressure at state 3 gets calculated. Temperature and volume at state 1 gets defined. Also, all the variables get printed using the \'print\' command.
8)Eigth Step - Isentropic expansion process (3 -> 4) variables are calculated by calling the function named engine_kinematics
Description:-
v_expansion %An array containing volume in m^3 during expansion process (crank angle varies from 0 to 180 degree) calculated by calling the function \'engine_kinematics\'
p_expansion %An array containing pressure values in N/m^2 for the corresponding values of volume during compression is calculated using values like v3, p3, gamma and array \'v_expansion\' and commands like \'np.array()\' and \'pow\'
Output of the step :- Volume at hundred different values of angle of crank and corresponding pressure during expansion process both gets calculated and presented in a row vector form. While calculating the volume, the crank angle varies from 0 degree to 180 degree. Also, all the variables get printed using the \'print\' command.
9)Ninth Step - State 4 variables are defined
Description:-
v4 = v1 %Volume (meter^3) at state 4 which is defined as equal to volume at state 1 as process 4 -> 1 is isochoric (i.e constant volume) process
t4 = t3*(pow((1/CR), (gamma - 1))) %Temperature (degree kelvin) at state 4 is calculated from the thermodynamic relations using values like t3, compression ratio and polytropic index and command \'pow\'
p4 = p3*(pow((1/CR), gamma)) %Pressure (N/m^2) at state 4 is calculated from the thermodynamic relations using values like p3, compression ratio and polytropic index and command \'pow\'
Output of the step:- Volume at state 4 gets defined. Temperature and pressure at state 4 gets calculated. Also, all the variables get printed using the \'print\' command.
10)Tenth Step - Calculation of output thermal efficiency and getting it print in the command window
Description:-
thermal_efficiency = (1 - (1/(pow(CR, (gamma - 1)))))*100 %Output thermal efficiency is calculated using the formula of thermal efficiency of the air standard cycle
print(\'Thermal Efficiency = \' + str(thermal_efficiency) + \'%\') %It is printed in the command window using \'print\' command and \'str\' command for converting numeric value to string data type and by passing thermal_efficiency as an input argument
Output of the step :- Output thermal efficiency gets calculated and printed in the command window.
11)Eleventh Step - Part 1 of the process of plotting of P-V diagram i.e creation aof an empty figure window and plotting of state variables and processes
Description:-
plt.figure() %Command used to create an empty figure window
plt.plot(v1, p1, marker = \'o\', markerfacecolor = \'w\', markeredgecolor = \'k\') %Command used to plot the state 1 variable
plt.plot(v2, p2, marker = \'o\', markerfacecolor = \'w\', markeredgecolor = \'k\') %Command used to plot the state 2 variable
plt.plot(v3, p3, marker = \'o\', markerfacecolor = \'w\', markeredgecolor = \'k\') %Command used to plot the state 3 variable
plt.plot(v4, p4, marker = \'o\', markerfacecolor = \'w\', markeredgecolor = \'k\') %Command used to plot the state 4 variable
plt.plot(v_compression, p_compression, color = \'r\', linewidth = 2) %Command used to plot an isentropic compression (1 -> 2) process
plt.plot([v2, v3], [p2, p3], color = \'y\', linewidth = 2) %Command used to plot an isochoric (2 -> 3) process
plt.plot(v_expansion, p_expansion, color = \'b\', linewidth = 2) %Command used to plot an isentropic expansion (3 -> 4) process
plt.plot([v4, v1], [p4, p1], color = \'m\', linewidth = 2) %Command used to plot an isochoric (4 -> 1) process
Output of the step :- An empty figure window gets created. All state variables and processes get plotted in the figure window.
12)Twelfth Step - Part 2 of the process of plotting of P-V diagram i.e giving title to the P-V diagram and labelling the x and y axes
Description:-
plt.title(\'P-V Diagram for Otto Cycle\', fontsize = 20, fontweight = \'bold\') %Command used to give title to the plot
plt.xlabel(\'Volume (m^3)\', fontsize = 12, fontweight = \'bold\') %Command used to give label to the x-axis
plt.ylabel(\'Pressure (N/m^2)\', fontsize = 12, fontweight = \'bold\') %Command used to give label to the y-axis
Output of the step :- Plot gets titled using the \'plt.title\' command. X and Y axes gets labelled using the \'plt.xlabel\' and \'plt.ylabel\' commands respectively.
13)Thirteenth Step - Part 3 of the process of plotting of P-V diagram i.e inserting the required text 1, 2, 3, 4 and the text Thermal Efficiency = 62.98928275128466% in the figure window
Description:-
plt.text(0.000857313, -127545, \'1\', fontweight = \'bold\') %Command used to insert text \'1\' with \'bold\' feature by passing the x and y coordinates and \" fontweight = \'bold\' \" as an input arguments
plt.text(0.0000637116, 3073080, \'2\', fontweight = \'bold\') %Command used to insert text \'2\' with \'bold\' feature by passing the x and y coordinates and \" fontweight = \'bold\' \" as an input arguments
plt.text(0.0000637116, 7350000, \'3\', fontweight = \'bold\') %Command used to insert text \'3\' with \'bold\' feature by passing the x and y coordinates and \" fontweight = \'bold\' \" as an input arguments
plt.text(0.000857313, 279000, \'4\', fontweight = \'bold\') %Command used to insert text \'4\' with \'bold\' feature by passing the x and y coordinates and \" fontweight = \'bold\' \" as an input arguments
plt.text(0.000331141, 4748000, \'Thermal Efficiency = \' + str(thermal_efficiency) + \'%\', color = \'g\', fontweight = \'bold\', fontsize = 15) %Command used to insert the text \" Thermal Efficiency = 62.98928275128466% \" and numeric data type is converted into string data type using the command \'str\'
Output of the step :- The required text gets inserted in the figure window.
14)Fourteenth Step - Saving the figure window and making the plots visible
Description:-
plt.savefig(\'P-V Diagram for Otto Cycle.png\') %Command used for saving the plot in PNG format
plt.show() %Command used for making the figure window visible once it is loaded and created
Output of the step :- The figure window gets saved in PNG format and the plot becomes visible.
Errors Encountered:-
1. Error 1 - Syntax error while printing the state variables occurred due to the use of an invalid syntax i.e not giving comma after single quotes while using the \'print\' command
Step taken - Use of correct syntax i.e giving comma after single quotes while using the \'print\' command. After this all state variables were printed in the command window and no error occurred.
2. Error 2 - Due to not defining an empty array \'volume\' before the for loop, a name error occurred
Step taken - An empty array is defined and created before the for loop, which further gets appended in the for loop. Array \'volume\' gets printed and no error occurred thereafter.
Output of the program:-
Variables (i.e volume, pressure and temperature) at all states and plot of PV diagram both are obtained. Thermal efficiency of the air standard cycle is also obtained.
Plot of PV diagram:-
Screenshot of Command Window :-
Output Thermal Efficiency:-
Thermal Efficiency = 62.98928275128466%