Creating Rankine cycle simulator using MATLAB

Aim – To create a Rankine cycle simulator using MATLAB.

Theory –

The Rankine cycle or Rankine Vapour Cycle is an idealized thermodynamic cycle of a heat engine that converts heat into mechanical work while undergoing phase change. The Rankine cycle is the process widely used by power plants such as coal-fired power plants or nuclear reactors. In this mechanism, a fuel is used to produce heat within a boiler, converting water into steam which then expands through a turbine producing useful work.

There are four processes in the Rankine cycle. The states are identified by numbers in the T–s diagram.

• Process 1–2: Isentropic Compression – The working fluid is pumped from low to high pressure. As the fluid is a liquid at this stage, the pump requires little input energy.
• Process 2–3: Constant Pressure Heat Addition – The high-pressure liquid enters a boiler, where it is heated at constant pressure by an external heat source to become superheated vapour.
• Process 3–4: Isentropic Expansion – The superheated vapour expands through a turbine, generating power. This decreases the temperature and pressure of the vapour, and some condensation may occur.
• Process 4–1: Constant Pressure Heat Rejection – The wet vapour then enters a condenser, where it is condensed at a constant pressure to become a saturated liquid.

Governing Equations –

1. Turbine work –

W_T = h₃ – h₄

2. Pump work –

W_P = h₁ – h₂

3. Net work –

W_net = W_T – W_P

4. Heat added by boiler –

Q_in = h₃ – h₂

5. Back work ratio –

BWR = W_P/W_T

6. Thermal efficiency –

ƞ_th = W_net / Q_in

7. Specific steam consumption –

(SSC) = 3600 / W_net

where,

W = work

Q = heat

h_i = enthalpy at point i

ƞ = efficiency.

Given inputs –

Turbine Inlet temperature = 400° C

Turbine Inlet Pressure = 30 bars

Turbine Outlet Pressure = 0.05 bars = Condenser Inlet Pressure

The link for the standard Xsteam library to get the steam table data for water, is given below. It is to be used in calculation of the points in between the state variables.

Click here to download the XSteam Library

Objective –

1. To write a code which creates a Rankine cycle simulator.
2. To calculate the state points of the Rankine Cycle based on inputs used.
3. To plot the corresponding T-s and h-s plots for the given set of inputs.
4. To calculate Net work output, Thermal efficiency, Back work ratio.

Code –

clear all
close all
clc
 
disp('Simulation of Rankine Cycle')
disp('Process 1-2 is Isentropic compression in Pump')
disp('Process 2-3 is constant pressure heat addition in Boiler')
disp('Process 3-4 is isentropic expansion in Turbine')
disp('Process 4-1 is constant pressure heat rejection in Condenser')
fprintf('\n');
 
% Inputs
P3 = input('\n Enter the Pressure at the Turbine Inlet(in bar): ');
T3 = input('Enter the Temperature at the Turbine Inlet(in Degree Celcius): ');
P4 = input('Enter the Pressure at the Condenser (in bar): ');
 
% Calculating the value of P,T,H,S at each point
% At state 3
H3 = XSteam('h_pT',P3,T3);
S3 = XSteam('s_pT',P3,T3);
H_L3 = XSteam('hL_p',P3);
H_G3 = XSteam('hV_p',P3);
S_L3 = XSteam('sL_p',P3);
S_G3 = XSteam('sV_p',P3);
T_sat3 = XSteam('Tsat_p',P3);
 
% At state 4
S4 = S3;
S_L4 = XSteam('sL_p',P4);
S_G4 = XSteam('sV_p',P4);
% Calculating the dryness fraction
D = (S4-S_L4)/(S_G4-S_L4);
H_L4 = XSteam('hL_p',P4);
H_G4 = XSteam('hV_p',P4);
% Calculating enthalpy at point 4
H4 = H_L4 + D*(H_G4 - H_L4);
T4 = XSteam('T_ph',P4,H4);
T_sat4 = XSteam('Tsat_p',P4);
 
% At state 1
P1 = P4;
T1 = T4;
H1 = H_L4;
S1 = S_L4;
 
% At state 2
P2 = P3;
S2 = S1;
T2 = XSteam('T_ps',P2,S2);
H2 = XSteam('H_ps',P2,S2);
 
 
% Turbine work
W_T = H3 - H4;
 
% Pump work
W_P = H2 - H1;
 
% Net work
W_net = W_T + W_P;
 
% Heat added by boiler
Qin = H3 - H2;
 
% Back work ratio
BWR = W_P/W_T;
 
% Thermal efficiency
Eff = (W_net/(Qin))*100 ;
 
