Dual cycle simulator using MATLAB-Final Independent Project

Aim

The aim of this program is to simulate a dual cycle and calculate the cycle's thermal efficiency, net work done and the mean effective pressure. The variation of the efficiency against cycle parameters such as compression ratio, cutoff ratio and alpha value are plotted.

Dual cycle

1-2 is isentropic compression

2-2.5 is heat addition at constant volume

2.5-3 is heat addition at constant pressure

3-4 is isentropic expansion

4-1 is heat rejection at constant volume

Compression ratio(r)=v1/v2

Cutoff ratio(r_c)=v3/v2

Alpha value=P2.5/P2

Cutoff ratio = 1 signifies that the cycle is an otto cycle and alpha value =1 signifies the cycle is a diesel cycle.

For this program compression ratio values between 8 and 23 are accepted. Value from 12 to 16 are considered for dual cycle, 8 to 11 are considered for otto cycle and 14 to 23 are considered for diesel cycle. Based on the cycle parameters inputted by the user the program decides which cycle it is simulating.

Formulae used

1. `P*v=R*T`
2. `P*v^gamma=const`
3. `eta=1-(1/r^(gamma-1))*[(alpha*r_c^gamma-1)/(alpha-1+alpha*gamma*(r_c-1))]`
4. `q_i=c_v*(T_2.5-T_2)+c_p*(T_3-T_2.5)`
5. `w=q_i*eta`
6. `mep=w/(v_1-v_2);`

Program

The temperature and pressure at state 1 are taken as inputs along with the alpha value, cutoff ratio and compression ratio. A sample input for dual cycle is given below.

A message dialogue box is created using 'msgbox' to display error messages to the user when invalid values are entered.

f=msgbox('Invalid input','Error','error');

The cycle to be simulated is based on the cycle parameters. The pressure, temperature and specific volume of each state are then calculated using equations 1 and 2. Pressure and specific volume values during the isentropic compression and expansion processes are also calculated.

The cycle is then plotted along with an appropriate legend. If the cycle is a dual cycle then the variation of its efficiency with compression ratio, cutoff ratio and alpha value are plotted in another figure. The thermal efficiency, net work done and mean effective pressure are calculated using the formulae given above and displayed using a message dialolgue box. 'sprintf' is used to display both the message and values.

In both plots 'uicontrol' command is used to create a push button in the plot figure to save the plot as a png file. The uicontrol object is stored in 'button'.

button=uicontrol(1,'String','Save','Callback',{@saveclicked,1});

The callback function 'saveclicked' is used to save the plot in the folder 'Dual cycle' in the current directory. The name of the image is entered by the user through a dialogue box using 'inputdlg' command and the name is stored in 'imgname'. Using 'saveas' command the plot is stored as a png file.

function saveclicked(source,event,num)
imgname=inputdlg('Enter image name','Image name',[1 50]);
pathname="C:\Users\Venkat JR\Documents\MATLAB\Results images";
mkdir(pathname,'Dual cycle')
newpathname=pathname+"\"+"Dual cycle";
if(isempty(imgname)||strcmp(imgname{1},''))
    f=msgbox('Invalid input','Error','error');
    return;
end
saveas(num,fullfile(newpathname,imgname{1}),'png');
end

Conclusion

The following inferences can be made from the efficiency comparison plot below:

1. A higher compression ratio is desirable for better efficiency although this is limited by the size of the engine.
2. A lower cutoff ratio is desirable for better efficiency
3. The efficiency does not vary much with alpha value after a point. An optimal efficiency can be obtained with a lower maximum pressure.

Plots and Results

Dual cycle

Otto cycle and Diesel cycle

Code

clear all
close all
clc
b=true;
%input
P1=input('Enter the pressure at state 1 in Pascal\n');
T1=input('Enter the temperature at state 1 in Kelvin\n');

disp('Alpha value is the ratio between pressure at state 2.5 and state 2');
disp('Alpha value equal to 1 implies its a diesel cycle');
while(b)
    alpha=input('Enter alpha value\n');
    if(alpha<1)
        f=msgbox('Invalid input','Error','error');
        disp('Please enter valid input');
    else
        b=0; 
    end
end
disp('Cutoff ratio(r_c) is the ratio between volume at state 3 and state 2');
disp('Cutoff ratio equal to 1 implies its an otto cycle');
while(~b)
    cutoffr=input('Enter cutoff ratio\n');
    if(cutoffr<1)
        f=msgbox('Invalid input','Error','error');
        disp('Please enter valid input');
    else
        b=1; 
    end
end

disp('Compression ratio(r) is the ratio between volume at state 1 and state 2');
disp('Compression ratio for dual cycle is between 12 and 16');
disp('Compression ratio for otto cycle is between 8 and 11');
disp('Compression ratio for diesel cycle is between 14 and 23');
ident=0; % to identify which cycle to simulate 1 is dual,2 is otto and 3 is diesel
while(b)
    compr=input('Enter compression ratio\n');
    if(compr<8 || compr>23)
        f=msgbox('Invalid input','Error','error');
        disp('Please enter valid input')
        clc
    elseif(compr<=11)
         cutoffr=1;
         b=0;
         titl='Otto cycle'; 
         ident=2;
    elseif(compr>=17)
        alpha=1;
        b=0;
        titl='Diesel cycle';
        ident=3;
    elseif(compr>=14 && compr<=16)
        b=0;
        if(alpha==1)
            titl='Diesel cycle';
            ident=3;
        else
            titl='Dual cycle';
            ident=1;
        end
    else
        b=0;
        titl='Dual cycle';
        ident=1;
    end
