Water1984.frink

View or download Water1984.frink in plain text format


/** This is an implementation of water properties based on the book
    Harr L., Kell, G.S., and Gallagher, J.S. NBS/NCR Steam Tables. United States: N. p., 1984.

    This distills at least 12 different studies of water properties, epecially
    vapor states, and contains detailed curve fits of these study results into
    unified equations in SI units.

    This follows the equations and tables in Appendix A.  Equation and constant
    numbers are the numbers defined in that appendix.

    UPDATE:  I have suspended development of this library in favor of the
    implementation of the International Association for the Properties of Water
    and Steam (IAPWS) formulation IAPWS-95.  That implementation can be found
    in water.frink.

    There are several reasons for this decision:

    * The Steam Tables book was *mostly* dimensionally-consistent and specified
      its units of measure, up until the points that it wasn't.  And these were
      crucial places, where it was difficult or impossible to proceed by
      guessing the units of measure.

      For example, all of the equations are described beginning in Appendix A.
      Only one page into this appendix, in equation A.5, they take the
      exponential function *of a unit with dimensions of density*.  You can't
      do that.  You can only take exponential of a dimensionless number.
      Attempts to guess what density to divide beta by in this equation were
      not helped by the book.  It could be g/cm^3, or the critical density of
      water, or something entirely different.  None of them made sense,
      especially due to the lack of worked examples.

    * Terminology is unspecific and vague.  The whole point of this
      formulation is to calculate an equation they just call the "Helmholtz
      Function".  That's the title of Appendix A where all the calculations
      are explained.  You have to research outside this book to understand
      that this actually means the "Helmholtz free energy" of a thermodynamic
      system.  Then you realize that the quantity computed in this book isn't
      actully the Helmholtz free energy, but a *dimensionless* version of the
      Helmholtz free energy, where the temperature is converted to a
      dimensionless number by "Tratio = T0/T" (where T0 is the critical
      temperature of water, (647.073 K)) and the density of something is maybe
      converted to the density of water (it's not specified if it's liquid or
      vapor phase in any place) which is maybe converted to a dimensionless
      form but the equations in the appendix don't describe how (see the point
      above.)  Then you discover that the result isn't the "Helmholtz free
      energy" nor a dimensionless version of the Helmholtz free energy, but it
      should be called something like the "specific dimensionless Helmholtz
      free energy" which maybe must be multiplied by the molar mass of water
      or something.  But you can't actually multiply it by the
      currently-best-known value of the molar mass of water because the
      formulation is just a bunch of curve-fits which assumed some old value
      of the molar mass of water, so you have to multiply by their old 6-digit
      estimate of the molar mass of water and hope that the results are still
      valid in today's SI.

    * The entire formulation of this book is a function that calculates what
      the authors call the "Helmholtz equation" but, as above, is probably
      better called something like the "specific dimensionless Helmholtz free
      energy" which is a single function that takes as its input [density,
      temperature] and only returns the "specific dimensionless Helmholtz free
      energy" (whatever that means to you) That's *all* this book shows you
      how to calculate.  If you know the temperature and pressure, you're out
      of luck.  You should have measured the density!  But is it the density
      of the liquid or vapor phase?!  LOL the book won't tell you.

    * Do you know the density of your pot of water to feed to this model?
      Get out your densitometer.  Do you not have a thing that measures
      density?  Why not?

    * In order to compute even the most basic properties of a system of water,
      like its pressure, you first need to know its density and temperature.
      LOL again, the book won't tell you if it's the density of liquid water or
      water vapor that it wants to know.

    * The book does note that if you want to calculate the pressure of a
      system, that you can use the equation "P = rho^2 dA/drho" where dA/drho
      is the partial derivative of A, (which is what the authors call the
      "Helmholtz equation" but see above to see what that means) with respect
      to rho, the density (but you probably don't know the density, and it's
      not clear if this is the vapor density or the liquid density.)  But this
      book doesn't give you the derivative functions of A, or, more precisely
      A[density, temperature] which is *the only thing this model calculates.*
      The authors don't specify the derivative of the equation, which is a sum
      of dozens of terms.  Are you supposed to derive the derivative yourself
      or obtain it numerically by taking small steps, or using the secant
      method, or what?  The IAPWS formulation is better because they
      explicitly enumerate the derivatives (and second derivatives) of the
      equation with respect to several variables.

    * The FORTRAN code is undocumented and inconsistent and mysterious.  There
      is an appendix supposedly containing FORTRAN code that generated the
      tables that make up the volume of the book.  But almost none of these
      calculations are addressed in the body of the book.  These will be
      broken down further:

    * The FORTRAN code is undocumented:  There are almost no comments anywhere
      in the FORTRAN code.

    * The FORTRAN code is inconsistent:  Where there are comments in the code,
      they refer to functions that don't actually exist.  For example, the
      function TSAT mentions that it uses the results of the function PSAT.
      There is no function PSAT.  Maybe they meant PS?

    * The FORTRAN code is mysterious:  Most of the functions in the FORTRAN
      code do things that aren't described in the appendix at all.  Even basic
      things like finding the saturation point of water at a specific
      temperature are deep dark mysteries hidden in the FORTRAN code with no
      equivalents in the text.  But the whole first chapter of the book is
      a giant table of saturation points of water.  But not a SINGLE WORD in
      the book mentions how they were calculated.

    * The code does not do any of the most basic things.  For example, the
      intent of implementing this library was to find the boiling point of
      water.  That's it.  I want a good estimate for the boiling point of
      water.  A definition of the "boiling point" of water is when the
      saturation pressure of the vapor phase of a gas equals the ambient
      atmospheric pressure.  The FORTRAN code to try and do this is utterly
      undocumented with lists of magical constants, no units of measure,
      nothing.

    * Parentheses are wrong, and give incorrect values if interpreted according
      to the rules of the SI, e.g.
      https://physics.nist.gov/cuu/Units/checklist.html
      For example, the discussion of equation A.6 says "T_R = T / 100 K".
      Under SI rules, or under the basic precedence rules of algebra that say
      that multiplication and division have the same precedence and are
      performed left to right, this is correctly interpreted as "(T/100) * K".
      But that's not what they mean, and you have to rewrite as T / (100 K) to
      get dimensionally-consistent results.

    * There are *no worked examples* showing partial results.  If you have an
      error somewhere in your calculations, all you know is that it doesn't
      match their undocumented results at the very end of a long chain of
      undocumented calculations.

    * Constants are old, and were probably not even the best-known constants
      when the book was written.  For example, the molar mass of water, the
      gas constant, the specific gas constant, the critical point of water,
      were probably not the best-known values at the time and are only
      specified to 6 decimal places.  If you try to use these calculations
      with modern constants, the results of the model are probably less
      correct.

    * There are no "inverse equations" for finding things like the boiling
      point of water.  (i.e. saturation vapor pressure = ambient pressure)
      Other formulations like IAPWS-95 and IAPWS-97 at least provide "good
      guess" equations for finding these quantities.  See water.frink for
      more.
*/


