Project 2 - Rankine cycle Simulator

Code for Rankine cycle

close all
clear all
clc

% Rankine Cycle Simulation
% Given Input
% Turbine Inlet temperature = 400° C
% Turbine Inlet Pressure = 30 bars
% Turbine Outlet Pressure = 0.05 bars = Condenser Inlet Pressure

% QH = heat transferred to the working fluid (kJ/kg)
% QL = heat rejected from the working fluid (kJ/kg)
% WT = work transferred from the working fluid (kJ/kg)
% WP = work transferred into the working fluid (kJ/kg) 
% x  = Dryness fraction
% Process inside the rankine cycle 
% Process 1-2: Constant pressure heat transfer process in the boiler
% Process 2-3: Reversible adiabatic expansion process in the steam turbine
% Process 3-4: Constant pressure heat transfer process in the condenser and
% Process 4-1: Reversible adiabatic compression process in the pump. 

% For calculating the W_net(Net work) value of H should be calculated

% Process 1-2: Constant pressure heat transfer process in the boiler
% State point 1 
P1=30; % Bar inlet pressure in Turbine
T1=400; % Degree Celcius inlet Temperature in Turbine
H1=XSteam('h_pt',P1,T1); % Enthalpy at point 1
S1=XSteam('s_pt',P1,T1); % Entropy at point 1

% Process 2-3: Reversible adiabatic expansion process in the steam turbine
% State point 2 
P2=0.05; % Bar Turbine Outlet Pressure
T2=XSteam('Tsat_p',P2);% Degree Celcius Turbine outlet Temperature 
S2=S1; % Entropy at point 2 (Reversible adiabatic process)                           
sf2= XSteam('sL_p',P2); % Saturated liquid entropy                       
sg2= XSteam('sV_p',P2); % Saturated vapour entropy
sfg2= sg2-sf2;  % sfg = sg − sf
x= (S2-sf2)/sfg2; % Dryness fraction
hf2=XSteam('hL_p',P2); % Saturated liquid enthalpy            
hg2=XSteam('hV_p',P2); % Saturated vapour enthalpy
hfg2 = hg2-hf2; % hfg = hg − hf
H2=hf2+x.*(hfg2);% Enthalpy at point 2

% Process 3-4: Constant pressure heat transfer process in the condenser
% State point 3
P3=P2; % Bar Condenser Inlet Pressure (Reversible Isobaric heat rejection process)
T3=T2; % Degree Celcius Condenser Inlet Temperature
H3=XSteam('hL_p',P3);% Enthalpy at point 3
S3=XSteam('sL_p',P3);% Entropy at point 3 
V3=XSteam('vL_p',P3);% Saturated liquid volume at point 3

% Process 4-1: Reversible adiabatic compression process in the pump.
% State point 4
P4=P1; % Bar Pump outlet Pressure(It will be equal to inlet pressure in Turbine)
S4=S3;% Entropy at point 4
Cp=XSteam('Cp_ps',P4,S4); % Specific heat capacity at constant pressure
Cv=XSteam('Cv_ps',P4,S4); % Specific heat capacity at constant volume
k= Cp/Cv; % heat capacity ratio
o=(k-1)/k; 

T4= T3*((P4/P3)^o); % Degree Celcius Boiler Inlet Temperature (The relation came from ideal gas law for isentropic process)