end
gamma=1.4;
cp=1.005;
cv=0.718;
%state 1
v1=0.287*T1/P1;

%state 2
P2=P1*((compr)^gamma);
v2=v1/compr;
T2=T1*(compr)^(gamma-1);

%state 2.5
v2_5=v2;
P2_5=P2*alpha;
T2_5=alpha*T2;

%state 3
T3=cutoffr*T2_5;
P3=P2_5;
v3=T3*v2_5/T2_5;

%state 4
v4=v1;
T4=T3/((compr/cutoffr)^(gamma-1));
P4=0.287*T4/v4;

% v and P values during isentropic compression and expansion
v_compression=linspace(v1,v2,500);
const=P1*(v1^gamma);
P_compression=const./((v_compression).^gamma);
v_expansion=linspace(v3,v4,500);
const=P3*v3^gamma;
P_expansion=const./((v_expansion).^gamma);

%plots
figure(1)
hold on

plot(v_compression,P_compression,'color','black');
if(ident==1||ident==2)
    plot([v2,v2_5],[P2,P2_5],'color','black');
end
if(ident==1||ident==3)
    plot([v2_5,v3],[P2_5,P3],'color','black');
end
plot(v_expansion,P_expansion,'color','black');
plot([v4,v1],[P4,P1],'color','black');

plot(v1,P1,'.','MarkerSize',8,'color','red');
plot(v2,P2,'.','MarkerSize',8,'color','red');
plot(v2_5,P2_5,'.','MarkerSize',8,'color','red');
plot(v3,P3,'.','MarkerSize',8,'color','red');
plot(v4,P4,'.','MarkerSize',8,'color','red');

title(titl);
xlabel('Specific volume(m^3/kg)');
ylabel('Pressure(Pa)'); 
%legends specific to the cycle being simulated
if(ident==1)
    legend('1-2 isentropic compression','2-2.5 constant volume heat addition','2.5-3 constant pressure heat addition','3-4 isentropic expansion','4-1 constant volume heat rejection');
elseif(ident==2)
    legend('1-2 isentropic compression','2-3 constant volume heat addition','3-4 isentropic expansion','4-1 constant volume heat rejection');
else
    legend('1-2 isentropic compression','2-3 constant pressure heat addition','3-4 isentropic expansion','4-1 constant volume heat rejection');
end
hold off
button=uicontrol(1,'String','Save','Callback',{@saveclicked,1});

%efficiency comparisons with cycle parameters if its a dual cycle
if(ident==1)
    figure(2)
    tcutoffr=1.2;
    talpha=1.5;
    tcompr=14;
    
    subplot(3,1,1);
    compr_range=linspace(12,16,10);
    eff1=1-(1./compr_range.^(gamma-1)).*((talpha*tcutoffr^gamma-1)/(talpha-1+talpha*gamma*(tcutoffr-1)));
    eff1=eff1.*100;
    plot(compr_range,eff1,'-o');
    xlabel('Compression ratio');
    ylabel('Efficiency(%)');
    title('Thermal efficiency vs Compression ratio(r_c=1.2,\alpha=1.5)');
    
    subplot(3,1,2);
    cutoffr_range=linspace(1.1,5,10);
    eff1=1-(1/tcompr^(gamma-1))*((talpha*cutoffr_range.^gamma-1)./(talpha-1+talpha*gamma*(cutoffr_range-1)));
    eff1=eff1.*100;
    plot(cutoffr_range,eff1,'-o');
    xlabel('Cutoff ratio');
    ylabel('Efficiency(%)');
    title('Thermal efficiency vs Cutoff ratio(r=14,\alpha=1.5)');
    
    subplot(3,1,3)
    alpha_range=linspace(1.1,5,10);
    eff1=1-(1/tcompr^(gamma-1))*((alpha_range.*tcutoffr^gamma-1)./(alpha_range-1+alpha_range.*gamma*(tcutoffr-1)));
    eff1=eff1.*100;
    plot(alpha_range,eff1,'-o');
    xlabel('Alpha value');
    ylabel('Efficiency(%)');
    title('Thermal efficiency vs Alpha value(r_c=1.2,r=14)');
    
    button=uicontrol(2,'String','Save','Callback',{@saveclicked,2});
end

%results
efficiency= 1-(1/compr^(gamma-1))*((alpha*cutoffr^gamma-1)/(alpha-1+alpha*gamma*(cutoffr-1)));
qin=cv*(T2_5-T2)+cp*(T3-T2_5);
wnet=qin*efficiency;
mep=wnet/(v1-v2);
f=msgbox({sprintf('Thermal efficiency=%g',efficiency*100);sprintf('Net work=%g kJ/kg',wnet);sprintf('Mean effective pressure=%g kPa',mep)},'Results');

%call back function
function saveclicked(source,event,num)
imgname=inputdlg('Enter image name','Image name',[1 50]);
pathname="C:\Users\Venkat JR\Documents\MATLAB\Results images";
mkdir(pathname,'Dual cycle')
newpathname=pathname+"\"+"Dual cycle";
if(isempty(imgname)||strcmp(imgname{1},''))
    f=msgbox('Invalid input','Error','error');
    return;
end
saveas(num,fullfile(newpathname,imgname{1}),'png');
end

References

1. Thermodynamics an engineering approach by Yunus A. Cengel and Michael A. Boles
2. https://in.mathworks.com/help/matlab/ref/uicontrol.html
3. https://in.mathworks.com/help/matlab/ref/msgbox.html