MATLAB - Rankine cycle Simulator

AIM

To create a Rankine Cycle Simulator using OCTAVE

CODE EXPLANATION

Rankine cycle

The thermodynamic cycle used in most of the turbines. It consists of Heat exchanger(HE) like boiler which heats the working fluid usually water, turbine where the steam expands, condenser where expanded steam condense, pump which circulates the water from condenser to HE.

In HE, the outlet will be superheated steam, consist of no moisture at all in order to protect turbine blades, whereas pump usually works well with subcooled liquid.

The following process happens

1. 1-2 Isentropic (entropy constant) Expansion in the Turbine
2. 2-3 Constant Pressure Heat Rejection by the Condenser
3. 3-2 Isentropic (entropy constant) Compression in the Pump
4. 4-1 Constant Pressure Heat Addition by the Boiler

CODE logic

1. Inputs are turbine inlet pressure and temperature and turbine outlet pressure. They are received using input command and stored in respective variables.
2. For first state point, h1 and s1 can be found only in superheated condition, using Xsteam they can be found.
3. For second state, we know s1=s2. using this we can find x2 which is nothing but degree of saturation(how much moisture it contains, 1=dry, 0=liquid) by the formula shown in code. from x2, h2 can be easily found as per formula shown. t2 can be found wrt pressure using Xsteam.
4. For third state, we know h3, s4 as they lie in left side of saturation line, so fluidtype at pressure=0.05bar. Also t2=t3.
5. For fourth state, we have to find pump work in order to find h4. Basically work=PdV (or) VdP and W=del(H), based on this Wpump found from which h4 found. We know s4=s3(isentropic process).
6. From these T-s and H-S plots can be plotted. For saturation line, for loop used in which Xsteam table values at saturation temperature ranfges stored.
7. Turbine work can be found using W=del(H)=h1-h2, similarly heat added q1=h1-h4. From these efficiency calculated, (Wt-Wpump)/q1.
8. specific steam consumption found by (3600/Wnet).

CODE

%Rankine cycle

clear all
close all
clc

%INTODUCTION
fprintf('                     RANKINE CYCLE SIMULATOR\n\n')                    
fprintf('1-2  Isentropic Expansion in the Turbine \n')
fprintf('2-3  Constant Pressure Heat Rejection by the Condenser \n')
fprintf('3-2  Isentropic Compression in the Pump \n')
fprintf('4-1  Constant Pressure Heat Addition by the Boiler \n\n')

%inputs
p1= input('Enter the Pressure at the Turbine Inlet (in bar): ')
t1= input('Enter the Temperature at the Turbine inlet (in Degree celsius): ')
p2= input('Enter the Pressure at the Condenser (in bar): ')
fprintf('\n                     RESULTS\n')

%state point 1 - b/w turbine inlet and HE
h1=XSteam('h_pt',p1,t1);
s1=XSteam('s_pt',p1,t1);
%display
fprintf('At state point 1: \n')
fprintf('P1 is : %.2f Bar\n',p1)
fprintf('T1 is : %.2f Deg celsius\n',t1)
fprintf('H1 is : %.2f kJ/kg\n',h1)
fprintf('S1 is : %.2f kJ/kgK\n\n',s1)

%state point 2 - b/w turbine outlet and condenser
t2=XSteam('Tsat_p',p2);
%isentropic process
s2=s1;
sf2=XSteam('sL_p',p2);
sg2=XSteam('sV_p',p2);
sfg2=sg2-sf2;
%s2=(sf2+(x2*sfg2))
x2=(s2-sf2)/sfg2;
hf2=XSteam('hL_p',p2);
hg2=XSteam('hV_p',p2);
hfg2=hg2-hf2;
h2=(hf2+(x2*hfg2));
%display
fprintf('At state point 2:\n')
fprintf('P2 is : %.2f Bar\n',p2)
fprintf('T2 is : %.2f Deg celsius\n',t2)
fprintf('H2 is : %.2f kJ/kg\n',h2)
fprintf('S2 is : %.2f kJ/kgK\n\n',s2)

