function OutVar = H2Oprop(query, Var1, Var2, Var3) % % H2Oprop.m v1.01 % % Developed by Bryan Ruddy (ruddy AT mit.edu), 9/2006 - 5/2009 % % The most recent version of this function can be found at http://web.mit.edu/ruddman/Public/Thermodynamics/H2Oprop.m % % This work is licensed under the Creative Commons Attribution 3.0 United % States License. To view a copy of this license, visit % http://creativecommons.org/licenses/by/3.0/us/ or send a letter to Creative % Commons, 171 Second Street, Suite 300, San Francisco, California, 94105, % USA. % % This script computes the thermodynamic and transport properties of water % using the latest (December 2008) IAPWS scientific formulations, with % rigorous bounds-checking. Information on the range of validity and % accuracy of these properties can be found in the IAPWS releases at % www.iapws.org. A total of 70 different combinations of input and output % parameters are supported, including 45 basic thermodynamic relationships, % 14 saturation relationships, and 11 other thermodynamic and transport properties. % % All quantities are represented in SI units. % % Usage: % The function should be called with a case-insensitive text string (query) % and up to three scalar numerical parameters. The identities of these % parameters depend on the choice of query, as described below. The % function outputs a single scalar value, whose meaning also depends on the % choice of query. Note that queries in terms of density and temperature % are the most computationally efficient; others use root-finding and/or % minimization routines. % % Query Options: % % t_dp Temperature from density and pressure (gives high-temperature point near density anomaly) % % t_dh Temperature from density and specific enthalpy % % t_ph Temperature from pressure and specific enthalpy % % t_du Temperature from density and specific internal energy % % t_pu Temperature from pressure and specific internal energy % % t_ds Temperature from density and specific entropy % % t_ps Temperature from pressure and specific entropy % % t_hs Temperature from specific enthalpy and specific entropy % % t_us Temperature from specific internal energy and specific entropy % % p_dt Pressure from density and temperature % % p_dh Pressure from density and specific enthalpy % % p_du Pressure from density and specific internal energy % % p_ds Pressure from density and specific entropy % % p_ts Pressure from temperature and specific entropy % % p_hs Pressure from specific enthalpy and specific entropy % % p_us Pressure from specific internal energy and specific entropy % % d_pt Density from pressure and temperature (gives liquid value if saturated) % % d_ph Density from pressure and specific enthalpy % % d_pu Density from pressure and specific internal energy % % d_ps Density from pressure and specific entropy % % d_ts Density from temperature and specific entropy % % d_hs Density from specific enthalpy and specific entropy % % d_us Density from specific internal energy and specific entropy % % h_dt Specific enthalpy from density and temperature % % h_pt Specific enthalpy from pressure and temperature (gives liquid value if saturated) % % h_du Specific enthalpy from density and specific internal energy % % h_ds Specific enthalpy from density and specific entropy % % h_pu Specific enthalpy from pressure and specific internal energy % % h_ps Specific enthalpy from pressure and specific entropy % % h_ts Specific enthalpy from temperature and specific entropy % % h_us Specific enthalpy from specific internal energy and specific entropy % % u_dt Specific internal energy from density and temperature % % u_dh Specific internal energy from density and specific enthalpy % % u_ds Specific internal energy from density and specific entropy % % u_pt Specific internal energy from pressure and temperature (gives liquid value if saturated) % % u_ph Specific internal energy from pressure and specific enthalpy % % u_ps Specific internal energy from pressure and specific entropy % % u_ts Specific internal energy from temperature and specific entropy % % u_hs Specific internal energy from specific enthalpy and specific entropy % % s_dt Specific entropy from density and temperature % % s_dh Specific entropy from density and specific enthalpy % % s_du Specific entropy from density and specific internal energy % % s_pt Specific entropy from pressure and temperature (gives liquid value if saturated) % % s_ph Specific entropy from pressure and specific enthalpy % % s_pu Specific entropy from pressure and specific internal energy % % psat_t Saturation pressure from temperature % % tsat_p Saturation temperature from pressure % % tsat_dl Saturation temperature from liquid density (gives high-temperature point near density anomaly) % % tsat_dv Saturation temperature from vapor density % % dlsat_t Saturated liquid density from temperature % % dlsat_p Saturated liquid density from pressure % % dvsat_t Saturated vapor density from temperature % % dvsat_p Saturated vapor density from pressure % % hlsat_t Saturated liquid specific enthalpy from temperature % % hlsat_p Saturated liquid specific enthalpy from pressure % % hvsat_t Saturated vapor specific enthalpy from temperature % % hvsat_p Saturated vapor specific enthalpy from pressure % % ulsat_t Saturated liquid specific internal energy from temperature % % ulsat_p Saturated liquid specific internal energy from pressure % % uvsat_t Saturated vapor specific internal energy from temperature % % uvsat_p Saturated vapor specific internal energy from pressure % % slsat_t Saturated liquid specific entropy from temperature % % slsat_p Saturated liquid specific entropy from pressure % % svsat_t Saturated vapor specific entropy from temperature % % svsat_p Saturated vapor specific entropy from pressure % % x_dt Steam quality from density and temperature % % cp_dt Constant pressure specific heat from density and temperature % % cv_dt Constant volume specific heat from density and temperature % % betat_dt Thermal expansion coefficient from density and temperature % % kappat_dt Isothermal compressibility from density and temperature % % w_dt Sound velocity from density and temperature % % mu_dt Viscosity from density and temperature % % k_dt Thermal conductivity from density and temperature % % sigma_t Surface tension from temperature % % epsilon_dt Dielectric constant from density and temperature % % n_dtl Refractive index from density, temperature, and wavelength % Programming Notes: % % Code known to be problematic, for efficiency and/or accuracy reasons, is % marked with a *** in the comments. This code should be fixed in future % versions. % % In general, a major objective in future revisions should be to reduce the % degree to which the code is recursive. This will increase performance % and reliability by eliminating unneccessary bounds checks and reducing % the use of nested fzero calls. % % Change Log: % % v1.01: Replaced recursive function references with references to this % function's current name, rather than the development name. %--------------------------------------------------------------------------------------- % Equation of State Constants % 1995 IAPWS Scientific Formulation global Tc Tt Pc rhoc R n0 gamma0 c d t a b B C D A alpha beta gamma epsilon n Tc = 647.096; %Critical Temperature, K rhoc = 322; %Critical Density, kg/m^3 R = 461.51805; %Water Gas Constant, J/kg-K n0 = [-8.32044648374976; 6.68321052759325; 3.00632; 0.012436; 0.97315; 1.27950; ... 0.96956; 0.24873]; gamma0 = [0; 0; 0; 1.28728967; 3.53734222; 7.74073708; 9.24437796; 27.5075105]; c = [zeros(7,1); ones(15,1); 2*ones(20,1); 3; 3; 3; 3; 4; 6; 6; 6; 6; zeros(5,1)]; d = [1; 1; 1; 2; 2; 3; 4; 1; 1; 1; 2; 2; 3; 4; 4; 5; 7; 9; 10; 11; 13; 15; 1; 2;... 2; 2; 3; 4; 4; 4; 5; 6; 6; 7; 9*ones(5,1); 10; 10; 12; 3; 4; 4; 5; 14; 3; 6;... 6; 6; 3; 3; 3; 0; 0]; t = [-0.5; 0.875; 1; 0.5; 0.75; 0.375; 1; 4; 6; 12; 1; 5; 4; 2; 13; 9; 3; 4; 11;... 4; 13; 1; 7; 1; 9; 10; 10; 3; 7; 10; 10; 6; 10; 10; 1; 2; 3; 4; 8; 6; 9; 8;... 16; 22; 23; 23; 10; 50; 44; 46; 50; 0; 1; 4; 0; 0]; a = [zeros(54,1); 3.5; 3.5]; b = [zeros(54,1); 0.85; 0.95]; B = [zeros(54,1); 0.2; 0.2]; C = [zeros(54,1); 28; 32]; D = [zeros(54,1); 700; 800]; A = [zeros(54,1); 0.32; 0.32]; alpha = [zeros(51,1); 20; 20; 20; 0; 0]; beta = [ones(51,1); 150; 150; 250; 0.3; 0.3]; gamma = [zeros(51,1); 1.21; 1.21; 1.25; 0; 0]; epsilon = [zeros(51,1); 1; 1; 1; 0; 0]; n = ... [0.12533547935523e-1 0.78957634722828e1 -0.87803203303561e1 0.31802509345418 ... -0.26145533859358 -0.78199751687981e-2 0.88089493102134e-2 -0.66856572307965 ... 0.20433810950965 -0.66212605039687e-4 -0.19232721156002 -0.25709043003438 ... 0.16074868486251 -0.40092828925807e-1 0.39343422603254e-6 -0.75941377088144e-5 ... 0.56250979351888e-3 -0.15608652257135e-4 0.11537996422951e-8 0.36582165144204e-6 ... -0.13251180074668e-11 -0.62639586912454e-9 -0.10793600908932 0.17611491008752e-1 ... 0.22132295167546 -0.40247669763528 0.58083399985759 0.49969146990806e-2 ... -0.31358700712549e-1 -0.74315929710341 0.47807329915480 0.20527940895948e-1 ... -0.13636435110343 0.14180634400617e-1 0.83326504880713e-2 -0.29052336009585e-1 ... 0.38615085574206e-1 -0.20393486513704e-1 -0.16554050063734e-2 0.19955571979541e-2 ... 0.15870308324157e-3 -0.16388568342530e-4 0.43613615723811e-1 0.34994005463765e-1 ... -0.76788197844621e-1 0.22446277332006e-1 -0.62689710414685e-4 -0.55711118565645e-9 ... -0.19905718354408 0.31777497330738 -0.11841182425981 -0.31306260323435e2 ... 0.31546140237781e2 -0.25213154341695e4 -0.14874640856724 0.31806110878444]'; %--------------------------------------------------------------------------------------- % Phase Boundary Constants % Used in a variety of bounds-checking, root-finding, and minimization % tasks throughout the code. Pc = 22.064e6; % Critical Pressure (IAPWS95) Pt = 611.654771007894; % Triple Point Pressure (IAPWS95) %Pmax = 1e9; % Maximum Pressure for Experimental Verification (IAPWS95) Picemax = 20617812820.4493; % Maximum Pressure on Melting Curve (715 K) (IAPWS08) Pmin = 1.9349587189864e-040; % Minimum Pressure on Sublimation curve (50 K) (IAPWS08) Tt = 273.16; % Triple Point Temperature (IAPWS95) Tmax = 1273; % Maximum Temperature for Experimental Verification (IAPWS95) Tmin = 251.165; % Minimum Liquid Temperature (Ice-Ice-Liquid Triple Point) (IAPWS08) Trhomax = 277.148195578572; % Highest Temperature of Density Maximum Anomaly Tabsmin = 50; % Minimum Temperature for Experimental Verification (IAPWS08) Ticemax = 715; % Maximum Freezing Temperature (IAPWS08) Tsmin = 694.274446246169; % Temperature of Local Minimum Freezing Entropy Tsmax = 603.867446392581; % Temperature of Maximum Freezing Entropy rhovt = 0.00485457572477878; % Triple Point Vapor Density rholt = 999.792520031622; % Triple Point Liquid Density rhoamax = 999.925131084086; % Saturation Density Anomaly Point rhomax = 1237.39112310848; % Maximum Density (1GPa on melting curve) rholTmin = 1091.68582714549; % Liquid Density at Tmin rhomin = 3.29347680703935e-46; % Minimum Density (Pmin at Tmax) rhoicemax = 1930.89627377429; % Maximum Known Liquid Density rhosatmin = 9.86827543571457e-045; % Minimum Phase Boundary Density (50 K) hlt = 0.611781702320426; % Saturated Liquid Enthalpy at Triple Point hvt = 2500915.19146571; % Saturated Vapor Enthalpy at Triple Point ha = 16805.6061030613; % Saturated Liquid Enthalpy at Density Anomaly hvsatmin = 2088366.77643276; % Saturated Vapor Enthalpy at 50 K uvsatmin = 2065290.87393276; % Saturated Vapor Energy at 50 K ult = -0.000236902622937808; % Saturated Liquid Energy at Triple Point uvt = 2374919.6757225; % Saturated Vapor Energy at Triple Point ua = 16804.792597053; % Saturated Liquid Energy at Density Anomaly svsatmin = 51104.228263276876532; % Saturated Vapor Entropy at 50 K slt = 2.72849374306335e-012; % Saturated Liquid Entropy at Triple Point svt = 9155.49340929858; % Saturated Vapor Entropy at Triple Point sa = 61.0753924234768; % Saturated Liquid Entropy at Density Anomaly spimax = 1300.31975929576; % Highest Liquid Entropy on Pressurized Ice Line spitmax = 1246.430489334724371; % Liquid Entropy at High End of Pressurized Ice Line spimin = -337.232036172874; % Lowest Liquid Entropy (at Tmin) spitmin = 1213.949641489883561; % Liquid Entropy Minimum near Top of Hot Ice Line scrit = 4407.116565421005362; % Entropy at Critical Point sdtmax = 3267.889828921342541; % Entropy at Tmax and rhoicemax %--------------------------------------------------------------------------------------- % Function Body query = lower(query); % Strip out case from input parameter switch query case 'p_dt' % Pressure from density, temperature rho = Var1; % Input density T = Var2; % Input temperature if (T > Tc) % Supercritical P = Water_P(rho, T); % Compute directly from correlation elseif (T >= Tt) % Liquid - Vapor Regime [Psat rholsat rhovsat] = GoodSatProp(T); % Find accurate saturation properties if (rho > rhovsat) && (rho < rholsat) % Two-Phase P = Psat; % Clearly the fluid is saturated... else % Liquid or Vapor P = Water_P(rho, T); % We aren't saturated, so directly compute pressure if P > MeltPress(T) % High-Pressure Ice P = NaN; % There's no correlation for this stuff... end end elseif (T >= Tmin) % Ice - Liquid - Vapor Regime [Pm Ps Pf] = MeltPress(T); % Find melting, freezing, and sublimation pressures rhovsub = fzero(@(x) Water_P(x, T) - Ps, Ps/R/T); % Find vapor density at sublimation point rholmin = fzero(@(x) Water_P(x, T) - Pm, rhomax); % Find liquid density at ice Ih melting point rholmax = fzero(@(x) Water_P(x, T) - Pf, rhomax); % Find liquid density at high-pressure ice melting point if (rho <= rhovsub) || ((rho >= rholmin) && (rho <= rholmax)) %Liquid or Gas P = Water_P(rho, T); % Single-phase, so directly compute pressure else % Ice P = NaN; % Would need the ice correlation to compute this... end else % Ice - Vapor Regime Ps = MeltPress(T); % Find sublimation pressure rhovsat = fzero(@(x) Water_P(x, T) - Ps, Ps/R/T); % Find saturated vapor density if (rho <= rhovsat) % Vapor P = Water_P(rho, T); % Single-phase, so directly compute pressure else % Ice P = NaN; % Would need the ice correlation to compute this... end end OutVar = P; % Output the pressure %--------------------------------------------------------------------------------------- case 'd_pt' % Density from pressure, temperature P = Var1; % Input pressure T = Var2; % Input temperature Tm = MeltTemp(P); % Compute solid-fluid phase transition temperature if (T < (Tm - 2*eps(Tm))) || isnan(Tm) % Ice or Absurd Pressure rho = NaN; % Would need the ice correlation to compute this... elseif T > Tc % Supercritical rho = fzero(@(x) Water_P(x, T) - P, [rhomin*.999 rhoicemax*1.01]); % Use root-finding method to compute density elseif T > Tt % Liquid - Vapor Regime [Psat rholsat rhovsat] = GoodSatProp(T); % Find accurate saturation properties if (P >= Psat) % Liquid rho = fzero(@(x) Water_P(x, T) - P, [rholsat*.999 rhoicemax]); % Use root-finding method, and give liquid density if saturated else % Vapor rho = fzero(@(x) Water_P(x, T) - P, [rhomin rhovsat]); % Use root-finding method to compute density end elseif T >= Tabsmin % Liquid - Ice - Vapor Regime if (P < Pt) % Thin Vapor rho = fzero(@(x) Water_P(exp(x), T) - P, log(P/R/T)); % Use root-finding method to compute density rho = exp(rho); % Due to the extreme pressure range, examine the logarithm of density else % Cold Liquid rho = fzero(@(x) Water_P(x, T) - P, rhomax); % Use root-finding method to compute density end else % Too Cold rho = NaN; % Vapor is effectively nonexistent below Tabsmin... end OutVar = rho; % Output density %--------------------------------------------------------------------------------------- case 't_dp' % Temperature from density, pressure, only valid above density anomaly rho = Var1; % Input density P = Var2; % Input pressure Tm = MeltTemp(P); % Melting temperature if (rho >= rholt) && (P < 3e7) % Density anomaly region [Trhomax1, rhomax1] = fminbnd(@(x) 1/H2Oprop('d_pt',P,x), ... Tm, Trhomax*1.001, optimset('TolX',1e-15)); % Find density anomaly rhomax1 = 1/rhomax1; % Correct the maximum density if (rho > rhomax1) % Non-Physical T = NaN; % At a given pressure, ice at any temperature floats... else % Liquid Tsat = H2Oprop('tsat_p',P); % Find the saturation temperature T = fzero(@(x) Water_P(rho,x) - P, [Trhomax1 Tsat]); % Use root-finding method to compute temperature end else % No density anomaly, skip the minimization Tsat = H2Oprop('tsat_p',P); % Find the saturation temperature from pressure [junk rhol rhov] = GoodSatProp(Tsat); % Find the saturation densities rhom = H2Oprop('d_pt',P,Tm); if (rho >= rhom) % Ice or Non-Physical T = NaN; elseif (rho >= rhol) % Liquid Tmin1 = H2Oprop('tsat_dl',rho); % Find the saturation temperature from density if isnan(Tmin1) % Check to see if that density is actually on the saturation line Tmin1 = Tm; % If not, use the melting temperature as the lower limit instead end T = fzero(@(x) Water_P(rho,x) - P, [Tmin1 Tsat]); % Use root-finding method to compute temperature elseif (rho <= rhov) % Vapor T = fzero(@(x) Water_P(rho,x) - P, [Tsat Tmax]); % Use root-finding method to compute temperature elseif (P > Pc) % Supercritical fluid if rho > rhoc % Liquid-like fluid Tmin1 = H2Oprop('tsat_dl',rho); % Find the saturation temperature from density else % Vapor-like fluid Tmin1 = H2Oprop('tsat_dv',rho); % Find the saturation temperature from density end if isnan(Tmin1) % Find the saturation temperature from density Tmin1 = Tm; % If not, use the melting temperature as the lower limit instead end T = fzero(@(x) Water_P(rho,x) - P, [Tmin1 Tmax]); % Use root-finding method to compute temperature else % Two-phase fluid T = Tsat; % Fluid is at saturation temperature end end OutVar = T; % Output temperature %--------------------------------------------------------------------------------------- case 'u_dt' % Specific energy from density, temperature rho = Var1; % Input density T = Var2; % Input temperature if (T < Tmin) % Ice or Vapor Pm = MeltPress(T); % Find the sublimation pressure from temperature rhos = H2Oprop('d_pt',Pm,T); % Compute the vapor density at the sublimation point if (rho > rhos) % Ice u = NaN; % Would need the ice correlation to compute this... else % Vapor u = Water_u(rho, T); % Directly compute the specific energy end elseif (T < Tt) % Ice, Liquid, or Vapor [Pm Ps Pf] = MeltPress(T); % Find all phase boundary pressures at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the Ice Ih melting point rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density at the high-pressure ice melting point if (rho < rhos) % Vapor u = Water_u(rho, T); % Directly compute the specific energy elseif (rho < rhom) % Ice Ih u = NaN; % Would need the ice correlation to compute this... elseif (rho < rhof) % Liquid u = Water_u(rho, T); % Directly compute the specific energy else % High-Pressure Ice u = NaN; % There's no correlation for this stuff... end else % Liquid or Vapor [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature if (rho > rhovsat) && (rho < rholsat) % Two-phase Mixture v = 1/rho; % Convert density to specific volume vl = 1/rholsat; % Convert saturated liquid density to specific volume vv = 1/rhovsat; % Convert saturated vapor density to specific volume X = (v - vl)/(vv - vl); % Compute the quality from specific volume ul = Water_u(rholsat, T); % Directly compute the saturated liquid specific energy uv = Water_u(rhovsat, T); % Directly compute the saturated vapor specific energy u = X*uv + (1-X)*ul; % Find the mixture specific energy else %Single-phase u = Water_u(rho, T); % Directly compute the specific energy end end OutVar = u; % Output specific internal energy %--------------------------------------------------------------------------------------- case 's_dt' % Specific entropy from density, temperature rho = Var1; % Input density T = Var2; % Input temperature if (T < Tmin) % Ice or Vapor rhos = H2Oprop('dvsat_t',T); % Compute the vapor density at the sublimation point if (rho > rhos) % Ice s = NaN; % Would need the ice correlation to compute this... else % Vapor s = Water_s(rho, T); % Directly compute the specific entropy end elseif (T < Tt) % Ice, Liquid, or Vapor [Pm Ps Pf] = MeltPress(T); % Find all phase boundary pressures at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the Ice Ih melting point rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density at the high-pressure ice melting point if (rho < (rhos+10*eps(rhos))) % Vapor s = Water_s(rho, T); % Directly compute the specific entropy elseif (rho < rhom) % Ice Ih s = NaN; % Would need the ice correlation to compute this... elseif (rho <= rhof) % Liquid s = Water_s(rho, T); % Directly compute the specific entropy else % High-Pressure Ice s = NaN; % There's no correlation for this stuff... end else % Liquid or Vapor [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature if (rho > rhovsat) && (rho < rholsat) % Two-phase Mixture v = 1/rho; % Convert density to specific volume vl = 1/rholsat; % Convert saturated liquid density to specific volume vv = 1/rhovsat; % Convert saturated vapor density to specific volume X = (v - vl)/(vv - vl); % Compute the quality from specific volume sl = Water_s(rholsat, T); % Directly compute the saturated liquid specific entropy sv = Water_s(rhovsat, T); % Directly compute the saturated vapor specific entropy s = X*sv + (1-X)*sl; % Find the mixture specific entropy else %Single-phase s = Water_s(rho, T); % Directly compute the specific entropy end end OutVar = s; % Output specific entropy %--------------------------------------------------------------------------------------- case 'h_dt' % Specific enthalpy from density, temperature rho = Var1; % Input density T = Var2; % Input temperature if (T < Tmin) % Ice or Vapor Pm = MeltPress(T); % Find the sublimation pressure from temperature rhos = H2Oprop('d_pt',Pm,T); % Compute the vapor density at the sublimation point if (rho > (rhos+10*eps(rhos))) % Ice h = NaN; % Would need the ice correlation to compute this... else % Vapor h = Water_h(rho, T); % Directly compute the specific enthalpy end elseif (T < Tt) % Ice, Liquid, or Vapor [Pm Ps Pf] = MeltPress(T); % Find all phase boundary pressures at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the Ice Ih melting point rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density at the high-pressure ice melting point if (rho < rhos) % Vapor h = Water_h(rho, T); % Directly compute the specific enthalpy elseif (rho < rhom) % Ice Ih h = NaN; % Would need the ice correlation to compute this... elseif (rho <= rhof) % Liquid h = Water_h(rho, T); % Directly compute the specific enthalpy else % High-Pressure Ice h = NaN; % There's no correlation for this stuff... end else % Liquid or Vapor [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature if (rho > rhovsat) && (rho < rholsat) % Two-phase Mixture v = 1/rho; % Convert density to specific volume vl = 1/rholsat; % Convert saturated liquid density to specific volume vv = 1/rhovsat; % Convert saturated vapor density to specific volume X = (v - vl)/(vv - vl); % Compute the quality from specific volume hl = Water_h(rholsat, T); % Directly compute the saturated liquid specific enthalpy hv = Water_h(rhovsat, T); % Directly compute the saturated vapor specific enthalpy h = X*hv + (1-X)*hl; % Find the mixture specific enthalpy else %Single-phase h = Water_h(rho, T); % Directly compute the specific enthalpy end end OutVar = h; % Output specific enthalpy %--------------------------------------------------------------------------------------- case 't_dh' % Temperature from density, enthalpy rho = Var1; % Input density h = Var2; % Input specific enthalpy if (rho < rhomin) % Too rarefied T = NaN; % Densities this low violate continuum approximations anyway... elseif (rho < rhosatmin) % Extremely Thin Vapor (less dense than 50 K saturated vapor) if (h >= hvsatmin) % Vapor T = fzero(@(x) Water_h(rho, x) - h, [Tabsmin Tmax]); % Use root-finding method to compute temperature else % Non-Physical T = NaN; % It's not ice at this density... end elseif (rho < rhovt) % Thin Vapor (less dense than triple point vapor) Tsat = H2Oprop('tsat_dv',rho); % Find saturation temperature from density hsat = Water_h(rho, Tsat); % Find saturated vapor enthalpy at this density if (h >= hsat) % Vapor T = fzero(@(x) Water_h(rho, x) - h, [Tsat Tmax]); % Use root-finding method to compute temperature else % Two-Phase with Ice T = NaN; % Would need the ice correlation to compute this... end elseif (rho < rhoc) % Vapor or Two-Phase with Liquid Tsat = H2Oprop('tsat_dv',rho); % Find saturation temperature from density hsat = Water_h(rho, Tsat); % Find saturated vapor enthalpy at this density Xmin = H2Oprop('x_dt',rho,Tt); % Find quality for this density at the triple point hmin = Xmin*hvt + (1-Xmin)*hlt; % Find enthalpy for two-phase mixture at triple point if (h >= hsat) % Vapor T = fzero(@(x) Water_h(rho, x) - h, [Tsat Tmax]); % Use root-finding method to compute temperature elseif (h >= hmin) % Two-Phase with Liquid T = fzero(@(x) H2Oprop('h_dt',rho,x) - h, [Tt Tsat]); % Use root-finding method to compute temperature else % Multi-Phase with Ice T = NaN; % Would need the ice correlation to compute this... end elseif (rho < rholt) % Liquid or Two-Phase with Vapor (No Density Anomaly) Tsat = H2Oprop('tsat_dl',rho); % Find saturation temperature from density hsat = Water_h(rho, Tsat); % Find saturated liquid enthalpy at this density Xmin = H2Oprop('x_dt',rho,Tt); % Find quality for this density at the triple point hmin = Xmin*hvt + (1-Xmin)*hlt; % Find enthalpy for two-phase mixture at triple point if (h >= hsat) % Liquid T = fzero(@(x) Water_h(rho, x) - h, [Tsat Tmax]); % Use root-finding method to compute temperature elseif (h >= hmin) % Two-Phase with Vapor T = fzero(@(x) H2Oprop('h_dt',rho,x) - h, [Tt Tsat]); % Use root-finding method to compute temperature else % Ice T = NaN; % Would need the ice correlation to compute this... end elseif (rho < rhoamax) % Liquid or Two-Phase with Vapor (Density Anomaly) if h >= ha % Warm Liquid or Two-Phase with Vapor Tsat = H2Oprop('tsat_dl',rho); % Find saturation temperature from density hsat = Water_h(rho, Tsat); % Find saturated liquid enthalpy at this density if (h >= hsat) % Liquid T = fzero(@(x) Water_h(rho, x) - h, [Tsat Tmax]); % Use root-finding method to compute temperature else % Two-Phase with Vapor T = fzero(@(x) H2Oprop('h_dt',rho,x) - h, [Tt Tsat]); % Use root-finding method to compute temperature end else % Cold Liquid or Two-Phase with Vapor Tsat = fzero(@(x) GoodSatProp(x,'rholsat') - rho, [Tt Trhomax]);% Find saturation temperature from density hsat = Water_h(rho, Tsat); % Find saturated liquid enthalpy at this density Tm = fzero(@(x) Water_P(rho, x) - MeltPress(x), [Tmin Tt]); % Find melting temperature at this density hm = Water_h(rho, Tm); % Find enthalpy of liquid at melting point with this density if (h >= hsat) % Two-Phase with Vapor T = fzero(@(x) H2Oprop('h_dt',rho,x) - h, [Tsat Tmax]); % Use root-finding method to compute temperature elseif (h >= hm) % Liquid T = fzero(@(x) Water_h(rho, x) - h, [Tm Tsat]); % Use root-finding method to compute temperature else % Ice T = NaN; % Would need the ice correlation to compute this... end end elseif (rho < rholTmin) % Pressurized Liquid (up to density of liquid at liquid-Ice Ih-Ice III triple point) Tm = fzero(@(x) H2Oprop('d_pt',MeltPress(x),x) - rho, [Tmin Tt]);% Find melting temperature from density hm = Water_h(rho,Tm); % Find enthalpy of liquid at melting point with this density if (h >= hm) % Liquid T = fzero(@(x) Water_h(rho, x) - h, [Tm Tmax]); % Use root-finding method to compute temperature else % Ice T = NaN; % Would need the ice correlation to compute this... end elseif (rho < rhoicemax) % Very Pressurized Liquid (Over 200 MPa) Tm = fzero(@(x) H2Oprop('d_pt',MeltPress(x,'hp'),x) - rho, ... [Tmin Ticemax]); % Find melting temperature from density hm = Water_h(rho,Tm); % Find enthalpy of liquid at melting point with this density if (h >= hm) % Liquid T = fzero(@(x) Water_h(rho, x) - h, [Tm Tmax]); % Use root-finding method to compute temperature else % High-Pressure Ice T = NaN; % There's no correlation for this stuff... end else % Too Dense (High-Pressure Ice or Over 20 GPa) T = NaN; % There's no correlation for this stuff... end OutVar = T; % Output Temperature %--------------------------------------------------------------------------------------- case 'd_ph' % Density from pressure, enthalpy P = Var1; % Input pressure h = Var2; % Input specific enthalpy if (P < Pmin) % Too rarefied rho = NaN; % Pressures this low likely violate continuum approximations... elseif (P < Pt) % Thin vapor (below triple point pressure) Tsat = MeltTemp(P); % Compute the sublimation temperature at this pressure rhosat = H2Oprop('d_pt',P,Tsat); % Compute the density of saturated vapor at this pressure hsat = Water_h(rhosat, Tsat); % Compute the enthalpy of saturated vapor at this pressure if (h >= hsat) % Vapor T = fzero(@(x) h - H2Oprop('h_pt',P,x), [Tsat Tmax]); % Use root-finding method to compute temperature rho = H2Oprop('d_pt',P,T); % Use the temperature and given pressure to find density else % Ice rho = NaN; % Would need the ice correlation to compute this... end elseif (P < Pc) % Liquid or vapor Tm = MeltTemp(P); % Compute the melting temperature at this pressure rhom = H2Oprop('d_pt',P,Tm); % Compute the density of freezing liquid at this pressure hm = Water_h(rhom,Tm); % Compute the enthalpy of freezing liquid at this pressure Tsat = H2Oprop('tsat_p',P); % Compute the saturation temperature at this pressure [junk rholsat rhovsat] = GoodSatProp(Tsat); % Compute the saturated liquid and vapor densities at this pressure hlsat = Water_h(rholsat,Tsat); % Compute the saturated liquid enthalpy at this pressure hvsat = Water_h(rhovsat,Tsat); % Compute the saturated vapor enthalpy at this pressure if (h <= hm) % Ice rho = NaN; % Would need the ice correlation to compute this... elseif (h <= hlsat) % Liquid T = fzero(@(x) h - H2Oprop('h_pt',P,x), [Tm Tsat]); % Use root-finding method to compute temperature rho = H2Oprop('d_pt',P,T); % Use the temperature and given pressure to find density elseif (h < hvsat) % Two-Phase rho = fzero(@(x) H2Oprop('h_dt',x,Tsat) - h, [rhovsat rholsat]); % Use root-finding method to compute density else % Vapor T = fzero(@(x) h - H2Oprop('h_pt',P,x), [Tsat Tmax]); % Use root-finding method to compute temperature rho = H2Oprop('d_pt',P,T); % Use the temperature and given pressure to find density end elseif (P < Picemax) % Supercritical fluid Tm = MeltTemp(P); % Compute the melting temperature at this pressure rhom = H2Oprop('d_pt',P,Tm); % Compute the density of freezing fluid at this pressure hm = Water_h(rhom,Tm); % Compute the enthalpy of freezing fluid at this pressure if (h >= hm) % Fluid T = fzero(@(x) h - H2Oprop('h_pt',P,x), [Tm Tmax]); % Use root-finding method to compute temperature rho = H2Oprop('d_pt',P,T); % Use the temperature and given pressure to find density else % Ice rho = NaN; % There's no correlation for this stuff... end else % Too pressurized (cannot verify fluidity) rho = NaN; % There's no correlation for this stuff... end OutVar = rho; % Output density %--------------------------------------------------------------------------------------- case 't_ph' % Temperature from pressure, enthalpy P = Var1; % Input pressure h = Var2; % Input specific enthalpy rho = H2Oprop('d_ph',P,h); % Compute density from given pressure and enthalpy if isnan(rho) % Do the given pressure and enthalpy lie in a valid region? T = NaN; % If not, propagate the NaN else % If so, we can compute temperature T = H2Oprop('t_dh',rho,h); % Compute temperature from density and enthalpy end OutVar = T; % Output temperature %--------------------------------------------------------------------------------------- case 't_du' % Temperature from density, energy rho = Var1; % Input density u = Var2; % Input specific internal energy if (rho < rhomin) % Too rarefied T = NaN; % Densities this low likely violate continuum approximations... elseif (rho < rhosatmin) % Extremely Thin Vapor (less dense than 50 K saturated vapor) if (u >= uvsatmin) % Vapor T = fzero(@(x) Water_u(rho, x) - u, [Tabsmin Tmax]); % Use root-finding method to compute temperature else % Non-Physical T = NaN; % It's not ice at this density... end elseif (rho < rhovt) % Thin Vapor (less dense than triple point vapor) Tsat = H2Oprop('tsat_dv',rho); % Find saturation temperature from density usat = Water_u(rho, Tsat); % Find saturated vapor energy at this density if (u >= usat) % Vapor T = fzero(@(x) Water_u(rho, x) - u, [Tsat Tmax]); % Use root-finding method to compute temperature else % Ice T = NaN; % Would need the ice correlation to compute this... end elseif (rho < rhoc) % Vapor or Two-Phase with Liquid Tsat = H2Oprop('tsat_dv',rho); % Find saturation temperature from density usat = Water_u(rho, Tsat); % Find saturated vapor energy at this density Xmin = H2Oprop('x_dt',rho,Tt); % Find quality for this density at the triple point umin = Xmin*uvt + (1-Xmin)*ult; % Find energy for two-phase mixture at triple point if (u >= usat) % Vapor T = fzero(@(x) Water_u(rho, x) - u, [Tsat Tmax]); % Use root-finding method to compute temperature elseif (u >= umin) % Two-Phase with Liquid T = fzero(@(x) H2Oprop('u_dt',rho,x) - u, [Tt Tsat]); % Use root-finding method to compute temperature else % Multi-Phase with Ice T = NaN; % Would need the ice correlation to compute this... end elseif (rho < rholt) % Liquid or Two-Phase with Vapor (No Density Anomaly) Tsat = H2Oprop('tsat_dl',rho); % Find saturation temperature from density usat = Water_u(rho, Tsat); % Find saturated liquid energy at this density Xmin = H2Oprop('x_dt',rho,Tt); % Find quality for this density at the triple point umin = Xmin*uvt + (1-Xmin)*ult; % Find energy for two-phase mixture at triple point if (u >= usat) % Liquid T = fzero(@(x) Water_u(rho, x) - u, [Tsat Tmax]); % Use root-finding method to compute temperature elseif (u >= umin) % Two-Phase with Vapor T = fzero(@(x) H2Oprop('u_dt',rho,x) - u, [Tt Tsat]); % Use root-finding method to compute temperature else % Ice T = NaN; % Would need the ice correlation to compute this... end elseif (rho < rhoamax) % Liquid or Two-Phase with Vapor (Density Anomaly) if u >= ua % Warm Liquid or Two-Phase with Vapor Tsat = H2Oprop('tsat_dl',rho); % Find saturation temperature from density (assume above Trhomax) usat = Water_u(rho, Tsat); % Find saturated liquid energy at this density above Trhomax if (u >= usat) % Liquid T = fzero(@(x) Water_u(rho, x) - u, [Tsat Tmax]); % Use root-finding method to compute temperature else % Two-Phase with Vapor T = fzero(@(x) H2Oprop('u_dt',rho,x) - u, [Tt Tsat]); % Use root-finding method to compute temperature end else % Cold Liquid or Two-Phase with Vapor Tsat = fzero(@(x) GoodSatProp(x,'rholsat') - rho, [Tt Trhomax]); % Find saturation temperature from density (assume below Trhomax) usat = Water_u(rho, Tsat); % Find saturated liquid energy at this density below Trhomax Tm = fzero(@(x) Water_P(rho, x) - MeltPress(x), [Tmin Tt]); % Find melting temperature at this density um = Water_u(rho, Tm); % Find energy of liquid at melting point with this density if (u >= usat) % Two-Phase with Vapor T = fzero(@(x) H2Oprop('u_dt',rho,x) - u, [Tsat Tmax]); % Use root-finding method to compute temperature elseif (u >= um) % Liquid T = fzero(@(x) Water_u(rho, x) - u, [Tm Tsat]); % Use root-finding method to compute temperature else % Ice T = NaN; % Would need the ice correlation to compute this... end end elseif (rho < rholTmin) % Pressurized Liquid (up to density of liquid at liquid-Ice Ih-Ice III triple point) Tm = fzero(@(x) H2Oprop('d_pt',MeltPress(x),x) - rho, [Tmin Tt]); % Find melting temperature from density um = Water_u(rho,Tm); % Find energy of liquid at melting point with this density if (u >= um) % Liquid T = fzero(@(x) Water_u(rho, x) - u, [Tm Tmax]); % Use root-finding method to compute temperature else % Ice T = NaN; % Would need the ice correlation to compute this... end elseif (rho < rhoicemax) % Very Pressurized Liquid (over 200 MPa) Tm = fzero(@(x) H2Oprop('d_pt',MeltPress(x,'hp'),x) - rho, ... [Tmin Ticemax]); % Find melting temperature from density um = Water_u(rho,Tm); % Find energy of liquid at melting point with this density if (u >= um) % Liquid T = fzero(@(x) Water_u(rho, x) - u, [Tm Tmax]); % Use root-finding method to compute temperature else % High-Pressure Ice T = NaN; % There's no correlation for this stuff... end else % Too Dense (high-pressure ice or over 20 GPa) T = NaN; % There's no correlation for this stuff, and we are well outside the region of validity! end OutVar = T; % Output temperature %--------------------------------------------------------------------------------------- case 'd_pu' % Density from pressure, energy P = Var1; % Input pressure u = Var2; % Input specific internal energy if (P < Pmin) % Too rarefied rho = NaN; % Pressures this low likely violate continuum approximations... elseif (P < Pt) % Thin vapor (below triple point pressure) Tsat = MeltTemp(P); % Compute the sublimation temperature at this pressure rhosat = H2Oprop('d_pt',P,Tsat); % Compute the density of saturated vapor at this pressure usat = Water_u(rhosat, Tsat); % Compute the energy of saturated vapor at this pressure if (u >= usat) % Vapor T = fzero(@(x) u - H2Oprop('u_pt',P,x), [Tsat Tmax]); % Use root-finding method to compute temperature rho = H2Oprop('d_pt',P,T); % Use the temperature and given pressure to find density else % Ice rho = NaN; % Would need the ice correlation to compute this... end elseif (P < Pc) % Liquid, vapor, or two-phase Tm = MeltTemp(P); % Compute the melting temperature at this pressure rhom = H2Oprop('d_pt',P,Tm); % Compute the density of freezing liquid at this pressure um = Water_u(rhom,Tm); % Compute the energy of freezing liquid at this pressure Tsat = H2Oprop('tsat_p',P); % Compute the saturation temperature at this pressure [junk rholsat rhovsat] = GoodSatProp(Tsat); % Compute the saturated liquid and vapor densities at this pressure ulsat = Water_u(rholsat,Tsat); % Compute the saturated liquid energy at this pressure uvsat = Water_u(rhovsat,Tsat); % Compute the saturated vapor energy at this pressure if (u <= um) % Ice rho = NaN; % Would need the ice correlation to compute this... elseif (u <= ulsat) % Liquid T = fzero(@(x) u - H2Oprop('u_pt',P,x), [Tm Tsat]); % Use root-finding method to compute temperature rho = H2Oprop('d_pt',P,T); % Use the temperature and given pressure to find density elseif (u < uvsat) % Two-Phase rho = fzero(@(x) H2Oprop('u_dt',x,Tsat) - u, [rhovsat rholsat]); % Use root-finding method to compute density else % Vapor T = fzero(@(x) u - H2Oprop('u_pt',P,x), [Tsat Tmax]); % Use root-finding method to compute temperature rho = H2Oprop('d_pt',P,T); % Use the temperature and given pressure to find density end elseif (P < Picemax) % Supercritical fluid Tm = MeltTemp(P); % Compute the melting temperature at this pressure rhom = H2Oprop('d_pt',P,Tm); % Compute the density of freezing fluid at this pressure um = Water_u(rhom,Tm); % Compute the energy of freezing fluid at this pressure if (u >= um) % Fluid T = fzero(@(x) u - H2Oprop('u_pt',P,x), [Tm Tmax]); % Use root-finding method to compute temperature rho = H2Oprop('d_pt',P,T); % Use the temperature and given pressure to find density else % Ice rho = NaN; % There's no correlation for this stuff... end else % Too pressurized (cannot verify fluidity) rho = NaN; % There's no correlation for this stuff... end OutVar = rho; % Output density %--------------------------------------------------------------------------------------- case 't_pu' % Temperature from pressure, energy P = Var1; % Input pressure u = Var2; % Input specific internal energy rho = H2Oprop('d_pu',P,u); % Find density from pressure and energy if isnan(rho) % Do the given pressure and energy lie in a valid region? T = NaN; % If not, propagate the NaN else % If so, we can compute temperature T = H2Oprop('t_du',rho,u); % Compute temperature from density and enthalpy end OutVar = T; % Output temperature %--------------------------------------------------------------------------------------- case 'd_ts' % Density from temperature, entropy T = Var1; % Input temperature s = Var2; % Input specific entropy if (T < Tmin) % Cold Vapor Psat = MeltPress(T); % Find sublimation pressure from temperature rhosat = H2Oprop('d_pt',Psat,T); % Find vapor density at sublimation point ssat = Water_s(rhosat,T); % Find vapor entropy at sublimation point if (s >= ssat*.999) % Vapor rho = fzero(@(x) s - Water_s(x,T),[rhomin*.999 rhosat*1.001]); % Use root-finding method to compute density else % Ice rho = NaN; % Would need the ice correlation to compute this... end elseif (T < Tt) % Cold Vapor or Pressurized Liquid [Pm Ps Pf] = MeltPress(T); % Find all phase boundary pressures at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find vapor density at sublimation point rhom = H2Oprop('d_pt',Pm,T); % Find liquid density at Ice Ih melting point rhof = H2Oprop('d_pt',Pf,T); % Find liquid density at high-pressure ice freezing point ss = Water_s(rhos,T); % Find vapor entropy at sublimation point sm = Water_s(rhom,T); % Find liquid entropy at Ice Ih melting point sf = Water_s(rhof,T); % Find liquid entropy at high-pressure ice freezing point if (s < sf) % High-pressure Ice rho = NaN; % There's no correlation for this stuff... elseif (s <= sm) % Pressurized Liquid rho = fzero(@(x) s - Water_s(x,T),[rhom rhof]); % Use root-finding method to compute density elseif (s < ss) % Ice Ih rho = NaN; % Would need the ice correlation to compute this... else % Vapor rho = fzero(@(x) s - Water_s(x,T),[rhomin rhos]); % Use root-finding method to compute density end elseif (T < Tc) % Vapor, Liquid, or Two-Phase [Psat rholsat rhovsat] = GoodSatProp(T); % Find saturation properties at this temperature sl = Water_s(rholsat,T); % Find entropy of saturated liquid sv = Water_s(rhovsat,T); % Find entropy of saturated vapor Pm = MeltPress(T); % Find melting pressure at this temperature rhom = H2Oprop('d_pt',Pm,T); % Find density of liquid at melting point sm = Water_s(rhom,T); % Find entropy of freezing liquid if (s < sm*.999) % High-Pressure Ice rho = NaN; % There's no correlation for this stuff... elseif (s <= sl) % Liquid rho = fzero(@(x) s - Water_s(x,T),[rholsat rhom*1.001]); % Use root-finding method to compute density elseif (s < sv) % Two-Phase Mixture rho = fzero(@(x) s - H2Oprop('s_dt',x,T),[rhovsat rholsat]); % Use root-finding method to compute density else % Vapor rho = fzero(@(x) s - Water_s(x,T),[rhomin rhovsat]); % Use root-finding method to compute density end elseif (T < Ticemax) % Supercritical Fluid Pm = MeltPress(T); % Find melting pressure at this temperature rhom = H2Oprop('d_pt',Pm,T); % Find density of liquid at melting point sm = Water_s(rhom,T); % Find entropy of freezing liquid if (s < (sm - 10*eps(sm))) % High-Pressure Ice rho = NaN; % There's no correlation for this stuff... else % Supercritical Fluid rho = fzero(@(x) s - Water_s(x,T),[rhomin rhom+5*eps(rhom)]); % Use root-finding method to compute density end else % Hot Supercritical Fluid try % Use error handling in case the true density is too large rho = fzero(@(x) s - Water_s(x,T),[rhomin rhoicemax]); % Use root-finding method to compute density catch % If that didn't work... [errmsg,errcode] = lasterr; % Fetch the error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error was the expected failure of fzero... rho = NaN; % ...fail gracefully with a NaN else rethrow(lasterror); % Otherwise, pass the error through end end end OutVar = rho; % Output density %--------------------------------------------------------------------------------------- case 't_ds' % Temperature from density, entropy rho = Var1; % Input density s = Var2; % Input specific entropy if (rho < rhomin) % Too rarefied T = NaN; % Pressures this low likely violate continuum approximations... elseif (rho < rhosatmin) % Extremely Thin Vapor (less dense than 50 K saturated vapor) if (s >= svsatmin) % Vapor T = fzero(@(x) Water_s(rho, x) - s, [Tabsmin Tmax]); % Use root-finding method to compute temperature else % Non-Physical T = NaN; % It's not ice at this density... end elseif (rho < rhovt) % Thin Vapor (less dense than triple point vapor) Tsat = H2Oprop('tsat_dv',rho); % Find sublimation temperature from density ssat = Water_s(rho, Tsat); % Find saturated vapor entropy if (s >= ssat) % Vapor T = fzero(@(x) Water_s(rho, x) - s, [Tsat Tmax]); % Use root-finding method to compute temperature else % Ice T = NaN; % Would need the ice correlation to compute this... end elseif (rho < rhoc) % Vapor or Two-Phase with Liquid Tsat = H2Oprop('tsat_dv',rho); % Find saturation temperature from density ssat = Water_s(rho, Tsat); % Find saturated vapor entropy Xmin = H2Oprop('x_dt',rho,Tt); % Find quality for this density at the triple point smin = Xmin*svt + (1-Xmin)*slt; % Find entropy for two-phase mixture at the triple point if (s >= ssat) % Vapor T = fzero(@(x) Water_s(rho, x) - s, [Tsat Tmax]); % Use root-finding method to compute temperature elseif (s >= smin) % Two-Phase with Liquid T = fzero(@(x) H2Oprop('s_dt',rho,x) - s, [Tt Tsat]); % Use root-finding method to compute temperature else % Ice T = NaN; % Would need the ice correlation to compute this... end elseif (rho < rholt) % Liquid or Two-Phase with Vapor (No Density Anomaly) Tsat = H2Oprop('tsat_dl',rho); % Find saturation temperature from density ssat = Water_s(rho, Tsat); % Find saturated liquid entropy Xmin = H2Oprop('x_dt',rho,Tt); % Find quality for this density at the triple point smin = Xmin*svt + (1-Xmin)*slt; % Find entropy for two-phase mixture at the triple point if (s >= ssat) % Liquid T = fzero(@(x) Water_s(rho, x) - s, [Tsat Tmax]); % Use root-finding method to compute temperature elseif (s >= smin) % Two-Phase T = fzero(@(x) H2Oprop('s_dt',rho,x) - s, [Tt Tsat]); % Use root-finding method to compute temperature else % Ice T = NaN; % Would need the ice correlation to compute this... end elseif (rho < rhoamax) % Liquid or Two-Phase with Vapor (Density Anomaly) if s >= sa % Warm Liquid or Two-Phase with Vapor Tsat = H2Oprop('tsat_dl',rho); % Find saturation temperature from density (assume above Trhomax) ssat = Water_s(rho, Tsat); % Find saturated liquid entropy if (s >= ssat) % Liquid T = fzero(@(x) Water_s(rho, x) - s, [Tsat Tmax]); % Use root-finding method to compute temperature else % Two-Phase with Vapor T = fzero(@(x) H2Oprop('s_dt',rho,x) - s, [Tt Tsat]); % Use root-finding method to compute temperature end else % Cold Liquid or Two-Phase Mixture Tsat = fzero(@(x) GoodSatProp(x,'rholsat') - rho, [Tt Trhomax]); % Find saturation temperature from density (assume below Trhomax) ssat = Water_s(rho, Tsat); % Find saturated liquid entropy Tm = fzero(@(x) Water_P(rho, x) - MeltPress(x), [Tmin Tt]); % Find melting temperature at this density sm = Water_s(rho, Tm); % Find entropy of liquid at freezing point if (s >= ssat) % Two-Phase with Vapor T = fzero(@(x) H2Oprop('s_dt',rho,x) - s, [Tsat Tmax]); % Use root-finding method to compute temperature elseif (s >= sm) % Liquid T = fzero(@(x) Water_s(rho, x) - s, [Tm Tsat]); % Use root-finding method to compute temperature else % Ice T = NaN; % Would need the ice correlation to compute this... end end elseif (rho < rholTmin) % Pressurized Liquid (up to density of liquid at liquid-Ice Ih-Ice III triple point) Tm = fzero(@(x) H2Oprop('d_pt',MeltPress(x),x) - rho, [Tmin Tt]); % Find melting temperature from density sm = Water_s(rho,Tm); % Find entropy of liquid at melting point if (s >= sm) % Liquid T = fzero(@(x) Water_s(rho, x) - s, [Tm Tmax]); % Use root-finding method to compute temperature else % Ice Ih T = NaN; % Would need the ice correlation to compute this... end elseif (rho <= rhoicemax) % Very Pressurized Liquid (over 200 MPa) Tm = fzero(@(x) H2Oprop('d_pt',MeltPress(x,'hp'),x) - rho, ... [Tmin Ticemax]); % Find melting temperature from density sm = Water_s(rho,Tm); % Find entropy of liquid at melting point if (s >= sm) % Liquid T = fzero(@(x) Water_s(rho, x) - s, [Tm Tmax]); % Use root-finding method to compute temperature else % High-Pressure Ice T = NaN; % There's no correlation for this stuff... end else % Too Dense (high-pressure ice or over 20 GPa) T = NaN; % There's no correlation for this stuff, and we are well outside the region of validity! end OutVar = T; % Output temperature %--------------------------------------------------------------------------------------- case 'p_du' % Pressure from density, energy rho = Var1; % Input density u = Var2; % Input specific internal energy T = H2Oprop('t_du',rho,u); % Calculate temperature from density and energy if isnan(T) % Was the temperature calculation successful? P = NaN; % If not, propagate the NaN else P = H2Oprop('p_dt',rho,T); % If so, use the temperature and density to find pressure end OutVar = P; % Output pressure %--------------------------------------------------------------------------------------- case 'p_dh' % Pressure from density, enthalpy rho = Var1; % Input density h = Var2; % Input specific enthalpy T = H2Oprop('t_dh',rho,h); % Calculate temperature from density and enthalpy if isnan(T) % Was the temperature calculation successful? P = NaN; % If not, propagate the NaN else P = H2Oprop('p_dt',rho,T); % If so, use the temperature and density to find pressure end OutVar = P; % Output pressure %--------------------------------------------------------------------------------------- case 'p_ds' % Pressure from density, entropy rho = Var1; % Input density s = Var2; % Input specific entropy T = H2Oprop('t_ds',rho,s); % Calculate temperature from density and entropy if isnan(T) % Was the temperature calculation successful? P = NaN; % If not, propagate the NaN else P = H2Oprop('p_dt',rho,T); % If so, use the temperature and density to find pressure end OutVar = P; % Output pressure %--------------------------------------------------------------------------------------- case 'p_ts' % Pressure from temperature, entropy T = Var1; % Input temperature s = Var2; % Input specific entropy rho = H2Oprop('d_ts',T,s); % Calculate density from temperature and entropy if isnan(rho) % Was the density calculation successful? P = NaN; % If not, propagate the NaN else P = H2Oprop('p_dt',rho,T); % If so, use the temperature and density to find pressure end OutVar = P; % Output pressure %--------------------------------------------------------------------------------------- case 'u_dh' % Specific energy from density, enthalpy rho = Var1; % Input density h = Var2; % Input specific enthalpy T = H2Oprop('t_dh',rho,h); % Calculate temperature from density and enthalpy if isnan(T) % Was the temperature calculation successful? u = NaN; % If not, propagate the NaN else u = H2Oprop('u_dt',rho,T); % If so, use the temperature and density to find specific internal energy end OutVar = u; % Output specific internal energy %--------------------------------------------------------------------------------------- case 'u_ds' % Specific energy from density, entropy rho = Var1; % Input density s = Var2; % Input specific entropy T = H2Oprop('t_ds',rho,s); % Calculate temperature from density and entropy if isnan(T) % Was the temperature calculation successful? u = NaN; % If not, propagate the NaN else u = H2Oprop('u_dt',rho,T); % If so, use the temperature and density to find specific internal energy end OutVar = u; % Output specific internal energy %--------------------------------------------------------------------------------------- case 'u_pt' % Specific energy from pressure, temperature P = Var1; % Input pressure T = Var2; % Input temperature rho = H2Oprop('d_pt',P,T); % Calculate density from pressure and temperature if isnan(rho) % Was the density calculation successful? u = NaN; % If not, propagate the NaN else u = H2Oprop('u_dt',rho,T); % If so, use the temperature and density to find specific internal energy end OutVar = u; % Output specific internal energy %--------------------------------------------------------------------------------------- case 'u_ph' % Specific energy from pressure, enthalpy P = Var1; % Input pressure h = Var2; % Input specific enthalpy T = H2Oprop('t_ph',P,h); % Calculate temperature from pressure and enthalpy rho = H2Oprop('d_ph',P,h); % Calculate density from pressure and enthalpy if isnan(T)||isnan(rho) % Were the temperature and density calculations successful? u = NaN; % If not, propagate the NaN else u = H2Oprop('u_dt',rho,T); % If so, use the temperature and density to find specific internal energy end OutVar = u; % Output specific internal energy %--------------------------------------------------------------------------------------- case 'u_ts' % Specific energy from entropy, temperature T = Var1; % Input temperature s = Var2; % Input specific entropy rho = H2Oprop('d_ts',T,s); % Calculate density from temperature and entropy if isnan(rho) % Was the density calculation successful? u = NaN; % If not, propagate the NaN else u = H2Oprop('u_dt',rho,T); % If so, use the temperature and density to find specific internal energy end OutVar = u; % Output specific internal energy %--------------------------------------------------------------------------------------- case 'h_du' % Specific enthalpy from density, energy rho = Var1; % Input density u = Var2; % Input specific internal energy T = H2Oprop('t_du',rho,u); % Calculate temperature from density and internal energy if isnan(T) % Was the temperature calculation successful? h = NaN; % If not, propagate the NaN else h = H2Oprop('h_dt',rho,T); % If so, use the temperature and density to find specific enthalpy end OutVar = h; % Output specific enthalpy %--------------------------------------------------------------------------------------- case 'h_ds' % Specific enthalpy from density, entropy rho = Var1; % Input density s = Var2; % Input specific entropy T = H2Oprop('t_ds',rho,s); % Calculate temperature from density and entropy if isnan(T) % Was the temperature calculation successful? h = NaN; % If not, propagate the NaN else h = H2Oprop('h_dt',rho,T); % If so, use the temperature and density to find specific enthalpy end OutVar = h; % Output specific enthalpy %--------------------------------------------------------------------------------------- case 'h_pt' % Specific enthalpy from pressure, temperature P = Var1; % Input pressure T = Var2; % Input temperature rho = H2Oprop('d_pt',P,T); % Calculate density from pressure and temperature if isnan(rho) % Was the density calculation successful? h = NaN; % If not, propagate the NaN else h = H2Oprop('h_dt',rho,T); % If so, use the temperature and density to find specific enthalpy end OutVar = h; % Output specific enthalpy %--------------------------------------------------------------------------------------- case 'h_pu' % Specific enthalpy from pressure, energy P = Var1; % Input pressure u = Var2; % Input specific internal energy T = H2Oprop('t_pu',P,u); % Calculate temperature from pressure and energy rho = H2Oprop('d_pu',P,u); % Calculate density from pressure and energy if isnan(T)||isnan(rho) % Were the density and temperature calculations successful? h = NaN; % If not, propagate the NaN else h = H2Oprop('h_dt',rho,T); % If so, use the temperature and density to find specific enthalpy end OutVar = h; % Output specific enthalpy %--------------------------------------------------------------------------------------- case 'h_ts' % Specific enthalpy from entropy, temperature T = Var1; % Input temperature s = Var2; % Input specific entropy rho = H2Oprop('d_ts',T,s); % Calculate density from temperature and entropy if isnan(rho) % Was the density calculation successful? h = NaN; % If not, propagate the NaN else h = H2Oprop('h_dt',rho,T); % If so, use the temperature and density to find specific enthalpy end OutVar = h; % Output specific enthalpy %--------------------------------------------------------------------------------------- case 's_dh' % Specific entropy from density, enthalpy rho = Var1; % Input density h = Var2; % Input specific enthalpy T = H2Oprop('t_dh',rho,h); % Calculate temperature from density and enthalpy if isnan(T) % Was the temperature calculation successful? s = NaN; % If not, propagate the NaN else s = H2Oprop('s_dt',rho,T); % If so, use the temperature and density to find specific entropy end OutVar = s; % Output specific entropy %--------------------------------------------------------------------------------------- case 's_du' % Specific entropy from density, energy rho = Var1; % Input density u = Var2; % Input specific internal energy T = H2Oprop('t_du',rho,u); % Calculate temperature from density and energy if isnan(T) % Was the temperature calculation successful? s = NaN; % If not, propagate the NaN else s = H2Oprop('s_dt',rho,T); % If so, use the temperature and density to find specific entropy end OutVar = s; % Output specific entropy %--------------------------------------------------------------------------------------- case 's_pt' % Specific entropy from pressure, temperature P = Var1; % Input pressure T = Var2; % Input temperature rho = H2Oprop('d_pt',P,T); % Calculate density from pressure and temperature if isnan(rho) % Was the density calculation successful? s = NaN; % If not, propagate the NaN else s = H2Oprop('s_dt',rho,T); % If so, use the temperature and density to find specific entropy end OutVar = s; % Output specific entropy %--------------------------------------------------------------------------------------- case 's_ph' % Specific entropy from pressure, enthalpy P = Var1; % Input pressure h = Var2; % Input specific enthalpy T = H2Oprop('t_ph',P,h); % Calculate temperature from pressure and enthalpy rho = H2Oprop('d_ph',P,h); % Calculate density from pressure and enthalpy if isnan(T)||isnan(rho) % Were the density and temperature calculations successful? s = NaN; % If not, propagate the NaN else s = H2Oprop('s_dt',rho,T); % If so, use the temperature and density to find specific entropy end OutVar = s; % Output specific entropy %--------------------------------------------------------------------------------------- case 's_pu' % Specific entropy from pressure, energy P = Var1; % Input pressure u = Var2; % Input specific internal energy T = H2Oprop('t_pu',P,u); % Calculate temperature from pressure and energy rho = H2Oprop('d_pu',P,u); % Calculate density from pressure and energy if isnan(T)||isnan(rho) % Were the density and temperature calculations successful? s = NaN; % If not, propagate the NaN else s = H2Oprop('s_dt',rho,T); % If so, use the temperature and density to find specific entropy end OutVar = s; % Output specific entropy %--------------------------------------------------------------------------------------- case 't_ps' % Temperature from pressure, entropy P = Var1; % Input pressure s = Var2; % Input specific entropy if (P < Pmin) % Too rarefied T = NaN; % Pressures this low likely violate continuum approximations... elseif (P < Pt) % Thin vapor (below triple point pressure) Tsat = MeltTemp(P); % Find sublimation temperature at this pressure rhosat = H2Oprop('d_pt',P,Tsat); % Find vapor density at sublimation point ssat = Water_s(rhosat, Tsat); % Find vapor entropy at sublimation point if (s >= ssat) % Vapor T = fzero(@(x) s - H2Oprop('s_pt',P,x), [Tsat Tmax]); % Use root-finding method to find temperature else % Ice T = NaN; % Would need the ice correlation to calculate this... end elseif (P < Pc) % Liquid or vapor Tm = MeltTemp(P); % Find Ice Ih melting point at this pressure rhom = H2Oprop('d_pt',P,Tm); % Find liquid density at melting point sm = Water_s(rhom,Tm); % Find liquid entropy at melting point Tsat = H2Oprop('tsat_p',P); % Find saturation temperature at this pressure [junk rholsat rhovsat] = GoodSatProp(Tsat); % Calculate saturated liquid and vapor densities at this pressure slsat = Water_s(rholsat,Tsat); % Find entropy of saturated liquid svsat = Water_s(rhovsat,Tsat); % Find entropy of saturated vapor if (s <= sm) % Ice T = NaN; % Would need the ice correlation to calculate this... elseif (s <= slsat) % Liquid T = fzero(@(x) s - H2Oprop('s_pt',P,x), [Tm Tsat]); % Use root-finding method to find temperature elseif (s < svsat) % Two-Phase Mixture T = Tsat; % Return the saturation temperature else % Vapor T = fzero(@(x) s - H2Oprop('s_pt',P,x), [Tsat Tmax]); % Use root-finding method to find temperature end elseif (P < Picemax) % Supercritical Fluid Tm = MeltTemp(P); % Find the melting point at this pressure rhom = H2Oprop('d_pt',P,Tm); % Find liquid density at the melting point sm = Water_s(rhom,Tm); % Find liquid entropy at the melting point if (s >= sm) % Fluid T = fzero(@(x) s - H2Oprop('s_pt',P,x), [Tm Tmax]); % Use root-finding method to find temperature else % Ice T = NaN; % Would need the ice correlation to calculate this... end else % Too pressurized (cannot verify fluidity) T = NaN; % There's no correlation for this stuff... end OutVar = T; % Output temperature %--------------------------------------------------------------------------------------- case 'd_ps' % Density from pressure, entropy P = Var1; % Input pressure s = Var2; % Input specific entropy T = H2Oprop('t_ps',P,s); % Calculate temperature from pressure and entropy if isnan(T) % Did the temperature calculation work? rho = NaN; % If not, propagate the NaN else rho = H2Oprop('d_ts',T,s); % If so, use the temperature and entropy to find the density end OutVar = rho; % Output density %--------------------------------------------------------------------------------------- case 'u_ps' % Specific energy from pressure, entropy P = Var1; % Input pressure s = Var2; % Input specific entropy T = H2Oprop('t_ps',P,s); % Calculate temperature from pressure and entropy rho = H2Oprop('d_ts',T,s); % Calculate density from temperature and entropy if isnan(T) % Did the temperature calculation work? u = NaN; % If not, propagate the NaN else u = H2Oprop('u_dt',rho,T); % If so, use the temperature and density to find the internal energy end OutVar = u; % Output specific internal energy %--------------------------------------------------------------------------------------- case 'h_ps' % Specific enthalpy from pressure, entropy P = Var1; % Input pressure s = Var2; % Input specific entropy T = H2Oprop('t_ps',P,s); % Calculate temperature from pressure and entropy rho = H2Oprop('d_ts',T,s); % Calculate density from temperature and entropy if isnan(T) % Did the temperature calculation work? h = NaN; % If not, propagate the NaN else h = H2Oprop('h_dt',rho,T); % If so, use the temperature and density to find the enthalpy end OutVar = h; % Output specific enthalpy %--------------------------------------------------------------------------------------- case 't_hs' % Temperature from enthalpy, entropy % Note that there are local maxima and minima in entropy along the high-pressure melting curve. % These are outside the official pressure range, and likely represent % model inaccuaracies rather than a physical phenomenon. This function % does some gymnastics around these local extrema for numerical % reasons, but one should be cautious in interpreting results for fluid % states above 8.9 GPa and 600 K. h = Var1; % Input specific enthalpy s = Var2; % Input specific entropy if (s < spimin) % Ice T = NaN; % Would need the ice correlation to compute this... elseif (s < slt) % Ice or Liquid Tm = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x),x), ... [Tmin+100*eps(Tmin) Tt]); % Use root-finding method to find Ice Ih melting point from entropy Tf = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ... [Tmin+100*eps(Tmin) Tsmax]); % Use root-finding method to find high-pressure ice melting point from entropy try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tm Tf]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end elseif (s < spitmin) % Ice, Liquid, or Two-Phase with Vapor (below anomalies near maximum of melting curve) Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy Tf = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), [Tt Tsmax]); % Use root-finding method to find high-pressure ice melting point from entropy hl = H2Oprop('hlsat_t',Tv); % Find enthalpy of saturated liquid if (h < hl) % Two-Phase or Ice try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end else % Liquid or High-Pressure Ice try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tv Tf]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end end elseif (s <= spitmax) % Ice, Liquid, or Two-Phase with Vapor (between entropy minimum and entropy at top of melting curve) Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy Tfl = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ... [Tt Tsmax]); % Use root-finding method to find high-pressure ice melting point below entropy maximum from entropy Tfh = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ... [Tsmax Tsmin]); % Use root-finding method to find Ice VII melting point between entropy maximum and entropy minimum from entropy Tfm = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ... [Tsmin Ticemax]); % Use root-finding method to find Ice VII melting point above entropy minimum from entropy hl = H2Oprop('hlsat_t',Tv); % Find enthalpy of saturated liquid hil = H2Oprop('h_pt',MeltPress(Tfl,'hp'),Tfl); % Find enthalpy of freezing liquid below entropy maximum hih = H2Oprop('h_pt',MeltPress(Tfh,'hp'),Tfh); % Find enthalpy of freezing liquid between entropy extrema if (h < hl) % Two-Phase with Vapor or Ice try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end elseif (h < hil) % Liquid (below anomalies) T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tv Tfl]); % Use root-finding method to find temperature from enthalpy and entropy elseif (h < hih) % Ice VII (between anomalies) T = NaN; % There's no correlation for this stuff... else % Hot Liquid (above anomalies) try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tfh Tfm]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end end elseif (s <= spimax) % Ice, Liquid, or Two-Phase with Vapor (between entropy at top of melting curve and entropy maximum) Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy Tfl = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ... [Tt Tsmax]); % Use root-finding method to find high-pressure ice melting point below entropy maximum from entropy Tfh = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ... [Tsmax Ticemax]); % Use root-finding method to find high-pressure ice melting point above entropy maximum from entropy hl = H2Oprop('hlsat_t',Tv); % Find enthalpy of saturated liquid hil = H2Oprop('h_pt',MeltPress(Tfl,'hp'),Tfl); % Find enthalpy of freezing liquid below entropy maximum hih = H2Oprop('h_pt',MeltPress(Tfh,'hp'),Tfh); % Find enthalpy of freezing liquid between entropy extrema if (h < hl) % Two-Phase with Vapor or Ice try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end elseif (h < hil) % Liquid (below anomalies) T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tv Tfl]); % Use root-finding method to find temperature from enthalpy and entropy elseif (h < hih) % Ice VII (between anomalies) T = NaN; % There's no correlation for this stuff... else % Hot Liquid (above anomalies) Tfm = H2Oprop('t_ds',rhoicemax,s); % Find temperature at which material at rhoicemax has this entropy try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tfh Tfm]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range, or if the state has density > rhoicemax [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN elseif isnan(Tfm) % If this h,s point is at a density above rhoicemax, Tfm would be NaN and fzero would also have failed T = NaN; % Density too great, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end end elseif (s <= sdtmax) % Liquid or Two-Phase with Vapor (below entropy at maximum density and temperature) Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy hl = H2Oprop('hlsat_t',Tv); % Find enthalpy of saturated liquid if (h > hl) % Liquid Tm = H2Oprop('t_ds',rhoicemax,s); % Find temperature at which material at rhoicemax has this entropy try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tv Tm-10*eps(Tm)]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range, or if the state has density > rhoicemax [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN elseif isnan(Tm) % If this h,s point is at a density above rhoicemax, Tm would be NaN and fzero would also have failed T = NaN; % Density too great, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end else % Two-Phase with Vapor try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end end elseif (s <= scrit) % Liquid or Two-Phase with Vapor (above entropy at maximum density and temperature) try Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy catch % The root-finding would fail if the actual entropy were very close to the critical entropy [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" boiling temperature was too close to Tc... Tv = Tc; % Just set the boiling point to equal the critical point else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end hl = H2Oprop('hlsat_t',Tv); % Find enthalpy of saturated liquid if (h > hl) % Liquid try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tt Tmax*.999]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end else % Two-Phase with Vapor try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end end elseif (s <= svt) % Vapor or Two-Phase with Liquid try Tv = fzero(@(x) s - H2Oprop('svsat_t',x), [Tabsmin Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy catch % The root-finding would fail if the actual entropy were very close to the critical entropy [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" boiling temperature was too close to Tc... Tv = Tc; % Just set the boiling point to equal the critical point else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end hv = H2Oprop('hvsat_t',Tv); % Find entropy of saturated vapor if (h < hv) % Two-Phase with Liquid try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end else % Vapor try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tv+1e5*eps(Tv) Tmax]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end end elseif (s <= svsatmin) % Vapor or Ice Ts = fzero(@(x) s - H2Oprop('svsat_t',x), [Tabsmin Tt]); % Use root-finding method to find sublimation point from entropy try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Ts+1e5*eps(Ts) Tmax]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end else % Only vapor allowed try T = fzero(@(x) h - H2Oprop('h_ts',x,s), [Tabsmin Tmax]); % Use root-finding method to find temperature from enthalpy and entropy catch % The root-finding would fail if enthalpy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end end OutVar = T; % Output temperature %--------------------------------------------------------------------------------------- case 'd_hs' % Density from enthalpy, entropy h = Var1; % Input specific enthalpy s = Var2; % Input specific entropy T = H2Oprop('t_hs',h,s); % Calculate temperature from enthalpy and entropy if isnan(T) % If that didn't work... rho = NaN; % ...propagate the NaN else rho = H2Oprop('d_ts',T,s); % Otherwise, calculate density from temperature and entropy end OutVar = rho; % Output density %--------------------------------------------------------------------------------------- case 'p_hs' % Pressure from enthalpy, entropy h = Var1; % Input specific enthalpy s = Var2; % Input specific entropy T = H2Oprop('t_hs',h,s); % Calculate temperature from enthalpy and entropy if isnan(T) % If that didn't work... P = NaN; % ...propagate the NaN else P = H2Oprop('p_ts',T,s); % Otherwise, calculate pressure from temperature and entropy end OutVar = P; % Output density %--------------------------------------------------------------------------------------- case 'u_hs' % Specific Energy from enthalpy, entropy h = Var1; % Input specific enthalpy s = Var2; % Input specific entropy T = H2Oprop('t_hs',h,s); % Calculate temperature from enthalpy and entropy if isnan(T) % If that didn't work... u = NaN; % ...propagate the NaN else u = H2Oprop('u_ts',T,s); % Otherwise, calculate specific internal energy from temperature and entropy end OutVar = u; % Output specific internal energy %--------------------------------------------------------------------------------------- case 't_us' % Temperature from energy, entropy u = Var1; % Input specific internal energy s = Var2; % Input specific entropy if (s < spimin) % Ice T = NaN; % Would need the ice correlation to compute this... elseif (s < slt) % Ice or Liquid Tm = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x),x), ... [Tmin+100*eps(Tmin) Tt]); % Use root-finding method to find Ice Ih melting point from entropy Tf = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ... [Tmin+100*eps(Tmin) Tsmax]); % Use root-finding method to find high-pressure ice melting point from entropy try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tm Tf]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end elseif (s < spitmin) % Ice, Liquid, or Two-Phase with Vapor (below anomalies near maximum of melting curve) Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy Tf = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), [Tt Tsmax]); % Use root-finding method to find high-pressure ice melting point from entropy ul = H2Oprop('ulsat_t',Tv); % Find energy of saturated liquid if (u < ul) % Two-Phase or Ice try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end else % Liquid or High-Pressure Ice try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tv Tf]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end end elseif (s <= spitmax) % Ice, Liquid, or Two-Phase with Vapor (between entropy minimum and entropy at top of melting curve) Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy Tfl = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ... [Tt Tsmax]); % Use root-finding method to find high-pressure ice melting point below entropy maximum from entropy Tfh = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ... [Tsmax Tsmin]); % Use root-finding method to find Ice VII melting point between entropy maximum and entropy minimum from entropy Tfm = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ... [Tsmin Ticemax]); % Use root-finding method to find Ice VII melting point above entropy minimum from entropy ul = H2Oprop('ulsat_t',Tv); % Find energy of saturated liquid uil = H2Oprop('u_pt',MeltPress(Tfl,'hp'),Tfl); % Find energy of freezing liquid below entropy maximum uih = H2Oprop('u_pt',MeltPress(Tfh,'hp'),Tfh); % Find energy of freezing liquid between entropy extrema if (u < ul) % Two-Phase with Vapor or Ice try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end elseif (u < uil) % Liquid (below anomalies) T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tv Tfl]); % Use root-finding method to find temperature from energy and entropy elseif (u < uih) % Ice VII (between anomalies) T = NaN; % There's no correlation for this stuff... else % Hot Liquid (above anomalies) try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tfh Tfm]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end end elseif (s <= spimax) % Ice, Liquid, or Two-Phase with Vapor (between entropy at top of melting curve and entropy maximum) Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy Tfl = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ... [Tt Tsmax]); % Use root-finding method to find high-pressure ice melting point below entropy maximum from entropy Tfh = fzero(@(x) s - H2Oprop('s_pt',MeltPress(x,'hp'),x), ... [Tsmax Ticemax]); % Use root-finding method to find high-pressure ice melting point above entropy maximum from entropy ul = H2Oprop('ulsat_t',Tv); % Find energy of saturated liquid uil = H2Oprop('u_pt',MeltPress(Tfl,'hp'),Tfl); % Find energy of freezing liquid below entropy maximum uih = H2Oprop('u_pt',MeltPress(Tfh,'hp'),Tfh); % Find energy of freezing liquid between entropy extrema if (u < ul) % Two-Phase with Vapor or Ice try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end elseif (u < uil) % Liquid (below anomalies) T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tv Tfl]); % Use root-finding method to find temperature from energy and entropy elseif (u < uih) % Ice VII (between anomalies) T = NaN; % There's no correlation for this stuff... else % Hot Liquid (above anomalies) Tfm = H2Oprop('t_ds',rhoicemax,s); % Find temperature at which material at rhoicemax has this entropy try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tfh Tfm]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range, or if the state has density > rhoicemax [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN elseif isnan(Tfm) % If this u,s point is at a density above rhoicemax, Tfm would be NaN and fzero would also have failed T = NaN; % Density too great, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end end elseif (s <= sdtmax) % Liquid or Two-Phase with Vapor (below entropy at maximum density and temperature) Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy ul = H2Oprop('ulsat_t',Tv); % Find energy of saturated liquid if (u > ul) % Liquid Tm = H2Oprop('t_ds',rhoicemax,s); % Find temperature at which material at rhoicemax has this entropy try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tv Tm-10*eps(Tm)]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range, or if the state has density > rhoicemax [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN elseif isnan(Tm) % If this u,s point is at a density above rhoicemax, Tm would be NaN and fzero would also have failed T = NaN; % Density too great, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end else % Two-Phase with Vapor try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end end elseif (s <= scrit) % Liquid or Two-Phase with Vapor (above entropy at maximum density and temperature) try Tv = fzero(@(x) s - H2Oprop('slsat_t',x), [Tt Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy catch % The root-finding would fail if the actual entropy were very close to the critical entropy [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" boiling temperature was too close to Tc... Tv = Tc; % Just set the boiling point to equal the critical point else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end ul = H2Oprop('ulsat_t',Tv); % Find energy of saturated liquid if (u > ul) % Liquid try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tt Tmax*.999]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end else % Two-Phase with Vapor try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end end elseif (s <= svt) % Vapor or Two-Phase with Liquid try Tv = fzero(@(x) s - H2Oprop('svsat_t',x), [Tabsmin Tc-1.01e-5]); % Use root-finding method to find boiling point from entropy catch % The root-finding would fail if the actual entropy were very close to the critical entropy [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" boiling temperature was too close to Tc... Tv = Tc; % Just set the boiling point to equal the critical point else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end uv = H2Oprop('uvsat_t',Tv); % Find energy of saturated vapor if (u < uv) % Two-Phase with Liquid try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tt Tv]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end else % Vapor try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tv+1e5*eps(Tv) Tmax]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end end elseif (s <= svsatmin) % Vapor or Ice Ts = fzero(@(x) s - H2Oprop('svsat_t',x), [Tabsmin Tt]); % Use root-finding method to find sublimation point from entropy try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Ts+1e5*eps(Ts) Tmax]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end else % Only vapor allowed try T = fzero(@(x) u - H2Oprop('u_ts',x,s), [Tabsmin Tmax]); % Use root-finding method to find temperature from energy and entropy catch % The root-finding would fail if energy were outside the physical range [errmsg,errcode] = lasterr; % Retrieve error code if strcmp(errcode, 'MATLAB:fzero:ValuesAtEndPtsSameSign') % If the error indicates that the "true" temperature was outside the fluid range... T = NaN; % Ice or Non-Physical, so return NaN else % Other kinds of errors indicate real problems... rethrow(lasterror); % ...so kick the error upstairs end end end OutVar = T; %--------------------------------------------------------------------------------------- case 'd_us' % Density from energy, entropy u = Var1; % Input specific internal energy s = Var2; % Input specific entropy T = H2Oprop('t_us',u,s); % Calculate temperature from energy and entropy if isnan(T) % If that didn't work... rho = NaN; % ...propagate the NaN else rho = H2Oprop('d_ts',T,s); % Otherwise, calculate density from temperature and entropy end OutVar = rho; % Output density %--------------------------------------------------------------------------------------- case 'p_us' % Pressure from energy, entropy u = Var1; % Input specific internal energy s = Var2; % Input specific entropy T = H2Oprop('t_us',u,s); % Calculate temperature from energy and entropy if isnan(T) % If that didn't work... P = NaN; % ...propagate the NaN else P = H2Oprop('p_ts',T,s); % Otherwise, calculate pressure from temperature and entropy end OutVar = P; % Output pressure %--------------------------------------------------------------------------------------- case 'h_us' % Specific Enthalpy from energy, entropy u = Var1; % Input specific internal energy s = Var2; % Input specific entropy T = H2Oprop('t_us',u,s); % Calculate temperature from energy and entropy if isnan(T) % If that didn't work... h = NaN; % ...propagate the NaN else h = H2Oprop('h_ts',T,s); % Otherwise, calculate enthalpy from temperature and entropy end OutVar = h; % Output specific enthalpy %--------------------------------------------------------------------------------------- case 'x_dt' % Quality from density, temperature rho = Var1; % Input density T = Var2; % Input temperature [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature if isnan(Psat) % Did the given temperature correspond to a point on the saturation line? X = NaN; % If not, return NaN elseif rho < rhovsat % Is the given state a pure vapor? X = 1; % If so, the quality is one elseif rho > rholsat % Is the given state a pure liquid? X = 0; % If so, the quality is zero else % We must actually have a two-phase mixture, then v = 1/rho; % Convert given density to specific volume vl = 1/rholsat; % Convert saturated liquid density to specific volume vv = 1/rhovsat; % Convert saturated vapor density to specific volume X = (v - vl)/(vv - vl); % Compute the quality from the specific volumes end OutVar = X; % Output quality %--------------------------------------------------------------------------------------- case 'cp_dt' % Constant Pressure Specific Heat from density, temperature rho = Var1; % Input density T = Var2; % Input temperature if (T > Ticemax) % Supercritical Fluid cp = Water_cp(rho,T); % Directly compute specific heat elseif (T > Tt) % Regular Fluid [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature Pm = MeltPress(T); % Find the melting pressure at this temperature rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point if (rho > rhom) % High-Pressure Ice cp = NaN; % There's no correlation for this stuff... elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase Mixture cp = Inf; % By definition, the constant-pressure specific heat at saturation is infinite else % Single Fluid Phase cp = Water_cp(rho,T); % Directly compute specific heat end elseif (T > Tmin) % Cold Fluid [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase cp = Water_cp(rho,T); % Directly compute specific heat else % Ice cp = NaN; % Would need an appropriate correlation... end elseif (T > Tabsmin) % Cold Vapor Ps = MeltPress(T); % Find the sublimation pressure at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point if (rho <= rhos) % Vapor cp = Water_cp(rho,T); % Directly compute specific heat else % Ice Ih cp = NaN; % Would need the ice correlation to compute this... end end OutVar = cp; % Output constant-pressure specific heat %--------------------------------------------------------------------------------------- case 'cv_dt' % Constant Volume Specific Heat from density, temperature rho = Var1; % Input density T = Var2; % Input temperature if (T > Ticemax) % Supercritical Fluid cv = Water_cv(rho,T); % Directly compute specific heat elseif (T > Tt) % Regular Fluid [Psat rholsat rhovsat] = GoodSatProp(T); % Find saturation properties at this temperature Pm = MeltPress(T); % Find melting pressure at this temperature rhom = H2Oprop('d_pt',Pm,T); % Find liquid density at melting point if (rho > rhom) % Ice cv = NaN; % Would need an appropriate correlation... elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase cv = (H2Oprop('u_dt',rho,T+1e6*eps(T)) - ... H2Oprop('u_dt',rho,T))/1e6/eps(T); % Use an approximate numerical derivative to find the specific heat else % Single Phase cv = Water_cv(rho,T); % Directly compute specific heat end elseif (T > Tmin) % Cold Fluid [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase cv = Water_cv(rho,T); % Directly compute specific heat else % Ice cv = NaN; % Would need an appropriate correlation... end elseif (T > Tabsmin) % Cold Vapor Ps = MeltPress(T); % Find the sublimation pressure at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point if (rho <= rhos) % Vapor cv = Water_cv(rho,T); % Directly compute specific heat else % Ice Ih cv = NaN; % Would need the ice correlation to compute this... end end OutVar = cv; % Output constant volume specific heat %--------------------------------------------------------------------------------------- case 'betat_dt' % Thermal Expansion Coefficient from density, temperature rho = Var1; % Input density T = Var2; % Input temperature if (T > Ticemax) % Supercritical Fluid betat = Water_betat(rho,T); % Directly compute thermal expansion coefficient elseif (T > Tt) % Ordinary Fluid Pm = MeltPress(T); % Find the melting pressure at this temperature rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature if (rho > rhom) % Ice betat = NaN; % There's no correlation for this stuff... elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase betat = NaN; % Thermal expansion is poorly defined for a saturated two-phase mixture else % Single Phase betat = Water_betat(rho,T); % Directly compute thermal expansion coefficient end elseif (T > Tmin) % Cold Fluid [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase betat = Water_betat(rho,T); % Directly compute thermal expansion coefficient else % Ice betat = NaN; % Would need an appropriate correlation... end elseif (T > Tabsmin) % Cold Vapor Ps = MeltPress(T); % Find the sublimation pressure at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point if (rho <= rhos) % Vapor betat = Water_betat(rho,T); % Directly compute thermal expansion coefficient else % Ice Ih betat = NaN; % Would need the ice correlation to compute this... end end OutVar = betat; % Output thermal expansion coefficient %--------------------------------------------------------------------------------------- case 'kappat_dt' % Isothermal Compressibility from density, temperature rho = Var1; % Input density T = Var2; % Input temperature if (T > Ticemax) % Supercritical Fluid kappat = Water_kappat(rho,T); % Directly compute compressibility elseif (T > Tt) % Ordinary Fluid Pm = MeltPress(T); % Find the melting pressure at this temperature rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature if (rho > rhom) % High-Pressure Ice kappat = NaN; % There's no correlation for this stuff... elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase kappat = NaN; % Compressibility is poorly defined for a saturated mixture else % Single Phase kappat = Water_kappat(rho,T); % Directly compute compressibility end elseif (T > Tmin) % Cold Fluid [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase kappat = Water_kappat(rho,T); % Directly compute compressibility else % Ice kappat = NaN; % Would need an appropriate correlation... end elseif (T > Tabsmin) % Cold Vapor Ps = MeltPress(T); % Find the sublimation pressure at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point if (rho <= rhos) % Vapor kappat = Water_kappat(rho,T); % Directly compute compressibility else % Ice Ih kappat = NaN; % Would need the ice correlation to compute this... end end OutVar = kappat; % Output isothermal compressibility %--------------------------------------------------------------------------------------- case 'w_dt' % Sound Velocity from density, temperature rho = Var1; % Input density T = Var2; % Input temperature if (T > Ticemax) % Supercritical Fluid w = Water_w(rho,T); % Directly compute sound velocity elseif (T > Tt) % Ordinary Fluid Pm = MeltPress(T); % Find the melting pressure at this temperature rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature if (rho > rhom) % High-Pressure Ice w = NaN; % There's no correlation for this stuff... elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase w = NaN; % The liquid and vapor components have independent sound velocities else % Single Phase Fluid w = Water_w(rho,T); % Directly compute sound velocity end elseif (T > Tmin) % Cold Fluid [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase w = Water_w(rho,T); % Directly compute sound velocity else % Ice w = NaN; % Would need an appropriate correlation... end elseif (T > Tabsmin) % Cold Vapor Ps = MeltPress(T); % Find the sublimation pressure at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point if (rho <= rhos) % Vapor w = Water_w(rho,T); % Directly compute sound velocity else % Ice Ih w = NaN; % Would need the ice correlation to compute this... end end OutVar = w; % Output sound velocity %--------------------------------------------------------------------------------------- case 'mu_dt' % Viscosity from density, temperature rho = Var1; % Input density T = Var2; % Input temperature if (T > Ticemax) % Supercritical Fluid mu = Water_mu(rho,T); % Directly compute viscosity elseif (T > Tt) % Ordinary Fluid Pm = MeltPress(T); % Find the melting pressure at this temperature rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature if (rho > rhom) % High-Pressure Ice mu = Inf; % Solid elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase mu = NaN; % The liquid and vapor components have independent viscosities else % Single Phase mu = Water_mu(rho,T); % Directly compute viscosity end elseif (T > Tmin) % Cold Fluid [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase mu = Water_mu(rho,T); % Directly compute viscosity else % Ice mu = Inf; % Solid end elseif (T > Tabsmin) % Cold Vapor Ps = MeltPress(T); % Find the sublimation pressure at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point if (rho <= rhos) % Vapor mu = Water_mu(rho,T); % Directly compute viscosity else % Ice Ih mu = Inf; % Solid end end OutVar = mu; % Output viscosity %--------------------------------------------------------------------------------------- case 'k_dt' % Thermal Conductivity from density, temperature rho = Var1; % Input density T = Var2; % Input temperature if (T > Ticemax) % Supercritical Fluid k = Water_k(rho,T); % Directly compute thermal conductivity elseif (T > Tt) % Ordinary Fluid Pm = MeltPress(T); % Find the melting pressure at this temperature rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature if (rho > rhom) % High-Pressure Ice k = NaN; % There's no correlation for this stuff elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase k = NaN; % The liquid and vapor components have independent thermal conductivities else % Single Phase k = Water_k(rho,T); % Directly compute thermal conductivity end elseif (T > Tmin) % Cold Fluid [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase k = Water_k(rho,T); % Directly compute thermal conductivity else % Ice k = NaN; % Would need an appropriate correlation... end elseif (T > Tabsmin) % Cold Vapor Ps = MeltPress(T); % Find the sublimation pressure at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point if (rho <= rhos) % Vapor k = Water_k(rho,T); % Directly compute thermal conductivity else % Ice Ih k = NaN; % Would need an appropriate correlation... end end OutVar = k; % Output thermal conductivity %--------------------------------------------------------------------------------------- case 'sigma_t' % Surface tension from temperature T = Var1; % Input temperature if (T < Tt) || (T > Tc) % Is this temperature off of the saturation line? sigma = NaN; % If so, there's no surface tension else sigma = Water_sigma(T); % Otherwise, directly compute surface tension of the saturated liquid end OutVar = sigma; % Output surface tension %--------------------------------------------------------------------------------------- case 'epsilon_dt' % Dielectric Constant from density, temperature rho = Var1; % Input density T = Var2; % Input temperature if (T > Ticemax) % Supercritical Fluid epsilon = Water_epsilon(rho,T); % Directly compute dielectric constant elseif (T > Tt) % Ordinary Fluid Pm = MeltPress(T); % Find the melting pressure at this temperature rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature if (rho > rhom) % High-Pressure Ice epsilon = NaN; % There's no correlation for this stuff elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase epsilon = NaN; % The liquid and vapor components have independent dielctric constants else % Single Phase epsilon = Water_epsilon(rho,T); % Directly compute dielectric constant end elseif (T > Tmin) % Cold Fluid [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase epsilon = Water_epsilon(rho,T); % Directly compute dielectric constant else % Ice epsilon = NaN; % There's no correlation for this stuff end elseif (T > Tabsmin) % Cold Vapor Ps = MeltPress(T); % Find the sublimation pressure at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point if (rho <= rhos) % Vapor epsilon = Water_epsilon(rho,T); % Directly compute dielectric constant else % Ice Ih epsilon = NaN; % There's no correlation for this stuff end end OutVar = epsilon; % Output dielectric constant %--------------------------------------------------------------------------------------- case 'n_dtl' % Refractive Index from density, temperature, wavelength rho = Var1; % Input density T = Var2; % Input temperature lambda = Var3; % Input wavelength if (T > Ticemax) % Supercritical Fluid n = Water_n(rho, T, lambda); % Directly compute refractive index elseif (T > Tt) % Ordinary Fluid Pm = MeltPress(T); % Find the melting pressure at this temperature rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the melting point [Psat rholsat rhovsat] = GoodSatProp(T); % Find the saturation properties at this temperature if (rho > rhom) % High-Pressure Ice n = NaN; % There's no correlation for this stuff elseif (rho < rholsat) && (rho > rhovsat) % Two-Phase n = NaN; % The liquid and vapor components have independent refractive indices else % Single Phase n = Water_n(rho, T, lambda); % Directly compute refractive index end elseif (T > Tmin) % Cold Fluid [Pm Ps Pf] = MeltPress(T); % Find all of the phase boundary pressures at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point rhom = H2Oprop('d_pt',Pm,T); % Find the liquid density at the ice Ih melting point rhof = H2Oprop('d_pt',Pf,T); % Find the liquid density ath the high-pressure ice melting point if (rho <= rhos) || ((rho >= rhom) && (rho <= rhof)) % Single Fluid Phase n = Water_n(rho, T, lambda); % Directly compute refractive index else % Ice n = NaN; % There's no correlation for this stuff end elseif (T > Tabsmin) % Cold Vapor Ps = MeltPress(T); % Find the sublimation pressure at this temperature rhos = H2Oprop('d_pt',Ps,T); % Find the vapor density at the sublimation point if (rho <= rhos) % Vapor n = Water_n(rho, T, lambda); % Directly compute refractive index else % Ice Ih n = NaN; % There's no correlation for this stuff end end OutVar = n; % Output refractive index %--------------------------------------------------------------------------------------- case 'tsat_p' % Saturation temperature from pressure P = Var1; % Input pressure if (P < Pt) % Are we on the sublimation curve? T = MeltTemp(P); % If so, find the sublimation temperature elseif (P < Pc) % Are we on the saturation curve? T = fzero(@(x) GoodSatProp(x) - P, [Tt Tc - .01]); % If so, use root finding to find the saturation temperature else % Otherwise, there's no vapor phase boundary... T = NaN; % ...so we return NaN end OutVar = T; % Output temperature %--------------------------------------------------------------------------------------- case 'tsat_dl' % Saturation temperature from liquid density, valid above 277.134 K density anomaly rhol = Var1; % Input saturated liquid density if (rhol > rhoc) && (rhol < rhoamax) % Are we on the saturation curve? T = fzero(@(x) GoodSatProp(x,'rholsat') - rhol, [Trhomax Tc - .01]); % If so, use root finding to find the saturation temperature else % Otherwise, the input was not a valid saturated liquid density... T = NaN; % ...so we return NaN end OutVar = T; % Output temperature %--------------------------------------------------------------------------------------- case 'tsat_dv' % Saturation temperature from vapor density rhov = Var1; % Input saturated vapor density if (rhov < rhovt) % Are we on the sublimation curve? T = fzero(@(x) Water_P(rhov, x) - MeltPress(x,'sub'), Tt); % If so, use root finding to find the sublimation temperature elseif (rhov < rhoc) % Are we on the saturation curve? T = fzero(@(x) GoodSatProp(x,'rhovsat') - rhov, [Tt Tc]); % If so, use root finding to find the saturation temperature else % Otherwise, the input was not a valid saturated vapor density... T = NaN; % ...so we return NaN end OutVar = T; % Output temperature %--------------------------------------------------------------------------------------- case 'psat_t' % Saturation pressure from temperature T = Var1; % Input temperature if (T < Tt) % Are we on the sublimation curve? P = MeltPress(T,'sub'); % If so, find the sublimation pressure elseif (T < Tc) % Are we on the saturation curve? P = GoodSatProp(T); % if so, find the saturation pressure else % Otherwise, there's no vapor phase boundary... P = NaN; % ...so we return NaN end OutVar = P; % Output pressure %--------------------------------------------------------------------------------------- case 'dlsat_t' % Saturated liquid density from temperature T = Var1; % Input temperature dl = GoodSatProp(T,'rholsat'); % Find the saturated liquid density OutVar = dl; % Output saturated liquid density %--------------------------------------------------------------------------------------- case 'dvsat_t' % Saturated vapor density from temperature T = Var1; % Input temperature if (T > Tt) % Are we on or above the saturation curve? dv = GoodSatProp(T,'rhovsat'); % If so, find the saturated vapor density else % Otherwise, we must be on the sublimation curve Ps = MeltPress(T,'sub'); % Find the sublimation pressure dv = fzero(@(x) Water_P(x, T) - Ps, rhomin); % Use a root-finding method to find the vapor density at the sublimation point end OutVar = dv; % Output saturated vapor density %--------------------------------------------------------------------------------------- case 'dlsat_p' % Saturated liquid density from pressure P = Var1; % Input pressure if (P < Pt) % Is the pressure too low for liquid to exist? P = NaN; % If so, set it to NaN and let it propagate end T = H2Oprop('tsat_p',P); % Find the saturation temperature if isnan(T) % Is the pressure supercritical, or was it NaN? dl = NaN; % If so, return NaN else dl = GoodSatProp(T,'rholsat'); % Otherwise, calculate the saturated liquid density end OutVar = dl; % Output saturated liquid density %--------------------------------------------------------------------------------------- case 'dvsat_p' % Saturated vapor density from pressure P = Var1; % Input pressure T = H2Oprop('tsat_p',P); % Find the saturation or sublimation temperature if isnan(T) % Does a saturation or sublimation temperature exist? dv = NaN; % If not, return NaN elseif (P < Pt) % Is the process of interest sublimation? dv = fzero(@(x) Water_P(x, T) - P, P/R/T); % If so, use a root-finding method to find the vapor density else % Otherwise, we must be looking at the saturation line dv = GoodSatProp(T,'rhovsat'); % And so, calculate the saturated vapor density end OutVar = dv; % Output saturated vapor density %--------------------------------------------------------------------------------------- case 'ulsat_t' % Saturated liquid specific energy from temperature T = Var1; % Input temperature dl = GoodSatProp(T,'rholsat'); % Find the saturated liquid density if isnan(dl) % Does saturated liquid exist? ul = NaN; % If not, return NaN else ul = Water_u(dl,T); % Otherwise, directly compute the saturated liquid specific internal energy end OutVar = ul; % Output saturated liquid specific internal energy %--------------------------------------------------------------------------------------- case 'ulsat_p' % Saturated liquid specific energy from pressure P = Var1; % Input pressure if (P < Pt) % Is the pressure too low for liquid to exist? P = NaN; % If so, set it to NaN and let it propagate end T = H2Oprop('tsat_p',P); % Find the saturation temperature if isnan(T) % Is the pressure supercritical, or was it NaN? ul = NaN; % If so, return NaN else dl = GoodSatProp(T,'rholsat'); % Otherwise, find the saturated liquid density ul = Water_u(dl,T); % Directly compute the saturated liquid specific internal energy end OutVar = ul; % Output saturated liquid specific internal energy %--------------------------------------------------------------------------------------- case 'uvsat_t' % Saturated vapor specific energy from temperature T = Var1; % Input temperature if (T < Tt) % Are we on the sublimation curve? Ps = MeltPress(T,'sub'); % If so, find the sublimation pressure dv = fzero(@(x) Water_P(x, T) - Ps, rhomin); % Use a root-finding method to compute the saturated vapor density uv = Water_u(dv,T); % Directly compute the saturated vapor specific internal energy elseif (T > Tc) % Are we supercritical? uv = NaN; % If so, return NaN else dv = GoodSatProp(T,'rhovsat'); % Otherwise, find the saturated vapor density uv = Water_u(dv,T); % Directly compute the saturated vapor specific internal energy end OutVar = uv; % Output saturated vapor specific internal energy %--------------------------------------------------------------------------------------- case 'uvsat_p' % Saturated vapor specific energy from pressure P = Var1; % Input pressure T = H2Oprop('tsat_p',P); % Find the saturation or sublimation temperature if isnan(T) % Does a saturation or sublimation temperature exist? uv = NaN; % If not, return NaN elseif (P < Pt) % Is the process of interest sublimation? dv = fzero(@(x) Water_P(x, T) - P, P/R/T); % If so, use a root-finding method to find the vapor density uv = Water_u(dv,T); % Directly compute the saturated vapor specific internal energy else dv = GoodSatProp(T,'rhovsat'); % Otherwise, find the saturated vapor density uv = Water_u(dv,T); % Directly compute the saturated vapor specific internal energy end OutVar = uv; % Output saturated vapor specific internal energy %--------------------------------------------------------------------------------------- case 'hlsat_t' % Saturated liquid specific enthalpy from temperature T = Var1; % Input temperature dl = GoodSatProp(T,'rholsat'); % Find the saturated liquid density if isnan(dl) % Does saturated liquid exist? hl = NaN; % If not, return NaN else hl = Water_h(dl,T); % Otherwise, directly compute the saturated liquid specific enthalpy end OutVar = hl; % Output saturated liquid specific enthalpy %--------------------------------------------------------------------------------------- case 'hlsat_p' % Saturated liquid specific enthalpy from pressure P = Var1; % Input pressure if (P < Pt) % Is the pressure too low for liquid to exist? P = NaN; % If so, set it to NaN and let it propagate end T = H2Oprop('tsat_p',P); % Find the saturation temperature if isnan(T) % Is the pressure supercritical, or was it NaN? hl = NaN; % If so, return NaN else dl = GoodSatProp(T,'rholsat'); % Otherwise, find the saturated liquid density hl = Water_h(dl,T); % Directly compute the saturated liquid specific enthalpy end OutVar = hl; % Output saturated liquid specific enthalpy %--------------------------------------------------------------------------------------- case 'hvsat_t' % Saturated vapor specific enthalpy from temperature T = Var1; % Input temperature if (T < Tt) % Are we on the sublimation curve? Ps = MeltPress(T,'sub'); % If so, find the sublimation pressure dv = fzero(@(x) Water_P(x, T) - Ps, rhomin); % Use a root-finding method to compute the saturated vapor density hv = Water_h(dv,T); % Directly compute the saturated vapor specific internal energy elseif (T > Tc) % Are we supercritical? hv = NaN; % If so, return NaN else dv = GoodSatProp(T,'rhovsat'); % Otherwise, find the saturated vapor density hv = Water_h(dv,T); % Directly compute the saturated vapor specific internal energy end OutVar = hv; % Output saturated vapor specific internal energy %--------------------------------------------------------------------------------------- case 'hvsat_p' % Saturated vapor specific enthalpy from pressure P = Var1; % Input pressure T = H2Oprop('tsat_p',P); % Find the saturation or sublimation temperature if isnan(T) % Does a saturation or sublimation temperature exist? hv = NaN; % If not, return NaN elseif (P < Pt) % Is the process of interest sublimation? dv = fzero(@(x) Water_P(x, T) - P, P/R/T); % If so, use a root-finding method to find the vapor density hv = Water_h(dv,T); % Directly compute the saturated vapor specific enthalpy else dv = GoodSatProp(T,'rhovsat'); % Otherwise, find the saturated vapor density hv = Water_h(dv,T); % Directly compute the saturated vapor specific enthalpy end OutVar = hv; % Output saturated vapor specific enthalpy %--------------------------------------------------------------------------------------- case 'slsat_t' % Saturated liquid specific entropy from temperature T = Var1; % Input temperature dl = GoodSatProp(T,'rholsat'); % Find the saturated liquid density if isnan(dl) % Does saturated liquid exist? sl = NaN; % If not, return NaN else sl = Water_s(dl,T); % Otherwise, directly compute the saturated liquid specific entropy end OutVar = sl; % Output saturated liquid specific entropy %--------------------------------------------------------------------------------------- case 'slsat_p' % Saturated liquid specific entropy from pressure P = Var1; % Input pressure if (P < Pt) % Is the pressure too low for liquid to exist? P = NaN; % If so, set it to NaN and let it propagate end T = H2Oprop('tsat_p',P); % Find the saturation temperature if isnan(T) % Is the pressure supercritical, or was it NaN? sl = NaN; % If so, return NaN else dl = GoodSatProp(T,'rholsat'); % Otherwise, find the saturated liquid density sl = Water_s(dl,T); % Directly compute the saturated liquid specific entropy end OutVar = sl; % Output saturated liquid specific entropy %--------------------------------------------------------------------------------------- case 'svsat_t' % Saturated vapor specific entropy from temperature T = Var1; % Input temperature dv = H2Oprop('dvsat_t',T); % Find the saturated vapor density if (T > Tc) % If the temperature is supercritical... sv = NaN; % ...return NaN else sv = Water_s(dv,T); % Otherwise, directly compute the saturated vapor specific entropy end OutVar = sv; % Output saturated vapor specific entropy %--------------------------------------------------------------------------------------- case 'svsat_p' % Saturated vapor specific entropy from pressure P = Var1; % Input pressure T = H2Oprop('tsat_p',P); % Find the saturation temperature dv = H2Oprop('dvsat_p',P); % Find the saturated vapor density if isnan(T) % If a saturation or sublimation temperature could not be found... sv = NaN; % ...return NaN else sv = Water_s(dv,T); % Otherwise, directly compute the saturated vapor specific entropy end OutVar = sv; % Output saturated vapor specific entropy %--------------------------------------------------------------------------------------- otherwise help H2Oprop % Display help text to assist with choice of valid query error('Invalid input string.') % Throw an exception to stop any code that might have called this function end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Internal Helper Functions function P = Water_P(rho, T) % Compute pressure using the equation of state global Tc rhoc R % Inherit a few basic physical constants tau = Tc/T; % Compute the reduced temperature delta = rho/rhoc; % Compute the reduced density phir_d = phirx_d(delta,tau); % Compute the partial derivative of the residual part of Helmholtz energy with respect to density P = rho*R*T*(1 + delta*phir_d); % Calculate the pressure %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function u = Water_u(rho, T) % Compute specific internal energy using the equation of state global Tc rhoc R n0 gamma0 % Inherit a few basic physical constants and equation of state parameters tau = Tc/T; % Compute the reduced temperature delta = rho/rhoc; % Compute the reduced density phi0_t = n0(2) + n0(3)/tau + sum(n0(4:8).*gamma0(4:8).*((1 -... exp(-gamma0(4:8)*tau)).^(-1) - 1)); % Compute the partial derivative of the ideal gas part of Helmholtz energy with respect to temperature phir_t = phirx_t(delta,tau); % Compute the partial derivative of the residual part of Helmholtz energy with respect to temperature u = R*T*tau*(phi0_t + phir_t); % Calculate the specific internal energy %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function s = Water_s(rho, T) % Compute specific entropy using the equation of state global Tc rhoc R n0 gamma0 % Inherit a few basic physical constants and equation of state parameters tau = Tc/T; % Compute the reduced temperature delta = rho/rhoc; % Compute the reduced density phi0 = log(delta) + n0(1) + n0(2)*tau + n0(3)*log(tau) + sum(n0(4:8).*log(1-... exp(-gamma0(4:8)*tau))); % Compute the ideal gas part of Helmholtz energy phi0_t = n0(2) + n0(3)/tau + sum(n0(4:8).*gamma0(4:8).*((1 -... exp(-gamma0(4:8)*tau)).^(-1) - 1)); % Compute the partial derivative of the ideal gas part of Helmholtz energy with respect to temperature phir = phirx(delta,tau); % Compute the residual part of Helmholtz energy phir_t = phirx_t(delta,tau); % Compute the partial derivative of the residual part of Helmholtz energy with respect to temperature s = R*tau*(phi0_t + phir_t) - R*phi0 - R*phir; % Calculate the specific entropy %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function h = Water_h(rho, T) % Compute specific enthalpy using the equation of state global Tc rhoc R n0 gamma0 % Inherit a few basic physical constants and equation of state parameters tau = Tc/T; % Compute the reduced temperature delta = rho/rhoc; % Compute the reduced density phi0_t = n0(2) + n0(3)/tau + sum(n0(4:8).*gamma0(4:8).*((1 -... exp(-gamma0(4:8)*tau)).^(-1) - 1)); % Compute the partial derivative of the ideal gas part of Helmholtz energy with respect to temperature phir_d = phirx_d(delta,tau); % Compute the partial derivative of the residual part of Helmholtz energy with respect to density phir_t = phirx_t(delta,tau); % Compute the partial derivative of the residual part of Helmholtz energy with respect to temperature h = R*T*(1 + tau*(phi0_t + phir_t) + delta*phir_d); % Calculate the specific enthalpy %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function cp = Water_cp(rho, T) % Compute the constant-pressure specific heat using the equation of state global Tc rhoc R n0 gamma0 % Inherit a few basic physical constants and equation of state parameters tau = Tc/T; % Compute the reduced temperature delta = rho/rhoc; % Compute the reduced density phi0_tt = -n0(3)/tau^2 - sum(n0(4:8).*gamma0(4:8).^2.*exp(-gamma0(4:8)*tau).*... (1 - exp(-gamma0(4:8)*tau)).^(-2)); % Compute the second partial derivative of the ideal gas part of Helmholtz energy with respect to temperature phir_d = phirx_d(delta,tau); % Compute the partial derivative of the residual part of Helmholtz energy with respect to density phir_dd = phirx_dd(delta,tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density phir_tt = phirx_tt(delta,tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to temperature phir_dt = phirx_dt(delta,tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density and temperature cp = -R*tau^2*(phi0_tt + phir_tt) + R*(1 + delta*phir_d - delta*tau*phir_dt)^2/... (1 + 2*delta*phir_d + delta^2*phir_dd); % Calculate the constant-pressure specific heat %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function betat = Water_betat(rho, T) % Compute the thermal expansion coefficient using the equation of state global Tc rhoc % Inherit a few basic physical constants tau = Tc/T; % Compute the reduced temperature delta = rho/rhoc; % Compute the reduced density phir_d = phirx_d(delta,tau); % Compute the partial derivative of the residual part of Helmholtz energy with respect to density phir_dd = phirx_dd(delta,tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density phir_dt = phirx_dt(delta,tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density and temperature betat = 1/T*(1+delta*phir_d-delta*tau*phir_dt)/(1+2*delta*phir_d+delta^2*phir_dd); % Calculate the thermal expansion coefficient %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function cv = Water_cv(rho, T) % Compute the constant-volume specific heat using the equation of state global Tc rhoc R n0 gamma0 % Inherit a few basic physical constants and equation of state parameters tau = Tc/T; % Compute the reduced temperature delta = rho/rhoc; % Compute the reduced density phi0_tt = -n0(3)/tau^2 - sum(n0(4:8).*gamma0(4:8).^2.*exp(-gamma0(4:8)*tau).*... (1 - exp(-gamma0(4:8)*tau)).^(-2)); % Compute the second partial derivative of the ideal gas part of Helmholtz energy with respect to temperature phir_tt = phirx_tt(delta,tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to temperature cv = -R*tau^2*(phi0_tt + phir_tt); % Calculate the constant-volume specific heat %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function kappat = Water_kappat(rho, T) % Compute the isothermal compressibility using the equation of state global Tc rhoc R n0 gamma0 % Inherit a few basic physical constants and equation of state parameters tau = Tc/T; % Compute the reduced temperature delta = rho/rhoc; % Compute the reduced density phi0_tt = -n0(3)/tau^2 - sum(n0(4:8).*gamma0(4:8).^2.*exp(-gamma0(4:8)*tau).*... (1 - exp(-gamma0(4:8)*tau)).^(-2)); % Compute the second partial derivative of the ideal gas part of Helmholtz energy with respect to temperature phir_d = phirx_d(delta,tau); % Compute the partial derivative of the residual part of Helmholtz energy with respect to density phir_dd = phirx_dd(delta,tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density phir_tt = phirx_tt(delta,tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to temperature phir_dt = phirx_dt(delta,tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density and temperature % Properties cv = -R*tau^2*(phi0_tt + phir_tt); % Isochoric Heat Capacity cp = -R*tau^2*(phi0_tt + phir_tt) + R*(1 + delta*phir_d - delta*tau*phir_dt)^2/... (1 + 2*delta*phir_d + delta^2*phir_dd); % Isobaric Heat Capacity w = sqrt(R*T*(1 + 2*delta*phir_d + delta^2*phir_dd - (1 + delta*phir_d - ... delta*tau*phir_dt)^2/tau^2/(phi0_tt + phir_tt))); % Speed of Sound kappat = cp/cv/rho/w^2; % Isothermal Compressibility %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function w = Water_w(rho, T) % Compute the speed of sound using the equation of state global Tc rhoc R n0 gamma0 % Inherit a few basic physical constants and equation of state parameters tau = Tc/T; % Compute the reduced temperature delta = rho/rhoc; % Compute the reduced density phi0_tt = -n0(3)/tau^2 - sum(n0(4:8).*gamma0(4:8).^2.*exp(-gamma0(4:8)*tau).*... (1 - exp(-gamma0(4:8)*tau)).^(-2)); % Compute the second partial derivative of the ideal gas part of Helmholtz energy with respect to temperature phir_d = phirx_d(delta,tau); % Compute the partial derivative of the residual part of Helmholtz energy with respect to density phir_dd = phirx_dd(delta,tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density phir_tt = phirx_tt(delta,tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to temperature phir_dt = phirx_dt(delta,tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density and temperature w = sqrt(R*T*(1 + 2*delta*phir_d + delta^2*phir_dd - (1 + delta*phir_d - ... delta*tau*phir_dt)^2/tau^2/(phi0_tt + phir_tt))); % Speed of Sound %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function sigma = Water_sigma(T) % Compute the surface tension using the IAPWS 1994 formulation global Tc % Inherit a few basic physical constants TAU = 1 - T/Tc; % Compute the (alternate) reduced temperature sigma = 0.2358*TAU^1.256*(1 - 0.625*TAU); % Use the IAPWS correlation to compute surface tension %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function mu = Water_mu(rho, T) % Compute the viscosity using the IAPWS 2008 formulation % *** NOTE: SOME DISAGREEMENT WITH DATA NEAR CRITICAL POINT % Probably a subtle error in isothermal compressibility, shared with % thermal conductivity formulation global Tc rhoc R % Inherit a few basic physical constants Tstar = 647.096; % Reference Temperature rhostar = 322.0; % Reference Density Pstar = 22.064e6; % Reference Pressure mustar = 1e-6; % Reference Viscosity % Viscosity correlation constants iv = [0 1 2 3 0 1 2 3 5 0 1 2 3 4 0 1 0 3 4 3 5]'; jv = [0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 4 4 5 6 6]'; Nv = ... [ 5.20094e-1 8.50895e-2 -1.08374 -2.89555e-1 2.22531e-1 9.99115e-1 ... 1.88797 1.26613 1.20573e-1 -2.81378e-1 -9.06851e-1 -7.72479e-1 ... -4.89837e-1 -2.57040e-1 1.61913e-1 2.57399e-1 -3.25372e-2 6.98452e-2 ... 8.72102e-3 -4.35673e-3 -5.93264e-4]'; xi0 = 0.13e-9; nu = 0.630; gamma0 = 1.239; Gamma0 = 0.06; xmu = 0.068; % Critical exponent for viscosity qC = 1/1.9e-9; qD = 1/1.1e-9; xi1 = 0.3817016416e-9; % Viscosity correlation computation Tbar = T/Tstar; % Compute the (alternative) reduced temperature Tbar1 = 1.5; % Define a reference reduced temperature T1 = Tbar1*Tstar; % Compute the reference absolute temperature rhobar = rho/rhostar; % Compute the (alternative) reduced density tau = Tc/T; % Compute the (conventional) reduced temperature tau1 = Tc/T1; % Compute the reference (conventional) reduced temperature delta = rho/rhoc; % Compute the (conventional) reduced density if (Tbar <= 1.01) && (Tbar >= 0.996) && (rhobar <= 1.36) && (rhobar >= 0.31) % Near the critical point? If so, do the expensive critical enhancement calculation phir_d = phirx_d(delta, tau); % Compute the partial derivative of the residual part of Helmholtz energy with respect to density phir_d1 = phirx_d(delta, tau1); % Compute the partial derivative of the residual part of Helmholtz energy with respect to density at the reference temperature phir_dd = phirx_dd(delta, tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density phir_dd1 = phirx_dd(delta, tau1); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density at the reference temperature chibar = rhobar*Pstar/rhostar/R/T/(1 + 2*delta*phir_d + delta^2*phir_dd); % Compute the compressibility chibar1 = rhobar*Tbar1/Tbar * ... Pstar/rhostar/R/T1/(1 + 2*delta*phir_d1 + delta^2*phir_dd1); % Compute the compressibility at the reference temperature dchibar = max((chibar - chibar1), 0); % Compute the compressibility difference with the reference temperature and discard mu2 if it's negative xi = xi0*(dchibar/Gamma0)^(nu/gamma0); % Compute the correlation length if xi <= xi1 % Use the below equation for Y if the correlation length is very short Y = .2*qC*xi*(qD*xi)^5*(1 - qC*xi + (qC*xi)^2 - 765/504*(qD*xi)^2); else % Otherwise, make the more complex calculation below for Y psiD = acos((1 + (qD*xi)^2)^-.5); w = abs((qC*xi - 1)/(qC*xi + 1))^.5 * tan(psiD/2); L = (qC*xi <= 1) * 2*atan(abs(w)) + (qC*xi > 1) * log((1 + w)/(1 - w)); Y = sin(3*psiD)/12 - sin(2*psiD)/4/qC/xi + ... (1 - 1.25*(qC*xi)^2)/(qC*xi)^2*sin(psiD) - ... ((1 - 1.5*(qC*xi)^2)*psiD - abs((qC*xi)^2 - 1)^1.5*L)/(qC*xi)^3; end mu2 = exp(xmu*Y); % Compute the critical enhancement for viscosity else mu2 = 1; % Neglect any critical enhancement if far from the critical point end mu0 = 100*sqrt(Tbar)/(1.67752 + 2.20462/Tbar + 0.6366564/Tbar^2 - 0.241605/Tbar^3); % Compute dilute gas viscosity mu1 = exp(rhobar*sum(Nv.*(1/Tbar - 1).^iv.*(rhobar - 1).^jv)); % Compute finite density viscosity enhancement mu = mu0*mu1*mu2*mustar; % Use the IAPWS 2008 formulation to compute viscosity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function k = Water_k(rho, T) % Compute the thermal conductivity using the IAPWS 2008 formulation % *** NOTE: SOME DISAGREEMENT WITH DATA NEAR CRITICAL POINT % Probably a subtle error in isothermal compressibility, shared with % viscosity formulation global Tc rhoc R % Inherit a few basic physical constants Tstar = 647.226; % Reference Temperature rhostar = 317.763; % Reference Density Pstar = 22.115e6; % Reference Pressure kstar = 0.4945; % Reference Thermal Conductivity % Thermal conductivity correlation constants ik = [0 0 0 0 0 0 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 4 4]'; jk = [0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 0 1]'; Nk = ... [ 1.3293046 -0.40452437 0.24409490 0.018660751 -0.12961068 0.044809953 ... 1.7018363 -2.2156845 1.6511057 -0.76736002 0.37283344 -0.11203160 ... 5.2246158 -10.124111 4.9874687 -0.27297694 -0.43083393 0.13333849 ... 8.7127675 -9.5000611 4.3786606 -0.91783782 -1.8525999 0.93404690]'; iv = [0 1 2 3 0 1 2 3 5 0 1 2 3 4 0 1 0 3 4 3 5]'; jv = [0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 4 4 5 6 6]'; Nv = ... [ 5.20094e-1 8.50895e-2 -1.08384 -2.89555e-1 2.22531e-1 9.99115e-1 ... 1.88797 1.26613 1.20573e-1 -2.81378e-1 -9.06851e-1 -7.72479e-1 ... -4.89837e-1 -2.57040e-1 1.61913e-1 2.57399e-1 -3.25372e-2 6.98452e-2 ... 8.72102e-3 -4.35673e-3 -5.93264e-4]'; % Thermal conductivity correlation computation Tbar = T/Tstar; % Compute the (alternative) reduced temperature rhobar = rho/rhostar; % Compute the (alternative) reduced density tau = Tc/T; % Compute the (conventional) reduced temperature delta = rho/rhoc; % Compute the (conventional) reduced density phir_d = phirx_d(delta, tau); % Compute the partial derivative of the residual part of Helmholtz energy with respect to density phir_dd = phirx_dd(delta, tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density phir_dt = phirx_dt(delta, tau); % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density and temperature chibar = rhobar*Pstar/rhostar/R/T/(1 + 2*delta*phir_d + delta^2*phir_dd); % Compute the compressibility dPbar = rho*R*Tstar/Pstar*(1 + delta*phir_d - delta*tau*phir_dt); % Compute the isochoric pressure-temperature coefficient mu0 = 100*sqrt(Tbar)/(1.67752 + 2.20462/Tbar + 0.6366564/Tbar^2 - 0.241605/Tbar^3); % Compute dilute gas viscosity mu1 = exp(rhobar*sum(Nv.*(1/Tbar - 1).^iv.*(rhobar - 1).^jv)); % Compute finite density viscosity enhancement k0 = sqrt(Tbar)/(1 + 6.978267/Tbar + 2.599096/Tbar^2 - 0.998254/Tbar^3); % Compute dilute gas thermal conductivity k1 = exp(rhobar*sum(Nk.*(1/Tbar - 1).^ik.*(rhobar - 1).^jk)); % Compute finite density thermal conductivity enhancement k2 = 55.071*0.0013848/mu0/mu1*(Tbar/rhobar)^2*dPbar^2*chibar^0.4678*rhobar^0.5 * ... exp(-18.66*(Tbar - 1)^2 - (rhobar - 1)^4); % Compute critical enhancement in thermal conductivity k = (k0*k1 + k2)*kstar; % Use the IAPWS 2008 formulation to compute thermal conductivity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function epsilon = Water_epsilon(rho, T) % Compute the dielectric constant using the IAPWS 1997 formulation global Tc rhoc % Inherit a few basic physical constants % Dielectric constant correlation constants Nh = ... [ 0.978224486826 -0.957771379375 0.237511794148 0.714692244396 ... -0.298217036956 -0.108863472196 0.949327488264e-1 -0.980469816509e-2 ... 0.165167634970e-4 0.937359795772e-4 -0.123179218720e-9 0.196096504426e-2]'; ih = [1 1 1 2 3 3 4 5 6 7 10]'; jh = [0.25 1 2.5 1.5 1.5 2.5 2 2 5 0.5 10]'; Na = 6.0221367e23; epsilon0 = 1/(4e-7*pi*299792458^2); alpha = 1.636e-40; Mw = 0.018015268; mu = 6.138e-30; kB = 1.380658e-23; % Dielectric constant correlation computation tau = Tc/T; % Compute the reduced temperature delta = rho/rhoc; % Compute the reduced density g = 1 + sum(Nh(1:11).*delta.^ih.*tau.^jh) + Nh(12)*delta*(T/228 - 1)^-1.2; % Compute the Harris-Alder g-factor Ah = Na*mu^2/epsilon0/kB * rho*g/T/Mw; Bh = Na*alpha/epsilon0/3 * rho/Mw; epsilon = (1 + Ah + 5*Bh + sqrt(9 + 2*Ah + 18*Bh + Ah^2 + 10*Ah*Bh + 9*Bh^2))/ ... (4 - 4*Bh); % Use the IAPWS 1997 formulation to compute the dielectric constant %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function n = Water_n(rho, T, lambda) % Compute the refractive index using the IAPWS 1997 formulation Tbar = T/273.15; % Calculate the reduced temperature rhobar = rho/1000; % Calculate the reduced density lambdabar = lambda/589e-9; % Calculate the reduced wavelength nx = rhobar*(0.244257733 + 9.74634476e-3*rhobar - 3.73234996e-3*Tbar + ... 2.68678472e-4*lambdabar^2*Tbar + 1.58920570e-3/lambdabar^2 + ... 2.45934259e-3/(lambdabar^2-0.2292020^2) + ... 0.900704920/(lambdabar^2-5.432937^2) - 1.66626219e-2*rhobar^2); n = sqrt((1 + 2*nx)/(1 - nx)); % Use the IAPWS 1997 formulation to compute the refractive index %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Phase Boundary Functions function [Psat, rholsat, rhovsat] = GoodSatProp(T, prop) % Use the Maxwell relation to find saturation properties from the equation of state global Tc rhoc R % Inherit a few basic physical constants tau = Tc/T; % Compute the reduced temperature % Use the 1992 IAPWS correlation for saturation properties to determine a % starting guess for the saturation properties [Psat rholsat rhovsat] = SatProp(T); % Compute the approximate saturation properties if isnan(Psat) || ((Tc - T) < 1e-5) % If the given temperature is not on the saturation line, or is very close to the critical point, keep the approximate values return; end % Refine the initial estimate by iterating the Maxwell relation tol1 = 1e-10; % Tolerance on iteration-to-iteration property change tol2 = 1e-12; % Tolerance on change in iteration-to-iteration property change err = 10; % Initial iteration-to-iteration property change err1 = 20; % Initial previous iteration-to-iteration property change kmax = 20; % Maximum iteration count k = 0; % Starting iteration while (err > tol1) && (k < kmax) && (abs(err - err1) > tol2) rholsat1 = fzero(@(x) (Psat - R*T*x*(1 + x/rhoc*phirx_d(x/rhoc,tau))), rholsat); % Find best saturated liquid density given Psat if (rholsat1 < rhoc) % Found wrong zero, which is possible near the critical point Guess = R*T*rhoc*(1 + phirx_d(1,tau)); % Use the directly calculated pressure for critical density and the current temperature in place of Psat rholsat1 = fzero(@(x) (Guess - R*T*x*(1 + x/rhoc*phirx_d(x/rhoc,tau))), ... [322.001 350]); % Use firm bounds to ensure that the correct zero is located end rhovsat1 = fzero(@(x) (Psat - R*T*x*(1 + x/rhoc*phirx_d(x/rhoc,tau))), rhovsat); % Find best saturated vapor density given Psat if (rhovsat1 > rhoc) % Found wrong zero, which is possible near the critical point Guess = R*T*rhoc*(1 + phirx_d(1,tau)); % Use the directly calculated pressure for critical density and the current temperature in place of Psat rhovsat1 = fzero(@(x) (Guess - R*T*x*(1 + x/rhoc*phirx_d(x/rhoc,tau))), ... [300 321.999]); % Use firm bounds to ensure that the correct zero is located end Psat1 = R*T/(1/rhovsat1 - 1/rholsat1)*(log(rholsat1/rhovsat1) + ... phirx(rholsat1/rhoc,tau) - phirx(rhovsat1/rhoc,tau)); % Use the new estimates of the saturation densities to improve the estimate of Psat errmat = ([Psat rholsat rhovsat] - [Psat1 rholsat1 rhovsat1])./... [Psat rholsat rhovsat]; % Construct a matrix of the relative change in saturation properties err1 = err; % Store the previous iteration-to-iteration property change err = max(abs(errmat)); % Note the largest relative iteration-to-iteration property change Psat = Psat1; % Update the saturation pressure with the new estimate rholsat = rholsat1; % Update the saturated liquid density with the new estimate rhovsat = rhovsat1; % Update the saturated vapor density with the new estimate k = k + 1; % Increment the iteration counter end if (nargin > 1) % If a specific property was requested in the first function output, put it there switch prop case 'rholsat' Psat = rholsat; % Return the saturated liquid density first case 'rhovsat' Psat = rhovsat; % Return the saturated vapor density first end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [Psat rholsat rhovsat] = SatProp(T, prop) % Estimate the saturation properties using the 1992 IAPWS supplementary release global Tc Tt Pc rhoc % Inherit a few basic physical constants % Saturation property correlation constants asat = [-7.85951783 1.84408259 -11.7866497 22.6807411 -15.9618719 1.80122502]'; Asat = [1 1.5 3 3.5 4 7.5]'; bsat = [1.99274064 1.09965342 -0.510839303 -1.75493479 -45.5170352 -6.74694450e5]'; Bsat = [1/3 2/3 5/3 16/3 43/3 110/3]'; csat = [-2.03150240 -2.68302940 -5.38626492 -17.2991605 -44.7586581 -63.9201063]'; Csat = [2/6 4/6 8/6 18/6 37/6 71/6]'; % Saturation propery computation TAU = 1 - T/Tc; % Calculate the reduced temperature if (T > Tc) || (T < Tt) % If we aren't on the saturation line, return NaN Psat = NaN; rholsat = NaN; rhovsat = NaN; else % Otherwise, use the correlations to calculate the saturation properties Psat = Pc*exp(Tc/T*sum(asat.*TAU.^Asat)); rholsat = rhoc*(1 + sum(bsat.*TAU.^Bsat)); rhovsat = rhoc*exp(sum(csat.*TAU.^Csat)); end if (nargin > 1) % If a specific property was requested in the first function output, put it there switch prop case 'rholsat' Psat = rholsat; % Return the saturated liquid density first case 'rhovsat' Psat = rhovsat; % Return the saturated vapor density first end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function T = MeltTemp(P) % 2008 IAPWS Formulation % Describes melting curve from 273.16 K to 715 K % Describes sublimation curve from 50 K to 273.16 K % Triple Point Constants P1 = 611.657; %Ice I - Liquid - Vapor T1 = 273.16; %Ice I - Liquid - Vapor P2 = 208.566566032976716757e6; %Ice I - Ice III - Liquid (Adjusted to match correlations; nominally 208.566 MPa) T2 = 251.165; %Ice I - Ice III - Liquid P3 = 350.1e6; %Ice III - Ice V - Liquid T3 = 256.164; %Ice III - Ice V - Liquid P4 = 632.4e6; %Ice V - Ice VI - Liquid T4 = 273.31; %Ice V - Ice VI - Liquid P5 = 2216e6; %Ice VI - Ice VII - Liquid T5 = 355; %Ice VI - Ice VII - Liquid Tmax = 715; %Limits of Validity Tmin = 50; Pmax = 20617812820.4493; % Highest pressure on melting line Pmin = 1.9349587189864e-040; % Lowest pressure on melting line % Computation if (P < Pmin) %Below valid range T = NaN; elseif (P < P1) %Sublimation fun = @(x) P1*exp(T1/x*(-0.212144006e2*(x/T1)^(1/300) + ... 0.273203819e2*(x/T1)^(181/150) - ... 0.610598130e1*(x/T1)^(511/300))); T = fzero(@(x) fun(x) - P, [Tmin T1]); % Use a root-finding technique to find the sublimation temperature elseif (P < P2) %Ice Ih fun = @(x) P1*(1 + 0.119539337e7*(1 - (x/T1)^3) + ... 0.808183159e5*(1 - (x/T1)^25.75) + ... 0.333826860e4*(1 - (x/T1)^103.75)); T = fzero(@(x) fun(x) - P, [T2 T1]); % Use a root-finding technique to find the Ice Ih melting point elseif (P < P3) %Ice III T = T2*((-0.700052 + P/P2)/0.299948)^(1/60); % Invert the correlation explicitly to find the Ice III melting point elseif (P < P4) %Ice V T = T3*((0.18721 + P/P3)/1.18721)^(1/8); % Invert the correlation explicitly to find the Ice V melting point elseif (P < P5) %Ice VI T = T4*((0.07476 + P/P4)/1.07476)^(1/4.6); % Invert the correlation explicitly to find the Ice VI melting point elseif (P < Pmax) %Ice VII fun = @(x) P5*exp(1.73683*(1-T5/x) - 0.0544606*(1-(x/T5)^5) + ... 8.06106e-8*(1-(x/T5)^22)); T = fzero(@(x) fun(x) - P, [T5 Tmax]); % Use a root-finding technique to find the Ice VII melting point else %Above valid range T = NaN; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [Pm Pm1 Pm2] = MeltPress(T, prop) % 2008 IAPWS Formulation % Gives sublimation curve from 130 K to 251.165 K % Gives Ice Ih curve from 251.165 K to 273.16 K (Plus other phases) % Gives Ice V curve from 273.16 K to 273.31 K % Gives Ice VI curve from 273.31 K to 355 K % Gives Ice VII curve from 355 K to 715 K % From 251.165 K to 273.16 K, Pm is the Ice Ih pressure, Pm1 is the % sublimation pressure, and Pm2 is the Ice III or Ice V pressure. This can % be overridden by supplying the 'prop' argument, wthich can have one of % two values: 'sub' returns sublimation pressure in Pm, and 'hp' returns % the high-pressure ice melting mressure in Pm. The 'prop' argument is % useful in some minimization and root-finding problems, where only the % first output parameter can be used. % Triple Point Constants P1 = 611.657; % Ice I - Liquid - Vapor T1 = 273.16; % Ice I - Liquid - Vapor P2 = 208.566566032976716757e6; % Ice I - Ice III - Liquid (Adjusted to match correlations; nominally 208.566 MPa) T2 = 251.165; % Ice I - Ice III - Liquid P3 = 350.1e6; % Ice III - Ice V - Liquid T3 = 256.164; % Ice III - Ice V - Liquid P4 = 632.4e6; % Ice V - Ice VI - Liquid T4 = 273.31; % Ice V - Ice VI - Liquid P5 = 2216e6; % Ice VI - Ice VII - Liquid T5 = 355; % Ice VI - Ice VII - Liquid Tabsmin = 50; % Lowest Sublimation Temperature Ticemax = 715; % Highest Freezing Temperature % Computation if (T < Tabsmin) % No Data Pm = NaN; Pm1 = NaN; Pm2 = NaN; elseif (T < T2) % Sublimation Only Pm = P1*exp(T1/T*(-0.212144006e2*(T/T1)^(1/300) + 0.273203819e2*(T/T1)^(181/150) - ... 0.610598130e1*(T/T1)^(511/300))); % Report sublimation pressure in Pm Pm1 = NaN; Pm2 = NaN; elseif (T < T3) % Ice Ih + Sublimation + Ice III Pm = P1*(1 + 0.119539337e7*(1 - (T/T1)^3) + 0.808183159e5*(1 - (T/T1)^25.75) + ... 0.333826860e4*(1 - (T/T1)^103.75)); % Report Ice Ih pressure in Pm Pm1 = P1*exp(T1/T*(-0.212144006e2*(T/T1)^(1/300) + ... 0.273203819e2*(T/T1)^(181/150) - ... 0.610598130e1*(T/T1)^(511/300))); % Report sublimation pressure in Pm1 Pm2 = P2*(1 - 0.299948*(1 - (T/T2)^60)); % Report Ice III pressure in Pm2 elseif (T <= T1) % Ice Ih + Sublimation + Ice V Pm = P1*(1 + 0.119539337e7*(1 - (T/T1)^3) + 0.808183159e5*(1 - (T/T1)^25.75) + ... 0.333826860e4*(1 - (T/T1)^103.75)); % Report Ice Ih pressure in Pm Pm1 = P1*exp(T1/T*(-0.212144006e2*(T/T1)^(1/300) + ... 0.273203819e2*(T/T1)^(181/150) - ... 0.610598130e1*(T/T1)^(511/300))); % Report sublimation pressure in Pm1 Pm2 = P3*(1 - 1.18721*(1 - (T/T3)^8)); % Report Ice V pressure in Pm2 elseif (T < T4) % Ice V Pm = P3*(1 - 1.18721*(1 - (T/T3)^8)); % Report Ice V pressure in Pm Pm1 = NaN; Pm2 = NaN; elseif (T < T5) % Ice VI Pm = P4*(1 - 1.07476*(1 - (T/T4)^4.6)); % Report Ice VI pressure in Pm Pm1 = NaN; Pm2 = NaN; elseif (T <= Ticemax) % Ice VII Pm = P5*exp(1.73683*(1 - T5/T) - 0.0544606*(1 - (T/T5)^5) + ... 8.06106e-8*(1 - (T/T5)^22)); % Report Ice VII pressure in Pm Pm1 = NaN; Pm2 = NaN; else % Temperature above experimental data; melting line unknown Pm = NaN; Pm1 = NaN; Pm2 = NaN; end if (nargin > 1) && ~isnan(Pm1) % Only listen to 'prop' parameter if it's relevant switch prop case 'sub' % Report sublimation pressure in Pm Pm = Pm1; case 'hp' % Report high-pressure ice melting pressure in Pm Pm = Pm2; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Low-level Equation of State Functions % Consult the IAPWS 1995 release for details of these functions function Q = phirx(delta, tau) % Compute the residual part of Helmholtz energy global c d t a b B C D A alpha beta gamma epsilon n psi = exp(-C*(delta - 1)^2 - D*(tau - 1)^2); theta = 1 - tau + A.*((delta - 1)^2).^(1/2./beta); DELTA = theta.^2 + B.*((delta - 1)^2).^a; DELTAb = DELTA.^b; Q = sum(n(1:7).*delta.^d(1:7).*tau.^t(1:7)) + ... sum(n(8:51).*delta.^d(8:51).*tau.^t(8:51).*exp(-delta.^c(8:51))) + ... sum(n(52:54).*delta.^d(52:54).*tau.^t(52:54).*exp(-alpha(52:54).*(delta - ... epsilon(52:54)).^2 - beta(52:54).*(tau - gamma(52:54)).^2)) + ... sum(n(55:56).*DELTAb(55:56).*delta.*psi(55:56)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function Q = phirx_d(delta, tau) % Compute the partial derivative of the residual part of Helmholtz energy with respect to density global c d t a b B C D A alpha beta gamma epsilon n psi = exp(-C*(delta - 1)^2 - D*(tau - 1)^2); psi_d = -2.*C*(delta - 1).*psi; theta = 1 - tau + A.*((delta - 1)^2).^(1/2./beta); DELTA = theta.^2 + B.*((delta - 1)^2).^a; DELTAb = DELTA.^b; DELTA_d = (delta - 1)*(2*A.*theta./beta.*((delta - 1)^2).^(1/2./beta - 1) + ... 2*B.*a.*((delta - 1)^2).^(a-1)); DELTAb_d = b.*DELTA.^(b - 1).*DELTA_d; Q = sum(n(1:7).*d(1:7).*delta.^(d(1:7) - 1).*tau.^t(1:7)) + ... sum(n(8:51).*exp(-delta.^c(8:51)).*delta.^(d(8:51) - 1).*tau.^t(8:51).*... (d(8:51)-c(8:51).*delta.^c(8:51))) + ... sum(n(52:54).*delta.^d(52:54).*tau.^t(52:54).*exp(-alpha(52:54).*(delta - ... epsilon(52:54)).^2 - beta(52:54).*(tau - gamma(52:54)).^2).*(d(52:54)/... delta - 2*alpha(52:54).*(delta - epsilon(52:54)))) + ... sum(n(55:56).*(DELTAb(55:56).*(psi(55:56) + delta*psi_d(55:56)) + ... delta*DELTAb_d(55:56).*psi(55:56))); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function Q = phirx_dd(delta, tau) % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density global c d t a b B C D A alpha beta gamma epsilon n psi = exp(-C*(delta - 1)^2 - D*(tau - 1)^2); psi_d = -2.*C*(delta - 1).*psi; psi_dd = 2*C.*psi.*(2*C*(delta - 1)^2 - 1); theta = 1 - tau + A.*((delta - 1)^2).^(1/2./beta); DELTA = theta.^2 + B.*((delta - 1)^2).^a; DELTAb = DELTA.^b; DELTA_d = (delta - 1)*(2*A.*theta./beta.*((delta - 1)^2).^(1/2./beta - 1) + ... 2*B.*a.*((delta - 1)^2).^(a-1)); DELTAb_d = b.*DELTA.^(b - 1).*DELTA_d; DELTA_dd = 2*A.*theta./beta.*((delta - 1)^2).^(1/2./beta - 1) + ... 2*B.*a.*((delta - 1)^2).^(a-1) + (delta - 1)^2*(4*B.*a.*(a - 1).*... ((delta - 1)^2).^(a - 2) + 2*A.^2.*beta.^-2.*(((delta - 1)^2).^... (1/2./beta - 1)).^2) + 4*A.*theta./beta.*(2./beta - 1).*... ((delta - 1)^2).^(1/2./beta); DELTAb_dd = b.*(DELTA.^(b - 1).*DELTA_dd + (b - 1).*DELTA.^(b - 2).*DELTA_d.^2); Q = sum(n(1:7).*d(1:7).*(d(1:7) - 1).*delta.^(d(1:7) - 2).*tau.^t(1:7)) + ... sum(n(8:51).*exp(-delta.^c(8:51)).*delta.^(d(8:51) - 2).*tau.^t(8:51).*... ((d(8:51) - c(8:51).*delta.^c(8:51)).*(d(8:51) - 1 - ... c(8:51).*delta.^c(8:51)) - c(8:51).^2.*delta.^c(8:51))) + ... sum(n(52:54).*tau.^t(52:54).*exp(-alpha(52:54).*(delta - ... epsilon(52:54)).^2 - beta(52:54).*(tau - gamma(52:54)).^2).*... (-2*alpha(52:54).*delta.^d(52:54) + 4*alpha(52:54).^2.*delta.^d(52:54).*... (delta - epsilon(52:54)).^2 - 4*d(52:54).*alpha(52:54).*... delta.^(d(52:54) - 1).*(delta - epsilon(52:54)) + d(52:54).*(d(52:54) - ... 1).*delta.^(d(52:54) - 2))) + ... sum(n(55:56).*(DELTAb(55:56).*(2*psi_d(55:56) + ... delta*psi_dd(55:56)) + 2*DELTAb_d(55:56).*(psi(55:56) + ... delta*psi_d(55:56)) + delta*DELTAb_dd(55:56).*psi(55:56))); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function Q = phirx_t(delta, tau) % Compute the partial derivative of the residual part of Helmholtz energy with respect to temperature global c d t a b B C D A alpha beta gamma epsilon n psi = exp(-C*(delta - 1)^2 - D*(tau - 1)^2); psi_t = -2*D.*psi*(tau - 1); theta = 1 - tau + A.*((delta - 1)^2).^(1/2./beta); DELTA = theta.^2 + B.*((delta - 1)^2).^a; DELTAb = DELTA.^b; DELTAb_t = -2*theta.*b.*DELTA.^(b - 1); Q = sum(n(1:7).*t(1:7).*delta.^d(1:7).*tau.^(t(1:7) - 1)) + ... sum(n(8:51).*t(8:51).*delta.^d(8:51).*tau.^(t(8:51) - 1).*... exp(-delta.^c(8:51))) + ... sum(n(52:54).*delta.^d(52:54).*tau.^t(52:54).*exp(-alpha(52:54).*(delta - ... epsilon(52:54)).^2 - beta(52:54).*(tau - gamma(52:54)).^2).*... (t(52:54)/tau - 2*beta(52:54).*(tau - gamma(52:54)))) + ... sum(n(55:56)*delta.*(DELTAb_t(55:56).*psi(55:56) + DELTAb(55:56).*... psi_t(55:56))); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function Q = phirx_tt(delta, tau) % Compute the second partial derivative of the residual part of Helmholtz energy with respect to temperature global c d t a b B C D A alpha beta gamma epsilon n psi = exp(-C*(delta - 1)^2 - D*(tau - 1)^2); psi_t = -2*D.*psi*(tau - 1); psi_tt = 2*D.*psi.*(2*D*(tau - 1)^2 - 1); theta = 1 - tau + A.*((delta - 1)^2).^(1/2./beta); DELTA = theta.^2 + B.*((delta - 1)^2).^a; DELTAb = DELTA.^b; DELTAb_t = -2*theta.*b.*DELTA.^(b - 1); DELTAb_tt = 2*b.*DELTA.^(b - 1) + 4*theta.^2.*b.*(b - 1).*DELTA.^(b - 2); Q = sum(n(1:7).*t(1:7).*(t(1:7) - 1).*delta.^d(1:7).*tau.^(t(1:7) - 2)) + ... sum(n(8:51).*t(8:51).*(t(8:51) - 1).*delta.^d(8:51).*tau.^(t(8:51) - 2).*... exp(-delta.^c(8:51))) + ... sum(n(52:54).*delta.^d(52:54).*tau.^t(52:54).*exp(-alpha(52:54).*(delta - ... epsilon(52:54)).^2 - beta(52:54).*(tau - gamma(52:54)).^2).*... ((t(52:54)/tau - 2*beta(52:54).*(tau - gamma(52:54))).^2 - ... t(52:54)/tau^2 - 2*beta(52:54))) + ... sum(n(55:56)*delta.*(DELTAb_tt(55:56).*psi(55:56) + 2*DELTAb_t(55:56).*... psi_t(55:56) + DELTAb(55:56).*psi_tt(55:56))); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function Q = phirx_dt(delta, tau) % Compute the second partial derivative of the residual part of Helmholtz energy with respect to density and temperature global c d t a b B C D A alpha beta gamma epsilon n psi = exp(-C*(delta - 1)^2 - D*(tau - 1)^2); psi_d = -2.*C*(delta - 1).*psi; psi_t = -2*D.*psi*(tau - 1); psi_dt = 4*C.*D.*psi*(delta - 1)*(tau - 1); theta = 1 - tau + A.*((delta - 1)^2).^(1/2./beta); DELTA = theta.^2 + B.*((delta - 1)^2).^a; DELTA_d = (delta - 1)*(2*A.*theta./beta.*((delta - 1)^2).^(1/2./beta - 1) + ... 2*B.*a.*((delta - 1)^2).^(a-1)); DELTAb = DELTA.^b; DELTAb_d = b.*DELTA.^(b - 1).*DELTA_d; DELTAb_t = -2*theta.*b.*DELTA.^(b - 1); DELTAb_dt = -2*(delta - 1)*A.*b./beta.*DELTA.^(b - 1).*((delta - 1)^2).^... (1/2./beta - 1) - 2*theta.*b.*(b-1).*DELTA.^(b - 2).*DELTA_d; Q = sum(n(1:7).*d(1:7).*t(1:7).*delta.^(d(1:7) - 1).*tau.^(t(1:7) - 1)) + ... sum(n(8:51).*t(8:51).*delta.^(d(8:51) - 1).*tau.^(t(8:51) - 1).*(d(8:51) - ... c(8:51).*delta.^c(8:51)).*exp(-delta.^c(8:51))) + ... sum(n(52:54).*delta.^d(52:54).*tau.^t(52:54).*exp(-alpha(52:54).*(delta - ... epsilon(52:54)).^2 - beta(52:54).*(tau - gamma(52:54)).^2).*... (d(52:54)/delta - 2*alpha(52:54).*(delta - epsilon(52:54))).*... (t(52:54)/tau - 2*beta(52:54).*(tau - gamma(52:54)))) + ... sum(n(55:56).*(DELTAb(55:56).*(psi_t(55:56) + delta*psi_dt(55:56)) + ... delta*DELTAb_d(55:56).*psi_t(55:56) + DELTAb_t(55:56).*(psi(55:56) + ... delta*psi_d(55:56)) + delta*DELTAb_dt(55:56).*psi(55:56)));