Plotting of an Otto Cycle graph in MATLAB

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:

  1. Function to compute adiabatic volume to eventually compute pressure in adiabatic processes
  2. Main code to call the above function and to plot the graph

1:

function [v] = ottoengine_kinm(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 = (pi/4)*(bore^2)*(stroke);

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

%crank angle theta
%Replacing 720 degrees below with startcrank and endcrank
theta = 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 - cosd(theta);
t3 = (R^2 - (sind(theta)).^2).^0.5;

v = vclear.*(1 + t1 * (t2-t3));

end

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.

 

2:

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;

%Swept volume (when at BDC)
vswept = (pi/4)*(bore^2)*stroke;

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

%Total volume = v1
v1 = vswept + vclear

%volume at top dead centre = v2
v2 = vclear

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

%Computing t2 using ideal gas law p1v1/t1 = p2v2/t2
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*v1^gamma;
v_comp = ottoengine_kinm(bore, stroke, conrod, cr, 180, 0);
p_comp = c_comp./(v_comp.^gamma);

%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
t3 %printing t3 for ease of reading

c_exp = p3*v3^gamma;
v_exp = ottoengine_kinm(bore, stroke, conrod, cr, 0, 180);
p_exp = c_exp./(v_exp.^gamma);

%State variables at state point 4
%v4 = v total
v4 = v1

%p4 using p3v3^gamma = p4v4^gamma
p4 = p3*(v3/v4)^gamma

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

figure(1)
hold on
plot(v2,p1,'o','color','r')%Zero point representing intake phase at atmospheric pressure
plot([v2 v1], [p1 p1], 'linewidth',2, 'color','b') %line showing intake stroke between point 0 and 1st point
plot(v1,p1,'*','color','r') %1st point
plot(v2,p2,'*','color','r') %2nd point
plot(v_comp, p_comp, 'color', 'k') %adiabatic curve from 1st to 2nd point
plot(v3,p3,'*','color','r') %3rd point
plot([v2 v3], [p2 p3], 'color', 'k') %line between 2nd and 3rd point
plot(v_exp, p_exp,'color', 'k') %adiabatic curve from 3rd point to 4th point
plot(v4,p4,'*','color','r') %4th point
plot([v4 v1], [p4 p1], 'color', 'k') %line between 4th point and 1st point
xlabel('Volume'), ylabel('Pressure')

figure(2)
plot([v1 v2 v3 v4 v1], [p1 p2 p3 p4 p1], 'color', 'k') %figure if adiabatic curves weren't plotted using array
xlabel('Volume'), ylabel('Pressure')

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

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

Incorrect Plot:

Incorrect PlotDirect line is drawn between points for adiabatic compressio and expansion. This is wrong as it doesn't show the right characteristics.

 

Efficiency of the Engine:

v1 =  0.00086239
v2 =  0.000095821
p1 =  101325
p2 =  2196120.29612
t1 =  500
t2 =  1204.1
v3 =  0.000095821
p3 =  4194855
t3 =  2300
v4 =  0.00086239
p4 =  193542.98743
t4 =  955.06

nth =  0.58476
nth1 =  0.58476

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:

Error 1= variables incorrectly typed: vtot by v1 and vtdc by v2

vtot =  0.00086239
vtdc =  0.000095821
p2 =  2196120.29612
error: 'v2' undefined near line 35 column 10
error: called from
    otto1 at line 35 column 4

Error 2= v = vclear * (1 + t1 * (t2-t3)).

>> otto_engine_kinematics_arrays

parse error near line 28 of file D:CoursesSkill-LyncMech EssentialsMATLABWeek 2Air Standard CycleOtto cycleotto_engine_kinematics_arrays.m

  syntax error

>>> v = vclear * (1 + t1 * (t2-t3)).

Usage of "." for array multiplication incorrect. Corrected in function above.

 

Error 2:

Plot command mistake

Changed to v2 v1, p1 p1.


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

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

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

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