Air standard Cycle - MATLAB

Objective : Write MATLAB Code for creating P-V diagram and Calculate the thermal efficiency of an Air standard cycle (Otto Cycle) of SI Internal Combustion engine.

Theory :

The air standard otto cycle is ideal cycle dor the spark-ignition (SI) internal combustion engines.

The otto cycle 1-2-3-4 cosists of four processes:

1. Process 1-2 : Isentropic compression of Air
2. Process 2-3 : constant volume Heat adition
3. Process 3-4 : Isenrtopic expansion
4. Process 4-1 : costant volume Heat rejection

The Cycle is shown below as P-V and T-s diagrams.

Thermal efficiency of Air standard cycle = Net Work done / Net Heat added

Geometric Inputs : Bore, Stroke, Connecting rod length, starting and ending crank angle

Volume as function of Theta:

Here

a = stroke/2

R = con_rod/a

MATLAB CODE :

clear all
close all
clc

%otto cycle
%inputs
p1 = 101325;                   % Pressure at stage 1 in Pascal
t1 = 560;                      % Temprature at stage 1 in Kelvin
t3 = 2400;                     % Temprature at stage 3 in Kelvin
gamma = 1.2;                   % value of gamma 


%engine parameter
bore = 0.1;                    % Bore diameter in meter
stroke = 0.1;                  % length of stroke in meter
con_rod = 0.15;                % length of connecting rod in meter
CR = 12;                       % compression ratio

v_swept = (pi/4)*bore^2*stroke;   % swept volume
v_clearance = v_swept/(CR-1);     % clearance volume

v1 = v_swept + v_clearance;       % volume of piston at stage 1
v2 = v_clearance;                 % volume of piston at stage 1

p2 = p1*(CR)^gamma;               % pressure at stage 2

t2 = (p2*v2*t1)/(p1*v1);          % Temperature at stage 2

constant = p1*v1^gamma;           % using Isentropic process formula
v_compression = otto_cycle_pv(0.1,0.1,0.15,12,180,0); % calling function for calculating compression volume as function of theta
p_compression = constant./v_compression.^gamma;       % compression pressure

p3 = p2*(t3/t2);                  % pressure at state 3 using gas law

v3 = v2;                          % Isobaric process
v4 = v1;                          % Isobaric process
p4 = p3*(v2/v1)^gamma;            % pressure at stage 4

constant = p3*v3^gamma;           % Isentropic process
v_expansion = otto_cycle_pv(0.1,0.1,0.15,12,180,0); % calling function for calculating expansion volume as function of theta
p_expansion = constant./v_expansion.^gamma;         % expansion pressure

t4 = (p4*v4*t3)/(p3*v3);          % Temperature at stage 4

%thermal efficiency
T_E = (1-(t4-t1)/(t3-t2))*100 ;   % Calculating Thermal efficiency
fprintf('Thermal efficiency is %f',T_E)    % Displaying Thermal efficiency
%plotting          

figure(1)    %Plotting PV Diagram of Otto cycle.
hold on
plot(v1,p1,'*','color','r')
plot(v_compression,p_compression)
plot(v2,p2,'*','color','r')
plot(v3,p3,'*','color','r')
plot(v4,p4,'*','color','r')
plot(v_expansion,p_expansion)
plot([v2 v3],[p2 p3],[v4 v1],[p4 p1],'linewidth',1)

xlabel('volume in meter')
ylabel('Pressure in pascal')
title('PV diagram for otto cycle')

% Function for calculating Volume of piston as function of theta
function [v] = otto_cycle_pv(bore,stroke,con_rod,CR,start_crank,end_crank)

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

v_s = (pi/4)*bore^2*stroke;  % calculating swept volume
v_c = (v_s)/(CR-1);          % v_c is clearance volume and CR is compression ratio

theta = linspace(start_crank,end_crank,100); % crank angle
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;   % volume of piston as function of theta

end

Code Expanation: P1, t1, t3, bore, stroke, con_rod, CR and gamma are inputs.Then v_swept and v_clearance is calculated using appropriate formulae.Other unknown parameters calculated using given inputs and Thermodynamics processes.A function is created for calculating the compression and expansion volume as function of theta.Then Thermal efficiency is calculated using all Temperature and plot of P-V diagram is created.

Result :

Thermal efficiency is 39.163566