Rankine Cycle Simulator

Introduction

Rankine cycle is a thermodynamic cycle derived from Carnot vapour power cycle for overcoming its limitations. A Matlab program is to be made to simulate the Rankine Cycle and plot state variables at the start of all of the four processes.

T – S, h – S and P – V representations are as shown below.

Rankine cycle has the following thermodynamic processes.

1 – 2 = Isobaric heat addition (in the boiler)

2 – 3 = Adiabatic expansion (in turbine)

3 – 4 = Isobaric heat release (in the condenser)

4 – 1 = Adiabatic pumping (in pump)

Here the values calculated by the steam table is calculated by a program named XSteam.

Code

1. Firstly three inputs are taken i.e. pressure at points 1 and 2 also the temperature at point 1.
2. After then the state variables for points 1,2,3 and 4 are calculated.
3. Turbine and pump work also heat added and removed is calculated along with efficiency.
4. Now two points are calculated which lies on both boiler and saturated lines.
5. Using these points the T-S diagram is plotted. For the saturation curve, a different matrix is filled.
6. These points are also used to plot the h-S diagram but the curve between these points is calculated using their thermodynamic relation.
7. All the state variables calculated are displayed along with Work, efficiency and heat (addition and released) is displayed.

close all
clear all
clc


fprintf("\n              RANKINE CYCLE SIMULATOR             ");
fprintf("\n\n1-2 is Isentropic Expansion in the Turbine.");
fprintf("\n2-3 is Constant Pressure Heat Rejection by the Condensor.");
fprintf("\n3-4 is Isentropic Compression in the Pump.");
fprintf("\n4-1 is Constant Pressure Heat Addition by the Boiler.\n");

% USER INPUTS
p1 = input("\nEnter the Pressure at the Turbine Inlet (in bar): ");
t1 = input("Enter the Temperature at the Turbine Inlet (in degree Celcius): ");
p2 = input("Enter the Pressure at the Condensor (in bar): ");

% CALCULATIONS


h1 = XSteam('h_pT',p1,t1);
s1 = XSteam('s_pT',p1,t1);


s2 = s1;
sf2 = XSteam('sL_p',p2);
sg2 = XSteam('sV_p',p2);
hf2 = XSteam('hL_p',p2);
hg2 = XSteam('hV_p',p2);
x2 = ((s2-sf2)/(sg2-sf2));
h2 = hf2+(x2*(hg2-hf2));
t2 = XSteam('Tsat_p',p2);

p3 = p2;
t3 = t2;
h3 = XSteam('hL_p',p3);
s3 = XSteam('sL_p',p3);
v3 = XSteam('vL_p',p3);

p4 = p1;
s4 = s3;
cp4 = XSteam('Cp_ps',p3,s4);
cv4 = XSteam('Cv_ps',p3,s4);
g4 = cp4/cv4;
o = (g4-1)/g4;
t4 = t3*((p4/p3)^o);
Wp = v3*(p4-p3)*100;
h4 = h3+Wp;

Wt = h1-h2;
Wnet = Wt-Wp;
Qin = h1-h4;
Qout = h2-h3;

%Efficiency
Ntherm = (Wnet/Qin)*100;

%specific steam consumption
SSC=3600/Wnet;
tf_sat=XSteam('Tsat_p',p1);
sf_sat=XSteam('sL_p',p1);
hf_sat=XSteam('hL_p',p1);
sg_sat=XSteam('sV_p',p1);
hg_sat=XSteam('hV_p',p1);

figure(1)
%plotting T-S diagram
hold on;

%for plottin saturation curve
tsat1=linspace(0,374,1000);

for i=1:length(tsat1)
  ss1(i)=XSteam('sL_T',tsat1(i));
end

for j=1:length(tsat1)
  ss2(j)=XSteam('sV_T',tsat1(j));
end


plot(ss1,tsat1,'b','lineWidth',2)

plot(ss2,tsat1,'m','lineWidth',2)

%for saturation temp at p1
temp=XSteam('Tsat_p',p1);
t=linspace(t4,temp,30);
for i=1:length(t)
    S(i)=XSteam('s_pt',p1,t(i));
    hold on
   
end

t0=linspace(temp,t1,300);
for i=1:length(t0)
    S0(i)=XSteam('s_pt',p1,t0(i));