% Specific steam consumption
SSC = 3600/W_T;
 
% Displaying the results
% Values of state 1
disp('Results');
disp('At State point 1');
fprintf('P1 is : %.2f Bar \n',P1);
fprintf('T1 is : %.2f Deg Celcius \n',T1);
fprintf('H1 is : %.2f KJ/kg \n',H1);
fprintf('S1 is : %.2f KJ/kgK \n',S1);
fprintf('\n');
 
% Values of state 2
disp('At State point 2');
fprintf('P2 is : %.2f Bar \n',P2);
fprintf('T2 is : %.2f Deg Celcius \n',T2);
fprintf('H2 is : %.2f KJ/kg \n',H2);
fprintf('S2 is : %.2f KJ/kgK \n',S2);
fprintf('\n');
 
% Values of state 3
disp('At State point 3');
fprintf('P3 is : %.2f Bar \n',P3);
fprintf('T3 is : %.2f Deg Celcius \n',T3);
fprintf('H3 is : %.2f KJ/kg \n',H3);
fprintf('S3 is : %.2f KJ/kgK \n',S3);
fprintf('\n');
 
% Values of state 4
disp('At State point 4');
fprintf('P4 is : %.2f Bar \n',P4);
fprintf('T4 is : %.2f Deg Celcius \n',T4);
fprintf('H4 is : %.2f KJ/kg \n',H4);
fprintf('S4 is : %.2f KJ/kgK \n',S4);
fprintf('\n');
 
% Turbine work
fprintf('Wt is : %.2f KJ/kg \n',W_T);
 
% Pump work
fprintf('Wp is : %.2f KJ/kg \n',W_P);
 
% Net work
fprintf('Wnet is : %.2f KJ/kg \n',W_net);
 
% Heat added by boiler
fprintf('Qin is : %2f KJ/kg \n',Qin);
 
% Back work ratio
fprintf('Back work ratio is : %f \n',BWR);
 
% Thermal efficiency
fprintf('Efficiency is : %.2f Percent \n',Eff);
 
% Specific steam Consumption
fprintf('S.S.C. is : %.2f kg/KWh \n',SSC);
 
 
 
% plotting T-s curve
t = linspace(0,373,374);
figure(1)
hold on
 
for i = 1:length(t)
plot(XSteam('sL_T',t(i)),t(i),'r.')
plot(XSteam('sV_T',t(i)),t(i),'r.')
end
 
plot([S3 S4],[T3 T4],'linewidth',2,'color','b')
text(S3,T3,'3','FontSize',10)
if D>1
T1 = T_sat4;
plot([S4 S_G4 S1 S2],[T4 T_sat4 T1 T2],'linewidth',2,'color','b')
T4 = T_sat4;
text(S4,T4,'4','FontSize',10)
text(S1,T1,'1','FontSize',10)
else
plot([S4 S1 S2],[T4 T1 T2],'linewidth',2,'color','b')
text(S4,T4,'4','FontSize',10)
text(S1,T1,'1','FontSize',10)
end
 
n = linspace(T3,T4,500);
for i = 1:length(n)
plot(XSteam('s_pT',P3,n(i)),n(i),'.','linewidth',2,'color','b')
end
 
plot([S_L3 S_G3],[T_sat3 T_sat3],'linewidth',2,'color','b')
text(S_L3,T_sat3,'2','FontSize',10)
xlabel('Entropy(KJ/kg-K)')
ylabel('Temperature(K)')
title('T-S Diagram')
 
 
%Plotting of H-S Curve
figure(2) 
hold on
for i = 1:length(t)
 plot(XSteam('sL_T',t(i)),XSteam('hL_T',t(i)),'r.')
 plot(XSteam('sV_T',t(i)),XSteam('hV_T',t(i)),'r.')
end
plot([S1 S2 S3 S4 S1],[H1 H2 H3 H4 H1],'linewidth',2,'color','b')
text(S1,H1,'1','FontSize',10)
text(S2,H2,'2','FontSize',10)
text(S3,H3,'3','FontSize',10)
text(S4,H4,'4','FontSize',10)
xlabel('Entropy(KJ/kg-K)')
ylabel('Enthalpy(KJ/kg)')
title('H-S Diagram')

Errors –

1. Many typing errors were found.
2. Multiple results were showing in one line in command window.

Forgot to use '\n' to change the line. After using the required command, error was rectified.

Result –

1. T-s diagram –

2. H-s diagram –

3. Command Window output –

References –

1. https://en.wikipedia.org/wiki/Rankine_cycle

2. http://sounak4u.weebly.com/vapour--combined-power-cycle.html