Rankine cycle simulator

The goal of this challenge is to take in input conditions at the turbine inlet plot the characteristic Rankine cycle T-s and h-s plots using the XSteam function.

A conventional Rankine cycle consists of 4 state points, with 1-2 being isentropic expansion in the turbine, 2-3 being isochoric heat removal in the condenser, 3-4 being isentropic compression in the pump and 4-1 being heat addition in the boiler. For the sake of the challenge, reheating and regeneration are not taken into consideration. Turbine and pump inefficiencies are not taken into consideration for this problem.

x in the figure refers to the dryness of the working fluid, with 0 being wet fluid and 1 being dry vapour. x above 1 refers to superheated steam. Superheated steam is preferred for the turbine inlet. In the region under the curve, it is in liquid-vapour equilibrium.

The code to simulate the rankine cycle is as below:

The code is written such that based on the 3 required input values that will be prompted at the beginning, the program will be able to find the state points and plot the necessary T-s and h-s plots.

clear all
close all
clc

disp(\'        RANKINE CYCLE SIMULATOR\');
disp(\'1-2 is Isentropic Expansion in the turbine\');
disp(\'2-3 is Isochoric heat removal in the condenser\');
disp(\'3-4 is Isentropic compression in the pump\');
disp(\'4-1 is Isochoric heat addition in the boiler\');

%Obtaining the input data
prompt_p = \'Value of pressure at turbine inlet, p1(bar) = \';
p1 = input(prompt_p);
prompt_T = \'Value of temperature at turbine inlet, T1(C) = \';
T1 = input(prompt_T);
prompt = \'Value of temperature at condenser inlet, p2(bar) = \';
p2 = input(prompt);

%At stage 1
h1 = XSteam(\'h_pT\', p1, T1);
s1 = XSteam(\'s_pT\', p1, T1);

%At stage 2
%s2 = sf2-x(sfg2)

s2 = s1;

sf2= XSteam(\'sL_p\', p2);
sg2 = XSteam(\'sV_p\', p2); 
x2 = (s2-sf2)/(sg2-sf2);

hf2 = XSteam(\'hL_p\', p2);
hg2 = XSteam(\'hV_p\', p2);
h2 = hf2 + (x2*(hg2-hf2));

%At stage 3

p3 = p2;

s3= XSteam(\'sL_p\', p3);

h3 = XSteam(\'hL_p\', p3);

%At stage 4

s4 = s3;
p4 = p1;

%Work done by turbine
Wt = h1 - h2;

%Work done to run pump
v3 = XSteam(\'vL_p\', p3);
Wp = v3*(p4 - p3)*100;
h4 = h3 + Wp;

%net work
Wnet = Wt - Wp;

%Heat into the system
Qin = h1 - h4;

%Heat out of the system
Qout = h2 - h3;

%Thermal efficiency
eff = (Wnet/Qin)*100;

%Specific Steam Consumption
%SSC = mass flow rate/power

SSC = (3600/Wnet);

%Calculating temperatures at the state points
T3 = XSteam(\'Tsat_p\', p3);
T2 = T3;
Cp4 = XSteam(\'Cp_ps\', p4, s4);
Cv4 = XSteam(\'Cv_ps\', p4, s4);
k = Cp4/Cv4;
T4 = T3/((p3/p4)^((k-1)/k));

%Calculating temperature, entropy and enthalpy at the vapour and liquid
%line for plotting
T6 = XSteam(\'Tsat_p\', p1);
s5 = XSteam(\'sL_p\', p1);
s6 = XSteam(\'sV_p\', p1);
h5 = XSteam(\'hL_p\', p1);
h6 = XSteam(\'hV_p\', p1);
T5 = T6;

%Plotting the saturation curve
T = linspace(1,375,1000);

for i = 1:length(T)
    
    sf(i) = XSteam(\'sL_T\', T(i));
    sg(i) = XSteam(\'sV_T\', T(i));
    hf(i) = XSteam(\'hL_T\', T(i));
    hg(i) = XSteam(\'hV_T\', T(i));
   
end

%Plotting T-s diagram
figure(1)
hold on
plot(sf, T, \'r\')
plot(sg, T, \'r\')
plot([s1 s2], [T1 T2], \'b\', \'linewidth\', 3)
plot([s2 s3], [T2 T3] , \'b\', \'linewidth\', 3)
plot([s3 s4], [T3 T4], \'b\', \'linewidth\', 3)
plot([s4 s5], [T4 T5], \'b\', \'linewidth\', 3)
plot([s5 s6], [T5 T6], \'b\', \'linewidth\', 3)

sc1 = linspace(s6, s1, 1000);
for l = 1:length(sc1)
    
    Tc1(l) = XSteam(\'T_ps\', p1, sc1(l));  
    
end
plot(sc1, Tc1, \'b\', \'linewidth\', 3)
text(s1+0.1, T1, \'1\')
text(s2+0.1, T2, \'2\')
text(s3, T3-20, \'3\')
text(s4, T4+20, \'4\')
title(\'T-s plot\')
xlabel(\'Entropy\')
ylabel(\'Temperature\')

%Plotting h-s diagram
figure(2)
hold on
plot(sf, hf, \'r\')
plot(sg, hg, \'r\')
plot([s1 s2], [h1 h2], \'b\', \'linewidth\', 3)
plot([s2 s3], [h2 h3], \'b\', \'linewidth\', 3)
plot([s3 s4], [h3 h4], \'b\', \'linewidth\', 3)

hc = linspace(h4, h5, 1000);
for k = 1:length(hc)
    
    sch(k) = XSteam(\'s_ph\', p1, hc(k));  
    
end
plot(sch,hc, \'b\', \'linewidth\', 3)

plot([s5 s6], [h5 h6], \'b\', \'linewidth\', 3)

for m = 1:length(sc1)
    
    hc1(m) = XSteam(\'h_ps\', p1, sc1(m));
    
end
plot(sc1, hc1, \'b\', \'linewidth\', 3)
text(s1+0.1, h1, \'1\')
text(s2+0.1, h2, \'2\')
text(s3, h3-100, \'3\')
text(s4, h4+200, \'4\')
title(\'h-s plot\')
xlabel(\'Entropy\')
ylabel(\'Enthalpy\')

disp(\'   RESULTS\');
fprintf(\'At State Point 1:\\n\');
fprintf(\'p1: %8.3f bar\\n\', p1);
fprintf(\'T1: %8.3f C\\n\', T1);
fprintf(\'h1: %8.3f KJ/kg\\n\', h1);
fprintf(\'s1: %8.3f KJ/kgK\\n\\n\', s1);

fprintf(\'At State Point 2:\\n\');
fprintf(\'p2: %8.3f bar\\n\', p2);
fprintf(\'T2: %8.3f C\\n\', T2);
fprintf(\'h2: %8.3f kJ/kg\\n\', h2);
fprintf(\'s2: %8.3f kJ/kgK\\n\', s2);
fprintf(\'x2: %8.3f\\n\\n\', x2);

fprintf(\'At State Point 3:\\n\');
fprintf(\'p3: %8.3f bar\\n\', p3);
fprintf(\'T3: %8.3f K\\n\', T3);
fprintf(\'h3: %8.3f kJ/kg\\n\', h3);
fprintf(\'s3: %8.3f kJ/kgK\\n\\n\', s3);

fprintf(\'At State Point 4:\\n\');
fprintf(\'p4: %8.3f bar\\n\', p4);
fprintf(\'T4: %8.3f K\\n\', T4);
fprintf(\'h4: %8.3f KJ/kg\\n\', h4);
fprintf(\'s4: %8.3f kJ/kgK\\n\\n\', s4);

fprintf(\'Turbine Work: %8.3f kJ/kg\\n\', Wt);
fprintf(\'Pump Work: %8.3f kJ/kg\\n\', Wp);
fprintf(\'Wnet: %8.3f kJ/kg\\n\', Wnet);
fprintf(\'Efficiency: %8.3f\\n\', eff);
fprintf(\'SSC: %8.3f kJ/kWh\', SSC); 

Initial conditions are taken as 44bar and 538 C at the turbine inlet and 5.83 bar at the turbine outlet for the sake of this simulation.

Using the initial pressure and temperature, the enthalpy and entropy at state point 1 are found using the XSteam function. Since 1-2 is an isentropic process, s1 = s2. In order to determine h2, we require the dryness fraction at state point 2, which we find out using the vapour and liquid saturation values at point 2 and the relation s = sf + x(sg-sf).

The conditions at state point 3 are assumed to lie on the liquid line and their values are obtained from the steam tables using the XSteam function. Finally, at state point 4, the conditions are determined by assuming s3=s4. The work done to run the pump is found out using the relation Wp = v(p4-p3) and Wp = h4-h3 is used to find out the value of h4.

The efficiency of the turbine is calculated next using the formula Efficiency = Wout/Qin.

Using the data obtained at all the state points, the Ts and hs plots are plot.

In order to accurately capture the physical conditions inside the boiler, the temperature, entropy and enthalpy are split into 3 parts: T4, T5 and T6. With T4 being the temperature at the inlet of the boiler and T6 being the temperature at the outlet. The salient enthalpy and entropy are found in this range and plot to obtain the accurate plot.

For the given data, the plots so obtained:

T-s plot:

h-s plot:

Both the characteristic curves are obtained with the profiles as expected.

The results at the state points are printed in the command window.