Project 2 - Rankine cycle Simulator

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

% Step 3
% Taking all inputs arguments
P1 = input('Enter the Pressure at the 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): ');

fprintf('          RESULT\n');

%% Calculation for General properties
% Step 4
% Calculation of Cp value at P1 & P2 for saturated liquid & vapour
CP1_L=XSteam('CpL_p',P1);
CP1_V=XSteam('CpV_p',P1);
CP2_L=XSteam('CpL_p',P2);
CP2_V=XSteam('CpV_p',P2);

% Step 5
% Calculation of Saturation Temperature at P1 & P2
Tsat_P1=XSteam('Tsat_p',P1);
Tsat_P2=XSteam('Tsat_p',P2);

%% Calculation of Enthalpy & Entropy at saturation points
% Step 6
% Calculation of saturation properties at state 1
Hg_1=XSteam('hV_p',P1);
Sg_1=XSteam('sV_p',P1);

% Step 7
% Calculation of saturation properties at state 2
Hg_2=XSteam('hV_p',P2);
Sg_2=XSteam('sV_p',P2);

% Step 8
% Calculation of saturation properties at state 3
Hf_3=XSteam('hL_p',P2);
Sf_3=XSteam('sL_p',P2);

% Step 9
% Calculation of saturation properties at state 4
Hf_4=XSteam('hL_p',P1);
Sf_4=XSteam('sL_p',P1);

%% Calculation of Process properties at state 1
% Step 10
% Calculation of Process properties at state 1
fprintf('At State point 1\n');

fprintf('P1 is :%f Bar\n',P1);

fprintf('T1 is :%f degC\n',T1);

H1=XSteam('h_pT',P1,T1);
fprintf('H1 is :%f kJ/kg\n',H1);

S1=XSteam('s_pT',P1,T1);
fprintf('S1 is :%f kJ/kg\n\n',S1);

%% Calculation of Process properties at state 2
% Step 11
% Calculation of Process properties at state 2

fprintf('At state point 2\n');

fprintf('P2 is : %f Bar\n', P2);

S2 = S1;
fprintf('S2 is : %f kJ/kgK\n', S2);

%Calculating value of X
if(S2<=Sg_2)
    X2 = (S2 - Sf_3) / (Sg_2-Sf_3);   % Dryness Fraction at state point 2
else
    X2=1;    
end
fprintf('X2 is :%f \n', X2);

%calculating value of T2
if(X2<1)
    T2=Tsat_P2;
else 
    %S2=Sg_2+CP2_V*(log(T2/Tsat_P2))
    T2=(Tsat_P2+273.15)*exp((S2-Sg_2)/CP2_V) ;
    T2=T2-273.15;
end    
fprintf('T2 is : %f degC \n', T2)

if(X2<1)
   H2 = Hf_3 + X2 * (Hf_3-Hg_2);
else
   H2=XSteam('h_pT',P2,T2);
end   
fprintf('H2 is : %f kJ/kg\n\n', H2)

%% Calculation of Process properties at state 3
% Step 12
% Calculation of Process properties at state 3
fprintf('At state point 3\n')

P3 = P2;
fprintf('P3 is : %f Bar\n', P3)

V3 = XSteam('vL_p', P3);

T3 = Tsat_P2;    % No subcooling assumed
fprintf('T3 is : %f degC\n', T3)

H3 =Hf_3;
fprintf('H3 is : %f kJ/kg\n', H3)

S3 =Sf_3;
fprintf('S3 is : %f kJ/kgK\n\n', S3)

%% Calculation of Process properties at state 4
% Step 13
% Calculation of Process properties at state 4

fprintf('At state point 4\n')

P4 = P1;
fprintf('P4 is : %f Bar\n', P4)

S4 = S3;
fprintf('S4 is : %f kJ/kgK\n', S4)

