To simulate Dual cycle using MATLAB.
AIM - To simulate Dual cycle and calculate the Dual cycle performance parameters using MATLAB.
Theory -
The dual combustion cycle (also known as the mixed cycle, Trinkler cycle, Seliger cycle or Sabathe cycle) is a thermal cycle that is a combination of the Otto cycle and the Diesel cycle, first introduced by Russian-German engineer Gustav Trinkler, who never claimed to have developed the cycle though. Heat is added partly at constant volume (isochoric) and partly at constant pressure (isobaric), the significance of which is that more time is available for the fuel to completely combust. Because of lagging characteristics of fuel this cycle is invariably used for Diesel and ignition engines. It consists of two adiabatic and two constant volume and one constant pressure processes.
The dual cycle consists of following operations:
- Process 1-2: Isentropic compression
- Process 2-3: Addition of heat at constant volume.
- Process 3-4: Addition of heat at constant pressure.
- Process 4-5: Isentropic expansion.
- Process 5-1: Rejection of heat at constant volume.
The EGR rate is defined as the percentage of exhaust gases from the total gas mass induced into the engine. For example, an EGR rate of 33 % means that a third of the gas which enters the cylinders is actually exhaust gas and 67 % is fresh air. The higher the EGR rate the lower the NOx emissions levels.
Assumptions in parameters -
1.Default set data-
V=4; Engine Displacement Volume(Litr)
N=4; Number of Cylinders
n=2300; Engine Speed(rpm)
`r_c`=16; Compression Ratio
X=1.27; Stroke-to-Bore Ratio S/B
`eta_m`=0.86; Mechanical Efficiency
EGR=2; Internal EGR Percentage
cvp=20; Percent of Heat Added During Constant-Volume Combustion Process
cpp=100-cvp; Percent of Heat Added During Constant-Pressure Combustion Process
x=12; number of Carbon atoms
y=23; number of Hydrogen atoms
AF=18; Air-Fuel Mass Ratio
`eta_(conv)`=1; Combustion Efficiency
K=1.35; Specific Heat Ratio of Air
R=0.287; Particular Gas Constant of Air (kJ/kg.K)
`C_v`=0.821; Specific Heat for Constant Volume of Air (kJ/kg)
`C_p`=R+C_v; Specific Heat for Constant Pressure of Air (kJ/kg)
P_0=100; Atmospheric Pressure (kPa)
T_0=60; Combustion Chamber Temperature at Start of Compression (Deg C)
`rho_0`=P_0/(R*(293.15)); Density of Air at Reference Condition (kg/cubic meter)
2.User input assumption -
Engine Displacement Volume (Litres): 3
Number of Cylinders: 4
Engine Speed(rpm): 3000
Compression Ratio: 16
Stroke-to-Bore Ratio S/B: 1.5
Efficiency:0.80
EGR Percentage: 3%
Heat Added During Constant-Volume Combustion Process: 20%
Fuel Composition
Carbon atoms: 13
Hydrogen atoms: 20
Import Air-Fuel Mass Ratio: 18
Import Combustion Efficiency: 1
Import Atmospheric Pressure (kPa): 102
Import Combustion Chamber Temperature at Start of Compression (Deg C): 65
Import Particular Gas Constant of Air (kJ/kg.K): 0.287
Import Specific Heat Ratio of Air: 1.30
Import Specific Heat for Constant Volume of Air: (kJ/kg) 0.821
Objectives -
1. To calculate the Dual cycle paramters.
2. To obtain P-V diagram of the Dual Cycle.
Procedure -
Main Program -
clear all
close all
clc
%To simulate Dual Cycle with set data and the user input using MATLAB.
%--------------------------INPUT DATA----------------------------
a=input('Please enter 1 for default input data, 0 for arbitrary ones:');
if a==1
%-------------------------Default Input Data---------------------
%----------Engine Inputs------
V=4; %Engine Displacement Volume(Litr)
N=4; %Number of Cylinders
n=2300; %Engine Speed(rpm)
r_c=16; %Compression Ratio
X=1.27; %Stroke-to-Bore Ratio S/B
eta_m=0.86; %Mechanical Efficiency
EGR=2; %Internal EGR Percentage
cvp=20; %Percent of Heat Added During Constant-Volume Combustion Process
cpp=100-cvp; %Percent of Heat Added During Constant-Pressure Combustion Process
%---------Fuel Inputs/Combustion Efficiency----
x=12; %number of Carbon atoms
y=23; %number of Hydrogen atoms
AF=18; %Air-Fuel Mass Ratio
eta_conv=1; %Combustion Efficiency
%Properties of Working Fluids----------------
%a constant average value between two exteremes ( SOC & EOCombustion)
K=1.35; %Specific Heat Ratio of Air
R=0.287; %Particular Gas Constant of Air (kJ/kg.K)
C_v=0.821; %Specific Heat for Constant Volume of Air (kJ/kg)
C_p=R+C_v; %Specific Heat for Constant Pressure of Air (kJ/kg)
%---------Atmospheric Inputs------------
P_0=100; %Atmospheric Pressure (kPa)
T_0=60; %Combustion Chamber Temperature at Start of Compression (Deg C)
row_0=P_0/(R*(293.15)); %Density of Air at Reference Condition (kg/cubic meter)
%-----------------------------------------------------------------
else
%-----------------------Arbitrary Input Data------------------
%----------Engine Inputs------
V=input('Import Engine Displacement Volume(Litr): ');
N=input('Import Number of Cylinders: ');
n=input('Import Engine Speed(rpm): ');
r_c=input('Import Compression Ratio: ');
X=input('Import Stroke-to-Bore Ratio S/B: ');
eta_m=input('Import Mechanical Efficiency:');
EGR=input('Import Internal EGR Percentage: % ');
cvp=input('Import Percent of Heat Added During Constant-Volume Combustion Process: % ');
cpp=100-cvp;
%---------Fuel Inputs/Efficiencies----
disp('---------------CxHy Fuel Composition------------ ');
x=input('Please import number of Carbon atoms: ');
y=input('Please import number of Hydrogen atoms: ');
b=x/y;
AF_s=(34.56*(b+4))/(12.011+(1.008*b)); %Stoichiometric Air/Fuel ratio (Heywood (5-3))
Z = ['For Complete Combustion Assumption The Air-Fuel Mass Ratio Must be ',num2str(AF_s),' or Higher'];
disp(Z)
AF=input('Import Air-Fuel Mass Ratio: ');
eta_conv=input('Import Combustion Efficiency: ');
%---------Atmospheric Inputs------------
P_0=input('Import Atmospheric Pressure (kPa): ');
T_0=input('Import Combustion Chamber Temperature at Start of Compression(Deg C): ');
R=input('Import Particular Gas Constant of Air (kJ/kg.K): ');
K=input('Import Specific Heat Ratio of Air: ');
C_v=input('Import Specific Heat for Constant Volume of Air: (kJ/kg) ');
C_p=R+C_v; %Specific Heat for Constant Pressure of Air (kJ/kg)
row_0=P_0/(R*(293.15)); %Density of Air at Reference Condition (kg/cubic meter)
%-----------------------------------------------------
end
%--------------------------Calculation---------------------------
%--------Fuel Heating Value & Combustion Efficiency----
M_C=(x*12)/((x*12)+(y)); %Mass Fraction of Carbon In Fuel
M_H=(y*0.68)/((x*12)+(y)); %MODIFIED Mass Fraction of Hydrogen In Fuel
FHV=((14600*M_C)+(62000*M_H))*2.326; %Fuel Heating Value (kJ/kg) [Dulong’s formula]
b=x/y;
AF_s=(34.56*(b+4))/(12.011+(1.008*b)); %Stoichiometric Air/Fuel ratio (Heywood (5-3))
lambda=AF/AF_s;
Th_Air=lambda*100; %Theoritical Air Percentage
EAP=Th_Air-100; %Exess Air Percentage
if lambda >= 1
eta_ch=1;
else
eta_ch=(1.3773*lambda)-0.3773;
end
eta_c=eta_conv*eta_ch;
%-----------------------------------------------------------------
%------for one cylinder
V_d=(V/N)/1000; %Displacement Volume (Cubic Meter)
V_c=V_d/(r_c-1); %Clearance Volume (Cubic Meter)
B=((4*V_d/pi*X)^(1/3))*100; %Bore (cm)
S=B*X; %Stroke (cm)
A_p=(pi*(B^2))/4 ; %Piston Area (square cm)
%----State 1 (Start of Compression)-----
T_1=(T_0)+273; %Temprature of Gas Mixture in Cylinder at State 1 (Deg K)
P_1=P_0; %Pressure of Gas Mixture in Cylinder at State 1 (kPa)
V_1=V_d+V_c; %Volume of Gas Mixture in Cylinder at State 1 (Cubic Meter)
%----Mass of Gases Calculation------
M_m=(P_1*V_1)/(R*T_1); %Mass of Gas Mixture in Cylinder (kg)
%the mass within the cylinder will then remain the same for the entire
%cycle
M_a=(AF/(AF+1))*((100-EGR)/100)*M_m; %Mass of Air in Cylinder (kg)
M_f=(1/(AF+1))*((100-EGR)/100)*M_m; %Mass of Fuel in Cylinder (kg)
M_ex=(EGR/100)*M_m; %Mass of Exhaust Residual in Cylinder at SOC (kg)
%------------State 2 (Start of Combustion)-----
T_2=(T_1)*((r_c)^(K-1)); %Temprature of Gas Mixture in Cylinder at State 2 (Deg K)
P_2=(P_1)*((r_c)^K); %Pressure of Gas Mixture in Cylinder at State 2 (kPa)
V_2=V_c; %Volume of Gas Mixture in Cylinder at State 2 (Cubic Meter)
%------Heat Added During Two Combustion Processes-------------
Q_in=M_f*FHV*eta_c; %Total Heat Added During Combustion
Q_cv=(cvp/100)*Q_in; %Heat Added During Constant-Volume Combustion
Q_cp=(cpp/100)*Q_in; %Heat Added During Constant-Volume Combustion
%------------State x (End of Constant-Volume Combustion)-----
T_x=(Q_cv/(M_m*C_v))+T_2; %Temprature of Gas Mixture in Cylinder at State x (Deg K)
P_x=P_2*(T_x/T_2); %Pressure of Gas Mixture in Cylinder at State x (kPa)
V_x=V_2; %Volume of Gas Mixture in Cylinder at State 3 (Cubic Meter)
%------------State 3 (End of Combustion)-----
T_3=((Q_cp)/(M_m*C_p))+T_x; %Temprature of Gas Mixture in Cylinder at State 3 (Deg K)
P_3=P_x; %Pressure of Gas Mixture in Cylinder at State 3 (kPa)
V_3=(T_3/T_x)*V_x; %Volume of Gas Mixture in Cylinder at State 3 (Cubic Meter)
%----Cutoff Ratio (Change in Volume During Combustion)----
Beta=V_3/V_x;
%----Pressure Ratio (Change in Pressure During Combustion)----
Alpha=P_x/P_2;
%------------State 4 (End of Expression)-----
V_4=V_1; %Volume of Gas Mixture in Cylinder at State 4 (Cubic Meter)
T_4=T_3*((V_3/V_4)^(K-1)); %Temprature of Gas Mixture in Cylinder at State 4 (Deg K)
P_4=P_3*((V_3/V_4)^K); %Pressure of Gas Mixture in Cylinder at State 4 (kPa)
%------------State 5
T_5=T_1; %Temprature of Gas Mixture in Cylinder at State 5 (Deg K)
P_5=P_1; %Pressure of Gas Mixture in Cylinder at State 5 (kPa)
V_5=V_1; %Volume of Gas Mixture in Cylinder at State 5 (Cubic Meter)
%------------State 6
T_6=T_0; %Temprature of Gas Mixture in Cylinder at State 6 (Deg K)
P_6=P_0; %Pressure of Gas Mixture in Cylinder at State 6 (kPa)
V_6=V_2; %Volume of Gas Mixture in Cylinder at State 6 (Cubic Meter)
%------------------------------------------------------------------------------------------------------
%---------------------------------Processes------------------------------
%--- 1-2 Process (Compression)
Q_12=0;
W_12=M_m*C_v*(T_1-T_2); %(kJ)
%--- 2-x Process (Constant-Volume Combustion)
Q_2x=M_m*C_v*(T_x-T_2); %(kJ)
W_2x=0;
%--- x-3 Process (Constant-Pressure Combustion)
Q_x3=M_m*C_p*(T_3-T_x); %(kJ)
W_x3=P_3*(V_3-V_x); %(kJ)
%--- 3-4 Process (Expression)
Q_34=0;
W_34=M_m*C_v*(T_3-T_4); %(kJ)
%--- 4-5 Process (Exhaust)
Q_45=M_m*C_v*(T_5-T_4); %(kJ)
W_45=0;
%------------------------------ Overal Engine Performance----------
U_p=2*S*n*(1/6000); %Mean Piston Speed (m/s)
W_i=W_12+W_x3+W_34; %Net Indicated Work for One Cylinder During One Cycle (kJ)
eta_t=W_i/Q_in; %Thermal Efficiency
imep=W_i/(V_d); %Indicated Mean Effective Pressure (kPa)
bmep=imep*eta_m; %Brake Mean Effective Pressure (kPa)
W_idat=W_i*N*(n/(60*2)); %Indicated Power at {n} RPM (KW)
W_b=eta_m*W_i; %Net Brake Work for One Cylinder During One Cycle (kJ)
W_bdat=W_b*N*(n/(60*2)); %Brake Power at {n} RPM (KW)
Taw=(W_bdat/(2*pi*(n/60)))*1000; %Torque (N-m)
W_fdat=W_idat-W_bdat; %Friction Power Lost (KW)
BSP=W_bdat/A_p; %Brake Specific Power (KW/square cm)
OPD=W_bdat/V; %Outpot per Displacement (KW/Litr)
M_fdat=M_f*(n/60)*(1/2)*N; %Rate of Fuel That Inject in All Cylinder During One Cycle (kg/sec)
bsfc=(M_fdat/W_bdat)*3600000; %Brake Specific Fuel Consumption (gr/KW-hr)
eta_v=M_a/(row_0*V_d); %Volumetric Efficiency
%--------------------Display Outputs--------------------------
disp(' ');
disp(' ');
disp('---------------------------Fuel & Combustion Properties-------------------------- ');
F=[x y];
disp(' C H')
disp(F)
A = ['The Actual Air/Fuel Ratio is: ',num2str(AF)];
disp(A)
Z = ['The Stoichiometric Air/Fuel Ratio is: ',num2str(AF_s)];
disp(Z)
T=['The Percent of Theoritical Air is: ',num2str(Th_Air)];
disp(T)
E=['The Percent of Excess Air is: ',num2str(EAP)];
disp(E)
Eff=['Total Combustion Efficiency is: ',num2str(eta_c)];
disp(Eff)
HV=['Fuel Heating Value is: ',num2str(FHV), ' kJ/kg'];
disp(HV)
disp(' ');
disp(' ');
disp('---------------------------Cycle Operating Point Properties-------------------------- ');
disp(' ');
disp(' ');
aa=['The Maximum Temprature of The Cycle is: ',num2str(T_3-273.15), ' Deg C'];
disp(aa)
ab=['The Maximum Pressure of The Cycle is : ',num2str(P_3), ' kPa'];
disp(ab)
disp(' ');
disp(' ');
disp('---------------------------Overal Engine Performances-------------------------- ');
disp(' ');
disp(' ');
a15=['The Mean Piston Speed is: ',num2str(U_p),' m/s'];
disp(a15)
a1=['The Thermal Efficiency is: ',num2str(eta_t)];
disp(a1)
a2=['The Indicated Mean Effective Pressure imep is : ',num2str(imep), ' kPa'];
disp(a2)
a3 = ['The Brake Mean Effective Pressure is : ',num2str(bmep), ' kPa'];
disp(a3)
a5 = ['The Indicated Power at ',num2str(n),' RPM is : ',num2str(W_idat), ' KW'];
disp(a5)
a6 = ['The Brake Power at ',num2str(n),' RPM is : ',num2str(W_bdat), ' KW'];
disp(a6)
a8=['The Friction Power Lost is: ',num2str(W_fdat), ' KW'];
disp(a8)
a7 = ['The Torque at ',num2str(n),' RPM is : ',num2str(Taw), ' N-m'];
disp(a7)
a9=['The Brake Specific Fuel Consumption bsfc is : ',num2str(bsfc), ' gr/KW-hr'];
disp(a9)
a12 = ['The Volumetric Efficiency is : ',num2str(eta_v)];
disp(a12)
a10 = ['The Outpot per Displacement OPD is : ',num2str(OPD), ' KW/Litr'];
disp(a10)
a11 = ['The Brake Specific Power BSP is : ',num2str(BSP), ' KW/cm^2'];
disp(a11)
a13 = ['The Pressure Increase Ratio During Contant-Volume Combustion is : ',num2str(Alpha)];
disp(a13)
a14 = ['The CutOff Ratio is : ',num2str(Beta)];
disp(a14)
%-------------------%Plot Statements----------------------
%---Plot 1-2 Process (Compressio)
V_12=linspace(V_1,V_2,100);
P_12=((V_1./V_12).^K)*P_1;
plot(V_12,P_12);
title('P-V Diagram','FontSize', 20 , 'FontName', 'Helvetica','FontWeight','bold');xlabel('Volume (m^3)','FontSize', 14 , 'FontName', 'Helvetica','FontWeight','bold');ylabel('Pressure (kPa)','FontSize', 14 , 'FontName', 'Helvetica','FontWeight','bold');
hold on
%---Plot 2-x Process (Constant-Volume Combustion)
V_2x=[V_2 V_x];
P_2x=[P_2 P_x];
plot(V_2x,P_2x);
%---Plot x-3 Process (Constant-Pressure Combustion)
V_x3=[V_x V_3];
P_x3=[P_x P_3];
plot(V_x3,P_x3);
%---Plot 3-4 Process (Expression)
V_34=linspace(V_3,V_4,100);
P_34=((V_3./V_34).^K)*P_3;
plot(V_34,P_34);
%---Plot 4-1 Process (Exhaust)
V_41=[V_4 V_1];
P_41=[P_4 P_1];
plot(V_41,P_41);
%---Plot 5-6 Process (Pumping-Expulsion)
V_56=[V_5 V_6];
P_56=[P_5 P_6];
plot(V_56,P_56);
hold off
text(V_1,P_1, '1,5' , 'FontSize' ,12)
text(V_2,P_2, '2' , 'FontSize' ,12)
text(V_x,P_x, 'x' , 'FontSize' ,12)
text(V_3,P_3, '3' , 'FontSize' ,12)
text(V_4,P_4, '4' , 'FontSize' ,12)
text(V_6,P_6, '6' , 'FontSize' ,12)
Results -
1.For default data set -
P-V diagram plot -
Command window -
2.For User inputs data -
P-V diagram plot -
Command window -
Conclusion -
The PV diagram plotted with calculation of engine performance parameters.