Processing of Rankine Cycle Simulation in MATLAB.

Aim: To calculate all the state variables of a rankine cycle and plot t-s and h-s diagram.

Introduction:

Rankine cycle: Rankine cycle is a thermodynamic cycle which is used to evaluate the performance of a steam turbine systems. Previously, it was used to study the performance of reciprocating steam engines. It is an idealized thermodynamic cycle of a heat engine that converts heat into mechanical work while undergoing phase change. 

As we know, steam turbine system has four basic components. Those are boiler, steam turbine, condenser/heat exchanger and feed water pump.

Following is the simple representation of steam turbine power plant working and t-s diagram of the cycle.

As we know that a simple rankine cycle consists of four thermodynamic processes. 

Process 1-2: Isentropic expansion in the turbine.

Process 2-3: Constant pressure heat rejection by the condenser. 

Process 3-4: Isentropic compression in the feed water pump.

Process 4-1: Constant pressure heat addition by the boiler.

Objectives:

1. To calculate all the state points and performance evaluating  parameters.

2. To plot t-s and h-s diagram of rankine cycle

 Formulae to evaluate performance of the rankine cycle:

1. Isentropic expansion:

T2/T1 = (P2/P1)^1-(1/gamma)

2. Workdone by the pump(W_p):

W_p = Specific volume*(P1-p2)*100

3. Workdone by the turbine(W_t):

W_t = h1 - h2

4. Net workdone(W_net):

W_net = W_t - W_p

5. Back work ratio(W_back):

W_back = W_p/W_t

6. Heat supplied in the boiler(Q_in):

Q_in = h1 - h4

7. Heat rejected in the condenser(Q_out):

Q_out = h2 - h3

8. Thermal efficiency of a rankine cycle (therm_eff) is given by,

therm_eff = W_net/Q_in

9. Specific steam consumption (SSC):

SSC = 3600/W_net

10. Dryness fraction (X):

X = (S - Sf)/Sfg 

Procedure and Description of the MATLAB code:

• First of all, we have to make a rough representation of our rankine cycle simulation and should refer normal t-s and h-s diagram of the rankine cycle to avoid unnecessary calculations. Then start the MATLAB program with "clear all", "close all"and "clc".
• Initially print the header and sub-header for our rankine cycle calculations. Name all the step-wise thermodynamic processes of a rankine cycle and try to run the program. Now we have our header and process structure ready.
• Secondly, insert the input conditions of our cycle on which further calculations are evaluated upon.

Note: We have a function file named as "XSteam", which contains sub-functions i.e.(formulae) required to calculate various thermodynamic properties of steam and water at various stages.

• For state 1, the dry superheated steam is at the inlet of steam turbine and certainly has some inlet pressure and temperature. So we have the given inlet pressure and temperature of steam, but we need to calculate enthalpy and entropy of the steam at state 1 condition. So will call the "XSteam" function and will calculate enthalpy and entropy at state 1.
• Similarly, we will calculate all the required variables (Temperature, pressure, enthalpy, entropy and dryness fraction) at each state respectively with appropriate comments to avoid baffling of the calculations.
• After calculating all the variables at various states, we have to calculate the work done at the turbine and pump as well as resulted net work and back work ratio.
• Similarly, we will be calculating thermal efficiency and specific steam consumption with the use of above variables calculated.
• However, to plot t-s and h-s diagram of the rankine cycle, we will need saturation curve in order to comprehend physical state of working fluid at various stages. So we will take standard saturation temperature range for water as 0 to 374 Deg. celsius by linspace command with 1000 values in between.
• So for saturation liquid and vapour curve, again we will call the "XSteam function" to calculate the enthalpy and entropy at the estimated temperature range. After resulting these, we will be able to plot the saturation curve in t-s and h-s diagram.
• Now we have to plot all those calculated variables and state points with different colors in t-s and h-s diagram by using "plot' command. Then we have to define and label the axis as well as provide a title and legends to the plot.
• Finally, we will have a comprehensive visualization of calculated data in the form of graph describing all the processes, state points and changes in their physical state accordingly.
• Additionally, we need to print all those calculations in the command window, because plots will provide visuallization and we must have detailed report of the results and the input conditions. So will print each variable with its state significance in the command window using "fprintf" command.

 MATLAB CODE:

clear all
close all
clc

fprintf('                                        RANKINE CYCLE SIMULATORrn')
fprintf(' 1-2 : Isentropic expansion in the turbinern')
fprintf(' 2-3 : Constant pressure heat rejection by the condenserrn')
fprintf(' 3-4 : Isentropic compression in the pumprn')
fprintf(' 4-1 : Constant pressure heat addition by the boilerrn')

% Inputs
p1 = input('Enter the pressure at 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) = ');


