BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

Rankine cycle Simulator

SSushanto Sarkar21 May 2020
MATLAB

PROJECT 2- RANKINE CYCLE SIMULATOR USING MATLAB CODING

Objective: Make a Rankine Cycle simulator on MATLAB

  • Understand the basic concept of the Rankine cycle.
  • Calculate the state points based on user inputs.
  • Plot T-S & H-S graphs as per the user inputs.

Introduction :

RANKINE CYCLE: THE IDEAL CYCLE FOR VAPOUR POWER CYCLE

Many of the impracticalities associated with the Carnot cycle can be eliminated by superheating the steam inside the boiler and condensing it completely in a condenser. The cycle that results is the ideal Rankine cycle, this cycle does not involve any internal reversibility.

Rankine cycle is the practical working cycle for the thermal power plant & consists of 4 processes.

Figure.1 A simple steam plant process 1-2-3-4

Process 1-2: Reversible adiabatic/isentropic expansion process of high enthalpy steam from the boiler to condenser pressure in the turbine.

Process 2-3: Constant pressure heat rejection process in the condenser to get saturated liquid.

Process 3-4: Reversible adiabatic/isentropic compression of liquid from the condenser to boiler pressure through the boiler feed pump.

Process 4-1: Constant pressure heat addition process to convert the liquid into superheated steam in a boiler.

Figure.2 Rankine cycle on a) P-V, b) T-S, and c) H-S diagram

For the purpose of analysis the ranking cycle is assumed to be carried out in a steady flow operation. applying the steady flow energy equation to each of the processes on the basis of a unit mass of fluid and neglecting change in K.E and P.E, the work and heat quantities can be evaluated in terms of fluid properties.

Applying the Steady flow energy equation (SFEE) across the turbine, condenser, boiler

Work done by the turbine `W_T= h_1-h_2`

work supplied to the feed pump `W_p= h_4-h_3`

Heat supplied by the boiler `Q_1= h_1-h_4`

Heat rejected by the condenser `Q_2= h_2- h_3`

The efficiency of a Rankine cycle is:

`eta_R=(W_(n et))/ Q_s=( W_T- W_P)/ Q_s= ((h_1-h_2)-(h_4- h_3))/(h_1-h_2)`

Back Work Ratio(BWR)

it is the part of work output from a turbine work is used to drive the feed pump which required work input, therefore, it is the ratio of pump work input /turbine work output

`BWR=W_c/W_T`

Pump Work

The pump handles liquid water which is incompressible therefore for isentropic compression `dh=vdp= h_4- h_3= W_(pump)`

The pump work is quite small compared to the turbine work and some times neglected as h4=h3. Therefore efforts are made to keep the specific volume of the fluid as small as possible during the compression process to minimize the work input and as large as possible during an expansion process to maximize the work output.

Steam property table

XSteam.m file provides accurate steam and water properties from 0 - 1000 bar and from 0 - 2000 deg C

First the wanted property then a _ then the wanted input properties.

Example.

T_ph is temperature as a function of pressure and enthalpy.

h_pT Enthalpy as a function of pressure and temperature.

sV_p Saturated vapour entropy
sL_p Saturated liquid entropy

hV_p Saturated vapour enthalpy
hL_p Saturated liquid enthalpy

Tsat_p Saturation temperature

XSteam takes 2 or 3 arguments. The first argument must always be the steam table function you want to use. The other arguments are the inputs to that function.

Example:

XSteam('h_pt',1,20) Returns the enthalpy of water at 1 bar and 20 degC
XSteam('TSat_p',1) Returns the saturation temperature of water at 1 bar

Methodology:

  1. Take inputs from the user
  2. Using the steam table data and user inputs calculation of other state points properties.
  3. Acquiring the data from the steam table for plotting the saturation curve.
  4. Plotting the state point on saturation curves T-Sor H-S.
  5. Display the data of the state points along with work relevant to a turbine or pump.

MATLAB CODES FOR RANKINE SIMULATOR

                        %     PROJECT 2-     %
                        % RANKINE SIMULATOR  %




close all
clear all
clc

disp('                       RANKINE CYCLE SIMULATOR                  ');
fprintf('\n');
disp('1-2 Isentropic Expansion In The Turbine');
disp('2-3 Constant Pressure Heat Rejection By The Condenser');
disp('3-4 Isentropic compression of pump');
disp('4-1 Constant Pressure Heat Addition'); 
fprintf('\n');

