THE EQUATION OF STATE FOR WATER IAPWS-95 PROGRAMMED IN MATHEMATICA 11.2

Mario-César Suárez-Arriaga, International Geothermal Association

Join the Conversation #WolframTechConf

Abstract

★	Together with air, water is the most important substance existing on Earth for living organisms, in science and in many industrial applications: computational fluid dynamics, steam turbines, high pressure devices and supercritical cycles.  

◆	The International Association for the Properties of Water and Steam (*IAPWS) published in 1995, a worldwide accepted mathematical formulation, which reproduces with maximum accuracy millions of experimental data obtained since 1824 up to present time. This formulation is presented here in the form of a computable Notebook written in the Wolfram programming Language of Mathematica version 11.2. 

◆ 	The water properties are calculated for single phase, liquid or vapor, as function of (ρ, T), and for two-phase water (liquid & steam) as a function of saturation temperature p ((Subscript[T, sat])). All water properties are defined in terms of the Helmholtz potential (Subscript[f, w]) (ρ, T) = e - T s : 

 p  =  (ρ^2)(Subscript[(((∂(Subscript[f, w]))/(∂ρ))), T]),  e=(Subscript[f, w])−T(Subscript[(((∂ (Subscript[f, w]))/(∂T))), ρ]),  h  =  (Subscript[f, w])−T (Subscript[(((∂ (Subscript[f, w]))/(∂T))), ρ])+ (Subscript[ρ (((∂ (Subscript[f, w]))/(∂ρ))), T]), s=−(Subscript[(((∂ (Subscript[f, w]))/(∂T))), ρ]), (Subscript[c, v])=(Subscript[(((∂ e)/(∂ T))),  v])

 (Subscript[c, p])=(Subscript[(((∂ h)/(∂T))), p]),  (Subscript[C, f])=(1/ρ)(Subscript[(((∂ρ)/(∂ p))), T]),  J=(Subscript[(((∂T)/(∂p))), h]) ,  (Subscript[v, s])=(Sqrt[(Subscript[(((∂ p)/(∂ρ))), s])]),  (Subscript[γ, f])=((− 1)/ρ)(Subscript[(((∂ρ)/(∂ T))), p]),  (Subscript[δ, f])=((Subscript[k, T])/(ρ (Subscript[c, p])))
 
UNITS: Helmholtz free energy (Subscript[f, W]) [kJ/kg], density ρ [kg/(m^3)], temperature T [K], pressure p [MPa], internal energy e [kJ/kg], enthalpy h [kJ/kg], entropy s [kJ/kg], specific isochoric heat cv [J/kg/K], specific isobaric heat cp [J/kg/K], compressibility [1/Pa], Joule coefficient [K/bar] speed of sound (Subscript[v, s]) [m/s], thermal expansivity (Subscript[γ, f])  [1/K], thermal diffusivity [(m^2)/s] and surface tension (Subscript[σ, w]) [N/m] (for two-phase water).

*(International Association for the Properties of Water and Steam, IAPWS R6-95,  Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use, 2016).

Introduction

●	In natural systems (aquifers, hydrocarbon & geothermal systems, atmosphere-ocean interactions), water is the most important acting fluid. 
★	Water in Geothermal Reservoirs covers a wide thermodynamic range, because the temperature gradient between the core and the surface of the Earth produces a continuous flow of natural heat, going from shallow depths ((10^2) m), to medium depths ((10^3) m), and to even deeper depths ((10^4) m), reaching the immense temperatures of molten rocks (T > 1000°C). Anomalous thermal gradient = 150°C/km.

          GS:  T ∈ [370, 5000] °C
Z ∈ [3, 6000] km

●	The fluids in deep reservoirs are at supercritical thermodynamic conditions (T >  375°C, p > 220 bar), with higher fluid density and more heat available. These systems could provide twenty times as much power, as the normal geothermal fluids used today (T <  360°C, p < 100 bar). This document provides a computational description of water, from very low temperatures and pressures up to supercritical conditions.

Molecular Structure of Water

♣	Water is a complex substance ((Subscript[H,   2])O). Two hydrogen atoms and one oxygen atom make an angle of approximately 104.5 degrees, giving an extraordinary stability to this molecular structure: each hydrogen atom is in line between the oxygen atom of its own molecule and the oxygen atom of a different molecule. 
     
										(Row[□,  Molecular structure of water ])

♣	These hydrogen bonds means extra attractive energy, causing many of the unusual properties of water : 

