# water.frink

``` /** This class contains equations about the properties of water.     It is primarily based on The International Association for the Properties     of Water and Steam (IAPWS, http://iapws.org/) document IAPWS R6-95(2018),     a PDF of which is usually available from:     http://www.iapws.org/relguide/IAPWS-95.html     Specifically, the file     http://www.iapws.org/relguide/IAPWS95-2018.pdf     All references to equation and table numbers are the equation numbers in     this document.     The development of this document is discussed in     https://aip.scitation.org/doi/10.1063/1.1461829     There is an online calculator that can be used to validate these equations     at:     http://twt.mpei.ac.ru/mcs/worksheets/iapws/IAPWS95.xmcd     IAPWS also produces a simpler formulation "for industry", called IAPWS-IF97     which is described at:     http://www.iapws.org/relguide/IF97-Rev.html     That page contains a PDF of IAPWS-IF97 and some "backward equations" that     make it easier to calculate without iterating inverse solutions.  These     may be used to find good "first guesses" for inverse solutions.     Simpler equations for vapor pressure and densities at the saturation point     are found here:     http://www.iapws.org/relguide/Supp-sat.html     especially in the PDF file:     http://www.iapws.org/relguide/Supp-sat.pdf     These can be used to obtain good initial guesses for finding     higher-accuracy results in the full IAPWS95 model.     A scientific paper describing the derivation of the saturation properties     can be found at:     https://aip.scitation.org/doi/10.1063/1.555926     All of the methods in this class are class-level methods so you don't     need an instance of the Water class to call them.  Just call them like:     Water.boilingPoint[1 atmosphere] */ class Water {    /** Constants as defined in IAPWS RS-95(2018).  These may not exactly match        the current best-known values of these constants in the SI, but they are        used as defined there for purposes of matching its output exactly.    */    /** The temperature of the critical point of water. (eq.1) */    class var Tc = 647.096 K    /** The critical density of water. (eq.2) */    class var rhoc = 322 kg m^-3    /** The specific gas constant.  Due to the use of the specific gas constant,        Eq. (4) corresponds to a mass-based formulation.  In order to convert        values of specific properties to molar properties, an appropriate value        for the molar mass should be used  (eq. 3) */    class var R = 0.46151805 kJ kg^-1 K^-1    /** The critical pressure of water from IAPWS SR1-86(1992) */    class var pc = 22.064 MPa    /** The triple point of water. */    class var TriplePoint = 273.16 K        /** IAPWS doesn't directly define the molar mass of water directly but        indicates that the best-known values of several of these constants have        changed since the model was published.  It sort of waves its hand at        George S. Kell, Effects of isotopic composition, temperature, pressure,        and dissolved gases on the density of liquid water Journal of Physical        and Chemical Reference Data 6, 1109 (1977);        https://doi.org/10.1063/1.555561        to hint that it contains a suitable definition of the molar mass.  Its        masses and isotope ratios are taken from "Atomic Weights of the        Elements 1973", Pure Appl. Chem 37, 589-603 (1974) Of course, none of        this matches the current definition of Avogadro's constant, the gas        constant, etc.  Anyway, this is the version of the molar mass of water        defined in that paper, and the one that should probably be used with        this model and its definition of R.  See the notes for R above.  */    class var molarMassH2O = 18.015242 g/mol                    /** Calculates the pressure at the specified density and temperature.                p = rho^2 (df/drho)_T                See Table 3 for this equation.    */    class pressure[rho is mass_density, T is temperature] :=    {       delta = rho / rhoc       return rho R T (1 + delta dphirdeltaDimensioned[rho, T])    }    /** Calculates the compressibility factor Z = p / (rho R T) */    class compressibilityFactor[rho is mass_density, T is temperature] :=    {       return pressure[rho, T] / (rho R T)    }        /** Calculates the specific internal energy of the system at the specified        density and temperature.  Multiply this by mass to get the energy.        u = f - T (df/dT)_rho        returns energy / mass    */    class specificInternalEnergy[rho is mass_density, T is temperature] :=    {       delta = rho / rhoc       tau = Tc / T       return tau R T dphidtau[delta, tau]    }    /** Calculates the specific entropy        s = -(df/dT)_rho        Multiply by mass to get the total entropy.            This has units of specific heat capacity (e.g. cal/g/degC)    */    class specificEntropy[rho is mass_density, T is temperature] :=    {       delta = rho / rhoc       tau = Tc / T       return R (tau dphidtau[delta, tau] - phi[delta, tau])    }            /** Calculates the specific enthalpy        h = f - T(df/dT)_rho + rho(df/drho)_T        Multiply by mass to get the total enthalpy.            This has units of energy/mass    */    class specificEnthalpy[rho is mass_density, T is temperature] :=    {       delta = rho / rhoc       tau = Tc / T       return R T (1 + tau dphidtau[delta, tau] + delta dphirdelta[delta, tau])    }    /** Calculates the specific isochoric heat capacity        c_v = (du/dT)_rho        This has units of specific heat capacity (e.g. cal/g/degC)    */    class specificIsochoricHeatCapacity[rho is mass_density, T is temperature] :=    {       delta = rho / rhoc       tau = Tc / T       return R (-tau^2 d2phitau[delta, tau])    }            /** Calculates the specific isobaric heat capacity        c_p = (dh/dT)_p            This has units of specific heat capacity (e.g. cal/g/degC)    */    class specificIsobaricHeatCapacity[rho is mass_density, T is temperature] :=    {       delta = rho / rhoc       tau = Tc / T       return R (-tau^2 d2phitau[delta, tau] + (1 + delta dphirdelta[delta, tau] - delta tau d2phirdeltatau[delta, tau])^2 / (1 + 2 delta dphirdelta[delta, tau] + delta^2 d2phirdelta[delta, tau]))    }        /** Calculates the speed of sound.        w = (dp/drho)^(1/2)            This has units of velocity    */    class speedOfSound[rho is mass_density, T is temperature] :=    {       delta = rho / rhoc       tau = Tc / T       return sqrt[R T (1 + 2 delta dphirdelta[delta, tau] + delta^2 d2phirdelta[delta,tau] - (1 + delta dphirdelta[delta, tau] - delta tau d2phirdeltatau[delta, tau])^2/(tau^2 (d2phitau[delta, tau])))]    }        /** The Joule-Thompson coefficient        mu = (dT/drho)_h        Results have units of temperature / pressure    */    class JouleThompsonCoefficient[rho is mass_density, T is temperature] :=    {       delta = rho / rhoc       tau = Tc / T       return (-(delta dphirdelta[delta, tau] + delta^2 d2phirdelta[delta,tau] + delta tau d2phirdeltatau[delta,tau]) / ((1 + delta dphirdelta[delta, tau] - delta tau d2phirdeltatau[delta,tau])^2 - tau^2 d2phitau[delta,tau] (1 + 2 delta dphirdelta[delta, tau] + delta^2 d2phirdelta[delta, tau]))) / (R rho)    }    /** The isentropic temperature-pressure coefficent        beta_s = (dT/dp)_s        Results have units of temperature / pressure    */    class isentropicTemperaturePressureCoefficient[rho is mass_density, T is temperature] :=    {       delta = rho / rhoc       tau = Tc / T       return  (1 + delta dphirdelta[delta, tau] - delta tau d2phirdeltatau[delta, tau]) / ((1 + delta dphirdelta[delta, tau] - delta tau d2phirdeltatau[delta,tau])^2 - tau^2 d2phitau[delta,tau] (1 + 2 delta dphirdelta[delta, tau] + delta^2 d2phirdelta[delta, tau])) / (R rho)    }        /** The isothermal throttling coefficient        delta_T = (dh/dp)_T        This has units of inverse density (e.g. m^3 / kg)    */    class isothermalThrottlingCoefficient[rho is mass_density, T is temperature] :=    {       delta = rho / rhoc       tau = Tc / T       return (1-(1 + delta dphirdelta[delta, tau] - delta tau d2phirdeltatau[delta, tau])/ (1 + 2 delta dphirdelta[delta, tau] + delta^2 d2phirdelta[delta,tau])) / rho    }    /** The second virial coefficient B(T).            "Even at low density, real gases don't quite obey the ideal gas law.         A systematic way to account for deviations from ideal behavior is the         virial expansion,         P V = n R T (1 + B(T)/(V/n) + C(T)/(V/n)^2 + ...)         where the functions B(T), C(T), and so on are called the "virial         coefficents".  When the density of the gas is fairly low, so that the         volume per mole is large, each term in the series is much smaller than         the one before.  In many situations it's sufficient to omit the third         term and concentrate on the second, whose coefficent B(T) is called the         second virial coefficient (the first coefficient being 1.)"         --Daniel V. Schroeder, An Introduction to Thermal Physics, p.9         B(tau) rhoc = lim[phirdelta(delta, tau), delta->0]         Units are m^3 kg^-1         This actually returns the *specific* version of B.         To get the dimensionless units that will work in the equation above,         you need to multiply by Water.molarMassH2O    */    class B[T is temperature] :=    {       tau = Tc / T       sum = 0              for i = [1,2,3,8,9,10,23]       {          [ci, di, ti, ni] = Table2@i          sum = sum + ni tau^ti       }       for i = 55 to 56       {          [ai, bi, Bi, ni, Ci, Di, Ai, betai] = Table2@i          thetai = Ai + 1 - tau       // See footnote b, p.12          sum = sum + ni (thetai^2 + Bi)^bi exp[-Ci - Di(tau-1)^2]       }              return sum/rhoc    }        /** The third virial coefficient C(T).  See the discussion for the second        virial coefficient B[T] above.        Units are m^6 kg^-2        This actually returns the *specific* version of C.        To get the dimensionless units for the equation documented in B[T],        you need to multiply by Water.molarMassH2O    */    class C[T is temperature] :=    {       tau = Tc / T       sum = 0              for i = [4,5,11,12,24,25,26]       {          [ci, di, ti, ni] = Table2@i          sum = sum + 2 ni tau^ti       }       for i = 8 to 10       {          [ci, di, ti, ni] = Table2@i          sum = sum - 2 ni tau^ti       }              for i = 55 to 56       {          [ai, bi, Bi, ni, Ci, Di, Ai, betai] = Table2@i          thetai = Ai + 1 - tau       // See footnote b, p.12          sum = sum + 4 ni (Ci (thetai^2 + Bi) - bi (Ai thetai / betai + Bi ai))* (thetai^2 + Bi)^(bi-1) exp[-Ci - Di(tau-1)^2]       }              return sum/rhoc^2    }        /** Calculates the specific internal energy of the system at the specified        atmospheric pressure equals the saturated vapor pressure.   You can use        StandardAtmosphereTest.frink which calls StandardAtmosphere.frink        to find the approximate atmospheric pressure for any altitude.    */    class boilingPoint[ps is pressure] :=    {       return saturatedTemperature[ps]    }        /** Returns the approximate saturated vapor pressure of water from        IAPWS SR1-86(1992).  This does not use the full IAPWS95 model, but it        can be used to get a very good first guess.   This is Section 3 */    class saturatedVaporPressure[T is temperature] :=    {       theta = T/Tc       tau = 1 - theta       p = pc exp[Tc/T (-7.85951783 tau +                         1.84408259 tau^(3/2) +                        -11.7866497 tau^3 +                         22.6807411 tau^(7/2) +                        -15.9618719 tau^4 +                         1.80122502 tau^(15/2))]       return p    }        /** Returns the approximate saturated vapor pressure of water from        IAPWS IF-97.  This does not use the full IAPWS95 model, but it        can be used to get a very good first guess.   This is Section 8 and        8.1 of IAPWS IF-97 */    class saturatedVaporPressure97[Ts is temperature] :=    {       Tstar = 1 K       Theta = Ts / Tstar + Table34@9 / ((Ts/Tstar) - Table34@10)  // eq. 29b       A = Theta^2 + Table34@1 Theta + Table34@2       B = Table34@3 Theta^2 + Table34@4 Theta + Table34@5       C = Table34@6 Theta^2 + Table34@7 Theta + Table34@8       p = 1 MPa ((2 C) / (-B + (B^2 - 4 A C)^(1/2)))^4   // eq. 30       return p    }        /** Returns the approximate saturated vapor temperature of water from        IAPWS97.  This does not use the full IAPWS95 model, but it        can be used to get a very good first guess.   This is Section 8.2 of        IAPWS97. */    class saturatedTemperature[ps is pressure] :=    {       // Equation 31       pstar = 1 MPa       beta = (ps / pstar)^(1/4)  // Eq. 29a              E = beta^2 + Table34@3 beta + Table34@6       F = Table34@1 beta^2 + Table34@4 beta + Table34@7       G = Table34@2 beta^2 + Table34@5 beta + Table34@8       D = (2 G)/(-F - (F^2 - 4 E G)^(1/2))              Ts = 1 K ((Table34@10 + D - ((Table34@10 + D)^2 - 4(Table34@9 + Table34@10 D))^(1/2)) / 2)       return Ts    }        /** Returns the approximate saturated liquid density of water from        IAPWS SR1-86(1992).  This does not use the full IAPWS95 model, but it        can be used to get a very good first guess.   This is Section 4.1 */    class saturatedLiquidDensity[T is temperature] :=    {       theta = T/Tc       tau = 1-theta       rho = rhoc (1 + 1.99274064 tau^(1/3)    + // b1                       1.09965342 tau^(2/3)    + // b2                     -0.510839303 tau^(5/3)    + // b3                      -1.75493479 tau^(16/3)   + // b4                       -45.5170352 tau^(43/3)   + // b5                    -6.74694450e5 tau^(110/3))   // b6       return rho    }        /** Returns the approximate saturated vapor density of water from        IAPWS SR1-86(1992).  This does not use the full IAPWS95 model, but it        can be used to get a very good first guess.   This is Section 4.2 */    class saturatedVaporDensity[T is temperature] :=    {       theta = T/Tc       tau = 1-theta       rho = rhoc exp[1 + -2.03150240 tau^(2/6)    + // c1                          -2.68302940 tau^(4/6)    + // c2                          -5.38626492 tau^(8/6)    + // c3                          -17.2991605 tau^(18/6)   + // c4                           -44.7586581 tau^(37/6)   + // c5                          -63.9201063 tau^(71/6)]    // c6       return rho    }    /** The specific Helmholtz free energy f (also sometimes called A).        The IAPWS formulation is implemented as a fundamental equation for the        specific Helmholtz free energy f.                This equation is expressed in dimensionless form, phi = f/(R T),        or, conversely f = phi R T        (see eq. 4 and following)    */    class f[rho is mass_density, T is temperature, debug=false] :=    {       return R T phiDimensioned[rho, T, debug]    }        /** The formulation of the phi equation is implemented as a dimensionless        version phi[delta, tau] where        delta = rho / rhoc        tau = Tc / T            This turns a dimensionally-correct call with density and temperature        into the internal dimensionless form.  (see eq.4 and following)    */    class phiDimensioned[rho is mass_density, T is temperature, debug=false] :=    {       return phi[rho / rhoc, Tc / T, debug]    }        /** The dimensionless phi function phi[delta, tau] is separated into two        parts, an ideal-gas part phi0[delta, tau] and a residual part        phir[delta, tau]   (eq. 4)    */    class phi[delta is dimensionless, tau is dimensionless, debug=false] :=    {       phi0 = phi0[delta, tau]       phir = phir[delta, tau]       if debug          println["phi0=\$phi0\nphir=\$phir\n"]              return phi0 + phir    }    /** Calculates the first derivative of phi with respect to tau.  This        includes both the ideal gas and residual parts. */    class dphidtau[delta is dimensionless, tau is dimensionless] :=    {       return dphi0dtau[delta, tau] + dphirtau[delta, tau]    }        /** Calculates the second derivative of phi with respect to tau.  This        includes both the ideal gas and residual parts. */    class d2phitau[delta is dimensionless, tau is dimensionless] :=    {       return d2phi0dtau[delta, tau] + d2phirtau[delta, tau]    }    /** Calculates the second derivative of phi with respect to delta.  This        includes both the ideal gas and residual parts. */    class d2phidelta[delta is dimensionless, tau is dimensionless] :=    {       return d2phi0ddelta[delta, tau] + d2phirdelta[delta, tau]    }    /** The ideal gas part phi0 of the dimensionless Helmholtz free energy.        This is obtained from an equation for the specific isobaric heat        capacity in the ideal-gas state developed by Cooper         This implements eq. 5    */    class phi0[delta is dimensionless, tau is dimensionless] :=    {       sum = ln[delta] + Table1@1  +       /* n1 0 */                         Table1@2 tau +    /* n2 0 */                         Table1@3 ln[tau]  /* n3 0 */       for i = 4 to 8       {          [ni, gammai] = Table1@i          sum = sum + ni ln[1 - exp[-gammai tau]]       }              return sum    }        /** The first partial derivative of the ideal gas part of the dimensionless        Helmholtz free energy with respect to delta. */    class dphi0ddelta[delta is dimensionless, tau is dimensionless] :=    {       return 1/delta    }        /** The second partial derivative of the ideal gas part of the dimensionless        Helmholtz free energy with respect to delta. */    class d2phi0ddelta[delta is dimensionless, tau is dimensionless] :=    {       return -1 / delta^2    }        /** The first partial derivative of the ideal gas part of the dimensionless        Helmholtz free energy with respect to tau. */    class dphi0dtau[delta is dimensionless, tau is dimensionless] :=    {       sum = Table1@2 + Table1@3 / tau       for i=4 to 8       {          [ni, gammai] = Table1@i          sum = sum + ni gammai ((1 - exp[-gammai tau])^-1 - 1)       }       return sum    }        /** The second partial derivative of the ideal gas part of the dimensionless        Helmholtz free energy with respect to tau. */    class d2phi0dtau[delta is dimensionless, tau is dimensionless] :=    {       sum = -Table1@3 / tau^2       for i=4 to 8       {          [ni, gammai] = Table1@i          sum = sum - ni gammai^2 exp[-gammai tau] (1 - exp[-gammai tau])^-2       }       return sum    }        /** The second partial derivative of the ideal gas part of the dimensionless        Helmholtz free energy with respect to delta then with respect to tau. */    class d2phi0ddeltadtau[delta is dimensionless, tau is dimensionless] :=    {       return 0    }        /** The residual part phir[delta, tau] of the dimensionless Helmholtz free        energy.        This implements eq. 6    */    class phir[delta is dimensionless, tau is dimensionless] :=    {       sum = 0       for i=1 to 7       {          [ci, di, ti, ni] = Table2@i          sum = sum + ni delta^di tau^ti       }       for i=8 to 51       {          [ci, di, ti, ni] = Table2@i          sum = sum + ni delta^di tau^ti e^(-delta^ci)       }       for i=52 to 54       {          [ci, di, ti, ni, alphai, betai, gammai, epsiloni] = Table2@i          sum = sum + ni delta^di tau^ti * exp[-alphai (delta - epsiloni)^2 - betai(tau - gammai)^2]       }       for i=55 to 56       {          [ai, bi, Bi, ni, Ci, Di, Ai, betai] = Table2@i          Psi = exp[-Ci (delta-1)^2 - Di (tau - 1)^2]          theta = (1 - tau) + Ai ((delta-1)^2)^(1/(2 betai))          Delta = theta^2 + Bi ((delta-1)^2)^ai          sum = sum + ni Delta^bi delta Psi       }       return sum    }        /** The partial derivative of phir[delta, tau] with respect to delta, with        correct dimensions. */    class dphirdeltaDimensioned[rho is mass_density, T is temperature] :=    {       return dphirdelta[rho / rhoc, Tc / T]    }        /** The partial derivative of phir[delta, tau] with respect to delta.        This implements eq. 6    */    class dphirdelta[delta is dimensionless, tau is dimensionless] :=    {       sum = 0       for i=1 to 7       {          [ci, di, ti, ni] = Table2@i          sum = sum + ni di delta^(di-1) tau^ti       }       for i=8 to 51       {          [ci, di, ti, ni] = Table2@i          sum = sum + ni exp[-delta^ci] (delta^(di-1) tau^ti(di - ci delta^ci))       }       for i=52 to 54       {          [ci, di, ti, ni, alphai, betai, gammai, epsiloni] = Table2@i          sum = sum + ni delta^di tau^ti exp[-alphai (delta-epsiloni)^2 - betai (tau-gammai)^2] (di/delta - 2 alphai (delta - epsiloni))       }              for i=55 to 56       {          [ai, bi, Bi, ni, Ci, Di, Ai, betai] = Table2@i          Psi = exp[-Ci (delta-1)^2 - Di (tau - 1)^2]          theta = (1 - tau) + Ai ((delta-1)^2)^(1/(2 betai))          Delta = theta^2 + Bi ((delta-1)^2)^ai          dPsiddelta = -2 Ci (delta-1) Psi          dDeltaddelta = (delta-1) (Ai theta 2 / betai ((delta-1)^2)^(1/(2 betai)-1) + 2 Bi ai((delta-1)^2)^(ai-1))          dDeltabiddelta = bi Delta^(bi-1) dDeltaddelta          sum = sum + ni (Delta^bi (Psi + delta dPsiddelta) + dDeltabiddelta delta Psi)       }       return sum    }        /** The second partial derivative of phir[delta, tau] with respect to delta,        with correct dimensions. */    class d2phirdeltaDimensioned[rho is mass_density, T is temperature] :=    {       return d2phirdelta[rho / rhoc, Tc / T]    }        /** The second partial derivative of phir[delta, tau] with respect to delta.    */    class d2phirdelta[delta is dimensionless, tau is dimensionless] :=    {       sum = 0       for i=1 to 7       {          [ci, di, ti, ni] = Table2@i          sum = sum + ni di (di-1) delta^(di-2) tau^ti       }       for i=8 to 51       {          [ci, di, ti, ni] = Table2@i          sum = sum + ni exp[-delta^ci] (delta^(di-2) tau^ti ((di - ci delta^ci)(di-1-ci delta^ci) - ci^2 delta^ci))       }       for i=52 to 54       {          [ci, di, ti, ni, alphai, betai, gammai, epsiloni] = Table2@i          sum = sum + ni tau^ti exp[-alphai (delta-epsiloni)^2 - betai(tau-gammai)^2] * (-2 alphai delta^di + 4 alphai^2 delta^di (delta-epsiloni)^2 - 4 di alphai delta^(di-1) (delta-epsiloni) + di (di-1) delta^(di-2))       }              for i=55 to 56       {          [ai, bi, Bi, ni, Ci, Di, Ai, betai] = Table2@i          Psi = exp[-Ci (delta-1)^2 - Di (tau - 1)^2]          theta = (1 - tau) + Ai ((delta-1)^2)^(1/(2 betai))          Delta = theta^2 + Bi ((delta-1)^2)^ai          dPsiddelta = -2 Ci (delta-1) Psi          dDeltaddelta = (delta-1) (Ai theta 2 / betai ((delta-1)^2)^(1/(2 betai)-1) + 2 Bi ai((delta-1)^2)^(ai-1))          d2Psiddelta = (2 Ci (delta-1)^2 -1) 2 Ci Psi          dDeltabiddelta = bi Delta^(bi-1) dDeltaddelta          d2Deltaddelta = 1 / (delta-1) dDeltaddelta + (delta-1)^2 (4 Bi ai (ai-1)((delta-1)^2)^(ai-2) + 2 Ai^2 (1/betai)^2 (((delta-1)^2)^(1/(2 betai) - 1))^2 + Ai theta 4/betai (1 / (2 betai) - 1)((delta-1)^2)^(1/(2 betai) - 2))          d2Deltabiddelta = bi (Delta^(bi-1) d2Deltaddelta + (bi - 1) Delta^(bi-2) dDeltaddelta^2)          sum = sum + ni (delta^bi (2 dPsiddelta + delta d2Psiddelta) + 2 dDeltabiddelta (Psi + delta dPsiddelta) + d2Deltabiddelta delta Psi)       }       return sum    }        /** The partial derivative of phir[delta, tau] with respect to tau, with        correct dimensions. */    class dphirtauDimensioned[rho is mass_density, T is temperature] :=    {       return dphirtau[rho / rhoc, Tc / T]    }        /** The partial derivative of phir[delta, tau] with respect to tau.    */    class dphirtau[delta is dimensionless, tau is dimensionless] :=    {       sum = 0       for i=1 to 7       {          [ci, di, ti, ni] = Table2@i          sum = sum + ni ti delta^di tau^(ti-1)       }       for i=8 to 51       {          [ci, di, ti, ni] = Table2@i          sum = sum + ni ti delta^di tau^(ti-1) exp[-delta^ci]       }       for i=52 to 54       {          [ci, di, ti, ni, alphai, betai, gammai, epsiloni] = Table2@i          sum = sum + ni delta^di tau^ti exp[-alphai (delta-epsiloni)^2 - betai (tau-gammai)^2] (ti/tau - 2 betai(tau - gammai))       }              for i=55 to 56       {          [ai, bi, Bi, ni, Ci, Di, Ai, betai] = Table2@i          Psi = exp[-Ci (delta-1)^2 - Di (tau - 1)^2]          theta = (1 - tau) + Ai ((delta-1)^2)^(1/(2 betai))          Delta = theta^2 + Bi ((delta-1)^2)^ai          dDeltabidtau = -2 theta bi Delta^(bi-1)          dPsidtau = -2 Di (tau-1) Psi          sum = sum + ni Delta (dDeltabidtau Psi + Delta^bi dPsidtau)       }       return sum    }        /** The partial second derivative of phir[delta, tau] with respect to tau,        with correct dimensions. */    class d2phirtauDimensioned[rho is mass_density, T is temperature] :=    {       return d2phirtau[rho / rhoc, Tc / T]    }        /** The partial second derivative of phir[delta, tau] with respect to tau.    */    class d2phirtau[delta is dimensionless, tau is dimensionless] :=    {       sum = 0       for i=1 to 7       {          [ci, di, ti, ni] = Table2@i          sum = sum + ni ti (ti-1) delta^di tau^(ti-2)       }       for i=8 to 51       {          [ci, di, ti, ni] = Table2@i          sum = sum + ni ti (ti-1) delta^di tau^(ti-2) exp[-delta^ci]       }       for i=52 to 54       {          [ci, di, ti, ni, alphai, betai, gammai, epsiloni] = Table2@i          sum = sum + ni delta^di tau^ti exp[-alphai (delta-epsiloni)^2 - betai (tau-gammai)^2] ((ti/tau - 2 betai(tau - gammai))^2 - ti/tau^2 - 2 betai)       }              for i=55 to 56       {          [ai, bi, Bi, ni, Ci, Di, Ai, betai] = Table2@i          Psi = exp[-Ci (delta-1)^2 - Di (tau - 1)^2]          theta = (1 - tau) + Ai ((delta-1)^2)^(1/(2 betai))          Delta = theta^2 + Bi ((delta-1)^2)^ai          dDeltabidtau = -2 theta bi Delta^(bi-1)          dPsidtau = -2 Di (tau-1) Psi          d2Psidtau = (2 Di (tau-1)^2 - 1) 2 Di Psi          d2Deltabidtau = 2 bi Delta^(bi-1) + 4 theta^2 bi (bi-1) Delta^(bi-2)          sum = sum + ni delta (d2Deltabidtau Psi + 2 dDeltabidtau dPsidtau + Delta^bi d2Psidtau)       }       return sum    }        /** The partial second derivative of phir[delta, tau] with respect to tau,        with correct dimensions. */    class d2phirdeltatauDimensioned[rho is mass_density, T is temperature] :=    {       return d2phirdeltatau[rho / rhoc, Tc / T]    }        /** The partial second derivative of phir[delta, tau] with respect to         delta then to tau.    */    class d2phirdeltatau[delta is dimensionless, tau is dimensionless] :=    {       sum = 0       for i=1 to 7       {          [ci, di, ti, ni] = Table2@i          sum = sum + ni di ti delta^(di-1) tau^(ti-1)       }       for i=8 to 51       {          [ci, di, ti, ni] = Table2@i          sum = sum + ni ti delta^(di-1) tau^(ti-1) (di - ci delta^ci) exp[-delta^ci]       }       for i=52 to 54       {          [ci, di, ti, ni, alphai, betai, gammai, epsiloni] = Table2@i          sum = sum + ni delta^di tau^ti exp[-alphai (delta-epsiloni)^2 - betai (tau-gammai)^2] (di/delta - 2 alphai (delta - epsiloni))(ti/tau - 2 betai(tau - gammai))       }              for i=55 to 56       {          [ai, bi, Bi, ni, Ci, Di, Ai, betai] = Table2@i          Psi = exp[-Ci (delta-1)^2 - Di (tau - 1)^2]          theta = (1 - tau) + Ai ((delta-1)^2)^(1/(2 betai))          Delta = theta^2 + Bi ((delta-1)^2)^ai          dDeltabidtau = -2 theta bi Delta^(bi-1)          dPsidtau = -2 Di (tau-1) Psi          d2Psidtau = (2 Di (tau-1)^2 - 1) 2 Di Psi          dPsiddelta = -2 Ci (delta-1) Psi          d2Deltabidtau = 2 bi Delta^(bi-1) + 4 theta^2 bi (bi-1) Delta^(bi-2)          dDeltaddelta = (delta-1) (Ai theta 2 / betai ((delta-1)^2)^(1/(2 betai)-1) + 2 Bi ai((delta-1)^2)^(ai-1))          dDeltabiddelta = bi Delta^(bi-1) dDeltaddelta          d2Psiddeltadtau = 4 Ci Di (delta-1)(tau-1) Psi          d2Deltabiddeltadtau = -Ai bi 2 / betai Delta^(bi-1) (delta-1)((delta-1)^2)^(1/(2 betai) - 1) - 2 theta bi (bi-1) Delta^(bi-2) dDeltaddelta          sum = sum + ni (Delta^(bi) (dPsidtau + delta d2Psiddeltadtau) + delta dDeltabiddelta dPsidtau + dDeltabidtau (Psi + delta dPsiddelta) + d2Deltabiddeltadtau delta Psi)       }       return sum    }       /** Numerical values of the coefficients and parameters of the ideal-gas        part of the dimensionless Helmholtz free energy, Eq. 5.  This is        table 1 in IAPWS-95.        The coluns are [ni0, gammai0]    */    class var Table1 = [undef,      // No element 0        -8.3204464837497,         6.6832105275932,         3.00632,        [0.012436, 1.28728967],        [0.97315,  3.53734222],        [1.27950,  7.74073708],        [0.96956,  9.24437796],        [0.24873, 27.5075105]]            /** Numerical values of the coefficients and parameters of the        residual part of the dimensionless Helmholtz free energy, Eq. (6) */    class var Table2 = [undef,    // No element 0     //[ci, di,  ti,           ni]            [x,  1, -1/2,   0.12533547935523e-1],        [x,  1, 7/8,    0.78957634722828e1],        [x,  1,  1,    -0.87803203303561e1],        [x,  2, 1/2,    0.31802509345418],        [x,  2, 3/4,   -0.26145533859358],        [x,  3, 3/8,   -0.78199751687981e-2],        [x,  4,  1,     0.88089493102134e-2],        [1,  1,  4,    -0.66856572307965],        [1,  1,  6,     0.20433810950965],        [1,  1, 12,    -0.66212605039687e-4],        [1,  2,  1,    -0.19232721156002],        [1,  2,  5,    -0.25709043003438],        [1,  3,  4,     0.16074868486251],        [1,  4,  2,    -0.40092828925807e-1],        [1,  4, 13,     0.39343422603254e-6],        [1,  5,  9,    -0.75941377088144e-5],        [1,  7,  3,     0.56250979351888e-3],        [1,  9,  4,    -0.15608652257135e-4],        [1, 10, 11,     0.11537996422951e-8],        [1, 11,  4,     0.36582165144204e-6],        [1, 13, 13,    -0.13251180074668e-11],        [1, 15,  1,    -0.62639586912454e-9],        [2,  1,  7,    -0.10793600908932],        [2,  2,  1,     0.17611491008752e-1],        [2,  2,  9,     0.22132295167546],        [2,  2, 10,    -0.40247669763528],        [2,  3, 10,     0.58083399985759],        [2,  4,  3,     0.49969146990806e-2],        [2,  4,  7,    -0.31358700712549e-1],        [2,  4, 10,    -0.74315929710341],        [2,  5, 10,     0.47807329915480],        [2,  6,  6,     0.20527940895948e-1],        [2,  6, 10,    -0.13636435110343],        [2,  7, 10,     0.14180634400617e-1],        [2,  9,  1,     0.83326504880713e-2],        [2,  9,  2,    -0.29052336009585e-1],        [2,  9,  3,     0.38615085574206e-1],        [2,  9,  4,    -0.20393486513704e-1],        [2,  9,  8,    -0.16554050063734e-2],        [2, 10,  6,     0.19955571979541e-2],        [2, 10,  9,     0.15870308324157e-3],        [2, 12,  8,    -0.16388568342530e-4],        [3,  3, 16,     0.43613615723811e-1],        [3,  4, 22,     0.34994005463765e-1],        [3,  4, 23,    -0.76788197844621e-1],        [3,  5, 23,     0.22446277332006e-1],        [4, 14, 10,    -0.62689710414685e-4],        [6,  3, 50,    -0.55711118565645e-9],        [6,  6, 44,    -0.19905718354408],        [6,  6, 46,     0.31777497330738],        [6,  6, 50,    -0.11841182425981],     // [ci, di, ti,         ni,          alphai, betai, gammai, epsiloni]        [x,  3,  0,  -0.31306260323435e2,   20,    150,   1.21,      1],        [x,  3,  1,   0.31546140237781e2,   20,    150,   1.21,      1],        [x,  3,  4,  -0.25213154341695e4,   20,    250,   1.25,      1],     // [ai,   bi,    Bi,         ni,         Ci, Di,  Ai,   betai]        [7/2, 85/100, 0.2, -0.14874640856724, 28, 700, 0.32, 3/10],        [7/2, 95/100, 0.2,  0.31806110878444, 32, 800, 0.32, 3/10]]            /** Numerical values for IAPWS97 table 34.  These are used to solve        saturation-pressure and saturation-temperature equations. */    class var Table34 = [undef,  // No n0         0.11670521452767e4, // n1        -0.72421316703206e6, // n2        -0.17073846940092e2, // n3         0.12020824702470e5, // n4        -0.32325550322333e7, // n5         0.14915108613530e2, // n6        -0.48232657361591e4, // n7         0.40511340542057e6, // n8        -0.23855557567849,   // n9         0.65017534844798e3] // n10 }```

This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 19032 days, 6 hours, 47 minutes ago.