Source code for iapws._iapws

#!/usr/bin/python
# -*- coding: utf-8 -*-
# pylint: disable=invalid-name
# pylint: disable=too-many-lines, too-many-locals, too-many-statements
# pylint: disable=too-many-branches, too-many-boolean-expressions

"""
Miscelaneous IAPWS standards. This module include:

    * :func:`_Ice`: Ice Ih state equation
    * :func:`_Liquid`: Properties of liquid water at 0.1 MPa
    * :func:`_Supercooled`: Thermodynamic properties of supercooled water
    * :func:`_Sublimation_Pressure`: Sublimation pressure correlation
    * :func:`_Melting_Pressure`: Melting pressure correlation
    * :func:`_Viscosity`: Viscosity correlation
    * :func:`_ThCond`: Themal conductivity correlation
    * :func:`_Tension`: Surface tension correlation
    * :func:`_Dielectric`: Dielectric constant correlation
    * :func:`_Refractive`: Refractive index correlation
    * :func:`_Kw`: Ionization constant correlation for ordinary water
    * :func:`_Conductivity`: Electrolytic conductivity correlation
    * :func:`_D2O_Viscosity`: Viscosity correlation for heavy water
    * :func:`_D2O_ThCond`: Thermal conductivity correlation for heavy water
    * :func:`_D2O_Tension`: Surface tension correlation for heavy water
    * :func:`_D2O_Sublimation_Pressure`: Sublimation Pressure correlation
      for heavy water
    * :func:`_D2O_Melting_Pressure`: Melting Pressure correlation for heavy
      water
    * :func:`_Henry`: Henry constant for liquid-gas equilibrium
    * :func:`_Kvalue`: Vapor-liquid distribution constant
"""

from __future__ import division

from cmath import log as log_c
from math import log, exp, tan, atan, acos, sin, pi, log10
import warnings

from scipy.optimize import newton


# Constants
M = 18.015268     # g/mol
R = 0.461526      # kJ/kg·K

# Table 1 from Release on the Values of Temperature, Pressure and Density of
# Ordinary and Heavy Water Substances at their Respective Critical Points
Tc = 647.096      # K
Pc = 22.064       # MPa
rhoc = 322.       # kg/m³
Tc_D2O = 643.847  # K
Pc_D2O = 21.6618   # MPa
rhoc_D2O = 355.9999698294    # kg/m³

Tt = 273.16       # K
Pt = 611.657e-6   # MPa
Tb = 373.1243     # K
f_acent = 0.3443

# IAPWS, Guideline on the Use of Fundamental Physical Constants and Basic
# Constants of Water, http://www.iapws.org/relguide/fundam.pdf
Dipole = 1.85498  # Debye


