MATLAB PROGRAM FOR SIMULATION OF AIR STANDARD DIESEL CYCLE

AIM:-

• Write a MATLAB program for simulation of Air Standard Diesel 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 Diesel cycle.
• The Diesel cycleis a combustion process of a reciprocating internal combustion engine. First proposed by Rudolph Diesel over 100 years ago. In it, fuel is ignited by heat generated during the compression of air in the combustion chamber, into which fuel is then injected. This is in contrast to igniting the fuel-air mixture with a spark plug as in the Otto cycle (four-stroke/petrol) engine. Diesel engines are used in aircraft, automobiles, power generation, diesel-electric locomotives, and both surface ships and submarines.
• In analysis of this cycle some air standard assumptions are taken which are:-

1. The working fluid is a fixed mass of air undergoing the complete cycle which is treated throughout as an ideal gas.
2. All processes are ideal, combustion is replaced by heat addition to the air, and exhaust is replaced by a heat rejection process which restores the air to the initial state.

GOVERNING EQUATION:-

• Diesel cycle consists four processes which are

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

Isentropic,

1) PV^ϒ = Constant

P1*V1 ^ϒ = P2*V2 ^ϒ

^{`P2 = P1*[(V1)/(V2)]^\gamma`}

^{2) `(P1*V1)/(T1) = (P2*V2)/(T2)`}

^{`T2 = (P2*V2*T1)/(P1*V1)}

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

1) P2 = P3

2) V3 = V2*co

3) `(P3*V3)/(T3) = (P2*V2)/(T2)`

`T3 = (V3*T2)/(V2)`

3) 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`

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

Where,

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

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

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

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

ϒ = ratio of specific heat (C_p/C_v)

Cr = compression ratio (V1/V2)

Co = cut-off ratio (V3/V2)

• 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/(Vc) = 1 + 1/2(cr-1)[R+1-cos(\theta)-\sqrt(R^2-sin^2(\theta))]`

Where,

V = Volume during compression and expansion

Vc = Clearance volume

Cr = compression ratio

R =`l/a`

a = Crank pin radius

l = connecting rod length

θ = crank angle.

• Thermal Efficiency of diesel cycle:-

`(\eta) = 1 - (1/(cr^(\gamma-1)))*((co^\gamma-1)/(\gamma*(co-1)))`

Where,

Cr = Compression ratio

Co = cut-off ratio

ϒ = ratio of specific heat (C_p/C_v)

OBJECTIVE:-

• Program to calculate parameters at all states in of Diesel cycle.
• Program to plot P-V diagram of Diesel cycle after calculating parameters.

PROGRAM:-

• Program for calculating parameters:-

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

% inputs for calculation of other parameters 
gamma = 1.4                       % Specific heat ratio


% state variables
P1 = 101325                       % Pressure at beginning og compression(N/m2)
T1 = 300                          % 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 = 18 ;                         % compression ratio(V1/V2)
co = 2 ;                          % cut off ratio (V3/V2)

% 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 Cilinder head & piston head when piston at TDC (m3)

V1 = V_swept + V_clearance        % Total Volume of Cyclinder (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 at conctant pressure (2-3)
P3 = P2                            % Pressure at beginning of expansion(N/m2)

V3 = V2*co
E =0.63

% P3*V3/T3 = P2*V2/T2 | T3 = T2*V3/V2
T3 = T2*V3/V2                   

V_heat_addition = linspace(4.62*10^-5,9.24*10^-5,180);
P_heat_addition = linspace(5795596,5795596,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 = P2*V3^gamma;
V_expansion = engine_kinematics_piston(D,L,l,cr,24.37389,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,267360,180);


% 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 diesel cycle','fontsize',25)
xlabel('Volume (m3)','fontsize',15)
ylabel('Pressure (N/m2)','fontsize',15)
grid on, axis square
% labeling the figure
text(V1,P1, 'State-1','fontsize',12)
text(V2,P2, 'State-2','fontsize',12)
text(V3,P3, 'State-3','fontsize',12)
text(V4,P4, 'State-4','fontsize',12)
text(V_compression(140),P_compression(140),'Compression','fontsize',12)
text(V_expansion(20),P_expansion(20),'Expansion','fontsize',12)

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

Output:-

• Program for plotting parameters:-

From this program we will get P-v diagram for diesel 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 Diesel cycle')

xlabel('Volume(m3)')

ylabel('Pressure(N/m2)')

grid on

Output:-

• 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:-

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 = linspace(4.62*10^-5,9.24*10^-5,180);

(During heat addition Volume changes from V2 to V3)

P_heat_addition = linspace(5795596,5795596,180);

(During heat addition pressure remains constant)

3) Equation for piston kinematics during expansion stroke

constant_c = P2*V3^gamma;

V_expansion = engine_kinematics_piston(D,L,l,cr,24.37389,180);

(Expansion process volume changes from V3 to V4)

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,267360,180);

(During heat rejection pressure changes from P3 to P4 )

Output:-

CONCLUSION:-

• The area bounded by the complete cycle path represents the total work that can be done during one cycle.
• 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.
• By this plot we can get an idea about performance of engine operated on diesel cycle.
• By this programming we can compare two different cycles for better performance and efficiency.