% Variables at the inlet of steam turbine.
% State 1 calculations
h1 = XSteam('h_pt',p1,t1);
s1 = XSteam('s_pt',p1,t1);
tsat = XSteam('tsat_p',p1);

% Variables at the inlet of the condenser.

% State 2 calculations
s2 = s1; %(Since it is isentropic process)

% Enthalpy in superheated region
sg2 = XSteam('sv_p',p2);

% Enthaply in subcooled region
sf2 = XSteam('sl_p',p2);

% Enthalpy in wet region.
sfg2 = sg2 - sf2;

% To calculate dryness fraction 
X2 = (s2 - sf2)/sfg2;

% To calculate the entalpy
hf2 = XSteam('hl_p',p2);
hg2 = XSteam('hv_p',p2);
hfg2 = hg2 - hf2;
h2 = hf2 + (X2*hfg2);
% To calculate the temperature.
t2 = XSteam('tsat_p',p2);

% Variables at the inlet of feed water pump.

% State 3 calculations
% To calculate the temperature.
t3 = t2; % (As the temperature is constant from process 2-3).

% To calculate the enthalpy.
p3 = p2; % (As it is constant pressure heat rejection in the condenser).

% To calculate the enthalpy.
h3 = XSteam('hl_p',p3);

% To calculate the entropy.
s3 = XSteam('sl_p',p3);

% To calculate the specific volume
sp_v = XSteam('vl_p',p3);

% To calculate the workdone by the pump.
W_p = sp_v*(p1 - p2)*100;

% To calculate the workdone by the turbine.
W_t = h1-h2;

% Variables at the inlet of boiler.

% State 4 calculations
% To calculate the pressure
p4 = p1; % (As there is constant pressure heat addition in the boiler.)

% To calculate the enthalpy
% (W_p = h4 - h3)
h4 = W_p + h3;

% To calculate the entropy
s4 = s3;

% To calculate the heat capacity at constant pressure
Cp = XSteam('cp_ps',p4,s4);
% To calculate the heat capacity at constant volume.
Cv = XSteam('cv_ps',p4,s4);
% To calculate the heat capacity ratio
gamma = Cp/Cv;

% To calculate the temperature
t4 = t3*((p4/p3)^(gamma - 1)/gamma);


% State 5 calculations (Sensible heating)
t5 = tsat;
h5 = XSteam('hL_p',p1);
s5 = XSteam('sL_p',p1);

% State 6 calculations (Latent heating)
t6 = t1;
h6 = XSteam('hV_p',p1);
s6 = XSteam('sV_p',p1);

% To calculate the net work produced in the cycle.
W_net = W_t - W_p;

% To calculate the back work ratio
W_back = W_p/W_t;

% To calculate the heat supplied by the boiler.
Q_in = h1 - h4;

% To calculate the heat rejected by the condenser.
Q_out = h2 - h3;

% To calculate the thermal efficiency of the cycle.
therm_eff = (W_net/Q_in);

% To calculate the specific steam consumption.
SSC = 3600/W_net;

% Calculations for T-S diagram
% To calculate saturation curve points.
t = linspace(0,374,1000);
for i = 1:length(t)
sf(i) = XSteam('sL_t',t(i));
hf(i)= XSteam('hL_t',t(i));
sg(i)= XSteam('sV_t',t(i));
hg(i)= XSteam('hV_t',t(i));
end

% Plotting T-S diagram
figure(1)
hold on
grid on
% Plotting the saturation curve
% Saturated liquid curve
plot(sf,t,'b');
% Saturated vapour curve
plot(sg,t,'r');

% Plotting work obtained in the turbine
plot([s1 s2],[t1 t2],'color','g');
text(s1,t1,'state 1')
text(s2,t2,'state 2')
% Plotting heat rejection in the condenser
plot([s2 s3],[t2 t3],'color','c');
text(s3,t3,'state 3')
% Plotting workdone by the pump.
plot([s3 s4],[t3 t4],'color','k');
text(s4,t4,'state 4')
% Plotting the sensible heating process (Increament in temperature but no
% phase change)
plot([s4 s5],[t4 t5],'color','r');
text(s5,t5,'state 5')
% Plotting latent heating process (phase changing process)
plot([s5 s6],[t5 tsat],'color','m');
text(s6,tsat,'state 6')
% Plotting the superheating process of steam.
plot([s6 s1],[tsat,t1],'color','c');
% Defining axis and labelling
axis([0 8 10 500])
xlabel('Enthalpy (KJ/Kg)')
ylabel('Temperature (Degree celsius)')
title ('Temperature vs Enthalpy (t-s) diagram')
legend('Saturated liquid curve','Saturated vapour curve','Turbine work', 'Const. pressure heat rejection in the condenser','Pump work','Sensible heating in the boiler','Latent heating in the boiler','Superheating in the boiler','location','northwest')


% Calculations for H-S diagram.
% Temperature for obtaining saturation curve of water.

