/** 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)
}
}