Project 2 - Rankine cycle Simulator

%Program to simulate the Rankine cycle
%Plot t-s and h-s diagram
clear all
close all
clc

%% Display the title and the process of Rankine cycle
disp('Rankine Cycle')
disp('1-2 The Isentropic Expansion in the Turbine')
disp('2-3 Constant Pressure Heat Rejection by the Condenser')
disp('3-4 Isentropic Compression in the Pump')
disp('4-1 Constant Pressure Heat Addition by the Boiler')

%% Input the initial value for Rankine cycle
p1 = input('Enter the pressure at turbine inlet, p1(in bar) = ');
T1 = input('Enter the temperature at turbine inlet, T1(in Degree Celcius) = ');
p2 = input('Enter the pressure at the condenser, p2(in bar) = ');

%% Extrating data and definind the variable
%State 1
h1 = XSteam('h_pT', p1, T1); % call XSteam function and find h1 with the help of input p1 and t1
s1 = XSteam('s_pT', p1, T1);
%State 2
s2 = s1;
sf2 = XSteam('sL_p', p2);
sg2 = XSteam('sV_p', p2);
x2 = (s2-sf2)/(sg2-sf2); %Dry ness fraction = mass of dry steam/Total mass of wet steam
hf2 = XSteam('hL_p', p2);
hg2 = XSteam('hV_p', p2);
h2 = hf2 + (x2 * (hg2-hf2));
%State3
p3 = p2; %Constant pressure heat rejection
s3 = XSteam('sL_p', p3);
h3 = XSteam('hL_p', p3);
%State4
s4 = s3;
p4 = p1;
v3 = XSteam('vL_p', p3);
wp = v3*(p4-p3)*100; % Pump work
h4 = h3+wp;

%% Work done by turbine
wt = h1-h2;
%% Net Work
wnet = wt-wp;
%% Heat into the system
Qin = h1-h4;
%% Heat out of the system
Qout = h2-h3;
%% Thermal efficiency
eff = (wnet/Qin)*100;
%% Specific heat consumption
SSC = (3600/wnet);

%% Temperature calculating at the state point
T3 = XSteam('Tsat_p', p3);
T2 = T3;
cp4 = XSteam('cp_ps', p4,s4);
cv4 = XSteam('cv_ps', p4,s4);
k = cp4/cv4;
T4 = T3/((p3/p4)^((k-1)/k));
T6 = XSteam('Tsat_p',p1);
T5 = T6;
s5 = XSteam('sL_p', p1);
s6 = XSteam('sV_p', p1);
h5 = XSteam('hL_p', p1);
h6 = XSteam('hV_p', p1);
%% Creating dome
T = linspace(1,375,500);
for i = 1:length(T)
    sf(i) = XSteam('sL_T', T(i));
    sg(i) = XSteam('sV_T', T(i));
    hf(i) = XSteam('hL_T', T(i));
    hg(i) = XSteam('hV_T', T(i));
    
end

%% Plot T-s diagram
figure(1)
hold on
grid on
plot(sf,T, 'b','linestyle','--','linewidth',3) %plotting a dome
plot(sg,T, 'b','linestyle','--','linewidth',3) %plotting a dome
plot([s1 s2], [T1 T2],'r','linewidth',2)
plot([s2 s3], [T2 T3],'r','linewidth',2)
plot([s3 s4], [T3 T4],'r','linewidth',2)
plot([s4 s5], [T4 T5],'r','linewidth',2)
plot([s5 s6], [T5 T6],'r','linewidth',2)
plot([s6 s1], [T6 T1],'r','linewidth',2)
text(s1,T1,'1')
text(s2,T2,'2')
text(s3,T3,'3','HorizontalAlignment','right','VerticalAlignment','top')
text(s4,T4,'4','HorizontalAlignment','right','VerticalAlignment','bottom')
text(s5,T5,'5')
text(s6,T6,'6')
xlabel('Entropy (KJ/KgK)')
ylabel('Temperature (C)')
title('Temperature VS Entropy')

%% Plot h-s diagram
figure(2)
hold on
grid on
plot(sf,hf, 'b','linestyle','--','linewidth',3) %plotting a dome
plot(sg,hg, 'b','linestyle','--','linewidth',3) %plotting a dome
plot([s1 s2], [h1 h2],'g','linewidth',3)
plot([s2 s3], [h2 h3],'g','linewidth',3)
plot([s3 s4], [h3 h4],'g','linewidth',3)
plot([s4 s5], [h4 h5],'g','linewidth',3)
plot([s5 s6], [h5 h6],'g','linewidth',3)
plot([s6 s1], [h6 h1],'g','linewidth',3)
text(s1,h1,'1')
text(s2,h2,'2')
text(s3,h3,'3','HorizontalAlignment','right','VerticalAlignment','top')
text(s4,h4,'4','HorizontalAlignment','right','VerticalAlignment','bottom')
text(s5,h5,'5')
text(s6,h6,'6')
xlabel('Entropy (KJ/KgK)')
ylabel('Enthalpy (KJ/Kg)')
title('Enthalpy VS Entropy')

%% Display result
disp('----------RESULT----------');
fprintf('\n At state point 1:\n')
fprintf('P1 is: %f Bar \n',p1)
fprintf('T1 is: %f Celcius \n',T1)
fprintf('H1 is: %f KJ/Kg \n',h1)
fprintf('S1 is: %f KJ/KgK \n',s1)

fprintf('\n At state point 2:\n')
fprintf('P2 is: %f Bar \n',p2)
fprintf('T2 is: %f Celcius \n',T2)
fprintf('H2 is: %f KJ/Kg \n',h2)
fprintf('S2 is: %f KJ/KgK \n',s2)

fprintf('\n At state point 3:\n')
fprintf('P3 is: %f Bar \n',p3)
fprintf('T3 is: %f Celcius \n',T3)
fprintf('H3 is: %f KJ/Kg \n',h3)
fprintf('S3 is: %f KJ/KgK \n',s3)

fprintf('\n At state point 4:\n')
fprintf('P4 is: %f Bar \n',p4)
fprintf('T4 is: %f Celcius \n',T4)
fprintf('H4 is: %f KJ/Kg \n',h4)
fprintf('S4 is: %f KJ/KgK \n',s4)

fprintf('\n Wt is: %f KJ/Kg \n',wt)
fprintf('Wp is: %f KJ/Kg \n',wp)
fprintf('Wnet is: %f KJ/Kg \n',wnet)
fprintf('S.S.C is: %f Kg/KWh \n',SSC)
fprintf('Thermal efficiency is: %f % \n',eff)