1. Large heat of vaporization, 2. Expansion upon freezing. 3. Low thermal conductivity. 4. Higher specific heat. 5. Very low liquid compressibility. 6. Very high steam compressibility. 7. Water is basic and acid at the same time (pKa/b = 15.74). 8. Water is a universal solvent. 9. Water reaches its maximum density at 4°C. 10. Below 4°C, water density declines and ice floats on liquid water.

The IAPWS-95 Formulation

●	The fundamental potential is the reduced specific Helmholtz free energy (dimensionless form of (Subscript[f, w]) =e - T s)

  	                    (((Subscript[f, w]) (ρ,T))/(R T))=ϕ (δ, τ)=(ϕ^0) (δ, τ) + (ϕ^r) (δ, τ) ; δ= (ρ/(Subscript[ρ, C])) , τ= ((Subscript[T, C])/T)  , Π= (p/(Subscript[P, C]))

(ϕ^0) is the ideal gas part and (ϕ^r) is the residual part of ϕ (δ, τ).

●	 Mass gas constant R = 0.46151805 [(kJ/(kg K))]; critical density (Subscript[ρ, C]) = 322 [(kg/(m^3))]; critical temperature (Subscript[T, C])= 647.096 K; critical pressure (Subscript[P, C]) = 220.64 bar. The Mathematica expressions of both functions are:

(ϕ^0)[δ_,τ_]:=Log[δ]+n0[1]+n0[2]τ+n0[3]Log[τ]+(∑_(i=4)^8)n0[i]Log[1-Exp[-γ0[[i]]τ]];
(ϕ^r)[δ_,τ_]:=(∑_(i=1)^7)n[i](δ^( d[i]))(τ^(t[i]))+(∑_(i=8)^51)n[i](δ^( d[i]))(τ^(t[i]))Exp[-(δ^( c[i]))]+ 
(∑_(i=52)^54)n[i](δ^( d[i]))(τ^(t[i]))Exp[-α[i](((δ-1))^2)-β[i](((τ-γ[i]))^2)]+(∑_(i=55)^56)n[i]((Δ[δ,τ])^(b[i]))δ\[Cross]ψ[δ,τ];  
(* this Cell needs numerical constants that are not included here *)

●	And the following derivatives:

(∂_(δ ))(ϕ^0)(δ,τ); (∂_(δ ))(ϕ^r)(δ,τ); (∂_(τ ))(ϕ^0)(δ,τ); (∂_(τ ))(ϕ^r)(δ,τ); (∂_(δδ ))(ϕ^0)(δ,τ); (∂_(δδ ))(ϕ^r)(δ,τ);
(∂_(ττ))(ϕ^0)(δ,τ); (∂_(ττ))(ϕ^r)(δ,τ); (∂_(δτ ))(ϕ^0) (δ, τ); (∂_(δτ ))(ϕ^r)(δ,τ); (∂_(δ ))Π; (∂_(Π ))δ

●	The two last derivatives are needed to calculate the backward inversion formulation of ϕ (δ, τ) ⟷ ϕ (Π, τ).

The Four Phases of Water in the IAPWS-95 Formulation

●	The IAPWS-95 covers the 4 water phases, in the entire stable fluid region for p ∈  [0.1, 10,000] bar, and for T ∈  [0, 1000] °C (1273 K), reproducing with high precision the experimental data available at the time the release 2016 was prepared, including  the critical point  (374°C, 221 bar, 322 kg/(m^3), 2085 kJ/kg). The interpolating equations of this formulation can be extrapolated for density and enthalpy beyond the already mentioned limits, giving reasonable accurate results for pressures up to (10^6) bar and temperatures up to 4727 °C (5000 K).

Helmholtz Free Energy Φ (δ,τ) for Water (IAPWS-95 Equations):

1. Pressure: p (δ,τ)=R (Subscript[ρ, C]) (Subscript[T, C]) (δ/τ)  (1 + δ  ((∂ (ϕ^r))/(∂ δ )))

2. Internal Energy:  e (δ,τ)=R (Subscript[T, C])   ( ((∂ (ϕ^0))/(∂ τ ))+((∂ (ϕ^r))/(∂ τ )))

3. Specific Enthalpy: h (δ,τ)=  R (Subscript[T, C]) ((1/τ) + ((∂ (ϕ^0))/(∂ τ ))+((∂ (ϕ^r))/(∂ τ )))