Wp = V3 * (P4 - P3) * 100;
H4 = Wp + H3;
fprintf('H4 is : %f kJ/kg\n', H4)

T4= (Tsat_P1+273.15)/(exp((Sf_4-S4)/CP1_L));
T4=T4-273.15;
fprintf('T4 is : %f Deg celsius \n\n', T4)

%% Calculation for plotting saturation Liquid and Vapour lines 
% Step 14
% Calculation & plotting saturation Liquid and Vapour lines

temp = linspace(1, T1, 1000);
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 .
% Step 15
% Plotting T-s Plot

figure(1)
hold on
plot([S1 S2], [T1 T2], 'r', 'linewidth', 2)
text(S1,T1,' 1')
text(S2,T2,' 2')
plot([S2 Sg_2], [T2 Tsat_P2], 'r', 'linewidth', 2)
text(Sg_2,Tsat_P2,' 2s')
plot([Sg_2 S3], [Tsat_P2 T3], 'r', 'linewidth', 2)
text(S3,T3,' 3')
plot([S3 S4], [T3 T4], 'r', 'linewidth', 2)
text(S4,T4,' 4')
plot([S4 Sf_4], [T4 Tsat_P1], 'r', 'linewidth', 2)
text(Sf_4,Tsat_P1,' 4s')
plot([Sf_4 Sg_1], [Tsat_P1 Tsat_P1], 'r', 'linewidth', 2)
text(Sg_1,Tsat_P1,' 1s')
plot([Sg_1 S1], [Tsat_P1 T1], 'r', 'linewidth', 2)
plot(S_L, temp, 'b', 'linestyle','--','linewidth', 1)
plot(S_V, temp, 'b', 'linestyle','--','linewidth', 1)

title('T-S Diagram of Ideal Rankine Cycle')
ylabel('Temprature (degC)')
xlabel('Specific Entropy (kJ/kgK)')

% Step 15
% Plotting H-S Plot
figure(2)
hold on
plot([S1 S2], [H1 H2], 'r', 'linewidth', 2)
text(S1,H1,' 1')
text(S2,H2,' 2')
plot([S2 Sg_2], [H2 Hg_2], 'r', 'linewidth', 2)
text(Sg_2,Hg_2,' 2s')
plot([Sg_2 S3], [Hg_2 H3], 'r', 'linewidth', 2)
text(S3,H3,' 3')
plot([S3 S4], [H3 H4], 'r', 'linewidth', 2)
text(S4,H4,' 4')
plot([S4 Sf_4], [H4 Hf_4], 'r', 'linewidth', 2)
text(Sf_4,Hf_4,' 4s')
plot([Sf_4 Sg_1], [Hf_4 Hg_1], 'r', 'linewidth', 2)
text(Sg_1,Hg_1,' 1s')
plot([Sg_1 S1], [Hg_1 H1], 'r', 'linewidth', 2)
plot(S_L, H_L, 'b','linestyle','--', 'linewidth', 1.5)
plot(S_V, H_V, 'b','linestyle','--', 'linewidth', 1.5)      
title('H-S Diagram of Ideal Rankine Cycle')
ylabel('Specific Enthalpy (kJ/kg)')
xlabel('Specific Entropy (kJ/kgK)')

%% Work done, heat, SSC and efficiency calculations
% Step 16
% Calculation of Work done, heat, SSC and efficiency
Wt = H1 - H2;
fprintf('Work by turbine is : %f kJ/kg\n', Wt)
fprintf('Work by compressor is : %f kJ/kg\n', Wp)
Wnet = Wt - Wp;
fprintf('Net Wok done is : %f kJ/kg\n', Wnet)
Qin = H1 - H4;
eff_thermal = Wnet*100 / Qin ;    % Effieciency
fprintf('Thermal efficiency is : %f', eff_thermal); disp('%')
SSC = 3600 / Wt;    % Specific steam consuption
fprintf('S.S.C is : %f kg/kWh ', SSC)