Ideal Rankine Cycle simulator using MATLAB
I. AIM:
The aim of this project is to develop an ideal Rankine cycle simulator using MATLAB programming.
A. Introduction:
The Rankine cycle is an ideal cycle for vapor power plants in which the working fluid is alternately vaporised and condensed. The term ideal implies that the cycle doesn’t involve any internal irreversibilities. The working fluid is generally water and its vapor (steam). The schematic diagram and the T-S diagram of a simple (ideal) Rankine cycle is shown in the below figure:
Schematic of Ideal Rankine cycle (Reference: Ron Hugo, https://www.youtube.com/channel/UCrRKDc1uHcif575CQexIA5Q)
It is a thermodynamic cycle consisting of four processes namely:
1-2 Isentropic compression in a pump
2-3 Constant pressure or Isobaric heat addition in a boiler
3-4 Isentropic expansion in a turbine
4-1 Constant pressure or Isobaric heat rejection in a condenser
1. Process 1-2:
Water enters the pump at state 1 as saturated liquid at condenser pressure and then it is isentropically compressed to the boiler operating pressure. During this process the pump draws the energy from a prime mover, which is the work input for the cycle. Since, water is an incompressible liquid the temperature and enthalpy rise during the compression is of negligible value.
2. Process 2-3:
The compressed water from the pump enters the boiler at state2, where heat is added to it and leaves the boiler at state3. The boiler is a large heat exchanger where heat from combustion gases/nuclear reactor or other sources is used to convert the liquid water to steam at constant pressure. The boiler, together with the section where the steam is superheated (the super heater), is often called the steam generator.
3. Process 3-4:
The superheated vapor at state 3 enters the turbine, where it expands isentropically and produces work by rotating the shaft connected to an electric generator. The pressure and the temperature of steam drop during this process to the values at state 4, where steam enters the condenser.
4.Process 4-1:
During this process, steam at saturated liquid–vapor mixture with a high quality at state 4 is condensed at constant pressure in the condenser, which is basically a large heat exchanger, by rejecting heat to a cooling medium such as a lake, a river, or the atmosphere. Steam leaves the condenser as saturated liquid and enters the pump at state1, completing the cycle.
The area under the curve 2-3 on the T-S represent the heat added to the system and the area under the curve 2-3 on the T-S represent the heat rejected by the system. The difference between these two (the area enclosed by the cycle curve) is the net work (turbine) produced during the cycle.
B. Energy Analysis of the Ideal Rankine Cycle
The components present in the Rankine cycle (the pump, boiler, turbine, and condenser) are all steady-flow devices. Hence the cycle can be analysed using steady flow analysis. The kinetic and potential energies changes of the working fluid are very small compared to the work and heat transfer during the cycle. Hence these can be neglected.
The steady flow energy per unit mass can be written as:
`( q_(In) − q_(Out) ) + ( w_(In) −w_(Out) ) = h_e − h_i (kJ)/(kg)`
There is no work done in the boiler and condenser. Also, the pump and turbine are assumed to be isentropic. Then the conservation of energy relation for each device can be expressed as follows:
Pump (q = 0):
`w_( pump,In) = h _(2 )− h_1 = v(p2-p1)`
`h _1 = h_(f) (p1) and v ≅ v 1 = v _(f) (p1)`
Boiler ( w = 0):
`q _(In) = h_( 3) − h_( 2)`
Turbine ( q = 0):
`w_(turb,out) = h_3 − h_4`
Condenser ( w = 0):
`q_(out) = h_4 − h_1`
The thermal efficiency of the Rankine cycle is determined from:
`η_(th) = w_(Net)/q _(In)`
where, `w_(Net) = q_(In)− q_(Out) = w_(turb,out) − w_(pump,In)`
The thermal efficiency can also be interpreted as the ratio of the area enclosed by the cycle on a T-s diagram to the area under the heat-addition process. The analysis of the Rankine cycle can be done using steam tables, which provide information about state properties.
II. Objectives:
- Create an ideal Rankine Cycle Simulator using MATLAB.
- Calculate the state points of the Rankine Cycle based on user inputs.
- Plot the corresponding T-s and h-s plots for the given set of inputs.
III. Program:
The program for ideal Rankine simulator is given in the code below:
%program to simulate the ideal rankine cycle of a steam turbine
clear all
close all
clc
fprintf(' nn IDEAL RANKINE CYCLE SIMULATOR')
fprintf(' n n Process 1-2 is isentropic compression in the pumpn')
fprintf(' Process 2-3 is constant pressure heat addition in the boilern')
fprintf(' Process 3-4 is isentropic expansion in the steam turbinen')
fprintf(' Process 4-1 is constant pressure heat rejection in condensernn')
%Inputs from the user
fprintf(' Inputs requiredn')
p3 = input('n 1. Enter the pressure at turbine inlet (in bar) = ')
T3 = input(' n 2. Enter the temperature at turbine inlet (in degree celcius) = ')
p4 = input('n 3. Enter the pressure at turbine outlet/condenser inlet(in bar) = ')
%computing all the state properties
%state1 properties
p1 = p4 ; %Process 4-1 is constant pressure heat rejection
T1 = XSteam('Tsat_p',p1); %temperature at state1 1 is equal to saturated liquid enthalpy at constant pressure state 1 using Xsteam function
h1 = XSteam('hL_p',p1); %enthalpy at state1 1 is equal to saturated liquid enthalpy at constant pressure state 1 using Xsteam function
s1 = XSteam('sL_p',p1); %entropy at state 1 is equal to saturated liquid enthalpy at constant pressure state 1 using Xsteam function
v1 = XSteam('vL_p',p1) ; %liquid volume at 1 is equal to saturated liquid volume at constant pressure state 1 using Xsteam function
sf1 = XSteam('sL_p',p1);
%state2 properties
p2 = p3 ; %Process 2-3 is constant pressure heat addition
s2 = s1 ; %process 1-2 is isentropic compression
T2 = XSteam('T_ps',p2,s2); %temperature at state1 2 is calculated using XSteam function with inputs arguments of p2 and s2
v2 = XSteam('v_pT',p2,s2) ; %specific volume at state 2 is calculated using XSteam function with inputs arguments of p2 and s2
h2 = XSteam('h_ps',p2,s2); %enthalpy at state 2 is calculated using XSteam function with inputs arguments of p2 and s2
Tsat_2 = XSteam('Tsat_p',p2); %Saturation temperature at p2
sf2 = XSteam('sL_T',Tsat_2); %Saturated liquid entropy at p2
hf2 = XSteam('hL_T',Tsat_2); %Saturated liquid Enthaply at p2
sg2 = XSteam('sV_T',Tsat_2); %Saturated vapor entropy at p2
hg2 = XSteam('hV_T',Tsat_2); %Saturated vapor Enthaply at p2
Cp = XSteam('Cp_ps',p2,s2); %specific heat at constant pressure at state 2
Cv = XSteam('Cv_ps',p2,s2); %specific heat at constant volume at state 2
C = Cp/Cv; %ratio of specific heats
%Pump work in (kJ/kg)
w_p = h2-h1 ; %pump work in terms of enthalpy change during process 1-2 is Wp = h2-h1
%state3 properties
h3 = XSteam('h_pt',p3,T3); %enthalpy at state 3 is calculated using XSteam function with inputs arguments of p3 and T2
s3 = XSteam('s_pT',p3,T3);
v3 = XSteam('v_pT',p3,T3); %specific volume at state3 is equal to vapour volume at p3 and T3 using XSteam function with inputs p3 ,T3
x3 = XSteam('vx_ph',p3,h3) ; %Vapour volume fraction at p3
%state4 properties
s4 = s3; %process 3-4 is isentropic expansion
x4 = XSteam('x_ps',p4,s4); %vapour fraction/dryness fraction at state4
h4 = XSteam('h_ps',p4,s4); %enthalpy at state 4
T4 = XSteam('T_ph',p4,h4); %temperature at state4
sg4 = XSteam('sV_p',p4); %Saturated vapor entropy at p4
%heat added in boiler
q_in = h3-h2;
%heat rejected in condenser
q_out = h4-h1;
%work output of turbine
w_t = h3 -h4;
%net work output
w_net = w_t-w_p;
%thermal efficiency of the cycle
Eta_th = (w_net/q_in)*100 ;
%back work ratio
BWR = w_p/w_t; %it is the ratio of pump work to turbine work
%vapor-dome values
t_s = linspace(5,400,500);
for i =1:length(t_s)
s(i) = XSteam('sL_T',t_s(i));
S(i) = XSteam('sV_T',t_s(i));
h(i) = XSteam('hL_T',t_s(i));
H(i) = XSteam('hV_T',t_s(i));
end
%Temperature-Entropy diagram
figure(1)
hold on
line([s1 s2],[T1 T2],'color','r','Marker','x','linewidth',2,'Markersize',8,'MarkerEdgeColor','k')
line([s2 sf2],[T2 Tsat_2],'color','b','Marker','x','linewidth',2,'Markersize',8,'MarkerEdgeColor','k')
line([s3 s4],[T3 T4],'color','g','Marker','x','linewidth',2,'Markersize',8,'MarkerEdgeColor','k')
line([s4 s1],[T4 T1],'color','m','Marker','x','linewidth',2,'Markersize',8,'MarkerEdgeColor','k')
line([sf2 sg2],[Tsat_2 Tsat_2],'color','b','Marker','x','linewidth',2,'Markersize',8,'MarkerEdgeColor','k')
line([sg2 s3],[Tsat_2 T3],'color','b','Marker','x','linewidth',2,'Markersize',8,'MarkerEdgeColor','k')
plot(s,t_s,'color','black');
plot(S,t_s,'color','black');
xlabel('Entropy [kJ/kg]')
ylabel('Temperature [deg-celcius]')
title('T-s Diagram of ideal Rankine cycle')
legend('isentropic compression in pump','isobaric heating in boiler','isentropic expansion in turbine','isobaric heat rejection in condenser')
%enthaply-entropy diagram
figure(2)
hold on
line([s1 s2],[h1 h2],'color','r','Marker','x','linewidth',2,'Markersize',8,'MarkerEdgeColor','k')
line([s2 sf2],[h2 hf2],'color','b','Marker','x','linewidth',2,'Markersize',8,'MarkerEdgeColor','k')
line([s3 s4],[h3 h4],'color','g','Marker','x','linewidth',2,'Markersize',8,'MarkerEdgeColor','k')
line([s4 s1],[h4 h1],'color','m','Marker','x','linewidth',2,'Markersize',8,'MarkerEdgeColor','k')
line([sf2 sg2],[hf2 hg2],'color','b','Marker','x','linewidth',2,'Markersize',8,'MarkerEdgeColor','k')
line([sg2 s3],[hg2 h3],'color','b','Marker','x','linewidth',2,'Markersize',8,'MarkerEdgeColor','k')
plot(s,h,'color','black');
plot(S,H,'color','black');
xlabel('Entropy [kJ/kg]')
ylabel('Enthalpy [kJ/kg]')
title('h-s Diagram of ideal Rankine cycle')
legend('isentropic compression in pump','isobaric heating in boiler','isentropic expansion in turbine','isobaric heat rejection in condenser')
%output from the simulator
fprintf('n Output:n')
fprintf('n A. State1(Pump inlet) properties:')
fprintf('n 1. Pressure= %d bar',p1)
fprintf('n 2. Temperature = %d K',T1)
fprintf('n 4. Enthalpy = %d kJ/kg',h1)
fprintf('n 5. Entropy = %d kJ/kg',s1)
fprintf('nn B. State2(Boiler inlet) properties:')
fprintf('n 1. Pressure= %d bar',p2)
fprintf('n 2. Temperature = %d K',T2)
fprintf('n 3. Enthalpy = %d kJ/kg',h2)
fprintf('n 4. Entropy = %d kJ/kg',s2)
fprintf('nn C. properties at Boiler pressure Saturation:')
fprintf('n 1. Saturation Temperature= %d K',Tsat_2)
fprintf('n 2. Saturation liquid enthalpy = %d kJ/kg',hf2)
fprintf('n 3. Saturation liquid entropy = %d kJ/kg',sf2)
fprintf('n 2. Saturation vapor enthalpy = %d kJ/kg',hg2)
fprintf('n 3. Saturation vapor entropy = %d kJ/kg',sg2)
fprintf('nn D. State3(Turbine inlet) properties:')
fprintf('n 1. Pressure= %d bar',p3)
fprintf('n 2. Temperature = %d K',T3)
fprintf('n 3. Enthalpy = %d kJ/kg',h3)
fprintf('n 4. Entropy = %d kJ/kg',s3)
fprintf('n 5. Vapor volume fraction = %d ',x3)
fprintf('nn E. State4(condenser inlet) properties:')
fprintf('n 1. Pressure= %d bar',p4)
fprintf('n 2. Temperature = %d K',T4)
fprintf('n 3. Enthalpy = %d kJ/kg',h4)
fprintf('n 4. Entropy = %d kJ/kg',s4)
fprintf('n 5. Vapor volume fraction = %d ',x4)
fprintf('nn F. Cycle Results:')
fprintf('n 1. Pump Work = %d kJ/kg',w_p)
fprintf('n 2. Heat_boiler = %d kJ/kg',q_in)
fprintf('n 3. Turbine_work = %d kJ/kg',w_t)
fprintf('n 4. Net work output = %d kJ/kg',w_net)
fprintf('n 5. Thermal Efficiency = %d',Eta_th)
fprintf('n 6. Back work ratio = %d',BWR)
fprintf('n ')
Code Explanation:
Step1: clear all – command to clear the workspace.
Step 2: close all – command to close all opened previous plots.
Step 3: clc – command to clear the command window.
Step4: The first step is to print the header content containing the title of the program and the four processes involved in the ideal Rankine cycle into the command window using ‘fprintf’ command.
Step 5: Inputs are taken from the user for different state properties available using the ‘input’ command.
Step 6:
Computing the state properties of pressure, temperature, enthalpy, entropy, specific volume, saturation properties, vapour fraction, specific heats at constant pressure and constant volume, ratio of specific heats using the steam tables, which are taken from the function code ‘XSteam’.
The state properties are calculated at
1: pump inlet
- Boiler inlet
- Turbine inlet
- Condenser inlet
The function code ‘XSteam’ containing the steam tables library for water is given below:
%h_prho behöver T_prho för samtliga regioner!!!!
%***********************************************************************************************************
%* Water and steam properties according to IAPWS IF-97 *
%* By Magnus Holmgren, www.x-eng.com *
%* The steam tables are free and provided as is. *
%* We take no responsibilities for any errors in the code or damage thereby. *
%* You are free to use, modify and distribute the code as long as authorship is properly acknowledged. *
%* Please notify me at [email protected] if the code is used in commercial applications *
%***********************************************************************************************************
%
% XSteam provides accurate steam and water properties from 0 - 1000 bar and from 0 - 2000 deg C according to
% the standard IAPWS IF-97. For accuracy of the functions in different regions see IF-97 (www.iapws.org)
%
% *** Using XSteam *****************************************************************************************
%XSteam take 2 or 3 arguments. The first argument must always be the steam table function you want to use.
%The other arguments are the inputs to that function.
%Example: XSteam('h_pt',1,20) Returns the enthalpy of water at 1 bar and 20 degC
%Example: XSteam('TSat_p',1) Returns the saturation temperature of water at 1 bar.
%For a list of valid Steam Table functions se bellow or the XSteam macros for MS Excel.
%
%*** Nomenclature ******************************************************************************************
% First the wanted property then a _ then the wanted input properties.
% Example. T_ph is temperature as a function of pressure and enthalpy.
% For a list of valid functions se bellow or XSteam for MS Excel.
% T Temperature (deg C)
% p Pressure (bar)
% h Enthalpy (kJ/kg)
% v Specific volume (m3/kg)
% rho Density
% s Specific entropy
% u Specific internal energy
% Cp Specific isobaric heat capacity
% Cv Specific isochoric heat capacity
% w Speed of sound
% my Viscosity
% tc Thermal Conductivity
% st Surface Tension
% x Vapour fraction
% vx Vapour Volume Fraction
%
%*** Valid Steam table functions. ****************************************************************************
%
%Temperature
%Tsat_p Saturation temperature
%T_ph Temperture as a function of pressure and enthalpy
%T_ps Temperture as a function of pressure and entropy
%T_hs Temperture as a function of enthalpy and entropy
%
%Pressure
%psat_T Saturation pressure
%p_hs Pressure as a function of h and s.
%p_hrho Pressure as a function of h and rho. Very unaccurate for solid water region since it's almost incompressible!
%
%Enthalpy
%hV_p Saturated vapour enthalpy
%hL_p Saturated liquid enthalpy
%hV_T Saturated vapour enthalpy
%hL_T Saturated liquid enthalpy
%h_pT Entalpy as a function of pressure and temperature.
%h_ps Entalpy as a function of pressure and entropy.
%h_px Entalpy as a function of pressure and vapour fraction
%h_prho Entalpy as a function of pressure and density. Observe for low temperatures (liquid) this equation has 2 solutions.
%h_Tx Entalpy as a function of temperature and vapour fraction
%
%Specific volume
%vV_p Saturated vapour volume
%vL_p Saturated liquid volume
%vV_T Saturated vapour volume
%vL_T Saturated liquid volume
%v_pT Specific volume as a function of pressure and temperature.
%v_ph Specific volume as a function of pressure and enthalpy
%v_ps Specific volume as a function of pressure and entropy.
%
%Density
%rhoV_p Saturated vapour density
%rhoL_p Saturated liquid density
%rhoV_T Saturated vapour density
%rhoL_T Saturated liquid density
%rho_pT Density as a function of pressure and temperature.
%rho_ph Density as a function of pressure and enthalpy
%rho_ps Density as a function of pressure and entropy.
%
%Specific entropy
%sV_p Saturated vapour entropy
%sL_p Saturated liquid entropy
%sV_T Saturated vapour entropy
%sL_T Saturated liquid entropy
%s_pT Specific entropy as a function of pressure and temperature (Returns saturated vapour entalpy if mixture.)
%s_ph Specific entropy as a function of pressure and enthalpy
%
%Specific internal energy
%uV_p Saturated vapour internal energy
%uL_p Saturated liquid internal energy
%uV_T Saturated vapour internal energy
%uL_T Saturated liquid internal energy
%u_pT Specific internal energy as a function of pressure and temperature.
%u_ph Specific internal energy as a function of pressure and enthalpy
%u_ps Specific internal energy as a function of pressure and entropy.
%
%Specific isobaric heat capacity
%CpV_p Saturated vapour heat capacity
%CpL_p Saturated liquid heat capacity
%CpV_T Saturated vapour heat capacity
%CpL_T Saturated liquid heat capacity
%Cp_pT Specific isobaric heat capacity as a function of pressure and temperature.
%Cp_ph Specific isobaric heat capacity as a function of pressure and enthalpy
%Cp_ps Specific isobaric heat capacity as a function of pressure and entropy.
%
%Specific isochoric heat capacity
%CvV_p Saturated vapour isochoric heat capacity
%CvL_p Saturated liquid isochoric heat capacity
%CvV_T Saturated vapour isochoric heat capacity
%CvL_T Saturated liquid isochoric heat capacity
%Cv_pT Specific isochoric heat capacity as a function of pressure and temperature.
%Cv_ph Specific isochoric heat capacity as a function of pressure and enthalpy
%Cv_ps Specific isochoric heat capacity as a function of pressure and entropy.
%
%Speed of sound
%wV_p Saturated vapour speed of sound
%wL_p Saturated liquid speed of sound
%wV_T Saturated vapour speed of sound
%wL_T Saturated liquid speed of sound
%w_pT Speed of sound as a function of pressure and temperature.
%w_ph Speed of sound as a function of pressure and enthalpy
%w_ps Speed of sound as a function of pressure and entropy.
%
%Viscosity
%Viscosity is not part of IAPWS Steam IF97. Equations from
%"Revised Release on the IAPWS Formulation 1985 for the Viscosity of Ordinary Water Substance", 2003 are used.
%Viscosity in the mixed region (4) is interpolated according to the density. This is not true since it will be two fases.
%my_pT Viscosity as a function of pressure and temperature.
%my_ph Viscosity as a function of pressure and enthalpy
%my_ps Viscosity as a function of pressure and entropy.
%
%Thermal Conductivity
%Revised release on the IAPS Formulation 1985 for the Thermal Conductivity of ordinary water substance (IAPWS 1998)
%tcL_p Saturated vapour thermal conductivity
%tcV_p Saturated liquid thermal conductivity
%tcL_T Saturated vapour thermal conductivity
%tcV_T Saturated liquid thermal conductivity
%tc_pT Thermal conductivity as a function of pressure and temperature.
%tc_ph Thermal conductivity as a function of pressure and enthalpy
%tc_hs Thermal conductivity as a function of enthalpy and entropy
%
%Surface tension
%st_T Surface tension for two phase water/steam as a function of T
%st_p Surface tension for two phase water/steam as a function of T
%Vapour fraction
%x_ph Vapour fraction as a function of pressure and enthalpy
%x_ps Vapour fraction as a function of pressure and entropy.
%
%Vapour volume fraction
%vx_ph Vapour volume fraction as a function of pressure and enthalpy
%vx_ps Vapour volume fraction as a function of pressure and entropy.
function Out=XSteam(fun,In1,In2)
%*Contents.
%*1 Calling functions
%*1.1
%*1.2 Temperature (T)
%*1.3 Pressure (p)
%*1.4 Enthalpy (h)
%*1.5 Specific Volume (v)
%*1.6 Density (rho)
%*1.7 Specific entropy (s)
%*1.8 Specific internal energy (u)
%*1.9 Specific isobaric heat capacity (Cp)
%*1.10 Specific isochoric heat capacity (Cv)
%*1.11 Speed of sound
%*1.12 Viscosity
%*1.13 Prandtl
%*1.14 Kappa
%*1.15 Surface tension
%*1.16 Heat conductivity
%*1.17 Vapour fraction
%*1.18 Vapour Volume Fraction
%
%*2 IAPWS IF 97 Calling functions
%*2.1 Functions for region 1
%*2.2 Functions for region 2
%*2.3 Functions for region 3
%*2.4 Functions for region 4
%*2.5 Functions for region 5
%
%*3 Region Selection
%*3.1 Regions as a function of pT
%*3.2 Regions as a function of ph
%*3.3 Regions as a function of ps
%*3.4 Regions as a function of hs
%
%4 Region Borders
%4.1 Boundary between region 1 and 3.
%4.2 Region 3. pSat_h and pSat_s
%4.3 Region boundary 1to3 and 3to2 as a functions of s
%
%5 Transport properties
%5.1 Viscosity (IAPWS formulation 1985)
%5.2 Thermal Conductivity (IAPWS formulation 1985)
%5.3 Surface Tension
%
%6 Units
%
%7 Verification
%7.1 Verifiy region 1
%7.2 Verifiy region 2
%7.3 Verifiy region 3
%7.4 Verifiy region 4
%7.5 Verifiy region 5
%***********************************************************************************************************
%*1 Calling functions *
%***********************************************************************************************************
%***********************************************************************************************************
%*1.1
fun=lower(fun);
switch fun
%***********************************************************************************************************
%*1.2 Temperature
case 'tsat_p'
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
Out = fromSIunit_T(T4_p(p));
else
Out = NaN;
end
case 'tsat_s'
s = toSIunit_s(In1);
if s > -0.0001545495919 && s < 9.155759395
ps = p4_s(s);
Out = fromSIunit_T(T4_p(ps));
else
Out = NaN;
end
case 't_ph'
p = toSIunit_p(In1);
h = toSIunit_h(In2);
Region = region_ph(p, h);
switch Region
case 1
Out = fromSIunit_T(T1_ph(p, h));
case 2
Out = fromSIunit_T(T2_ph(p, h));
case 3
Out = fromSIunit_T(T3_ph(p, h));
case 4
Out = fromSIunit_T(T4_p(p));
case 5
Out = fromSIunit_T(T5_ph(p, h));
otherwise
Out = NaN;
end
case 't_ps'
p = toSIunit_p(In1);
s = toSIunit_s(In2);
Region = region_ps(p, s);
switch Region
case 1
Out = fromSIunit_T(T1_ps(p, s));
case 2
Out = fromSIunit_T(T2_ps(p, s));
case 3
Out = fromSIunit_T(T3_ps(p, s));
case 4
Out = fromSIunit_T(T4_p(p));
case 5
Out = fromSIunit_T(T5_ps(p, s));
otherwise
Out = NaN;
end
case 't_hs'
h = toSIunit_h(In1);
s = toSIunit_s(In2);
Region = region_hs(h, s);
switch Region
case 1
p1 = p1_hs(h, s);
Out = fromSIunit_T(T1_ph(p1, h));
case 2
p2 = p2_hs(h, s);
Out = fromSIunit_T(T2_ph(p2, h));
case 3
p3 = p3_hs(h, s);
Out = fromSIunit_T(T3_ph(p3, h));
case 4
Out = fromSIunit_T(T4_hs(h, s));
case 5
error('functions of hs is not avlaible in region 5');
otherwise
Out = NaN;
end
%***********************************************************************************************************
%*1.3 Pressure (p)
case 'psat_t'
T = toSIunit_T(In1);
if T < 647.096 && T > 273.15
Out = fromSIunit_p(p4_T(T));
else
Out = NaN;
end
case 'psat_s'
s = toSIunit_s(In1);
if s > -0.0001545495919 && s < 9.155759395
Out = fromSIunit_p(p4_s(s));
else
Out = NaN;
end
case 'p_hs'
h = toSIunit_h(In1);
s = toSIunit_s(In2);
Region = region_hs(h, s);
switch Region
case 1
Out = fromSIunit_p(p1_hs(h, s));
case 2
Out = fromSIunit_p(p2_hs(h, s));
case 3
Out = fromSIunit_p(p3_hs(h, s));
case 4
tSat = T4_hs(h, s);
Out = fromSIunit_p(p4_T(tSat));
case 5
error('functions of hs is not avlaible in region 5');
otherwise
Out = NaN;
end
case 'p_hrho'
h=In1;
rho=In2;
%Not valid for water or sumpercritical since water rho does not change very much with p.
%Uses iteration to find p.
High_Bound = fromSIunit_p(100);
Low_Bound = fromSIunit_p(0.000611657);
ps = fromSIunit_p(10);
rhos = 1 / XSteam('v_ph',ps, h);
while abs(rho - rhos) > 0.0000001
rhos = 1 / XSteam('v_ph',ps, h);
if rhos >= rho
High_Bound = ps;
else
Low_Bound = ps;
end
ps = (Low_Bound + High_Bound) / 2;
end
Out = ps;
%***********************************************************************************************************
%*1.4 Enthalpy (h)
case 'hv_p'
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
Out = fromSIunit_h(h4V_p(p));
else
Out = NaN;
end
case 'hl_p'
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
Out = fromSIunit_h(h4L_p(p));
else
Out = NaN;
end
case 'hv_t'
T = toSIunit_T(In1);
if T > 273.15 && T < 647.096
p = p4_T(T);
Out = fromSIunit_h(h4V_p(p));
else
Out = NaN;
end
case 'hl_t'
T = toSIunit_T(In1);
if T > 273.15 && T < 647.096
p = p4_T(T);
Out = fromSIunit_h(h4L_p(p));
else
Out = NaN;
end
case 'h_pt'
p = toSIunit_p(In1);
T = toSIunit_T(In2);
Region = region_pT(p, T);
switch Region
case 1
Out = fromSIunit_h(h1_pT(p, T));
case 2
Out = fromSIunit_h(h2_pT(p, T));
case 3
Out = fromSIunit_h(h3_pT(p, T));
case 4
Out = NaN;
case 5
Out = fromSIunit_h(h5_pT(p, T));
otherwise
Out = NaN;
end
case 'h_ps'
p = toSIunit_p(In1);
s = toSIunit_s(In2);
Region = region_ps(p, s);
switch Region
case 1
Out = fromSIunit_h(h1_pT(p, T1_ps(p, s)));
case 2
Out = fromSIunit_h(h2_pT(p, T2_ps(p, s)));
case 3
Out = fromSIunit_h(h3_rhoT(1 / v3_ps(p, s), T3_ps(p, s)));
case 4
xs = x4_ps(p, s);
Out = fromSIunit_h(xs * h4V_p(p) + (1 - xs) * h4L_p(p));
case 5
Out = fromSIunit_h(h5_pT(p, T5_ps(p, s)));
otherwise
Out = NaN;
end
case 'h_px'
p = toSIunit_p(In1);
x = toSIunit_x(In2);
if x > 1 || x < 0 || p >= 22.064
Out = NaN;
return
end
hL = h4L_p(p);
hV = h4V_p(p);
Out = hL + x * (hV - hL);
case 'h_prho'
p = toSIunit_p(In1);
rho = 1 / toSIunit_v(1 / In2);
Region = Region_prho(p, rho);
switch Region
case 1
Out = fromSIunit_h(h1_pT(p, T1_prho(p, rho)));
case 2
Out = fromSIunit_h(h2_pT(p, T2_prho(p, rho)));
case 3
Out = fromSIunit_h(h3_rhoT(rho, T3_prho(p, rho)));
case 4
if p < 16.529
vV = v2_pT(p, T4_p(p));
vL = v1_pT(p, T4_p(p));
else
vV = v3_ph(p, h4V_p(p));
vL = v3_ph(p, h4L_p(p));
end
hV = h4V_p(p);
hL = h4L_p(p);
x = (1 / rho - vL) / (vV - vL);
Out = fromSIunit_h((1 - x) * hL + x * hV);
case 5
Out = fromSIunit_h(h5_pT(p, T5_prho(p, rho)));
otherwise
Out = NaN;
end
case 'h_tx'
T = toSIunit_T(In1);
x = toSIunit_x(In2);
if x > 1 || x < 0 || T >= 647.096
Out = NaN;
return
end
p = p4_T(T);
hL = h4L_p(p);
hV = h4V_p(p);
Out = hL + x * (hV - hL);
%***********************************************************************************************************
%*1.5 Specific Volume (v)
case {'vv_p','rhov_p'}
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
if p < 16.529
Out = fromSIunit_v(v2_pT(p, T4_p(p)));
else
Out = fromSIunit_v(v3_ph(p, h4V_p(p)));
end
else
Out = NaN;
end
if fun(1)=='r'
Out=1/Out;
end
case {'vl_p','rhol_p'}
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
if p < 16.529
Out = fromSIunit_v(v1_pT(p, T4_p(p)));
else
Out = fromSIunit_v(v3_ph(p, h4L_p(p)));
end
else
Out = NaN;
end
if fun(1)=='r'
Out=1/Out;
end
case {'vv_t','rhov_t'}
T = toSIunit_T(In1);
if T > 273.15 && T < 647.096
if T <= 623.15
Out = fromSIunit_v(v2_pT(p4_T(T), T));
else
Out = fromSIunit_v(v3_ph(p4_T(T), h4V_p(p4_T(T))));
end
else
Out = NaN;
end
if fun(1)=='r'
Out=1/Out;
end
case {'vl_t','rhol_t'}
T = toSIunit_T(In1);
if T > 273.15 && T < 647.096
if T <= 623.15
Out = fromSIunit_v(v1_pT(p4_T(T), T));
else
Out = fromSIunit_v(v3_ph(p4_T(T), h4L_p(p4_T(T))));
end
else
Out = NaN;
end
if fun(1)=='r'
Out=1/Out;
end
case {'v_pt','rho_pt'}
p = toSIunit_p(In1);
T = toSIunit_T(In2);
Region = region_pT(p, T);
switch Region
case 1
Out = fromSIunit_v(v1_pT(p, T));
case 2
Out = fromSIunit_v(v2_pT(p, T));
case 3
Out = fromSIunit_v(v3_ph(p, h3_pT(p, T)));
case 4
Out = NaN;
case 5
Out = fromSIunit_v(v5_pT(p, T));
otherwise
Out = NaN;
end
if fun(1)=='r'
Out=1/Out;
end
case {'v_ph','rho_ph'}
p = toSIunit_p(In1);
h = toSIunit_h(In2);
Region = region_ph(p, h);
switch Region
case 1
Out = fromSIunit_v(v1_pT(p, T1_ph(p, h)));
case 2
Out = fromSIunit_v(v2_pT(p, T2_ph(p, h)));
case 3
Out = fromSIunit_v(v3_ph(p, h));
case 4
xs = x4_ph(p, h);
if p < 16.529
v4v = v2_pT(p, T4_p(p));
v4L = v1_pT(p, T4_p(p));
else
v4v = v3_ph(p, h4V_p(p));
v4L = v3_ph(p, h4L_p(p));
end
Out = fromSIunit_v((xs * v4v + (1 - xs) * v4L));
case 5
Ts = T5_ph(p, h);
Out = fromSIunit_v(v5_pT(p, Ts));
otherwise
Out = NaN;
end
if fun(1)=='r'
Out=1/Out;
end
case {'v_ps','rho_ps'}
p = toSIunit_p(In1);
s = toSIunit_s(In2);
Region = region_ps(p, s);
switch Region
case 1
Out = fromSIunit_v(v1_pT(p, T1_ps(p, s)));
case 2
Out = fromSIunit_v(v2_pT(p, T2_ps(p, s)));
case 3
Out = fromSIunit_v(v3_ps(p, s));
case 4
xs = x4_ps(p, s);
if p < 16.529
v4v = v2_pT(p, T4_p(p));
v4L = v1_pT(p, T4_p(p));
else
v4v = v3_ph(p, h4V_p(p));
v4L = v3_ph(p, h4L_p(p));
end
Out = fromSIunit_v((xs * v4v + (1 - xs) * v4L));
case 5
Ts = T5_ps(p, s);
Out = fromSIunit_v(v5_pT(p, Ts));
otherwise
Out = NaN;
end
if fun(1)=='r'
Out=1/Out;
end
%***********************************************************************************************************
%*1.6 Density (rho)
% Density is calculated as 1/v. Se section 1.5 Volume
%***********************************************************************************************************
%*1.7 Specific entropy (s)
case 'sv_p'
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
if p < 16.529
Out = fromSIunit_s(s2_pT(p, T4_p(p)));
else
Out = fromSIunit_s(s3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)));
end
else
Out = NaN;
end
case 'sl_p'
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
if p < 16.529
Out = fromSIunit_s(s1_pT(p, T4_p(p)));
else
Out = fromSIunit_s(s3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)));
end
else
Out = NaN;
end
case 'sv_t'
T = toSIunit_T(In1);
if T > 273.15 && T < 647.096
if T <= 623.15
Out = fromSIunit_s(s2_pT(p4_T(T), T));
else
Out = fromSIunit_s(s3_rhoT(1 / (v3_ph(p4_T(T), h4V_p(p4_T(T)))), T));
end
else
Out = NaN;
end
case 'sl_t'
T = toSIunit_T(In1);
if T > 273.15 && T < 647.096
if T <= 623.15
Out = fromSIunit_s(s1_pT(p4_T(T), T));
else
Out = fromSIunit_s(s3_rhoT(1 / (v3_ph(p4_T(T), h4L_p(p4_T(T)))), T));
end
else
Out = NaN;
end
case 's_pt'
p = toSIunit_p(In1);
T = toSIunit_T(In2);
Region = region_pT(p, T);
switch Region
case 1
Out = fromSIunit_s(s1_pT(p, T));
case 2
Out = fromSIunit_s(s2_pT(p, T));
case 3
hs = h3_pT(p, T);
rhos = 1 / v3_ph(p, hs);
Out = fromSIunit_s(s3_rhoT(rhos, T));
case 4
Out = NaN;
case 5
Out = fromSIunit_s(s5_pT(p, T));
otherwise
Out = NaN;
end
case 's_ph'
p = toSIunit_p(In1);
h = toSIunit_h(In2);
Region = region_ph(p, h);
switch Region
case 1
T = T1_ph(p, h);
Out = fromSIunit_s(s1_pT(p, T));
case 2
T = T2_ph(p, h);
Out = fromSIunit_s(s2_pT(p, T));
case 3
rhos = 1 / v3_ph(p, h);
Ts = T3_ph(p, h);
Out = fromSIunit_s(s3_rhoT(rhos, Ts));
case 4
Ts = T4_p(p);
xs = x4_ph(p, h);
if p < 16.529
s4v = s2_pT(p, Ts);
s4L = s1_pT(p, Ts);
else
v4v = v3_ph(p, h4V_p(p));
s4v = s3_rhoT(1 / v4v, Ts);
v4L = v3_ph(p, h4L_p(p));
s4L = s3_rhoT(1 / v4L, Ts);
end
Out = fromSIunit_s((xs * s4v + (1 - xs) * s4L));
case 5
T = T5_ph(p, h);
Out = fromSIunit_s(s5_pT(p, T));
otherwise
Out = NaN;
end
%***********************************************************************************************************
%*1.8 Specific internal energy (u)
case 'uv_p'
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
if p < 16.529
Out = fromSIunit_u(u2_pT(p, T4_p(p)));
else
Out = fromSIunit_u(u3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)));
end
else
Out = NaN;
end
case 'ul_p'
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
if p < 16.529
Out = fromSIunit_u(u1_pT(p, T4_p(p)));
else
Out = fromSIunit_u(u3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)));
end
else
Out = NaN;
end
case 'uv_t'
T = toSIunit_T(In1);
if T > 273.15 && T < 647.096
if T <= 623.15
Out = fromSIunit_u(u2_pT(p4_T(T), T));
else
Out = fromSIunit_u(u3_rhoT(1 / (v3_ph(p4_T(T), h4V_p(p4_T(T)))), T));
end
else
Out = NaN;
end
case 'ul_t'
T = toSIunit_T(In1);
if T > 273.15 && T < 647.096
if T <= 623.15
Out = fromSIunit_u(u1_pT(p4_T(T), T));
else
Out = fromSIunit_u(u3_rhoT(1 / (v3_ph(p4_T(T), h4L_p(p4_T(T)))), T));
end
else
Out = NaN;
end
case 'u_pt'
p = toSIunit_p(In1);
T = toSIunit_T(In2);
Region = region_pT(p, T);
switch Region
case 1
Out = fromSIunit_u(u1_pT(p, T));
case 2
Out = fromSIunit_u(u2_pT(p, T));
case 3
hs = h3_pT(p, T);
rhos = 1 / v3_ph(p, hs);
Out = fromSIunit_u(u3_rhoT(rhos, T));
case 4
Out = NaN;
case 5
Out = fromSIunit_u(u5_pT(p, T));
otherwise
Out = NaN;
end
case 'u_ph'
p = toSIunit_p(In1);
h = toSIunit_h(In2);
Region = region_ph(p, h);
switch Region
case 1
Ts = T1_ph(p, h);
Out = fromSIunit_u(u1_pT(p, Ts));
case 2
Ts = T2_ph(p, h);
Out = fromSIunit_u(u2_pT(p, Ts));
case 3
rhos = 1 / v3_ph(p, h);
Ts = T3_ph(p, h);
Out = fromSIunit_u(u3_rhoT(rhos, Ts));
case 4
Ts = T4_p(p);
xs = x4_ph(p, h);
if p < 16.529
u4v = u2_pT(p, Ts);
u4L = u1_pT(p, Ts);
else
v4v = v3_ph(p, h4V_p(p));
u4v = u3_rhoT(1 / v4v, Ts);
v4L = v3_ph(p, h4L_p(p));
u4L = u3_rhoT(1 / v4L, Ts);
end
Out = fromSIunit_u((xs * u4v + (1 - xs) * u4L));
case 5
Ts = T5_ph(p, h);
Out = fromSIunit_u(u5_pT(p, Ts));
otherwise
Out = NaN;
end
case 'u_ps'
p = toSIunit_p(In1);
s = toSIunit_s(In2);
Region = region_ps(p, s);
switch Region
case 1
Ts = T1_ps(p, s);
Out = fromSIunit_u(u1_pT(p, Ts));
case 2
Ts = T2_ps(p, s);
Out = fromSIunit_u(u2_pT(p, Ts));
case 3
rhos = 1 / v3_ps(p, s);
Ts = T3_ps(p, s);
Out = fromSIunit_u(u3_rhoT(rhos, Ts));
case 4
if p < 16.529
uLp = u1_pT(p, T4_p(p));
uVp = u2_pT(p, T4_p(p));
else
uLp = u3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p));
uVp = u3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p));
end
xs = x4_ps(p, s);
Out = fromSIunit_u((xs * uVp + (1 - xs) * uLp));
case 5
Ts = T5_ps(p, s);
Out = fromSIunit_u(u5_pT(p, Ts));
otherwise
Out = NaN;
end
%***********************************************************************************************************
%*1.9 Specific isobaric heat capacity (Cp)
case 'cpv_p'
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
if p < 16.529
Out = fromSIunit_Cp(Cp2_pT(p, T4_p(p)));
else
Out = fromSIunit_Cp(Cp3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)));
end
else
Out = NaN;
end
case 'cpl_p'
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
if p < 16.529
Out = fromSIunit_Cp(Cp1_pT(p, T4_p(p)));
else
Out = fromSIunit_Cp(Cp3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)));
end
else
Out = NaN;
end
case 'cpv_t'
T = toSIunit_T(In1);
if T > 273.15 && T < 647.096
if T <= 623.15
Out = fromSIunit_Cp(Cp2_pT(p4_T(T), T));
else
Out = fromSIunit_Cp(Cp3_rhoT(1 / (v3_ph(p4_T(T), h4V_p(p4_T(T)))), T));
end
else
Out = NaN;
end
case 'cpl_t'
T = toSIunit_T(In1);
if T > 273.15 && T < 647.096
if T <= 623.15
Out = fromSIunit_Cp(Cp1_pT(p4_T(T), T));
else
Out = fromSIunit_Cp(Cp3_rhoT(1 / (v3_ph(p4_T(T), h4L_p(p4_T(T)))), T));
end
else
Out = NaN;
end
case 'cp_pt'
p = toSIunit_p(In1);
T = toSIunit_T(In2);
Region = region_pT(p, T);
switch Region
case 1
Out = fromSIunit_Cp(Cp1_pT(p, T));
case 2
Out = fromSIunit_Cp(Cp2_pT(p, T));
case 3
hs = h3_pT(p, T);
rhos = 1 / v3_ph(p, hs);
Out = fromSIunit_Cp(Cp3_rhoT(rhos, T));
case 4
Out = NaN;
case 5
Out = fromSIunit_Cp(Cp5_pT(p, T));
otherwise
Out = NaN;
end
case 'cp_ph'
p = toSIunit_p(In1);
h = toSIunit_h(In2);
Region = region_ph(p, h);
switch Region
case 1
Ts = T1_ph(p, h);
Out = fromSIunit_Cp(Cp1_pT(p, Ts));
case 2
Ts = T2_ph(p, h);
Out = fromSIunit_Cp(Cp2_pT(p, Ts));
case 3
rhos = 1 / v3_ph(p, h);
Ts = T3_ph(p, h);
Out = fromSIunit_Cp(Cp3_rhoT(rhos, Ts));
case 4
Out = NaN;
case 5
Ts = T5_ph(p, h);
Out = fromSIunit_Cp(Cp5_pT(p, Ts));
otherwise
Out = NaN;
end
case 'cp_ps'
p = toSIunit_p(In1);
s = toSIunit_s(In2);
Region = region_ps(p, s);
switch Region
case 1
Ts = T1_ps(p, s);
Out = fromSIunit_Cp(Cp1_pT(p, Ts));
case 2
Ts = T2_ps(p, s);
Out = fromSIunit_Cp(Cp2_pT(p, Ts));
case 3
rhos = 1 / v3_ps(p, s);
Ts = T3_ps(p, s);
Out = fromSIunit_Cp(Cp3_rhoT(rhos, Ts));
case 4
Out = NaN;
case 5
Ts = T5_ps(p, s);
Out = fromSIunit_Cp(Cp5_pT(p, Ts));
otherwise
Out = NaN;
end
%***********************************************************************************************************
%*1.10 Specific isochoric heat capacity (Cv)
case 'cvv_p'
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
if p < 16.529
Out = fromSIunit_Cv(Cv2_pT(p, T4_p(p)));
else
Out = fromSIunit_Cv(Cv3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)));
end
else
Out = NaN;
end
case 'cvl_p'
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
if p < 16.529
Out = fromSIunit_Cv(Cv1_pT(p, T4_p(p)));
else
Out = fromSIunit_Cv(Cv3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)));
end
else
Out = NaN;
end
case 'cvv_t'
T = toSIunit_T(In1);
if T > 273.15 && T < 647.096
if T <= 623.15
Out = fromSIunit_Cv(Cv2_pT(p4_T(T), T));
else
Out = fromSIunit_Cv(Cv3_rhoT(1 / (v3_ph(p4_T(T), h4V_p(p4_T(T)))), T));
end
else
Out = NaN;
end
case 'cvl_t'
T = toSIunit_T(In1);
if T > 273.15 && T < 647.096
if T <= 623.15
Out = fromSIunit_Cv(Cv1_pT(p4_T(T), T));
else
Out = fromSIunit_Cv(Cv3_rhoT(1 / (v3_ph(p4_T(T), h4L_p(p4_T(T)))), T));
end
else
Out = NaN;
end
case 'cv_pt'
p = toSIunit_p(In1);
T = toSIunit_T(In2);
Region = region_pT(p, T);
switch Region
case 1
Out = fromSIunit_Cv(Cv1_pT(p, T));
case 2
Out = fromSIunit_Cv(Cv2_pT(p, T));
case 3
hs = h3_pT(p, T);
rhos = 1 / v3_ph(p, hs);
Out = fromSIunit_Cv(Cv3_rhoT(rhos, T));
case 4
Out = NaN;
case 5
Out = fromSIunit_Cv(Cv5_pT(p, T));
otherwise
Out = NaN;
end
case 'cv_ph'
p = toSIunit_p(In1);
h = toSIunit_h(In2);
Region = region_ph(p, h);
switch Region
case 1
Ts = T1_ph(p, h);
Out = fromSIunit_Cv(Cv1_pT(p, Ts));
case 2
Ts = T2_ph(p, h);
Out = fromSIunit_Cv(Cv2_pT(p, Ts));
case 3
rhos = 1 / v3_ph(p, h);
Ts = T3_ph(p, h);
Out = fromSIunit_Cv(Cv3_rhoT(rhos, Ts));
case 4
Out = NaN;
case 5
Ts = T5_ph(p, h);
Out = fromSIunit_Cv(Cv5_pT(p, Ts));
otherwise
Out = NaN;
end
case 'cv_ps'
p = toSIunit_p(In1);
s = toSIunit_s(In2);
Region = region_ps(p, s);
switch Region
case 1
Ts = T1_ps(p, s);
Out = fromSIunit_Cv(Cv1_pT(p, Ts));
case 2
Ts = T2_ps(p, s);
Out = fromSIunit_Cv(Cv2_pT(p, Ts));
case 3
rhos = 1 / v3_ps(p, s);
Ts = T3_ps(p, s);
Out = fromSIunit_Cv(Cv3_rhoT(rhos, Ts));
case 4
Out = NaN; %(xs * CvVp + (1 - xs) * CvLp) / Cv_scale - Cv_offset
case 5
Ts = T5_ps(p, s);
Out = fromSIunit_Cv(Cv5_pT(p, Ts));
otherwise
Out = CVErr(xlErrValue);
end
%***********************************************************************************************************
%*1.11 Speed of sound
case 'wv_p'
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
if p < 16.529
Out = fromSIunit_w(w2_pT(p, T4_p(p)));
else
Out = fromSIunit_w(w3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p)));
end
else
Out = NaN;
end
case 'wl_p'
p = toSIunit_p(In1);
if p > 0.000611657 && p < 22.06395
if p < 16.529
Out = fromSIunit_w(w1_pT(p, T4_p(p)));
else
Out = fromSIunit_w(w3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p)));
end
else
Out = NaN;
end
case 'wv_t'
T = toSIunit_T(In1);
if T > 273.15 && T < 647.096
if T <= 623.15
Out = fromSIunit_w(w2_pT(p4_T(T), T));
else
Out = fromSIunit_w(w3_rhoT(1 / (v3_ph(p4_T(T), h4V_p(p4_T(T)))), T));
end
else
Out = NaN;
end
case 'wl_t'
T = toSIunit_T(In1);
if T > 273.15 && T < 647.096
if T <= 623.15
Out = fromSIunit_w(w1_pT(p4_T(T), T));
else
Out = fromSIunit_w(w3_rhoT(1 / (v3_ph(p4_T(T), h4L_p(p4_T(T)))), T));
end
else
Out = NaN;
end
case 'w_pt'
p = toSIunit_p(In1);
T = toSIunit_T(In2);
Region = region_pT(p, T);
switch Region
case 1
Out = fromSIunit_w(w1_pT(p, T));
case 2
Out = fromSIunit_w(w2_pT(p, T));
case 3
hs = h3_pT(p, T);
rhos = 1 / v3_ph(p, hs);
Out = fromSIunit_w(w3_rhoT(rhos, T));
case 4
Out = NaN;
case 5
Out = fromSIunit_w(w5_pT(p, T));
otherwise
Out = NaN;
end
case 'w_ph'
p = toSIunit_p(In1);
h = toSIunit_h(In2);
Region = region_ph(p, h);
switch Region
case 1
Ts = T1_ph(p, h);
Out = fromSIunit_w(w1_pT(p, Ts));
case 2
Ts = T2_ph(p, h);
Out = fromSIunit_w(w2_pT(p, Ts));
case 3
rhos = 1 / v3_ph(p, h);
Ts = T3_ph(p, h);
Out = fromSIunit_w(w3_rhoT(rhos, Ts));
case 4
Out = NaN;
case 5
Ts = T5_ph(p, h);
Out = fromSIunit_w(w5_pT(p, Ts));
otherwise
Out = NaN;
end
case 'w_ps'
p = toSIunit_p(In1);
s = toSIunit_s(In2);
Region = region_ps(p, s);
switch Region
case 1
Ts = T1_ps(p, s);
Out = fromSIunit_w(w1_pT(p, Ts));
case 2
Ts = T2_ps(p, s);
Out = fromSIunit_w(w2_pT(p, Ts));
case 3
rhos = 1 / v3_ps(p, s);
Ts = T3_ps(p, s);
Out = fromSIunit_w(w3_rhoT(rhos, Ts));
case 4
Out = NaN; %(xs * wVp + (1 - xs) * wLp) / w_scale - w_offset
case 5
Ts = T5_ps(p, s);
Out = fromSIunit_w(w5_pT(p, Ts));
otherwise
Out = NaN;
end
%***********************************************************************************************************
%*1.12 Viscosity
case 'my_pt'
p = toSIunit_p(In1);
T = toSIunit_T(In2);
Region = region_pT(p, T);
switch Region
case 4
Out = NaN;
case {1, 2, 3, 5}
Out = fromSIunit_my(my_AllRegions_pT(p, T));
otherwise
Out = NaN;
end
case {'my_ph'}
p = toSIunit_p(In1);
h = toSIunit_h(In2);
Region = region_ph(p, h);
switch Region
case {1, 2, 3, 5}
Out = fromSIunit_my(my_AllRegions_ph(p, h));
case {4}
Out = NaN;
otherwise
Out = NaN;
end
case 'my_ps'
h = XSteam('h_ps',In1, In2);
Out = XSteam('my_ph',In1, h);
%***********************************************************************************************************
%*1.13 Prandtl
case 'pr_pt'
Cp = toSIunit_Cp(XSteam('Cp_pT',In1, In2));
my = toSIunit_my(XSteam('my_pT',In1,In2));
tc = toSIunit_tc(XSteam('tc_pT',In1,In2));
Out = Cp * 1000 * my / tc;
case 'pr_ph'
Cp = toSIunit_Cp(XSteam('Cp_ph',In1, In2));
my = toSIunit_my(XSteam('my_ph',In1,In2));
tc = toSIunit_tc(XSteam('tc_ph',In1,In2));
Out = Cp * 1000 * my / tc;
%***********************************************************************************************************
%*1.14 Kappa
%***********************************************************************************************************
%***********************************************************************************************************
%*1.15 Surface tension
case 'st_t'
T = toSIunit_T(In1);
Out = fromSIunit_st(Surface_Tension_T(T));
case 'st_p'
T = XSteam('Tsat_p',In1);
T = toSIunit_T(T);
Out = fromSIunit_st(Surface_Tension_T(T));
%***********************************************************************************************************
%*1.16 Thermal conductivity
case 'tcl_p'
T = XSteam('Tsat_p',In1);
v = XSteam('vL_p',In1);
p = toSIunit_p(In1);
T = toSIunit_T(T);
v = toSIunit_v(v);
rho = 1 / v;
Out = fromSIunit_tc(tc_ptrho(p, T, rho));
case 'tcv_p'
ps = In1;
T = XSteam('Tsat_p',ps);
v = XSteam('vV_p',ps);
p = toSIunit_p(In1);
T = toSIunit_T(T);
v = toSIunit_v(v);
rho = 1 / v;
Out = fromSIunit_tc(tc_ptrho(p, T, rho));
case 'tcl_t'
Ts = In1;
p = XSteam('psat_T',Ts);
v = XSteam('vL_T',Ts);
p = toSIunit_p(p);
T = toSIunit_T(Ts);
v = toSIunit_v(v);
rho = 1 / v;
Out = fromSIunit_tc(tc_ptrho(p, T, rho));
case 'tcv_t'
Ts = In1;
p = XSteam('psat_T',Ts);
v = XSteam('vV_T',Ts);
p = toSIunit_p(p);
T = toSIunit_T(Ts);
v = toSIunit_v(v);
rho = 1 / v;
Out = fromSIunit_tc(tc_ptrho(p, T, rho));
case 'tc_pt'
Ts = In2;
ps = In1;
v = XSteam('v_pT',ps, Ts);
p = toSIunit_p(ps);
T = toSIunit_T(Ts);
v = toSIunit_v(v);
rho = 1 / v;
Out = fromSIunit_tc(tc_ptrho(p, T, rho));
case 'tc_ph'
hs = In2;
ps = In1;
v = XSteam('v_ph',ps, hs);
T = XSteam('T_ph',ps, hs);
p = toSIunit_p(ps);
T = toSIunit_T(T);
v = toSIunit_v(v);
rho = 1 / v;
Out = fromSIunit_tc(tc_ptrho(p, T, rho));
case 'tc_hs'
hs = In1;
p = XSteam('p_hs',hs, In2);
ps = p;
v = XSteam('v_ph',ps, hs);
T = XSteam('T_ph',ps, hs);
p = toSIunit_p(p);
T = toSIunit_T(T);
v = toSIunit_v(v);
rho = 1 / v;
Out = fromSIunit_tc(tc_ptrho(p, T, rho));
%***********************************************************************************************************
%*1.17 Vapour fraction
case 'x_ph'
p = toSIunit_p(In1);
h = toSIunit_h(In2);
if p > 0.000611657 && p < 22.06395
Out = fromSIunit_x(x4_ph(p, h));
else
Out = NaN;
end
case 'x_ps'
p = toSIunit_p(In1);
s = toSIunit_s(In2);
if p > 0.000611657 && p < 22.06395
Out = fromSIunit_x(x4_ps(p, s));
else
Out = NaN;
end
%***********************************************************************************************************
%*1.18 Vapour Volume Fraction
case 'vx_ph'
p = toSIunit_p(In1);
h = toSIunit_h(In2);
if p > 0.000611657 && p < 22.06395
if p < 16.529
vL = v1_pT(p, T4_p(p));
vV = v2_pT(p, T4_p(p));
else
vL = v3_ph(p, h4L_p(p));
vV = v3_ph(p, h4V_p(p));
end
xs = x4_ph(p, h);
Out = fromSIunit_vx((xs * vV / (xs * vV + (1 - xs) * vL)));
else
Out = NaN;
end
case 'vx_ps'
p = toSIunit_p(In1);
s = toSIunit_s(In2);
if p > 0.000611657 && p < 22.06395
if p < 16.529
vL = v1_pT(p, T4_p(p));
vV = v2_pT(p, T4_p(p));
else
vL = v3_ph(p, h4L_p(p));
vV = v3_ph(p, h4V_p(p));
end
xs = x4_ps(p, s);
Out = fromSIunit_vx((xs * vV / (xs * vV + (1 - xs) * vL)));
else
Out = NaN;
end
case 'check'
err=check();
otherwise
error(['Unknown calling function to XSteam, ',fun, ' See help XSteam for valid calling functions']);
end
%***********************************************************************************************************
%*2 IAPWS IF 97 Calling functions *
%***********************************************************************************************************
%
%***********************************************************************************************************
%*2.1 Functions for region 1
function v1_pT = v1_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation
%Eqution 7, Table 3, Page 6
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32];
J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41];
n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26];
R = 0.461526; %kJ/(kg K)
Pi = p / 16.53;
tau = 1386 / T;
gamma_der_pi = 0;
for i = 1 : 34
gamma_der_pi = gamma_der_pi - n1(i) * I1(i) * (7.1 - Pi) ^ (I1(i) - 1) * (tau - 1.222) ^ J1(i);
end
v1_pT = R * T / p * Pi * gamma_der_pi / 1000;
function h1_pT = h1_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation
%Eqution 7, Table 3, Page 6
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32];
J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41];
n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26];
R = 0.461526; %kJ/(kg K)
Pi = p / 16.53;
tau = 1386 / T;
gamma_der_tau = 0;
for i = 1 : 34
gamma_der_tau = gamma_der_tau + (n1(i) * (7.1 - Pi) ^ I1(i) * J1(i) * (tau - 1.222) ^ (J1(i) - 1));
end
h1_pT = R * T * tau * gamma_der_tau;
function u1_pT = u1_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation
%Eqution 7, Table 3, Page 6
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32];
J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41];
n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26];
R = 0.461526; %kJ/(kg K)
Pi = p / 16.53;
tau = 1386 / T;
gamma_der_tau = 0;
gamma_der_pi = 0;
for i = 1 : 34
gamma_der_pi = gamma_der_pi - n1(i) * I1(i) * (7.1 - Pi) ^ (I1(i) - 1) * (tau - 1.222) ^ J1(i);
gamma_der_tau = gamma_der_tau + (n1(i) * (7.1 - Pi) ^ I1(i) * J1(i) * (tau - 1.222) ^ (J1(i) - 1));
end
u1_pT = R * T * (tau * gamma_der_tau - Pi * gamma_der_pi);
function s1_pT = s1_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation
%Eqution 7, Table 3, Page 6
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32];
J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41];
n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26];
R = 0.461526; %kJ/(kg K)
Pi = p / 16.53;
tau = 1386 / T;
gamma = 0;
gamma_der_tau = 0;
for i = 1 : 34
gamma_der_tau = gamma_der_tau + (n1(i) * (7.1 - Pi) ^ I1(i) * J1(i) * (tau - 1.222) ^ (J1(i) - 1));
gamma = gamma + n1(i) * (7.1 - Pi) ^ I1(i) * (tau - 1.222) ^ J1(i);
end
s1_pT = R * tau * gamma_der_tau - R * gamma;
function Cp1_pT = Cp1_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation
%Eqution 7, Table 3, Page 6
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32];
J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41];
n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26];
R = 0.461526; %kJ/(kg K)
Pi = p / 16.53;
tau = 1386 / T;
gamma_der_tautau = 0;
for i = 1 :34
gamma_der_tautau = gamma_der_tautau + (n1(i) * (7.1 - Pi) ^ I1(i) * J1(i) * (J1(i) - 1) * (tau - 1.222) ^ (J1(i) - 2));
end
Cp1_pT = -R * tau ^ 2 * gamma_der_tautau;
function Cv1_pT = Cv1_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation
%Eqution 7, Table 3, Page 6
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32];
J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41];
n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26];
R = 0.461526; %kJ/(kg K)
Pi = p / 16.53;
tau = 1386 / T;
gamma_der_pi = 0;
gamma_der_pipi = 0;
gamma_der_pitau = 0;
gamma_der_tautau = 0;
for i = 1 : 34
gamma_der_pi = gamma_der_pi - n1(i) * I1(i) * (7.1 - Pi) ^ (I1(i) - 1) * (tau - 1.222) ^ J1(i);
gamma_der_pipi = gamma_der_pipi + n1(i) * I1(i) * (I1(i) - 1) * (7.1 - Pi) ^ (I1(i) - 2) * (tau - 1.222) ^ J1(i);
gamma_der_pitau = gamma_der_pitau - n1(i) * I1(i) * (7.1 - Pi) ^ (I1(i) - 1) * J1(i) * (tau - 1.222) ^ (J1(i) - 1);
gamma_der_tautau = gamma_der_tautau + n1(i) * (7.1 - Pi) ^ I1(i) * J1(i) * (J1(i) - 1) * (tau - 1.222) ^ (J1(i) - 2);
end
Cv1_pT = R * (-tau ^ 2 * gamma_der_tautau + (gamma_der_pi - tau * gamma_der_pitau) ^ 2 / gamma_der_pipi);
function w1_pT = w1_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation
%Eqution 7, Table 3, Page 6
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32];
J1 = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41];
n1 = [0.14632971213167, -0.84548187169114, -3.756360367204, 3.3855169168385, -0.95791963387872, 0.15772038513228, -0.016616417199501, 8.1214629983568E-04, 2.8319080123804E-04, -6.0706301565874E-04, -0.018990068218419, -0.032529748770505, -0.021841717175414, -5.283835796993E-05, -4.7184321073267E-04, -3.0001780793026E-04, 4.7661393906987E-05, -4.4141845330846E-06, -7.2694996297594E-16, -3.1679644845054E-05, -2.8270797985312E-06, -8.5205128120103E-10, -2.2425281908E-06, -6.5171222895601E-07, -1.4341729937924E-13, -4.0516996860117E-07, -1.2734301741641E-09, -1.7424871230634E-10, -6.8762131295531E-19, 1.4478307828521E-20, 2.6335781662795E-23, -1.1947622640071E-23, 1.8228094581404E-24, -9.3537087292458E-26];
R = 0.461526; %kJ/(kg K)
Pi = p / 16.53;
tau = 1386 / T;
gamma_der_pi = 0;
gamma_der_pipi = 0;
gamma_der_pitau = 0;
gamma_der_tautau = 0;
for i = 1 : 34
gamma_der_pi = gamma_der_pi - n1(i) * I1(i) * (7.1 - Pi) ^ (I1(i) - 1) * (tau - 1.222) ^ J1(i);
gamma_der_pipi = gamma_der_pipi + n1(i) * I1(i) * (I1(i) - 1) * (7.1 - Pi) ^ (I1(i) - 2) * (tau - 1.222) ^ J1(i);
gamma_der_pitau = gamma_der_pitau - n1(i) * I1(i) * (7.1 - Pi) ^ (I1(i) - 1) * J1(i) * (tau - 1.222) ^ (J1(i) - 1);
gamma_der_tautau = gamma_der_tautau + n1(i) * (7.1 - Pi) ^ I1(i) * J1(i) * (J1(i) - 1) * (tau - 1.222) ^ (J1(i) - 2);
end
w1_pT = (1000 * R * T * gamma_der_pi ^ 2 / ((gamma_der_pi - tau * gamma_der_pitau) ^ 2 / (tau ^ 2 * gamma_der_tautau) - gamma_der_pipi)) ^ 0.5;
function T1_ph = T1_ph(p, h)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation, 5.2.1 The Backward Equation T ( p,h )
%Eqution 11, Table 6, Page 10
I1 = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 4, 5, 6];
J1 = [0, 1, 2, 6, 22, 32, 0, 1, 2, 3, 4, 10, 32, 10, 32, 10, 32, 32, 32, 32];
n1 = [-238.72489924521, 404.21188637945, 113.49746881718, -5.8457616048039, -1.528548241314E-04, -1.0866707695377E-06, -13.391744872602, 43.211039183559, -54.010067170506, 30.535892203916, -6.5964749423638, 9.3965400878363E-03, 1.157364750534E-07, -2.5858641282073E-05, -4.0644363084799E-09, 6.6456186191635E-08, 8.0670734103027E-11, -9.3477771213947E-13, 5.8265442020601E-15, -1.5020185953503E-17];
Pi = p / 1;
eta = h / 2500;
T = 0;
for i = 1 : 20
T = T + n1(i) * Pi ^ I1(i) * (eta + 1) ^ J1(i);
end
T1_ph = T;
function T1_ps = T1_ps(p, s)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%5 Equations for Region 1, Section. 5.1 Basic Equation, 5.2.2 The Backward Equation T ( p, s )
%Eqution 13, Table 8, Page 11
I1 = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 4];
J1 = [0, 1, 2, 3, 11, 31, 0, 1, 2, 3, 12, 31, 0, 1, 2, 9, 31, 10, 32, 32];
n1 = [174.78268058307, 34.806930892873, 6.5292584978455, 0.33039981775489, -1.9281382923196E-07, -2.4909197244573E-23, -0.26107636489332, 0.22592965981586, -0.064256463395226, 7.8876289270526E-03, 3.5672110607366E-10, 1.7332496994895E-24, 5.6608900654837E-04, -3.2635483139717E-04, 4.4778286690632E-05, -5.1322156908507E-10, -4.2522657042207E-26, 2.6400441360689E-13, 7.8124600459723E-29, -3.0732199903668E-31];
Pi = p / 1;
Sigma = s / 1;
T = 0;
for i = 1 : 20
T = T + n1(i) * Pi ^ I1(i) * (Sigma + 2) ^ J1(i);
end
T1_ps = T;
function p1_hs = p1_hs(h, s)
%Supplementary Release on Backward Equations for Pressure as a Function of Enthalpy and Entropy p(h,s) to the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%5 Backward Equation p(h,s) for Region 1
%Eqution 1, Table 2, Page 5
I1 = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 5];
J1 = [0, 1, 2, 4, 5, 6, 8, 14, 0, 1, 4, 6, 0, 1, 10, 4, 1, 4, 0];
n1 = [-0.691997014660582, -18.361254878756, -9.28332409297335, 65.9639569909906, -16.2060388912024, 450.620017338667, 854.68067822417, 6075.23214001162, 32.6487682621856, -26.9408844582931, -319.9478483343, -928.35430704332, 30.3634537455249, -65.0540422444146, -4309.9131651613, -747.512324096068, 730.000345529245, 1142.84032569021, -436.407041874559];
eta = h / 3400;
Sigma = s / 7.6;
p = 0;
for i = 1 : 19
p = p + n1(i) * (eta + 0.05) ^ I1(i) * (Sigma + 0.05) ^ J1(i);
end
p1_hs = p * 100;
function T1_prho = T1_prho(p ,rho)
%Solve by iteration. Observe that for low temperatures this equation has 2 solutions.
%Solve with half interval method
Low_Bound = 273.15;
High_Bound = T4_p(p);
rhos=-1000;
while abs(rho - rhos) > 0.00001
Ts = (Low_Bound + High_Bound) / 2;
rhos = 1 / v1_pT(p, Ts);
if rhos < rho
High_Bound = Ts;
else
Low_Bound = Ts;
end
end
T1_prho = Ts;
%***********************************************************************************************************
%*2.2 functions for region 2
function v2_pT = v2_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2, Section. 6.1 Basic Equation
%Table 11 and 12, Page 14 and 15
J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3];
n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307];
Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24];
Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58];
nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07];
R = 0.461526; %kJ/(kg K)
Pi = p;
tau = 540 / T;
g0_pi = 1 / Pi;
gr_pi = 0;
for i = 1 : 43
gr_pi = gr_pi + nr(i) * Ir(i) * Pi ^ (Ir(i) - 1) * (tau - 0.5) ^ Jr(i);
end
v2_pT = R * T / p * Pi * (g0_pi + gr_pi) / 1000;
function h2_pT = h2_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2, Section. 6.1 Basic Equation
%Table 11 and 12, Page 14 and 15
J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3];
n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307];
Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24];
Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58];
nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07];
R = 0.461526; %kJ/(kg K)
Pi = p;
tau = 540 / T;
g0_tau = 0;
for i = 1 : 9
g0_tau = g0_tau + n0(i) * J0(i) * tau ^ (J0(i) - 1);
end
gr_tau = 0;
for i = 1 : 43
gr_tau = gr_tau + nr(i) * Pi ^ Ir(i) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1);
end
h2_pT = R * T * tau * (g0_tau + gr_tau);
function u2_pT = u2_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2, Section. 6.1 Basic Equation
%Table 11 and 12, Page 14 and 15
J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3];
n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307];
Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24];
Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58];
nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07];
R = 0.461526; %kJ/(kg K)
Pi = p;
tau = 540 / T;
g0_pi = 1 / Pi;
g0_tau = 0;
for i = 1 : 9
g0_tau = g0_tau + n0(i) * J0(i) * tau ^ (J0(i) - 1);
end
gr_pi = 0;
gr_tau = 0;
for i = 1 : 43
gr_pi = gr_pi + nr(i) * Ir(i) * Pi ^ (Ir(i) - 1) * (tau - 0.5) ^ Jr(i);
gr_tau = gr_tau + nr(i) * Pi ^ Ir(i) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1);
end
u2_pT = R * T * (tau * (g0_tau + gr_tau) - Pi * (g0_pi + gr_pi));
function s2_pT = s2_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2, Section. 6.1 Basic Equation
%Table 11 and 12, Page 14 and 15
J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3];
n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307];
Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24];
Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58];
nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07];
R = 0.461526; %kJ/(kg K)
Pi = p;
tau = 540 / T;
g0 = log(Pi);
g0_tau = 0;
for i = 1 : 9
g0 = g0 + n0(i) * tau ^ J0(i);
g0_tau = g0_tau + n0(i) * J0(i) * tau ^ (J0(i) - 1);
end
gr = 0;
gr_tau = 0;
for i = 1 : 43
gr = gr + nr(i) * Pi ^ Ir(i) * (tau - 0.5) ^ Jr(i);
gr_tau = gr_tau + nr(i) * Pi ^ Ir(i) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1);
end
s2_pT = R * (tau * (g0_tau + gr_tau) - (g0 + gr));
function Cp2_pT = Cp2_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2, Section. 6.1 Basic Equation
%Table 11 and 12, Page 14 and 15
J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3];
n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307];
Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24];
Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58];
nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07];
R = 0.461526; %kJ/(kg K)
Pi = p;
tau = 540 / T;
g0_tautau = 0;
for i = 1 : 9
g0_tautau = g0_tautau + n0(i) * J0(i) * (J0(i) - 1) * tau ^ (J0(i) - 2);
end
gr_tautau = 0;
for i = 1 : 43
gr_tautau = gr_tautau + nr(i) * Pi ^ Ir(i) * Jr(i) * (Jr(i) - 1) * (tau - 0.5) ^ (Jr(i) - 2);
end
Cp2_pT = -R * tau ^ 2 * (g0_tautau + gr_tautau);
function Cv2_pT = Cv2_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2, Section. 6.1 Basic Equation
%Table 11 and 12, Page 14 and 15
J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3];
n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307];
Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24];
Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58];
nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07];
R = 0.461526; %kJ/(kg K)
Pi = p;
tau = 540 / T;
g0_tautau = 0;
for i = 1 : 9
g0_tautau = g0_tautau + n0(i) * J0(i) * (J0(i) - 1) * tau ^ (J0(i) - 2);
end
gr_pi = 0;
gr_pitau = 0;
gr_pipi = 0;
gr_tautau = 0;
for i = 1 : 43
gr_pi = gr_pi + nr(i) * Ir(i) * Pi ^ (Ir(i) - 1) * (tau - 0.5) ^ Jr(i);
gr_pipi = gr_pipi + nr(i) * Ir(i) * (Ir(i) - 1) * Pi ^ (Ir(i) - 2) * (tau - 0.5) ^ Jr(i);
gr_pitau = gr_pitau + nr(i) * Ir(i) * Pi ^ (Ir(i) - 1) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1);
gr_tautau = gr_tautau + nr(i) * Pi ^ Ir(i) * Jr(i) * (Jr(i) - 1) * (tau - 0.5) ^ (Jr(i) - 2);
end
Cv2_pT = R * (-tau ^ 2 * (g0_tautau + gr_tautau) - (1 + Pi * gr_pi - tau * Pi * gr_pitau) ^ 2 / (1 - Pi ^ 2 * gr_pipi));
function w2_pT = w2_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2, Section. 6.1 Basic Equation
%Table 11 and 12, Page 14 and 15
J0 = [0, 1, -5, -4, -3, -2, -1, 2, 3];
n0 = [-9.6927686500217, 10.086655968018, -0.005608791128302, 0.071452738081455, -0.40710498223928, 1.4240819171444, -4.383951131945, -0.28408632460772, 0.021268463753307];
Ir = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24];
Jr = [0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58];
nr = [-1.7731742473213E-03, -0.017834862292358, -0.045996013696365, -0.057581259083432, -0.05032527872793, -3.3032641670203E-05, -1.8948987516315E-04, -3.9392777243355E-03, -0.043797295650573, -2.6674547914087E-05, 2.0481737692309E-08, 4.3870667284435E-07, -3.227767723857E-05, -1.5033924542148E-03, -0.040668253562649, -7.8847309559367E-10, 1.2790717852285E-08, 4.8225372718507E-07, 2.2922076337661E-06, -1.6714766451061E-11, -2.1171472321355E-03, -23.895741934104, -5.905956432427E-18, -1.2621808899101E-06, -0.038946842435739, 1.1256211360459E-11, -8.2311340897998, 1.9809712802088E-08, 1.0406965210174E-19, -1.0234747095929E-13, -1.0018179379511E-09, -8.0882908646985E-11, 0.10693031879409, -0.33662250574171, 8.9185845355421E-25, 3.0629316876232E-13, -4.2002467698208E-06, -5.9056029685639E-26, 3.7826947613457E-06, -1.2768608934681E-15, 7.3087610595061E-29, 5.5414715350778E-17, -9.436970724121E-07];
R = 0.461526; %kJ/(kg K)
Pi = p;
tau = 540 / T;
g0_tautau = 0;
for i = 1 : 9
g0_tautau = g0_tautau + n0(i) * J0(i) * (J0(i) - 1) * tau ^ (J0(i) - 2);
end
gr_pi = 0;
gr_pitau = 0;
gr_pipi = 0;
gr_tautau = 0;
for i = 1 : 43
gr_pi = gr_pi + nr(i) * Ir(i) * Pi ^ (Ir(i) - 1) * (tau - 0.5) ^ Jr(i);
gr_pipi = gr_pipi + nr(i) * Ir(i) * (Ir(i) - 1) * Pi ^ (Ir(i) - 2) * (tau - 0.5) ^ Jr(i);
gr_pitau = gr_pitau + nr(i) * Ir(i) * Pi ^ (Ir(i) - 1) * Jr(i) * (tau - 0.5) ^ (Jr(i) - 1);
gr_tautau = gr_tautau + nr(i) * Pi ^ Ir(i) * Jr(i) * (Jr(i) - 1) * (tau - 0.5) ^ (Jr(i) - 2);
end
w2_pT = (1000 * R * T * (1 + 2 * Pi * gr_pi + Pi ^ 2 * gr_pi ^ 2) / ((1 - Pi ^ 2 * gr_pipi) + (1 + Pi * gr_pi - tau * Pi * gr_pitau) ^ 2 / (tau ^ 2 * (g0_tautau + gr_tautau)))) ^ 0.5;
function T2_ph = T2_ph(p, h)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2,6.3.1 The Backward Equations T( p, h ) for Subregions 2a, 2b, and 2c
if p < 4
sub_reg = 1;
else
if p < (905.84278514723 - 0.67955786399241 * h + 1.2809002730136E-04 * h ^ 2)
sub_reg = 2;
else
sub_reg = 3;
end
end
switch sub_reg
case 1
%Subregion A
%Table 20, Eq 22, page 22
Ji = [0, 1, 2, 3, 7, 20, 0, 1, 2, 3, 7, 9, 11, 18, 44, 0, 2, 7, 36, 38, 40, 42, 44, 24, 44, 12, 32, 44, 32, 36, 42, 34, 44, 28];
Ii = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 7];
ni = [1089.8952318288, 849.51654495535, -107.81748091826, 33.153654801263, -7.4232016790248, 11.765048724356, 1.844574935579, -4.1792700549624, 6.2478196935812, -17.344563108114, -200.58176862096, 271.96065473796, -455.11318285818, 3091.9688604755, 252266.40357872, -6.1707422868339E-03, -0.31078046629583, 11.670873077107, 128127984.04046, -985549096.23276, 2822454697.3002, -3594897141.0703, 1722734991.3197, -13551.334240775, 12848734.66465, 1.3865724283226, 235988.32556514, -13105236.545054, 7399.9835474766, -551966.9703006, 3715408.5996233, 19127.72923966, -415351.64835634, -62.459855192507];
Ts = 0;
hs = h / 2000;
for i = 1 : 34
Ts = Ts + ni(i) * p ^ (Ii(i)) * (hs - 2.1) ^ Ji(i);
end
T2_ph = Ts;
case 2
%Subregion B
%Table 21, Eq 23, page 23
Ji = [0, 1, 2, 12, 18, 24, 28, 40, 0, 2, 6, 12, 18, 24, 28, 40, 2, 8, 18, 40, 1, 2, 12, 24, 2, 12, 18, 24, 28, 40, 18, 24, 40, 28, 2, 28, 1, 40];
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 6, 7, 7, 9, 9];
ni = [1489.5041079516, 743.07798314034, -97.708318797837, 2.4742464705674, -0.63281320016026, 1.1385952129658, -0.47811863648625, 8.5208123431544E-03, 0.93747147377932, 3.3593118604916, 3.3809355601454, 0.16844539671904, 0.73875745236695, -0.47128737436186, 0.15020273139707, -0.002176411421975, -0.021810755324761, -0.10829784403677, -0.046333324635812, 7.1280351959551E-05, 1.1032831789999E-04, 1.8955248387902E-04, 3.0891541160537E-03, 1.3555504554949E-03, 2.8640237477456E-07, -1.0779857357512E-05, -7.6462712454814E-05, 1.4052392818316E-05, -3.1083814331434E-05, -1.0302738212103E-06, 2.821728163504E-07, 1.2704902271945E-06, 7.3803353468292E-08, -1.1030139238909E-08, -8.1456365207833E-14, -2.5180545682962E-11, -1.7565233969407E-18, 8.6934156344163E-15];
Ts = 0;
hs = h / 2000;
for i = 1 : 38
Ts = Ts + ni(i) * (p - 2) ^ (Ii(i)) * (hs - 2.6) ^ Ji(i);
end
T2_ph = Ts;
otherwise
%Subregion C
%Table 22, Eq 24, page 24
Ji = [0, 4, 0, 2, 0, 2, 0, 1, 0, 2, 0, 1, 4, 8, 4, 0, 1, 4, 10, 12, 16, 20, 22];
Ii = [-7, -7, -6, -6, -5, -5, -2, -2, -1, -1, 0, 0, 1, 1, 2, 6, 6, 6, 6, 6, 6, 6, 6];
ni = [-3236839855524.2, 7326335090218.1, 358250899454.47, -583401318515.9, -10783068217.47, 20825544563.171, 610747.83564516, 859777.2253558, -25745.72360417, 31081.088422714, 1208.2315865936, 482.19755109255, 3.7966001272486, -10.842984880077, -0.04536417267666, 1.4559115658698E-13, 1.126159740723E-12, -1.7804982240686E-11, 1.2324579690832E-07, -1.1606921130984E-06, 2.7846367088554E-05, -5.9270038474176E-04, 1.2918582991878E-03];
Ts = 0;
hs = h / 2000;
for i = 1 : 23
Ts = Ts + ni(i) * (p + 25) ^ (Ii(i)) * (hs - 1.8) ^ Ji(i);
end
T2_ph = Ts;
end
function T2_ps = T2_ps(p, s)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%6 Equations for Region 2,6.3.2 The Backward Equations T( p, s ) for Subregions 2a, 2b, and 2c
%Page 26
if p < 4
sub_reg = 1;
else
if s < 5.85
sub_reg = 3;
else
sub_reg = 2;
end
end
switch sub_reg
case 1
%Subregion A
%Table 25, Eq 25, page 26
Ii = [-1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.25, -1.25, -1.25, -1, -1, -1, -1, -1, -1, -0.75, -0.75, -0.5, -0.5, -0.5, -0.5, -0.25, -0.25, -0.25, -0.25, 0.25, 0.25, 0.25, 0.25, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.75, 0.75, 0.75, 0.75, 1, 1, 1.25, 1.25, 1.5, 1.5];
Ji = [-24, -23, -19, -13, -11, -10, -19, -15, -6, -26, -21, -17, -16, -9, -8, -15, -14, -26, -13, -9, -7, -27, -25, -11, -6, 1, 4, 8, 11, 0, 1, 5, 6, 10, 14, 16, 0, 4, 9, 17, 7, 18, 3, 15, 5, 18];
ni = [-392359.83861984, 515265.7382727, 40482.443161048, -321.93790923902, 96.961424218694, -22.867846371773, -449429.14124357, -5011.8336020166, 0.35684463560015, 44235.33584819, -13673.388811708, 421632.60207864, 22516.925837475, 474.42144865646, -149.31130797647, -197811.26320452, -23554.39947076, -19070.616302076, 55375.669883164, 3829.3691437363, -603.91860580567, 1936.3102620331, 4266.064369861, -5978.0638872718, -704.01463926862, 338.36784107553, 20.862786635187, 0.033834172656196, -4.3124428414893E-05, 166.53791356412, -139.86292055898, -0.78849547999872, 0.072132411753872, -5.9754839398283E-03, -1.2141358953904E-05, 2.3227096733871E-07, -10.538463566194, 2.0718925496502, -0.072193155260427, 2.074988708112E-07, -0.018340657911379, 2.9036272348696E-07, 0.21037527893619, 2.5681239729999E-04, -0.012799002933781, -8.2198102652018E-06];
Pi = p;
Sigma = s / 2;
teta = 0;
for i = 1 : 46
teta = teta + ni(i) * Pi ^ Ii(i) * (Sigma - 2) ^ Ji(i);
end
T2_ps = teta;
case 2
%Subregion B
%Table 26, Eq 26, page 27
Ii = [-6, -6, -5, -5, -4, -4, -4, -3, -3, -3, -3, -2, -2, -2, -2, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 5];
Ji = [0, 11, 0, 11, 0, 1, 11, 0, 1, 11, 12, 0, 1, 6, 10, 0, 1, 5, 8, 9, 0, 1, 2, 4, 5, 6, 9, 0, 1, 2, 3, 7, 8, 0, 1, 5, 0, 1, 3, 0, 1, 0, 1, 2];
ni = [316876.65083497, 20.864175881858, -398593.99803599, -21.816058518877, 223697.85194242, -2784.1703445817, 9.920743607148, -75197.512299157, 2970.8605951158, -3.4406878548526, 0.38815564249115, 17511.29508575, -1423.7112854449, 1.0943803364167, 0.89971619308495, -3375.9740098958, 471.62885818355, -1.9188241993679, 0.41078580492196, -0.33465378172097, 1387.0034777505, -406.63326195838, 41.72734715961, 2.1932549434532, -1.0320050009077, 0.35882943516703, 5.2511453726066E-03, 12.838916450705, -2.8642437219381, 0.56912683664855, -0.099962954584931, -3.2632037778459E-03, 2.3320922576723E-04, -0.1533480985745, 0.029072288239902, 3.7534702741167E-04, 1.7296691702411E-03, -3.8556050844504E-04, -3.5017712292608E-05, -1.4566393631492E-05, 5.6420857267269E-06, 4.1286150074605E-08, -2.0684671118824E-08, 1.6409393674725E-09];
Pi = p;
Sigma = s / 0.7853;
teta = 0;
for i = 1 : 44
teta = teta + ni(i) * Pi ^ Ii(i) * (10 - Sigma) ^ Ji(i);
end
T2_ps = teta;
otherwise
%Subregion C
%Table 27, Eq 27, page 28
Ii = [-2, -2, -1, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 7, 7, 7, 7, 7];
Ji = [0, 1, 0, 0, 1, 2, 3, 0, 1, 3, 4, 0, 1, 2, 0, 1, 5, 0, 1, 4, 0, 1, 2, 0, 1, 0, 1, 3, 4, 5];
ni = [909.68501005365, 2404.566708842, -591.6232638713, 541.45404128074, -270.98308411192, 979.76525097926, -469.66772959435, 14.399274604723, -19.104204230429, 5.3299167111971, -21.252975375934, -0.3114733441376, 0.60334840894623, -0.042764839702509, 5.8185597255259E-03, -0.014597008284753, 5.6631175631027E-03, -7.6155864584577E-05, 2.2440342919332E-04, -1.2561095013413E-05, 6.3323132660934E-07, -2.0541989675375E-06, 3.6405370390082E-08, -2.9759897789215E-09, 1.0136618529763E-08, 5.9925719692351E-12, -2.0677870105164E-11, -2.0874278181886E-11, 1.0162166825089E-10, -1.6429828281347E-10];
Pi = p;
Sigma = s / 2.9251;
teta = 0;
for i = 1 : 30
teta = teta + ni(i) * Pi ^ Ii(i) * (2 - Sigma) ^ Ji(i);
end
T2_ps = teta;
end
function p2_hs = p2_hs(h, s)
%Supplementary Release on Backward Equations for Pressure as a function of Enthalpy and Entropy p(h,s) to the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%Chapter 6:Backward Equations p(h,s) for Region 2
if h < (-3498.98083432139 + 2575.60716905876 * s - 421.073558227969 * s ^ 2 + 27.6349063799944 * s ^ 3)
sub_reg = 1;
else
if s < 5.85
sub_reg = 3;
else
sub_reg = 2;
end
end
switch sub_reg
case 1
%Subregion A
%Table 6, Eq 3, page 8
Ii = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 4, 5, 5, 6, 7];
Ji = [1, 3, 6, 16, 20, 22, 0, 1, 2, 3, 5, 6, 10, 16, 20, 22, 3, 16, 20, 0, 2, 3, 6, 16, 16, 3, 16, 3, 1];
ni = [-1.82575361923032E-02, -0.125229548799536, 0.592290437320145, 6.04769706185122, 238.624965444474, -298.639090222922, 0.051225081304075, -0.437266515606486, 0.413336902999504, -5.16468254574773, -5.57014838445711, 12.8555037824478, 11.414410895329, -119.504225652714, -2847.7798596156, 4317.57846408006, 1.1289404080265, 1974.09186206319, 1516.12444706087, 1.41324451421235E-02, 0.585501282219601, -2.97258075863012, 5.94567314847319, -6236.56565798905, 9659.86235133332, 6.81500934948134, -6332.07286824489, -5.5891922446576, 4.00645798472063E-02];
eta = h / 4200;
Sigma = s / 12;
Pi = 0;
for i = 1 : 29
Pi = Pi + ni(i) * (eta - 0.5) ^ Ii(i) * (Sigma - 1.2) ^ Ji(i);
end
p2_hs = Pi ^ 4 * 4;
case 2
%Subregion B
%Table 7, Eq 4, page 9
Ii = [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 8, 8, 12, 14];
Ji = [0, 1, 2, 4, 8, 0, 1, 2, 3, 5, 12, 1, 6, 18, 0, 1, 7, 12, 1, 16, 1, 12, 1, 8, 18, 1, 16, 1, 3, 14, 18, 10, 16];
ni = [8.01496989929495E-02, -0.543862807146111, 0.337455597421283, 8.9055545115745, 313.840736431485, 0.797367065977789, -1.2161697355624, 8.72803386937477, -16.9769781757602, -186.552827328416, 95115.9274344237, -18.9168510120494, -4334.0703719484, 543212633.012715, 0.144793408386013, 128.024559637516, -67230.9534071268, 33697238.0095287, -586.63419676272, -22140322476.9889, 1716.06668708389, -570817595.806302, -3121.09693178482, -2078413.8463301, 3056059461577.86, 3221.57004314333, 326810259797.295, -1441.04158934487, 410.694867802691, 109077066873.024, -24796465425889.3, 1888019068.65134, -123651009018773];
eta = h / 4100;
Sigma = s / 7.9;
Pi = 0;
for i = 1 : 33
Pi = Pi + ni(i) * (eta - 0.6) ^ Ii(i) * (Sigma - 1.01) ^ Ji(i);
end
p2_hs = Pi ^ 4 * 100;
otherwise
%Subregion C
%Table 8, Eq 5, page 10
Ii = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 5, 5, 5, 5, 6, 6, 10, 12, 16];
Ji = [0, 1, 2, 3, 4, 8, 0, 2, 5, 8, 14, 2, 3, 7, 10, 18, 0, 5, 8, 16, 18, 18, 1, 4, 6, 14, 8, 18, 7, 7, 10];
ni = [0.112225607199012, -3.39005953606712, -32.0503911730094, -197.5973051049, -407.693861553446, 13294.3775222331, 1.70846839774007, 37.3694198142245, 3581.44365815434, 423014.446424664, -751071025.760063, 52.3446127607898, -228.351290812417, -960652.417056937, -80705929.2526074, 1626980172256.69, 0.772465073604171, 46392.9973837746, -13731788.5134128, 1704703926305.12, -25110462818730.8, 31774883083552, 53.8685623675312, -55308.9094625169, -1028615.22421405, 2042494187562.34, 273918446.626977, -2.63963146312685E+15, -1078908541.08088, -29649262098.0124, -1.11754907323424E+15];
eta = h / 3500;
Sigma = s / 5.9;
Pi = 0;
for i = 1 : 31
Pi = Pi + ni(i) * (eta - 0.7) ^ Ii(i) * (Sigma - 1.1) ^ Ji(i);
end
p2_hs = Pi ^ 4 * 100;
end
function T2_prho=T2_prho(p,rho)
%Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
%Solve with half interval method
if p < 16.5292
Low_Bound = T4_p(p);
else
Low_Bound = B23T_p(p);
end
High_Bound = 1073.15;
rhos=-1000;
while abs(rho - rhos) > 0.000001
Ts = (Low_Bound + High_Bound) / 2;
rhos = 1 / v2_pT(p, Ts);
if rhos < rho
High_Bound = Ts;
else
Low_Bound = Ts;
end
end
T2_prho = Ts;
%***********************************************************************************************************
%*2.3 functions for region 3
function p3_rhoT = p3_rhoT(rho, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
%Table 30 and 31, Page 30 and 31
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11];
Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26];
ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05];
R = 0.461526; %kJ/(KgK)
tc = 647.096; %K
pc = 22.064; %MPa
rhoc = 322; %kg/m3
delta = rho / rhoc;
tau = tc / T;
fidelta = 0;
for i = 2 : 40
fidelta = fidelta + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * tau ^ Ji(i);
end
fidelta = fidelta + ni(1) / delta;
p3_rhoT = rho * R * T * delta * fidelta / 1000;
function u3_rhoT = u3_rhoT(rho, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
%Table 30 and 31, Page 30 and 31
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11];
Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26];
ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05];
R = 0.461526; %kJ/(KgK)
tc = 647.096; %K
pc = 22.064; %MPa
rhoc = 322; %kg/m3
delta = rho / rhoc;
tau = tc / T;
fitau = 0;
for i = 2 : 40
fitau = fitau + ni(i) * delta ^ Ii(i) * Ji(i) * tau ^ (Ji(i) - 1);
end
u3_rhoT = R * T * (tau * fitau);
function h3_rhoT = h3_rhoT(rho, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
%Table 30 and 31, Page 30 and 31
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11];
Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26];
ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05];
R = 0.461526; %kJ/(KgK)
tc = 647.096; %K
pc = 22.064; %MPa
rhoc = 322; %kg/m3
delta = rho / rhoc;
tau = tc / T;
fidelta = 0;
fitau = 0;
for i = 2 : 40
fidelta = fidelta + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * tau ^ Ji(i);
fitau = fitau + ni(i) * delta ^ Ii(i) * Ji(i) * tau ^ (Ji(i) - 1);
end
fidelta = fidelta + ni(1) / delta;
h3_rhoT = R * T * (tau * fitau + delta * fidelta);
function s3_rhoT = s3_rhoT(rho, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
%Table 30 and 31, Page 30 and 31
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11];
Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26];
ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05];
R = 0.461526; %kJ/(KgK)
tc = 647.096; %K
pc = 22.064; %MPa
rhoc = 322; %kg/m3
delta = rho / rhoc;
tau = tc / T;
fi = 0;
fitau = 0;
for i = 2 : 40
fi = fi + ni(i) * delta ^ Ii(i) * tau ^ Ji(i);
fitau = fitau + ni(i) * delta ^ Ii(i) * Ji(i) * tau ^ (Ji(i) - 1);
end
fi = fi + ni(1) * log(delta);
s3_rhoT = R * (tau * fitau - fi);
function Cp3_rhoT = Cp3_rhoT(rho, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
%Table 30 and 31, Page 30 and 31
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11];
Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26];
ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05];
R = 0.461526; %kJ/(KgK)
tc = 647.096; %K
pc = 22.064; %MPa
rhoc = 322; %kg/m3
delta = rho / rhoc;
tau = tc / T;
fitautau = 0;
fidelta = 0;
fideltatau = 0;
fideltadelta = 0;
for i = 2 : 40
fitautau = fitautau + ni(i) * delta ^ Ii(i) * Ji(i) * (Ji(i) - 1) * tau ^ (Ji(i) - 2);
fidelta = fidelta + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * tau ^ Ji(i);
fideltatau = fideltatau + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * Ji(i) * tau ^ (Ji(i) - 1);
fideltadelta = fideltadelta + ni(i) * Ii(i) * (Ii(i) - 1) * delta ^ (Ii(i) - 2) * tau ^ Ji(i);
end
fidelta = fidelta + ni(1) / delta;
fideltadelta = fideltadelta - ni(1) / (delta ^ 2);
Cp3_rhoT = R * (-tau ^ 2 * fitautau + (delta * fidelta - delta * tau * fideltatau) ^ 2 / (2 * delta * fidelta + delta ^ 2 * fideltadelta));
function Cv3_rhoT = Cv3_rhoT(rho, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
%Table 30 and 31, Page 30 and 31
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11];
Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26];
ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05];
R = 0.461526; %kJ/(KgK)
tc = 647.096; %K
pc = 22.064; %MPa
rhoc = 322; %kg/m3
delta = rho / rhoc;
tau = tc / T;
fitautau = 0;
for i = 1 : 40
fitautau = fitautau + ni(i) * delta ^ Ii(i) * Ji(i) * (Ji(i) - 1) * tau ^ (Ji(i) - 2);
end
Cv3_rhoT = R * -(tau * tau * fitautau);
function w3_rhoT = w3_rhoT(rho, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%7 Basic Equation for Region 3, Section. 6.1 Basic Equation
%Table 30 and 31, Page 30 and 31
Ii = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 8, 9, 9, 10, 10, 11];
Ji = [0, 0, 1, 2, 7, 10, 12, 23, 2, 6, 15, 17, 0, 2, 6, 7, 22, 26, 0, 2, 4, 16, 26, 0, 2, 4, 26, 1, 3, 26, 0, 2, 26, 2, 26, 2, 26, 0, 1, 26];
ni = [1.0658070028513, -15.732845290239, 20.944396974307, -7.6867707878716, 2.6185947787954, -2.808078114862, 1.2053369696517, -8.4566812812502E-03, -1.2654315477714, -1.1524407806681, 0.88521043984318, -0.64207765181607, 0.38493460186671, -0.85214708824206, 4.8972281541877, -3.0502617256965, 0.039420536879154, 0.12558408424308, -0.2799932969871, 1.389979956946, -2.018991502357, -8.2147637173963E-03, -0.47596035734923, 0.0439840744735, -0.44476435428739, 0.90572070719733, 0.70522450087967, 0.10770512626332, -0.32913623258954, -0.50871062041158, -0.022175400873096, 0.094260751665092, 0.16436278447961, -0.013503372241348, -0.014834345352472, 5.7922953628084E-04, 3.2308904703711E-03, 8.0964802996215E-05, -1.6557679795037E-04, -4.4923899061815E-05];
R = 0.461526; %kJ/(KgK)
tc = 647.096; %K
pc = 22.064; %MPa
rhoc = 322; %kg/m3
delta = rho / rhoc;
tau = tc / T;
fitautau = 0;
fidelta = 0;
fideltatau = 0;
fideltadelta = 0;
for i = 2 : 40
fitautau = fitautau + ni(i) * delta ^ Ii(i) * Ji(i) * (Ji(i) - 1) * tau ^ (Ji(i) - 2);
fidelta = fidelta + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * tau ^ Ji(i);
fideltatau = fideltatau + ni(i) * Ii(i) * delta ^ (Ii(i) - 1) * Ji(i) * tau ^ (Ji(i) - 1);
fideltadelta = fideltadelta + ni(i) * Ii(i) * (Ii(i) - 1) * delta ^ (Ii(i) - 2) * tau ^ Ji(i);
end
fidelta = fidelta + ni(1) / delta;
fideltadelta = fideltadelta - ni(1) / (delta ^ 2);
w3_rhoT = (1000 * R * T * (2 * delta * fidelta + delta ^ 2 * fideltadelta - (delta * fidelta - delta * tau * fideltatau) ^ 2 / (tau ^ 2 * fitautau))) ^ 0.5;
function T3_ph = T3_ph(p, h)
%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%2004
%Section 3.3 Backward Equations T(p,h) and v(p,h) for Subregions 3a and 3b
%Boundary equation, Eq 1 Page 5
h3ab = 2014.64004206875 + 3.74696550136983 * p - 2.19921901054187E-02 * p ^ 2 + 8.7513168600995E-05 * p ^ 3;
if h < h3ab
%Subregion 3a
%Eq 2, Table 3, Page 7
Ii = [-12, -12, -12, -12, -12, -12, -12, -12, -10, -10, -10, -8, -8, -8, -8, -5, -3, -2, -2, -2, -1, -1, 0, 0, 1, 3, 3, 4, 4, 10, 12];
Ji = [0, 1, 2, 6, 14, 16, 20, 22, 1, 5, 12, 0, 2, 4, 10, 2, 0, 1, 3, 4, 0, 2, 0, 1, 1, 0, 1, 0, 3, 4, 5];
ni = [-1.33645667811215E-07, 4.55912656802978E-06, -1.46294640700979E-05, 6.3934131297008E-03, 372.783927268847, -7186.54377460447, 573494.7521034, -2675693.29111439, -3.34066283302614E-05, -2.45479214069597E-02, 47.8087847764996, 7.64664131818904E-06, 1.28350627676972E-03, 1.71219081377331E-02, -8.51007304583213, -1.36513461629781E-02, -3.84460997596657E-06, 3.37423807911655E-03, -0.551624873066791, 0.72920227710747, -9.92522757376041E-03, -0.119308831407288, 0.793929190615421, 0.454270731799386, 0.20999859125991, -6.42109823904738E-03, -0.023515586860454, 2.52233108341612E-03, -7.64885133368119E-03, 1.36176427574291E-02, -1.33027883575669E-02];
ps = p / 100;
hs = h / 2300;
Ts = 0;
for i = 1 : 31
Ts = Ts + ni(i) * (ps + 0.24) ^ Ii(i) * (hs - 0.615) ^ Ji(i);
end
T3_ph = Ts * 760;
else
%Subregion 3b
%Eq 3, Table 4, Page 7,8
Ii = [-12, -12, -10, -10, -10, -10, -10, -8, -8, -8, -8, -8, -6, -6, -6, -4, -4, -3, -2, -2, -1, -1, -1, -1, -1, -1, 0, 0, 1, 3, 5, 6, 8];
Ji = [0, 1, 0, 1, 5, 10, 12, 0, 1, 2, 4, 10, 0, 1, 2, 0, 1, 5, 0, 4, 2, 4, 6, 10, 14, 16, 0, 2, 1, 1, 1, 1, 1];
ni = [3.2325457364492E-05, -1.27575556587181E-04, -4.75851877356068E-04, 1.56183014181602E-03, 0.105724860113781, -85.8514221132534, 724.140095480911, 2.96475810273257E-03, -5.92721983365988E-03, -1.26305422818666E-02, -0.115716196364853, 84.9000969739595, -1.08602260086615E-02, 1.54304475328851E-02, 7.50455441524466E-02, 2.52520973612982E-02, -6.02507901232996E-02, -3.07622221350501, -5.74011959864879E-02, 5.03471360939849, -0.925081888584834, 3.91733882917546, -77.314600713019, 9493.08762098587, -1410437.19679409, 8491662.30819026, 0.861095729446704, 0.32334644281172, 0.873281936020439, -0.436653048526683, 0.286596714529479, -0.131778331276228, 6.76682064330275E-03];
hs = h / 2800;
ps = p / 100;
Ts = 0;
for i = 1 : 33
Ts = Ts + ni(i) * (ps + 0.298) ^ Ii(i) * (hs - 0.72) ^ Ji(i);
end
T3_ph = Ts * 860;
end
function v3_ph = v3_ph(p, h)
%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%2004
%Section 3.3 Backward Equations T(p,h) and v(p,h) for Subregions 3a and 3b
%Boundary equation, Eq 1 Page 5
h3ab = 2014.64004206875 + 3.74696550136983 * p - 2.19921901054187E-02 * p ^ 2 + 8.7513168600995E-05 * p ^ 3;
if h < h3ab
%Subregion 3a
%Eq 4, Table 6, Page 9
Ii = [-12, -12, -12, -12, -10, -10, -10, -8, -8, -6, -6, -6, -4, -4, -3, -2, -2, -1, -1, -1, -1, 0, 0, 1, 1, 1, 2, 2, 3, 4, 5, 8];
Ji = [6, 8, 12, 18, 4, 7, 10, 5, 12, 3, 4, 22, 2, 3, 7, 3, 16, 0, 1, 2, 3, 0, 1, 0, 1, 2, 0, 2, 0, 2, 2, 2];
ni = [5.29944062966028E-03, -0.170099690234461, 11.1323814312927, -2178.98123145125, -5.06061827980875E-04, 0.556495239685324, -9.43672726094016, -0.297856807561527, 93.9353943717186, 1.92944939465981E-02, 0.421740664704763, -3689141.2628233, -7.37566847600639E-03, -0.354753242424366, -1.99768169338727, 1.15456297059049, 5683.6687581596, 8.08169540124668E-03, 0.172416341519307, 1.04270175292927, -0.297691372792847, 0.560394465163593, 0.275234661176914, -0.148347894866012, -6.51142513478515E-02, -2.92468715386302, 6.64876096952665E-02, 3.52335014263844, -1.46340792313332E-02, -2.24503486668184, 1.10533464706142, -4.08757344495612E-02];
ps = p / 100;
hs = h / 2100;
vs = 0;
for i = 1 : 32
vs = vs + ni(i) * (ps + 0.128) ^ Ii(i) * (hs - 0.727) ^ Ji(i);
end
v3_ph = vs * 0.0028;
else
%Subregion 3b
%Eq 5, Table 7, Page 9
Ii = [-12, -12, -8, -8, -8, -8, -8, -8, -6, -6, -6, -6, -6, -6, -4, -4, -4, -3, -3, -2, -2, -1, -1, -1, -1, 0, 1, 1, 2, 2];
Ji = [0, 1, 0, 1, 3, 6, 7, 8, 0, 1, 2, 5, 6, 10, 3, 6, 10, 0, 2, 1, 2, 0, 1, 4, 5, 0, 0, 1, 2, 6];
ni = [-2.25196934336318E-09, 1.40674363313486E-08, 2.3378408528056E-06, -3.31833715229001E-05, 1.07956778514318E-03, -0.271382067378863, 1.07202262490333, -0.853821329075382, -2.15214194340526E-05, 7.6965608822273E-04, -4.31136580433864E-03, 0.453342167309331, -0.507749535873652, -100.475154528389, -0.219201924648793, -3.21087965668917, 607.567815637771, 5.57686450685932E-04, 0.18749904002955, 9.05368030448107E-03, 0.285417173048685, 3.29924030996098E-02, 0.239897419685483, 4.82754995951394, -11.8035753702231, 0.169490044091791, -1.79967222507787E-02, 3.71810116332674E-02, -5.36288335065096E-02, 1.6069710109252];
ps = p / 100;
hs = h / 2800;
vs = 0;
for i = 1 : 30
vs = vs + ni(i) * (ps + 0.0661) ^ Ii(i) * (hs - 0.72) ^ Ji(i);
end
v3_ph = vs * 0.0088;
end
function T3_ps = T3_ps(p, s)
%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%2004
%3.4 Backward Equations T(p,s) and v(p,s) for Subregions 3a and 3b
%Boundary equation, Eq 6 Page 11
if s <= 4.41202148223476
%Subregion 3a
%Eq 6, Table 10, Page 11
Ii = [-12, -12, -10, -10, -10, -10, -8, -8, -8, -8, -6, -6, -6, -5, -5, -5, -4, -4, -4, -2, -2, -1, -1, 0, 0, 0, 1, 2, 2, 3, 8, 8, 10];
Ji = [28, 32, 4, 10, 12, 14, 5, 7, 8, 28, 2, 6, 32, 0, 14, 32, 6, 10, 36, 1, 4, 1, 6, 0, 1, 4, 0, 0, 3, 2, 0, 1, 2];
ni = [1500420082.63875, -159397258480.424, 5.02181140217975E-04, -67.2057767855466, 1450.58545404456, -8238.8953488889, -0.154852214233853, 11.2305046746695, -29.7000213482822, 43856513263.5495, 1.37837838635464E-03, -2.97478527157462, 9717779473494.13, -5.71527767052398E-05, 28830.794977842, -74442828926270.3, 12.8017324848921, -368.275545889071, 6.64768904779177E+15, 0.044935925195888, -4.22897836099655, -0.240614376434179, -4.74341365254924, 0.72409399912611, 0.923874349695897, 3.99043655281015, 3.84066651868009E-02, -3.59344365571848E-03, -0.735196448821653, 0.188367048396131, 1.41064266818704E-04, -2.57418501496337E-03, 1.23220024851555E-03];
Sigma = s / 4.4;
Pi = p / 100;
teta = 0;
for i = 1 : 33
teta = teta + ni(i) * (Pi + 0.24) ^ Ii(i) * (Sigma - 0.703) ^ Ji(i);
end
T3_ps = teta * 760;
else
%Subregion 3b
%Eq 7, Table 11, Page 11
Ii = [-12, -12, -12, -12, -8, -8, -8, -6, -6, -6, -5, -5, -5, -5, -5, -4, -3, -3, -2, 0, 2, 3, 4, 5, 6, 8, 12, 14];
Ji = [1, 3, 4, 7, 0, 1, 3, 0, 2, 4, 0, 1, 2, 4, 6, 12, 1, 6, 2, 0, 1, 1, 0, 24, 0, 3, 1, 2];
ni = [0.52711170160166, -40.1317830052742, 153.020073134484, -2247.99398218827, -0.193993484669048, -1.40467557893768, 42.6799878114024, 0.752810643416743, 22.6657238616417, -622.873556909932, -0.660823667935396, 0.841267087271658, -25.3717501764397, 485.708963532948, 880.531517490555, 2650155.92794626, -0.359287150025783, -656.991567673753, 2.41768149185367, 0.856873461222588, 0.655143675313458, -0.213535213206406, 5.62974957606348E-03, -316955725450471, -6.99997000152457E-04, 1.19845803210767E-02, 1.93848122022095E-05, -2.15095749182309E-05];
Sigma = s / 5.3;
Pi = p / 100;
teta = 0;
for i = 1 : 28
teta = teta + ni(i) * (Pi + 0.76) ^ Ii(i) * (Sigma - 0.818) ^ Ji(i);
end
T3_ps = teta * 860;
end
function v3_ps = v3_ps(p, s)
%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h) and T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%2004
%3.4 Backward Equations T(p,s) and v(p,s) for Subregions 3a and 3b
%Boundary equation, Eq 6 Page 11
if s <= 4.41202148223476
%Subregion 3a
%Eq 8, Table 13, Page 14
Ii = [-12, -12, -12, -10, -10, -10, -10, -8, -8, -8, -8, -6, -5, -4, -3, -3, -2, -2, -1, -1, 0, 0, 0, 1, 2, 4, 5, 6];
Ji = [10, 12, 14, 4, 8, 10, 20, 5, 6, 14, 16, 28, 1, 5, 2, 4, 3, 8, 1, 2, 0, 1, 3, 0, 0, 2, 2, 0];
ni = [79.5544074093975, -2382.6124298459, 17681.3100617787, -1.10524727080379E-03, -15.3213833655326, 297.544599376982, -35031520.6871242, 0.277513761062119, -0.523964271036888, -148011.182995403, 1600148.99374266, 1708023226634.27, 2.46866996006494E-04, 1.6532608479798, -0.118008384666987, 2.537986423559, 0.965127704669424, -28.2172420532826, 0.203224612353823, 1.10648186063513, 0.52612794845128, 0.277000018736321, 1.08153340501132, -7.44127885357893E-02, 1.64094443541384E-02, -6.80468275301065E-02, 0.025798857610164, -1.45749861944416E-04];
Pi = p / 100;
Sigma = s / 4.4;
omega = 0;
for i = 1 : 28
omega = omega + ni(i) * (Pi + 0.187) ^ Ii(i) * (Sigma - 0.755) ^ Ji(i);
end
v3_ps = omega * 0.0028;
else
%Subregion 3b
%Eq 9, Table 14, Page 14
Ii = [-12, -12, -12, -12, -12, -12, -10, -10, -10, -10, -8, -5, -5, -5, -4, -4, -4, -4, -3, -2, -2, -2, -2, -2, -2, 0, 0, 0, 1, 1, 2];
Ji = [0, 1, 2, 3, 5, 6, 0, 1, 2, 4, 0, 1, 2, 3, 0, 1, 2, 3, 1, 0, 1, 2, 3, 4, 12, 0, 1, 2, 0, 2, 2];
ni = [5.91599780322238E-05, -1.85465997137856E-03, 1.04190510480013E-02, 5.9864730203859E-03, -0.771391189901699, 1.72549765557036, -4.67076079846526E-04, 1.34533823384439E-02, -8.08094336805495E-02, 0.508139374365767, 1.28584643361683E-03, -1.63899353915435, 5.86938199318063, -2.92466667918613, -6.14076301499537E-03, 5.76199014049172, -12.1613320606788, 1.67637540957944, -7.44135838773463, 3.78168091437659E-02, 4.01432203027688, 16.0279837479185, 3.17848779347728, -3.58362310304853, -1159952.60446827, 0.199256573577909, -0.122270624794624, -19.1449143716586, -1.50448002905284E-02, 14.6407900162154, -3.2747778718823];
Pi = p / 100;
Sigma = s / 5.3;
omega = 0;
for i = 1 : 31
omega = omega + ni(i) * (Pi + 0.298) ^ Ii(i) * (Sigma - 0.816) ^ Ji(i);
end
v3_ps = omega * 0.0088;
end
function p3_hs = p3_hs(h, s)
%Supplementary Release on Backward Equations ( ) , p h s for Region 3,
%Equations as a function of h and s for the Region Boundaries, and an Equation
%( ) sat , T hs for Region 4 of the IAPWS Industrial formulation 1997 for the
%Thermodynamic Properties of Water and Steam
%2004
%Section 3 Backward functions p(h,s), T(h,s), and v(h,s) for Region 3
if s < 4.41202148223476
%Subregion 3a
%Eq 1, Table 3, Page 8
Ii = [0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 6, 7, 8, 10, 10, 14, 18, 20, 22, 22, 24, 28, 28, 32, 32];
Ji = [0, 1, 5, 0, 3, 4, 8, 14, 6, 16, 0, 2, 3, 0, 1, 4, 5, 28, 28, 24, 1, 32, 36, 22, 28, 36, 16, 28, 36, 16, 36, 10, 28];
ni = [7.70889828326934, -26.0835009128688, 267.416218930389, 17.2221089496844, -293.54233214597, 614.135601882478, -61056.2757725674, -65127225.1118219, 73591.9313521937, -11664650591.4191, 35.5267086434461, -596.144543825955, -475.842430145708, 69.6781965359503, 335.674250377312, 25052.6809130882, 146997.380630766, 5.38069315091534E+19, 1.43619827291346E+21, 3.64985866165994E+19, -2547.41561156775, 2.40120197096563E+27, -3.93847464679496E+29, 1.47073407024852E+24, -4.26391250432059E+31, 1.94509340621077E+38, 6.66212132114896E+23, 7.06777016552858E+33, 1.75563621975576E+41, 1.08408607429124E+28, 7.30872705175151E+43, 1.5914584739887E+24, 3.77121605943324E+40];
Sigma = s / 4.4;
eta = h / 2300;
Pi = 0;
for i = 1 : 33
Pi = Pi + ni(i) * (eta - 1.01) ^ Ii(i) * (Sigma - 0.75) ^ Ji(i);
end
p3_hs = Pi * 99;
else
%Subregion 3b
%Eq 2, Table 4, Page 8
Ii = [-12, -12, -12, -12, -12, -10, -10, -10, -10, -8, -8, -6, -6, -6, -6, -5, -4, -4, -4, -3, -3, -3, -3, -2, -2, -1, 0, 2, 2, 5, 6, 8, 10, 14, 14];
Ji = [2, 10, 12, 14, 20, 2, 10, 14, 18, 2, 8, 2, 6, 7, 8, 10, 4, 5, 8, 1, 3, 5, 6, 0, 1, 0, 3, 0, 1, 0, 1, 1, 1, 3, 7];
ni = [1.25244360717979E-13, -1.26599322553713E-02, 5.06878030140626, 31.7847171154202, -391041.161399932, -9.75733406392044E-11, -18.6312419488279, 510.973543414101, 373847.005822362, 2.99804024666572E-08, 20.0544393820342, -4.98030487662829E-06, -10.230180636003, 55.2819126990325, -206.211367510878, -7940.12232324823, 7.82248472028153, -58.6544326902468, 3550.73647696481, -1.15303107290162E-04, -1.75092403171802, 257.98168774816, -727.048374179467, 1.21644822609198E-04, 3.93137871762692E-02, 7.04181005909296E-03, -82.910820069811, -0.26517881813125, 13.7531682453991, -52.2394090753046, 2405.56298941048, -22736.1631268929, 89074.6343932567, -23923456.5822486, 5687958081.29714];
Sigma = s / 5.3;
eta = h / 2800;
Pi = 0;
for i = 1 : 35
Pi = Pi + ni(i) * (eta - 0.681) ^ Ii(i) * (Sigma - 0.792) ^ Ji(i);
end
p3_hs = 16.6 / Pi;
end
function h3_pT = h3_pT(p, T)
%Not avalible with if 97
%Solve function T3_ph-T=0 with half interval method.
%ver2.6 Start corrected bug
if p < 22.06395 %Bellow tripple point
Ts = T4_p(p); %Saturation temperature
if T <= Ts %Liquid side
High_Bound = h4L_p(p); %Max h är liauid h.
Low_Bound = h1_pT(p, 623.15);
else
Low_Bound = h4V_p(p); %Min h är Vapour h.
High_Bound = h2_pT(p, B23T_p(p));
end
else %Above tripple point. R3 from R2 till R3.
Low_Bound = h1_pT(p, 623.15);
High_Bound = h2_pT(p, B23T_p(p));
end
%ver2.6 End corrected bug
Ts = T+1;
while abs(T - Ts) > 0.00001
hs = (Low_Bound + High_Bound) / 2;
Ts = T3_ph(p, hs);
if Ts > T
High_Bound = hs;
else
Low_Bound = hs;
end
end
h3_pT = hs;
function T3_prho = T3_prho(p, rho)
%Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
%Solve with half interval method
Low_Bound = 623.15;
High_Bound = 1073.15;
ps=-1000;
while abs(p - ps) > 0.00000001
Ts = (Low_Bound + High_Bound) / 2;
ps = p3_rhoT(rho, Ts);
if ps > p
High_Bound = Ts;
else
Low_Bound = Ts;
end
end
T3_prho = Ts;
%***********************************************************************************************************
%*2.4 functions for region 4
function p4_T = p4_T(T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Section 8.1 The Saturation-Pressure Equation
%Eq 30, Page 33
teta = T - 0.23855557567849 / (T - 650.17534844798);
a = teta ^ 2 + 1167.0521452767 * teta - 724213.16703206;
B = -17.073846940092 * teta ^ 2 + 12020.82470247 * teta - 3232555.0322333;
C = 14.91510861353 * teta ^ 2 - 4823.2657361591 * teta + 405113.40542057;
p4_T = (2 * C / (-B + (B ^ 2 - 4 * a * C) ^ 0.5)) ^ 4;
function T4_p = T4_p(p)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Section 8.2 The Saturation-Temperature Equation
%Eq 31, Page 34
beta = p ^ 0.25;
E = beta ^ 2 - 17.073846940092 * beta + 14.91510861353;
f = 1167.0521452767 * beta ^ 2 + 12020.82470247 * beta - 4823.2657361591;
G = -724213.16703206 * beta ^ 2 - 3232555.0322333 * beta + 405113.40542057;
D = 2 * G / (-f - (f ^ 2 - 4 * E * G) ^ 0.5);
T4_p = (650.17534844798 + D - ((650.17534844798 + D) ^ 2 - 4 * (-0.23855557567849 + 650.17534844798 * D)) ^ 0.5) / 2;
function h4_s = h4_s(s)
%Supplementary Release on Backward Equations ( ) , p h s for Region 3,Equations as a function of h and s for the Region Boundaries, and an Equation( ) sat , T hs for Region 4 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%4 Equations for Region Boundaries Given Enthalpy and Entropy
% Se picture page 14
if (s > -0.0001545495919 & s <= 3.77828134)==1
%hL1_s
%Eq 3,Table 9,Page 16
Ii = [0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 7, 8, 12, 12, 14, 14, 16, 20, 20, 22, 24, 28, 32, 32];
Ji = [14, 36, 3, 16, 0, 5, 4, 36, 4, 16, 24, 18, 24, 1, 4, 2, 4, 1, 22, 10, 12, 28, 8, 3, 0, 6, 8];
ni = [0.332171191705237, 6.11217706323496E-04, -8.82092478906822, -0.45562819254325, -2.63483840850452E-05, -22.3949661148062, -4.28398660164013, -0.616679338856916, -14.682303110404, 284.523138727299, -113.398503195444, 1156.71380760859, 395.551267359325, -1.54891257229285, 19.4486637751291, -3.57915139457043, -3.35369414148819, -0.66442679633246, 32332.1885383934, 3317.66744667084, -22350.1257931087, 5739538.75852936, 173.226193407919, -3.63968822121321E-02, 8.34596332878346E-07, 5.03611916682674, 65.5444787064505];
Sigma = s / 3.8;
eta = 0;
for i = 1 : 27
eta = eta + ni(i) * (Sigma - 1.09) ^ Ii(i) * (Sigma + 0.0000366) ^ Ji(i);
end
h4_s = eta * 1700;
elseif (s > 3.77828134 & s <= 4.41202148223476 )==1
%hL3_s
%Eq 4,Table 10,Page 16
Ii = [0, 0, 0, 0, 2, 3, 4, 4, 5, 5, 6, 7, 7, 7, 10, 10, 10, 32, 32];
Ji = [1, 4, 10, 16, 1, 36, 3, 16, 20, 36, 4, 2, 28, 32, 14, 32, 36, 0, 6];
ni = [0.822673364673336, 0.181977213534479, -0.011200026031362, -7.46778287048033E-04, -0.179046263257381, 4.24220110836657E-02, -0.341355823438768, -2.09881740853565, -8.22477343323596, -4.99684082076008, 0.191413958471069, 5.81062241093136E-02, -1655.05498701029, 1588.70443421201, -85.0623535172818, -31771.4386511207, -94589.0406632871, -1.3927384708869E-06, 0.63105253224098];
Sigma = s / 3.8;
eta = 0;
for i = 1 : 19
eta = eta + ni(i) * (Sigma - 1.09) ^ Ii(i) * (Sigma + 0.0000366) ^ Ji(i);
end
h4_s = eta * 1700;
elseif (s > 4.41202148223476 & s <= 5.85 )==1
%Section 4.4 Equations ( ) 2ab " h s and ( ) 2c3b "h s for the Saturated Vapor Line
%Page 19, Eq 5
%hV2c3b_s(s)
Ii = [0, 0, 0, 1, 1, 5, 6, 7, 8, 8, 12, 16, 22, 22, 24, 36];
Ji = [0, 3, 4, 0, 12, 36, 12, 16, 2, 20, 32, 36, 2, 32, 7, 20];
ni = [1.04351280732769, -2.27807912708513, 1.80535256723202, 0.420440834792042, -105721.24483466, 4.36911607493884E+24, -328032702839.753, -6.7868676080427E+15, 7439.57464645363, -3.56896445355761E+19, 1.67590585186801E+31, -3.55028625419105E+37, 396611982166.538, -4.14716268484468E+40, 3.59080103867382E+18, -1.16994334851995E+40];
Sigma = s / 5.9;
eta = 0;
for i = 1 : 16
eta = eta + ni(i) * (Sigma - 1.02) ^ Ii(i) * (Sigma - 0.726) ^ Ji(i);
end
h4_s = eta ^ 4 * 2800;
elseif (s > 5.85 & s < 9.155759395)==1
%Section 4.4 Equations ( ) 2ab " h s and ( ) 2c3b "h s for the Saturated Vapor Line
%Page 20, Eq 6
Ii = [1, 1, 2, 2, 4, 4, 7, 8, 8, 10, 12, 12, 18, 20, 24, 28, 28, 28, 28, 28, 32, 32, 32, 32, 32, 36, 36, 36, 36, 36];
Ji = [8, 24, 4, 32, 1, 2, 7, 5, 12, 1, 0, 7, 10, 12, 32, 8, 12, 20, 22, 24, 2, 7, 12, 14, 24, 10, 12, 20, 22, 28];
ni = [-524.581170928788, -9269472.18142218, -237.385107491666, 21077015581.2776, -23.9494562010986, 221.802480294197, -5104725.33393438, 1249813.96109147, 2000084369.96201, -815.158509791035, -157.612685637523, -11420042233.2791, 6.62364680776872E+15, -2.27622818296144E+18, -1.71048081348406E+31, 6.60788766938091E+15, 1.66320055886021E+22, -2.18003784381501E+29, -7.87276140295618E+29, 1.51062329700346E+31, 7957321.70300541, 1.31957647355347E+15, -3.2509706829914E+23, -4.18600611419248E+25, 2.97478906557467E+34, -9.53588761745473E+19, 1.66957699620939E+24, -1.75407764869978E+32, 3.47581490626396E+34, -7.10971318427851E+38];
Sigma1 = s / 5.21;
Sigma2 = s / 9.2;
eta = 0;
for i = 1 : 30
eta = eta + ni(i) * (1 / Sigma1 - 0.513) ^ Ii(i) * (Sigma2 - 0.524) ^ Ji(i);
end
h4_s = exp(eta) * 2800;
else
h4_s = -99999;
end
function p4_s = p4_s(s)
%Uses h4_s and p_hs for the diffrent regions to determine p4_s
h_sat = h4_s(s);
if (s > -0.0001545495919 & s <= 3.77828134)==1
p4_s = p1_hs(h_sat, s);
elseif (s > 3.77828134 & s <= 5.210887663)==1
p4_s = p3_hs(h_sat, s);
elseif (s > 5.210887663 & s < 9.155759395)==1
p4_s = p2_hs(h_sat, s);
else
p4_s = -99999;
end
function h4L_p = h4L_p(p)
if (p > 0.000611657 & p < 22.06395)==1
Ts = T4_p(p);
if p < 16.529
h4L_p = h1_pT(p, Ts);
else
%Iterate to find the the backward solution of p3sat_h
Low_Bound = 1670.858218;
High_Bound = 2087.23500164864;
ps=-1000;
while abs(p - ps) > 0.00001
hs = (Low_Bound + High_Bound) / 2;
ps = p3sat_h(hs);
if ps > p
High_Bound = hs;
else
Low_Bound = hs;
end
end
h4L_p = hs;
end
else
h4L_p = -99999;
end
function h4V_p = h4V_p(p)
if (p > 0.000611657 & p < 22.06395)==1
Ts = T4_p(p);
if p < 16.529
h4V_p = h2_pT(p, Ts);
else
%Iterate to find the the backward solution of p3sat_h
Low_Bound = 2087.23500164864;
High_Bound = 2563.592004+5;
ps=-1000;
while abs(p - ps) > 0.000001
hs = (Low_Bound + High_Bound) / 2;
ps = p3sat_h(hs);
if ps < p
High_Bound = hs;
else
Low_Bound = hs;
end
end
h4V_p = hs;
end
else
h4V_p = -99999;
end
function x4_ph = x4_ph(p, h)
%Calculate vapour fraction from hL and hV for given p
h4v = h4V_p(p);
h4L = h4L_p(p);
if h > h4v
x4_ph = 1;
elseif h < h4L
x4_ph = 0;
else
x4_ph = (h - h4L) / (h4v - h4L);
end
function x4_ps = x4_ps(p, s)
if p < 16.529
ssv = s2_pT(p, T4_p(p));
ssL = s1_pT(p, T4_p(p));
else
ssv = s3_rhoT(1 / (v3_ph(p, h4V_p(p))), T4_p(p));
ssL = s3_rhoT(1 / (v3_ph(p, h4L_p(p))), T4_p(p));
end
if s < ssL
x4_ps = 0;
elseif s > ssv
x4_ps = 1;
else
x4_ps = (s - ssL) / (ssv - ssL);
end
function T4_hs = T4_hs(h, s)
%Supplementary Release on Backward Equations ( ) , p h s for Region 3,
%Chapter 5.3 page 30.
%The if 97 function is only valid for part of region4. Use iteration outsida.
Ii = [0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 5, 5, 5, 5, 6, 6, 6, 8, 10, 10, 12, 14, 14, 16, 16, 18, 18, 18, 20, 28];
Ji = [0, 3, 12, 0, 1, 2, 5, 0, 5, 8, 0, 2, 3, 4, 0, 1, 1, 2, 4, 16, 6, 8, 22, 1, 20, 36, 24, 1, 28, 12, 32, 14, 22, 36, 24, 36];
ni = [0.179882673606601, -0.267507455199603, 1.162767226126, 0.147545428713616, -0.512871635973248, 0.421333567697984, 0.56374952218987, 0.429274443819153, -3.3570455214214, 10.8890916499278, -0.248483390456012, 0.30415322190639, -0.494819763939905, 1.07551674933261, 7.33888415457688E-02, 1.40170545411085E-02, -0.106110975998808, 1.68324361811875E-02, 1.25028363714877, 1013.16840309509, -1.51791558000712, 52.4277865990866, 23049.5545563912, 2.49459806365456E-02, 2107964.67412137, 366836848.613065, -144814105.365163, -1.7927637300359E-03, 4899556021.00459, 471.262212070518, -82929439019.8652, -1715.45662263191, 3557776.82973575, 586062760258.436, -12988763.5078195, 31724744937.1057];
if (s > 5.210887825 & s < 9.15546555571324)==1
Sigma = s / 9.2;
eta = h / 2800;
teta = 0;
for i = 1 : 36
teta = teta + ni(i) * (eta - 0.119) ^ Ii(i) * (Sigma - 1.07) ^ Ji(i);
end
T4_hs = teta * 550;
else
%function psat_h
if (s > -0.0001545495919 && s <= 3.77828134)==1
Low_Bound = 0.000611;
High_Bound = 165.291642526045;
hL=-1000;
while abs(hL - h) > 0.00001 && abs(High_Bound - Low_Bound) > 0.0001
PL = (Low_Bound + High_Bound) / 2;
Ts = T4_p(PL);
hL = h1_pT(PL, Ts);
if hL > h
High_Bound = PL;
else
Low_Bound = PL;
end
end
elseif s > 3.77828134 && s <= 4.41202148223476
PL = p3sat_h(h);
elseif s > 4.41202148223476 && s <= 5.210887663
PL = p3sat_h(h);
end
Low_Bound = 0.000611;
High_Bound = PL;
sss=-1000;
while (abs(s - sss) > 0.000001 & abs(High_Bound - Low_Bound) > 0.0000001)==1
p = (Low_Bound + High_Bound) / 2;
%Calculate s4_ph
Ts = T4_p(p);
xs = x4_ph(p, h);
if p < 16.529
s4v = s2_pT(p, Ts);
s4L = s1_pT(p, Ts);
else
v4v = v3_ph(p, h4V_p(p));
s4v = s3_rhoT(1 / v4v, Ts);
v4L = v3_ph(p, h4L_p(p));
s4L = s3_rhoT(1 / v4L, Ts);
end
sss = (xs * s4v + (1 - xs) * s4L);
if sss < s
High_Bound = p;
else
Low_Bound = p;
end
end
T4_hs = T4_p(p);
end
%***********************************************************************************************************
%*2.5 functions for region 5
function h5_pT = h5_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Basic Equation for Region 5
%Eq 32,33, Page 36, Tables 37-41
Ji0 = [0, 1, -3, -2, -1, 2];
ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917];
Iir = [1, 1, 1, 2, 3];
Jir = [0, 1, 3, 9, 3];
nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07];
R = 0.461526; %kJ/(kg K)
tau = 1000 / T;
Pi = p;
gamma0_tau = 0;
for i = 1 : 6
gamma0_tau = gamma0_tau + ni0(i) * Ji0(i) * tau ^ (Ji0(i) - 1);
end
gammar_tau = 0;
for i = 1 : 5
gammar_tau = gammar_tau + nir(i) * Pi ^ Iir(i) * Jir(i) * tau ^ (Jir(i) - 1);
end
h5_pT = R * T * tau * (gamma0_tau + gammar_tau);
function v5_pT = v5_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Basic Equation for Region 5
%Eq 32,33, Page 36, Tables 37-41
Ji0 = [0, 1, -3, -2, -1, 2];
ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917];
Iir = [1, 1, 1, 2, 3];
Jir = [0, 1, 3, 9, 3];
nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07];
R = 0.461526; %kJ/(kg K)
tau = 1000 / T;
Pi = p;
gamma0_pi = 1 / Pi;
gammar_pi = 0;
for i = 1 : 5
gammar_pi = gammar_pi + nir(i) * Iir(i) * Pi ^ (Iir(i) - 1) * tau ^ Jir(i);
end
v5_pT = R * T / p * Pi * (gamma0_pi + gammar_pi) / 1000;
function u5_pT = u5_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Basic Equation for Region 5
%Eq 32,33, Page 36, Tables 37-41
Ji0 = [0, 1, -3, -2, -1, 2];
ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917];
Iir = [1, 1, 1, 2, 3];
Jir = [0, 1, 3, 9, 3];
nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07];
R = 0.461526; %kJ/(kg K)
tau = 1000 / T;
Pi = p;
gamma0_pi = 1 / Pi;
gamma0_tau = 0;
for i = 1 : 6
gamma0_tau = gamma0_tau + ni0(i) * Ji0(i) * tau ^ (Ji0(i) - 1);
end
gammar_pi = 0;
gammar_tau = 0;
for i = 1 : 5
gammar_pi = gammar_pi + nir(i) * Iir(i) * Pi ^ (Iir(i) - 1) * tau ^ Jir(i);
gammar_tau = gammar_tau + nir(i) * Pi ^ Iir(i) * Jir(i) * tau ^ (Jir(i) - 1);
end
u5_pT = R * T * (tau * (gamma0_tau + gammar_tau) - Pi * (gamma0_pi + gammar_pi));
function Cp5_pT = Cp5_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Basic Equation for Region 5
%Eq 32,33, Page 36, Tables 37-41
Ji0 = [0, 1, -3, -2, -1, 2];
ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917];
Iir = [1, 1, 1, 2, 3];
Jir = [0, 1, 3, 9, 3];
nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07];
R = 0.461526; %kJ/(kg K)
tau = 1000 / T;
Pi = p;
gamma0_tautau = 0;
for i = 1 : 6
gamma0_tautau = gamma0_tautau + ni0(i) * Ji0(i) * (Ji0(i) - 1) * tau ^ (Ji0(i) - 2);
end
gammar_tautau = 0;
for i = 1 : 5
gammar_tautau = gammar_tautau + nir(i) * Pi ^ Iir(i) * Jir(i) * (Jir(i) - 1) * tau ^ (Jir(i) - 2);
end
Cp5_pT = -R * tau ^ 2 * (gamma0_tautau + gammar_tautau);
function s5_pT = s5_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Basic Equation for Region 5
%Eq 32,33, Page 36, Tables 37-41
Ji0 = [0, 1, -3, -2, -1, 2];
ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917];
Iir = [1, 1, 1, 2, 3];
Jir = [0, 1, 3, 9, 3];
nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07];
R = 0.461526; %kJ/(kg K)
tau = 1000 / T;
Pi = p;
gamma0 = log(Pi);
gamma0_tau = 0;
for i = 1 : 6
gamma0_tau = gamma0_tau + ni0(i) * Ji0(i) * tau ^ (Ji0(i) - 1);
gamma0 = gamma0 + ni0(i) * tau ^ Ji0(i);
end
gammar = 0;
gammar_tau = 0;
for i = 1 : 5
gammar = gammar + nir(i) * Pi ^ Iir(i) * tau ^ Jir(i);
gammar_tau = gammar_tau + nir(i) * Pi ^ Iir(i) * Jir(i) * tau ^ (Jir(i) - 1);
end
s5_pT = R * (tau * (gamma0_tau + gammar_tau) - (gamma0 + gammar));
function Cv5_pT = Cv5_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Basic Equation for Region 5
%Eq 32,33, Page 36, Tables 37-41
Ji0 = [0, 1, -3, -2, -1, 2];
ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917];
Iir = [1, 1, 1, 2, 3];
Jir = [0, 1, 3, 9, 3];
nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07];
R = 0.461526; %kJ/(kg K)
tau = 1000 / T;
Pi = p;
gamma0_tautau = 0;
for i = 1 : 6
gamma0_tautau = gamma0_tautau + ni0(i) * (Ji0(i) - 1) * Ji0(i) * tau ^ (Ji0(i) - 2);
end
gammar_pi = 0;
gammar_pitau = 0;
gammar_pipi = 0;
gammar_tautau = 0;
for i = 1 : 5
gammar_pi = gammar_pi + nir(i) * Iir(i) * Pi ^ (Iir(i) - 1) * tau ^ Jir(i);
gammar_pitau = gammar_pitau + nir(i) * Iir(i) * Pi ^ (Iir(i) - 1) * Jir(i) * tau ^ (Jir(i) - 1);
gammar_pipi = gammar_pipi + nir(i) * Iir(i) * (Iir(i) - 1) * Pi ^ (Iir(i) - 2) * tau ^ Jir(i);
gammar_tautau = gammar_tautau + nir(i) * Pi ^ Iir(i) * Jir(i) * (Jir(i) - 1) * tau ^ (Jir(i) - 2);
end
Cv5_pT = R * (-(tau ^ 2 *(gamma0_tautau + gammar_tautau)) - (1 + Pi * gammar_pi - tau * Pi * gammar_pitau)^2 / (1 - Pi ^ 2 * gammar_pipi));
function w5_pT = w5_pT(p, T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam, September 1997
%Basic Equation for Region 5
%Eq 32,33, Page 36, Tables 37-41
Ji0 = [0, 1, -3, -2, -1, 2];
ni0 = [-13.179983674201, 6.8540841634434, -0.024805148933466, 0.36901534980333, -3.1161318213925, -0.32961626538917];
Iir = [1, 1, 1, 2, 3];
Jir = [0, 1, 3, 9, 3];
nir = [-1.2563183589592E-04, 2.1774678714571E-03, -0.004594282089991, -3.9724828359569E-06, 1.2919228289784E-07];
R = 0.461526; %kJ/(kg K)
tau = 1000 / T;
Pi = p;
gamma0_tautau = 0;
for i = 1 : 6
gamma0_tautau = gamma0_tautau + ni0(i) * (Ji0(i) - 1) * Ji0(i) * tau ^ (Ji0(i) - 2);
end
gammar_pi = 0;
gammar_pitau = 0;
gammar_pipi = 0;
gammar_tautau = 0;
for i = 1 : 5
gammar_pi = gammar_pi + nir(i) * Iir(i) * Pi ^ (Iir(i) - 1) * tau ^ Jir(i);
gammar_pitau = gammar_pitau + nir(i) * Iir(i) * Pi ^ (Iir(i) - 1) * Jir(i) * tau ^ (Jir(i) - 1);
gammar_pipi = gammar_pipi + nir(i) * Iir(i) * (Iir(i) - 1) * Pi ^ (Iir(i) - 2) * tau ^ Jir(i);
gammar_tautau = gammar_tautau + nir(i) * Pi ^ Iir(i) * Jir(i) * (Jir(i) - 1) * tau ^ (Jir(i) - 2);
end
w5_pT = (1000 * R * T * (1 + 2 * Pi * gammar_pi + Pi ^ 2 * gammar_pi ^ 2) / ((1 - Pi ^ 2 * gammar_pipi) + (1 + Pi * gammar_pi - tau * Pi * gammar_pitau) ^ 2 / (tau ^ 2 * (gamma0_tautau + gammar_tautau)))) ^ 0.5;
function T5_ph = T5_ph(p, h)
%Solve with half interval method
Low_Bound = 1073.15;
High_Bound = 2273.15;
hs=h-1;
while abs(h - hs) > 0.00001
Ts = (Low_Bound + High_Bound) / 2;
hs = h5_pT(p, Ts);
if hs > h
High_Bound = Ts;
else
Low_Bound = Ts;
end
end
T5_ph = Ts;
function T5_ps = T5_ps(p, s)
%Solve with half interval method
Low_Bound = 1073.15;
High_Bound = 2273.15;
ss=s-1;
while abs(s - ss) > 0.00001
Ts = (Low_Bound + High_Bound) / 2;
ss = s5_pT(p, Ts);
if ss > s
High_Bound = Ts;
else
Low_Bound = Ts;
end
end
T5_ps = Ts;
function T5_prho=T5_prho(p,rho)
%Solve by iteration. Observe that fo low temperatures this equation has 2 solutions.
%Solve with half interval method
Low_Bound = 1073.15;
High_Bound = 2073.15;
rhos=-1000;
while abs(rho - rhos) > 0.000001
Ts = (Low_Bound + High_Bound) / 2;
rhos = 1 / v2_pT(p, Ts);
if rhos < rho
High_Bound = Ts;
else
Low_Bound = Ts;
end
end
T5_prho = Ts;
%***********************************************************************************************************
%*3 Region Selection
%***********************************************************************************************************
%*3.1 Regions as a function of pT
function region_pT = region_pT(p, T)
if T > 1073.15 && p < 10 && T < 2273.15 && p > 0.000611
region_pT = 5;
elseif T <= 1073.15 && T > 273.15 && p <= 100 && p > 0.000611
if T > 623.15
if p > B23p_T(T)
region_pT = 3;
if T < 647.096
ps = p4_T(T);
if abs(p - ps) < 0.00001
region_pT = 4;
end
end
else
region_pT = 2;
end
else
ps = p4_T(T);
if abs(p - ps) < 0.00001
region_pT = 4;
elseif p > ps
region_pT = 1;
else
region_pT = 2;
end
end
else
region_pT = 0; %**Error, Outside valid area
end
%***********************************************************************************************************
%*3.2 Regions as a function of ph
function region_ph = region_ph( p, h)
%Check if outside pressure limits
if p < 0.000611657 || p > 100
region_ph = 0;
return
end
%Check if outside low h.
if h < 0.963 * p + 2.2 %Linear adaption to h1_pt()+2 to speed up calcualations.
if h < h1_pT(p, 273.15)
region_ph = 0;
return
end
end
if p < 16.5292 %Bellow region 3,Check region 1,4,2,5
%Check Region 1
Ts = T4_p(p);
hL = 109.6635 * log(p) + 40.3481 * p + 734.58; %Approximate function for hL_p
if abs(h - hL) < 100 %if approximate is not god enough use real function
hL = h1_pT(p, Ts);
end
if h <= hL
region_ph = 1;
return
end
%Check Region 4
hV = 45.1768 * log(p) - 20.158 * p + 2804.4; %Approximate function for hV_p
if abs(h - hV) < 50 %if approximate is not god enough use real function
hV = h2_pT(p, Ts);
end
if h < hV
region_ph = 4;
return
end
%Check upper limit of region 2 Quick Test
if h < 4000
region_ph = 2;
return
end
%Check region 2 (Real value)
h_45 = h2_pT(p, 1073.15);
if h <= h_45
region_ph = 2;
return
end
%Check region 5
if p > 10
region_ph = 0;
return
end
h_5u = h5_pT(p, 2273.15);
if h < h_5u
region_ph = 5;
return
end
region_ph = 0;
return
else %for p>16.5292
%Check if in region1
if h < h1_pT(p, 623.15)
region_ph = 1;
return
end
%Check if in region 3 or 4 (Bellow Reg 2)
if h < h2_pT(p, B23T_p(p))
%Region 3 or 4
if p > p3sat_h(h)
region_ph = 3;
return
else
region_ph = 4;
return
end
end
%Check if region 2
if h < h2_pT(p, 1073.15)
region_ph = 2;
return
end
end
region_ph = 0;
%***********************************************************************************************************
%*3.3 Regions as a function of ps
function region_ps = region_ps( p, s)
if p < 0.000611657 || p > 100 || s < 0 || s > s5_pT(p, 2273.15)
region_ps = 0;
return
end
%Check region 5
if s > s2_pT(p, 1073.15)
if p <= 10
region_ps = 5;
return
else
region_ps = 0;
return
end
end
%Check region 2
if p > 16.529
ss = s2_pT(p, B23T_p(p)); %Between 5.047 & 5.261. Use to speed up!
else
ss = s2_pT(p, T4_p(p));
end
if s > ss
region_ps = 2;
return
end
%Check region 3
ss = s1_pT(p, 623.15);
if p > 16.529 && s > ss
if p > p3sat_s(s)
region_ps = 3;
return
else
region_ps = 4;
return
end
end
%Check region 4 (Not inside region 3)
if p < 16.529 && s > s1_pT(p, T4_p(p))
region_ps = 4;
return
end
%Check region 1
if p > 0.000611657 && s > s1_pT(p, 273.15)
region_ps = 1;
return
end
region_ps = 1;
%***********************************************************************************************************
%*3.4 Regions as a function of hs
function region_hs = region_hs( h, s)
if s < -0.0001545495919
region_hs = 0;
return
end
%Check linear adaption to p=0.000611. if bellow region 4.
hMin = (((-0.0415878 - 2500.89262) / (-0.00015455 - 9.155759)) * s);
if s < 9.155759395 && h < hMin
region_hs = 0;
return
end
%******Kolla 1 eller 4. (+liten bit över B13)
if s >= -0.0001545495919 && s <= 3.77828134
if h < h4_s(s)
region_hs = 4;
return
elseif s < 3.397782955 %100MPa line is limiting
TMax = T1_ps(100, s);
hMax = h1_pT(100, TMax);
if h < hMax
region_hs = 1;
return
else
region_hs = 0;
return
end
else %The point is either in region 4,1,3. Check B23
hB = hB13_s(s);
if h < hB
region_hs = 1;
return
end
TMax = T3_ps(100, s);
vmax = v3_ps(100, s);
hMax = h3_rhoT(1 / vmax, TMax);
if h < hMax
region_hs = 3;
return
else
region_hs = 0;
return
end
end
end
%******Kolla region 2 eller 4. (Övre delen av område b23-> max)
if s >= 5.260578707 && s <= 11.9212156897728
if s > 9.155759395 %Above region 4
Tmin = T2_ps(0.000611, s);
hMin = h2_pT(0.000611, Tmin);
%function adapted to h(1073.15,s)
hMax = -0.07554022 * s ^ 4 + 3.341571 * s ^ 3 - 55.42151 * s ^ 2 + 408.515 * s + 3031.338;
if h > hMin && h < hMax
region_hs = 2;
return
else
region_hs = 0;
return
end
end
hV = h4_s(s);
if h < hV %Region 4. Under region 3.
region_hs = 4;
return
end
if s < 6.04048367171238
TMax = T2_ps(100, s);
hMax = h2_pT(100, TMax);
else
%function adapted to h(1073.15,s)
hMax = -2.988734 * s ^ 4 + 121.4015 * s ^ 3 - 1805.15 * s ^ 2 + 11720.16 * s - 23998.33;
end
if h < hMax %Region 2. Över region 4.
region_hs = 2;
return
else
region_hs = 0;
return
end
end
%Kolla region 3 eller 4. Under kritiska punkten.
if s >= 3.77828134 && s <= 4.41202148223476
hL = h4_s(s);
if h < hL
region_hs = 4;
return
end
TMax = T3_ps(100, s);
vmax = v3_ps(100, s);
hMax = h3_rhoT(1 / vmax, TMax);
if h < hMax
region_hs = 3;
return
else
region_hs = 0;
return
end
end
%Kolla region 3 eller 4 från kritiska punkten till övre delen av b23
if s >= 4.41202148223476 && s <= 5.260578707
hV = h4_s(s);
if h < hV
region_hs = 4;
return
end
%Kolla om vi är under b23 giltighetsområde.
if s <= 5.048096828
TMax = T3_ps(100, s);
vmax = v3_ps(100, s);
hMax = h3_rhoT(1 / vmax, TMax);
if h < hMax
region_hs = 3;
return
else
region_hs = 0;
return
end
else %Inom området för B23 i s led.
if (h > 2812.942061) %Ovanför B23 i h_led
if s > 5.09796573397125
TMax = T2_ps(100, s);
hMax = h2_pT(100, TMax);
if h < hMax
region_hs = 2;
return
else
region_hs = 0;
return
end
else
region_hs = 0;
return
end
end
if (h < 2563.592004) %Nedanför B23 i h_led men vi har redan kollat ovanför hV2c3b
region_hs = 3;
return
end
%Vi är inom b23 området i både s och h led.
Tact = TB23_hs(h, s);
pact = p2_hs(h, s);
pBound = B23p_T(Tact);
if pact > pBound
region_hs = 3;
return
else
region_hs = 2;
return
end
end
end
region_hs = 0;
%***********************************************************************************************************
%*3.5 Regions as a function of p and rho
function Region_prho = Region_prho(p,rho)
v = 1 / rho;
if p < 0.000611657 || p > 100
Region_prho = 0;
return
end
if p < 16.5292 %Bellow region 3, Check region 1,4,2
if v < v1_pT(p, 273.15) %Observe that this is not actually min of v. Not valid Water of 4°C is ligther.
Region_prho = 0;
return
end
if v <= v1_pT(p, T4_p(p))
Region_prho = 1;
return
end
if v < v2_pT(p, T4_p(p))
Region_prho = 4;
return
end
if v <= v2_pT(p, 1073.15)
Region_prho = 2;
return
end
if p > 10 %Above region 5
Region_prho = 0;
return
end
if v <= v5_pT(p, 2073.15)
Region_prho = 5;
return
end
else %Check region 1,3,4,3,2 (Above the lowest point of region 3.)
if v < v1_pT(p, 273.15) %Observe that this is not actually min of v. Not valid Water of 4°C is ligther.
Region_prho = 0;
return
end
if v < v1_pT(p, 623.15)
Region_prho = 1;
return
end
%Check if in region 3 or 4 (Bellow Reg 2)
if v < v2_pT(p, B23T_p(p))
%Region 3 or 4
if p > 22.064 %Above region 4
Region_prho = 3;
return
end
if v < v3_ph(p, h4L_p(p)) || v > v3_ph(p, h4V_p(p)) %Uses iteration!!
Region_prho = 3;
return
else
Region_prho = 4;
return
end
end
%Check if region 2
if v < v2_pT(p, 1073.15)
Region_prho = 2;
return
end
end
Region_prho = 0;
%***********************************************************************************************************
%*4 Region Borders
%***********************************************************************************************************
%***********************************************************************************************************
%*4.1 Boundary between region 2 and 3.
function B23p_T = B23p_T(T)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%1997
%Section 4 Auxiliary Equation for the Boundary between Regions 2 and 3
%Eq 5, Page 5
B23p_T = 348.05185628969 - 1.1671859879975 * T + 1.0192970039326E-03 * T ^ 2;
function B23T_p = B23T_p(p)
%Release on the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water and Steam
%1997
%Section 4 Auxiliary Equation for the Boundary between Regions 2 and 3
%Eq 6, Page 6
B23T_p = 572.54459862746 + ((p - 13.91883977887) / 1.0192970039326E-03) ^ 0.5;
%***********************************************************************************************************
%*4.2 Region 3. pSat_h & pSat_s
function p3sat_h = p3sat_h(h)
%Revised Supplementary Release on Backward Equations for the functions T(p,h), v(p,h) & T(p,s), v(p,s) for Region 3 of the IAPWS Industrial formulation 1997 for the Thermodynamic Properties of Water & Steam
%2004
%Section 4 Boundary Equations psat(h) & psat(s) for the Saturation Lines of Region 3
%Se pictures Page 17, Eq 10, Table 17, Page 18
Ii = [0, 1, 1, 1, 1, 5, 7, 8, 14, 20, 22, 24, 28, 36];
Ji = [0, 1, 3, 4, 36, 3, 0, 24, 16, 16, 3, 18, 8, 24];
ni = [0.600073641753024, -9.36203654849857, 24.6590798594147, -107.014222858224, -91582131580576.8, -8623.32011700662, -23.5837344740032, 2.52304969384128E+17, -3.89718771997719E+18, -3.33775713645296E+22, 35649946963.6328, -1.48547544720641E+26, 3.30611514838798E+18, 8.13641294467829E+37];
hs = h / 2600;
ps = 0;
for i = 1:14
ps = ps + ni(i) * (hs - 1.02) ^ Ii(i) * (hs - 0.608) ^ Ji(i);
end
p3sat_h = ps * 22;
function p3sat_s = p3sat_s(s)
Ii = [0, 1, 1, 4, 12, 12, 16, 24, 28, 32];
Ji = [0, 1, 32, 7, 4, 14, 36, 10, 0, 18];
ni = [0.639767553612785, -12.9727445396014, -2.24595125848403E+15, 1774667.41801846, 7170793495.71538, -3.78829107169011E+17, -9.55586736431328E+34, 1.87269814676188E+23, 119254746466.473, 1.10649277244882E+36];
Sigma = s / 5.2;
Pi = 0;
for i = 1:10
Pi = Pi + ni(i) * (Sigma - 1.03) ^ Ii(i) * (Sigma - 0.699) ^ Ji(i);
end
p3sat_s = Pi * 22;
%***********************************************************************************************************
%4.3 Region boundary 1to3 & 3to2 as a functions of s
function hB13_s = hB13_s(s)
%Supplementary Release on Backward Equations ( ) , p h s for Region 3,
%Chapter 4.5 page 23.
Ii = [0, 1, 1, 3, 5, 6];
Ji = [0, -2, 2, -12, -4, -3];
ni = [0.913965547600543, -4.30944856041991E-05, 60.3235694765419, 1.17518273082168E-18, 0.220000904781292, -69.0815545851641];
Sigma = s / 3.8;
eta = 0;
for i = 1 : 6
eta = eta + ni(i) * (Sigma - 0.884) ^ Ii(i) * (Sigma - 0.864) ^ Ji(i);
end
hB13_s = eta * 1700;
function TB23_hs = TB23_hs(h, s)
%Supplementary Release on Backward Equations ( ) , p h s for Region 3,
%Chapter 4.6 page 25.
Ii = [-12, -10, -8, -4, -3, -2, -2, -2, -2, 0, 1, 1, 1, 3, 3, 5, 6, 6, 8, 8, 8, 12, 12, 14, 14];
Ji = [10, 8, 3, 4, 3, -6, 2, 3, 4, 0, -3, -2, 10, -2, -1, -5, -6, -3, -8, -2, -1, -12, -1, -12, 1];
ni = [6.2909626082981E-04, -8.23453502583165E-04, 5.15446951519474E-08, -1.17565945784945, 3.48519684726192, -5.07837382408313E-12, -2.84637670005479, -2.36092263939673, 6.01492324973779, 1.48039650824546, 3.60075182221907E-04, -1.26700045009952E-02, -1221843.32521413, 0.149276502463272, 0.698733471798484, -2.52207040114321E-02, 1.47151930985213E-02, -1.08618917681849, -9.36875039816322E-04, 81.9877897570217, -182.041861521835, 2.61907376402688E-06, -29162.6417025961, 1.40660774926165E-05, 7832370.62349385];
Sigma = s / 5.3;
eta = h / 3000;
teta = 0;
for i = 1 : 25
teta = teta + ni(i) * (eta - 0.727) ^ Ii(i) * (Sigma - 0.864) ^ Ji(i);
end
TB23_hs = teta * 900;
%***********************************************************************************************************
%*5 Transport properties
%***********************************************************************************************************
%*5.1 Viscosity (IAPWS formulation 1985, Revised 2003)
%***********************************************************************************************************
function my_AllRegions_pT = my_AllRegions_pT( p, T)
h0 = [0.5132047, 0.3205656, 0, 0, -0.7782567, 0.1885447];
h1 = [0.2151778, 0.7317883, 1.241044, 1.476783, 0, 0];
h2 = [-0.2818107, -1.070786, -1.263184, 0, 0, 0];
h3 = [0.1778064, 0.460504, 0.2340379, -0.4924179, 0, 0];
h4 = [-0.0417661, 0, 0, 0.1600435, 0, 0];
h5 = [0, -0.01578386, 0, 0, 0, 0];
h6 = [0, 0, 0, -0.003629481, 0, 0];
%Calcualte density.
switch region_pT(p, T)
case 1
rho = 1 / v1_pT(p, T);
case 2
rho = 1 / v2_pT(p, T);
case 3
hs = h3_pT(p, T);
rho = 1 / v3_ph(p, hs);
case 4
rho = NaN;
case 5
rho = 1 / v5_pT(p, T);
otherwise
my_AllRegions_pT = NaN;
return
end
rhos = rho / 317.763;
Ts = T / 647.226;
ps = p / 22.115;
%Check valid area
if T > 900 + 273.15 || (T > 600 + 273.15 && p > 300) || (T > 150 + 273.15 && p > 350) || p > 500
my_AllRegions_pT = NaN;
return
end
my0 = Ts ^ 0.5 / (1 + 0.978197 / Ts + 0.579829 / (Ts ^ 2) - 0.202354 / (Ts ^ 3));
Sum = 0;
for i = 0 : 5
Sum = Sum + h0(i+1) * (1 / Ts - 1) ^ i + h1(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 1 + h2(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 2 + h3(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 3 + h4(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 4 + h5(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 5 + h6(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 6;
end
my1 = exp(rhos * Sum);
mys = my0 * my1;
my_AllRegions_pT = mys * 0.000055071;
function my_AllRegions_ph = my_AllRegions_ph( p, h)
h0 = [0.5132047, 0.3205656, 0, 0, -0.7782567, 0.1885447];
h1 = [0.2151778, 0.7317883, 1.241044, 1.476783, 0, 0];
h2 = [-0.2818107, -1.070786, -1.263184, 0, 0, 0];
h3 = [0.1778064, 0.460504, 0.2340379, -0.4924179, 0, 0];
h4 = [-0.0417661, 0, 0, 0.1600435, 0, 0];
h5 = [0, -0.01578386, 0, 0, 0, 0];
h6 = [0, 0, 0, -0.003629481, 0, 0];
%Calcualte density.
switch region_ph(p, h)
case 1
Ts = T1_ph(p, h);
T = Ts;
rho = 1 / v1_pT(p, Ts);
case 2
Ts = T2_ph(p, h);
T = Ts;
rho = 1 / v2_pT(p, Ts);
case 3
rho = 1 / v3_ph(p, h);
T = T3_ph(p, h);
case 4
xs = x4_ph(p, h);
if p < 16.529
v4v = v2_pT(p, T4_p(p));
v4L = v1_pT(p, T4_p(p));
else
v4v = v3_ph(p, h4V_p(p));
v4L = v3_ph(p, h4L_p(p));
end
rho = 1 / (xs * v4v + (1 - xs) * v4L);
T = T4_p(p);
case 5
Ts = T5_ph(p, h);
T = Ts;
rho = 1 / v5_pT(p, Ts);
otherwise
my_AllRegions_ph = NaN;
return
end
rhos = rho / 317.763;
Ts = T / 647.226;
ps = p / 22.115;
%Check valid area
if T > 900 + 273.15 || (T > 600 + 273.15 && p > 300) || (T > 150 + 273.15 && p > 350) || p > 500
my_AllRegions_ph = NaN;
return
end
my0 = Ts ^ 0.5 / (1 + 0.978197 / Ts + 0.579829 / (Ts ^ 2) - 0.202354 / (Ts ^ 3));
Sum = 0;
for i = 0 : 5
Sum = Sum + h0(i+1) * (1 / Ts - 1) ^ i + h1(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 1 + h2(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 2 + h3(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 3 + h4(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 4 + h5(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 5 + h6(i+1) * (1 / Ts - 1) ^ i * (rhos - 1) ^ 6;
end
my1 = exp(rhos * Sum);
mys = my0 * my1;
my_AllRegions_ph = mys * 0.000055071;
%***********************************************************************************************************
%*5.2 Thermal Conductivity (IAPWS formulation 1985)
function tc_ptrho = tc_ptrho( p, T, rho)
%Revised release on the IAPWS formulation 1985 for the Thermal Conductivity of ordinary water
%IAPWS September 1998
%Page 8
%ver2.6 Start corrected bug
if T < 273.15
tc_ptrho = NaN; %Out of range of validity (para. B4)
return
elseif T < 500 + 273.15
if p > 100
tc_ptrho = NaN; %Out of range of validity (para. B4)
return
end
elseif T <= 650 + 273.15
if p > 70
tc_ptrho = NaN; %Out of range of validity (para. B4)
return
end
else T <= 800 + 273.15
if p > 40
tc_ptrho = NaN; %Out of range of validity (para. B4)
return
end
end
%ver2.6 End corrected bug
T = T / 647.26;
rho = rho / 317.7;
tc0 = T ^ 0.5 * (0.0102811 + 0.0299621 * T + 0.0156146 * T ^ 2 - 0.00422464 * T ^ 3);
tc1 = -0.39707 + 0.400302 * rho + 1.06 * exp(-0.171587 * (rho + 2.39219) ^ 2);
dT = abs(T - 1) + 0.00308976;
Q = 2 + 0.0822994 / dT ^ (3 / 5);
if T >= 1
s = 1 / dT;
else
s = 10.0932 / dT ^ (3 / 5);
end
tc2 = (0.0701309 / T ^ 10 + 0.011852) * rho ^ (9 / 5) * exp(0.642857 * (1 - rho ^ (14 / 5))) + 0.00169937 * s * rho ^ Q * exp((Q / (1 + Q)) * (1 - rho ^ (1 + Q))) - 1.02 * exp(-4.11717 * T ^ (3 / 2) - 6.17937 / rho ^ 5);
tc_ptrho = tc0 + tc1 + tc2;
%***********************************************************************************************************
%5.3 Surface Tension
function Surface_Tension_T = Surface_Tension_T( T)
%IAPWS Release on Surface Tension of Ordinary Water Substance,
%September 1994
tc = 647.096; %K
B = 0.2358; %N/m
bb = -0.625;
my = 1.256;
if T < 0.01 || T > tc
Surface_Tension_T = NaN; %"Out of valid region"
return
end
tau = 1 - T / tc;
Surface_Tension_T = B * tau ^ my * (1 + bb * tau);
%***********************************************************************************************************
%*6 Units *
%***********************************************************************************************************
function toSIunit_p = toSIunit_p( Ins )
%Translate bar to MPa
toSIunit_p = Ins / 10;
function fromSIunit_p = fromSIunit_p( Ins )
%Translate bar to MPa
fromSIunit_p = Ins * 10;
function toSIunit_T = toSIunit_T( Ins )
%Translate degC to Kelvon
toSIunit_T = Ins + 273.15;
function fromSIunit_T = fromSIunit_T( Ins )
%Translate Kelvin to degC
fromSIunit_T = Ins - 273.15;
function toSIunit_h = toSIunit_h( Ins )
toSIunit_h = Ins;
function fromSIunit_h = fromSIunit_h( Ins )
fromSIunit_h = Ins;
function toSIunit_v = toSIunit_v( Ins )
toSIunit_v = Ins;
function fromSIunit_v = fromSIunit_v( Ins )
fromSIunit_v = Ins;
function toSIunit_s = toSIunit_s( Ins )
toSIunit_s = Ins;
function fromSIunit_s = fromSIunit_s( Ins )
fromSIunit_s = Ins;
function toSIunit_u = toSIunit_u( Ins )
toSIunit_u = Ins;
function fromSIunit_u = fromSIunit_u( Ins )
fromSIunit_u = Ins;
function toSIunit_Cp = toSIunit_Cp( Ins )
toSIunit_Cp = Ins;
function fromSIunit_Cp = fromSIunit_Cp( Ins )
fromSIunit_Cp = Ins;
function toSIunit_Cv = toSIunit_Cv( Ins )
toSIunit_Cv = Ins;
function fromSIunit_Cv = fromSIunit_Cv( Ins )
fromSIunit_Cv = Ins;
function toSIunit_w = toSIunit_w( Ins )
toSIunit_w = Ins;
function fromSIunit_w = fromSIunit_w( Ins )
fromSIunit_w = Ins;
function toSIunit_tc = toSIunit_tc( Ins )
toSIunit_tc = Ins;
function fromSIunit_tc = fromSIunit_tc( Ins )
fromSIunit_tc = Ins;
function toSIunit_st = toSIunit_st( Ins )
toSIunit_st = Ins;
function fromSIunit_st = fromSIunit_st( Ins )
fromSIunit_st = Ins;
function toSIunit_x = toSIunit_x( Ins )
toSIunit_x = Ins;
function fromSIunit_x = fromSIunit_x( Ins )
fromSIunit_x = Ins;
function toSIunit_vx = toSIunit_vx( Ins )
toSIunit_vx = Ins;
function fromSIunit_vx = fromSIunit_vx( Ins )
fromSIunit_vx = Ins;
function toSIunit_my = toSIunit_my( Ins )
toSIunit_my = Ins;
function fromSIunit_my = fromSIunit_my( Ins )
fromSIunit_my = Ins;
%***********************************************************************************************************
%*7 Verification *
%***********************************************************************************************************
function err = check()
err=0;
%*********************************************************************************************************
%* 7.1 Verifiy region 1
%IF-97 Table 5, Page 9
p=[30/10,800/10,30/10];
T=[300,300,500];
Fun={'v1_pT','h1_pT','u1_pT','s1_pT','Cp1_pT','w1_pT'};
IF97=[0.00100215168,0.000971180894,0.001202418;...
115.331273,184.142828,975.542239;...
112.324818,106.448356,971.934985;...
0.392294792,0.368563852,2.58041912;...
4.17301218,4.01008987,4.65580682;...
1507.73921,1634.69054,1240.71337];
R1=zeros(6,3);
for i=1:3
for j=1:6
R1(j,i)=eval([char(Fun(j)),'(',num2str(p(i)),',',num2str(T(i)),');']);
end
end
Region1_error=abs((R1-IF97)./IF97)
err=err+sum(sum(Region1_error>1E-8))
%IF-97 Table 7, Page 11
p=[30/10,800/10,800/10];
h=[500,500,1500];
IF97=[391.798509,378.108626,611.041229];
R1=zeros(1,3);
for i=1:3
R1(i)=T1_ph(p(i),h(i));
end
T1_ph_error=abs((R1-IF97)./IF97)
err=err+sum(sum(T1_ph_error>1E-8))
%IF-97 Table 9, Page 12
p=[30/10,800/10,800/10];
s=[0.5,0.5,3];
IF97=[307.842258,309.979785,565.899909];
R1=zeros(1,3);
for i=1:3
R1(i)=T1_ps(p(i),s(i));
end
T1_ps_error=abs((R1-IF97)./IF97)
err=err+sum(sum(T1_ps_error>1E-8))
%Supplementary Release on Backward Equations
%for Pressure as a Function of Enthalpy and Entropy p(h,s)
%Table 3, Page 6
h=[0.001,90,1500];
s=[0,0,3.4];
IF97=[0.0009800980612,91.929547272,58.68294423];
R1=zeros(1,3);
for i=1:3
R1(i)=p1_hs(h(i),s(i));
end
p1_hs_error=abs((R1-IF97)./IF97)
err=err+sum(sum(p1_hs_error>1E-8))
%*********************************************************************************************************
%* 7.2 Verifiy region 2
% IF-97 Table 15, Page 17
p=[0.035/10,0.035/10,300/10];
T=[300,700,700];
Fun={'v2_pT','h2_pT','u2_pT','s2_pT','Cp2_pT','w2_pT'};
IF97=[39.4913866,92.3015898,0.00542946619;...
2549.91145,3335.68375,2631.49474;...
2411.6916,3012.62819,2468.61076;...
8.52238967,10.1749996,5.17540298;...
1.91300162,2.08141274,10.3505092;...
427.920172,644.289068,480.386523];
R2=zeros(6,3);
for i=1:3
for j=1:6
R2(j,i)=eval([char(Fun(j)),'(',num2str(p(i)),',',num2str(T(i)),');']);
end
end
Region2_error=abs((R2-IF97)./IF97)
err=err+sum(sum(Region2_error>1E-8))
%IF-97 Table 24, Page 25
p=[0.01/10,30/10,30/10,50/10,50/10,250/10,400/10,600/10,600/10];
h=[3000,3000,4000,3500,4000,3500,2700,2700,3200];
IF97=[534.433241,575.37337,1010.77577,801.299102,1015.31583,875.279054,743.056411,791.137067,882.75686];
R2=zeros(1,9);
for i=1:9
R2(i)=T2_ph(p(i),h(i));
end
T2_ph_error=abs((R2-IF97)./IF97)
err=err+sum(sum(T2_ph_error>1E-8))
%IF-97 Table 29, Page 29
p=[1/10,1/10,25/10,80/10,80/10,900/10,200/10,800/10,800/10];
s=[7.5,8,8,6,7.5,6,5.75,5.25,5.75];
IF97=[399.517097,514.127081,1039.84917,600.48404,1064.95556,1038.01126,697.992849,854.011484,949.017998];
R2=zeros(1,9);
for i=1:9
R2(i)=T2_ps(p(i),s(i));
end
T2_ps_error=abs((R2-IF97)./IF97)
err=err+sum(sum(T2_ps_error>1E-8))
%Supplementary Release on Backward Equations for Pressure as a Function of Enthalpy and Entropy p(h,s)
%Table 3, Page 6
h=[2800,2800,4100,2800,3600,3600,2800,2800,3400];
s=[6.5,9.5,9.5,6,6,7,5.1,5.8,5.8];
IF97=[1.371012767,0.001879743844,0.1024788997,4.793911442,83.95519209,7.527161441,94.3920206,8.414574124,83.76903879];
R2=zeros(1,9);
for i=1:9
R2(i)=p2_hs(h(i),s(i));
end
p2_hs_error=abs((R2-IF97)./IF97)
err=err+sum(sum(p2_hs_error>1E-8))
%*********************************************************************************************************
%* 7.3 Verifiy region 3
% IF-97 Table 33, Page 32
T=[650,650,750];
rho=[500,200,500];
Fun={'p3_rhoT','h3_rhoT','u3_rhoT','s3_rhoT','Cp3_rhoT','w3_rhoT'};
IF97=[25.5837018,22.2930643,78.3095639;...
1863.43019,2375.12401,2258.68845;...
1812.26279,2263.65868,2102.06932;...
4.05427273,4.85438792,4.46971906;...
13.8935717,44.6579342,6.34165359;...
502.005554,383.444594,760.696041];
R3=zeros(6,3);
for i=1:3
for j=1:6
R3(j,i)=eval([char(Fun(j)),'(',num2str(rho(i)),',',num2str(T(i)),');']);
end
end
Region3_error=abs((R3-IF97)./IF97)
err=err+sum(sum(Region3_error>1E-8))
%T3_ph
p=[200/10,500/10,1000/10,200/10,500/10,1000/10];
h=[1700,2000,2100,2500,2400,2700];
IF97=[629.3083892,690.5718338,733.6163014,641.8418053,735.1848618,842.0460876];
R3=zeros(1,6);
for i=1:6
R3(i)=T3_ph(p(i),h(i));
end
T3_ph_error=abs((R3-IF97)./IF97)
err=err+sum(sum(T3_ph_error>1E-8))
%v3_ph
p=[200/10,500/10,1000/10,200/10,500/10,1000/10];
h=[1700,2000,2100,2500,2400,2700];
IF97=[0.001749903962,0.001908139035,0.001676229776,0.006670547043,0.0028012445,0.002404234998];
R3=zeros(1,6);
for i=1:6
R3(i)=v3_ph(p(i),h(i));
end
v3_ph_error=abs((R3-IF97)./IF97)
err=err+sum(sum(v3_ph_error>1E-7))
%T3_ps
p=[200/10,500/10,1000/10,200/10,500/10,1000/10];
s=[3.7,3.5,4,5,4.5,5];
IF97=[620.8841563,618.1549029,705.6880237,640.1176443,716.3687517,847.4332825];
R3=zeros(1,6);
for i=1:6
R3(i)=T3_ps(p(i),s(i));
end
T3_ps_error=abs((R3-IF97)./IF97)
err=err+sum(sum(T3_ps_error>1E-8))
%v3_ps
p=[200/10,500/10,1000/10,200/10,500/10,1000/10];
s=[3.7,3.5,4,5,4.5,5];
IF97=[0.001639890984,0.001423030205,0.001555893131,0.006262101987,0.002332634294,0.002449610757];
R3=zeros(1,6);
for i=1:6
R3(i)=v3_ps(p(i),s(i));
end
v3_ps_error=abs((R3-IF97)./IF97)
err=err+sum(sum(v3_ps_error>1E-8))
%p3_hs
h=[1700,2000,2100,2500,2400,2700];
s=[3.8,4.2,4.3,5.1,4.7,5];
IF97=[25.55703246,45.40873468,60.7812334,17.20612413,63.63924887,88.39043281];
R3=zeros(1,6);
for i=1:6
R3(i)=p3_hs(h(i),s(i));
end
p3_hs_error=abs((R3-IF97)./IF97)
err=err+sum(sum(p3_hs_error>1E-8))
%h3_pT (Iteration)
p=[255.83702,222.93064,783.09564]./10;
T=[650,650,750];
IF97=[1863.271389,2375.696155,2258.626582];
R3=zeros(1,3);
for i=1:3
R3(i)=h3_pT(p(i),T(i));
end
h3_pT_error=abs((R3-IF97)./IF97)
err=err+sum(sum(h3_pT_error>1E-6)) %Decimals in IF97
%*********************************************************************************************************
%* 7.4 Verifiy region 4
%Saturation pressure, If97, Table 35, Page 34
T=[300,500,600];
IF97=[0.0353658941, 26.3889776, 123.443146]/10;
R3=zeros(1,3);
for i=1:3
R4(i)=p4_T(T(i));
end
p4_t_error=abs((R4-IF97)./IF97)
err=err+sum(sum( p4_t_error>1E-7))
T=[1,10,100]/10;
IF97=[372.755919,453.035632,584.149488];
R3=zeros(1,3);
for i=1:3
R4(i)=T4_p(T(i));
end
T4_p_error=abs((R4-IF97)./IF97)
err=err+sum(sum( T4_p_error>1E-7))
s=[1,2,3,3.8,4,4.2,7,8,9,5.5,5,4.5];
IF97=[308.5509647,700.6304472,1198.359754,1685.025565,1816.891476,1949.352563,2723.729985,2599.04721,2511.861477,2687.69385,2451.623609,2144.360448];
R3=zeros(1,12);
for i=1:12
R4(i)=h4_s(s(i));
end
h4_s_error=abs((R4-IF97)./IF97)
err=err+sum(sum( h4_s_error>1E-7))
%*********************************************************************************************************
%* 7.5 Verifiy region 5
% IF-97 Table 42, Page 39
T=[1500,1500,2000];
p=[5,80,80]/10;
Fun={'v5_pT','h5_pT','u5_pT','s5_pT','Cp5_pT','w5_pT'};
IF97=[1.38455354,0.0865156616,0.115743146;...
5219.76332,5206.09634,6583.80291;...
4527.48654,4513.97105,5657.85774;...
9.65408431,8.36546724,9.15671044;...
2.61610228,2.64453866,2.8530675;...
917.071933,919.708859,1054.35806];
R5=zeros(6,3);
for i=1:3
for j=1:6
R5(j,i)=eval([char(Fun(j)),'(',num2str(p(i)),',',num2str(T(i)),');']);
end
end
Region5_error=abs((R5-IF97)./IF97)
err=err+sum(sum(Region5_error>1E-8))
%T5_ph (Iteration)
p=[0.5,8,8];
h=[5219.76331549428,5206.09634477373,6583.80290533381];
IF97=[1500,1500,2000];
R5=zeros(1,3);
for i=1:3
R5(i)=T5_ph(p(i),h(i));
end
T5_ph_error=abs((R5-IF97)./IF97)
err=err+sum(sum(T5_ph_error>1E-7)) %Decimals in IF97
%T5_ps (Iteration)
p=[0.5,8,8];
s=[9.65408430982588,8.36546724495503,9.15671044273249];
IF97=[1500,1500,2000];
R5=zeros(1,3);
for i=1:3
R5(i)=T5_ps(p(i),s(i));
end
T5_ps_error=abs((R5-IF97)./IF97)
err=err+sum(sum(T5_ps_error>1E-4)) %Decimals in IF97
%*********************************************************************************************************
%* 7.6 Verifiy calling funtions
fun={'Tsat_p','T_ph','T_ps','T_hs','psat_T','p_hs','hV_p','hL_p','hV_T','hL_T','h_pT','h_ps','h_px','h_prho','h_Tx','vV_p','vL_p','vV_T','vL_T','v_pT','v_ph','v_ps','rhoV_p','rhoL_p','rhoV_T','rhoL_T','rho_pT','rho_ph','rho_ps','sV_p','sL_p','sV_T','sL_T','s_pT','s_ph','uV_p','uL_p','uV_T','uL_T','u_pT','u_ph','u_ps','CpV_p','CpL_p','CpV_T','CpL_T','Cp_pT','Cp_ph','Cp_ps','CvV_p','CvL_p','CvV_T','CvL_T','Cv_pT','Cv_ph','Cv_ps','wV_p','wL_p','wV_T','wL_T','w_pT','w_ph','w_ps','my_pT','my_ph','my_ps','tcL_p','tcV_p','tcL_T','tcV_T','tc_pT','tc_ph','tc_hs','st_T','st_p','x_ph','x_ps','vx_ph','vx_ps'};
In1={'1','1','1','100','100','84','1','1','100','100','1','1','1','1','100','1','1','100','100','1','1','1','1','1','100','100','1','1','1','0.006117','0.0061171','0.0001','100','1','1','1','1','100','100','1','1','1','1','1','100','100','1','1','1','1','1','100','100','1','1','1','1','1','100','100','1','1','1','1','1','1','1','1','25','25','1','1','100','100','1','1','1','1','1'};
In2={'0','100','1','0.2','0','0.296','0','0','0','0','20','1','0.5','2','0.5','0','0','0','0','100','1000','5','0','0','0','0','100','1000','1','0','0','0','0','20','84.01181117','0','0','0','0','100','1000','1','0','0','0','0','100','200','1','0','0','0','0','100','200','1','0','0','0','0','100','200','1','100','100','1','0','0','0','0','25','100','0.34','0','0','1000','4','418','4'};
Control={'99.60591861','23.84481908','73.70859421','13.84933511','1.014179779','2.295498269','2674.949641','417.4364858','2675.572029','419.099155','84.01181117','308.6107171','1546.193063','1082.773391','1547.33559210927','1.694022523','0.001043148','1.671860601','0.001043455','1.695959407','0.437925658','1.03463539','0.590310924','958.6368897','0.598135993','958.3542773','0.589636754','2.283492601','975.6236788','9.155465556','1.8359E-05','9.155756716','1.307014328','0.296482921','0.296813845','2505.547389','417.332171','2506.015308','418.9933299','2506.171426','956.2074342','308.5082185','2.075938025','4.216149431','2.077491868','4.216645119','2.074108555','4.17913573168802','4.190607038','1.552696979','3.769699683','1.553698696','3.76770022','1.551397249','4.035176364','3.902919468','472.0541571','1545.451948','472.2559492','1545.092249','472.3375235','1542.682475','1557.8585','1.22704E-05','0.000914003770302108','0.000384222','0.677593822','0.024753668','0.607458162','0.018326723','0.607509806','0.605710062','0.606283124','0.0589118685876641','0.058987784','0.258055424','0.445397961','0.288493093','0.999233827'};
for i=1:length(fun)
T=['XSteam(''',char(fun(i)),''',',char(In1(i)),',',char(In2(i)),')'];
Res=eval(T);
Error=(Res-(str2num(char(Control(i)))))/str2num(char(Control(i)));
Check=[T,'=',num2str(Res),' - (Control)',char(Control(i)),'=',num2str(Error)]
if Error>1E-5
err=err+1;
error('To large error')
end
end
Step7: The pump work is calculated from the steady flow energy equation using the enthalpy change from state 2 and state 3.
Step8: The heat added in boiler and heat rejected in condenser are also calculated from the steady flow energy equation using the enthalpy change during process (2-3) and process (4-1) respectively.
Step9: The turbine output work is calculated from the steady flow energy equation using the enthalpy change from state 3 and state 4.
Step10: The net output work is calculated as the difference between turbine work and pump work.
Step11: The thermal efficiency is calculated as the ratio of net work output to heat added in boiler.
Step12: The back work ratio is calculated as the ratio of pump work to turbine work.
IV Errors:
- While computing the vapor dome values using the for loop command an error is prompted as shown in the below figure:
The part of the code with error is:
t_s = linspace(5,400,500);
for i =1:length(t_s)
s(i) = XSteam('sL_T',t_s);
S(i) = XSteam('sV_T',t_s(i));
h(i) = XSteam('hL_T',t_s(i));
H(i) = XSteam('hV_T',t_s(i));
end
It is due to syntax error in indexing the for loop. It is rectified by properly indexing the variables in for loop.
t_s = linspace(5,400,500);
for i =1:length(t_s)
s(i) = XSteam('sL_T',t_s(i));
S(i) = XSteam('sV_T',t_s(i));
h(i) = XSteam('hL_T',t_s(i));
H(i) = XSteam('hV_T',t_s(i));
end
V. Output:
The outputs from this code for simulating ideal Rankine cycle is the properties at different states in the cycle and energy conversion from heat added in boiler to work generated in the turbine. Also, the pump work required, the thermal efficiency of the cycle, the back work ratio for the given input parameters can also be calculated using this simulator. The details of the outputs from the simulator code are printed in the command window and are given in the below figures:
The thermal efficiency of the ideal Rankine cycle is about 36% for the given input parameters.
The output T-s diagram and h-s diagram for the cycle are shown in the below figures:
The T-s diagram shown above represents the temperature and entropy changes during the cycle. The blue line represents the isobaric heating in the boiler(conversion of water into steam) from state2 to state3. The green line represents the isentropic expansion in the turbine. The pink line represents the constant pressure heat rejection in the condenser. The zoomed view on the left of the figure is the isentropic compression(red line) during the pumping of water from condenser to boiler pressure shown for better clarity. The black curve is the saturated liquid vapour-dome.
Similarly, the h-s diagram is plotted and shown in the below figure:
VI. Conclusion:
In this project, a program is executed in MATLAB for simulating the ideal Rankine cycle for the given input parameters. A function code containing a steam table library is used for calculating the various state properties, turbine work and thermal efficiency of the cycle. From the input taken from the user end and the properties calculated T-s and h-s diagrams are plotted for the Ideal Rankine cycle.
References:
1.https://www.instructables.com/Simple-Rankine-Cycle-Thermodynamics/
2. ENGINEERING THERMODYNAMICS, R.K. Rajput
3. https://www.thermopedia.com/content/1072/
4. https://in.mathworks.com/help/matlab/ref/for.html
5. Thermodynamics: An Engineering Approach, YUNUS. TURNER CENGEL (ROBERT. CIMBALA, JOHN.), Michael A. Boles