# IAPWS-06 for Ice
[docs]def _Ice(T, P): """Basic state equation for Ice Ih Parameters ---------- T : float Temperature, [K] P : float Pressure, [MPa] Returns ------- prop : dict Dict with calculated properties of ice. The available properties are: * rho: Density, [kg/m³] * h: Specific enthalpy, [kJ/kg] * u: Specific internal energy, [kJ/kg] * a: Specific Helmholtz energy, [kJ/kg] * g: Specific Gibbs energy, [kJ/kg] * s: Specific entropy, [kJ/kgK] * cp: Specific isobaric heat capacity, [kJ/kgK] * alfav: Cubic expansion coefficient, [1/K] * beta: Pressure coefficient, [MPa/K] * xkappa: Isothermal compressibility, [1/MPa] * ks: Isentropic compressibility, [1/MPa] * gt: [∂g/∂T]P * gtt: [∂²g/∂T²]P * gp: [∂g/∂P]T * gpp: [∂²g/∂P²]T * gtp: [∂²g/∂T∂P] Notes ----- Raise :class:`NotImplementedError` if input isn't in limit: * T ≤ 273.16 * P ≤ 208.566 * State below the melting and sublimation lines Examples -------- >>> st1 = _Ice(100, 100) >>> st1["rho"], st1["h"], st1["s"] 941.678203297 -483.491635676 -2.61195122589 >>> st2 = _Ice(273.152519,0.101325) >>> st2["a"], st2["u"], st2["cp"] -0.00918701567 -333.465403393 2.09671391024 >>> st3 = _Ice(273.16,611.657e-6) >>> st3["alfav"], st3["beta"], st3["xkappa"], st3["ks"] 0.000159863102566 1.35714764659 1.17793449348e-04 1.14161597779e-04 References ---------- IAPWS, Revised Release on the Equation of State 2006 for H2O Ice Ih September 2009, http://iapws.org/relguide/Ice-2009.html """ # Check input in range of validity if T > 273.16: # No Ice Ih stable warnings.warn("Metastable ice") elif P > 208.566: # Ice Ih limit upper pressure raise NotImplementedError("Incoming out of bound") elif P < Pt: Psub = _Sublimation_Pressure(T) if Psub > P: # Zone Gas warnings.warn("Metastable ice in vapor region") elif T > 251.165: Pmel = _Melting_Pressure(T) if Pmel < P: # Zone Liquid warnings.warn("Metastable ice in liquid region") Tr = T/Tt Pr = P/Pt P0 = 101325e-6/Pt s0 = -0.332733756492168e4*1e-3 # Express in kJ/kgK gok = [-0.632020233335886e6, 0.655022213658955, -0.189369929326131e-7, 0.339746123271053e-14, -0.556464869058991e-21] r2k = [complex(-0.725974574329220e2, -0.781008427112870e2)*1e-3, complex(-0.557107698030123e-4, 0.464578634580806e-4)*1e-3, complex(0.234801409215913e-10, -0.285651142904972e-10)*1e-3] t1 = complex(0.368017112855051e-1, 0.510878114959572e-1) t2 = complex(0.337315741065416, 0.335449415919309) r1 = complex(0.447050716285388e2, 0.656876847463481e2)*1e-3 go = gop = gopp = 0 for k in range(5): go += gok[k]*1e-3*(Pr-P0)**k for k in range(1, 5): gop += gok[k]*1e-3*k/Pt*(Pr-P0)**(k-1) for k in range(2, 5): gopp += gok[k]*1e-3*k*(k-1)/Pt**2*(Pr-P0)**(k-2) r2 = r2p = 0 for k in range(3): r2 += r2k[k]*(Pr-P0)**k for k in range(1, 3): r2p += r2k[k]*k/Pt*(Pr-P0)**(k-1) r2pp = r2k[2]*2/Pt**2 c = r1*((t1-Tr)*log_c(t1-Tr)+(t1+Tr)*log_c(t1+Tr)-2*t1*log_c( t1)-Tr**2/t1)+r2*((t2-Tr)*log_c(t2-Tr)+(t2+Tr)*log_c( t2+Tr)-2*t2*log_c(t2)-Tr**2/t2) ct = r1*(-log_c(t1-Tr)+log_c(t1+Tr)-2*Tr/t1)+r2*( -log_c(t2-Tr)+log_c(t2+Tr)-2*Tr/t2) ctt = r1*(1/(t1-Tr)+1/(t1+Tr)-2/t1) + r2*(1/(t2-Tr)+1/(t2+Tr)-2/t2) cp = r2p*((t2-Tr)*log_c(t2-Tr)+(t2+Tr)*log_c( t2+Tr)-2*t2*log_c(t2)-Tr**2/t2) ctp = r2p*(-log_c(t2-Tr)+log_c(t2+Tr)-2*Tr/t2) cpp = r2pp*((t2-Tr)*log_c(t2-Tr)+(t2+Tr)*log_c( t2+Tr)-2*t2*log_c(t2)-Tr**2/t2) g = go-s0*Tt*Tr+Tt*c.real gt = -s0+ct.real gp = gop+Tt*cp.real gtt = ctt.real/Tt gtp = ctp.real gpp = gopp+Tt*cpp.real propiedades = {} propiedades["gt"] = gt propiedades["gp"] = gp propiedades["gtt"] = gtt propiedades["gpp"] = gpp propiedades["gtp"] = gtp propiedades["T"] = T propiedades["P"] = P propiedades["v"] = gp/1000 propiedades["rho"] = 1000./gp propiedades["h"] = g-T*gt propiedades["s"] = -gt propiedades["cp"] = -T*gtt propiedades["u"] = g-T*gt-P*gp propiedades["g"] = g propiedades["a"] = g-P*gp propiedades["alfav"] = gtp/gp propiedades["beta"] = -gtp/gpp propiedades["xkappa"] = -gpp/gp propiedades["ks"] = (gtp**2-gtt*gpp)/gp/gtt return propiedades
# IAPWS-08 for Liquid water at 0.1 MPa
[docs]def _Liquid(T, P=0.1): """Supplementary release on properties of liquid water at 0.1 MPa Parameters ---------- T : float Temperature, [K] P : float Pressure, [MPa] Although this relation is for P=0.1MPa, can be extrapoled at pressure 0.3 MPa Returns ------- prop : dict Dict with calculated properties of water. The available properties are: * h: Specific enthalpy, [kJ/kg] * u: Specific internal energy, [kJ/kg] * a: Specific Helmholtz energy, [kJ/kg] * g: Specific Gibbs energy, [kJ/kg] * s: Specific entropy, [kJ/kgK] * cp: Specific isobaric heat capacity, [kJ/kgK] * cv: Specific isochoric heat capacity, [kJ/kgK] * w: Speed of sound, [m/s²] * rho: Density, [kg/m³] * v: Specific volume, [m³/kg] * vt: [∂v/∂T]P, [m³/kgK] * vtt: [∂²v/∂T²]P, [m³/kgK²] * vp: [∂v/∂P]T, [m³/kg/MPa] * vtp: [∂²v/∂T∂P], [m³/kg/MPa] * alfav: Cubic expansion coefficient, [1/K] * xkappa : Isothermal compressibility, [1/MPa] * ks: Isentropic compressibility, [1/MPa] * mu: Viscosity, [Pas] * k: Thermal conductivity, [W/mK] * epsilon: Dielectric constant, [-] Notes ----- Raise :class:`NotImplementedError` if input isn't in limit: * 253.15 ≤ T ≤ 383.15 * 0.1 ≤ P ≤ 0.3 Examples -------- >>> st1 = _Liquid(260) >>> st1["rho"], st1["h"], st1["s"] 997.0683602710492 -55.86223174460868 -0.20998554842619535 References ---------- IAPWS, Revised Supplementary Release on Properties of Liquid Water at 0.1 MPa, http://www.iapws.org/relguide/LiquidWater.html """ # Check input in range of validity if T <= 253.15 or T >= 383.15 or P < 0.1 or P > 0.3: raise NotImplementedError("Incoming out of bound") if P != 0.1: # Raise a warning if the P value is extrapolated warnings.warn("Using extrapolated values") Rg = 0.46151805 # kJ/kgK Po = 0.1 Tr = 10 tau = T/Tr alfa = Tr/(593-T) beta = Tr/(T-232) a = [None, -1.661470539e5, 2.708781640e6, -1.557191544e8, None, 1.93763157e-2, 6.74458446e3, -2.22521604e5, 1.00231247e8, -1.63552118e9, 8.32299658e9, -7.5245878e-6, -1.3767418e-2, 1.0627293e1, -2.0457795e2, 1.2037414e3] b = [None, -8.237426256e-1, 1.908956353, -2.017597384, 8.546361348e-1, 5.78545292e-3, -1.53195665E-2, 3.11337859e-2, -4.23546241e-2, 3.38713507e-2, -1.19946761e-2, -3.1091470e-6, 2.8964919e-5, -1.3112763e-4, 3.0410453e-4, -3.9034594e-4, 2.3403117e-4, -4.8510101e-5] c = [None, -2.452093414e2, 3.869269598e1, -8.983025854] n = [None, 4, 5, 7, None, None, 4, 5, 7, 8, 9, 1, 3, 5, 6, 7] m = [None, 2, 3, 4, 5, 1, 2, 3, 4, 5, 6, 1, 3, 4, 5, 6, 7, 9] suma1 = sum(a[i]*alfa**n[i] for i in range(1, 4)) suma2 = sum(b[i]*beta**m[i] for i in range(1, 5)) go = Rg*Tr*(c[1]+c[2]*tau+c[3]*tau*log(tau)+suma1+suma2) suma1 = sum(a[i]*alfa**n[i] for i in range(6, 11)) suma2 = sum(b[i]*beta**m[i] for i in range(5, 11)) vo = Rg*Tr/Po/1000*(a[5]+suma1+suma2) suma1 = sum(a[i]*alfa**n[i] for i in range(11, 16)) suma2 = sum(b[i]*beta**m[i] for i in range(11, 18)) vpo = Rg*Tr/Po**2/1000*(suma1+suma2) suma1 = sum(n[i]*a[i]*alfa**(n[i]+1) for i in range(1, 4)) suma2 = sum(m[i]*b[i]*beta**(m[i]+1) for i in range(1, 5)) so = -Rg*(c[2]+c[3]*(1+log(tau))+suma1-suma2) suma1 = sum(n[i]*(n[i]+1)*a[i]*alfa**(n[i]+2) for i in range(1, 4)) suma2 = sum(m[i]*(m[i]+1)*b[i]*beta**(m[i]+2) for i in range(1, 5)) cpo = -Rg*(c[3]+tau*suma1+tau*suma2) suma1 = sum(n[i]*a[i]*alfa**(n[i]+1) for i in range(6, 11)) suma2 = sum(m[i]*b[i]*beta**(m[i]+1) for i in range(5, 11)) vto = Rg/Po/1000*(suma1-suma2) # This properties are only neccessary for computing thermodynamic # properties at pressures different from 0.1 MPa suma1 = sum(n[i]*(n[i]+1)*a[i]*alfa**(n[i]+2) for i in range(6, 11)) suma2 = sum(m[i]*(m[i]+1)*b[i]*beta**(m[i]+2) for i in range(5, 11)) vtto = Rg/Tr/Po/1000*(suma1+suma2) suma1 = sum(n[i]*a[i]*alfa**(n[i]+1) for i in range(11, 16)) suma2 = sum(m[i]*b[i]*beta**(m[i]+1) for i in range(11, 18)) vpto = Rg/Po**2/1000*(suma1-suma2) if P != 0.1: go += vo*(P-0.1) so -= vto*(P-0.1) cpo -= T*vtto*(P-0.1) vo -= vpo*(P-0.1) vto += vpto*(P-0.1) vppo = 3.24e-10*Rg*Tr/0.1**3 vpo += vppo*(P-0.1) h = go+T*so u = h-P*vo a = go-P*vo cv = cpo+T*vto**2/vpo xkappa = -vpo/vo alfa = vto/vo ks = -(T*vto**2/cpo+vpo)/vo w = (-vo**2*1e9/(vpo*1e3+T*vto**2*1e6/cpo))**0.5 propiedades = {} propiedades["g"] = go propiedades["T"] = T propiedades["P"] = P propiedades["v"] = vo propiedades["vt"] = vto propiedades["vp"] = vpo propiedades["vpt"] = vpto propiedades["vtt"] = vtto propiedades["rho"] = 1/vo propiedades["h"] = h propiedades["s"] = so propiedades["cp"] = cpo propiedades["cv"] = cv propiedades["u"] = u propiedades["a"] = a propiedades["xkappa"] = xkappa propiedades["alfav"] = vto/vo propiedades["ks"] = ks propiedades["w"] = w # Viscosity correlation, Eq 7 a = [None, 280.68, 511.45, 61.131, 0.45903] b = [None, -1.9, -7.7, -19.6, -40] T_ = T/300 mu = sum(a[i]*T_**b[i] for i in range(1, 5))/1e6 propiedades["mu"] = mu # Thermal conductivity correlation, Eq 8 c = [None, 1.6630, -1.7781, 1.1567, -0.432115] d = [None, -1.15, -3.4, -6.0, -7.6] k = sum(c[i]*T_**d[i] for i in range(1, 5)) propiedades["k"] = k # Dielectric constant correlation, Eq 9 e = [None, -43.7527, 299.504, -399.364, 221.327] f = [None, -0.05, -1.47, -2.11, -2.31] epsilon = sum(e[i]*T_**f[i] for i in range(1, 5)) propiedades["epsilon"] = epsilon return propiedades
# IAPWS-15 for supercooled liquid water
[docs]def _Supercooled(T, P): """Guideline on thermodynamic properties of supercooled water Parameters ---------- T : float Temperature, [K] P : float Pressure, [MPa] Returns ------- prop : dict Dict with calculated properties of water. The available properties are: * L: Ordering field, [-] * x: Mole fraction of low-density structure, [-] * rho: Density, [kg/m³] * s: Specific entropy, [kJ/kgK] * h: Specific enthalpy, [kJ/kg] * u: Specific internal energy, [kJ/kg] * a: Specific Helmholtz energy, [kJ/kg] * g: Specific Gibbs energy, [kJ/kg] * alfap: Thermal expansion coefficient, [1/K] * xkappa : Isothermal compressibility, [1/MPa] * cp: Specific isobaric heat capacity, [kJ/kgK] * cv: Specific isochoric heat capacity, [kJ/kgK] * w: Speed of sound, [m/s²] Notes ----- Raise :class:`NotImplementedError` if input isn't in limit: * Tm ≤ T ≤ 300 * 0 < P ≤ 1000 The minimum temperature in range of validity is the melting temperature, it depend of pressure Raise :class:`RuntimeError` if solution isn't founded Examples -------- >>> liq = _supercooled(235.15, 0.101325) >>> liq["rho"], liq["cp"], liq["w"] 968.09999 5.997563 1134.5855 References ---------- iapws, guideline on thermodynamic properties of supercooled water, http://iapws.org/relguide/Supercooled.html """ # Check input in range of validity if P < 198.9: Tita = T/235.15 Ph = 0.1+228.27*(1-Tita**6.243)+15.724*(1-Tita**79.81) if P < Ph or T > 300: raise NotImplementedError("Incoming out of bound") else: Th = 172.82+0.03718*P+3.403e-5*P**2-1.573e-8*P**3 if T < Th or T > 300 or P > 1000: raise NotImplementedError("Incoming out of bound") # Parameters, Table 1 Tll = 228.2 rho0 = 1081.6482 Rg = 0.461523087 pi0 = 300e3/rho0/Rg/Tll omega0 = 0.5212269 L0 = 0.76317954 k0 = 0.072158686 k1 = -0.31569232 k2 = 5.2992608 # Reducing parameters, Eq 2 tau = T/Tll-1 p = P*1000/rho0/Rg/Tll tau_ = tau+1 p_ = p+pi0 # Eq 3 ci = [-8.1570681381655, 1.2875032, 7.0901673598012, -3.2779161e-2, 7.3703949e-1, -2.1628622e-1, -5.1782479, 4.2293517e-4, 2.3592109e-2, 4.3773754, -2.9967770e-3, -9.6558018e-1, 3.7595286, 1.2632441, 2.8542697e-1, -8.5994947e-1, -3.2916153e-1, 9.0019616e-2, 8.1149726e-2, -3.2788213] ai = [0, 0, 1, -0.2555, 1.5762, 1.6400, 3.6385, -0.3828, 1.6219, 4.3287, 3.4763, 5.1556, -0.3593, 5.0361, 2.9786, 6.2373, 4.0460, 5.3558, 9.0157, 1.2194] bi = [0, 1, 0, 2.1051, 1.1422, 0.9510, 0, 3.6402, 2.0760, -0.0016, 2.2769, 0.0008, 0.3706, -0.3975, 2.9730, -0.3180, 2.9805, 2.9265, 0.4456, 0.1298] di = [0, 0, 0, -0.0016, 0.6894, 0.0130, 0.0002, 0.0435, 0.0500, 0.0004, 0.0528, 0.0147, 0.8584, 0.9924, 1.0041, 1.0961, 1.0228, 1.0303, 1.6180, 0.5213] phir = phirt = phirp = phirtt = phirtp = phirpp = 0 for c, a, b, d in zip(ci, ai, bi, di): phir += c*tau_**a*p_**b*exp(-d*p_) phirt += c*a*tau_**(a-1)*p_**b*exp(-d*p_) phirp += c*tau_**a*p_**(b-1)*(b-d*p_)*exp(-d*p_) phirtt += c*a*(a-1)*tau_**(a-2)*p_**b*exp(-d*p_) phirtp += c*a*tau_**(a-1)*p_**(b-1)*(b-d*p_)*exp(-d*p_) phirpp += c*tau_**a*p_**(b-2)*((d*p_-b)**2-b)*exp(-d*p_) # Eq 5 K1 = ((1+k0*k2+k1*(p-k2*tau))**2-4*k0*k1*k2*(p-k2*tau))**0.5 K2 = (1+k2**2)**0.5 # Eq 6 omega = 2+omega0*p # Eq 4 L = L0*K2/2/k1/k2*(1+k0*k2+k1*(p+k2*tau)-K1) # Define interval of solution, Table 4 if omega < 10/9*(log(19)-L): xmin = 0.049 xmax = 0.5 elif 10/9*(log(19)-L) <= omega < 50/49*(log(99)-L): xmin = 0.0099 xmax = 0.051 else: xmin = 0.99*exp(-50/49*L-omega) xmax = min(1.1*exp(-L-omega), 0.0101) # Eq 8 def f(x): "Function for iterative calculation" if x < xmin: x = xmin if x > xmax: x = xmax return L+log(x/(1-x))+omega*(1-2*x) x = None for xo in (xmin, xmax, (xmin+xmax)/2): try: x, sol = newton(f, xo, full_output=True) except RuntimeError: pass else: if sol.converged: break # Exit when solution don't found if not x: raise RuntimeError("Solution don't found") # Eq 12 fi = 2*x-1 Xi = 1/(2/(1-fi**2)-omega) # Derivatives, Table 3 Lt = L0*K2/2*(1+(1-k0*k2+k1*(p-k2*tau))/K1) Lp = L0*K2*(K1+k0*k2-k1*p+k1*k2*tau-1)/2/k2/K1 Ltt = -2*L0*K2*k0*k1*k2**2/K1**3 Ltp = 2*L0*K2*k0*k1*k2/K1**3 Lpp = -2*L0*K2*k0*k1/K1**3 prop = {} prop["L"] = L prop["x"] = x # Eq 13 prop["rho"] = rho0/((tau+1)/2*(omega0/2*(1-fi**2)+Lp*(fi+1))+phirp) # Eq 1 prop["g"] = phir+(tau+1)*(x*L+x*log(x)+(1-x)*log(1-x)+omega*x*(1-x)) # Eq 14 prop["s"] = -Rg*((tau+1)/2*Lt*(fi+1) + (x*L+x*log(x)+(1-x)*log(1-x)+omega*x*(1-x))+phirt) # Basic derived state properties prop["h"] = prop["g"]+T*prop["s"] prop["u"] = prop["h"]+P/prop["rho"] prop["a"] = prop["u"]-T*prop["s"] # Eq 15 prop["xkappa"] = prop["rho"]/rho0**2/Rg*1000/Tll*( (tau+1)/2*(Xi*(Lp-omega0*fi)**2-(fi+1)*Lpp)-phirpp) prop["alfap"] = prop["rho"]/rho0/Tll*( Ltp/2*(tau+1)*(fi+1) + (omega0*(1-fi**2)/2+Lp*(fi+1))/2 - (tau+1)*Lt/2*Xi*(Lp-omega0*fi) + phirtp) prop["cp"] = -Rg*(tau+1)*(Lt*(fi+1)+(tau+1)/2*(Ltt*(fi+1)-Lt**2*Xi)+phirtt) # Eq 16 prop["cv"] = prop["cp"]-T*prop["alfap"]**2/prop["rho"]/prop["xkappa"]*1e3 # Eq 17 prop["w"] = (prop["rho"]*prop["xkappa"]*1e-6*prop["cv"]/prop["cp"])**-0.5 return prop
[docs]def _Sublimation_Pressure(T): """Sublimation Pressure correlation Parameters ---------- T : float Temperature, [K] Returns ------- P : float Pressure at sublimation line, [MPa] Notes ----- Raise :class:`NotImplementedError` if input isn't in limit: * 50 ≤ T ≤ 273.16 Examples -------- >>> _Sublimation_Pressure(230) 8.947352740189152e-06 References ---------- IAPWS, Revised Release on the Pressure along the Melting and Sublimation Curves of Ordinary Water Substance, http://iapws.org/relguide/MeltSub.html. """ if 50 <= T <= 273.16: Tita = T/Tt suma = 0 a = [-0.212144006e2, 0.273203819e2, -0.61059813e1] expo = [0.333333333e-2, 1.20666667, 1.70333333] for ai, expi in zip(a, expo): suma += ai*Tita**expi return exp(suma/Tita)*Pt raise NotImplementedError("Incoming out of bound")
[docs]def _Melting_Pressure(T, ice="Ih"): """Melting Pressure correlation Parameters ---------- T : float Temperature, [K] ice: string Type of ice: Ih, III, V, VI, VII. Below 273.15 is a mandatory input, the ice Ih is the default value. Above 273.15, the ice type is unnecesary. Returns ------- P : float Pressure at sublimation line, [MPa] Notes ----- Raise :class:`NotImplementedError` if input isn't in limit: * 251.165 ≤ T ≤ 715 Examples -------- >>> _Melting_Pressure(260) 8.947352740189152e-06 >>> _Melting_Pressure(254, "III") 268.6846466336108 References ---------- IAPWS, Revised Release on the Pressure along the Melting and Sublimation Curves of Ordinary Water Substance, http://iapws.org/relguide/MeltSub.html. """ if ice == "Ih" and 251.165 <= T <= 273.16: # Ice Ih Tref = Tt Pref = Pt Tita = T/Tref a = [0.119539337e7, 0.808183159e5, 0.33382686e4] expo = [3., 0.2575e2, 0.10375e3] suma = 1 for ai, expi in zip(a, expo): suma += ai*(1-Tita**expi) P = suma*Pref elif ice == "III" and 251.165 < T <= 256.164: # Ice III Tref = 251.165 Pref = 208.566 Tita = T/Tref P = Pref*(1-0.299948*(1-Tita**60.)) elif (ice == "V" and 256.164 < T <= 273.15) or 273.15 < T <= 273.31: # Ice V Tref = 256.164 Pref = 350.100 Tita = T/Tref P = Pref*(1-1.18721*(1-Tita**8.)) elif 273.31 < T <= 355: # Ice VI Tref = 273.31 Pref = 632.400 Tita = T/Tref P = Pref*(1-1.07476*(1-Tita**4.6)) elif 355. < T <= 715: # Ice VII Tref = 355 Pref = 2216.000 Tita = T/Tref P = Pref*exp(1.73683*(1-1./Tita)-0.544606e-1*(1-Tita**5) + 0.806106e-7*(1-Tita**22)) else: raise NotImplementedError("Incoming out of bound") return P
# Transport properties
[docs]def _Viscosity(rho, T, fase=None, drho=None): """Equation for the Viscosity Parameters ---------- rho : float Density, [kg/m³] T : float Temperature, [K] fase: dict, optional for calculate critical enhancement phase properties drho: float, optional for calculate critical enhancement [∂ρ/∂P]T at reference state, Returns ------- μ : float Viscosity, [Pa·s] Examples -------- >>> _Viscosity(998, 298.15) 0.0008897351001498108 >>> _Viscosity(600, 873.15) 7.743019522728247e-05 References ---------- IAPWS, Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance, http://www.iapws.org/relguide/viscosity.html """ Tr = T/Tc Dr = rho/rhoc # Eq 11 H = [1.67752, 2.20462, 0.6366564, -0.241605] mu0 = 100*Tr**0.5/sum(Hi/Tr**i for i, Hi in enumerate(H)) # Eq 12 li = [0, 1, 2, 3, 0, 1, 2, 3, 5, 0, 1, 2, 3, 4, 0, 1, 0, 3, 4, 3, 5] lj = [0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 4, 4, 5, 6, 6] Hij = [0.520094, 0.850895e-1, -0.108374e1, -0.289555, 0.222531, 0.999115, 0.188797e1, 0.126613e1, 0.120573, -0.281378, -0.906851, -0.772479, -0.489837, -0.257040, 0.161913, 0.257399, -0.325372e-1, 0.698452e-1, 0.872102e-2, -0.435673e-2, -0.593264e-3] mu1 = exp(Dr*sum((1/Tr-1)**i*h*(Dr-1)**j for i, j, h in zip(li, lj, Hij))) # Critical enhancement if fase and drho: qc = 1/1.9 qd = 1/1.1 # Eq 21 DeltaX = Pc*Dr**2*(fase.drhodP_T/rho-drho/rho*1.5/Tr) if DeltaX < 0: DeltaX = 0 # Eq 20 X = 0.13*(DeltaX/0.06)**(0.63/1.239) if X <= 0.3817016416: # Eq 15 Y = qc/5*X*(qd*X)**5*(1-qc*X+(qc*X)**2-765./504*(qd*X)**2) else: Fid = acos((1+qd**2*X**2)**-0.5) # Eq 17 w = abs((qc*X-1)/(qc*X+1))**0.5*tan(Fid/2) # Eq 19 # Eq 18 if qc*X > 1: Lw = log((1+w)/(1-w)) else: Lw = 2*atan(abs(w)) # Eq 16 Y = sin(3*Fid)/12-sin(2*Fid)/4/qc/X+(1-5/4*(qc*X)**2)/( qc*X)**2*sin(Fid)-((1-3/2*(qc*X)**2)*Fid-abs(( qc*X)**2-1)**1.5*Lw)/(qc*X)**3 # Eq 14 mu2 = exp(0.068*Y) else: mu2 = 1 # Eq 10 mu = mu0*mu1*mu2 return mu*1e-6
[docs]def _ThCond(rho, T, fase=None, drho=None): """Equation for the thermal conductivity Parameters ---------- rho : float Density, [kg/m³] T : float Temperature, [K] fase: dict, optional for calculate critical enhancement phase properties drho: float, optional for calculate critical enhancement [∂ρ/∂P]T at reference state, Returns ------- k : float Thermal conductivity, [W/mK] Examples -------- >>> _ThCond(998, 298.15) 0.6077128675880629 >>> _ThCond(0, 873.15) 0.07910346589648833 References ---------- IAPWS, Release on the IAPWS Formulation 2011 for the Thermal Conductivity of Ordinary Water Substance, http://www.iapws.org/relguide/ThCond.html """ d = rho/rhoc Tr = T/Tc # Eq 16 no = [2.443221e-3, 1.323095e-2, 6.770357e-3, -3.454586e-3, 4.096266e-4] k0 = Tr**0.5/sum(n/Tr**i for i, n in enumerate(no)) # Eq 17 li = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4] lj = [0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 0, 1, 2, 3, 4, 5] nij = [1.60397357, -0.646013523, 0.111443906, 0.102997357, -0.0504123634, 0.00609859258, 2.33771842, -2.78843778, 1.53616167, -0.463045512, 0.0832827019, -0.00719201245, 2.19650529, -4.54580785, 3.55777244, -1.40944978, 0.275418278, -0.0205938816, -1.21051378, 1.60812989, -0.621178141, 0.0716373224, -2.7203370, 4.57586331, -3.18369245, 1.1168348, -0.19268305, 0.012913842] k1 = exp(d*sum((1/Tr-1)**i*n*(d-1)**j for i, j, n in zip(li, lj, nij))) # Critical enhancement if fase: Rg = 0.46151805 if not drho: # Industrial formulation # Eq 25 if d <= 0.310559006: ai = [6.53786807199516, -5.61149954923348, 3.39624167361325, -2.27492629730878, 10.2631854662709, 1.97815050331519] elif d <= 0.776397516: ai = [6.52717759281799, -6.30816983387575, 8.08379285492595, -9.82240510197603, 12.1358413791395, -5.54349664571295] elif d <= 1.242236025: ai = [5.35500529896124, -3.96415689925446, 8.91990208918795, -12.0338729505790, 9.19494865194302, -2.16866274479712] elif d <= 1.863354037: ai = [1.55225959906681, 0.464621290821181, 8.93237374861479, -11.0321960061126, 6.16780999933360, -0.965458722086812] else: ai = [1.11999926419994, 0.595748562571649, 9.88952565078920, -10.3255051147040, 4.66861294457414, -0.503243546373828] drho = 1/sum(a*d**i for i, a in enumerate(ai))*rhoc/Pc DeltaX = d*(Pc/rhoc*fase.drhodP_T-Pc/rhoc*drho*1.5/Tr) if DeltaX < 0: DeltaX = 0 X = 0.13*(DeltaX/0.06)**(0.63/1.239) # Eq 22 y = X/0.4 # Eq 20 # Eq 19 if y < 1.2e-7: Z = 0 else: Z = 2/pi/y*(((1-1/fase.cp_cv)*atan(y)+y/fase.cp_cv)-( 1-exp(-1/(1/y+y**2/3/d**2)))) # Eq 18 k2 = 177.8514*d*fase.cp/Rg*Tr/fase.mu*1e-6*Z else: # No critical enhancement k2 = 0 # Eq 10 k = k0*k1+k2 return 1e-3*k
[docs]def _Tension(T): """Equation for the surface tension Parameters ---------- T : float Temperature, [K] Returns ------- σ : float Surface tension, [N/m] Notes ----- Raise :class:`NotImplementedError` if input isn't in limit: * 248.15 ≤ T ≤ 647 * Estrapolate to -25ºC in supercooled liquid metastable state Examples -------- >>> _Tension(300) 0.0716859625 >>> _Tension(450) 0.0428914992 References ---------- IAPWS, Revised Release on Surface Tension of Ordinary Water Substance June 2014, http://www.iapws.org/relguide/Surf-H2O.html """ if 248.15 <= T <= Tc: tau = 1-T/Tc sigma = 235.8 * tau**1.256 * (1-0.625*tau) # The equation give surface tension in mN/m², converted to N/m² return 1e-3*sigma raise NotImplementedError("Incoming out of bound")
[docs]def _Dielectric(rho, T): """Equation for the Dielectric constant Parameters ---------- rho : float Density, [kg/m³] T : float Temperature, [K] Returns ------- epsilon : float Dielectric constant, [-] Notes ----- Raise :class:`NotImplementedError` if input isn't in limit: * 238 ≤ T ≤ 1200 Examples -------- >>> _Dielectric(999.242866, 298.15) 78.5907250 >>> _Dielectric(26.0569558, 873.15) 1.12620970 References ---------- IAPWS, Release on the Static Dielectric Constant of Ordinary Water Substance for Temperatures from 238 K to 873 K and Pressures up to 1000 MPa, http://www.iapws.org/relguide/Dielec.html """ # Check input parameters if T < 238 or T > 1200: raise NotImplementedError("Incoming out of bound") k = 1.380658e-23 Na = 6.0221367e23 alfa = 1.636e-40 epsilon0 = 8.854187817e-12 mu = 6.138e-30 d = rho/rhoc Tr = Tc/T li = [1, 1, 1, 2, 3, 3, 4, 5, 6, 7, 10] lj = [0.25, 1, 2.5, 1.5, 1.5, 2.5, 2, 2, 5, 0.5, 10] ni = [0.978224486826, -0.957771379375, 0.237511794148, 0.714692244396, -0.298217036956, -0.108863472196, 0.949327488264e-1, -.980469816509e-2, 0.165167634970e-4, 0.937359795772e-4, -0.12317921872e-9, 0.196096504426e-2] g = 1+ni[11]*d/(Tc/228/Tr-1)**1.2 for n, i, j in zip(ni, li, lj): g += n * d**i * Tr**j A = Na*mu**2*rho*g/M*1000/epsilon0/k/T B = Na*alfa*rho/3/M*1000/epsilon0 e = (1+A+5*B+(9+2*A+18*B+A**2+10*A*B+9*B**2)**0.5)/4/(1-B) return e
[docs]def _Refractive(rho, T, lr=0.5893): """Equation for the refractive index Parameters ---------- rho : float Density, [kg/m³] T : float Temperature, [K] lr : float, optional Light Wavelength, [μm] Returns ------- n : float Refractive index, [-] Notes ----- Raise :class:`NotImplementedError` if input isn't in limit: * 0 ≤ ρ ≤ 1060 * 261.15 ≤ T ≤ 773.15 * 0.2 ≤ λ ≤ 1.1 Examples -------- >>> _Refractive(997.047435, 298.15, 0.2265) 1.39277824 >>> _Refractive(30.4758534, 773.15, 0.5893) 1.00949307 References ---------- IAPWS, Release on the Refractive Index of Ordinary Water Substance as a Function of Wavelength, Temperature and Pressure, http://www.iapws.org/relguide/rindex.pdf """ # Check input parameters if rho < 0 or rho > 1060 or \ T < 261.15 or T > 773.15 or \ lr < 0.2 or lr > 1.1: raise NotImplementedError("Incoming out of bound") Lir = 5.432937 Luv = 0.229202 d = rho/1000. Tr = T/273.15 L = lr/0.589 a = [0.244257733, 0.974634476e-2, -0.373234996e-2, 0.268678472e-3, 0.158920570e-2, 0.245934259e-2, 0.900704920, -0.166626219e-1] A = d*(a[0]+a[1]*d+a[2]*Tr+a[3]*L**2*Tr+a[4]/L**2+a[5]/(L**2-Luv**2)+a[6]/( L**2-Lir**2)+a[7]*d**2) return ((2*A+1)/(1-A))**0.5
[docs]def _Kw(rho, T): """Equation for the ionization constant of ordinary water Parameters ---------- rho : float Density, [kg/m³] T : float Temperature, [K] Returns ------- pKw : float Ionization constant in -log10(kw), [-] Notes ----- Raise :class:`NotImplementedError` if input isn't in limit: * 0 ≤ ρ ≤ 1250 * 273.15 ≤ T ≤ 1073.15 Examples -------- >>> _Kw(1000, 300) 13.906565 References ---------- IAPWS, Release on the Ionization Constant of H2O, http://www.iapws.org/relguide/Ionization.pdf """ # Check input parameters if rho < 0 or rho > 1250 or T < 273.15 or T > 1073.15: raise NotImplementedError("Incoming out of bound") # The internal method of calculation use rho in g/cm³ d = rho/1000. # Water molecular weight different Mw = 18.015268 gamma = [6.1415e-1, 4.825133e4, -6.770793e4, 1.01021e7] pKg = 0 for i, g in enumerate(gamma): pKg += g/T**i Q = d*exp(-0.864671+8659.19/T-22786.2/T**2*d**(2./3)) pKw = -12*(log10(1+Q)-Q/(Q+1)*d*(0.642044-56.8534/T-0.375754*d)) + \ pKg+2*log10(Mw/1000) return pKw
[docs]def _Conductivity(rho, T): """Equation for the electrolytic conductivity of liquid and dense supercrítical water Parameters ---------- rho : float Density, [kg/m³] T : float Temperature, [K] Returns ------- K : float Electrolytic conductivity, [S/m] Notes ----- Raise :class:`NotImplementedError` if input isn't in limit: * 600 ≤ ρ ≤ 1200 * 273.15 ≤ T ≤ 1073.15 Examples -------- >>> _Conductivity(1000, 373.15) 1.13 References ---------- IAPWS, Electrolytic Conductivity (Specific Conductance) of Liquid and Dense Supercritical Water from 0°C to 800°C and Pressures up to 1000 MPa, http://www.iapws.org/relguide/conduct.pdf """ # density in g/l rho_ = rho/1000 # This guideline predates the current standard on the ionization constant, # therefore the standard accepted at that time must be used in order to # obtain the values of the tables for testing. # Marshall, W.L., Franck, E.U. # Ion product of water substance, 0-1000ºC, 1-10,000 bars New International # Formulation and its background # J. Phys. Chem. Ref. Data 10(2) (1981) 295-304 # doi: 10.1063/1.555643 # Eq 4 kw = 10**(-4.098 - 3245.2/T + 2.2362e5/T**2 - 3.984e7/T**3 + (13.957 - 1262.3/T + 8.5641e5/T**2)*log10(rho_)) # kw = 10**-_Kw(rho, T) A = [1850., 1410., 2.16417e-6, 1.81609e-7, -1.75297e-9, 7.20708e-12] B = [16., 11.6, 3.26e-4, -2.3e-6, 1.1e-8] t = T-273.15 Loo = A[0]-1/(1/A[1] + A[2]*t + A[3]*t**2 + A[4]*t**3 + A[5]*t**4) # Eq 5 rho_h = B[0]-1/(1/B[1] + B[2]*t + B[3]*t**2 + B[4]*t**3) # Eq 6 # Eq 4 L_o = (rho_h-rho_)*Loo/rho_h # Eq 1 k = 1e-3*L_o*kw**0.5*rho_ return k*1e2
# Heavy water transport properties
[docs]def _D2O_Viscosity(rho, T, fase=None, drho=None): """Equation for the Viscosity of heavy water Parameters ---------- rho : float Density, [kg/m³] T : float Temperature, [K] fase: dict, optional for calculate critical enhancement phase properties drho: float, optional for calculate critical enhancement [∂ρ/∂P]T at reference state, Returns ------- μ : float Viscosity, [Pa·s] Examples -------- >>> _D2O_Viscosity(998, 298.15) 0.0008897351001498108 >>> _D2O_Viscosity(600, 873.15) 7.743019522728247e-05 References ---------- IAPWS, Release on the IAPWS Formulation 2020 for the Viscosity of Heavy Water, http://iapws.org/relguide/D2Ovisc.pdf """ Tr = T/Tc_D2O rhor = rho/356 # Eq 11, viscosity in the dilute-gas limit no = 0.889754+61.22217*Tr-44.8866*Tr**2+111.5812*Tr**3+3.547412*Tr**4 do = 0.79637+2.38127*Tr-0.33463*Tr**2+2.669*Tr**3+0.000211366*Tr**4 mu0 = Tr**0.5 * no/do # Eq 12 hi = [0, 2, 3, 4, 5, 6, 0, 1, 3, 4, 6, 0, 1, 5, 0, 1, 2, 5, 6, 0, 2, 3, 5, 2, 2] hj = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 6] Hij = [0.510953, -0.558947, -2.718820, 0.480990, 2.404510, -1.824320, 0.275847, 0.762957, 1.760340, 0.0819086, 1.417750, -0.228148, -0.321497, -2.302500, 0.0661035, 0.0449393, 1.466670, 0.938984, -0.108354, -0.00481265, -1.545710, -0.0570938, -0.0753783, 0.553080, -0.0650201] arr = [(1/Tr-1)**i*(rhor-1)**j*hij for i, j, hij in zip(hi, hj, Hij)] mu1 = exp(rhor*sum(arr)) # Critical enhancement if fase and drho: qc = 1/1.9 qd = 1/0.4 # Eq 21 DeltaX = Pc_D2O*rhor**2*(fase.drhodP_T/rho-drho/rho*1.5/Tr) if DeltaX < 0: DeltaX = 0 # Eq 20 X = 0.13*(DeltaX/0.06)**(0.63/1.239) if X <= 0.03021806692: # Eq 15 Y = qc/5*X*(qd*X)**5*(1-qc*X+(qc*X)**2-765/504*(qd*X)**2) else: Fid = acos((1+qd**2*X**2)**-0.5) # Eq 17 w = abs((qc*X-1)/(qc*X+1))**0.5*tan(Fid/2) # Eq 19 # Eq 18 if qc*X > 1: Lw = log((1+w)/(1-w)) else: Lw = 2*atan(abs(w)) # Eq 16 Y = sin(3*Fid)/12-sin(2*Fid)/4/qc/X+(1-5/4*(qc*X)**2)/( qc*X)**2*sin(Fid)-((1-3/2*(qc*X)**2)*Fid-abs(( qc*X)**2-1)**1.5*Lw)/(qc*X)**3 # Eq 14 mu2 = exp(0.068*Y) else: mu2 = 1 return mu0*mu1*mu2*1e-6
[docs]def _D2O_ThCond(rho, T, fase=None, drho=None): """Equation for the thermal conductivity of heavy water Parameters ---------- rho : float Density, [kg/m³] T : float Temperature, [K] fase: dict, optional for calculate critical enhancement phase properties drho: float, optional for calculate critical enhancement [∂ρ/∂P]T at reference state, Returns ------- k : float Thermal conductivity, [W/mK] Examples -------- >>> _D2O_ThCond(998, 298.15) 0.6077128675880629 >>> _D2O_ThCond(0, 873.15) 0.07910346589648833 References ---------- IAPWS, Release on the IAPWS Formulation 2021 for the Thermal Conductivity of Heavy Water, http://iapws.org/relguide/D2OThCond.pdf """ d = rho/356 Tr = T/643.847 # Eq 16 no = [1, 3.3620798, -1.0191198, 2.8518117] do = [0.10779213, -0.034637234, 0.036603464, 0.0091018912] k0 = Tr**0.5 * sum(Li*Tr**i for i, Li in enumerate(no)) / \ sum(Li*Tr**i for i, Li in enumerate(do)) # Eq 17 li = [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4] lj = [0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5] nij = [1.50933576, -0.65831078, 0.111174263, 0.140185152, -0.0656227722, 0.00785155213, 2.8414715, -2.9826577, 1.34357932, -0.599233641, 0.28116337, -0.0533292833, 4.86095723, -6.19784468, 2.20941867, 0.224691518, -0.322191265, 0.0596204654, 2.06156007, -3.48612456, 1.47962309, 0.625101458, -0.56123225, 0.0974446139, -2.06105687, 0.416240028, 2.92524513, -2.81703583, 1.00551476, -0.127884416] k1 = exp(d*sum((1/Tr-1)**i * n*(d-1)**j for i, j, n in zip(li, lj, nij))) # Critical enhancement if fase and drho: Rg = 0.415151994 DeltaX = d*(Pc_D2O/rhoc_D2O*fase.drhodP_T-Pc_D2O/rhoc_D2O*drho*1.5/Tr) if DeltaX < 0: DeltaX = 0 X = 0.13*(DeltaX/0.06)**(0.63/1.239) # Eq 22 y = X/0.36 # Eq 20 # Eq 19 if y < 1.2e-7: Z = 0 else: Z = 2/pi/y*(((1-1/fase.cp_cv)*atan(y)+y/fase.cp_cv)-( 1-exp(-1/(1/y+y**2/3/d**2)))) # Eq 18 k2 = 175.987*d*fase.cp/Rg*Tr/fase.mu*1e-6*Z else: # No critical enhancement k2 = 0 # Eq 15 k = k0*k1+k2 return 1e-3*k
[docs]def _D2O_Tension(T): """Equation for the surface tension of heavy water Parameters ---------- T : float Temperature, [K] Returns ------- σ : float Surface tension, [N/m] Notes ----- Raise :class:`NotImplementedError` if input isn't in limit: * 269.65 ≤ T ≤ 643.847 Examples -------- >>> _D2O_Tension(298.15) 0.07186 >>> _D2O_Tension(573.15) 0.01399 References ---------- IAPWS, Release on Surface Tension of Heavy Water Substance, http://www.iapws.org/relguide/surfd2o.pdf """ Tr = T/643.847 if 269.65 <= T < 643.847: return 1e-3*(238*(1-Tr)**1.25*(1-0.639*(1-Tr))) raise NotImplementedError("Incoming out of bound")
[docs]def _D2O_Sublimation_Pressure(T): """Sublimation Pressure correlation for heavy water Parameters ---------- T : float Temperature, [K] Returns ------- P : float Pressure at sublimation line, [MPa] Notes ----- Raise :class:`NotImplementedError` if input isn't in limit: * 210 ≤ T ≤ 276.969 Examples -------- >>> _Sublimation_Pressure(245) 3.27390934e-5 References ---------- IAPWS, Revised Release on the IAPWS Formulation 2017 for the Thermodynamic Properties of Heavy Water, http://www.iapws.org/relguide/Heavy.html. """ if 210 <= T <= 276.969: Tita = T/276.969 suma = 0 ai = [-0.1314226e2, 0.3212969e2] ti = [-1.73, -1.42] for a, t in zip(ai, ti): suma += a*(1-Tita**t) return exp(suma)*0.00066159 raise NotImplementedError("Incoming out of bound")
[docs]def _D2O_Melting_Pressure(T, ice="Ih"): """Melting Pressure correlation for heavy water Parameters ---------- T : float Temperature, [K] ice: string Type of ice: Ih, III, V, VI, VII. Below 276.969 is a mandatory input, the ice Ih is the default value. Above 276.969, the ice type is unnecesary. Returns ------- P : float Pressure at melting line, [MPa] Notes ----- Raise :class:`NotImplementedError` if input isn't in limit: * 254.415 ≤ T ≤ 315 Examples -------- >>> _D2O__Melting_Pressure(260) 8.947352740189152e-06 >>> _D2O__Melting_Pressure(254, "III") 268.6846466336108 References ---------- IAPWS, Revised Release on the Pressure along the Melting and Sublimation Curves of Ordinary Water Substance, http://iapws.org/relguide/MeltSub.html. """ if ice == "Ih" and 254.415 <= T <= 276.969: # Ice Ih, Eq 9 Tita = T/276.969 ai = [-0.30153e5, 0.692503e6] ti = [5.5, 8.2] suma = 1 for a, t in zip(ai, ti): suma += a*(1-Tita**t) P = suma*0.00066159 elif ice == "III" and 254.415 < T <= 258.661: # Ice III, Eq 10 Tita = T/254.415 P = 222.41*(1-0.802871*(1-Tita**33)) elif ice == "V" and 258.661 < T <= 275.748: # Ice V, Eq 11 Tita = T/258.661 P = 352.19*(1-1.280388*(1-Tita**7.6)) elif (ice == "VI" and 275.748 < T <= 276.969) or 276.969 < T <= 315: # Ice VI Tita = T/275.748 P = 634.53*(1-1.276026*(1-Tita**4)) else: raise NotImplementedError("Incoming out of bound") return P
[docs]def _Henry(T, gas, liquid="H2O"): """Equation for the calculation of Henry's constant Parameters ---------- T : float Temperature, [K] gas : string Name of gas to calculate solubility liquid : string Name of liquid solvent, can be H20 (default) or D2O Returns ------- kw : float Henry's constant, [MPa] Notes ----- The gas availables for H2O solvent are He, Ne, Ar, Kr, Xe, H2, N2, O2, CO, CO2, H2S, CH4, C2H6, SF6 For D2O as solvent He, Ne, Ar, Kr, Xe, D2, CH4 Raise :class:`NotImplementedError` if input gas or liquid are unsupported Examples -------- >>> _Henry(500, "He") 1.1973 >>> _Henry(300, "D2", "D2O") 1.6594 References ---------- IAPWS, Guideline on the Henry's Constant and Vapor-Liquid Distribution Constant for Gases in H2O and D2O at High Temperatures, http://www.iapws.org/relguide/HenGuide.html """ if liquid == "D2O": gas += "(D2O)" limit = { "He": (273.21, 553.18), "Ne": (273.20, 543.36), "Ar": (273.19, 568.36), "Kr": (273.19, 525.56), "Xe": (273.22, 574.85), "H2": (273.15, 636.09), "N2": (278.12, 636.46), "O2": (274.15, 616.52), "CO": (278.15, 588.67), "CO2": (274.19, 642.66), "H2S": (273.15, 533.09), "CH4": (275.46, 633.11), "C2H6": (275.44, 473.46), "SF6": (283.14, 505.55), "He(D2O)": (288.15, 553.18), "Ne(D2O)": (288.18, 549.96), "Ar(D2O)": (288.30, 583.76), "Kr(D2O)": (288.19, 523.06), "Xe(D2O)": (295.39, 574.85), "D2(D2O)": (288.17, 581.00), "CH4(D2O)": (288.16, 517.46)} # Check input parameters if liquid not in ("D2O", "H2O"): raise NotImplementedError("Solvent liquid unsupported") if gas not in limit: raise NotImplementedError("Gas unsupported") Tmin, Tmax = limit[gas] if T < Tmin or T > Tmax: warnings.warn("Temperature out of data of correlation") if liquid == "D2O": Tc_ = Tc_D2O Pc_ = 21.671 else: Tc_ = Tc Pc_ = Pc Tr = T/Tc_ tau = 1-Tr # Eq 4 if liquid == "H2O": ai = [-7.85951783, 1.84408259, -11.7866497, 22.6807411, -15.9618719, 1.80122502] bi = [1, 1.5, 3, 3.5, 4, 7.5] else: ai = [-7.896657, 24.73308, -27.81128, 9.355913, -9.220083] bi = [1, 1.89, 2, 3, 3.6] ps = Pc_*exp(1/Tr*sum(a*tau**b for a, b in zip(ai, bi))) # Select values from Table 2 par = { "He": (-3.52839, 7.12983, 4.47770), "Ne": (-3.18301, 5.31448, 5.43774), "Ar": (-8.40954, 4.29587, 10.52779), "Kr": (-8.97358, 3.61508, 11.29963), "Xe": (-14.21635, 4.00041, 15.60999), "H2": (-4.73284, 6.08954, 6.06066), "N2": (-9.67578, 4.72162, 11.70585), "O2": (-9.44833, 4.43822, 11.42005), "CO": (-10.52862, 5.13259, 12.01421), "CO2": (-8.55445, 4.01195, 9.52345), "H2S": (-4.51499, 5.23538, 4.42126), "CH4": (-10.44708, 4.66491, 12.12986), "C2H6": (-19.67563, 4.51222, 20.62567), "SF6": (-16.56118, 2.15289, 20.35440), "He(D2O)": (-0.72643, 7.02134, 2.04433), "Ne(D2O)": (-0.91999, 5.65327, 3.17247), "Ar(D2O)": (-7.17725, 4.48177, 9.31509), "Kr(D2O)": (-8.47059, 3.91580, 10.69433), "Xe(D2O)": (-14.46485, 4.42330, 15.60919), "D2(D2O)": (-5.33843, 6.15723, 6.53046), "CH4(D2O)": (-10.01915, 4.73368, 11.75711)} A, B, C = par[gas] # Eq 3 kh = ps*exp(A/Tr+B*tau**0.355/Tr+C*Tr**-0.41*exp(tau)) return kh
[docs]def _Kvalue(T, gas, liquid="H2O"): """Equation for the vapor-liquid distribution constant Parameters ---------- T : float Temperature, [K] gas : string Name of gas to calculate solubility liquid : string Name of liquid solvent, can be H20 (default) or D2O Returns ------- kd : float Vapor-liquid distribution constant, [-] Notes ----- The gas availables for H2O solvent are He, Ne, Ar, Kr, Xe, H2, N2, O2, CO, CO2, H2S, CH4, C2H6, SF6 For D2O as solvent He, Ne, Ar, Kr, Xe, D2, CH4 Raise :class:`NotImplementedError` if input gas or liquid are unsupported Examples -------- >>> _Kvalue(600, "He") 3.8019 >>> _Kvalue(300, "D2", "D2O") 14.3520 References ---------- IAPWS, Guideline on the Henry's Constant and Vapor-Liquid Distribution Constant for Gases in H2O and D2O at High Temperatures, http://www.iapws.org/relguide/HenGuide.html """ if liquid == "D2O": gas += "(D2O)" limit = { "He": (273.21, 553.18), "Ne": (273.20, 543.36), "Ar": (273.19, 568.36), "Kr": (273.19, 525.56), "Xe": (273.22, 574.85), "H2": (273.15, 636.09), "N2": (278.12, 636.46), "O2": (274.15, 616.52), "CO": (278.15, 588.67), "CO2": (274.19, 642.66), "H2S": (273.15, 533.09), "CH4": (275.46, 633.11), "C2H6": (275.44, 473.46), "SF6": (283.14, 505.55), "He(D2O)": (288.15, 553.18), "Ne(D2O)": (288.18, 549.96), "Ar(D2O)": (288.30, 583.76), "Kr(D2O)": (288.19, 523.06), "Xe(D2O)": (295.39, 574.85), "D2(D2O)": (288.17, 581.00), "CH4(D2O)": (288.16, 517.46)} # Check input parameters if liquid not in ("D2O", "H2O"): raise NotImplementedError("Solvent liquid unsupported") if gas not in limit: raise NotImplementedError("Gas unsupported") Tmin, Tmax = limit[gas] if T < Tmin or T > Tmax: warnings.warn("Temperature out of data of correlation") if liquid == "D2O": Tc_ = Tc_D2O else: Tc_ = Tc Tr = T/Tc_ tau = 1-Tr # Eq 6 if liquid == "H2O": ci = [1.99274064, 1.09965342, -0.510839303, -1.75493479, -45.5170352, -6.7469445e5] di = [1/3, 2/3, 5/3, 16/3, 43/3, 110/3] q = -0.023767 else: ci = [2.7072, 0.58662, -1.3069, -45.663] di = [0.374, 1.45, 2.6, 12.3] q = -0.024552 f = sum(c*tau**d for c, d in zip(ci, di)) # Select values from Table 2 par = {"He": (2267.4082, -2.9616, -3.2604, 7.8819), "Ne": (2507.3022, -38.6955, 110.3992, -71.9096), "Ar": (2310.5463, -46.7034, 160.4066, -118.3043), "Kr": (2276.9722, -61.1494, 214.0117, -159.0407), "Xe": (2022.8375, 16.7913, -61.2401, 41.9236), "H2": (2286.4159, 11.3397, -70.7279, 63.0631), "N2": (2388.8777, -14.9593, 42.0179, -29.4396), "O2": (2305.0674, -11.3240, 25.3224, -15.6449), "CO": (2346.2291, -57.6317, 204.5324, -152.6377), "CO2": (1672.9376, 28.1751, -112.4619, 85.3807), "H2S": (1319.1205, 14.1571, -46.8361, 33.2266), "CH4": (2215.6977, -0.1089, -6.6240, 4.6789), "C2H6": (2143.8121, 6.8859, -12.6084, 0), "SF6": (2871.7265, -66.7556, 229.7191, -172.7400), "He(D2O)": (2293.2474, -54.7707, 194.2924, -142.1257), "Ne(D2O)": (2439.6677, -93.4934, 330.7783, -243.0100), "Ar(D2O)": (2269.2352, -53.6321, 191.8421, -143.7659), "Kr(D2O)": (2250.3857, -42.0835, 140.7656, -102.7592), "Xe(D2O)": (2038.3656, 68.1228, -271.3390, 207.7984), "D2(D2O)": (2141.3214, -1.9696, 1.6136, 0), "CH4(D2O)": (2216.0181, -40.7666, 152.5778, -117.7430)} E, F, G, H = par[gas] # Eq 5 kd = exp(q*F+E/T*f+(F+G*tau**(2./3)+H*tau)*exp((273.15-T)/100)) return kd