% user input for the cycle
t1= input('Enter The Turbine Input Temperature in Degree Celcius :  ');
p1= input('Enter The Turbine Input Pressure in Bar :  ');
p2= input('Enter the turbine exit Pressure in Bar:  ');
fprintf('\n');
gamma=1.4;

% ploting the saturated curve
t=linspace(0,373,200);

for i=1:length(t)

    sf_sat(i)=XSteam('sL_T',t(i));
    sg_sat(i)=XSteam('sV_T',t(i));
    hf_sat(i)=XSteam('hL_T',t(i));
    hg_sat(i)=XSteam('hV_T',t(i));
end
% Temp & pressure vaulue at critical point 
t_crt= 373.946;
p_crt=220.064;
s_crt= XSteam('s_pT',p_crt,t_crt);
h_crt= XSteam('h_pT',p_crt,t_crt);

% State point properties of Rankine cycle 
% State 1 (Turbine inlet) 
h1 =XSteam('h_pt',p1,t1);
s1=XSteam('s_pT',p1,t1);

%process 1-2 Isentropic Expansion In The Turbine s1=s2
% For Isentropic process PV^gamma=constant P1/P2= (V2/V1)^gamma
% with PV=NRT, [ T2/T1=(p2/p1)^(gamma-1)/gamma) where T1,T2 is in Kelvin
 %t2= ((t1+273)*(p2/p1)^((gamma-1)/gamma))-273;

t2= XSteam('Tsat_p',p2);
% State 2(Condensor Inlet)
s2=s1;

% Dyrness fraction 'x' calculation 
sf2_sat=XSteam('sL_p',p2);
sg2_sat=XSteam('sV_p',p2);
hf2_sat=XSteam('hL_p',p2);
hg2_sat=XSteam('hV_p',p2);

% s2= sf_sat_p2+x*(sv_sat_p2-sf_sat_p2)where x is the dryness fraction 
x=(s2-sf2_sat)/(sg2_sat-sf2_sat);
h2 = hf2_sat+x*(hg2_sat-hf2_sat);
% t2=XSteam('T_ph', p2, h2)

%State3(Pump Inlet)
p3=p2;
s3=sf2_sat;
t3=t2;
h3=XSteam('hL_p',p2);

%state 4(boiler inlet)
s4=s3;
p4=p1;
%t4= ((t3+273)*(p3/p4)^((1-gamma)/gamma))-273;
t4= XSteam('T_ps',p4,s4);
h4=XSteam('h_pT', p4,t4);
tsat1 = XSteam('Tsat_p',p1);
sf1_sat=XSteam('sL_p',p1);
sg1_sat=XSteam('sV_p',p1);
hf1_sat=XSteam('hL_p',p1);
hg1_sat=XSteam('hV_p',p1);
% plot graphs 
% T-S plot for rankine cycle 
hold on
figure(1)
plot (sf_sat, t, 'color','b');
plot (sg_sat, t,'color','r');
plot([s1 s2],[t1 t2],'marker', '*','linewidth',2);
plot([s2 s3], [t2 t3],'marker', '*','linewidth',2);
plot([s3 s4],[t3 t4],'marker', '*','linewidth',2);
plot([s4 sf1_sat],[t4 tsat1],'marker', '*','linewidth',2,'color','c');
plot([sf1_sat sg1_sat],[tsat1 tsat1],'marker', '*','linewidth',2,'color','c');
plot([sg1_sat s1],[tsat1 t1],'marker', '*','linewidth',2,'color','c');

xlabel('Entropy-(S) /(KJ/kgK)');
ylabel(' Temperature-(T) /( degree C)');
title('T-S plot for Rankine Cycle)');
legend('Liquid Saturation curve',' Vapour Saturation curve','1-2 Isentropic Expansion', '2-3 Constant Pr. Heat Rejection','3-4 Isentropic Compression ','4-1 Constant Pr. Heat Addition '   )

