Solve an Otto cycle through MATLAB

AIM : Write a code and solve an otto cycle and make plots for it and create a PV diagram through Matlab

OBJECTIVE : create a PV diagram of otto cycle and output the thermal efficiency of the engine.

EXPLANATION : 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 processes are described by:

• Process 0–1 (intake stroke)

A mass of air is drawn into piston/cylinder arrangement at constant pressure.

• Process 1–2 (compression stroke)

It 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 (ignition phase)

It `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 (expansion stroke)

It is an adiabatic (isentropic) expansion (power stroke).

• Process 4–1(idealized heat rejection)

It completes the cycle by a constant-volume process in which heat is rejected from the air while the piston is at bottom dead centre.

• Process 1–0 (exhaust stroke)

The mass of air is released to the atmosphere in a constant pressure process.

INPUTS :

• Volume at state one V₁
• Temperature at state one T_1’
• With these two we can calculate pressure at state one P₁
• Compression ratio (CR) with these we can calculate V₂ and P₂ .With these two we can calculate temperature T₂
• Either T₃ or P₃ are assumed to calculate V₃ ,because V₁ and V₂ are same volume addition
• When T₃ and V₃ is known P₃ is calculated
• We will able to known hoe volume changes from one to two. This is depend on engine kinematics
• In this particular problem instead of starting from V₁ . We are going to join from the engine geometric parameters

Engine geometric parameters :

• Bore
• Stroke
• Connecting rod length
• Compression ratio
• We will able to solve the kinematics equation for the piston and we will able to calculate swept volume and clearance volume which information required to calculate volume at each state.

PROGRAM :

clear all
close all
clc
 
 

 %Inputs
gamma = 1.4
t3 = 2300
 
%state variables
p1 = 101325
t1 = 500
 
%Engine geometric parameters
bore = 0.1;
stroke = 0.1;
con_rod = 0.15;
Cr = 12;
%calculating the swept volume and the clearance volume
v_swept = (pi/4)*bore^2*stroke;
v_clearance = v_swept/(cr-1);
 
v1 = v_swept + v_clearance
v2 = V_clearance
 
% State variable at state point 2
%p2v2^gamma = p1v1^gamma
% p2 = p1*(v1/v2)^gamma = p1*cr^gamma
 p2 = p1*cr^gamma;
 
 %p1v1/t1 = p2v2/t2 | t2 = p2*v2*t1/(p1*v1)
 t2 = p2*v2*t1/(p1*v1);
 
 % state variables at state point 3
 v3 = v2
 
 %p3v3/t3 = p2v2/t2 | p3 = p2*t3/t2
 p3 = p2*t3/t2
 
 %state variables at state point 4
 v4 = v1
 
 %p3v3^gamma = p4v4^gamma | p4 = p3(v3/v4)^gamma
 p4 = p3*(v3/v4)^gamma
 
 figure(1)
 hold on 
 plot(v1,p1,'*','color','r')
 plot(v2,p2,'*','color','r')
 plot(v3,p3,'*','color','r')
 plot(v4,p4,'*','color','r')
 figure (2)
 plot([v1 v2 v3 v4 v1],[p1 p2 p3 p4 p1],'color','b')
 xlabel('volume (m^3)')
 ylabel('pressure (pa)')
grid on

ERRORS:

UPDATED PROGRAM :

clear all
close all
clc


%Inputs
gamma = 1.4
t3 = 2300

%state variables
p1 = 1.01325
t1 = 500

%Engine geometric parameters
bore = 0.1;
stroke = 0.1;
con_rod = 0.15;
cr = 7;
%calculating the swept volume and the clearance volume
v_swept = (pi/4)*bore^2*stroke;
v_clearance = v_swept/(cr-1);

v1 = v_swept + v_clearance
v2 = v_clearance

% State variable at state point 2
%p2v2^gamma = p1v1^gamma
% p2 = p1*(v1/v2)^gamma = p1*cr^gamma
 p2 = p1*cr^gamma;
 
 %p1v1/t1 = p2v2/t2 | t2 = p2*v2*t1/(p1*v1)
 t2 = p2*v2*t1/(p1*v1);
 

 
 
 % state variables at state point 3
 v3 = v2
 
 %p3v3/t3 = p2v2/t2 | p3 = p2*t3/t2
 p3 = p2*t3/t2

 
 %state variables at state point 4
 v4 = v1
 
 %p3v3^gamma = p4v4^gamma | p4 = p3(v3/v4)^gamma
 p4 = p3*(v3/v4)^gamma
 
 
 figure(1)
 hold on 
 plot(v1,p1,'*','color','r')
 plot(v2,p2,'*','color','r')
 plot(v3,p3,'*','color','r')
 plot(v4,p4,'*','color','r')
 figure (2)
 plot([v1 v2 v3 v4 v1],[p1 p2 p3 p4 p1],'color','b')
 xlabel('volume (m^3)')
 ylabel('pressure (pa)')
 grid on
 
 
 

OUTPUT :Which is completly incorrect, expansion process represented by wrong path curve.

• In order to determine the correctly we have to make two things

1. We have to know themodynamics relationship dring compression process
2. We need to use piston kinematics and determine how the volume changes from Top dead center(TDC) and Bottom dead center(BDC)

UPDATED PROGRAM 2:

clear all
close all
clc


%Inputs
gamma = 1.4
t3 = 2300

%state variables
p1 = 1.01325
t1 = 500

%Engine geometric parameters
bore = 0.1;
stroke = 0.1;
con_rod = 0.15;
cr = 8;
%calculating the swept volume and the clearance volume
v_swept = (pi/4)*bore^2*stroke;
v_clearance = v_swept/(cr-1);

v1 = v_swept + v_clearance
v2 = v_clearance

% State variable at state point 2
%p2v2^gamma = p1v1^gamma
% p2 = p1*(v1/v2)^gamma = p1*cr^gamma
 p2 = p1*cr^gamma;
 
 %p1v1/t1 = p2v2/t2 | t2 = p2*v2*t1/(p1*v1)
 t2 = p2*v2*t1/(p1*v1);
 

 constant_c = p1*v1^gamma;
 V_compression = engine_kinematics(bore, stroke, con_rod,cr,180,0); 
 
 P_compression = constant_c./V_compression.^gamma;
 
 
 % state variables at state point 3
 v3 = v2
 
 %p3v3/t3 = p2v2/t2 | p3 = p2*t3/t2
 p3 = p2*t3/t2
 
 constant_c = p3*v3^gamma;
 V_expansion = engine_kinematics(bore, stroke,con_rod,cr,0,180);
 
 P_expansion = constant_c./V_expansion.^gamma;
 
 %state variables at state point 4
 v4 = v1
 
 %p3v3^gamma = p4v4^gamma | p4 = p3(v3/v4)^gamma
 p4 = p3*(v3/v4)^gamma
 
 
 figure(1)
 hold on 
 plot(v1,p1,'*','color','r')
 plot(V_compression,P_compression)
 plot(v2,p2,'*','color','r')
 plot([v2 v3],[p2 p3],'color','b')
 plot(v3,p3,'*','color','r')
 plot(V_expansion,P_expansion)
 plot(v4,p4,'*','color','r')
 plot([v4 v1],[p4 p1],'color','b')
 xlabel('volume (m^3)')
 ylabel('pressure (pa)')
 grid on
 
 
 

Function used in this program: Funtion is nothing but peace of code that u can use just by calling its name and passing in right parameters.

function [V] = engine_kinematics(bore, stroke, con_rod, cr, start_crank, end_crank)

a = stroke/2;
R =con_rod/a;

v_s = pi/4*bore^2*stroke;
v_c = v_s/(cr-1);

theta = linspace(start_crank, end_crank,100);
 term1 = 0.5*(cr-1);
 term2 = R+1 - cosd(theta);
 term3 = (R^2 - sind(theta).^2).^0.5;
 
 V = (1 + term1*(term2 - term3)).*v_c;
 
end

OUTPUT: P V diagram

Conclusion:

1. The Thermal effieciency is calculated to be 0.5804 or 58.04%.

2. The PV diagram clearly say that shows changes in isentropic processes are a curve and not straight line.