Air standard cycle graph plotting in Python

Project involves plotting an Otto Cycle graph in MATLAB

The goal is to plot all strokes for a petrol operating by the otto cycle from intake, compression, power to exhaust, including both adiabatic curves.

Description:

The Otto Cycle is the idealised cycle that describes the functioning of a typical spark ignition piston engine.

Processes:

  1. 0-1: Constant Pressure intake stroke for air intake
  2. 1-2: Adiabatic (isentropic) compression of charge as piston moves from BDC to TDC
  3. 2-3: Constant volume heat transfer from external source to gas when piston is at TDC; represents ignition
  4. 3-4: Adiabatic (isentropic) expansion; power stroke
  5. 4-1: Constant volume heat rejection when piston is at BDC
  6. 1-0: exhaust stroke; mass of air released

Code:

The code is split into 2 sections, where one if a function defined inside the code (and minimised) to compute adiabatic volume to eventually compute pressure in adiabatic processes, and the other is the entre code block to plot the graph, which is shown below.

The function above takes bore dia, stroke, connecting rod length, compressionr ratio, and start and end crank angles to compute volume using the formula:

v/vc = 1 + 0.5*(cr-1)[R + 1 - cos(theta) - (R^2 - (sin(theta))^2)^0.5]

This "v" is used to compute volume and pressure change in the adiabaic processes in the function below.

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

# Inputs
gamma = 1.4
t3 = 2300
# State variables- operaing conditions
p1 = 101325
t1 = 500

# Engine geometry; rough dimensions: assumed GM LS9 engine that went into the 2013 Corvette C6 ZR1
bore = 0.103
stroke = 0.092
conrod = 0.15
cr =  9

a = stroke/2 #crank pin radius
R = conrod/a #ratio of conrod length to crank pin radius

# Swept volume (when at BDC)
vswept = (math.pi/4)*pow(bore,2)*stroke

# clearance volume
vclear = vswept/(cr-1)

# Total volume = v1
v1 = vswept + vclear

#volume at top dead centre = v2
v2 = vclear
print('v2 = ', v2)

#State varaiables at state point 2
#p1v1^gamma = p2v2^gamma
#cr = vtot/vtdc
print('p1 = ', p1) #printing p1 for ease of reading
p2 = p1*pow(cr, gamma)
print('p2 = ', p2)

# Computing t2 using ideal gas law p1v1/t1 = p2v2/t2
print('t1 = ', t1) #printing t1 for ease of reading
t2 = (p2*v2*t1)/(p1*v1)

# Adiabatic constants denoted by c_comp and c_exp
c_comp = p1*pow(v1,gamma)
#-------------------------------------------------------
#-------------------------------------------------------
def enginekin(bore, stroke, conrod, cr, startcrank, endcrank):

	a = stroke/2; # crank pin radius
	R = conrod/a; # ratio of conrod length to crank pin radius

	# Swept volume
	vswept = (math.pi/4)*pow(bore,2)*(stroke);

	# clearance volume
	vclear = vswept/(cr-1);

	# crank angle theta
	# Replacing 720 degrees below with 180 and 0
	theta = nm.radians(nm.linspace(startcrank, endcrank, 100))

	# formula => v/vc = 1 + 0.5*(cr-1)[R + 1 - cos(theta) - (R^2 - (sin(theta))^2)^0.5]
	T1 = 0.5*(cr-1);
	T2 = R + 1 - nm.cos(theta);
	# t3 = (pow(R,2) - (nm.sin(theta))^2)^0.5;
	T3 = pow((pow(R,2)-pow(nm.sin(theta),2)), 0.5)

	v = vclear*(1 + T1 * (T2-T3))

	return v
#-------------------------------------------------------
#-------------------------------------------------------

v_comp = enginekin(bore, stroke, conrod, cr, 180, 0)
p_comp = c_comp/pow(v_comp, gamma)

print('p_comp = ', p_comp, 'c_comp', c_comp)

# State variables at point 3; constant volume
v3 = v2

# computing p3 using ideal gas law
# p3v3/t3 = p2v2/t2; v3 = v2
p3 = p2*t3/t2
print('t3 = ', t3) # printing t3 for ease of reading
print('p3 > p2?', p3 > p2)

c_exp = p3*pow(v3, gamma);
v_exp = enginekin(bore, stroke, conrod, cr, 180, 0);
p_exp = c_exp/pow(v_exp,gamma);

print('p_exp = ', p_exp, 'c_exp', c_exp)
# State variables at state point 4
# v4 = v total
v4 = v1

# p4 using p3v3^gamma = p4v4^gamma
p4 = p3*pow((v3/v4), gamma)

