Rankine cycle simulator using MATLAB
Aim
The aim of this program is to simulate a rankine cycle and plot its T-s diagram and h-s diagrams. It also calculates the efficiency and back work ratio of the cycle based on the users inputs.
In this rankine cycle:
1-2 is isentropic compression in a pump
2-3 is constant pressure heat addition in a boiler
3-4 is isentropic expansion in a turbine
4-1 is constant pressure heat rejection in a condenser
Program
The program takes the users input for pressure and temperature at turbine inlet. Since the water enters the turbine as a superheated vapour the program checks if the inputted temperature is greater or equal to the saturation temperature at that pressure. If not the user is asked to input again. If it is greater it means state 3 is a superheated vapour and if it's equal then its a saturated vapour. In both cases the specific enthalpy and specific entropy for state 3 are calculated.
The user is then asked to enter the pressure at the turbine exit. Since water leaving the turbine is not a superheated vapour but a saturated mixture, the pressure value is checked if it is valid.
The values for pressure(P), temperature(T), specific enthalpy(h), specific entropy(s) are calculated for states 1,2,3 and 4. The work produced by the turbine, consumed by the pump and net work is calculated along with the back work ratio and efficiency of the cylce.
The cycle is plotted in T-s diagram and h-s diagram with the saturated liquid and saturated vapour lines(red).
The results are displayed in the command window.
Code
clear all
close all
clc
disp('1-2 is isentropic compression in pump');
disp('2-3 is constant pressure heat addition in boiler');
disp('3-4 is isentropic expansion in turbine');
disp('4-1 is constant pressure heat rejection in condenser');
P3=input('Enter pressure at turbine inlet in bar\n');
b=1;
while(b~=0)
T3=input('Enter temperature at turbine input in degree celsius\n');
if(T3>XSteam('Tsat_p',P3))
s3=XSteam('s_pT',P3,T3);
h3=XSteam('h_pT',P3,T3);
b=0;
elseif(T3==XSteam('Tsat_p',P3))
s3=XSteam('sV_p',P3);
h3=XSteam('hV_p',P3);
b=0;
else
fprintf('Temperature should be greater than or equal to saturation temperature: %f\n',XSteam('Tsat_p',P3));
end
end
while(b==0)
s4=s3;
P4=input('Enter pressure at turbine exit or condenser inlet\n');
htemp=XSteam('h_ps',P4,s4);
if(htemp>XSteam('hV_p',P4))
disp('Enter a lower pressure as steam exiting turbine cannot be superheated');
else
b=1;
end
end
P2=P3;
P1=P4;
%state 1 into pump
v1=XSteam('vL_p',P1);
w_pumpin=v1*(P2-P1)*100;
h1=XSteam('hL_p',P1);
s1=XSteam('sL_p',P1);
T1=XSteam('T_ph',P1,h1);
%state2 into boiler
h2=h1+w_pumpin;
s2=s1;
T2=XSteam('T_ph',P2,h2);
%state 2'
Tsat2=XSteam('Tsat_p',P2);
ssat2=XSteam('sL_p',P2);
%state3 into turbine out of boiler
h3;
s3;
%state3'
Tsat3=XSteam('Tsat_p',P3);
ssat3=XSteam('sV_p',P3);
%state4 into condenser out of turbine
s4=s3;
s4f=XSteam('sL_p',P4);
s4g=XSteam('sV_p',P4);
s4fg=s4g-s4f;
x4=(s4-s4f)/s4fg;
h4=XSteam('h_px',P4,x4);
T4=XSteam('T_ph',P4,h4);
%efficiency and back work ratio
w_turbout=h3-h4;
w_net=w_turbout-w_pumpin;
qin=h3-h2;
eff=(w_net/qin)*100;
bwr=w_pumpin/w_turbout;
%data for plotting saturated liquid and saturated vapour lines
tcr=linspace(0,374,1000);
for i=1:length(tcr)
sf(i)=XSteam('sL_T',tcr(i));
sg(i)=XSteam('sV_T',tcr(i));
hf(i)=XSteam('hL_T',tcr(i));
hg(i)=XSteam('hV_T',tcr(i));
end
%data for plotting the isobaric process
trange1=linspace(T2,Tsat2,500);
trange2=linspace(Tsat3,T3,500);
for i=1:length(trange1)
srange1(i)=XSteam('s_pT',P2,trange1(i));
srange2(i)=XSteam('s_pT',P2,trange2(i));
end
%T-s graph
figure(1)
hold on
%plotting points
plot(s1,T1,'o','color','black');
text(s1,T1,'1');
plot(s2,T2,'o','color','black');
text(s2,T2,'2');
plot(s3,T3,'o','color','black');
text(s3,T3,'3');
plot(s4,T4,'o','color','black');
text(s4,T4,'4');
plot(ssat2,Tsat2,'o','color','black');
text(ssat2,Tsat2,"2'");
plot(ssat3,Tsat3,'o','color','black');
text(ssat3,Tsat3,"3'");
plot(sf,tcr,'color','red');
plot(sg,tcr,'color','red');
plot([s1,s2],[T1,T2],'color','black');
plot([s3,s4],[T3,T4],'color','black');
plot([s4,s1],[T4,T1],'color','black');
plot([ssat2,ssat3],[Tsat2,Tsat3],'color','black');
plot(srange1,trange1,'color','black');
plot(srange2,trange2,'color','black');
hold off
xlabel('Entropy(kJ/kg-K)');
ylabel('Temperature(C)');
title('Rankine cycle T-s graph')
%h-s graph
figure(2)
hold on
plot(s1,h1,'o','color','black');
text(s1,h1,'1');
plot(s2,h2,'o','color','black');
text(s2,h2,'2');
plot(s3,h3,'o','color','black');
text(s3,h3,'3');
plot(s4,h4,'o','color','black');
text(s4,h4,'4');
plot([s1,s2],[h1,h2],'color','black');
plot([s2,s3],[h2,h3],'color','black');
plot([s3,s4],[h3,h4],'color','black');
plot([s4,s1],[h4,h1],'color','black');
plot(sf,hf,'color','red');
plot(sg,hf,'color','red');
hold off
xlabel('Entropy(kJ/kg-K)');
ylabel('Enthalpy(kJ/kg)');
title('Rankine cycle h-s graph');
%Results
disp('Results');
disp('State 1');
fprintf('P1 is: %f bar\n',P1);
fprintf('T1 is: %f degree celsius\n',T1);
fprintf('s1 is: %f kJ/kg-K\n',s1);
fprintf('h1 is: %f kJ/kg\n',h1);
disp('State 2');
fprintf('P2 is: %f bar\n',P2);
fprintf('T2 is: %f degree celsius\n',T2);
fprintf('s2 is: %f kJ/kg-K\n',s2);
fprintf('h2 is: %f kJ/kg\n',h2);
disp('State 3');
fprintf('P3 is: %f bar\n',P3);
fprintf('T3 is: %f degree celsius\n',T3);
fprintf('s3 is: %f kJ/kg-K\n',s3);
fprintf('h3 is: %f kJ/kg\n',h3);
disp('State 4');
fprintf('P4 is: %f bar\n',P4);
fprintf('T4 is: %f degree celsius\n',T4);
fprintf('s4 is: %f kJ/kg-K\n',s4);
fprintf('h4 is: %f kJ/kg\n',h4);
disp('Efficiency');
fprintf('Work turbine: %f kJ/kg\n',w_turbout);
fprintf('Work pump: %f kJ/kg\n',w_pumpin);
fprintf('Work net: %f kJ/kg\n',w_net);
fprintf('Efficiency: %f percentage\n',eff);
fprintf('Back work ratio: %f\n',bwr)Results
Plots
References
- https://www.learnthermo.com/T1-tutorial/ch09/lesson-B/pg07.php
- Thermodynamics an engineering approach by Yunus A. Cengel and Michael A. Boles