Project 2 - Rankine cycle Simulator

clear all;
close all;
clc

disp(' RANKINE CYCLE SIMULATOR')
disp('1-2 Isentropic Expansion in the Turbine')
disp('2-3 Constant Pressure Heat Rejection by the Condenser')
disp('3-2 Isentropic Compression in the Pump')
disp('4-1 Constant Pressure Heat Addition by the Boiler')

% Required Inputs
P1 = 30
T1 = 400
P2 = 0.5
disp (' RESULTS')

% getting Cp and Cv values from Xsteam
% Input detail for Vapour
high_CpV= XSteam('CpV_p', P1);
low_CpV=XSteam('CpV_p', P2);
high_CvV=XSteam('CvV_p', P1);
low_CvV=XSteam('CvV_p', P2);

% Input detail for Liquid
high_CpL= XSteam('CpL_p', P1);
low_CpL=XSteam('CpL_p', P2);
high_CvL=XSteam('CvL_p', P1);
low_CvL=XSteam('CvL_p', P2);

% Reading from Xctream file for important parameters
% Input detail for Pressure
T1s = XSteam('Tsat_p', P1);
S1g = XSteam('sV_p', P1);
S1f = XSteam('SL_p', P1);
S1fg = S1g - S1f;
H1f = XSteam('hL_p', P1);
H1g = XSteam('hV_p', P1);
H1fg = H1g - H1f;

% Input detail for Low Pressure
S2g = XSteam('sV_p', P2);
S2f = XSteam('SL_p', P2);
S2fg = S2g - S2f;
T2s = XSteam('Tsat_p', P2);
H2f = XSteam('hL_p', P2);
H2g = XSteam('hV_p', P2);
H2fg = H2g - H2f;

% State point 1
disp ('At state point 1:')
fprintf ('P1 is : %f Barn', P1)
fprintf ('T1 is : %f Deg celsiusn', T1)
H1 = XSteam('h_pt', P1, T1);
fprintf ('H1 is : %f kJ/kgn', H1)
S1 = XSteam('s_pt', P1, T1);
fprintf('S1 is : %f kJ/kgKnn', S1)

% State point 2
disp('At state point 2:')
fprintf('P2 is : %f Barn', P2)
S2 = S1;
fprintf ('S2 is : %f kJ/kgKn', S2)
% Dryness Fraction of point 2 in T-s
if(S2<=S2g)
X2 = (S2 - S2f) / (S2fg);
else
X2=1;
end
fprintf('X2 is :%f n', X2)

% Calculating T2 for superheat, dry, wet
if(X2<1)
T2=T2s;
else
% Now calculating T2 by s2=s2g+low_cpV*(log(T2/T2s))
T2=(T2s+273.15)*exp((S2-S2g)/low_CpV) ;
T2=T2-273.15;
end
fprintf('T2 is : %f Deg celsius n', T2)
if(X2<1)
H2 = H2f + X2 * H2fg;
else
H2=H2g+ low_CpV*(T2-T2s);
end
fprintf('H2 is : %f kJ/kgn', H2)

% State point 3 where assuming no subcooling
disp('At state point 3:')
P3 = P2;
fprintf('P3 is : %f Barn', P3)
V3 = XSteam('vL_p', P3);
T3 = T2s;
fprintf('T3 is : %f Deg celsiusn', T3)
H3 =H2f;
fprintf('H3 is : %f kJ/kgn', H3)
S3 =S2f;
fprintf('S3 is : %f kJ/kgKnn', S3)

% State point 4
disp('At state point 4:')
P4 = P1;
fprintf('P4 is : %f Barn', P4)
S4 = S3;
fprintf('S4 is : %f kJ/kgKn', S4)
Wp = V3 * (P4 - P3) * 100;
H4 = Wp + H3;
fprintf('H4 is : %f kJ/kgn', H4)

% For calculating T4 s4=s3=sf1-high_cpL*ln(t1s/t4)
term=exp((S1f-S4)/high_CpL);
T4= T1s/term;
fprintf('T4 is : %f Deg celsius nn', T4)

% Work done, heat, SSC and efficiency calculations
Wt = H1 - H2;
fprintf('Work by turbine is : %f kJ/kgn', Wt)
fprintf('Work by compressor is : %f kJ/kgn', Wp)
Wnet = Wt - Wp;
fprintf('Net Wok done is : %f kJ/kgn', Wnet)
Qin = H1 - H4;
eff_thermal = Wnet / Qin ; % Effieciency
fprintf('Thermal efficiency is : %f', eff_thermal); disp('%')
SSC = 3600 / Wt;
fprintf('S.S.C is : %f kg/kWh ', SSC)
% Calculation for plotting Liquid and Vapour lines
temp = linspace(1, 373, 500);
for i = 1 : length(temp)
H_L(i) = XSteam('hL_T', temp(i));
H_V(i) = XSteam('hV_T', temp(i));
S_L(i) = XSteam('sL_T', temp(i));
S_V(i) = XSteam('sV_T', temp(i));
end
% Ploting T - S and H - S Diagrams
% T-S diagram of Ideal Rankine Cycle
figure(1)
title('T-S Diagram of Ideal Rankine Cylce')
ylabel('Temprature (C)')
xlabel('Specific Entropy (kJ/kgK)')
grid on
hold on
plot([S1 S2], [T1 T2], 'r', 'linewidth', 2)
text(S1,T1,' 1')
text(S2,T2,' 2')
plot([S2 S2g], [T2 T2s], 'g', 'linewidth', 2)
text(S2g,T2s,' 2s')
plot([S2g S3], [T2s T3], 'g', 'linewidth', 2)
text(S3,T3,' 3')
plot([S3 S4], [T3 T4], 'c', 'linewidth', 2)
text(S4,T4,' 4')
plot([S4 S1f], [T4 T1s], 'y', 'linewidth', 2)
text(S1f,T1s,' 4s')
plot([S1f S1g], [T1s T1s], 'k', 'linewidth', 2)
text(S1g,T1s,' 1s')
plot([S1g S1], [T1s T1], 'm', 'linewidth', 2)
plot(S_L, temp, 'b', 'linestyle','--','linewidth', 2)
plot(S_V, temp, 'b', 'linestyle','--','linewidth', 2)
xlim([0.00 10.00])
ylim([0 400])
% H-S Diagram
figure(2)
title('H-S Diagram', 'FontSize', 15)
ylabel('Specific Enthalpy (kJ/kg)')
xlabel('Specific Entropy (kJ/kgK)')
grid on
hold on
plot([S1 S2], [H1 H2], 'r', 'linewidth', 2);
text(S1,H1,' 1')
plot([S2 S3], [H2 H3], 'g', 'linewidth', 2);
text(S2,H2,' 2')
plot([S3 S4], [H3 H4], 'm', 'linewidth', 2);
text(S3,H3,' 3')
plot([S4 S1], [H4 H1], 'k', 'linewidth', 2);
text(S4,H4,' 4')
plot(S_L, H_L, 'b','linestyle','--', 'linewidth', 2.5)
plot(S_V, H_V, 'b','linestyle','--', 'linewidth', 2.5)
xlim([0.00 10.00])
ylim([0 3500])