class Water1984
{
   /** Constants as defined in Steam Tables.  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.
   */

   
   /** Molar mass of water  [1] */
   class var M = 18.0152 g/mol

   /** Gas constant R (with overbar) [2] */ 
   class var Rbar = 8.31441 J / (mol K)

   /** Specific gas constant R = Rbar / M */
   class var R = Rbar / M

   /** A reference temperature used in several calculations. */
   class var  T0 = 647.073 K
   
   /** Tables used for the residual Helmholtz equations.  These are table
       A.2 in Steam Tables, so that's what we call the variables.

       The columns are k(i), l(i), G(i)
   */

   class var A2 = [[undef, undef, undef],         // no i=0 term
       [1, 1, -530.62968529023 J/g],  // i=1
       [1, 2, 2274.4901424408 J/g],   // i=2
       [1, 4, 2274.4901424408 J/g],   // i=3
       [1, 6, -69.830527374994 J/g],  // i=4
       [2, 1, 17863.832875422 J/g],   // i=5
       [2, 2, -39514.731563338 J/g],  // i=6
       [2, 4, 33803.884280753 J/g],   // i=7
       [2, 6, -13855.050202703 J/g],  // i=8
       [3, 1, -256374.36613260 J/g],  // i=9
       [3, 2, 482125.75981415 J/g],   // i=10
       [3, 4, -341830.16969660 J/g],  // i=11
       [3, 6, 122231.56417448 J/g],   // i=12
       [4, 1, 1179743.3655832 J/g],   // i=13
       [4, 2, -2173481.0110373 J/g],  // i=14
       [4, 4, 1082995.2168620  J/g],  // i=15
       [4, 6, -254419.98064049 J/g],  // i=16
       [5, 1, -3137777.4947767 J/g],  // i=17
       [5, 2, 5291191.0757704 J/g],   // i=18
       [5, 4, -1380257.7177877 J/g],  // i=19
       [5, 6, -251099.14369001 J/g],  // i=20
       [6, 1, 4656182.6115608 J/g],   // i=21
       [6, 2, -7275277.3275387 J/g],  // i=22
       [6, 4, 417742.46148294 J/g],   // i=23
       [6, 6, 1401635.8244614 J/g],   // i=24
       [7, 1, -3155523.1392127 J/g],  // i=25
       [7, 2, 4792966.6384584 J/g],   // i=26
       [7, 4, 409126.64781209 J/g],   // i=27
       [7, 6, -1362636.9388386 J/g],  // i=28
       [9, 1, 696252.20862664 J/g],   // i=29
       [9, 2, -1083490.0096447 J/g],  // i=30
       [9, 4, -227228.27401688 J/g],  // i=31
       [9, 6, 383654.86000660 J/g],   // i=32
       [3, 0, 6883.3257944332 J/g],   // i=33
       [3, 3, 21757.245522644 J/g],   // i=34
       [1, 3, -2662.7944829770 J/g],  // i=35
       [5, 3, -70730.418082074 J/g]]  // i=36