#t4 using ideal gas law p4v4/t4 = p3v3/t3
t4 = (p4*v4*t3)/(p3*v3)

print('v4 = ', v4, 'p4 = ', p4, 't4 = ', t4)


plt.figure(1)
plt.plot([v2, v3], [p2, p3])
plt.plot(v_comp, p_comp)
plt.plot(v_exp, p_exp)
plt.plot([v4, v1], [p4, p1])
plt.xlabel('Volume')
plt.ylabel('Pressure')
plt.show()

# Thermal Efficiency
nth = 1 - (1/pow(cr,(gamma-1))) # eff = 1 - 1/compression ratio^gamma
nth1 = 1 - ((t4 - t1)/(t3 - t2)) # eff = 1 - (t4 - t1)/(t3 - t2)

print('Thermal Efficienty = ', nth)

The code above starts with assumption of state variables and constants related to the engine we are working on (here, a GM LS9).

  1. Swept, total and clearance volume are calculated.
  2. Adiabatic equations are used to calculate p2 and ideal gas law for t2.
  3. Step repeated for 3rd to 4th state, and t4 calculated using ideal gas law again.

Figure 1 and 2 are two versions of the plot. Figure 1 indicates the actual plot where the curve of the adiabatic processes are plotted. Figure 2 is a simple straight lines closed quadrilateral graph, which is incorrect since it doesn't correctly show the curve of the adiabatic processes.

The function returns an array of volumes using compression ratio, crank pin radius and ratio of conrod to crank pin. This volume, depending on input variables for start crank and end crank (180 to 0 for compression and 0 to 180 for expansion) gives the respective values of volume, which we use along wih the respective pressure values calculated to make the plot.

Correct Plot:

Correct Plot

Efficiency of the Engine:

Thermal Efficienty [eff = 1 - 1/compression ratio^gamma] =  0.5847563534614941
Thermal Efficienty [eff = 1 - (t4 - t1)/(t3 - t2)] =  0.5847563534614941

Above block shows the output screen. Thermal efficiency is calculated using:

  1. 1 - 1/(compression ratio)^gamma
  2. 1 - (T4 - T1)/(T3 - T2)

Both values result in an efficiency of 0.58476, or 58.48%.

 

Errors:

Error 1: Attempting to define function in a separate code and call; didn't work. Had to define and call in main code


Projects by Arjun Bhat

This project explores the simulation of the motion of a pendulum for the equation of damped vibration of a pendulum. Below is the equation (taken from "Challenges" page): The above equation is solved using ODE function, and then odeint is used to get an ar Read more

Objective: The project involves simulating a 2-joint robotic arm.   Code and Description: import math import matplotlib.pyplot as plt import numpy as nm # Lengths l1 = 1 l2 = 0.5 # Angles ang1 = nm.radians(nm.linspace(0, 90, 19)) ang2 = nm.radians Read more

Flow over bicycle in Python
Arjun Bhat · 2019-11-18 16:09:30

Objective: To calculate the drag force for various drag coefficient and velocity values for a certain geometry.   Introduction: Drag force is the force exerted by the air as the vehicle moves through it. The force exerted is directly proportional to the square Read more

Aim: To parse the NASA Thermodynamc Properties data file, calculate thermodynamic properties and plot the characteristics. Introduction: NASA's thermodynamic properties file consists of multiple species of gases, their respective operating temperature ranges, and Read more

The project aims at optimization of a stalagmite function and finding the global maxima. For this, we will be using the concept of genetic algorithm. Following are the steps used to optimise the function: Code the function given in a separate function code. This is Read more

Curve Fitting using MATLAB
Arjun Bhat · 2019-08-19 20:20:53

Project aims at finding the best fit curve for a given data and finding the error between the original and said curve. The data being used is that of Specific heat (Cp) vs Temeperature. Plot will be made to show variation of Cp with Temp first.  Code: % Data cp Read more

This project explores the simulation of the motion of a pendulum for the equation of damped vibration of a pendulum. Below is the equation (taken from "Challenges" page): An attempt is made to solve the above equation, and then call it using ode45 to get an array of Read more

Project involves plotting an Otto Cycle graph in MATLAB The goal is to plot all strokes for a petrol operating by the otto cycle from intake, compression, power to exhaust, including both adiabatic curves. Description: The Otto Cycle is the idealised cycle that descr Read more

Code and Description: The project involves simulating a 2-joint robotic arm. % Lengths l1 = 1; l2 = 0.5; % Angles ang1 = linspace (0, 90, 20); ang2 = linspace (0, 90, 20); count = 1; % loop for i = 1:length(ang1) ANG1 = ang1(i); for j = 1:length(ang2 Read more


Loading...

The End