Project 2 - Rankine cycle Simulator

Program:Rankine Cycle Simulator using Matlab

1.Aim: To simulate Rankine cycle using Matlab and calculate Net Work Output and Back Work Ratio by making use of the standard Xsteam Library.

2.Objectives:

1. To create a Rankine Cycle Simulator using Matlab.
2. To calculate the state points of the Rankine Cycle based on user inputs and display them in the Command Window.
3. To calculate Net Work Output and Back Work Ratio.
4. To plot T-s and h-s diagrams for the given set of inputs.

3.Theory:

A Rankine cycle power system is a heat engine converting thermal energy into work. Similar to the vapor compression heat pump, it comprises four main components: a boiler (sometimes called evaporator), a turbine, a condenser and a pump.

The working fluid in a low-pressure slightly subcooled liquid state is brought to high pressure by the pump. The pump consumes power. The high-pressure liquid enters the boiler, where it is heated-up, vaporized and potentially heated again to leave the boiler in superheated vapor state. Thermal power must be supplied to the boiler by an external heat source. At the outlet of the boiler, the high-pressure vapor enters the turbine where it is expanded down to the low pressure. The turbine delivers mechanical power.

As a function of the nature of the working fluid and turbine inlet state, the fluid at the end of the expansion can be in superheated vapour state or in two-phase state. At the outlet of the turbine, the fluid enters the condenser where it is condensed to liquid state and, then, slightly subcooled. Thermal power rejected at the condenser is dumped into a heat sink. The low-pressure subcooled liquid enters the pump again.
The net power delivered by the Rankine cycle system is the power delivered by the turbine minus that consumed by the pump.

• Process 1-2 Isentropic expansion: The vapour is expanded in the turbine, thus producing work which may be converted to electricity. In practice, the expansion is limited by the temperature of the cooling medium and by the erosion of the turbine blades by liquid entrainment in the vapor stream as the process moves further into the two-phase region. Exit vapor qualities should be greater than 90%.
• Process 2-3 Isobaric Heat Rejection: The vapour-liquid mixture leaving the turbine (2) is condensed at low pressure, usually in a surface condenser using cooling water. In well designed and maintained condensers, the pressure of the vapor is well below atmospheric pressure, approaching the saturation pressure of the operating fluid at the cooling water temperature.
• Process 3-4 Isentropic Compression: The pressure of the condensate is raised in the feed pump. Because of the low specific volume of liquids, the pump work is relatively small and often neglected in thermodynamic calculations.
• Process 4-1 Isobaric Heat Addition: High pressure liquid enters the boiler from the feed pump (4) and is heated to the saturation temperature (4s). Further addition of energy causes evaporation of the liquid until it is fully converted to saturated steam (1s). Then it will be superheated to state 1.

XSteam:

XSteam is an implementation of the IAPWS IF97 standard formulation. It provides accurate data for water and steam and mixtures of water and steam properties from 0 - 1000 bar and from 0 - 2000°C. It is available for MATLAB, Excel, OpenOffice, and as a DLL.

Provided thermodynamic properties are:

• Temperature
• Pressure
• Enthalpy
• Specific volume
• Density
• Specific entropy
• Specific internal energy
• Specific isobaric heat capacity
• Specific isochoric heat capacity
• Speed of sound
• Viscosity
• Vapour fraction

All properties can be calculated with the inputs, p and T known, p and h known, h and s known and some with pressure and density known.

XSteam is a full implementation of the IF-97 formulation including all regions and all backward functions for good calculation speed. The code is speed optimized with pressure and enthalpy as inputs for dynamic simulations.

4.Governing Equations:

1. Heat addition, `q_"in"=h1-h4`
2. Heat rejection, `q_(out)=h2-h3`
3. Work done by the Turbine, `W_T=h1-h2`
4. Work done on the Pump ,`W_P=h3-h4=v_3*(P3-P4)`
5. Net Work Output ,`W_"net"=W_T-W_P`
6. Thermal Efficiency, `eta_"thermal"=W_"net"/q_"in"`
7. Back Work Ratio, `BWR=W_"-ve"/W_"+ve"=W_P/W_T`
8. Specific Steam Consumption, `SSC=3600/W_"net"`

All the parameters mentioned above are specific parameters,i.e., per unit kilogram.

5.Code:

clear all;
close all;
clc;

%Rankine Cycle Simulator using Matlab