% for calculating the enthalpy at point 4 equate pump work=change in enthalpy
H4=H3+(100*V3*(P4-P3)); % Enthalpy at point 4 (note 100 is multiplied because to balance the unit on both side in KJ/Kg pressure is change from bar to pascal and the to (10^5 * N/m^2)

% State point 5 and state point 6
P5=P1;
P6=P5;
S5=XSteam('sL_p',P1);% Entropy at point 5
S6=XSteam('sV_p',P1);% Entropy at point 6
H5=XSteam('hL_p',P1);% Enthalpy at point 5
H6=XSteam('hV_p',P1);% Enthalpy at point 6
T6=XSteam('Tsat_p',P1);% Degree Celcius
T5=T6;

% Calculation of Net work done and Back work ratio
% Net work done Wnet=Wt-Wp (Wt is work done by Turbine and Wp is work done by pump)

% Work done by Turbine
Wt= H1-H2;

% Work done by pump

Wp=H4-H3;

% Net work done 

Wnet=Wt-Wp;

% Heat Supplied in turbine

Qadd= H1-H4;

% Heat rejected in condenser

Qrej= H2-H3;

% Thermal efficiency of Rankine cycle

Ntherm= (Wnet/Qadd)*100;

% Back work ratio (Pump work input to work developed by turbine)

W_back=Wp/Wt;

% Plotting the saturation line 

T=linspace(1,800,1000);

for i=1:length(T)
sf(i)=XSteam('sL_T',T(i)); % Saturated liquid entropy
sg(i)=XSteam('sV_T',T(i)); % Saturated vapour entropy
hf(i)=XSteam('hL_T',T(i)); % Saturated liquid enthalpy
hg(i)=XSteam('hV_T',T(i)); % Saturated vapour enthalpy
end

% Plot the T-S curve for Rankine cycle
figure(1)
hold on

plot(sf,T,'color','r','LineStyle','--','LineWidth',2);
plot(sg,T,'color','r','LineStyle','--','LineWidth',2);
plot([S1 S2],[T1 T2],'*','color','b','LineStyle','-','LineWidth',1);
plot([S2 S3],[T2 T3],'*','color','g','LineStyle','-','LineWidth',1);
plot([S3 S4],[T3 T4],'*','color','r','LineStyle','-','LineWidth',1);
plot([S4 S5],[T4 T5],'*','color','k','LineStyle','-','LineWidth',1);
plot([S5 S6],[T5 T6],'*','color','y','LineStyle','-','LineWidth',1);
plot([S6 S1],[T6 T1],'*','color','c','LineStyle','-','LineWidth',1);
xlabel('Entropy (J⋅K−1⋅kg−1)');
ylabel('Temperature(° C)');
title('Temperature vs Entropy (T-S) curve');

figure(2)
hold on
plot(sf,hf,'color','r','LineStyle','--','LineWidth',2);
plot(sg,hg,'color','r','LineStyle','--','LineWidth',2);
plot([S1 S2],[H1 H2],'*','color','b','LineStyle','-','LineWidth',1);
plot([S2 S3],[H2 H3],'*','color','g','LineStyle','-','LineWidth',1);
plot([S3 S4],[H3 H4],'*','color','r','LineStyle','-','LineWidth',1);
plot([S4 S5],[H4 H5],'*','color','k','LineStyle','-','LineWidth',1);
plot([S5 S6],[H5 H6],'*','color','y','LineStyle','-','LineWidth',1);
plot([S6 S1],[H6 H1],'*','color','c','LineStyle','-','LineWidth',1);
xlabel('Entropy (J⋅K−1⋅kg−1)');
ylabel('Enthalpy(KJ/Kg)');
title('Enthalpy vs Entropy (H-S) curve');

% Printing all the Outputs
Pressure=[P1,P2,P3,P4,P5,P6];
Temperature=[T1,T2,T3,T4,T5,T6];
Enthalpy=[H1,H2,H3,H4,H5,H6];
Entropy=[S1,S2,S3,S4,S5,S6];

for i=1:6
    fprintf('State Point %1.0f \n',i);
    fprintf('P%1.0f: %0.3f (Bar)  \n',i,Pressure(i));
    fprintf('T%1.0f: %0.3f (° C) \n',i,Temperature(i));
    fprintf('H%1.0f: %0.3f (KJ/Kg) \n',i,Enthalpy(i));
    fprintf('S%1.0f: %0.3f (J⋅K−1⋅kg−1) \n',i,Entropy(i));
end

 fprintf('Work done by turbine: %0.3f (KJ/Kg)',Wt);
 fprintf('Work done by Pump: %0.3f (KJ/Kg) \n',Wp);
 fprintf('Heat Supplied in turbine: %0.3f (KJ) \n',Qadd);
 fprintf('Heat rejected in condenser: %0.3f (KJ) \n',Qrej);
 fprintf('Net Work output: %0.3f (KJ/Kg) \n',Wnet);
 fprintf('Thermal Efficiency: %0.2f (KJ/Kg) \n',Ntherm);
 fprintf('Back Work Ratio: %0.3f (KJ/Kg) \n',W_back)