% H-S plot for rankine cycle
figure(2)
hold on 
plot ( sf_sat,hf_sat, 'color','b');
plot (sg_sat,hg_sat,'color','r');
plot([s1 s2],[h1 h2],'marker', '*','linewidth',2)
plot([s2 s3], [h2 h3],'marker', '*','linewidth',2)
plot([s3 s4],[h3 h4],'marker', '*','linewidth',2)
plot([s4 s1],[h4 h1],'marker', '*','linewidth',2)
xlabel('Entropy-(S) /(KJ/kgK)');
ylabel(' Enthalpy -(H) /(KJ/kg');
title('H-S plot for Rankine Cycle');
legend('Liquid Saturation curve',' Vapour Saturation curve','1-2 Isentropic Expansion', '2-3 Constant Pr. Heat Rejection','3-4 Isentropic Compression ','4-1 Constant Pr. Heat Addition ')

%work calculation 

wt= h1-h2;
wp= h4-h3;
wnet= wt-wp;
ntherm= wnet/(h1-h4)*100;
back_work= wp/wt;
ssc=3600/wnet;

% Simulator Display

disp('                                RESULT                           ');
fprintf('\n');
disp('AT STATE POINT 1');
P1=sprintf('P1 is :  %.4f Bar', p1);
disp(P1);
T1=sprintf('T1 is :  %.4f Degree celcius', t1);
disp(T1);
H1=sprintf('H1 is :  %.4f KJ/kg', h1);
disp(H1);
S1=sprintf('S1 is :  %.4f KJ/kgK', s1);
disp(S1);

fprintf('\n');

disp('AT STATE POINT 2');
P2=sprintf('P2 is :  %.4f Bar', p2);
disp(P2);
T2=sprintf('T2 is :  %.4f Degree celcius', t2);
disp(T2);
H2=sprintf('H2 is :  %.4f KJ/kg', h2);
disp(H2);
S2=sprintf('S2 is :  %.4f KJ/kgK', s2);
disp(S2);

fprintf('\n');

disp('AT STATE POINT 3');
P3=sprintf('P1 is :  %.4f Bar', p3);
disp(P3);
T3=sprintf('T3 is :  %.4f Degree celcius', t3);
disp(T3);
H3=sprintf('H3 is :  %.4f KJ/kg', h3);
disp(H3);
S3=sprintf('S3 is :  %.4f KJ/kgK', s3);
disp(S3);

fprintf('\n');

disp('AT STATE POINT 4');
P4=sprintf('P1 is :  %.4f Bar', p4);
disp(P4);
T4=sprintf('T4 is :  %.4f Degree celcius', t4);
disp(T4);
H4=sprintf('H4 is :  %.4f KJ/kg', h4);
disp(H4);
S4=sprintf('S4 is :  %.4f KJ/kgK', s4);
disp(S4);

fprintf('\n');

WT=sprintf('Wt is :  %.4f KJ/kg', wt);
disp(WT);
WP=sprintf('Wp is :  %.4f KJ/kg', wp);
disp(WP);
WNET=sprintf('Wnet is :  %.4f KJ/kg', wnet);
disp(WNET);
NTHERM=sprintf('Ntherm is :  %.4f percent',ntherm);
disp(NTHERM);
SSC=sprintf('S.S.C is :  %.4f kg/kWh', ssc);
disp(SSC);

Code Explanation

1) First, we start by displaying the process details.

disp('                       RANKINE CYCLE SIMULATOR                  ');
fprintf('\n');
disp('1-2 Isentropic Expansion In The Turbine');
disp('2-3 Constant Pressure Heat Rejection By The Condenser');
disp('3-4 Isentropic compression of pump');
disp('4-1 Constant Pressure Heat Addition'); 
fprintf('\n');

2) Then we will take inputs from the user 

% user input for the cycle
t1= input('Enter The Turbine Input Temperature in Degree Celcius :  ');
p1= input('Enter The Turbine Input Pressure in Bar :  ');
p2= input('Enter the turbine exit Pressure in Bar:  ');
fprintf('\n');

3) Saturation curve data acquiring and curve plotting using the steam table function 

t=linspace(0,373,200);

for i=1:length(t)

    sf_sat(i)=XSteam('sL_T',t(i));
    sg_sat(i)=XSteam('sV_T',t(i));
    hf_sat(i)=XSteam('hL_T',t(i));
    hg_sat(i)=XSteam('hV_T',t(i));
end
% Temp & pressure vaulue at critical point 
t_crt= 373.946;
p_crt=220.064;
s_crt= XSteam('s_pT',p_crt,t_crt);
h_crt= XSteam('h_pT',p_crt,t_crt);