%Thermodynamic processes in Rankine cycle
fprintf('%65s\n', '<strong>RANKINE CYCLE SIMULATOR</strong>')
disp(' ')
disp('1-2 is Isentropic Expansion in the Turbine')
disp('2-3 is Constant Pressure Heat Rejection by the Condensor')
disp('3-4 is Isentropic Compression in the Pump')
disp('4-1 is Constant Pressure Heat Addition by the Boiler')
disp(' ')

%Inputs
p1=input('Enter the Pressure at the Turbine inlet(in bar):\n');
T1=input('Enter the Temperature at the Turbine inlet(in Degree Celcius):\n');
p2=input('Enter the Pressure at the Condensor inlet(in bar):\n');
disp(' ')

%checking whether the given rankine cycle is suoperheated or not
T1_s=XSteam('Tsat_p',p1); %if T1>T1_s, it is superheated rankine cycle

%Calculation of Thermodynamic Properties at state points

%State point1 || Inlet of Turbine
h1=XSteam('h_pT',p1,T1); %h is the specific enthalpy in KJ/kg
s1=XSteam('s_pT',p1,T1) ;%s is the specific entropy in KJ/kgK
sg_1=XSteam('sV_p',p1); %Specific entropy of saturated vapour at state point 1
sf_1=XSteam('sL_p',p1); %Specific entropy of saturated liquid at state point 1

%State point 2 || Inlet of condensor
s2=s1; %since process 1-2 is isentropic expansion
sg_2=XSteam('sV_p',p2); %Specific entropy of saturated vapour at state point 2
sf_2=XSteam('sL_p',p2); %Specific entropy of saturated liquid at state point 2

if(s2<sg_2)
x2=(s2-sf_2)/(sg_2-sf_2); %where x2 is the Dryness fraction at state point 2|| s2=sf_2+x2(sg_2-sf_2)
hg_2=XSteam('hV_p',p2); %Specific enthalpy of saturated vapour at state point 2
hf_2=XSteam('hL_p',p2); %Specific enthalpy of saturated liquid at state point 2
h2=hf_2+x2*(hg_2-hf_2);
T2_s=XSteam('Tsat_p',p2);
T2=T2_s;

else x2=1;
Cp2 = XSteam('CpL_p',p2);
T2=(T2_s+273)*exp((s2-sg_2)/Cp2); %Since for solids/liquids ds=mCpln(T(final)/T(initial))
T2=T2-273;

end

%State point 3 || inlet of compressor
T3=T2;
p3=p2;
hf_3=XSteam('hL_T',T3);  %Specific enthalpy of saturated liquid at state point 3
h3=hf_3; 
s3=XSteam('sL_T',T3);

%State point 4 || %inlet of boiler
s4=s3;%process 3-4 is isentropic process
p4=p1; %process 1-4 is constant pressure heat addition process
sf_4=XSteam('sL_p',p4);%Specific entropy of saturated liquid at state point 4
h4=XSteam('h_ps',p4,s4);
Cp4 = XSteam('cpl_p',p4);
T4_s=XSteam('Tsat_p',p4);  
T4=T1_s/exp((sf_4-s4)/Cp4); %Since for solids/liquids ds=mCpln(T(final)/T(initial))

%Calculation of Work,Efficiency,SSC and Back Work Ratio
w_t=h1-h2; %Work done by the Turbine in KJ/kg
v3=XSteam('v_ps',p3,s3); %specific volume at inlet of pump
w_p=v3*(p3-p4)*100; %Work done on the pump in KJ/kg || 1bar=100Kpa
w_p=-w_p;
w_net=w_t-w_p; %Net work required in Rankine cycle in KJ/kg
BWR=w_p/w_t; %Back Work Ratio(BWR)=Negative work/Positive work
SSC=3600/w_net; %in kg/KWh
qin=h1-h4; %Heat input in KJ/kg
Ntherm=(w_net/qin)*100; %Thermal efficiency of Rankine cycle in percentage

%Calculating Liquid and Vapour line for water
temp=linspace(1,374,500);
for i=1:length(temp)
    hL(i)=XSteam('hL_T',temp(i)); %Enthalpy along the liquid line
    hv(i)=XSteam('hv_T',temp(i)); %Enthalpy along the vapour line
    sL(i)=XSteam('sL_T',temp(i)); %Entropy along the liquid line
    sv(i)=XSteam('sv_T',temp(i)); %Entropy along the vapour line
end