   /** This is terms 37-40 of Table A.2 which has more columns.

       The columns are
       [k(i), l(i), rho(i), T(i), alpha(i), beta(i), g(i)]
   */

   class var A2prime = [[2, 0, 0.319 g/cm^3,  640. K,   34, 20000, -0.225 J/g],   // i=37
    [2, 2, 0.319 g/cm^3,  640. K,   40, 20000,  -1.68 J/g],   // i=38
    [2, 0, 0.319 g/cm^3, 641.6 K,   30, 40000,  0.055 J/g],   // i=39
    [4, 0, 1.55  g/cm^3,  270. K, 1050,    25,  -93.0 J/g]]   // i=40


   /** This is table A.3, the coefficients for the ideal gas function. */
   class var C =  [undef, // No coefficient 0
                   1.97302710180e1, // C1
                   2.09662681977e1, // C2
                   -0.483429455355, // C3
                   6.05743189245,   // C4
                   2.256023885e1,   // C5
                   -9.875324420,    // C6
                   -4.3135538513,   // C7
                   4.581557810e-1,  // C8
                   -4.7754901883e-2,// C9
                   4.1238460633e-3, // C10
                   -2.7929052852e-4,// C11
                   1.4481695261e-5, // C12
                   -5.6473658748e-7,// C13
                   1.620044600e-8,  // C14
                  -3.3038227960e-10,// C15
                  4.51916067368e-12,// C16
                 -3.70734122708e-14,// C17
                  1.37546068238e-16]// C18
   
   /** This is the Helmholtz function (also known as the Helmholtz free energy)
       for fluid water, which adds contributions
       from the base equation, the residual equation, and the ideal gas
       equation.  This is equation A.1.
   */

   class A[rho is mass_density, T is temperature, debug = false] :=
   {
      ab = Abase[rho, T]
      ar = Aresidual[rho, T]
      ai = Aidealgas[T]

      if debug
         println["Abase=$ab\nAresidual=$ar\nAideal=$ai\n"]
      
      return ab + ar + ai
   }

   /** The base Helmholtz free energy equation.
       "The base function was derived from the Ursell-Mayer [4] virial theory
       for the equation of state for a fluid.   However, for real fluids, that
       theory is useful only for the dilute gas.  To obtain a theoretical
       equation  that is accurate at high values of density, the virial series
       had been transformed [5,6] into a perturbation expansion in density,
       the reference term for which is the equation of state for a fluid
       composed of molecules that are (nearly) hard ellipsoids.  The equation
       of state so obtained is a function of only two molecular parameters; it
       is easily integrated to yield the Helmholtz function."
   */

   class Abase[rho is mass_density, T is temperature] :=
   {
      T0ratio = T0/T