%state point 3 - b/w pump and condenser outlet
t3=t2;
p3=p2;
s3=XSteam('sL_p',p3);            
h3=XSteam('hL_p',p3);
v3=XSteam('vL_p',p3);
%display
fprintf('At state point 3:\n')
fprintf('P3 is : %.2f Bar\n',p3)
fprintf('T3 is : %.2f Deg celsius\n',t3)
fprintf('H3 is : %.2f kJ/kg\n',h3)
fprintf('S3 is : %.2f kJ/kgK\n\n',s3)

%state point 4 - b/w pump and HE
p4=p1;
%isentropic process
s4=s3;
%to find t4
cp4=XSteam('Cp_ps',p4,s4);
cv4=XSteam('Cv_ps',p4,s4);
gama=cp4/cv4;
r=(gama-1)/gama;
t4=t3*(p4/p3)^r;
%pump work:W=PdV
wpump=v3*(p1-p2)*100;
%wpump=h4-h3
h4=wpump+h3;
%display
fprintf('At state point 4:\n')
fprintf('P4 is : %.2f Bar\n',p4)
fprintf('T4 is : %.2f Deg celsius\n',t4)
fprintf('H4 is : %.2f kJ/kg\n',h4)
fprintf('S4 is : %.2f kJ/kgK\n\n',s4)

%heat added
q1=h1-h4;
%turbne work
wt=h1-h2;

%efficiency
nth=(wt-wpump)/q1;

%specific steam consumption
ssc=3600/(wt-wpump);

%display
fprintf('Wt (turbine work) is : %.2f kJ/kg\n',wt)
fprintf('Wp (pump work) is : %.2f kJ/kg\n',wpump)
fprintf('Wnet (net work output) is : %.2f kJ/kg\n',(wt-wpump))
fprintf('Ntherm (thermal efficiency) is : %.2f percentage \n',(nth*100))
fprintf('S.S.C is : %.2f kJ/kWh\n',ssc)

t5=XSteam('Tsat_p',p1);
t6=t5;
s5=XSteam('sL_p',p1);
s6=XSteam('sV_p',p1);

%for saturation line/dome
T_curve = linspace(1,400,500);
for i=1:length(T_curve)
    sL(i) = XSteam('sL_T',T_curve(i));
    sV(i) = XSteam('sV_T',T_curve(i));
    hL(i) = XSteam('hL_T',T_curve(i));
    hV(i) = XSteam('hV_T',T_curve(i));
end

%plots
%T-S diagram
figure(1)
hold on
plot(sV,T_curve,'--','color','r')
plot(sL,T_curve,'--','color','r')
plot([s1 s2 s3 s4 s5 s6 s1],[t1 t2 t3 t4 t5 t6 t1],'color','b','linewidth',2);
text(s1,t1,' 1','fontsize',10)
text(s2,t2,' 2','fontsize',10)
text(s3,t3,'3','fontsize',10,'HorizontalAlignment','right','VerticalAlignment','top')
text(s4,t4,'4','fontsize',10,'HorizontalAlignment','right','VerticalAlignment','bottom')
xlabel('Entropy in kJ/KgK')
ylabel('Temperature in degree C')
title('T-S diagram')

%H-S diagram
figure(2)
hold on
plot(sV,hV,'--','color','r')
plot(sL,hL,'--','color','r')
plot([s1 s2 s3 s4 s1],[h1 h2 h3 h4 h1],'color','b','linewidth',2);
text(s1,h1,' 1','fontsize',10)
text(s2,h2,' 2','fontsize',10)
text(s3,h3,'3','fontsize',10,'HorizontalAlignment','right','VerticalAlignment','top')
text(s4,h4,'4','fontsize',10,'HorizontalAlignment','right','VerticalAlignment','bottom')
xlabel('Entropy in kJ/KgK')
ylabel('Enthalpy in kJ/kg')
title('H-S diagram')

RESULTS

OUTPUT

HS diagram

From the diagram we can see 1 is at right side of saturation dome reffering it is in superheated condition. Whereas 4 lies left of saturation dome telling it is at subcooled condition. 2 lies inside saturation dom and not on the dome, so it has some degree of moisture.

Efficiency here refers, out 100% head given to plant (coal burning, nuclear fission etc) only ~36% is converted to useful work. This can be altered by T1, T2 temperatures but since there is limitations in material available for HE it is restricted to certain value.

SSC refers how much steam required to produce 1kWh of work.

ERRORS FACED

1. Few brackets missed but later added appropriately.
2. In fprintf '%' was used which produced error, later replaced with 'percentage'.