Project 2 - Rankine cycle Simulator

%Project-2
% Rankine Cycle Simulator
clear all;
close all;
clc;

%Display According the Output Window 
disp('         RANKINE CYCLE SIMULATOR          ');
disp('1-2 is Isentropic Expansion in the Turbine');
disp('2-3 is Constant Pressure Heat Rejection by the Condenser');
disp('3-4 is Isentropic Compression in the Pump');
disp('4-1 is Constant Pressure Heat Addition by the Boiler');

% Initial parameters
t  = linspace(0,373.15,400);
P1 = input('Enter the pressure at the Turbine inlet(in bar)   : ');
T1 = input(['Enter the temperature at the Turbine inlet(in ' char(176) 'C) : ']);
P2 = input('Enter the pressure at the Condenser(in bar)       : ');
disp(' ')

% State 1
h1    = XSteam('h_pT',P1,T1);               
s1    = XSteam('s_pT',P1,T1);               
hl1   = XSteam('hL_p',P1);                 
hg1   = XSteam('hV_p',P1);                 
sl1   = XSteam('sL_p',P1);                  
sg1   = XSteam('sV_p',P1);                  
Tsat1 = XSteam('Tsat_p',P1);             

% State 2
s2    = s1;              
sl2   = XSteam('sL_p',P2);                 
sg2   = XSteam('sV_p',P2);                  
d     = (s2-sl2)/(sg2-sl2);                 
hl2   = XSteam('hL_p',P2);
hg2   = XSteam('hV_p',P2);
h2    = hl2 + d*(hg2-hl2);
T2    = XSteam('T_ph',P2,h2);
Tsat2 = XSteam('Tsat_p',P2); 

% State 3
h3 = hl2;
s3 = sl2;
P3 = P2;
T3 = T2;

% State 4
s4 = s3;
P4 = P1;
h4 = XSteam('h_ps',P4,s4);
T4 = XSteam('T_ps',P4,s4);


% Plot the graph
figure(1), clf
hold on

for i = 1:length(t)
    plot(XSteam('sL_T',t(i)),t(i),'b.')
    plot(XSteam('sV_T',t(i)),t(i),'b.')
end

% Plot Command
plot([s1 s2],[T1 T2],'linew',2)
text(s1,T1,'1','FontSize',15)
if d > 1
    T3 = Tsat2;
    plot([s2 sg2 s3 s4],[T2 Tsat2 T3 T4],'linew',2)
    T2 = Tsat2;
    text(s2,T2,'2','FontSize',15)
    text(s3,T3,'3','FontSize',15)
else
    plot([s2 s3 s4],[T2 T3 T4],'linew',2)
    text(s2,T2,'2','FontSize',15)
    text(s3,T3,'3','FontSize',15)
end

n = linspace(T1,T2,500);
for i = 1:length(n)
    plot(XSteam('s_pT',P1,n(i)),n(i),'.','linew',2)
end

plot([sl1 sg1],[Tsat1 Tsat1],'linew',2)
text(sl1,Tsat1,'4','FontSize',15)
xlabel('Entropy')
ylabel('Temperature')
title('T-S')

% Next Graph
figure(2), clf
hold on

for i = 1:length(t)
    plot(XSteam('sL_T',t(i)),XSteam('hL_T',t(i)),'g.')
    plot(XSteam('sV_T',t(i)),XSteam('hV_T',t(i)),'g.')
end
%Plotting 
plot([s1 s2 s3 s4 s1],[h1 h2 h3 h4 h1],'linew',2)
text(s1,h1,'1','FontSize',15)
text(s2,h2,'2','FontSize',15)
text(s3,h3,'3','FontSize',15)
text(s4,h4+90,'4','FontSize',15)

xlabel('Entropy')
ylabel('Enthalpy')
title('H-S')

%Work Done by Various Sources
% Turbine
Wt = h1 - h2

% Pump
Wp = h4 - h3

% Dome
Wnet = Wt - Wp

% Calculate the Thermal efficiency
Ntherm = (1 - ((h2 - h3)/(h1 - h4)))*100

% Back work ratio
BWR = Wp/Wt

% Specific Steam Conductivity
S.S.C = 3600/Wt

% Displaying commands for all
disp('RESULTS :')
disp('At state point 1: ')                                                   
fprintf('P1 = %.2f bar n',P1)
fprintf(['T1 = %.2f' char(176) 'C n'],T1)
fprintf('h1 = %.2f kJ/kg n',h1)
fprintf('s1 = %.2f kJ/kgK n',s1)
disp(' ')

disp('At state point 2: ')
fprintf('P2 = %.2f bar /n',P2)
fprintf(['T2 = %.2f' char(176) 'C n'],T2)
fprintf('h2 = %.2f kJ/kg n',h2)
fprintf('s2 = %.2f kJ/kgK n',s2)
disp(' ')

disp('At state point 3: ')
fprintf('P3 = %.2f bar n',P3)
fprintf(['T3 = %.2f' char(176) 'C n'],T3)
fprintf('h3 = %.2f kJ/kg n',h3)
fprintf('s3 = %.2f kJ/kgK n',s3)
disp(' ')

disp('At state point 4: ')
fprintf('P4 = %.2f bar n',P4)
fprintf(['T4 = %.2f' char(176) 'C n'],T4)
fprintf('h4 = %.2f kJ/kg n',h4)
fprintf('s4 = %.2f kJ/kgK n',s4)
disp(' ')

fprintf('Work done by turbine        = %.2f kJ/kg n', Wt)
fprintf('Work done by pump           = %.2f kJ/kg n', Wp)
fprintf('Net work done               = %.2f kJ/kg n', Wnet)
fprintf('Thermal efficiency          = %.2f %% n', Ntherm)
fprintf('Back-work ratio             = %f n', BWR)
fprintf('Specific Steam Conductivity = %.2f kg/kWh n', S.S.C)