      // Calculate molecular parameters which are functions of the temperature
      // eq. A.3
      b = -0.3540782 cm^3 g^-1 ln[T/T0] +      // b1 term, from Table A.1
           0.7478629 cm^3 g^-1          +      // b0 term
           0.007159876 cm^3 g^-1 T0ratio^3  +  // b3 term
          -0.003528426 cm^3 g^-1 T0ratio^5     // b5 term

      // eq. A.4
      Bbar =  1.1278334 cm^3 g^-1           +  // B0 term, from Table A.1
             -0.5944001 cm^3 g^-1 T0ratio   +  // B1 term
             -5.010996  cm^3 g^-1 T0ratio^2 +  // B2 term
              0.63684256 cm^3 g^-1 T0ratio^4   // B4 term

      // Geometric constants
      alpha = 11
      beta  = 133/3
      gamma = 7/2
      
      y = b rho / 4

      // Equation A.2
      Abase = R T (-ln[1-y] - (beta-1)/(1-y) + (alpha+beta+1)/(2(1-y)^2) +
                   4y (Bbar/b - gamma) - (alpha-beta+3)/2 + ln[rho R T / atm])
      
      return Abase
   }


   /** The residual Helmholtz free energy function.
       "The residual function consists of 40 terms.  The first 36 of these are
       used in a global, least-squares fit to experimental data.  Their
       functional form, illustrated in Fig. A.1, is that of a standard growth
       curve: thus the contributions to the Helmholtz function asymptotically
       approach zero at the zero of density and are independent of density at
       high values of density; but they contribute importantly over the range
       of density values for which data exist.  Following this, three terms
       were added that contribute only in the immediate neighborhood of the
       critical point, (t_k - 5 degC) <= t <= (t_k + 5 degC),
       0.20 g/cm^3 <= rho <= 0.44 g/cm^3, and a single term that contributes
       only in the region of high values of pressure and low values of
       temperature, t < C[75] , P > 3 kbar (Fig A.2).  Except in these very
       limited regions, the residual function is given by the first 36 terms.
       A discussion of the Helmholtz function obtained with the residual
       function so restricted (terms 1-36) is given in [7]."
   */

   class Aresidual[rho is mass_density, T is temperature] :=
   {
      sum = 0 m^2/s^2

      for i = 1 to 36
      {
         [ki, li, gi] = A2@i  // Get the specified row from table A2

         // Eq. A.5, part 1
         // In this equation, there is a term e^-rho but this doesn't make
         // dimensional sense because rho has dimensions of density.  So
         // what factor do we divide by?  Water at the triple point?
         // g/cm^3?  Something else?
         sum = sum + gi/ki (T0 / T)^li (1-e^( -rho/(g/cm^3) ))^ki
      }

      // Eq. A.5, part 2
      for [ki, li, rhoi, Ti, alphai, betai, gi] = A2prime
      {
         deltai = (rho - rhoi) / rhoi
         taui = (T - Ti) / Ti
         sum = sum + gi deltai^li exp[-alphai deltai^ki - betai taui^2]
      }

      return sum
   }


   /** The ideal gas Helmholtz free energy function.
       "Values for thermodynamic properties for water in the ideal gas state
       have been reported recently by Wooley [8] as part of a detailed
       analysis of the rotation-vibration structure of the water molecule.
       Reference [8] very closely approximates a term by term summation of
       contributions from each of the vibration and rotation states.  That
       result is (this function, A.6) where T_R = T / (100 K).  The
       coefficients Ci are listed in Table A.3.  Equation (A.6) can be used
       for values of temperature below 3000 K.  It is accurate to within the
       tolerance of the gas constant over the range 50 K <= T <= 2000 K."

       This is equation A.6
   */

   class Aidealgas[T is temperature] :=
   {
      TR = T / (100 K)
      sum = 1 + (C@1 / TR + C@2) ln[TR]

      for i = 3 to 18
         sum = sum + C@i TR^(i-6)

      return -R T sum
   }


   class dAdrho[rho is mass_density, T is temperature, debug=false] :=
   {
      rho2 = rho * 1.000001
      Aa = A[rho, T, debug]
      Ab = A[rho2, T, debug]

      return (Aa-Ab) / (rho2-rho)
   }
}


View or download Water1984.frink in plain text format


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 18659 days, 15 hours, 12 minutes ago.