MATLAB PROGRAM FOR SIMULATING AIR STANDARD OTTO CYCLE

AIM:-

Ø Write a MATLAB program for simulation of Air Standard Otto cycle.

INTRODUCTION:-

Ø This program is for simulating theoretical limit of particular air standard cycle, by which we can get an idea about performance of real world systems like IC engine & Gas turbines.

Ø This MATLAB program is for simulate the Otto cycle.

Ø The Air Standard Otto cycle is the ideal cycle for four stroke Spark ignition internal combustion engines, first proposed by Nikolause Otto over 130 years ago, which is currently used in most of the motor vehicles.

Ø In analysis of Otto cycle we will use ideal “air standards” assumptions which are:-

1. Gas and air mixture are modeled as air and an ideal gas, which continuously circulates in a closed cycle. Thus, there are no intake and exhaust processes.
2. All the processes making up the cycle are internally reversible.
3. The combustion process is replaced by a heat-addition process from an external source.
4. The exhaust process is replaced by a heat-rejection process and the gas returns to its initial state.

Ø The ideal Otto cycle the fuel-air mixture is introduced during the induction stroke and compressed to a much lower compression ratio (around 8:1) and is then ignited by a spark. The combustion results in a sudden jump in pressure while the volume remains essentially constant. The continuation of the cycle including the expansion and exhaust processes.

GOVERNING EQUATION:-

Ø Otto cycle consists four processes which are

· Isentropic Compression Process (1-2):-Piston moves from bottom dead centre (BDC) to top dead centre (TDC).

Isentropic,

1) P*Vϒ = Constant

P1*V1 ϒ = P2*V2 ϒ

P2 = P1*(V1/V2)^(gamma)P2=P1∗(V1/V2)(ϒ)

P2 = P1*(cr)^(gamma)P2=P1∗(cr)(ϒ)

2) (P1*V1)/(T1) = (P2*V2)/(T2)(P1∗V1)/(T1)=(P2∗V2)/(T2)

(T2) = (P2*V2*T1) / (P1*V1)(T2)=(P2∗V2∗T1)/(P1∗V1)

· Constant Volume Heat addition Process (2-3):- Heat transfer to the working gas from an external source while the piston is at top dead center.

1) V2 = V3

2) (P3*V3)/(T3) = (P2*V2)/(T2)(P3∗V3)/(T3)=(P2∗V2)/(T2)

(P3) = (P2*T3)/(T2)(P3)=(P2∗T3)/(T2)

· Isentropic Expansion Process (3-4):-Piston moves from top dead center (TDC) to bottom dead center (BDC).

Isentropic,

1) V1 = V4

2) P3*V3ϒ = P4*V4 ϒ

P4 = (P3)*(V3/V4)^(gamma)P4=(P3)∗(V3/V4)(ϒ)

· Constant Volume Heat rejection Process (4-1):- Heat rejected from the working gas while the piston is at bottom dead center.

Where,

P1, V1, T1 = Pressure (N/m2), Volume (m3) & Temperature (K) at point 1 (Starting of compression)

P2, V2, T2 = Pressure (N/m2), Volume (m3) & Temperature (K) at point 2 (End of compression)

P3, V3, T3 = Pressure (N/m2), Volume (m3) & Temperature (K) at point 3 (Starting of expansion)

P4, V4, T4 = Pressure (N/m2), Volume (m3) & Temperature (K) at point 4 (End of expansion)

ϒ = ratio of specific heat (Cp/Cv)

cr = compression ratio

Ø Equations Used for calculating Swept volume:-

· Swept Volume:-

It is the volume between top dead center (TDC) and bottom dead center (BDC)

V_swept = (pi/4)*D^2*L

Where,

D = Diameter of cylinder (m)

L = Stroke of the cylinder (m)

· Clearance Volume:-

Clearance volume is a volume between the cylinder head and the piston top when the piston is at top dead center (TDC)

V_clearance = V_swept/ (cr-1)

V1 = V_swept + V_clearance (Total Volume of Cyclinder)

V2 = V_clearance (Volume of cylinder when piston at TDC)

Ø Equation for change in volume at changing crank angle:

