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