Air Standard Cycle Using MATLAB

Aim :

To write the code for PV diagram of the otto cycle and to calculate the thermal efficiency using MATLAB

Theory :

The otto cycle is a description of what happens to a mass of gas as it is subjected to changes of pressure, temperature,volume,addition of heat and removal of heat. It is the thermodynamic cycle most commonly found in automobile engines.

Otto cycle:

• process (0-1) air is drawn into pistion/cylinder at constant pressure
• process (1-2) isentropic compression
• process (2-3) constant volume heat addition
• process (3-4) isentropic expansion
• process (4-1) constant volume heat rejection

Code for PV diagram of otto cycle :

clear all
close all
clc

% input

p1=101325
t1=500
gamma=1.4;
t3=2300;

% engine parameters

bore=.1;
stroke=.1;
con_rod=.15;
cr=12;

% state 1

v_swept=(pi/4)*bore^2*stroke;  % formula for swept volume

v_clr=v_swept/(cr-1);          % formula for clearance volume


v1=v_swept+v_clr;
v2=v_clr;


% state 2
% p1*v1^gamma=p2*v2^gamma (isotropic compression process)
% but v1/v2=cr so p2=p1*cr^gamma, also p*v/t is constant


p2=p1*cr^gamma;

% to compute temp at state 2 p1*v1/t1=p2*v2/t2 (ideal gas law)

% so t2=p2*v2*t1/(p1*v1)

t2=p2*v2*t1/(p1*v1);

constant=p1*v1^gamma;
v_compression= engine_k(bore,stroke,con_rod,cr,180,0);

p_compression=constant./v_compression.^gamma


%state 3
%v3=v2 becasue its the end of combustion process

v3=v2

%p2*v2/t1=p3*v3/t3 ,so p3=p2*t3/t2
p3=p2*t3/t2

constant=p3*v3^gamma;
v_expansion= engine_k(bore,stroke,con_rod,cr,0,180);

p_expansion=constant./v_expansion.^gamma;

%state 4 
%v4=v1 ,p3*v3^gamma=p4*v4^gamma , so p4=p3*v3/v4^gamma

v4=v1

p4=p3*(v3/v4)^gamma

%to find temp at state 4 p3*v3/t3=p4v4/t4
%so t4=p4*v4*t3/(p3*v3)

t4=p4*v4*t3/(p3*v3);

%plotting

figure(1)
hold on 
plot(v1,p1,'*','color','r') 
plot(v_compression,p_compression)
plot(v2,p2,'*','color','r')
plot(v3,p3,'*','color','r')
plot(v_expansion,p_expansion)
plot(v4,p4,'*','color','r')
plot([v1,v4],[p1,p4])
plot([v2,v3],[p2,p3])
xlabel('volume')
ylabel('pressure')

%thermal efficiency 

eff =(1-(1/(cr)^(gamma-1)))
percentage_eff = eff*100

Function code :

function [v]= engine_k(bore,stroke,con_rod,cr,start_crank,end_crank)

 a=stroke/2;
 r=con_rod/a;
 
  v_swept=pi/4*stroke^2*bore;
  v_clr=v_swept/(cr-1);
 
  theta=linspace(start_crank,end_crank,70);
 
   term1=0.5*(cr-1);
   term2=r+1-cosd(theta);
   term3=[r^2-sind(theta).^2].^0.5;

   v=(1+term1*(term2-term3)).*v_clr;
end

PV Diagram:

Thermal efficiency output :

eff =

    0.6299


percentage_eff =

   62.9893