4. Specific Entropy: s (δ,τ)=  R  (τ (((∂ (ϕ^0))/(∂ τ ))+((∂ (ϕ^r))/(∂ τ )))-(ϕ^0)-(ϕ^r))

5. Isochoric Heat: (Subscript[C, v]) (δ,τ)=- R (τ^( 2)) ((((∂^( 2)) (ϕ^0))/(∂ (τ^( 2)) )) + (((∂^( 2)) (ϕ^r))/(∂ (τ^( 2)) )))

6. Isobaric Heat: (Subscript[C, p]) (δ,τ)=- R (τ^( 2)) ((((∂^( 2)) (ϕ^0))/(∂ (τ^( 2)) )) + (((∂^( 2)) (ϕ^r))/(∂ (τ^( 2)) )))+R ((((1 + δ(∂_δ)(ϕ^r)-  δ τ (∂_(δ τ))(ϕ^r)))^( 2))/(1 + 2δ(∂_δ)(ϕ^r)+ (δ^2) (∂_(δδ))(ϕ^r)(δ, τ)))

7. Bulk Modulus: (Subscript[K,  f]) (δ,τ)=((R (Subscript[ρ, C]) (Subscript[T, C]) )/τ) (δ+2(δ^2) ((∂ (ϕ^r))/(∂ δ ))(δ, τ)+(δ^( 3)) (((∂^( 2)) (ϕ^r))/(∂ (δ^( 2)) ))(δ, τ))

8. Isothermal Compressibility: (Subscript[C,  f]) (δ,τ)=((1 )/((Subscript[P, C])δ)) ((∂ δ)/(∂  Π )) (δ, τ) = (1/((Subscript[K,  f]) (δ,τ)))

9. Thermal Expansivity: (Subscript[γ,  f]) (δ,τ)=((τ  )/(Subscript[T, C])) (((1 + δ  (∂_δ)(ϕ^r)-  δ τ  (∂_(δ τ))(ϕ^r))/(1 +  2δ  (∂_δ)(ϕ^r)+  (δ^( 2))(∂_(δδ)) (ϕ^r))))

10. Thermal Conductivity: (Subscript[K,  T]) (δ,τ)=(Subscript[λ, 0 ])(τ)×(Subscript[λ, 1 ])(δ,τ)+(Subscript[λ,  2 ])(δ,τ)

11. Dynamic Viscosity: μ (δ,τ)=(Subscript[μ, 0 ])(δ,τ)×(Subscript[μ, 1 ])(δ,τ)

12. Speed of Sound: 
(Subscript[v,  s]) (δ,τ)=(Sqrt[R T (1 +  2δ ((∂ (ϕ^r))/(∂ δ ))+ (δ^( 2))(((∂^( 2)) (ϕ^r))/(∂ (δ^( 2)) ))-((((1 + δ  (∂_δ)(ϕ^r)-  δ τ  (∂_(δ τ))(ϕ^r)))^2)/((τ^( 2))((∂_(ττ)) (ϕ^0)+  (∂_(ττ)) (ϕ^r)))))(10^3)])

Extended IAPWS-95 and Backward Formulation Φ (Π,τ)

●	All functions and derivatives are computed as functions of (δ, τ) in IAPWS-95. In geothermal fields, for practical reasons, it is more efficient to use the field variables (p, T) because water density is not measured. A backward formulation entirely based on the IAPWS-95 was created using a two-dimensional projection of the Newton-Raphson algorithm in the following way (0 ⩽ error  ⩽ (10^(-9)), convergence in 2 or 3 iterations):

(Subscript[δ, n+1])= (Subscript[δ, n])-(((R (Subscript[ρ, C]) (Subscript[T, C]))/(Subscript[P, C]))((Subscript[δ, n]) +(Subscript[δ, n]^( 2))   ((∂ (ϕ^( r)))/(∂ (Subscript[δ, n]) )))-((Subscript[p, 0])/(Subscript[P, C])))  ((∂ (Subscript[δ, n]))/(∂ Π )) ⟶   ρ ((Subscript[p, 0]), (Subscript[T, 0])) ≈ (Subscript[ρ, C]) (Subscript[δ, n+1]) 