4) State point properties calculation by the steam table function 

% State point properties of Rankine cycle 
% State 1 (Turbine inlet) 
h1 =XSteam('h_pt',p1,t1);
s1=XSteam('s_pT',p1,t1);

%process 1-2 Isentropic Expansion In The Turbine s1=s2
% For Isentropic process PV^gamma=constant P1/P2= (V2/V1)^gamma
% with PV=NRT, [ T2/T1=(p2/p1)^(gamma-1)/gamma) where T1,T2 is in Kelvin
 %t2= ((t1+273)*(p2/p1)^((gamma-1)/gamma))-273;

t2= XSteam('Tsat_p',p2);
% State 2(Condensor Inlet)
s2=s1;

% Dyrness fraction 'x' calculation 
sf2_sat=XSteam('sL_p',p2);
sg2_sat=XSteam('sV_p',p2);
hf2_sat=XSteam('hL_p',p2);
hg2_sat=XSteam('hV_p',p2);

% s2= sf_sat_p2+x*(sv_sat_p2-sf_sat_p2)where x is the dryness fraction 
x=(s2-sf2_sat)/(sg2_sat-sf2_sat);
h2 = hf2_sat+x*(hg2_sat-hf2_sat);
% t2=XSteam('T_ph', p2, h2)

%State3(Pump Inlet)
p3=p2;
s3=sf2_sat;
t3=t2;
h3=XSteam('hL_p',p2);

%state 4(boiler inlet)
s4=s3;
p4=p1;
%t4= ((t3+273)*(p3/p4)^((1-gamma)/gamma))-273;
t4= XSteam('T_ps',p4,s4);
h4=XSteam('h_pT', p4,t4);
tsat1 = XSteam('Tsat_p',p1);
sf1_sat=XSteam('sL_p',p1);
sg1_sat=XSteam('sV_p',p1);
hf1_sat=XSteam('hL_p',p1);
hg1_sat=XSteam('hV_p',p1);

5) T-S and H-S curve plotting 

% plot graphs 
% T-S plot for rankine cycle 
hold on
figure(1)
plot (sf_sat, t, 'color','b');
plot (sg_sat, t,'color','r');
plot([s1 s2],[t1 t2],'marker', '*','linewidth',2);
plot([s2 s3], [t2 t3],'marker', '*','linewidth',2);
plot([s3 s4],[t3 t4],'marker', '*','linewidth',2);
plot([s4 sf1_sat],[t4 tsat1],'marker', '*','linewidth',2,'color','c');
plot([sf1_sat sg1_sat],[tsat1 tsat1],'marker', '*','linewidth',2,'color','c');
plot([sg1_sat s1],[tsat1 t1],'marker', '*','linewidth',2,'color','c');

xlabel('Entropy-(S) /(KJ/kgK)');
ylabel(' Temperature-(T) /( degree C)');
title('T-S plot for Rankine Cycle)');
legend('Liquid Saturation curve',' Vapour Saturation curve','1-2 Isentropic Expansion', '2-3 Constant Pr. Heat Rejection','3-4 Isentropic Compression ','4-1 Constant Pr. Heat Addition '   )

% H-S plot for rankine cycle
figure(2)
hold on 
plot ( sf_sat,hf_sat, 'color','b');
plot (sg_sat,hg_sat,'color','r');
plot([s1 s2],[h1 h2],'marker', '*','linewidth',2)
plot([s2 s3], [h2 h3],'marker', '*','linewidth',2)
plot([s3 s4],[h3 h4],'marker', '*','linewidth',2)
plot([s4 s1],[h4 h1],'marker', '*','linewidth',2)
xlabel('Entropy-(S) /(KJ/kgK)');
ylabel(' Enthalpy -(H) /(KJ/kg');
title('H-S plot for Rankine Cycle');
legend('Liquid Saturation curve',' Vapour Saturation curve','1-2 Isentropic Expansion', '2-3 Constant Pr. Heat Rejection','3-4 Isentropic Compression ','4-1 Constant Pr. Heat Addition ')

6) Work and efficiency calculation 

%work calculation 

wt= h1-h2;
wp= h4-h3;
wnet= wt-wp;
ntherm= wnet/(h1-h4)*100;
back_work= wp/wt;
ssc=3600/wnet;