V/VcV/VcV/Vc=1+(1/2)∗(cr−1)∗[R+1−cos(θ)−(R2−sin2(θ))]V/Vc = 1 + (1/2)*(cr-1)*[R + 1 - cos(\theta)-\sqrt(R^2-sin^2(\theta))]V/Vc=1+(1/2)∗(cr−1)∗[R+1−cos(θ)−(​R2−sin2(θ))]

Where,

V = Volume during compression and expansion

Vc = Clearance volume

Cr = compression ratio

R = l/al/a

a = Crank pin radius

l = connecting rod length

θ = crank angle.

Ø Thermal Efficiency of Otto cycle:-

(η)(\eta)

(eta) = 1 - (1)/(cr^(gamma-1))(η)=1−(1)/(cr(ϒ−1))

Where,

Cr = Compression ratio

ϒ = ratio of specific heat (Cp/Cv)

OBJECTIVE:-

Ø Program to calculate parameters at all states in of Otto cycle.

Ø Program to plot P-V diagram of Otto cycle after calculating parameters.

PROGRAM:-

1) Program for calculating parameters:-

% A MATLAB program for simulation of Otto cycle 
clear all
close all
clc

% inputs for calculation of other parameters 
gamma = 1.4                       % Specific heat ratio(Cp/Cv)
T3 = 2300                         % Temperature at beginning of expansion (K)

% state variables
P1 = 101325                       % Pressure at beginning of compression (N/m2)
T1 = 500                          % Temperature at beginning of compression (K)

% Engine geometric parameters
D = 0.1;                          % Bore of Cylinder (m)
L = 0.1;                          % Length of the cylinder or stroke (m)
l = 0.15;                         % Connecting rod length (m)
cr = 12 ;                         % compression ratio

% Calculating swept_volume & Clearance_volume
V_swept = (pi/4)*D^2*L ;          % Volume travel by piston during one stroke (m3)
V_clearance = V_swept/(cr-1);     % Volume between Cylinder head & piston head when piston at TDC (m3)

V1 = V_swept + V_clearance        % Total Volume of Cylinder (m3)
V2 = V_clearance                  % Volume of cylinder when piston at TDC (m3)

% State variables for compression process(1 - 2)
% P2V2^gamma = P1V1^gamma
% P2 = P1(V1/V2)^gamma = P1*cr^gamma
P2 = P1*cr^gamma                  % Pressure at end of compression (N/m2)

% P1V1/T1 = P2V2/T2 | T2 = P2*V2*T1/(P1*V1)
T2 = P2*V2*T1/(P1*V1);            % Temperature at end of compression (K)

% Equation for piston kinematics during compression stroke
constant_c = P1*V1^gamma;          
V_compression = engine_kinematics_piston(D,L,l,cr,180,0);
P_compression = constant_c./V_compression.^gamma;

% State variables for Heat addition process (2-3)
V3 = V2                           

% P3*V3/T3 = P2*V2/T2 | P3 = P2*T3/T2
P3 = P2*T3/T2                    % Pressure at beginning of expansion (N/m2)
V_heat_addition = engine_kinematics_piston(D,L,l,cr,0,0);
P_heat_addition = linspace(3285264,5593140,180);

% State variable for expansion process (3-4)
V4 = V1

% P3V3^gamma = P4V4^gamma | P4 = P3*(V3/V4)^gamma
P4 = P3*(V3/V4)^gamma             % Pressure at the end of expansion (N/m2)

% Equation for piston kinematics during expansion stroke
constant_c = P3*V3^gamma;
V_expansion = engine_kinematics_piston(D,L,l,cr,0,180);
P_expansion = constant_c./V_expansion.^gamma;

% State variable for heat rejection process (4-1)
V_heat_rejection = engine_kinematics_piston(D,L,l,cr,180,180);
P_heat_rejection = linspace(101325,172510,180);

% Thermal Efficiency of otto Cycle
E = 1 - 1/(cr^(gamma-1))

