#!/usr/bin/python
# -*- coding: utf-8 -*-
# pylint: disable=invalid-name
# pylint: disable=too-many-lines, too-many-locals, too-many-instance-attributes
# pylint: disable=too-many-branches, too-many-statements
"""
Implemented multiparameter equation of state as a Helmholtz free energy:
* :class:`MEoS`: Base class of multiparameter equation of state
* :class:`IAPWS95`: 2016 revision of 1995 formulation for ordinaty water
* :class:`D2O`: 2017 formulation for heavy water.
"""
from __future__ import division
try:
from concurrent.futures import ProcessPoolExecutor
except ImportError:
# python2
pass
from itertools import repeat
import json
import os
import platform
import warnings
from numpy import exp, log, ndarray
from numpy.polynomial import chebyshev
from scipy.optimize import fsolve
from .iapws97 import _TSat_P, IAPWS97
from ._iapws import M, Tc, Pc, rhoc, Tc_D2O, Pc_D2O, rhoc_D2O
from ._iapws import _Viscosity, _ThCond, _Dielectric, _Refractive, _Tension
from ._iapws import _D2O_Viscosity, _D2O_ThCond, _D2O_Tension
from ._utils import _fase, getphase, deriv_H, _instanceBuilder
[docs]
def _phir(tau, delta, coef):
"""Residual contribution to the adimensional free Helmholtz energy
Parameters
----------
tau : float
Inverse reduced temperature Tc/T, [-]
delta : float
Reduced density rho/rhoc, [-]
coef : dict
Dictionary with equation of state parameters
Returns
-------
fir : float
Adimensional free Helmholtz energy
References
----------
IAPWS, Revised Release on the IAPWS Formulation 1995 for the
Thermodynamic Properties of Ordinary Water Substance for General and
Scientific Use, September 2016, Table 5
http://www.iapws.org/relguide/IAPWS-95.html
"""
fir = 0
# Polinomial terms
nr1 = coef.get("nr1", [])
d1 = coef.get("d1", [])
t1 = coef.get("t1", [])
for n, d, t in zip(nr1, d1, t1):
fir += n*delta**d*tau**t
# Exponential terms
nr2 = coef.get("nr2", [])
d2 = coef.get("d2", [])
g2 = coef.get("gamma2", [])
t2 = coef.get("t2", [])
c2 = coef.get("c2", [])
for n, d, g, t, c in zip(nr2, d2, g2, t2, c2):
fir += n*delta**d*tau**t*exp(-g*delta**c)
# Gaussian terms
nr3 = coef.get("nr3", [])
d3 = coef.get("d3", [])
t3 = coef.get("t3", [])
a3 = coef.get("alfa3", [])
e3 = coef.get("epsilon3", [])
b3 = coef.get("beta3", [])
g3 = coef.get("gamma3", [])
for n, d, t, a, e, b, g in zip(nr3, d3, t3, a3, e3, b3, g3):
fir += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)
# Non analitic terms
nr4 = coef.get("nr4", [])
a4 = coef.get("a4", [])
b4 = coef.get("b4", [])
Ai = coef.get("A", [])
Bi = coef.get("B", [])
Ci = coef.get("C", [])
Di = coef.get("D", [])
bt4 = coef.get("beta4", [])
for n, a, b, A, B, C, D, bt in zip(nr4, a4, b4, Ai, Bi, Ci, Di, bt4):
Tita = (1-tau)+A*((delta-1)**2)**(0.5/bt)
F = exp(-C*(delta-1)**2-D*(tau-1)**2)
Delta = Tita**2+B*((delta-1)**2)**a
fir += n*Delta**b*delta*F
return fir
[docs]
def _phird(tau, delta, coef):
r"""Residual contribution to the adimensional free Helmholtz energy, delta
derivative
Parameters
----------
tau : float
Inverse reduced temperature Tc/T, [-]
delta : float
Reduced density rho/rhoc, [-]
coef : dict
Dictionary with equation of state parameters
Returns
-------
fird : float
.. math::
\left.\frac{\partial \phi^r_{\delta}}{\partial \delta}\right|_{\tau}
References
----------
IAPWS, Revised Release on the IAPWS Formulation 1995 for the
Thermodynamic Properties of Ordinary Water Substance for General and
Scientific Use, September 2016, Table 5
http://www.iapws.org/relguide/IAPWS-95.html
"""
fird = 0
# Polinomial terms
nr1 = coef.get("nr1", [])
d1 = coef.get("d1", [])
t1 = coef.get("t1", [])
for n, d, t in zip(nr1, d1, t1):
fird += n*d*delta**(d-1)*tau**t
# Exponential terms
nr2 = coef.get("nr2", [])
d2 = coef.get("d2", [])
g2 = coef.get("gamma2", [])
t2 = coef.get("t2", [])
c2 = coef.get("c2", [])
for n, d, g, t, c in zip(nr2, d2, g2, t2, c2):
fird += n*exp(-g*delta**c)*delta**(d-1)*tau**t*(d-g*c*delta**c)
# Gaussian terms
nr3 = coef.get("nr3", [])
d3 = coef.get("d3", [])
t3 = coef.get("t3", [])
a3 = coef.get("alfa3", [])
e3 = coef.get("epsilon3", [])
b3 = coef.get("beta3", [])
g3 = coef.get("gamma3", [])
for n, d, t, a, e, b, g in zip(nr3, d3, t3, a3, e3, b3, g3):
fird += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
d/delta-2*a*(delta-e))
# Non analitic terms
nr4 = coef.get("nr4", [])
a4 = coef.get("a4", [])
b4 = coef.get("b4", [])
Ai = coef.get("A", [])
Bi = coef.get("B", [])
Ci = coef.get("C", [])
Di = coef.get("D", [])
bt4 = coef.get("beta4", [])
for n, a, b, A, B, C, D, bt in zip(nr4, a4, b4, Ai, Bi, Ci, Di, bt4):
Tita = (1-tau)+A*((delta-1)**2)**(0.5/bt)
F = exp(-C*(delta-1)**2-D*(tau-1)**2)
Fd = -2*C*F*(delta-1)
Delta = Tita**2+B*((delta-1)**2)**a
Deltad = (delta-1)*(A*Tita*2/bt*((delta-1)**2)**(0.5/bt-1)
+ 2*B*a*((delta-1)**2)**(a-1))
if Delta == 0:
DeltaBd = 0
else:
DeltaBd = b*Delta**(b-1)*Deltad
fird += n*(Delta**b*(F+delta*Fd)+DeltaBd*delta*F)
return fird
[docs]
def _phirt(tau, delta, coef):
r"""Residual contribution to the adimensional free Helmholtz energy, tau
derivative
Parameters
----------
tau : float
Inverse reduced temperature Tc/T, [-]
delta : float
Reduced density rho/rhoc, [-]
coef : dict
Dictionary with equation of state parameters
Returns
-------
firt : float
.. math::
\left.\frac{\partial \phi^r_{\tau}}{\partial \tau}\right|_{\delta}
References
----------
IAPWS, Revised Release on the IAPWS Formulation 1995 for the
Thermodynamic Properties of Ordinary Water Substance for General and
Scientific Use, September 2016, Table 5
http://www.iapws.org/relguide/IAPWS-95.html
"""
firt = 0
# Polinomial terms
nr1 = coef.get("nr1", [])
d1 = coef.get("d1", [])
t1 = coef.get("t1", [])
for n, d, t in zip(nr1, d1, t1):
firt += n*t*delta**d*tau**(t-1)
# Exponential terms
nr2 = coef.get("nr2", [])
d2 = coef.get("d2", [])
g2 = coef.get("gamma2", [])
t2 = coef.get("t2", [])
c2 = coef.get("c2", [])
for n, d, g, t, c in zip(nr2, d2, g2, t2, c2):
firt += n*t*delta**d*tau**(t-1)*exp(-g*delta**c)
# Gaussian terms
nr3 = coef.get("nr3", [])
d3 = coef.get("d3", [])
t3 = coef.get("t3", [])
a3 = coef.get("alfa3", [])
e3 = coef.get("epsilon3", [])
b3 = coef.get("beta3", [])
g3 = coef.get("gamma3", [])
for n, d, t, a, e, b, g in zip(nr3, d3, t3, a3, e3, b3, g3):
firt += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
t/tau-2*b*(tau-g))
# Non analitic terms
nr4 = coef.get("nr4", [])
a4 = coef.get("a4", [])
b4 = coef.get("b4", [])
Ai = coef.get("A", [])
Bi = coef.get("B", [])
Ci = coef.get("C", [])
Di = coef.get("D", [])
bt4 = coef.get("beta4", [])
for n, a, b, A, B, C, D, bt in zip(nr4, a4, b4, Ai, Bi, Ci, Di, bt4):
Tita = (1-tau)+A*((delta-1)**2)**(0.5/bt)
F = exp(-C*(delta-1)**2-D*(tau-1)**2)
Ft = -2*D*F*(tau-1)
Delta = Tita**2+B*((delta-1)**2)**a
if Delta == 0:
DeltaBt = 0
else:
DeltaBt = -2*Tita*b*Delta**(b-1)
firt += n*delta*(DeltaBt*F+Delta**b*Ft)
return firt
[docs]
class MEoS(_fase):
r"""
General implementation of multiparameter equation of state. From this
derived all child class specified per individual compounds
Parameters
----------
T : float
Temperature, [K]
P : float
Pressure, [MPa]
rho : float
Density, [kg/m³]
v : float
Specific volume, [m³/kg]
h : float
Specific enthalpy, [kJ/kg]
s : float
Specific entropy, [kJ/kgK]
u : float
Specific internal energy, [kJ/kg]
x : float
Vapor quality, [-]
l : float, optional
Wavelength of light, for refractive index, [μm]
rho0 : float, optional
Initial value of density, to improve iteration, [kg/m³]
T0 : float, optional
Initial value of temperature, to improve iteration, [K]
x0 : Initial value of vapor quality, necessary in bad input pair definition
where there are two valid solution (T-h, T-s)
Notes
-----
* It needs two incoming properties of T, P, rho, h, s, u.
* v as a alternate input parameter to rho
* T-x, P-x, preferred input pair to specified a point in two phases region
The calculated instance has the following properties:
* P: Pressure, [MPa]
* T: Temperature, [K]
* x: Vapor quality, [-]
* g: Specific Gibbs free energy, [kJ/kg]
* a: Specific Helmholtz free energy, [kJ/kg]
* v: Specific volume, [m³/kg]
* r: Density, [kg/m³]
* h: Specific enthalpy, [kJ/kg]
* u: Specific internal energy, [kJ/kg]
* s: Specific entropy, [kJ/kg·K]
* cp: Specific isobaric heat capacity, [kJ/kg·K]
* cv: Specific isochoric heat capacity, [kJ/kg·K]
* cp_cv: Heat capacity ratio, [-]
* Z: Compression factor, [-]
* fi: Fugacity coefficient, [-]
* f: Fugacity, [MPa]
* gamma: Isoentropic exponent, [-]
* alfav: Isobaric cubic expansion coefficient, [1/K]
* kappa: Isothermal compressibility, [1/MPa]
* kappas: Adiabatic compresibility, [1/MPa]
* alfap: Relative pressure coefficient, [1/K]
* betap: Isothermal stress coefficient, [kg/m³]
* joule: Joule-Thomson coefficient, [K/MPa]
* betas: Isoentropic temperature-pressure coefficient, [-]
* Gruneisen: Gruneisen parameter, [-]
* virialB: Second virial coefficient, [m³/kg]
* virialC: Third virial coefficient, [m⁶/kg²]
* dpdT_rho: Derivatives, dp/dT at constant rho, [MPa/K]
* dpdrho_T: Derivatives, dp/drho at constant T, [MPa·m³/kg]
* drhodT_P: Derivatives, drho/dT at constant P, [kg/m³·K]
* drhodP_T: Derivatives, drho/dP at constant T, [kg/m³·MPa]
* dhdT_rho: Derivatives, dh/dT at constant rho, [kJ/kg·K]
* dhdP_T: Isothermal throttling coefficient, [kJ/kg·MPa]
* dhdT_P: Derivatives, dh/dT at constant P, [kJ/kg·K]
* dhdrho_T: Derivatives, dh/drho at constant T, [kJ·m³/kg²]
* dhdrho_P: Derivatives, dh/drho at constant P, [kJ·m³/kg²]
* dhdP_rho: Derivatives, dh/dP at constant rho, [kJ/kg·MPa]
* kt: Isothermal Expansion Coefficient, [-]
* ks: Adiabatic Compressibility, [1/MPa]
* Ks: Adiabatic bulk modulus, [MPa]
* Kt: Isothermal bulk modulus, [MPa]
* v0: Ideal specific volume, [m³/kg]
* rho0: Ideal gas density, [kg/m³]
* u0: Ideal specific internal energy, [kJ/kg]
* h0: Ideal specific enthalpy, [kJ/kg]
* s0: Ideal specific entropy, [kJ/kg·K]
* a0: Ideal specific Helmholtz free energy, [kJ/kg]
* g0: Ideal specific Gibbs free energy, [kJ/kg]
* cp0: Ideal specific isobaric heat capacity, [kJ/kg·K]
* cv0: Ideal specific isochoric heat capacity, [kJ/kg·K]
* w0: Ideal speed of sound, [m/s]
* gamma0: Ideal isoentropic exponent, [-]
* w: Speed of sound, [m/s]
* mu: Dynamic viscosity, [Pa·s]
* nu: Kinematic viscosity, [m²/s]
* k: Thermal conductivity, [W/m·K]
* alfa: Thermal diffusivity, [m²/s]
* sigma: Surface tension, [N/m]
* epsilon: Dielectric constant, [-]
* n: Refractive index, [-]
* Prandt: Prandtl number, [-]
* Pr: Reduced Pressure, [-]
* Tr: Reduced Temperature, [-]
* Hvap: Vaporization heat, [kJ/kg]
* Svap: Vaporization entropy, [kJ/kg·K]
* Z_rho: :math:`(Z-1)/\rho`, [m³/kg]
* IntP: Internal pressure, [MPa]
* invT: Negative reciprocal temperature, [1/K]
* hInput: Specific heat input, [kJ/kg]
"""
Fi0 = {}
_constants = {}
_Pv = {}
_rhoL = {}
_rhoG = {}
_surf = {}
kwargs = {"T": 0.0,
"P": 0.0,
"rho": 0.0,
"v": 0.0,
"h": None,
"s": None,
"u": None,
"x": None,
"l": 0.5893,
"rho0": None,
"T0": None,
"x0": 0.5}
name = None
M = None
Tc = None
Pc = None
rhoc = None
Tt = None
status = 0
msg = "Undefined"
_mode = None
Liquid = None
Gas = None
Vapor = None
T = None
Tr = None
P = None
Pr = None
x = None
phase = None
sigma = None
virialB = None
virialC = None
Hvap = None
Svap = None
invT = None
v0 = None
rho0 = None
h0 = None
u0 = None
s0 = None
a0 = None
g0 = None
cp0 = None
cv0 = None
cp0_cv = None
gamma0 = None
[docs]
@classmethod
def from_list(cls, p1name, p1val, p2name, p2val):
"""Speed up method using multiprocessing for multiple point calculation
with a fixed input and changing other input parameter
Parameters
----------
p1name : str
string with name of fixed input parameter
p1val : float
fixed input parameter value
p2name : str
string with name of changing input parameter
p2val : list
iterable with values of changing input parameter
Returns
-------
states : list
list with calculated states
"""
py_version = platform.python_version_tuple()[0]
lst = []
if py_version == "3":
with ProcessPoolExecutor() as executor:
states = executor.map(
_instanceBuilder, *(p2val, repeat(cls), repeat(p1name),
repeat(p1val), repeat(p2name)))
for state in states:
lst.append(state)
elif py_version == "2":
# Disable multithreading in python2
def fi(x):
return cls(**{p1name: p1val, p2name: x})
for xi in p2val:
lst.append(fi(xi))
return lst
def __init__(self, **kwargs):
"""Constructor, define common constant and initinialice kwargs"""
self.R = self._constants["R"]/self._constants.get("M", self.M)
self.Zc = self.Pc/self.rhoc/self.R/self.Tc
self.kwargs = MEoS.kwargs.copy()
self.__call__(**kwargs)
def __call__(self, **kwargs):
"""Make instance callable to can add input parameter one to one"""
# Alternative rho input
if "rhom" in kwargs:
kwargs["rho"] = kwargs["rhom"]*self.M
del kwargs["rhom"]
elif kwargs.get("v", 0):
kwargs["rho"] = 1./kwargs["v"]
del kwargs["v"]
elif kwargs.get("vm", 0):
kwargs["rho"] = self.M/kwargs["vm"]
del kwargs["vm"]
self.kwargs.update(kwargs)
if self.calculable:
try:
self.status = 1
self.calculo()
self.msg = ""
except RuntimeError as err:
self.status = 0
self.msg = err.args[0]
raise err
# Add msg for extrapolation state
if self.name == "water" and 130 <= self.T < 273.15:
self.msg = "Extrapolated state"
self.status = 3
warnings.warn("Using extrapolated values")
elif self.name == "water" and 50 <= self.T < 130:
self.msg = "Extrapolated state using Low-Temperature extension"
self.status = 3
warnings.warn("Using extrapolated values and Low-Temperature"
"extension")
@property
def calculable(self):
"""Check if inputs are enough to define state"""
self._mode = ""
if self.kwargs["T"] and self.kwargs["P"]:
self._mode = "TP"
elif self.kwargs["T"] and self.kwargs["rho"]:
self._mode = "Trho"
elif self.kwargs["T"] and self.kwargs["h"] is not None:
self._mode = "Th"
elif self.kwargs["T"] and self.kwargs["s"] is not None:
self._mode = "Ts"
elif self.kwargs["T"] and self.kwargs["u"] is not None:
self._mode = "Tu"
elif self.kwargs["P"] and self.kwargs["rho"]:
self._mode = "Prho"
elif self.kwargs["P"] and self.kwargs["h"] is not None:
self._mode = "Ph"
elif self.kwargs["P"] and self.kwargs["s"] is not None:
self._mode = "Ps"
elif self.kwargs["P"] and self.kwargs["u"] is not None:
self._mode = "Pu"
elif self.kwargs["rho"] and self.kwargs["h"] is not None:
self._mode = "rhoh"
elif self.kwargs["rho"] and self.kwargs["s"] is not None:
self._mode = "rhos"
elif self.kwargs["rho"] and self.kwargs["u"] is not None:
self._mode = "rhou"
elif self.kwargs["h"] is not None and self.kwargs["s"] is not None:
self._mode = "hs"
elif self.kwargs["h"] is not None and self.kwargs["u"] is not None:
self._mode = "hu"
elif self.kwargs["s"] is not None and self.kwargs["u"] is not None:
self._mode = "su"
elif self.kwargs["T"] and self.kwargs["x"] is not None:
self._mode = "Tx"
elif self.kwargs["P"] and self.kwargs["x"] is not None:
self._mode = "Px"
return bool(self._mode)
[docs]
def calculo(self):
"""Calculate procedure"""
T = self.kwargs["T"]
rho = self.kwargs["rho"]
P = self.kwargs["P"]
s = self.kwargs["s"]
h = self.kwargs["h"]
u = self.kwargs["u"]
x = self.kwargs["x"]
# Initial values
T0 = self.kwargs["T0"]
rho0 = self.kwargs["rho0"]
if T0 or rho0:
To = T0
rhoo = rho0
elif self.name == "air":
To = 300
rhoo = 1e-3
else:
try:
st0 = IAPWS97(**self.kwargs)
except NotImplementedError:
To = 300
rhoo = 900
else:
if st0.status:
To = st0.T
rhoo = st0.rho
else:
To = 300
rhoo = 900
self.R = self._constants["R"]/self._constants.get("M", self.M)
rho_c = self._constants.get("rhoref", self.rhoc)
T_c = self._constants.get("Tref", self.Tc)
propiedades = None
if self._mode not in ("Tx", "Px"):
# Method with iteration necessary to get x
if self._mode == "TP":
try:
if self.name != "water":
raise NotImplementedError
st0 = IAPWS97(**self.kwargs)
rhoo = st0.rho
except NotImplementedError:
if rho0:
rhoo = rho0
elif T < self.Tc and \
self._Vapor_Pressure(T) < P < self.Pc:
rhoo = self._Liquid_Density(T)
elif T < self.Tc and P < self.Pc:
rhoo = self._Vapor_Density(T)
else:
rhoo = self.rhoc*3
def f(rho):
delta = rho/rho_c
tau = T_c/T
fird = _phird(tau, delta, self._constants)
Po = (1+delta*fird)*self.R*T*rho
return Po-P*1000
rho = fsolve(f, rhoo)[0]
# Calculate quality
if T > self.Tc:
x = 1
else:
Ps = self._Vapor_Pressure(T)
if Ps*0.95 < P < Ps*1.05:
rhol, rhov, Ps = self._saturation(T)
Ps *= 1e-3
if Ps > P:
x = 1
else:
x = 0
elif self._mode == "Th":
tau = T_c/T
ideal = self._phi0(tau, 1)
fiot = ideal["fiot"]
def f(rho):
delta = rho/rho_c
fird = _phird(tau, delta, self._constants)
firt = _phirt(tau, delta, self._constants)
ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
return ho-h
if T >= self.Tc:
rhoo = self.rhoc
rho = fsolve(f, rhoo)[0]
else:
x0 = self.kwargs["x0"]
rhov = self._Vapor_Density(T)
rhol = self._Liquid_Density(T)
deltaL = rhol/rho_c
deltaG = rhov/rho_c
firdL = _phird(tau, deltaL, self._constants)
firtL = _phirt(tau, deltaL, self._constants)
firdG = _phird(tau, deltaG, self._constants)
firtG = _phirt(tau, deltaG, self._constants)
hl = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
hv = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)
if x0 not in (0, 1) and hl <= h <= hv:
rhol, rhov, Ps = self._saturation(T)
vapor = self._Helmholtz(rhov, T)
liquido = self._Helmholtz(rhol, T)
hv = vapor["h"]
hl = liquido["h"]
x = (h-hl)/(hv-hl)
rho = 1/(x/rhov+(1-x)/rhol)
P = Ps/1000
else:
if h > hv:
rhoo = rhov
else:
rhoo = rhol
rho = fsolve(f, rhoo)[0]
elif self._mode == "Ts":
tau = T_c/T
def f(rho):
if rho < 0:
rho = 1e-20
delta = rho/rho_c
ideal = self._phi0(tau, delta)
fio = ideal["fio"]
fiot = ideal["fiot"]
fir = _phir(tau, delta, self._constants)
firt = _phirt(tau, delta, self._constants)
so = self.R*(tau*(fiot+firt)-fio-fir)
return so-s
if T >= self.Tc:
rhoo = self.rhoc
rho = fsolve(f, rhoo)[0]
else:
rhov = self._Vapor_Density(T)
rhol = self._Liquid_Density(T)
deltaL = rhol/rho_c
deltaG = rhov/rho_c
idealL = self._phi0(tau, deltaL)
idealG = self._phi0(tau, deltaG)
fioL = idealL["fio"]
fioG = idealG["fio"]
fiot = idealL["fiot"]
firL = _phir(tau, deltaL, self._constants)
firtL = _phirt(tau, deltaL, self._constants)
sl = self.R*(tau*(fiot+firtL)-fioL-firL)
firG = _phir(tau, deltaG, self._constants)
firtG = _phirt(tau, deltaG, self._constants)
sv = self.R*(tau*(fiot+firtG)-fioG-firG)
if sl <= s <= sv:
rhol, rhov, Ps = self._saturation(T)
vapor = self._Helmholtz(rhov, T)
liquido = self._Helmholtz(rhol, T)
sv = vapor["s"]
sl = liquido["s"]
x = (s-sl)/(sv-sl)
rho = 1/(x/rhov+(1-x)/rhol)
P = Ps/1000
else:
if s > sv:
rhoo = rhov
else:
rhoo = rhol
rho = fsolve(f, rhoo)[0]
elif self._mode == "Tu":
tau = T_c/T
ideal = self._phi0(tau, 1)
fiot = ideal["fiot"]
def f(rho):
delta = rho/rho_c
fird = _phird(tau, delta, self._constants)
firt = _phirt(tau, delta, self._constants)
Po = (1+delta*fird)*self.R*T*rho
ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
return ho-Po/rho-u
if T >= self.Tc:
rhoo = self.rhoc
rho = fsolve(f, rhoo)[0]
else:
rhov = self._Vapor_Density(T)
rhol = self._Liquid_Density(T)
deltaL = rhol/rho_c
deltaG = rhov/rho_c
firdL = _phird(tau, deltaL, self._constants)
firtL = _phirt(tau, deltaL, self._constants)
firdG = _phird(tau, deltaG, self._constants)
firtG = _phirt(tau, deltaG, self._constants)
PoL = (1+deltaL*firdL)*self.R*T*rhol
PoG = (1+deltaG*firdG)*self.R*T*rhov
hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)
uv = hoG-PoG/rhov
ul = hoL-PoL/rhol
if ul <= u <= uv:
rhol, rhov, Ps = self._saturation(T)
vapor = self._Helmholtz(rhov, T)
liquido = self._Helmholtz(rhol, T)
uv = vapor["h"]-vapor["P"]/rhov
ul = liquido["h"]-liquido["P"]/rhol
x = (u-ul)/(uv-ul)
rho = 1/(x/rhov-(1-x)/rhol)
P = Ps/1000
else:
if u > uv:
rhoo = rhov
else:
rhoo = rhol
rho = fsolve(f, rhoo)[0]
elif self._mode == "Prho":
delta = rho/rho_c
def f(T):
tau = T_c/T
fird = _phird(tau, delta, self._constants)
Po = (1+delta*fird)*self.R*T*rho
return Po-P*1000
T = fsolve(f, To)[0]
rhol = self._Liquid_Density(T)
rhov = self._Vapor_Density(T)
if T == To or rhov <= rho <= rhol:
def f(parr):
T, rhol, rhog = parr
tau = T_c/T
deltaL = rhol/self.rhoc
deltaG = rhog/self.rhoc
firL = _phir(tau, deltaL, self._constants)
firdL = _phird(tau, deltaL, self._constants)
firG = _phir(tau, deltaG, self._constants)
firdG = _phird(tau, deltaG, self._constants)
Jl = rhol*(1+deltaL*firdL)
Jv = rhog*(1+deltaG*firdG)
K = firL-firG
Ps = self.R*T*rhol*rhog/(rhol-rhog)*(K+log(rhol/rhog))
return (Jl-Jv,
Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
Ps - P*1000)
for to in [To, 300, 400, 500, 600]:
rhoLo = self._Liquid_Density(to)
rhoGo = self._Vapor_Density(to)
sol = fsolve(f, [to, rhoLo, rhoGo], full_output=True)
T, rhoL, rhoG = sol[0]
x = (1./rho-1/rhoL)/(1/rhoG-1/rhoL)
if sol[2] == 1 and 0 <= x <= 1 and \
sum(abs(sol[1]["fvec"])) < 1e-5:
break
if sum(abs(sol[1]["fvec"])) > 1e-5:
raise RuntimeError(sol[3])
liquido = self._Helmholtz(rhoL, T)
vapor = self._Helmholtz(rhoG, T)
P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000
elif self._mode == "Ph":
def funcion(parr):
rho, T = parr
delta = rho/rho_c
tau = T_c/T
ideal = self._phi0(tau, delta)
fiot = ideal["fiot"]
fird = _phird(tau, delta, self._constants)
firt = _phirt(tau, delta, self._constants)
Po = (1+delta*fird)*self.R*T*rho
ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
return Po-P*1000, ho-h
rho, T = fsolve(funcion, [rhoo, To])
rhol = self._Liquid_Density(T)
rhov = self._Vapor_Density(T)
if rho == rhoo or rhov <= rho <= rhol:
def f(parr):
T, rhol, rhog, x = parr
tau = T_c/T
deltaL = rhol/self.rhoc
deltaG = rhog/self.rhoc
ideal = self._phi0(tau, deltaL)
fiot = ideal["fiot"]
firL = _phir(tau, deltaL, self._constants)
firdL = _phird(tau, deltaL, self._constants)
firtL = _phirt(tau, deltaL, self._constants)
hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
firG = _phir(tau, deltaG, self._constants)
firdG = _phird(tau, deltaG, self._constants)
firtG = _phirt(tau, deltaG, self._constants)
hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)
Jl = rhol*(1+deltaL*firdL)
Jv = rhog*(1+deltaG*firdG)
K = firL-firG
Ps = self.R*T*rhol*rhog/(rhol-rhog)*(K+log(rhol/rhog))
return (Jl-Jv,
Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
hoL*(1-x)+hoG*x - h,
Ps - P*1000)
for to in [To, 300, 400, 500, 600]:
rLo = self._Liquid_Density(to)
rGo = self._Vapor_Density(to)
sol = fsolve(f, [to, rLo, rGo, 0.5], full_output=True)
T, rhoL, rhoG, x = sol[0]
if sol[2] == 1 and 0 <= x <= 1 and \
sum(abs(sol[1]["fvec"])) < 1e-5:
break
if sum(abs(sol[1]["fvec"])) > 1e-5:
raise RuntimeError(sol[3])
liquido = self._Helmholtz(rhoL, T)
vapor = self._Helmholtz(rhoG, T)
P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000
elif self._mode == "Ps":
try:
x0 = st0.x
except NameError:
x0 = None
if x0 is None or x0 == 0 or x0 == 1:
def f(parr):
rho, T = parr
delta = rho/rho_c
tau = T_c/T
ideal = self._phi0(tau, delta)
fio = ideal["fio"]
fiot = ideal["fiot"]
fird = _phird(tau, delta, self._constants)
fir = _phir(tau, delta, self._constants)
firt = _phirt(tau, delta, self._constants)
Po = (1+delta*fird)*self.R*T*rho
so = self.R*(tau*(fiot+firt)-fio-fir)
return Po-P*1000, so-s
rho, T = fsolve(f, [rhoo, To])
else:
def funcion(parr):
rho, T = parr
rhol, rhov, Ps = self._saturation(T)
vapor = self._Helmholtz(rhov, T)
liquido = self._Helmholtz(rhol, T)
x = (1./rho-1/rhol)/(1/rhov-1/rhol)
return Ps-P*1000, vapor["s"]*x+liquido["s"]*(1-x)-s
rho, T = fsolve(funcion, [2., 500.])
rhol, rhov, Ps = self._saturation(T)
vapor = self._Helmholtz(rhov, T)
liquido = self._Helmholtz(rhol, T)
sv = vapor["s"]
sl = liquido["s"]
x = (s-sl)/(sv-sl)
elif self._mode == "Pu":
def f(parr):
rho, T = parr
delta = rho/rho_c
tau = T_c/T
ideal = self._phi0(tau, delta)
fiot = ideal["fiot"]
fird = _phird(tau, delta, self._constants)
firt = _phirt(tau, delta, self._constants)
Po = (1+delta*fird)*self.R*T*rho
ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
return ho-Po/rho-u, Po-P*1000
sol = fsolve(f, [rhoo, To], full_output=True)
rho, T = sol[0]
rhol = self._Liquid_Density(T)
rhov = self._Vapor_Density(T)
if rho == rhoo or sol[2] != 1:
def f(parr):
T, rhol, rhog, x = parr
tau = T_c/T
deltaL = rhol/self.rhoc
deltaG = rhog/self.rhoc
ideal = self._phi0(tau, deltaL)
fiot = ideal["fiot"]
firL = _phir(tau, deltaL, self._constants)
firdL = _phird(tau, deltaL, self._constants)
firtL = _phirt(tau, deltaL, self._constants)
hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
firG = _phir(tau, deltaG, self._constants)
firdG = _phird(tau, deltaG, self._constants)
firtG = _phirt(tau, deltaG, self._constants)
hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)
Jl = rhol*(1+deltaL*firdL)
Jv = rhog*(1+deltaG*firdG)
K = firL-firG
Ps = self.R*T*rhol*rhog/(rhol-rhog)*(K+log(rhol/rhog))
vu = hoG-Ps/rhog
lu = hoL-Ps/rhol
return (Jl-Jv,
Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
lu*(1-x)+vu*x - u,
Ps - P*1000)
for to in [To, 300, 400, 500, 600]:
rLo = self._Liquid_Density(to)
rGo = self._Vapor_Density(to)
sol = fsolve(f, [to, rLo, rGo, 0.5], full_output=True)
T, rhoL, rhoG, x = sol[0]
if sol[2] == 1 and 0 <= x <= 1 and \
sum(abs(sol[1]["fvec"])) < 1e-5:
break
if sum(abs(sol[1]["fvec"])) > 1e-5:
raise RuntimeError(sol[3])
liquido = self._Helmholtz(rhoL, T)
vapor = self._Helmholtz(rhoG, T)
P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000
elif self._mode == "rhoh":
delta = rho/rho_c
def f(T):
tau = T_c/T
ideal = self._phi0(tau, delta)
fiot = ideal["fiot"]
fird = _phird(tau, delta, self._constants)
firt = _phirt(tau, delta, self._constants)
ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
return ho-h
T = fsolve(f, To)[0]
rhol = self._Liquid_Density(T)
rhov = self._Vapor_Density(T)
if T == To or rhov <= rho <= rhol:
def f(parr):
T, rhol, rhog = parr
tau = T_c/T
deltaL = rhol/self.rhoc
deltaG = rhog/self.rhoc
ideal = self._phi0(tau, deltaL)
fiot = ideal["fiot"]
firL = _phir(tau, deltaL, self._constants)
firdL = _phird(tau, deltaL, self._constants)
firtL = _phirt(tau, deltaL, self._constants)
hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
firG = _phir(tau, deltaG, self._constants)
firdG = _phird(tau, deltaG, self._constants)
firtG = _phirt(tau, deltaG, self._constants)
hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)
Jl = rhol*(1+deltaL*firdL)
Jv = rhog*(1+deltaG*firdG)
K = firL-firG
x = (1./rho-1/rhol)/(1/rhog-1/rhol)
return (Jl-Jv,
Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
hoL*(1-x)+hoG*x - h)
for to in [To, 300, 400, 500, 600]:
rhoLo = self._Liquid_Density(to)
rhoGo = self._Vapor_Density(to)
sol = fsolve(f, [to, rhoLo, rhoGo], full_output=True)
T, rhoL, rhoG = sol[0]
x = (1./rho-1/rhoL)/(1/rhoG-1/rhoL)
if sol[2] == 1 and 0 <= x <= 1 and \
sum(abs(sol[1]["fvec"])) < 1e-5:
break
if sum(abs(sol[1]["fvec"])) > 1e-5:
raise RuntimeError(sol[3])
liquido = self._Helmholtz(rhoL, T)
vapor = self._Helmholtz(rhoG, T)
P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000
elif self._mode == "rhos":
delta = rho/rho_c
def f(T):
tau = T_c/T
ideal = self._phi0(tau, delta)
fio = ideal["fio"]
fiot = ideal["fiot"]
fir = _phir(tau, delta, self._constants)
firt = _phirt(tau, delta, self._constants)
so = self.R*(tau*(fiot+firt)-fio-fir)
return so-s
T = fsolve(f, To)[0]
rhol = self._Liquid_Density(T)
rhov = self._Vapor_Density(T)
if T == To or rhov <= rho <= rhol:
def f(parr):
T, rhol, rhog = parr
tau = T_c/T
deltaL = rhol/self.rhoc
deltaG = rhog/self.rhoc
idealL = self._phi0(tau, deltaL)
fioL = idealL["fio"]
fiot = idealL["fiot"]
idealG = self._phi0(tau, deltaG)
fioG = idealG["fio"]
firL = _phir(tau, deltaL, self._constants)
firdL = _phird(tau, deltaL, self._constants)
firtL = _phirt(tau, deltaL, self._constants)
soL = self.R*(tau*(fiot+firtL)-fioL-firL)
firG = _phir(tau, deltaG, self._constants)
firdG = _phird(tau, deltaG, self._constants)
firtG = _phirt(tau, deltaG, self._constants)
soG = self.R*(tau*(fiot+firtG)-fioG-firG)
Jl = rhol*(1+deltaL*firdL)
Jv = rhog*(1+deltaG*firdG)
K = firL-firG
x = (1./rho-1/rhol)/(1/rhog-1/rhol)
return (Jl-Jv,
Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
soL*(1-x)+soG*x - s)
for to in [To, 300, 400, 500, 600]:
rhoLo = self._Liquid_Density(to)
rhoGo = self._Vapor_Density(to)
sol = fsolve(f, [to, rhoLo, rhoGo], full_output=True)
T, rhoL, rhoG = sol[0]
x = (1./rho-1/rhoL)/(1/rhoG-1/rhoL)
if sol[2] == 1 and 0 <= x <= 1 and \
sum(abs(sol[1]["fvec"])) < 1e-5:
break
if sum(abs(sol[1]["fvec"])) > 1e-5:
raise RuntimeError(sol[3])
liquido = self._Helmholtz(rhoL, T)
vapor = self._Helmholtz(rhoG, T)
P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000
elif self._mode == "rhou":
delta = rho/rho_c
def f(T):
tau = T_c/T
ideal = self._phi0(tau, delta)
fiot = ideal["fiot"]
fird = _phird(tau, delta, self._constants)
firt = _phirt(tau, delta, self._constants)
Po = (1+delta*fird)*self.R*T*rho
ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
return ho-Po/rho-u
T = fsolve(f, To)[0]
rhol = self._Liquid_Density(T)
rhov = self._Vapor_Density(T)
if T == To or rhov <= rho <= rhol:
def f(parr):
T, rhol, rhog = parr
tau = T_c/T
deltaL = rhol/self.rhoc
deltaG = rhog/self.rhoc
ideal = self._phi0(tau, deltaL)
fiot = ideal["fiot"]
firL = _phir(tau, deltaL, self._constants)
firdL = _phird(tau, deltaL, self._constants)
firtL = _phirt(tau, deltaL, self._constants)
hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
firG = _phir(tau, deltaG, self._constants)
firdG = _phird(tau, deltaG, self._constants)
firtG = _phirt(tau, deltaG, self._constants)
hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)
Jl = rhol*(1+deltaL*firdL)
Jv = rhog*(1+deltaG*firdG)
K = firL-firG
x = (1./rho-1/rhol)/(1/rhog-1/rhol)
Ps = self.R*T*rhol*rhog/(rhol-rhog)*(K+log(rhol/rhog))
vu = hoG-Ps/rhog
lu = hoL-Ps/rhol
return (Jl-Jv,
Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
lu*(1-x)+vu*x - u)
for to in [To, 300, 400, 500, 600]:
rhoLo = self._Liquid_Density(to)
rhoGo = self._Vapor_Density(to)
sol = fsolve(f, [to, rhoLo, rhoGo], full_output=True)
T, rhoL, rhoG = sol[0]
x = (1./rho-1/rhoL)/(1/rhoG-1/rhoL)
if sol[2] == 1 and 0 <= x <= 1 and \
sum(abs(sol[1]["fvec"])) < 1e-5:
break
if sum(abs(sol[1]["fvec"])) > 1e-5:
raise RuntimeError(sol[3])
liquido = self._Helmholtz(rhoL, T)
vapor = self._Helmholtz(rhoG, T)
P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000
elif self._mode == "hs":
def f(parr):
rho, T = parr
delta = rho/rho_c
tau = T_c/T
ideal = self._phi0(tau, delta)
fio = ideal["fio"]
fiot = ideal["fiot"]
fird = _phird(tau, delta, self._constants)
fir = _phir(tau, delta, self._constants)
firt = _phirt(tau, delta, self._constants)
ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
so = self.R*(tau*(fiot+firt)-fio-fir)
return ho-h, so-s
rho, T = fsolve(f, [rhoo, To])
rhol = self._Liquid_Density(T)
rhov = self._Vapor_Density(T)
if rhov <= rho <= rhol:
def f(parr):
T, rhol, rhog, x = parr
tau = T_c/T
deltaL = rhol/self.rhoc
deltaG = rhog/self.rhoc
idealL = self._phi0(tau, deltaL)
fiot = idealL["fiot"]
fioL = idealL["fio"]
idealG = self._phi0(tau, deltaG)
fioG = idealG["fio"]
firL = _phir(tau, deltaL, self._constants)
firdL = _phird(tau, deltaL, self._constants)
firtL = _phirt(tau, deltaL, self._constants)
hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
soL = self.R*(tau*(fiot+firtL)-fioL-firL)
firG = _phir(tau, deltaG, self._constants)
firdG = _phird(tau, deltaG, self._constants)
firtG = _phirt(tau, deltaG, self._constants)
hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)
soG = self.R*(tau*(fiot+firtG)-fioG-firG)
Jl = rhol*(1+deltaL*firdL)
Jv = rhog*(1+deltaG*firdG)
K = firL-firG
return (Jl-Jv,
Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
hoL*(1-x)+hoG*x - h,
soL*(1-x)+soG*x - s)
for to in [To, 300, 400, 500, 600]:
rLo = self._Liquid_Density(to)
rGo = self._Vapor_Density(to)
sol = fsolve(f, [to, rLo, rGo, 0.5], full_output=True)
T, rhoL, rhoG, x = sol[0]
if sol[2] == 1 and 0 <= x <= 1 and \
sum(abs(sol[1]["fvec"])) < 1e-5:
break
if sum(abs(sol[1]["fvec"])) > 1e-5:
raise RuntimeError(sol[3])
liquido = self._Helmholtz(rhoL, T)
vapor = self._Helmholtz(rhoG, T)
P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000
elif self._mode == "hu":
def f(parr):
rho, T = parr
delta = rho/rho_c
tau = T_c/T
ideal = self._phi0(tau, delta)
fiot = ideal["fiot"]
fird = _phird(tau, delta, self._constants)
firt = _phirt(tau, delta, self._constants)
Po = (1+delta*fird)*self.R*T*rho
ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
return ho-Po/rho-u, ho-h
sol = fsolve(f, [rhoo, To], full_output=True)
rho, T = sol[0]
rhol = self._Liquid_Density(T)
rhov = self._Vapor_Density(T)
if sol[2] != 1 or rhov <= rho <= rhol:
def f(parr):
T, rhol, rhog, x = parr
tau = T_c/T
deltaL = rhol/self.rhoc
deltaG = rhog/self.rhoc
ideal = self._phi0(tau, deltaL)
fiot = ideal["fiot"]
firL = _phir(tau, deltaL, self._constants)
firdL = _phird(tau, deltaL, self._constants)
firtL = _phirt(tau, deltaL, self._constants)
hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
firG = _phir(tau, deltaG, self._constants)
firdG = _phird(tau, deltaG, self._constants)
firtG = _phirt(tau, deltaG, self._constants)
hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)
Jl = rhol*(1+deltaL*firdL)
Jv = rhog*(1+deltaG*firdG)
K = firL-firG
Ps = self.R*T*rhol*rhog/(rhol-rhog)*(K+log(rhol/rhog))
vu = hoG-Ps/rhog
lu = hoL-Ps/rhol
return (Jl-Jv,
Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
hoL*(1-x)+hoG*x - h,
lu*(1-x)+vu*x - u)
for to in [To, 300, 400, 500, 600]:
rLo = self._Liquid_Density(to)
rGo = self._Vapor_Density(to)
sol = fsolve(f, [to, rLo, rGo, 0.5], full_output=True)
T, rhoL, rhoG, x = sol[0]
if sol[2] == 1 and 0 <= x <= 1 and \
sum(abs(sol[1]["fvec"])) < 1e-5:
break
if sum(abs(sol[1]["fvec"])) > 1e-5:
raise RuntimeError(sol[3])
liquido = self._Helmholtz(rhoL, T)
vapor = self._Helmholtz(rhoG, T)
P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000
elif self._mode == "su":
def f(parr):
rho, T = parr
delta = rho/rho_c
tau = T_c/T
ideal = self._phi0(tau, delta)
fio = ideal["fio"]
fiot = ideal["fiot"]
fird = _phird(tau, delta, self._constants)
fir = _phir(tau, delta, self._constants)
firt = _phirt(tau, delta, self._constants)
ho = self.R*T*(1+tau*(fiot+firt)+delta*fird)
so = self.R*(tau*(fiot+firt)-fio-fir)
Po = (1+delta*fird)*self.R*T*rho
return ho-Po/rho-u, so-s
sol = fsolve(f, [rhoo, To], full_output=True)
rho, T = sol[0]
rhol = self._Liquid_Density(T)
rhov = self._Vapor_Density(T)
if sol[2] != 1 or rhov <= rho <= rhol:
def f(parr):
T, rhol, rhog, x = parr
tau = T_c/T
deltaL = rhol/self.rhoc
deltaG = rhog/self.rhoc
idealL = self._phi0(tau, deltaL)
fiot = idealL["fiot"]
fioL = idealL["fio"]
idealG = self._phi0(tau, deltaG)
fioG = idealG["fio"]
firL = _phir(tau, deltaL, self._constants)
firdL = _phird(tau, deltaL, self._constants)
firtL = _phirt(tau, deltaL, self._constants)
hoL = self.R*T*(1+tau*(fiot+firtL)+deltaL*firdL)
soL = self.R*(tau*(fiot+firtL)-fioL-firL)
firG = _phir(tau, deltaG, self._constants)
firdG = _phird(tau, deltaG, self._constants)
firtG = _phirt(tau, deltaG, self._constants)
hoG = self.R*T*(1+tau*(fiot+firtG)+deltaG*firdG)
soG = self.R*(tau*(fiot+firtG)-fioG-firG)
Jl = rhol*(1+deltaL*firdL)
Jv = rhog*(1+deltaG*firdG)
K = firL-firG
Ps = self.R*T*rhol*rhog/(rhol-rhog)*(K+log(rhol/rhog))
vu = hoG-Ps/rhog
lu = hoL-Ps/rhol
return (Jl-Jv,
Jl*(1/rhog-1/rhol)-log(rhol/rhog)-K,
soL*(1-x)+soG*x - s,
lu*(1-x)+vu*x - u)
for to in [To, 300, 400, 500, 600]:
rLo = self._Liquid_Density(to)
rGo = self._Vapor_Density(to)
sol = fsolve(f, [to, rLo, rGo, 0.5], full_output=True)
T, rhoL, rhoG, x = sol[0]
if sol[2] == 1 and 0 <= x <= 1 and \
sum(abs(sol[1]["fvec"])) < 1e-5:
break
if sum(abs(sol[1]["fvec"])) > 1e-5:
raise RuntimeError(sol[3])
liquido = self._Helmholtz(rhoL, T)
vapor = self._Helmholtz(rhoG, T)
P = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(
liquido["fir"]-vapor["fir"]+log(rhoL/rhoG))/1000
elif self._mode == "Trho":
if T < self.Tc:
rhov = self._Vapor_Density(T)
rhol = self._Liquid_Density(T)
if rhol > rho > rhov:
rhol, rhov, Ps = self._saturation(T)
if rhol > rho > rhov:
vapor = self._Helmholtz(rhov, T)
liquido = self._Helmholtz(rhol, T)
x = (1/rho-1/rhol)/(1/rhov-1/rhol)
P = Ps/1000
rho = float(rho)
T = float(T)
propiedades = self._Helmholtz(rho, T)
if T > self.Tc:
x = 1
elif x is None:
x = 0
if not P:
P = propiedades["P"]/1000.
elif self._mode == "Tx":
# Check input T in saturation range
if self.Tt > T or self.Tc < T or x > 1 or x < 0:
raise NotImplementedError("Incoming out of bound")
rhol, rhov, Ps = self._saturation(T)
vapor = self._Helmholtz(rhov, T)
liquido = self._Helmholtz(rhol, T)
if x == 0:
propiedades = liquido
elif x == 1:
propiedades = vapor
P = Ps/1000.
elif self._mode == "Px":
# Check input P in saturation range
if self.Pc < P or x > 1 or x < 0:
raise NotImplementedError("Incoming out of bound")
# Iterate over saturation routine to get T
def f(T):
rhol = self._Liquid_Density(T)
rhog = self._Vapor_Density(T)
deltaL = rhol/self.rhoc
deltaG = rhog/self.rhoc
tau = T_c/T
firL = _phir(tau, deltaL, self._constants)
firG = _phir(tau, deltaG, self._constants)
Ps = self.R*T*rhol*rhog/(rhol-rhog)*(
firL-firG+log(deltaL/deltaG))
return Ps/1000-P
if T0:
To = T0
elif self.name == "water":
To = _TSat_P(P)
else:
To = (self.Tc+self.Tt)/2
T = fsolve(f, To)[0]
rhol, rhov, Ps = self._saturation(T)
vapor = self._Helmholtz(rhov, T)
liquido = self._Helmholtz(rhol, T)
if x == 0:
propiedades = liquido
elif x == 1:
propiedades = vapor
self.T = T
self.Tr = T/self.Tc
self.P = P
self.Pr = self.P/self.Pc
self.x = x
if self._mode in ["Tx", "Px"] or 0 < x < 1:
region = 4
else:
region = 0
self.phase = getphase(self.Tc, self.Pc, self.T, self.P, self.x, region)
if 0 < x < 1:
rho = vapor["rho"]
else:
rho = propiedades["rho"]
self.Liquid = _fase()
self.Gas = _fase()
self.Vapor = self.Gas
if x == 0:
# liquid phase
self.fill(self.Liquid, propiedades)
self.fill(self, propiedades)
elif x == 1:
# vapor phase
self.fill(self.Gas, propiedades)
self.fill(self, propiedades)
else:
self.fill(self.Liquid, liquido)
self.fill(self.Gas, vapor)
self.v = x*self.Gas.v+(1-x)*self.Liquid.v
self.rho = 1./self.v
self.h = x*self.Gas.h+(1-x)*self.Liquid.h
self.s = x*self.Gas.s+(1-x)*self.Liquid.s
self.u = x*self.Gas.u+(1-x)*self.Liquid.u
self.a = x*self.Gas.a+(1-x)*self.Liquid.a
self.g = x*self.Gas.g+(1-x)*self.Liquid.g
self.Z = x*self.Gas.Z+(1-x)*self.Liquid.Z
self.f = x*self.Gas.f+(1-x)*self.Liquid.f
self.Z_rho = x*self.Gas.Z_rho+(1-x)*self.Liquid.Z_rho
self.IntP = x*self.Gas.IntP+(1-x)*self.Liquid.IntP
# Calculate special properties useful only for one phase
if self._mode in ("Px", "Tx") or (x < 1 and self.Tt <= T <= self.Tc):
self.sigma = self._surface(T)
else:
self.sigma = None
vir = self._virial(T)
self.virialB = vir["B"]/self.rhoc
self.virialC = vir["C"]/self.rhoc**2
if 0 < x < 1:
self.Hvap = vapor["h"]-liquido["h"]
self.Svap = vapor["s"]-liquido["s"]
else:
self.Hvap = None
self.Svap = None
self.invT = -1/self.T
# Ideal properties
self.v0 = self.R*self.T/self.P/1000
self.rho0 = 1./self.v0
cp0 = self._prop0(self.rho0, self.T)
self.h0 = cp0.h
self.u0 = self.h0-self.P*self.v0
self.s0 = cp0.s
self.a0 = self.u0-self.T*self.s0
self.g0 = self.h0-self.T*self.s0
self.cp0 = cp0.cp
self.cv0 = cp0.cv
self.cp0_cv = self.cp0/self.cv0
cp0.v = self.v0
self.gamma0 = -self.v0/self.P/1000*self.derivative("P", "v", "s", cp0)
[docs]
def fill(self, fase, estado):
"""Fill phase properties"""
fase.rho = estado["rho"]
fase.v = 1/fase.rho
fase.h = estado["h"]
fase.s = estado["s"]
fase.u = fase.h-self.P*1000*fase.v
fase.a = fase.u-self.T*fase.s
fase.g = fase.h-self.T*fase.s
fase.Z = self.P*fase.v/self.T/self.R*1e3
fase.fi = exp(estado["fir"]+estado["delta"]*estado["fird"]
- log(1+estado["delta"]*estado["fird"]))
fase.f = fase.fi*self.P
fase.cv = estado["cv"]
fase.rhoM = fase.rho/self.M
fase.hM = fase.h*self.M
fase.sM = fase.s*self.M
fase.uM = fase.u*self.M
fase.aM = fase.a*self.M
fase.gM = fase.g*self.M
fase.alfap = estado["alfap"]
fase.betap = estado["betap"]
fase.cp = self.derivative("h", "T", "P", fase)
fase.cp_cv = fase.cp/fase.cv
fase.w = (self.derivative("P", "rho", "s", fase)*1000)**0.5
fase.cvM = fase.cv*self.M
fase.cpM = fase.cp*self.M
fase.joule = self.derivative("T", "P", "h", fase)*1e3
fase.Gruneisen = fase.v/fase.cv*self.derivative("P", "T", "v", fase)
fase.alfav = self.derivative("v", "T", "P", fase)/fase.v
fase.kappa = -self.derivative("v", "P", "T", fase)/fase.v*1e3
fase.betas = self.derivative("T", "P", "s", fase)
fase.gamma = -fase.v/self.P \
* self.derivative("P", "v", "T", fase)*fase.cp_cv*1e-3
fase.kt = -fase.v/self.P*self.derivative("P", "v", "T", fase)*1e-3
fase.ks = -self.derivative("v", "P", "s", fase)/fase.v*1e3
fase.Kt = -fase.v*self.derivative("P", "v", "s", fase)*1e-3
fase.Ks = -fase.v*self.derivative("P", "v", "T", fase)*1e-3
fase.dhdT_rho = self.derivative("h", "T", "rho", fase)
fase.dhdT_P = self.derivative("h", "T", "P", fase)
fase.dhdP_T = self.derivative("h", "P", "T", fase)*1e3
fase.dhdP_rho = self.derivative("h", "P", "rho", fase)*1e3
fase.dhdrho_T = self.derivative("h", "rho", "T", fase)
fase.dhdrho_P = self.derivative("h", "rho", "P", fase)
fase.dpdT_rho = self.derivative("P", "T", "rho", fase)*1e-3
fase.dpdrho_T = self.derivative("P", "rho", "T", fase)*1e-3
fase.drhodP_T = self.derivative("rho", "P", "T", fase)*1e3
fase.drhodT_P = self.derivative("rho", "T", "P", fase)
fase.Z_rho = (fase.Z-1)/fase.rho
fase.IntP = self.T*self.derivative("P", "T", "rho", fase)*1e-3-self.P
fase.hInput = fase.v*self.derivative("h", "v", "P", fase)
fase.mu = self._visco(fase.rho, self.T, fase)
fase.k = self._thermo(fase.rho, self.T, fase)
fase.nu = fase.mu/fase.rho
fase.alfa = fase.k/1000/fase.rho/fase.cp
fase.Prandt = fase.mu*fase.cp*1000/fase.k
if self.name == "water":
try:
fase.epsilon = _Dielectric(fase.rho, self.T)
except NotImplementedError:
fase.epsilon = None
try:
fase.n = _Refractive(fase.rho, self.T, self.kwargs["l"])
except NotImplementedError:
fase.n = None
else:
fase.epsilon = None
fase.n = None
[docs]
def derivative(self, z, x, y, fase):
"""
Wrapper derivative for custom derived properties
where x, y, z can be: P, T, v, rho, u, h, s, g, a
"""
return deriv_H(self, z, x, y, fase)
[docs]
def _saturation(self, T):
"""Saturation calculation for two phase search"""
# Added support for superancillary equation to try speed up saturation
# iteration
# Using generated data from https://github.com/usnistgov/fastchebpure
Ps, rhoL, rhoG = None, None, None
anc_file = os.path.join(
os.path.dirname(__file__), "%s_anc.json" %self.__class__.__name__)
if os.path.isfile(anc_file):
with open(anc_file) as anc_file:
dat = json.load(anc_file)
for coefs in dat["jexpansions_p"]:
coef = coefs["coef"]
xmin = coefs["xmin"]
xmax = coefs["xmax"]
if xmin <= T <= xmax:
c = chebyshev.Chebyshev(coef, domain=(xmin, xmax))
Ps = c(T)/1e3
break
for coefs in dat["jexpansions_rhoL"]:
coef = coefs["coef"]
xmin = coefs["xmin"]
xmax = coefs["xmax"]
if xmin <= T <= xmax:
c = chebyshev.Chebyshev(coef, domain=(xmin, xmax))
rhoL = c(T)*self.M/1e3
break
for coefs in dat["jexpansions_rhoV"]:
coef = coefs["coef"]
xmin = coefs["xmin"]
xmax = coefs["xmax"]
if xmin <= T <= xmax:
c = chebyshev.Chebyshev(coef, domain=(xmin, xmax))
rhoG = c(T)*self.M/1e3
break
if Ps is None:
rho_c = self._constants.get("rhoref", self.rhoc)
T_c = self._constants.get("Tref", self.Tc)
T = min(T, T_c)
tau = T_c/T
rhoLo = self._Liquid_Density(T)
rhoGo = self._Vapor_Density(T)
def f(parr):
rhol, rhog = parr
deltaL = rhol/rho_c
deltaG = rhog/rho_c
phirL = _phir(tau, deltaL, self._constants)
phirG = _phir(tau, deltaG, self._constants)
phirdL = _phird(tau, deltaL, self._constants)
phirdG = _phird(tau, deltaG, self._constants)
Jl = deltaL*(1+deltaL*phirdL)
Jv = deltaG*(1+deltaG*phirdG)
Kl = deltaL*phirdL+phirL+log(deltaL)
Kv = deltaG*phirdG+phirG+log(deltaG)
return Kv-Kl, Jv-Jl
rhoL, rhoG = fsolve(f, [rhoLo, rhoGo])
if rhoL == rhoG:
Ps = self.Pc*1e3
else:
deltaL = rhoL/self.rhoc
deltaG = rhoG/self.rhoc
firL = _phir(tau, deltaL, self._constants)
firG = _phir(tau, deltaG, self._constants)
Ps = self.R*T*rhoL*rhoG/(rhoL-rhoG)*(firL-firG+log(deltaL/deltaG))
# print(rhoL, rhoG, Ps)
return rhoL, rhoG, Ps
[docs]
def _Helmholtz(self, rho, T):
"""Calculated properties from helmholtz free energy and derivatives
Parameters
----------
rho : float
Density, [kg/m³]
T : float
Temperature, [K]
Returns
-------
prop : dict
Dictionary with calculated properties:
* fir: [-]
* fird: ∂fir/∂δ|τ
* firdd: ∂²fir/∂δ²|τ
* delta: Reducen density rho/rhoc, [-]
* P: Pressure, [kPa]
* h: Enthalpy, [kJ/kg]
* s: Entropy, [kJ/kgK]
* cv: Isochoric specific heat, [kJ/kgK]
* alfav: Thermal expansion coefficient, [1/K]
* betap: Isothermal compressibility, [1/kPa]
References
----------
IAPWS, Revised Release on the IAPWS Formulation 1995 for the
Thermodynamic Properties of Ordinary Water Substance for General and
Scientific Use, September 2016, Table 3
http://www.iapws.org/relguide/IAPWS-95.html
"""
if isinstance(rho, ndarray):
rho = rho[0]
if isinstance(T, ndarray):
T = T[0]
if rho < 0:
rho = 1e-20
T = max(T, 50)
rho_c = self._constants.get("rhoref", self.rhoc)
T_c = self._constants.get("Tref", self.Tc)
delta = rho/rho_c
tau = T_c/T
ideal = self._phi0(tau, delta)
fio = ideal["fio"]
fiot = ideal["fiot"]
fiott = ideal["fiott"]
res = self._phir(tau, delta)
fir = res["fir"]
firt = res["firt"]
firtt = res["firtt"]
fird = res["fird"]
firdd = res["firdd"]
firdt = res["firdt"]
propiedades = {}
propiedades["fir"] = fir
propiedades["fird"] = fird
propiedades["firdd"] = firdd
propiedades["delta"] = delta
propiedades["rho"] = rho
propiedades["P"] = (1+delta*fird)*self.R*T*rho
propiedades["h"] = self.R*T*(1+tau*(fiot+firt)+delta*fird)
propiedades["s"] = self.R*(tau*(fiot+firt)-fio-fir)
propiedades["cv"] = -self.R*tau**2*(fiott+firtt)
propiedades["alfap"] = (1-delta*tau*firdt/(1+delta*fird))/T
propiedades["betap"] = rho*(
1+(delta*fird+delta**2*firdd)/(1+delta*fird))
return propiedades
[docs]
def _prop0(self, rho, T):
"""Ideal gas properties"""
rho_c = self._constants.get("rhoref", self.rhoc)
T_c = self._constants.get("Tref", self.Tc)
delta = rho/rho_c
tau = T_c/T
ideal = self._phi0(tau, delta)
fio = ideal["fio"]
fiot = ideal["fiot"]
fiott = ideal["fiott"]
propiedades = _fase()
propiedades.h = self.R*T*(1+tau*fiot)
propiedades.s = self.R*(tau*fiot-fio)
propiedades.cv = -self.R*tau**2*fiott
propiedades.cp = self.R*(-tau**2*fiott+1)
propiedades.alfap = 1/T
propiedades.betap = rho
return propiedades
[docs]
def _phi0(self, tau, delta):
"""Ideal gas Helmholtz free energy and derivatives
Parameters
----------
tau : float
Inverse reduced temperature Tc/T, [-]
delta : float
Reduced density rho/rhoc, [-]
Returns
-------
prop : dictionary with ideal adimensional helmholtz energy and deriv
fio, [-]
fiot: ∂fio/∂τ|δ
fiod: ∂fio/∂δ|τ
fiott: ∂²fio/∂τ²|δ
fiodt: ∂²fio/∂τ∂δ
fiodd: ∂²fio/∂δ²|τ
References
----------
IAPWS, Revised Release on the IAPWS Formulation 1995 for the
Thermodynamic Properties of Ordinary Water Substance for General and
Scientific Use, September 2016, Table 4
http://www.iapws.org/relguide/IAPWS-95.html
"""
Fi0 = self.Fi0
fio = Fi0["ao_log"][0]*log(delta)+Fi0["ao_log"][1]*log(tau)
fiot = +Fi0["ao_log"][1]/tau
fiott = -Fi0["ao_log"][1]/tau**2
fiod = 1/delta
fiodd = -1/delta**2
fiodt = 0
for n, t in zip(Fi0["ao_pow"], Fi0["pow"]):
fio += n*tau**t
if t != 0:
fiot += t*n*tau**(t-1)
if t not in [0, 1]:
fiott += n*t*(t-1)*tau**(t-2)
for n, t in zip(Fi0["ao_exp"], Fi0["titao"]):
fio += n*log(1-exp(-tau*t))
fiot += n*t*((1-exp(-t*tau))**-1-1)
fiott -= n*t**2*exp(-t*tau)*(1-exp(-t*tau))**-2
# Extension to especial terms of air
if "ao_exp2" in Fi0:
for n, g, C in zip(Fi0["ao_exp2"], Fi0["titao2"], Fi0["sum2"]):
fio += n*log(C+exp(g*tau))
fiot += n*g/(C*exp(-g*tau)+1)
fiott += C*n*g**2*exp(-g*tau)/(C*exp(-g*tau)+1)**2
prop = {}
prop["fio"] = fio
prop["fiot"] = fiot
prop["fiott"] = fiott
prop["fiod"] = fiod
prop["fiodd"] = fiodd
prop["fiodt"] = fiodt
return prop
[docs]
def _phir(self, tau, delta):
"""Residual contribution to the free Helmholtz energy
Parameters
----------
tau : float
Inverse reduced temperature Tc/T, [-]
delta : float
Reduced density rho/rhoc, [-]
Returns
-------
prop : dict
Dictionary with residual adimensional helmholtz energy and deriv:
* fir
* firt: ∂fir/∂τ|δ,x
* fird: ∂fir/∂δ|τ,x
* firtt: ∂²fir/∂τ²|δ,x
* firdt: ∂²fir/∂τ∂δ|x
* firdd: ∂²fir/∂δ²|τ,x
References
----------
IAPWS, Revised Release on the IAPWS Formulation 1995 for the
Thermodynamic Properties of Ordinary Water Substance for General and
Scientific Use, September 2016, Table 5
http://www.iapws.org/relguide/IAPWS-95.html
"""
fir = fird = firdd = firt = firtt = firdt = 0
# Polinomial terms
nr1 = self._constants.get("nr1", [])
d1 = self._constants.get("d1", [])
t1 = self._constants.get("t1", [])
for n, d, t in zip(nr1, d1, t1):
fir += n*delta**d*tau**t
fird += n*d*delta**(d-1)*tau**t
firdd += n*d*(d-1)*delta**(d-2)*tau**t
firt += n*t*delta**d*tau**(t-1)
firtt += n*t*(t-1)*delta**d*tau**(t-2)
firdt += n*t*d*delta**(d-1)*tau**(t-1)
# Exponential terms
nr2 = self._constants.get("nr2", [])
d2 = self._constants.get("d2", [])
g2 = self._constants.get("gamma2", [])
t2 = self._constants.get("t2", [])
c2 = self._constants.get("c2", [])
for n, d, g, t, c in zip(nr2, d2, g2, t2, c2):
fir += n*delta**d*tau**t*exp(-g*delta**c)
fird += n*exp(-g*delta**c)*delta**(d-1)*tau**t*(d-g*c*delta**c)
firdd += n*exp(-g*delta**c)*delta**(d-2)*tau**t * \
((d-g*c*delta**c)*(d-1-g*c*delta**c)-g**2*c**2*delta**c)
firt += n*t*delta**d*tau**(t-1)*exp(-g*delta**c)
firtt += n*t*(t-1)*delta**d*tau**(t-2)*exp(-g*delta**c)
firdt += n*t*delta**(d-1)*tau**(t-1)*(d-g*c*delta**c)*exp(
-g*delta**c)
# Gaussian terms
nr3 = self._constants.get("nr3", [])
d3 = self._constants.get("d3", [])
t3 = self._constants.get("t3", [])
a3 = self._constants.get("alfa3", [])
e3 = self._constants.get("epsilon3", [])
b3 = self._constants.get("beta3", [])
g3 = self._constants.get("gamma3", [])
for n, d, t, a, e, b, g in zip(nr3, d3, t3, a3, e3, b3, g3):
fir += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)
fird += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
d/delta-2*a*(delta-e))
firdd += n*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
-2*a*delta**d + 4*a**2*delta**d*(delta-e)**2
- 4*d*a*delta**(d-1)*(delta-e) + d*(d-1)*delta**(d-2))
firt += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
t/tau-2*b*(tau-g))
firtt += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
(t/tau-2*b*(tau-g))**2-t/tau**2-2*b)
firdt += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
t/tau-2*b*(tau-g))*(d/delta-2*a*(delta-e))
# Non analitic terms
nr4 = self._constants.get("nr4", [])
a4 = self._constants.get("a4", [])
b4 = self._constants.get("b4", [])
Ai = self._constants.get("A", [])
Bi = self._constants.get("B", [])
Ci = self._constants.get("C", [])
Di = self._constants.get("D", [])
bt4 = self._constants.get("beta4", [])
for n, a, b, A, B, C, D, bt in zip(nr4, a4, b4, Ai, Bi, Ci, Di, bt4):
Tita = (1-tau)+A*((delta-1)**2)**(0.5/bt)
F = exp(-C*(delta-1)**2-D*(tau-1)**2)
Fd = -2*C*F*(delta-1)
Fdd = 2*C*F*(2*C*(delta-1)**2-1)
Ft = -2*D*F*(tau-1)
Ftt = 2*D*F*(2*D*(tau-1)**2-1)
Fdt = 4*C*D*F*(delta-1)*(tau-1)
Delta = Tita**2+B*((delta-1)**2)**a
Deltad = (delta-1)*(A*Tita*2/bt*((delta-1)**2)**(0.5/bt-1)
+ 2*B*a*((delta-1)**2)**(a-1))
if delta == 1:
Deltadd = 0
else:
Deltadd = Deltad/(delta-1)+(delta-1)**2*(
4*B*a*(a-1)*((delta-1)**2)**(a-2)
+ 2*A**2/bt**2*(((delta-1)**2)**(0.5/bt-1))**2
+ A*Tita*4/bt*(0.5/bt-1)*((delta-1)**2)**(0.5/bt-2))
# Avoid ZeroDivisionError at critical point
if Delta == 0:
DeltaBd = 0
DeltaBdd = 0
DeltaBt = 0
DeltaBtt = 0
DeltaBdt = 0
else:
DeltaBd = b*Delta**(b-1)*Deltad
DeltaBdd = b*(Delta**(b-1)*Deltadd+(b-1)*Delta**(b-2)*Deltad**2)
DeltaBt = -2*Tita*b*Delta**(b-1)
DeltaBtt = 2*b*Delta**(b-1)+4*Tita**2*b*(b-1)*Delta**(b-2)
DeltaBdt = -A*b*2/bt*Delta**(b-1)*(delta-1)*((delta-1)**2)**(
0.5/bt-1)-2*Tita*b*(b-1)*Delta**(b-2)*Deltad
fir += n*Delta**b*delta*F
fird += n*(Delta**b*(F+delta*Fd)+DeltaBd*delta*F)
firdd += n*(Delta**b*(2*Fd+delta*Fdd) + 2*DeltaBd*(F+delta*Fd)
+ DeltaBdd*delta*F)
firt += n*delta*(DeltaBt*F+Delta**b*Ft)
firtt += n*delta*(DeltaBtt*F+2*DeltaBt*Ft+Delta**b*Ftt)
firdt += n*(Delta**b*(Ft+delta*Fdt)+delta*DeltaBd*Ft
+ DeltaBt*(F+delta*Fd)+DeltaBdt*delta*F)
prop = {}
prop["fir"] = fir
prop["firt"] = firt
prop["firtt"] = firtt
prop["fird"] = fird
prop["firdd"] = firdd
prop["firdt"] = firdt
return prop
[docs]
def _virial(self, T):
"""Virial coefficient
Parameters
----------
T : float
Temperature [K]
Returns
-------
prop : dict
Dictionary with residual adimensional helmholtz energy:
* B: ∂fir/∂δ|δ->0
* C: ∂²fir/∂δ²|δ->0
"""
T_c = self._constants.get("Tref", self.Tc)
tau = T_c/T
B = C = 0
delta = 1e-200
# Polinomial terms
nr1 = self._constants.get("nr1", [])
d1 = self._constants.get("d1", [])
t1 = self._constants.get("t1", [])
for n, d, t in zip(nr1, d1, t1):
B += n*d*delta**(d-1)*tau**t
C += n*d*(d-1)*delta**(d-2)*tau**t
# Exponential terms
nr2 = self._constants.get("nr2", [])
d2 = self._constants.get("d2", [])
g2 = self._constants.get("gamma2", [])
t2 = self._constants.get("t2", [])
c2 = self._constants.get("c2", [])
for n, d, g, t, c in zip(nr2, d2, g2, t2, c2):
B += n*exp(-g*delta**c)*delta**(d-1)*tau**t*(d-g*c*delta**c)
C += n*exp(-g*delta**c)*(delta**(d-2)*tau**t*(
(d-g*c*delta**c)*(d-1-g*c*delta**c)-g**2*c**2*delta**c))
# Gaussian terms
nr3 = self._constants.get("nr3", [])
d3 = self._constants.get("d3", [])
t3 = self._constants.get("t3", [])
a3 = self._constants.get("alfa3", [])
e3 = self._constants.get("epsilon3", [])
b3 = self._constants.get("beta3", [])
g3 = self._constants.get("gamma3", [])
for n, d, t, a, e, b, g in zip(nr3, d3, t3, a3, e3, b3, g3):
B += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
d/delta-2*a*(delta-e))
C += n*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
-2*a*delta**d+4*a**2*delta**d*(
delta-e)**2-4*d*a*delta**2*(
delta-e)+d*2*delta)
# Non analitic terms
nr4 = self._constants.get("nr4", [])
a4 = self._constants.get("a4", [])
b4 = self._constants.get("b4", [])
Ai = self._constants.get("A", [])
Bi = self._constants.get("B", [])
Ci = self._constants.get("C", [])
Di = self._constants.get("D", [])
bt4 = self._constants.get("beta4", [])
for n, a, b, A, B_, C_, D, bt in zip(nr4, a4, b4, Ai, Bi, Ci, Di, bt4):
Tita = (1-tau)+A*((delta-1)**2)**(0.5/bt)
Delta = Tita**2+B_*((delta-1)**2)**a
Deltad = (delta-1)*(A*Tita*2/bt*((delta-1)**2)**(
0.5/bt-1)+2*B_*a*((delta-1)**2)**(a-1))
Deltadd = Deltad/(delta-1) + (delta-1)**2*(
4*B_*a*(a-1)*((delta-1)**2)**(a-2)
+ 2*A**2/bt**2*(((delta-1)**2)**(0.5/bt-1))**2
+ A*Tita*4/bt*(0.5/bt-1)*((delta-1)**2)**(0.5/bt-2))
DeltaBd = b*Delta**(b-1)*Deltad
DeltaBdd = b*(Delta**(b-1)*Deltadd+(b-1)*Delta**(b-2)*Deltad**2)
F = exp(-C_*(delta-1)**2-D*(tau-1)**2)
Fd = -2*C_*F*(delta-1)
Fdd = 2*C_*F*(2*C_*(delta-1)**2-1)
B += n*(Delta**b*(F+delta*Fd)+DeltaBd*delta*F)
C += n*(Delta**b*(2*Fd+delta*Fdd)+2*DeltaBd*(F+delta*Fd)
+ DeltaBdd*delta*F)
prop = {}
prop["B"] = B
prop["C"] = C
return prop
[docs]
def _derivDimensional(self, rho, T):
"""Calcule the dimensional form or Helmholtz free energy derivatives
Parameters
----------
rho : float
Density, [kg/m³]
T : float
Temperature, [K]
Returns
-------
prop : dict
Dictionary with residual helmholtz energy and derivatives:
* fir, [kJ/kg]
* firt: ∂fir/∂T|ρ, [kJ/kgK]
* fird: ∂fir/∂ρ|T, [kJ/m³kg²]
* firtt: ∂²fir/∂T²|ρ, [kJ/kgK²]
* firdt: ∂²fir/∂T∂ρ, [kJ/m³kg²K]
* firdd: ∂²fir/∂ρ²|T, [kJ/m⁶kg]
References
----------
IAPWS, Guideline on an Equation of State for Humid Air in Contact with
Seawater and Ice, Consistent with the IAPWS Formulation 2008 for the
Thermodynamic Properties of Seawater, Table 7,
http://www.iapws.org/relguide/SeaAir.html
"""
if not rho:
prop = {}
prop["fir"] = 0
prop["firt"] = 0
prop["fird"] = 0
prop["firtt"] = 0
prop["firdt"] = 0
prop["firdd"] = 0
return prop
R = self._constants.get("R")/self._constants.get("M", self.M)
rho_c = self._constants.get("rhoref", self.rhoc)
T_c = self._constants.get("Tref", self.Tc)
delta = rho/rho_c
tau = T_c/T
ideal = self._phi0(tau, delta)
fio = ideal["fio"]
fiot = ideal["fiot"]
fiott = ideal["fiott"]
fiod = ideal["fiod"]
fiodd = ideal["fiodd"]
res = self._phir(tau, delta)
fir = res["fir"]
firt = res["firt"]
firtt = res["firtt"]
fird = res["fird"]
firdd = res["firdd"]
firdt = res["firdt"]
prop = {}
prop["fir"] = R*T*(fio+fir)
prop["firt"] = R*(fio+fir-(fiot+firt)*tau)
prop["fird"] = R*T/rho_c*(fiod+fird)
prop["firtt"] = R*tau**2/T*(fiott+firtt)
prop["firdt"] = R/rho_c*(fiod+fird-firdt*tau)
prop["firdd"] = R*T/rho_c**2*(fiodd+firdd)
return prop
[docs]
def _surface(self, T):
"""Generic equation for the surface tension
Parameters
----------
T : float
Temperature, [K]
Returns
-------
σ : float
Surface tension, [N/m]
Notes
-----
Need a _surf dict in the derived class with the parameters keys:
sigma: coefficient
exp: exponent
"""
tau = 1-T/self.Tc
sigma = 0
for n, t in zip(self._surf["sigma"], self._surf["exp"]):
sigma += n*tau**t
return sigma
[docs]
@classmethod
def _Vapor_Pressure(cls, T):
"""Auxiliary equation for the vapour pressure
Parameters
----------
T : float
Temperature, [K]
Returns
-------
Pv : float
Vapour pressure, [Pa]
References
----------
IAPWS, Revised Supplementary Release on Saturation Properties of
Ordinary Water Substance September 1992,
http://www.iapws.org/relguide/Supp-sat.html, Eq.1
"""
if T < cls.Tt:
T = cls.Tt
elif T > cls.Tc:
T = cls.Tc
Tita = 1-T/cls.Tc
suma = 0
for n, x in zip(cls._Pv["ao"], cls._Pv["exp"]):
suma += n*Tita**x
Pr = exp(cls.Tc/T*suma)
Pv = Pr*cls.Pc
return Pv
[docs]
@classmethod
def _Liquid_Density(cls, T):
"""Auxiliary equation for the density of saturated liquid
Parameters
----------
T : float
Temperature, [K]
Returns
-------
rho : float
Saturated liquid density, [kg/m³]
References
----------
IAPWS, Revised Supplementary Release on Saturation Properties of
Ordinary Water Substance September 1992,
http://www.iapws.org/relguide/Supp-sat.html, Eq.2
"""
if T < cls.Tt:
T = cls.Tt
elif T > cls.Tc:
T = cls.Tc
eq = cls._rhoL["eq"]
Tita = 1-T/cls.Tc
if eq == 2:
Tita = Tita**(1./3)
suma = 0
for n, x in zip(cls._rhoL["ao"], cls._rhoL["exp"]):
suma += n*Tita**x
Pr = suma+1
rho = Pr*cls.rhoc
return rho
[docs]
@classmethod
def _Vapor_Density(cls, T):
"""Auxiliary equation for the density of saturated vapor
Parameters
----------
T : float
Temperature, [K]
Returns
-------
rho : float
Saturated vapor density, [kg/m³]
References
----------
IAPWS, Revised Supplementary Release on Saturation Properties of
Ordinary Water Substance September 1992,
http://www.iapws.org/relguide/Supp-sat.html, Eq.3
"""
if T < cls.Tt:
T = cls.Tt
elif T > cls.Tc:
T = cls.Tc
eq = cls._rhoG["eq"]
Tita = 1-T/cls.Tc
if eq == 4:
Tita = Tita**(1./3)
suma = 0
for n, x in zip(cls._rhoG["ao"], cls._rhoG["exp"]):
suma += n*Tita**x
Pr = exp(suma)
rho = Pr*cls.rhoc
return rho
[docs]
@classmethod
def _dPdT_sat(cls, T):
"""Auxiliary equation for the dP/dT along saturation line
Parameters
----------
T : float
Temperature, [K]
Returns
-------
dPdT : float
dPdT, [MPa/K]
References
----------
IAPWS, Revised Supplementary Release on Saturation Properties of
Ordinary Water Substance September 1992,
http://www.iapws.org/relguide/Supp-sat.html, derived from Eq.1
"""
Tita = 1-T/cls.Tc
suma1 = 0
suma2 = 0
for n, x in zip(cls._Pv["ao"], cls._Pv["exp"]):
suma1 -= n*x*Tita**(x-1)/cls.Tc
suma2 += n*Tita**x
Pr = (cls.Tc*suma1/T-cls.Tc/T**2*suma2)*exp(cls.Tc/T*suma2)
dPdT = Pr*cls.Pc
return dPdT
[docs]
def mainClassDoc():
"""
Function decorator used to automatic adiction of base class MEoS in
subclass __doc__
"""
def decorator(f):
# __doc__ is only writable in python3.
# The doc build must be done with python3 so this snnippet do the work
py_version = platform.python_version_tuple()[0]
if py_version == 3:
doc = f.__doc__.split(os.linesep)
try:
ind = doc.index("")
except ValueError:
ind = 1
doc1 = os.linesep.join(doc[:ind])
doc3 = os.linesep.join(doc[ind:])
doc2 = os.linesep.join(MEoS.__doc__.split(os.linesep)[3:])
f.__doc__ = doc1 + os.linesep + os.linesep + \
doc2 + os.linesep + os.linesep + doc3
return f
return decorator
[docs]
@mainClassDoc()
class IAPWS95(MEoS):
"""Implementation of IAPWS Formulation 1995 for ordinary water substance,
(revised release of 2016), for internal procedures, see MEoS base class
Examples
--------
>>> water=IAPWS95(T=300, rho=996.5560)
>>> water.P, water.cv, water.w, water.s
0.0992418350 4.13018112 1501.51914 0.393062643
>>> water=IAPWS95(T=500, rho=0.435)
>>> water.P, water.cv, water.w, water.s
0.0999679423 1.50817541 548.31425 7.944882714
>>> water=IAPWS95(T=900., P=700)
>>> water.rho, water.cv, water.w, water.s
870.7690 2.66422350 2019.33608 4.17223802
>>> water=IAPWS95(T=300., P=0.1)
>>> water.P, water.rho, water.h, water.s, water.cp, water.w, water.virialB
0.10000 996.56 112.65 0.39306 4.1806 1501.5 -0.066682
>>> water=IAPWS95(T=500., P=0.1)
>>> water.P, water.rho, water.h, water.s, water.cp, water.w, water.virialB
0.10000 0.43514 2928.6 7.9447 1.9813 548.31 -0.0094137
>>> water=IAPWS95(T=450., x=0.5)
>>> water.T, water.P, water.rho, water.h, water.s, water.virialB
450.00 0.93220 9.5723 1761.8 4.3589 -0.013028
>>> water=IAPWS95(P=1.5, rho=1000.)
>>> water.T, water.rho, water.h, water.s, water.cp, water.w, water.virialB
286.44 1000.0 57.253 0.19931 4.1855 1462.1 -0.085566
>>> water=IAPWS95(h=3000, s=8.)
>>> water.T, water.P, water.h, water.s, water.cp, water.w, water.virialB
536.24 0.11970 3000.0 8.0000 1.9984 567.04 -0.0076606
>>> water=IAPWS95(h=150, s=0.4)
>>> water.T, water.P, water.rho, water.h, water.s, water.cp, water.w
301.27 35.50549 1011.48 150.00 0.40000 4.0932 1564.1
>>> water=IAPWS95(T=450., rho=300)
>>> water.T, water.P, water.rho, water.h, water.s, water.x, water.virialB
450.00 0.93220 300.00 770.82 2.1568 0.010693 -0.013028
>>> water=IAPWS95(rho=300., P=0.1)
>>> water.T, water.P, water.rho, water.h, water.s, water.x, water.virialB
372.76 0.10000 300.00 420.56 1.3110 0.0013528 -0.025144
>>> water=IAPWS95(h=1500., P=0.1)
>>> water.T, water.P, water.rho, water.h, water.s, water.x, water.virialB
372.76 0.10000 1.2303 1500.0 4.2068 0.47952 -0.025144
>>> water=IAPWS95(s=5., P=3.5)
>>> water.T, water.P, water.rho, water.h, water.s, water.x, water.virialB
515.71 3.5000 25.912 2222.8 5.0000 0.66921 -0.0085877
>>> water=IAPWS95(T=500., u=900)
>>> water.P, water.rho, water.u, water.h, water.s, water.cp, water.w
108.21 903.62 900.00 1019.8 2.4271 4.1751 1576.0
>>> water=IAPWS95(P=0.3, u=1550.)
>>> water.T, water.P, water.rho, water.u, water.h, water.s, water.x
406.67 0.30000 3.3029 1550.0 1640.8 4.3260 0.49893
>>> water=IAPWS95(rho=300, h=1000.)
>>> water.T, water.P, water.rho, water.u, water.h, water.s, water.x
494.92 2.3991 300.00 992.00 1000.0 2.6315 0.026071
>>> water=IAPWS95(rho=30, s=8.)
>>> water.T, water.P, water.rho, water.u, water.h, water.s, water.cp
1562.42 21.671 30.000 4628.5 5350.9 8.0000 2.7190
>>> water=IAPWS95(rho=30, s=4.)
>>> water.T, water.P, water.rho, water.u, water.h, water.s, water.x
495.00 2.4029 30.000 1597.3 1677.4 4.0000 0.39218
>>> water=IAPWS95(rho=300, u=1000.)
>>> water.T, water.P, water.rho, water.u, water.h, water.s, water.x
496.44 2.4691 300.000 1000.0 1008.2 2.6476 0.02680
>>> water=IAPWS95(s=3., h=1000.)
>>> water.T, water.P, water.rho, water.u, water.h, water.s, water.x
345.73 0.034850 0.73526 952.60 1000.0 3.0000 0.29920
>>> water=IAPWS95(u=995., h=1000.)
>>> water.T, water.P, water.rho, water.u, water.h, water.s, water.x
501.89 2.7329 546.58 995.00 1000.0 2.6298 0.00866
>>> water=IAPWS95(u=1000., s=3.)
>>> water.T, water.P, water.rho, water.u, water.h, water.s, water.x
371.24 0.094712 1.99072 1000.00 1047.6 3.0000 0.28144
References
----------
IAPWS, Revised Release on the IAPWS Formulation 1995 for the Thermodynamic
Properties of Ordinary Water Substance for General and Scientific Use,
September 2016, http://www.iapws.org/relguide/IAPWS-95.html
IAPWS, Revised Supplementary Release on Saturation Properties of Ordinary
Water Substance September 1992, http://www.iapws.org/relguide/Supp-sat.html
IAPWS, Guideline on a Low-Temperature Extension of the IAPWS-95 Formulation
for Water Vapor, http://www.iapws.org/relguide/LowT.html
IAPWS, Revised Advisory Note No. 3: Thermodynamic Derivatives from IAPWS
Formulations, http://www.iapws.org/relguide/Advise3.pdf
"""
name = "water"
CASNumber = "7732-18-5"
formula = "H2O"
synonym = "R-718"
Tc = Tc
rhoc = rhoc
Pc = Pc
M = M
Tt = 273.16
Tb = 373.1243
f_acent = 0.3443
momentoDipolar = 1.855
Fi0 = {"ao_log": [1, 3.00632],
"pow": [0, 1],
"ao_pow": [-8.3204464837497, 6.6832105275932],
"ao_exp": [0.012436, 0.97315, 1.2795, 0.96956, 0.24873],
"titao": [1.28728967, 3.53734222, 7.74073708, 9.24437796,
27.5075105]}
_constants = {
"R": 8.314371357587,
"nr1": [0.12533547935523e-1, 0.78957634722828e1, -0.87803203303561e1,
0.31802509345418, -0.26145533859358, -0.78199751687981e-2,
0.88089493102134e-2],
"d1": [1, 1, 1, 2, 2, 3, 4],
"t1": [-0.5, 0.875, 1, 0.5, 0.75, 0.375, 1],
"nr2": [-0.66856572307965, 0.20433810950965, -0.66212605039687e-4,
-0.19232721156002, -0.25709043003438, 0.16074868486251,
-0.40092828925807e-1, .39343422603254e-6, -0.75941377088144e-5,
0.56250979351888e-3, -0.15608652257135e-4, 0.11537996422951e-8,
.36582165144204e-6, -.13251180074668e-11, -.62639586912454e-9,
-0.10793600908932, 0.17611491008752e-1, 0.22132295167546,
-0.40247669763528, 0.58083399985759, 0.49969146990806e-2,
-0.31358700712549e-1, -0.74315929710341, 0.47807329915480,
0.20527940895948e-1, -0.13636435110343, 0.14180634400617e-1,
0.83326504880713e-2, -0.29052336009585e-1, 0.38615085574206e-1,
-0.20393486513704e-1, -0.16554050063734e-2, .19955571979541e-2,
0.15870308324157e-3, -0.16388568342530e-4, 0.43613615723811e-1,
0.34994005463765e-1, -0.76788197844621e-1, 0.22446277332006e-1,
-0.62689710414685e-4, -0.55711118565645e-9, -0.19905718354408,
0.31777497330738, -0.11841182425981],
"c2": [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 6, 6,
6, 6],
"d2": [1, 1, 1, 2, 2, 3, 4, 4, 5, 7, 9, 10, 11, 13, 15, 1, 2, 2, 2, 3,
4, 4, 4, 5, 6, 6, 7, 9, 9, 9, 9, 9, 10, 10, 12, 3, 4, 4, 5, 14,
3, 6, 6, 6],
"t2": [4, 6, 12, 1, 5, 4, 2, 13, 9, 3, 4, 11, 4, 13, 1, 7, 1, 9, 10,
10, 3, 7, 10, 10, 6, 10, 10, 1, 2, 3, 4, 8, 6, 9, 8, 16, 22, 23,
23, 10, 50, 44, 46, 50],
"gamma2": [1]*44,
"nr3": [-0.31306260323435e2, 0.31546140237781e2, -0.25213154341695e4],
"d3": [3]*3,
"t3": [0, 1, 4],
"alfa3": [20]*3,
"beta3": [150, 150, 250],
"gamma3": [1.21, 1.21, 1.25],
"epsilon3": [1.]*3,
"nr4": [-0.14874640856724, 0.31806110878444],
"a4": [3.5, 3.5],
"b4": [0.85, 0.95],
"B": [0.2, 0.2],
"C": [28, 32],
"D": [700, 800],
"A": [0.32, .32],
"beta4": [0.3, 0.3]}
_Pv = {
"ao": [-7.85951783, 1.84408259, -11.7866497, 22.6807411, -15.9618719,
1.80122502],
"exp": [1, 1.5, 3, 3.5, 4, 7.5]}
_rhoL = {
"eq": 2,
"ao": [1.99274064, 1.09965342, -0.510839303, -1.75493479, -45.5170352,
-6.74694450e5],
"exp": [1, 2, 5, 16, 43, 110]}
_rhoG = {
"eq": 4,
"ao": [-2.0315024, -2.6830294, -5.38626492, -17.2991605, -44.7586581,
-63.9201063],
"exp": [1, 2, 4, 9, 18.5, 35.5]}
[docs]
def _phi0(self, tau, delta):
"""Low temperature extension of the IAPWS-95"""
prop = MEoS._phi0(self, tau, delta)
T = self.Tc/tau
if 50 <= T < 130:
fex, fext, fextt = self._phiex(T)
prop["fio"] += fex
prop["fiot"] += fext
prop["fiott"] += fextt
return prop
[docs]
def _phiex(self, T):
"""Low temperature extension"""
tau = self.Tc/T
E = 0.278296458178592
ep = self.Tc/130
fex = E*(-1/2/tau-3/ep**2*(tau+ep)*log(tau/ep)-9/2/ep+9*tau/2/ep**2
+ tau**2/2/ep**3)
fext = E*(1/2/tau**2-3/tau/ep-3/ep**2*log(tau/ep)+3/2/ep**2+tau/ep**3)
fextt = E*(-1/tau+1/ep)**3
return fex, fext, fextt
[docs]
@classmethod
def _alfa_sat(cls, T):
"""Auxiliary equation for the alfa coefficient for calculate the
enthalpy along the saturation line
Parameters
----------
T : float
Temperature, [K]
Returns
-------
alfa : float
alfa coefficient, [kJ/kg]
References
----------
IAPWS, Revised Supplementary Release on Saturation Properties of
Ordinary Water Substance September 1992,
http://www.iapws.org/relguide/Supp-sat.html, Eq.4
"""
di = [-1135.905627715, -5.65134998e-8, 2690.66631, 127.287297,
-135.003439, 0.981825814]
expi = [0, -19, 1, 4.5, 5, 54.5]
Tita = T/cls.Tc
alfa = 0
for d, x in zip(di, expi):
alfa += d*Tita**x
return alfa
[docs]
@classmethod
def _phi_sat(cls, T):
"""Auxiliary equation for the phi coefficient for calculate the
entropy along the saturation line
Parameters
----------
T : float
Temperature, [K]
Returns
-------
phi : float
phi coefficient, [kJ/kgK]
References
----------
IAPWS, Revised Supplementary Release on Saturation Properties of
Ordinary Water Substance September 1992,
http://www.iapws.org/relguide/Supp-sat.html, Eq.5
"""
di = [2319.5246, -5.65134998e-8*19/20, 2690.66631, 127.287297*9/7,
-135.003439*5/4, 0.981825814*109/107]
expi = [0, -20, None, 3.5, 4, 53.5]
Tita = T/cls.Tc
suma = 0
for d, x in zip(di, expi):
if x is None:
suma += d*log(Tita)
else:
suma += d*Tita**x
phi = suma/cls.Tc
return phi
[docs]
@classmethod
def _Liquid_Enthalpy(cls, T):
"""Auxiliary equation for the specific enthalpy for saturated liquid
Parameters
----------
T : float
Temperature, [K]
Returns
-------
h : float
Saturated liquid enthalpy, [kJ/kg]
References
----------
IAPWS, Revised Supplementary Release on Saturation Properties of
Ordinary Water Substance September 1992,
http://www.iapws.org/relguide/Supp-sat.html, Eq.6
"""
alfa = cls._alfa_sat(T)
rho = cls._Liquid_Density(T)
dpdT = cls._dPdT_sat(T)
h = alfa+T/rho*dpdT*1000
return h
[docs]
@classmethod
def _Vapor_Enthalpy(cls, T):
"""Auxiliary equation for the specific enthalpy for saturated vapor
Parameters
----------
T : float
Temperature, [K]
Returns
-------
h : float
Saturated vapor enthalpy, [kJ/kg]
References
----------
IAPWS, Revised Supplementary Release on Saturation Properties of
Ordinary Water Substance September 1992,
http://www.iapws.org/relguide/Supp-sat.html, Eq.7
"""
alfa = cls._alfa_sat(T)
rho = cls._Vapor_Density(T)
dpdT = cls._dPdT_sat(T)
h = alfa+T/rho*dpdT*1000
return h
[docs]
@classmethod
def _Liquid_Entropy(cls, T):
"""Auxiliary equation for the specific entropy for saturated liquid
Parameters
----------
T : float
Temperature, [K]
Returns
-------
s : float
Saturated liquid entropy, [kJ/kgK]
References
----------
IAPWS, Revised Supplementary Release on Saturation Properties of
Ordinary Water Substance September 1992,
http://www.iapws.org/relguide/Supp-sat.html, Eq.8
"""
phi = cls._phi_sat(T)
rho = cls._Liquid_Density(T)
dpdT = cls._dPdT_sat(T)
s = phi+dpdT/rho*1000
return s
[docs]
@classmethod
def _Vapor_Entropy(cls, T):
"""Auxiliary equation for the specific entropy for saturated vapor
Parameters
----------
T : float
Temperature, [K]
Returns
-------
s : float
Saturated liquid entropy, [kJ/kgK]
References
----------
IAPWS, Revised Supplementary Release on Saturation Properties of
Ordinary Water Substance September 1992,
http://www.iapws.org/relguide/Supp-sat.html, Eq.9
"""
phi = cls._phi_sat(T)
rho = cls._Vapor_Density(T)
dpdT = cls._dPdT_sat(T)
s = phi+dpdT/rho*1000
return s
[docs]
def _visco(self, rho, T, fase):
ref = IAPWS95()
st = ref._Helmholtz(rho, 1.5*Tc)
delta = rho/rhoc
drho = 1e3/self.R/1.5/Tc/(1+2*delta*st["fird"]+delta**2*st["firdd"])
return _Viscosity(rho, T, fase, drho)
[docs]
def _thermo(self, rho, T, fase):
ref = IAPWS95()
st = ref._Helmholtz(rho, 1.5*Tc)
delta = rho/rhoc
drho = 1e3/self.R/1.5/Tc/(1+2*delta*st["fird"]+delta**2*st["firdd"])
return _ThCond(rho, T, fase, drho)
[docs]
def _surface(self, T):
s = _Tension(T)
return s
[docs]
class IAPWS95_PT(IAPWS95):
"""Derivated class for direct P and T input"""
def __init__(self, P, T):
IAPWS95.__init__(self, T=T, P=P)
[docs]
class IAPWS95_Ph(IAPWS95):
"""Derivated class for direct P and h input"""
def __init__(self, P, h):
IAPWS95.__init__(self, P=P, h=h)
[docs]
class IAPWS95_Ps(IAPWS95):
"""Derivated class for direct P and s input"""
def __init__(self, P, s):
IAPWS95.__init__(self, P=P, s=s)
[docs]
class IAPWS95_Px(IAPWS95):
"""Derivated class for direct P and v input"""
def __init__(self, P, x):
IAPWS95.__init__(self, P=P, x=x)
[docs]
class IAPWS95_Tx(IAPWS95):
"""Derivated class for direct T and x input"""
def __init__(self, T, x):
IAPWS95.__init__(self, T=T, x=x)
[docs]
@mainClassDoc()
class D2O(MEoS):
"""Implementation of IAPWS Formulation for heavy water substance,
for internal procedures, see MEoS base class
Examples
--------
>>> hwater=D2O(T=300, rho=996.5560)
>>> hwater.P, hwater.Liquid.cv, hwater.Liquid.w
0.0030675947 4.21191157 5332.04871
References
----------
IAPWS, Release on the IAPWS Formulation 2017 for the Thermodynamic
Properties of Heavy Water, http://www.iapws.org/relguide/Heavy-2017.pdf
IAPWS, Revised Advisory Note No. 3: Thermodynamic Derivatives from IAPWS
Formulations, http://www.iapws.org/relguide/Advise3.pdf
"""
name = "heavy water"
CASNumber = "7789-20-0"
formula = "D2O"
synonym = "deuterium oxide"
Tc = Tc_D2O
rhoc = rhoc_D2O
Pc = Pc_D2O
M = 20.027508 # g/mol
Tt = 276.97
Tb = 374.563
f_acent = 0.364
momentoDipolar = 1.9
Fi0 = {"ao_log": [1, 3],
"pow": [0, 1],
"ao_pow": [-8.670994022646, 6.96033578458778],
"ao_exp": [0.010633, 0.99787, 2.1483, 0.3549],
"titao": [308/Tc, 1695/Tc, 3949/Tc, 10317/Tc],
"ao_hyp": [], "hyp": []}
_constants = {
"R": 8.3144598,
"nr1": [0.122082060e-1, 0.296956870e1, -0.379004540e1, 0.941089600,
-0.922466250, -0.139604190e-1],
"d1": [4, 1, 1, 2, 2, 3],
"t1": [1.0000, 0.6555, 0.9369, 0.5610, 0.7017, 1.0672],
"nr2": [-0.125203570, -0.555391500e1, -0.493009740e1, -0.359470240e-1,
-0.936172870e1, -0.691835150],
"c2": [1, 2, 2, 1, 2, 2],
"d2": [1, 1, 3, 2, 2, 1],
"t2": [3.9515, 4.6000, 5.1590, 0.2000, 5.4644, 2.3660],
"gamma2": [1]*6,
"nr3": [-0.456110600e-1, -0.224513300e1, 0.860006070e1, -0.248410420e1,
0.164476900e2, 0.270393360e1, 0.375637470e2, -0.177607760e1,
0.220924640e1, 0.519652000e1, 0.421097400, -0.391921100],
"t3": [3.4553, 1.4150, 1.5745, 3.4540, 3.8106, 4.8950, 1.4300, 1.5870,
3.7900, 2.6200, 1.9000, 4.3200],
"d3": [1, 3, 1, 3, 1, 1, 2, 2, 2, 1, 1, 1],
"alfa3": [0.6014, 1.4723, 1.5305, 2.4297, 1.3086, 1.3528, 3.4456,
1.2645, 2.5547, 1.2148, 18.738, 18.677],
"beta3": [0.4200, 2.4318, 1.2888, 8.2710, 0.3673, 0.9504, 7.8318,
3.3281, 7.1753, 0.9465, 1177.0, 1167.0],
"epsilon3": [1.8663, 0.2895, 0.5803, 0.2236, 0.6815, 0.9495, 1.1158,
0.1607, 0.4144, 0.9683, 0.9488, 0.9487],
"gamma3": [1.5414, 1.3794, 1.7385, 1.3045, 2.7242, 3.5321, 2.4552,
0.8319, 1.3500, 2.5617, 1.0491, 1.0486]}
_Pv = {
"ao": [-0.80236e1, 0.23957e1, -0.42639e2, 0.99569e2, -0.62135e2],
"exp": [1.0, 1.5, 2.75, 3.0, 3.2]}
_rhoL = {
"eq": 1,
"ao": [0.26406e1, 0.97090e1, -0.18058e2, 0.87202e1, -0.74487e1],
"exp": [0.3678, 1.9, 2.2, 2.63, 7.3]}
_rhoG = {
"eq": 3,
"ao": [-0.37651e1, -0.38673e2, 0.73024e2, -0.13251e3, 0.75235e2,
-0.70412e2],
"exp": [0.409, 1.766, 2.24, 3.04, 3.42, 6.9]}
[docs]
def _visco(self, rho, T, fase):
ref = D2O()
s = ref._Helmholtz(rho, 1.5*Tc_D2O)
delta = rho/rhoc_D2O
drho = 1e3/self.R/1.5/Tc_D2O/(1+2*delta*s["fird"]+delta**2*s["firdd"])
return _D2O_Viscosity(rho, T, fase, drho)
[docs]
def _thermo(self, rho, T, fase):
ref = D2O()
s = ref._Helmholtz(rho, 1.5*Tc_D2O)
delta = rho/rhoc_D2O
drho = 1e3/self.R/1.5/Tc_D2O/(1+2*delta*s["fird"]+delta**2*s["firdd"])
return _D2O_ThCond(rho, T, fase, drho)
[docs]
def _surface(self, T):
s = _D2O_Tension(T)
return s