where:  ((∂ (Subscript[δ, n]))/(∂ Π )) =(((((∂ Π)/(∂ (Subscript[δ, n]) ))))^(-1)),and: ((∂ Π)/(∂ (Subscript[δ, n]) )) =((R (Subscript[ρ, C]) (Subscript[T, C]))/((Subscript[P, C]) τ))  (1+2(Subscript[δ, n])  ((∂ (ϕ^( r)))/(∂ (Subscript[δ, n]) ))+(Subscript[δ, n]^( 2))  (((∂^( 2)) (ϕ^( r)))/(∂ (Subscript[δ, n]^( 2)) )))

●	Both formulations were programmed within a Module in a single cell written in Mathematica 11.2 using the Wolfram Language: The accuracy of this code was verified using the numerical values of the corresponding Helmholtz function for water properties Φ = (ϕ^0) + (ϕ^( r)) and its partial derivatives, to test the results at T = 500 K, and ρ = 838.025 kg (m^(-3)). 

●	VERIFICATION: The corresponding code presented herein reproduces exactly the proposed values published by the IAPWS-95:

Wolfram code for the IAPWS-95 Formulation

●	The Mathematica code was created using the Wolfram Programming Language 11.2 on the basis of the Helmholtz potential previously described. The main module has the general form:

IAPWS95[ρ_,T_]:= Module[{δ=(ρ/322),τ=(647.096/T)},(* Program *) ... ] Displays the options in English or Spanish:

 PROGRAM OVERVIEW:

1. Approximate computation of density (IFC67)

2. Module IAPWS-95 (body of main code).

3. Definition of constants, (ϕ^0),(ϕ^r) and their derivatives.

4. Computation of single-phase thermodynamic properties.

5. Computation of two-phase thermodynamic properties. 

6. Display of results in tabular or in graphical form.

Wolfram 4D Graphics: IAPWS-95 Helmholtz potential (ρ, T, p)

●	The Mathematica code has 2D and 3D graphic capabilities. Actually, the numerical values of the Helmholtz potential are in 4D, because it depends on density, temperature and pressure  (Subscript[f, W]) (ρ, T, p). The next graphics show the 3D projection of the Helmholtz hypersurface in two pressure hyperplanes at 1 bar and 500 bar respectively, in the range ρ ∈ [0,1000] kg/(m^(3 ))and T ∈  [0, 1000] °C.

3D Graphics: The three regions of Geothermal Water T (h, p)

●	This figure is the 3D surface of temperature (°C) as function of enthalpy and pressure for T (h, p) ∈  [0, 800] °C, covering the most common geothermal range for practical use. The curved region below the critical point in this 3D surface corresponds to the two-phase zone of water. The liquid region is located to the left of this zone, the superheated steam region is above it, and the supercritical zone is up to the right.

●	The red ball shows the approximate location of the critical point in the plane (h, T).


                    T (°C)

CONCLUSIONS

●	The IAPWS-95 mathematical formulation for water is programmed within a Mathematica-11.2 code, written in the Wolfram Programming Language. The computable properties are: p, h, e, s, (Subscript[c, p]), (Subscript[c, v]), (Subscript[v, s]), functions of (Subscript[f, W])(ρ, T), and (Subscript[σ, w]) (only for two-phase p((Subscript[T, sat])) fluid).

●	This formulation has been verified experimentally and is strictly valid in the range p ∈ [0.01,1000] MPa and T ∈  [0, 1000] °C, covering geothermal water under magmatic conditions.

●	The IAPWS-95 formulation can be extrapolated reasonably well to cover an extra range of p ∈ [(10^3), (10^5)] MPa and T ∈  [1000, 4727] °C. Outside these intervals, the IAPWS-95 formulation is no longer valid.

●	This Mathematica code is an extended version of the IAPWS-95 original formulation, and includes (Subscript[k, T]), μ, (Subscript[C, f]), (Subscript[K, f]) and (Subscript[γ, f]) of water. It contains efficient and accurate algorithms to perform backward inversions of internal correlations to compute single-phase properties as functions of pressure (bar) and temperature (°C), and two-phase properties as functions of saturation pressure T((Subscript[p, sat])). It includes 2D and 3D graphic capabilities.

●	This is a very powerful computational program that uses many of the advanced features of the Wolfram programming language. The code is of fundamental importance in Earth Sciences, for the steam industry and geothermal technology, where the main working fluid is water. 

●	The code can be used accurately for all kind of geothermal reservoirs of low, medium, and high enthalpy, as well as in supercritical systems.

Section head (*Section*)

Lorem ipsum dolor sit amet, consectetuer adipiscing elit. Pellentesque enim dolor, euismod sit amet, rutrum nec, fermentum id, risus. (*Text*)