end
    
%for plotting turbine work
plot([s1 s2],[t1 t2],'g','lineWidth',1.5)
%for heat rejection in condenser
plot([s2 s3],[t2 t3],'c','lineWidth',1.5)
%for pump work
plot([s3 s4],[t3 t4],'r','lineWidth',1.5)
  
%plotting curved line(constant pressure line in wet region) from outlet of
%the pump to the point from where phase change starts
plot([s4 sf_sat],[t4 tf_sat],'y','lineWidth',1.5)

%plotting the phase change region
plot([sf_sat sg_sat],[tf_sat tf_sat],'y','lineWidth',1.5)
 
%plotting the curved constant pressure line from where the phase change
%ends to the superheated point
plot(S0,t0,'y','lineWidth',1.5)

legend('Saturation Curve L','Saturation Curve V','Turbine','condenser','pump','boiler','location','northwest')
title('T-S Diagram')
xlabel('Entropy')
ylabel('Temperature')
grid on;
hold off;

%plotting H-S diagram
figure(2)
hold on


%for plotting the saturation curve
tsat1=linspace(0,374,1000);

for i=1:length(tsat1)
  hh1(i)=XSteam('hL_T',tsat1(i));
end

for j=1:length(tsat1)
  hh2(j)=XSteam('hV_T',tsat1(j));
end




plot(ss1,hh1,'b','lineWidth',1.5)
plot(ss2,hh2,'r','lineWidth',1.5)

%for plotting tubine

plot([s1 s2],[h1 h2],'g','lineWidth',1.5)

%for condenser
plot([s2 s3],[h2 h3],'k','lineWidth',1.5)

%for plotting pump
plot([s3 s4],[h3 h4],'r','lineWidth',1.5)
s=linspace(s3,s1,1000);
for i=1:length(s)
  h(i)=XSteam('h_ps',p1,s(i));
end
% for plotting boiler%
plot([s4 s1],[h4 h1],'y','lineWidth',1.5)
plot(s,h,'y','lineWidth',1.5)

legend('Saturation Curve L','Saturation Curve V','Turbine','condenser','pump','boiler','location','northwest')
title('H-S Diagram')
xlabel('Entropy')
ylabel('Enthalpy')
grid on;
hold off;

%outputs
%at state 1
fprintf('\n\nAt State Point 1:')
fprintf('\nP1 is : %f Bar',p1)
fprintf('\nT1 is : %f Deg Celcius',t1)
fprintf('\nH1 is : %f KJ/Kg',h1)
fprintf('\nS1 is : %f KJ/KgK',s1)

%at state 2
fprintf('\n\nAt State Point 2:')
fprintf('\nP2 is : %f Bar',p2)
fprintf('\nT2 is : %f Deg Celcius',t2)
fprintf('\nH2 is : %f KJ/Kg',h2)
fprintf('\nS2 is : %f KJ/KgK',s2)
fprintf('\nX is  : %f %',x2)

%at state 3
fprintf('\n\nAt State Point 3:')
fprintf('\nP3 is : %f Bar',p3)
fprintf('\nT3 is : %f Deg Celcius',t3)
fprintf('\nH3 is : %f KJ/Kg',h3)
fprintf('\nS3 is : %f KJ/KgK',s3)

%state 4
fprintf('\n\nAt State Point 4:')
fprintf('\nP4 is : %f Bar',p4)
fprintf('\nT4 is : %f Deg Celcius',t4)
fprintf('\nH4 is : %f KJ/Kg',h4)
fprintf('\nS4 is : %f KJ/KgK',s4)

%calculations
fprintf('\n\nWork Done By Turbine : %f kJ/kg',Wt)
fprintf('\nWork Done By Pump : %f kJ/kg',Wp)
fprintf('\nHeat Added in The Boiler: %f KJ/kg',Qin)
fprintf('\nHeat Rejected in The Condenser: %f KJ/kg',Qout)
fprintf('\nThermal Efficiency is : %f %',Ntherm*100)
fprintf('\nSpecific Steam Consumption : %f KJ/kWh\n',SSC)

Inputs

`P1=44 "bar"`

`P2=5.83 "bar"`

`T1=538 ^oC`

Output