#!/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
import os
import platform
import warnings
from numpy import exp, log, ndarray
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
[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))
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
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
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
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()
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"""
rho_c = self._constants.get("rhoref", self.rhoc)
T_c = self._constants.get("Tref", self.Tc)
if T > T_c:
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
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))
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
if T < 50:
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))
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()
if py_version[0] == "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