Project 2 - Rankine cycle Simulator

close all
clear all
clc
%-----------------------------------user inputs---------------------%
Temp1 = "What is Temperature 1 in celsius? ";
T1 = input(Temp1); %turbine inlet temperature
Pressure1 = "What is pressure 1 in Bar? ";
P1 = input(Pressure1); %turbine inlet temperature
Pressure2 = "What is pressure 2 in Bar? ";
P2 = input(Pressure2); %turbine outlet temperature
%----------------------------------constant - Values-----------------%
P4 = P1;                      
P3 = P2;    
%----------------------------Steam-Table Data------------------------%
Tsat = XSteam('tsat_p',P1);
hf1 = XSteam('hl_T',Tsat);      % Specific enthalpy - liquid phase
sg1 = XSteam('sv_T',Tsat);      % Specific entropy  - Gas phase
hf2 = XSteam('hl_p',P2);        % Specific enthalpy - liquid phase
hg2 = XSteam('hv_p',P2);        % Specific enthalpy - Gas phase
sf2 = XSteam('sl_p',P2);        % Specific entropy - liquid phase
sg2 = XSteam('sv_p',P2);        % Specific entropy - Gas phase
hf4 = XSteam('hl_p',P4);        % Specific enthalpy - liquid phase
sf4 = XSteam('sl_p',P4);        % Specific entropy - liquid phase
hg4 = XSteam('hv_p',P4);        % Specific enthalpy - Gas phase
sg4 = XSteam('sv_p',P4);        % Specific entropy - Gas phase
%----------------------------Entropy values s1 t0 s4-----------------%
s1 = XSteam('s_pT',P1,T1);     
s2 = s1;                        
s3 = XSteam('sl_p',P3);         
s4 = s3;     
%-----------------------------Temp values T1 to T4--------------------%
T2 = XSteam('tsat_p',P2);       % Saturation temperature a
T3 = T2;                        %same temperature in both state 
% calculating Temp4 by using following equation : T4 = T3*(P4/P3)^((gamma - 1)/gamma)
cp4 = XSteam('cp_ps',P4,s4);    % Specific heat for p4 and s4
cv4 = XSteam('cv_ps',P4,s4);    % Specific heat for p4 and s4
gamma = cp4/cv4;                %ratio of Cp and Cv 
T4 = T3*(P4/P3)^((gamma - 1)/gamma);
%--------------------------- Enthalphy (h) calculation----------------%
h1 = XSteam('h_pT',P1,T1);      % Enthalpy from Steam table 
x2 = (s2 - sf2)/(sg2 - sf2);    % from s1 = sf2 + x2.sfg2  equation
h2 = hf2 + x2*(hg2 - hf2);      % from h2 = hf2 + x2.hfg2  equation
h3 = XSteam('hl_p',P3);         % Enthalpy from Steam table 
vf3 = XSteam('vl_p',P3);        % Specific volume of steam
h4 = h3 + 100*vf3*(P1 - P2);    % fidning h4 value
%------------------------- calculating Qin and out--------------------%
Qin = h1 - h4;                  % Heat to turbine
Qout = h2 - h3;                 % Heat from compressor
%---------------------work done calculation --------------------------%
Wt = h1 - h2;                   % Work_turbine
Wp = h4 - h3;                   % Work_pump
Wnet = Wt - Wp;                 % Net_workDone
%------------------------ cycle characterstic-------------------------%
n_eff = (Wnet/Qin)*100;
S_S_C = 3600/Wnet;
%------------------------- Display output ----------------------------%
fprintf("\n Stage 1-2: Isoentropic Expansion:- \n")
fprintf("\n     - P1 = %0.3f bar",P1)
fprintf("\n     - T1 = %0.2f degree °C",T1)
fprintf("\n     - h1 = %0.2f KJ/kg",h1)
fprintf("\n     - s1 = %0.2f KJ/Kg-K",s1)
fprintf("\n     - Tsat = %0.2f degree °C \n",Tsat)
%%
fprintf("\n Stage 2-3: Constant Pressure Heat:- \n")
fprintf("\n     - P2 = %0.3f bar",P2)
fprintf("\n     - T2 = %0.2f degree °C",T2)
fprintf("\n     - h1 = %0.2f KJ/kg",h2)
fprintf("\n     - s1 = %0.2f KJ/Kg-K",s2)
%%
fprintf("\n Stage 3-4: Isoentropic Compression:- \n")
fprintf("\n     - P3 = %0.3f bar",P3)
fprintf("\n     - T3 = %0.3f degree celsius",T3)
fprintf("\n     - h3 = %0.3f KJ/kg",h3)
fprintf("\n     - s3 = %0.3f KJ/Kg-K \n",s3)
%%
fprintf("\n Stage 4-1: Constant Pressure Heat Addition:- \n")
fprintf("\n     - P4 = %0.3f bar",P4)
fprintf("\n     - T4 = %0.3f degree celsius",T4)
fprintf("\n     - h4 = %0.3f KJ/kg",h4)
fprintf("\n     - s4 = %0.3f KJ/Kg-K \n",s4)
%%
fprintf("\n   - work_done by turbine, Wt = %0.3f KJ/Kg",Wt)
fprintf("\n   - work_done by pump, Wp = %0.3f KJ/Kg",Wp)
fprintf("\n   - Net_work_done, Wnet = %0.3f KJ/Kg",Wnet)
fprintf("\n   - Thermal efficiency (η) = %0.3f",n_eff)
fprintf("\n   - specific_steam_consumption, SSC = %0.3f Kg/KJ\n",S_S_C)
%----------------------------Plotting T-S ------------------------------%
figure(1)                     
hold on  
title(" Temperature T v/s Entropy s ")
xlabel('Entropy S , KJ/Kg-K')
ylabel('Temperature , °C')
plot([s1 s2],[T1 T2],'linewidth',3,'color','b')
plot([s2 s3],[T2 T3],'linewidth',3,'color','m')
plot([s4 sf4],[T4 Tsat],'linewidth',3,'color','r')
plot([sf4 sg1],[Tsat Tsat],'linewidth',3,'color','g')
plot([sg1 s1],[Tsat T1],'linewidth',3,'color','y')
axis([0 8 0 400])
legend('To Condenser', 'To pump','Pump to Boiler','To turbine', 'To turbine');
%----------------------------Plotting h-S ------------------------------%
figure(2)                     
hold on  
title(" Enthalpy-h v/s Entropy-s ")
xlabel('Entropy : S (KJ/Kg-K)')
ylabel('Enthalpy : H (KJ/Kg)')
plot([s1 s2],[h1 h2],'linewidth',3,'color','r')
plot([s2 s3],[h2 h3],'linewidth',3,'color','m')
plot([s3 s4],[h3 h4],'linewidth',3,'color','b')
plot([s4 sf4],[h4 hf4],'linewidth',3,'color','r')
plot([sf4 sg1],[hf4 hg4],'linewidth',3,'color','g')
plot([sg1 s1],[hg4 h1],'linewidth',3,'color','y')
axis([0 10 0 4000])
legend('To Condenser', 'To pump','Pump to Boiler','To turbine','To turbine', 'Expansion in turbine');
%------------------------------------END--------------------------------------------------------------%

Output:

Output :