Project 2 - Rankine cycle Simulator
% A program to Simulate the Rankine Cycle
close all;
clear;
clc
%{
The typical Rankine Cycle processes are as follows
Pump: Compression of the fluid to high pressure using a pump Steps 1 - 2
Boiler: The compressed fluid is heated to the final temperature Steps 2 - 3
Turbine: Expansion of the vapor in the turbine. Steps 3 - 4)
Condenser: Condensation of the vapor in the condenser Steps 4 to 1
%}
%{
Isobaric Heat transfer or constant pressure heat addition in the boiler
qIn = h3 - h2
Isentropic expansion in the turbine
Wtout = h3 - h4
Isobaric Heat rejection in the condenser
qout = h4 - h1
Isentropic compression in the pump
Workdone on pump
Wpump = h2 − h1
h1 = hf at P1 and v1 = vf at P1
Thermal Efficiency of the Rankine Cycle
η = (qin - qout)/qin = ((h3 - h2) - (h4 - h1))/(h3 - h2)
η = Wnet/qin = (wtout - wpump)/qin = (W34 - W12)/q23
%}
%{
However, for this project, the turbine inlet shall be taken as state point
1 and outlet shall be state point 2 since the conditions such as inlet temperature (400 degree celsius),
inlet pressure (30 bars) and outlet pressure (0.05 bars) were given.
%}
% Conditions at state 1 are as follows:
p1 = 30; % Pressure at Turbine inlet
T1 = 400; % Temperature at Turbine inlet
p2 = 0.05; % Pressure at state 2; Condenser inlet
% Obtaining relevant data for state 1
disp('Isentropic Expansion - State 1')
fprintf('p1 is: %f bar \n',p1)
fprintf('t1 is: %f deg celsius',T1)
h1 = XSteam('h_pt',p1,T1); % Enthalpy at state 1 for given pressure and temperature
fprintf('h1 is: %f kJ/kg \n',h1)
s1 = XSteam('s_pt',p1,T1); % Entropy at state 1 for given pressure and temperature
fprintf('s1 is: %f kJ/kg \n',s1)
tsat_1 = XSteam('tsat_p',p1); % Saturation temperature for corresponding pressure at state 1
sg1 = XSteam('sV_p',p1); % Saturated vapour entropy for given pressure
sf1 = XSteam('sL_p',p1); % Saturated liquid entropy for given pressure
sfg1 = sg1-sf1;
hf1 = XSteam('hL_p',p1); % Saturated liquid enthalpy for given pressure
hg1 = XSteam('hV_p',p1); % Saturated vapour enthalpy for given pressure
hfg1 = hg1 - hf1;
% State 2 Data
s2 = s1; % Constant Entropy between states 1 & 2
disp('Isobaric Heat Rejection - state 2')
fprintf('p2 is: %f bar \n',p2)
fprintf('s2 is: %f kJ/kgk \n',s2)
% Evaluation parameters for state 2
sg2 = XSteam('sV_p',p2); % Saturated vapour entropy for given pressure
sf2 = XSteam('sL_p',p2); % Saturated liquid entropy for given pressure
hf2 = XSteam('hL_p',p2); % Saturated liquid enthalpy for given pressure
hg2 = XSteam('hV_p',p2); % Saturated vapour enthalpy for given pressure
x2 = (s2 - sf2)/(sg2-sf2); % Vapour fraction at state 2
if x2 <= 1 % Liquid - Vapour mixture
h2 = hf2 + x2*(hg2-hf2); % Enthalpy at state 2
T2 = XSteam('T_ph', p2, h2);
else % Working fluid is superheated
Tsat_2 = XSteam('Tsat_p',p2); % Saturated temperature for given pressure
T2 = XSteam('T_ps', p2, s2);
cpV_2 = XSteam('cpV_p', p2);
cvV_2 = XSteam('cvV_p', p2);
h2 = hg2 + cpV*(T2 - Tsat_2);
end
% State 3
p3 = p2; % Constant Pressure between states 2 & 3 (Isobaric Heat Rejection)
disp('Isentropic Compression - state 3')
fprintf('p3 is: %f bar \n',p3)
% Evaluation parameters for state 3
s3 = XSteam('sL_p', p3); % Saturated liquid entropy at state 3 for given pressure
h3 = XSteam('hL_p', p3); % Saturated liquid enthalpy at state 3 for given pressure
v3 = XSteam('vL_p', p3); % Specific volume at state 3
T3 = XSteam('Tsat_p',p3); % Saturation temperature for corresponding pressure at state 3
%T3 = Tsat_3;
% State 4
s4 = s3; % Constant Entropy between states 3 & 4 (Isentropic Compression)
p4 = p1; %
disp('Isobaric Heat Addition - state 4')
fprintf('s4 is: %f kJ/kgk \n',s4)
h4 = XSteam('h_ps', p4, s4); % Enthalpy at state 4
T4 = XSteam('T_ph', p4, h4);
% Relevant parameters of Boiler
% Saturation temperature of Boiler
Tsat_4 = XSteam('Tsat_p',p4);
% Relevant parameters of Liquid
h_sL4 = XSteam('hL_p',p4); % Saturation Liquid Enthalpy
s_sL4 = XSteam('sL_p',p4); % Saturation Liquid Entropy
% Relevant parameters of vapour
h_sV4 = XSteam('hV_p',p4); % Saturation Vapour Enthalpy
s_sV4 = XSteam('sV_p',p4); % Saturation Vapour Entropy
Wpump = h4 - h3; % Work done by pump
wtout = h1 - h2; % Work output from Turbine
% Thermal Efficiency (n) of the Rankine Cycle
qin = h1 - h4;
qout = h3 - h2;
%n = (qin - qout)/qin;
Wnet = (wtout - Wpump);
n = (Wnet/qin)*100;
fprintf('Net Work Wnet is: %.2f kJ/kg \n',Wnet);
fprintf('The Thermal Efficiency n of the Rankine Cycle is: %.2f\n',n);
%{
On a T - V or P - V Diagrams, region to the left of the saturated liquid
line is the compressed liquid region and the region to the right of the
saturated vapour line is the superheated region
%}
% Evaluating Critical parameters at Critical Pressure and Temperature of Water
% Critical Pressure of Water is 220.6 bar
% Critical Temperature of Water is 373.95°C
p_Critical = 220.6;
T_Critical = 373.95;
v_Critical = XSteam('v_pT',p_Critical, T_Critical);
s_Critical = XSteam('s_pT',p_Critical, T_Critical);
%{
s1 > s_Critical, hence fluid is superheated at state 1
For the given Turbine inlet conditions, s1 is evaluated as 6.9233, s4 is evaluated as 0.4763
%}
% Assumptions
% Defining Entropy range from evaluated s1 (vapour region) and s4 (liquid region)
s_sub = 0; % Corresponding to compressed liquid or sub-cooled region
s_sup = 8.2; % Corresponding to super heated region
% Evaluating enthalpies and pressures along Saturation liquid curve
s_sL = linspace(s_sub,s_Critical,500);
for i=1:length(s_sL)
p_sL(i) = XSteam('Psat_s',s_sL(i));
h_sL(i) = XSteam('hL_p',p_sL(i));
end
% Evaluating compressed liquid properties
s_comp = linspace(s4,s_sL4,500);
for ii=1:length(s_comp)
T_comp(ii) = XSteam('T_ps',p4,s_comp(ii));
h_comp(ii) = XSteam('h_ps',p4,s_comp(ii));
end
% Evaluating enthalpies and pressures along Saturation vapour curve
s_sV = linspace(s_Critical,s_sup,500);
for i=1:length(s_sL)
p_sV(i) = XSteam('Psat_s',s_sV(i));
h_sV(i) = XSteam('hL_p',p_sV(i));
end
% Evaluating super heated vapour properties
s_sheat = linspace(s_sV4,s1,500);
for ii=1:length(s_sheat)
T_sheat(ii) = XSteam('T_ps',p1,s_sheat(ii));
h_sheat(ii) = XSteam('h_ps',p1,s_sheat(ii));
end
% Evaluating Saturation Temperature for the assumed entropy range 0 - 8.2 (lines 183 & 184)
s = linspace(s_sub,s_sup,500);
for i=1:length(s)
T_s(i) = XSteam('Tsat_s',s(i));
end
% Graphics & Visualization
% PLOTTING THE T -S Diagram
h = figure;
plot(s,T_s,'-','linewidth',2,'color',[0 0 0])
xlabel('$Entropy(s)$','interpreter','LaTex','FontSize',22)
ylabel('$Temperature(T)$','interpreter','LaTex','FontSize',22,'Rotation',0,'HorizontalAlignment','Right')
title('T-S Diagram','interpreter','LaTex','FontSize',22)
hold on
% Turbine Isentropic Expansion Process 1 - 2
plot([s1 s2],[T1 T2],'linewidth',2,'color',[0 1 0])
% Condenser Heat Rejection Process 2-3
plot([s2 s3],[T2 T3],'linewidth',2,'color',[1, 0, 0])
% Pump Isentropic Compression Process 3 - 4
plot([s3 s4],[T3 T4],'linewidth',2,'color',[0 1 1])
% Plot in compressed liquid region
plot(s_comp,T_comp,'linewidth',2,'color',[0, 1, 1])
% Plot in saturated liquid - vapour region
plot([s_sL4 s_sV4],[Tsat_4 Tsat_4],'linewidth',2,'color',[0.75 0.75 1])
%% Plot in super heated vapour region
plot(s_sheat,T_sheat,'linewidth',2,'color',[0, 0.5, 0.5])
savefig(h,'TS.fig')
% References
% Online Mentoring Session with Arhan Basha
% Thermal Fluid Sciences 5ed by YUNUS A. ÇENGEL, JOHN M. CIMBALA and ROBERT H. TURNER
% https://www.engineeringtoolbox.com/water-properties-d_1508.html