7) Command window display codes 

% Simulator Display

disp('                                RESULT                           ');
fprintf('\n');
disp('AT STATE POINT 1');
P1=sprintf('P1 is :  %.4f Bar', p1);
disp(P1);
T1=sprintf('T1 is :  %.4f Degree celcius', t1);
disp(T1);
H1=sprintf('H1 is :  %.4f KJ/kg', h1);
disp(H1);
S1=sprintf('S1 is :  %.4f KJ/kgK', s1);
disp(S1);

fprintf('\n');

disp('AT STATE POINT 2');
P2=sprintf('P2 is :  %.4f Bar', p2);
disp(P2);
T2=sprintf('T2 is :  %.4f Degree celcius', t2);
disp(T2);
H2=sprintf('H2 is :  %.4f KJ/kg', h2);
disp(H2);
S2=sprintf('S2 is :  %.4f KJ/kgK', s2);
disp(S2);

fprintf('\n');

disp('AT STATE POINT 3');
P3=sprintf('P1 is :  %.4f Bar', p3);
disp(P3);
T3=sprintf('T3 is :  %.4f Degree celcius', t3);
disp(T3);
H3=sprintf('H3 is :  %.4f KJ/kg', h3);
disp(H3);
S3=sprintf('S3 is :  %.4f KJ/kgK', s3);
disp(S3);

fprintf('\n');

disp('AT STATE POINT 4');
P4=sprintf('P1 is :  %.4f Bar', p4);
disp(P4);
T4=sprintf('T4 is :  %.4f Degree celcius', t4);
disp(T4);
H4=sprintf('H4 is :  %.4f KJ/kg', h4);
disp(H4);
S4=sprintf('S4 is :  %.4f KJ/kgK', s4);
disp(S4);

fprintf('\n');

WT=sprintf('Wt is :  %.4f KJ/kg', wt);
disp(WT);
WP=sprintf('Wp is :  %.4f KJ/kg', wp);
disp(WP);
WNET=sprintf('Wnet is :  %.4f KJ/kg', wnet);
disp(WNET);
NTHERM=sprintf('Ntherm is :  %.4f percent',ntherm);
disp(NTHERM);
SSC=sprintf('S.S.C is :  %.4f kg/kWh', ssc);
disp(SSC);

% State point number ploting
T=[t1 t2 t3 t4];
S=[s1 s2 s3 s4];
H=[h1 h2 h3 h4];
delta =0.05;
figure(1)
sp={'1','2','3','4'};
text(S+delta, T+delta,sp,'color','r','fontsize',14)

Command window display

Results

Below we have a T-S plot for a Rankine cycle with its user-defined inputs. The process 1-2-3-4 can be very well understood from the plot, as the turbine input(state 1) takes place in the superheated region and due to an isentropic process, we get a straight vertical line from pressure state 1 to 2and hence we can extract turbine work WT from process 1-2. As it can be seen that the turbine outlet state 2 is into a Wet region i.e liquid-vapour region at pressure p2 i.e 0.05 bar. so from the state 2 to 3, we need a condenser that converts vapour-liquid state to liquid state so that pump can isentropically increase the pressure from p3 to p4( pump only work for the liquid state) and this conversion can only be achieved by rejecting heat (process 2-3) till saturated liquid state. The pump work(3-4) is very small as compared to the turbine work and it is desired so, therefore, it not visible in the plot. Also, the specific volume is kept as small as possible to extract maximum Wnet(Wt-Wp). Final process 4-1 i.e heat addition process at constant pressure till the supersaturated state of the turbine inlet. Similarly, in the case of the H-S plot can be seen where instead of temperature on Y-Axis enthalpy is used.

FIGURE.3 T-S PLOT FOR RANKINE CYCLE

FIGURE.4 H-S PLOT FOR RANKINE CYCLE

Conclusion: From the above-explained Matlab codes for the Rankine cycle simulator it can be concluded that by just taking the input from the user for turbine inlet temperature(T1), pressure(P1), and outlet pressure (P2) we can simulate the complete cycle. This program can be used for 'nth' no. of cycles with any user-defined parameter. so we can say that the program is very flexible and reduces time in state properties calculation for different cycles. The Results are much convincing and data are much relevant.

References:

1) Engineering Thermodynamics by P.K Nag.

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.