% Plotting H-S diagram
figure(2)
hold on
grid on
% Plotting saturation curve
% Saturated liquid curve
plot(sf,hf,'b');
% Saturated vapour curve
plot(sg,hg,'r');

% Plotting all the state variables.
% Process 1-2
plot([s1 s2],[h1,h2],'g');
text(s1,h1,'state 1')
text(s2,h2,'state 2')
% Process 2-3
plot([s2 s3],[h2 h3],'y');
text(s3,h3,'state 3')
% Process 3-4
plot([s3 s4],[h3 h4],'m');
text(s4,h4,'state 4')
% Process 4-1
plot([s4 s1],[h4 h1],'c');
% Defining axis and labelling
axis([0 8 0 4000])
xlabel('Entropy (KJ/Kg K)')
ylabel('Enthalpy (KJ/Kg)')
title ('Enthalpy vs Entropy (h-s) diagram')
legend('Saturated liquid curve','saturated vapour curve','Isentropic Expansion','Heat Rejection at p = constant','Isentropic Compression','Heat Addition at p = constant','location','northwest')

% To print output
% State 1
fprintf('At state point 1:rn')
fprintf('Pressure at turbine inlet (p1) = %f barrn',p1)
fprintf('Temperature of steam at the turbine inlet (t1) = %f Deg.celsiusrn',t1)
fprintf('Enthalpy of steam at the turbine inlet (h1) = %f KJ/Kgrn',h1)
fprintf('Entropy of steam at the turbine inlet (s1) = %f KJ/Kg Krn',s1)

% State 2
fprintf('At state point 2:rn')
fprintf('Pressure at condenser inlet (p2) = %f barrn',p2)
fprintf('Temperature of steam at condenser inlet (t2) = %f Deg.celsiusrn',t2)
fprintf('Enthalpy of steam at condenser inlet (h2) = %f KJ/Kgrn',h2)
fprintf('Entropy of steam at condenser inlet (s2) = %f KJ/Kg Krn',s2)
fprintf('Dryness fraction of steam (X2) = %frn',X2)

% State 3 
fprintf('At state point 3:rn')
fprintf('Pressure at the feed water pump inlet (p3) = %f barrn',p3)
fprintf('Temperature at the feed water pump inlet (t3) = %f Deg.celsiusrn',t3)
fprintf('Enthalpy of water at feed pump inlet (h3) = %f KJ/Kgrn',h3)
fprintf('Entropy of water at feed pump inlet (s3) = %f KJ/Kg Krn',s3)

% State 4
fprintf('At state point 4:rn')
fprintf('Pressure of water at the boiler inlet (p4) = %f barrn',p4)
fprintf('Temperature of water at the boiler inlet (t4) = %f Deg.celsiusrn',t4)
fprintf('Enthalpy of water at the boiler inlet (h4) = %f KJ/Kgrn',h4)
fprintf('Entropy of water at the boiler inlet (s4) = %f KJ/Kg Krn',s4)

% Outputs
fprintf('Thermal efficiency = %frn',therm_eff)
fprintf('Net work output (W_net) = %f KJ/Kgrn',W_net)
fprintf('Back work ratio = %frn',W_back)
fprintf('Specific steam consumption = %f Kg/KJrn',SSC)

 Results:

A] T-S diagram of the rankine cycle.

 B] H-S diagram of the rankine cycle:

C] Detailed results in the command window:

Performance evaluation:

• The resulted thermal efficiency of rankine cycle is 36.17% .
• The net work output achieved is 1117.851878 KJ/Kg.
• The back work ratio is calculated as 0.2686%.
• Specific steam consumption is 3.2204 Kg/KJ.

Errors encountered:

Mostly I faced errors related to invalid calculation results or mismatched results. However, errors pertaining to MATLAB tools and commands were few.

Still I will share one of them.

1.  Error due to invalid expression:

Actually, while scripting I forgot to call the function and still gave my inputs. 

Command window:

Editor window:

Solution:

I called the "XSteam" function and gave my inputs then the program executed my inputs. After facing this error I carefully called function everytime when I calculated state variables.

Conclusion: In this project, all the thermodynamic properties pertaining to the rankine cycle were successfully calculated and their respective plots were resulted. The approach used in this project is to initially comprehend the basics of rankine cycle, its t-s and h-s diagram and to understand the variables that are used for evaluation. After that, start with MATLAB coding in step by step manner and run the code at specific intervals to validate working of the code. By using the given input conditions, we evaluated the efficiency of the cycle and work obtained. However, we can calculate the actual cycle efficiency and work by modifying code and accounting for loss at various stages.

References:

1. https://en.wikipedia.org/wiki/Rankine_cycle#:~:text=The%20Rankine%20cycle%20is%20a,work%20while%20undergoing%20phase%20change.

2. http://www.thermopedia.com/content/1072/