%Plotting
%plot of T-S diagram
figure(1)
plot([s1,s2],[T1,T2],'Linewidth',2,'color','g')
hold on
plot([s2,s3],[T2,T3],'LineWidth',2,'Color','b')
plot([s3,s4],[T3,T4],'LineWidth',2,'Color','[0.50,0.24,0.33]')
plot([s4,sf_4],[T4,T4_s],'LineWidth',2,'Color','r')
plot(sL,temp,'K','LineStyle','--','LineWidth',2)
plot(sv,temp,'K','LineStyle','--','LineWidth',2)
plot([sf_4,sg_1],[T4_s,T1_s],'LineWidth',2,'Color','r')
plot([sg_1,s1],[T1_s,T1],'Linewidth',2,'Color','r')
text(s1,T1,'1')
text(s2,T2,'2')
text(s3,T3,'3')
text(s4,T4,'4')
text(sf_4,T4_s,'4s')
text(sg_1,T1_s,'1s')
title('T-s Daigram of Rankine Cycle','FontSize',22)
xlabel('Entropy(KJ/kgK)')
ylabel('Temperature(Degree Celcius)')
legend('Isentropic Expansion in the Turbine','Constant Pressure Heat Rejection in the Condensor','Isentropic Compression in the Pump','Constant Pressure Heat Addition in the Boiler','Saturation Curve')
axis([0 12 0 600])


%plot of h-s Diagram
figure(2)
hold on 
plot([s1,s2],[h1,h2],'LineWidth',2,'Color','g')
plot([s2,s3],[h2,h3],'LineWidth',2,'Color','b')
plot([s3,s4],[h3,h4],'LineWidth',2,'Color','[0.50,0.24,0.33]')
plot([s4,s1],[h4,h1],'LineWidth',2,'Color','r')
plot(sL,hL,'K','LineStyle','--','LineWidth',2)
plot(sv,hv,'K','LineStyle','--','LineWidth',2)
text(s1,h1,'1')
text(s2,h2,'2')
text(s3,h3,'3')
text(s4,h4,'4')
title('h-s Diagram for Rankine Cycle','FontSize',22)
xlabel('Entropy(KJ/kgK')
ylabel('Enthalpy(KJ/kg)')
legend('Isentropic Expansion in the Turbine','Constant Pressure Heat Rejection in the Condensor','Isentropic Compression in the Pump','Constant Pressure Heat Addition in the Boiler','Saturation Curve')
axis([0 12 0 4000])


%Displaying the results in the Command Window
fprintf('%57s\n','<strong>RESULTS</strong>')
disp(' ')
fprintf('At State point1:\n')
fprintf('P1 is : %f Bar\n',p1)
fprintf('T1 is : %f Deg Celcius\n',T1)
fprintf('h1 is : %f KJ/kg\n',h1)
fprintf('s1 is : %f KJ/kgK\n',s1)
disp(' ')
fprintf('At State point2:\n')
fprintf('P2 is : %f Bar\n',p2)
fprintf('T2 is : %f Deg Celcius\n',T2)
fprintf('h2 is : %f KJ/kg\n',h2)
fprintf('s2 is : %f KJ/kgK\n',s2)
fprintf('x2 is : %f\n',x2)
disp(' ')
fprintf('At State point3:\n')
fprintf('P3 is : %f Bar\n',p3)
fprintf('T3 is : %f Deg Celcius\n',T3)
fprintf('h3 is : %f KJ/kg\n',h3)
fprintf('s3 is : %f KJ/kgK\n',s3)
disp(' ')
fprintf('At State point4:\n')
fprintf('P4 is : %f Bar\n',p4)
fprintf('T4 is : %f Deg Celcius\n',T4)
fprintf('h4 is : %f KJ/kg\n',h4)
fprintf('s4 is : %f KJ/kgK\n',s4)
disp(' ')
fprintf('wt is : %f KJ/kg\n',w_t)
fprintf('wp is : %f KJ/kg\n',w_p)
fprintf('wnet is : %f KJ/kg\n',w_net)
fprintf('Ntherm is : %f Percent\n',Ntherm)
fprintf('S.S.C is : %f kg/kWh\n',SSC)
fprintf('Back Work Ratio is : %f\n',BWR)
disp(' ')

Explanation:

1. Firstly, we start by clear all, clc, close all to clear the Workspace, Command Window and close all the plots.
2. Then we use fprintf() and disp() commands to display the Heading as RANKINE CYCLE SIMULATOR followed by the four processes. Here '<strong>RANKINE CYCLE SIMULATOR</strong>' displays the text 'RANKINE CYCLE SIMULATOR' in Bold letters.
3. Now we enter the inputs using the commands
   
   p1=input('Enter the Pressure at the Turbine inlet(in bar):\n');
   T1=input('Enter the Temperature at the Turbine inlet(in Degree Celcius):\n');

   p2=input('Enter the Pressure at the Condensor inlet(in bar):\n');

   disp(' ')

   [Given input values are:

   Turbine Inlet Pressure = 30 bars
   Turbine Inlet temperature = 400° C
   Turbine Outlet Pressure = 0.05 bars = Condenser Inlet Pressure]

   By doing so, what happens is each time Matlab reads an input command, the content within the paranthesis appears in the command window to which we have to assign the corresponding variable.
4. Then we check whether the given cycle is a superheated Rankine cycle or not, by using the XSteam.m function as
   
   T1_s=XSteam('Tsat_p',p1); %if T1>T1_s, it is superheated rankine cycle
5. Then we evaluate the specific thermodynamic properties using the Xsteam function at each of the four state points as well as the saturation curves. While stating the properties at the state point 2, we use the if condition as:
   
   
   if(s2<sg_2)

   x2=(s2-sf_2)/(sg_2-sf_2); %where x2 is the Dryness fraction at state point 2|| s2=sf_2+x2(sg_2-sf_2)

   hg_2=XSteam('hV_p',p2); %Specific enthalpy of saturated vapour at state point 2

   hf_2=XSteam('hL_p',p2); %Specific enthalpy of saturated liquid at state point 2

   h2=hf_2+x2*(hg_2-hf_2);

   T2_s=XSteam('Tsat_p',p2);

   T2=T2_s;

   else x2=1;

   Cp2 = XSteam('CpL_p',p2);

   T2=(T2_s+273)*exp((s2-sg_2)/Cp2); %Since for solids/liquids ds=mCpln(T(final)/T(initial))

   T2=T2-273;

   end

   If the if condition is satisfied,the state point 2 will be having a dryness fraction less than 1, else it will be 1, i.e., saturated vapour state.

   All the enthalpy,entropy,Cp,Dryness fraction values can be found out using the XSteam library.
6. In a similar fashion, we can find the properties at state points 3(inlet of pump) and 4(inlet of boiler).
7. In the next step, we calculate Work done by the Turbine, Work done on the Pump,Net Work Done,Specific Steam Consumption, Back Work Ratio, Thermal Efficiency, Heat added etc. using the Governing Equations mentioned in section4 of this report.
8. Then we calculate Liquid and Vapour line using for loop as shown below.
   temp=linspace(1,374,500);

   for i=1:length(temp)

   hL(i)=XSteam('hL_T',temp(i)); %Enthalpy along the liquid line

   hv(i)=XSteam('hv_T',temp(i)); %Enthalpy along the vapour line

   sL(i)=XSteam('sL_T',temp(i)); %Entropy along the liquid line

   sv(i)=XSteam('sv_T',temp(i)); %Entropy along the vapour line

   end
9. In the next step, we plot the T-s and h-s diagrams. We plot the figures using plot() command,hold on command is used so that all the plots can be plotted in the same diagram, text() command is used to name the state points, title() command is used to give the title, xlabel() and ylabel() are used to label the axes. Finally we use the legend() command and define the axes limit using the axis([xi xf yi yf]) command as shown in the code.
10. Now we display the 'RESULTS' in the Command Window using the disp() and printf() commands.

6.Output:

6.1Output Window:

6.2Command Window:

6.3Workspace:

7.Graph:

7.1 T-s Diagram:

7.2 h-s Diagram:

This is enlarged view of process 3-4 of h-s diagram.

8.Results:

1. Simulated Rankine Cycle using Matlab.
2. Calculated the state points of Rankine Cycle based on user inputs and displayed them in the Command Window.
3. Calculated Net Work Output and Back Work Ratio.
   Net Work Output=1117.9 KJ/kg
   Back Work Ratio=0.0027
4. Created a plot of T-s and h-s diagrams for the given set of inputs.

9.Reference:

1. Basic and Applied Thermodynamics, P K Nag.
2. https://www.thermopedia.com/content/1072/
3. https://in.mathworks.com/matlabcentral/fileexchange/9817-x-steam-thermodynamic-properties-of-water-and-steam