% Program for plotting parameters on P-V diagram
figure(1)
hold on
plot(V1,P1,'*','color','r');
plot(V2,P2,'*','color','r');
plot(V_compression,P_compression)
plot(V3,P3,'*','color','r');
plot(V_expansion,P_expansion)
plot(V4,P4,'*','color','r');
plot(V_heat_addition,P_heat_addition);
plot(V_heat_rejection,P_heat_rejection);
title('P-V diagram of otto cycle','fontsize',23)
xlabel('Volume (m3)','fontsize',15)
ylabel('Pressure (N/m2)','fontsize',15)
text(V1,P1, 'State-1','fontsize',10)
text(V2,P2, 'State-2','fontsize',10)
text(V3,P3, 'State-3','fontsize',10)
text(V4,P4, 'State-4','fontsize',10)
text(V_compression(140),P_compression(140),'Compression','fontsize',12)
text(V_expansion(40),P_expansion(40),'Expansion','fontsize',12)
grid on

figure(2)
plot([V1 V2 V3 V4 V1],[P1 P2 P3 P4 P1],'color','r');
title('P-V diagram of otto cycle')
xlabel('Volume(m3)')
ylabel('Pressure(N/m2)')
grid on

Output:-

2) Program for plotting parameters:-

From this program we will get P-v diagram for Otto cycle but, as per plot we can see all point is directly connected this is not the right way to connect the parameters. As we know in isentropic process there should be curve instead straight line. So , this graph is not showing correct results.

In this case to get correct result we plot isentropic process for changing crank angles, for these we required other function program.

figure(2)

plot([V1 V2 V3 V4 V1],[P1 P2 P3 P4 P1],'color','r');

title('P-V diagram of otto cycle')

xlabel('Volume(m3)')

ylabel('Pressure(N/m2)')

grid on

Output:-

3) Program for solve piston kinematics equation:-

This program basically tells us how piston moves from TDC to BDC & how combustion chamber volume changes as function of crank angle.

function [V] = engine_kinematics_piston(D,L,l,cr,start_crank,end_crank)
 
 % Engine geometric parameters
 a = L/2                    % crank pin radius
 R = l/a 
 
 % caclulating swept_volume & clearance_volume
 V_swept = (pi/4)*D^2*L ;          % Volume travel by piston during one stroke
 V_clearance = V_swept/(cr-1);
 
 theta = linspace(start_crank,end_crank,180);
 
 term1 = 0.5*(cr-1);
 term2 = R + 1 - cosd(theta);
 term3 = (R^2 - sind(theta).^2).^0.5;
 
 V = (1 + term1*(term2 - term3)).*V_clearance;
 
 
 end

Output:-

V/Vc = 1 + (1/2)*(cr-1) * [ R + 1 - cos(theta) - sqrt{R^2 - (sin*(theta)^2)}]

By calling this function program and inserting value as per change in crank angle during process we can get piston kinematics.

1) Equation for piston kinematics during compression stroke

constant_c = P1*V1^gamma;

V_compression = engine_kinematics_piston(D,L,l,cr,180,0); ( During compression stroke angle decreases from 180 to 0.)

P_compression = constant_c./V_compression.^gamma;

2) Equation for piston kinematics during heat addition process.

V_heat_addition = engine_kinematics_piston(D,L,l,cr,0,0); (During heat addition Volume remains constant)

P_heat_addition = linspace(3285264,5593140,180); (During heat addition pressure changes from P2 to P3)

3) Equation for piston kinematics during expansion stroke

constant_c = P3*V3^gamma;

V_expansion = engine_kinematics_piston(D,L,l,cr,0,180); ( During expansion stroke angle increases from 180 to 0.)

P_expansion = constant_c./V_expansion.^gamma;

4) Equation for piston kinematics during heat rejection process

V_heat_rejection = engine_kinematics_piston(D,L,l,cr,180,180); (During heat rejection Volume remains constant)

P_heat_rejection = linspace(101325,172510,180); (During heat rejection pressure changes from P3 to P4 )

After Using these equation if we plot P-V diagram it look as follows:-

Errors:-

1) I made a typo, wrote incorrect spelling of clearance.

Output:-

2) I gave function name incorrect, there should be function [V] = … , to avoid this kind of mistakes ensure function name and Equation name should be same.

Output:-

CONCLUSION:-

1. The area bounded by the complete cycle path represents the total work that can be done during one cycle.
2. In first plot we can see area under P-v diagram is much more as compare to second plot which is incorrect. Second plot represent correct total work done during the cycle.
3. By this plot we can get an idea about performance of engine operated on Otto cycle.
4. By this programming we can compare two different cycles for better performance and efficiency.