OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Improved H2 property calculations.

Browse Source 
Most H2 properties are calculated using Helmholtz free energy EOS.
Moved Helmholtz equations from Brine_H2 to H2 class.
New file with simple H2 property calculations based on ideal gas.
This commit is contained in:
Svenn Tveit
2022-03-14 17:09:12 +01:00
parent ef94fb8896
commit 066601ae52
3 changed files with 802 additions and 265 deletions
Download Patch File Download Diff File Expand all files Collapse all files

199
opm/material/binarycoefficients/Brine_H2.hpp
View File

@@ -28,7 +28,6 @@
#ifndef OPM_BINARY_COEFF_BRINE_H2_HPP
#define OPM_BINARY_COEFF_BRINE_H2_HPP
#include <opm/material/IdealGas.hpp>
#include <opm/material/binarycoefficients/FullerMethod.hpp>
namespace Opm {
@@ -40,7 +39,6 @@ namespace BinaryCoeff {
*/
template<class Scalar, class H2O, class H2, bool verbose = true>
class Brine_H2 {
using IdealGas = Opm::IdealGas<Scalar>;
static const int liquidPhaseIdx = 0; // index of the liquid phase
static const int gasPhaseIdx = 1; // index of the gas phase
@@ -149,7 +147,7 @@ public:
/*!
* \brief Calculate fugacity coefficient for H2 which is needed in calculation of H2 solubility in Li et al (2018).
* The equation used is based on Helmoltz free energy EOS. The formulas here are taken from Span et al., J. Phys.
* Chem. Ref. Data 29, 2000 and adapted to H2 in Li et al (2018).
* Chem. Ref. Data 29, 2000 and Leachman et al., J. Phys. Chem. Ref. Data 38, 2009, and Li et al. (2018).
*
* \param temperature temperature [K]
* \param pg gas phase pressure [Pa]
@@ -158,208 +156,21 @@ public:
static Evaluation fugacityCoefficientH2(const Evaluation& temperature, const Evaluation& pg)
{
// Convert pressure to reduced density and temperature to reduced temperature
Evaluation rho_red = convertPgToReducedRho_(temperature, pg);
Evaluation rho_red = H2::reducedMolarDensity(temperature, pg);
Evaluation T_red = H2::criticalTemperature() / temperature;
// Residual Helmholtz energy, Eq. (7) in Li et al. (2018)
Evaluation resHelm = residualHelmholtz_(T_red, rho_red);
Evaluation resHelm = H2::residualPartHelmholtz(T_red, rho_red);
// Derivative of residual Helmholtz energy wrt to reduced density, Eq. (73) in Span et al. (2018)
Evaluation dResdHelm = derivResidualHelmholtz_(T_red, rho_red);
Evaluation dResHelm_dRedRho = H2::derivResHelmholtzWrtRedRho(T_red, rho_red);
// Fugacity coefficient, Eq. (8) in Li et al. (2018)
Evaluation lnPhiH2 = resHelm + rho_red * dResdHelm - log(rho_red * dResdHelm + 1);
Evaluation lnPhiH2 = resHelm + rho_red * dResHelm_dRedRho - log(rho_red * dResHelm_dRedRho + 1);
return lnPhiH2;
}
/*!
* \brief Convert pressure to reduced density (rho/rho_crit) for further calculation of fugacity coefficient in Li et
* al. (2018) and Span et al. (2000). The conversion is done using the simplest root-finding algorithm, i.e. the
* bisection method.
*
* \param pg gas phase pressure [Pa]
* \param temperature temperature [K]
*/
template <class Evaluation>
static Evaluation convertPgToReducedRho_(const Evaluation& temperature, const Evaluation& pg)
{
// Interval for search. xmin = 0.0 is good since it guaranties negative fmin.
Evaluation rho_red_min = 0.0;
Evaluation fmin = -pg / 1.0e6; // at 0.0 we don't need to envoke function (see eq.(9) in Li et al. (2018))
// xmax must be high enough to give positive xmax. Check fmax for xmax given. If fmax is not positive we
// increase xmax
Evaluation rho_red_max = 2.0; // start
Evaluation fmax = rootFindingObj_(rho_red_max, temperature, pg);
if (Opm::getValue(fmax) <= 0.0) {
// increase 10 times, and if that does not give fmax > 0, we abort
for (int i=0; i < 10; ++i) {
rho_red_max *= 1.25;
fmax = rootFindingObj_(rho_red_max, temperature, pg);
if (Opm::getValue(fmax) < 0.0) {
break;
}
}
throw std::runtime_error("No max reduced density could be found for current pressure for bisection!");
}
// Bisection loop
for (int iteration=1; iteration<100; ++iteration) {
// New midpoint and its obj. value
Evaluation rho_red = (rho_red_min + rho_red_max) / 2;
Evaluation fmid = rootFindingObj_(rho_red, temperature, pg);
// Check if midpoint fulfills f=0 or x-xmin is sufficiently small
if (Opm::abs(fmid) < 1e-8 || Opm::abs((rho_red_max - rho_red_min) / 2) < 1e-8) {
return rho_red;
}
// Else we repeat with midpoint being either xmin or xmax (depending on the signs)
else if ((Opm::getValue(fmid) > 0.0 && Opm::getValue(fmin) < 0.0) ||
(Opm::getValue(fmid) < 0.0 && Opm::getValue(fmin) > 0.0)) {
// fmid has same sign as fmax so we set xmid as the new xmax
rho_red_max = rho_red;
}
else {
// fmid has same sign as fmin so we set xmid as the new xmin
rho_red_min = rho_red;
fmin = fmid;
}
}
throw std::runtime_error("No reduced density could be found for current pressure using bisection!");
}
/*!
* \brief Objective function in root-finding done in convertPgToReducedRho_ taken from Li et al. (2018).
*
* \param rho_red reduced density [-]
* \param pg gas phase pressure [Pa]
* \param temperature temperature [K]
*/
template <class Evaluation>
static Evaluation rootFindingObj_(const Evaluation& rho_red, const Evaluation& temperature, const Evaluation& pg)
{
// Temporary calculations
Evaluation T_red = H2::criticalTemperature() / temperature; // reciprocal reduced temp.
Evaluation p_MPa = pg / 1.0e6; // Pa --> MPa
Scalar R = IdealGas::R;
Evaluation rho_cRT = H2::criticalDensity() * R * temperature;
// Eq. (9)
Evaluation dResdH = derivResidualHelmholtz_(T_red, rho_red);
Evaluation obj = rho_red * rho_cRT * (1 + rho_red * dResdH) - p_MPa;
return obj;
}
/*!
* \brief Derivative of the residual part of Helmholtz energy wrt. reduced density. Used primarily to calculate
* fugacity coefficient for H2.
*
* \param T_red reduced temperature [-]
* \param rho_red reduced density [-]
*/
template <class Evaluation>
static Evaluation derivResidualHelmholtz_(const Evaluation& T_red, const Evaluation& rho_red)
{
// Various parameter values needed in calculations (Table 1 in Li et al. (2018))
static const Scalar N[14] = {-6.93643, 0.01, 2.1101, 4.52059, 0.732564, -1.34086, 0.130985, -0.777414,
0.351944, -0.0211716, 0.0226312, 0.032187, -0.0231752, 0.0557346};
static const Scalar t[14] = {0.6844, 1.0, 0.989, 0.489, 0.803, 1.1444, 1.409, 1.754, 1.311, 4.187, 5.646,
0.791, 7.249, 2.986};
static const int d[14] = {1, 4, 1, 1, 2, 2, 3, 1, 3, 2, 1, 3, 1, 1};
static const int p[2] = {1, 1};
static const Scalar phi[5] = {-1.685, -0.489, -0.103, -2.506, -1.607};
static const Scalar beta[5] = {-0.1710, -0.2245, -0.1304, -0.2785, -0.3967};
static const Scalar gamma[5] = {0.7164, 1.3444, 1.4517, 0.7204, 1.5445};
static const Scalar D[5] = {1.506, 0.156, 1.736, 0.670, 1.662};
// Derivative of Eq. (7) in Li et al. (2018), which can be compared with Eq. (73) in Span et al. (2000)
// First sum term
Evaluation s1 = 0.0;
for (int i = 0; i < 7; ++i) {
s1 += d[i] * N[i] * pow(rho_red, d[i]-1) * pow(T_red, t[i]);
}
// Second sum term
Evaluation s2 = 0.0;
for (int i = 7; i < 9; ++i) {
s2 += N[i] * pow(T_red, t[i]) * pow(rho_red, d[i]-1) * exp(-pow(rho_red, p[i-7])) *
(d[i] - p[i-7]*pow(rho_red, p[i-7]));
}
// Third, and last, sum term
Evaluation s3 = 0.0;
for (int i = 9; i < 15; ++i) {
s3 += N[i] * pow(T_red, t[i]) * pow(rho_red, d[i]-1) *
exp(phi[i-9] * pow(rho_red - D[i-9], 2) + beta[i-9] * pow(T_red - gamma[i-9], 2)) *
(d[i] + 2 * phi[i-9] * rho_red * (rho_red - D[i-9]));
}
// Return total sum
Evaluation s = s1 + s2 + s3;
return s;
}
/*!
* \brief The residual part of Helmholtz energy. Used primarily to calculate fugacity coefficient in Li et al (2018).
*
* \param T_red reduced temperature [-]
* \param rho_red reduced density [-]
*/
template <class Evaluation>
static Evaluation residualHelmholtz_(const Evaluation& T_red, const Evaluation& rho_red)
{
// Various parameter values needed in calculations (Table 1 in Li et al. (2018))
static const Scalar N[14] = {-6.93643, 0.01, 2.1101, 4.52059, 0.732564, -1.34086, 0.130985, -0.777414,
0.351944, -0.0211716, 0.0226312, 0.032187, -0.0231752, 0.0557346};
static const Scalar t[14] = {0.6844, 1.0, 0.989, 0.489, 0.803, 1.1444, 1.409, 1.754, 1.311, 4.187, 5.646,
0.791, 7.249, 2.986};
static const int d[14] = {1, 4, 1, 1, 2, 2, 3, 1, 3, 2, 1, 3, 1, 1};
static const int p[2] = {1, 1};
static const Scalar phi[5] = {-1.685, -0.489, -0.103, -2.506, -1.607};
static const Scalar beta[5] = {-0.1710, -0.2245, -0.1304, -0.2785, -0.3967};
static const Scalar gamma[5] = {0.7164, 1.3444, 1.4517, 0.7204, 1.5445};
static const Scalar D[5] = {1.506, 0.156, 1.736, 0.670, 1.662};
// Eq. (7) in Li et al. (2018), which can be compared with Eq. (55) in Span et al. (2000)
// First sum term
Evaluation s1 = 0.0;
for (int i = 0; i < 7; ++i) {
s1 += N[i] * pow(rho_red, d[i]) * pow(T_red, t[i]);
}
// Second sum term
Evaluation s2 = 0.0;
for (int i = 7; i < 9; ++i) {
s2 += N[i] * pow(T_red, t[i]) * pow(rho_red, d[i]) * exp(-pow(rho_red, p[i-7]));
}
// Third, and last, sum term
Evaluation s3 = 0.0;
for (int i = 9; i < 15; ++i) {
s3 += N[i] * pow(T_red, t[i]) * pow(rho_red, d[i]) *
exp(phi[i-9] * pow(rho_red - D[i-9], 2) + beta[i-9] * pow(T_red - gamma[i-9], 2));
}
// Return total sum
Evaluation s = s1 + s2 + s3;
return s;
}
/*!
* \brief Binary diffusion coefficent [m^2/s] for molecular water and H2 as an approximation for brine-H2 diffusion.
*

603
opm/material/components/H2.hpp
View File

@@ -21,10 +21,13 @@
copyright holders.
*/
/*!
* \file
* \copydoc Opm:H2
* \brief Properties of pure molecular hydrogen \f$H_2\f$.
*/
* \file
*
* \ingroup Components
*
* \copydoc Opm::H2
*
*/
#ifndef OPM_H2_HPP
#define OPM_H2_HPP
@@ -39,8 +42,15 @@ namespace Opm {
/*!
* \ingroup Components
* \brief Properties of pure molecular hydrogen \f$H_2\f$.
*
* \brief Properties of pure molecular hydrogen \f$H_2\f$. Most of the properties are calculated based on Leachman JW,
* Jacobsen RT, Penoncello SG, Lemmon EW. Fundamental equations of state for parahydrogen, normal hydrogen, and
* orthohydrogen. J Phys Chem Reference Data 2009;38:721e48. Their approach uses the fundamental Helmholtz free energy
* thermodynamic EOS, which is better suited for many gases such as H2. See also Span R, Lemmon EW, Jacobsen RT, Wagner
* W, Yokozeki A. A Reference Equation of State for the Thermodynamic Properties of Nitrogen for Temperatures from
* 63.151 to 1000 K and Pressures to 2200 MPa for explicit equations derived from the fundamental EOS formula.
*
* OBS: All equation and table references here are taken from Leachman et al. (2009) unless otherwise stated!
*
* \tparam Scalar The type used for scalar values
*/
template <class Scalar>
@@ -50,40 +60,52 @@ class H2 : public Component<Scalar, H2<Scalar> >
public:
/*!
* \brief A human readable name for the \f$H_2\f$.
*/
* \brief A human readable name for the \f$H_2\f$.
*/
static std::string name()
{ return "H2"; }
/*!
* \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of molecular hydrogen.
*/
* \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of molecular hydrogen.
*/
static constexpr Scalar molarMass()
{ return 2.01588e-3; }
/*!
* \brief Returns the critical temperature \f$\mathrm{[K]}\f$ of molecular hydrogen.
*/
* \brief Returns the critical temperature \f$\mathrm{[K]}\f$ of molecular hydrogen.
*/
static Scalar criticalTemperature()
{ return 33.2; /* [K] */ }
{ return 33.145; /* [K] */ }
/*!
* \brief Returns the critical pressure \f$\mathrm{[Pa]}\f$ of molecular hydrogen.
*/
* \brief Returns the critical pressure \f$\mathrm{[Pa]}\f$ of molecular hydrogen.
*/
static Scalar criticalPressure()
{ return 13.0e5; /* [N/m^2] */ }
{ return 1.2964e6; /* [N/m^2] */ }
/*!
* \brief Returns the critical density \f$\mathrm{[mol/cm^3]}\f$ of molecular hydrogen.
*/
* \brief Returns the critical density \f$\mathrm{[mol/cm^3]}\f$ of molecular hydrogen.
*/
static Scalar criticalDensity()
{ return 15.508e-3; /* [mol/cm^3] */ }
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ at molecular hydrogen's triple point.
*/
* \brief Returns the temperature \f$\mathrm{[K]}\f$ at molecular hydrogen's triple point.
*/
static Scalar tripleTemperature()
{ return 14.0; /* [K] */ }
{ return 13.957; /* [K] */ }
/*!
* \brief Returns the pressure \f$\mathrm{[Pa]}\f$ of molecular hydrogen's triple point.
*/
static Scalar triplePressure()
{ return 0.00736e6; /* [N/m^2] */ }
/*!
* \brief Returns the density \f$\mathrm{[mol/cm^3]}\f$ of molecular hydrogen's triple point.
*/
static Scalar tripleDensity()
{ return 38.2e-3; /* [mol/cm^3] */ }
/*!
* \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of pure molecular hydrogen
@@ -91,11 +113,6 @@ public:
*
*\param temperature temperature of component in \f$\mathrm{[K]}\f$
*
* Taken from:
*
* See: R. Reid, et al. (1987, pp 208-209, 669) \cite reid1987
*
* \todo implement the Gomez-Thodos approach...
*/
template <class Evaluation>
static Evaluation vaporPressure(Evaluation temperature)
@@ -106,12 +123,22 @@ public:
return 0; // H2 is solid: We don't take sublimation into
// account
// antoine equation
const Scalar A = -7.76451;
const Scalar B = 1.45838;
const Scalar C = -2.77580;
// Intermediate calculations involving temperature
Evaluation sigma = 1 - temperature/criticalTemperature();
Evaluation T_recp = criticalTemperature() / temperature;
return 1e5 * exp(A - B/(temperature + C));
// Parameters for normal hydrogen in Table 8
static const Scalar N[4] = {-4.89789, 0.988558, 0.349689, 0.499356};
static const Scalar k[4] = {1.0, 1.5, 2.0, 2.85};
// Eq. (33)
Evaluation s = 0.0; // sum calculation
for (int i = 0; i < 4; ++i) {
s += N[i] * pow(sigma, k[i]);
}
Evaluation lnPsigmaPc = T_recp * s;
return exp(lnPsigmaPc) * criticalPressure();
}
/*!
@@ -123,8 +150,11 @@ public:
template <class Evaluation>
static Evaluation gasDensity(Evaluation temperature, Evaluation pressure)
{
// Assume an ideal gas
return IdealGas::density(Evaluation(molarMass()), temperature, pressure);
// Use Eq. (56) in Span et al. (2000) to get molar density [mol/m3]
Evaluation rho_molar = gasMolarDensity(temperature, pressure);
// return density in [kg/m3]
return rho_molar * molarMass();
}
/*!
@@ -135,7 +165,7 @@ public:
*/
template <class Evaluation>
static Evaluation gasMolarDensity(Evaluation temperature, Evaluation pressure)
{ return IdealGas::molarDensity(temperature, pressure); }
{ return reducedMolarDensity(temperature, pressure) * criticalDensity() * 1e6; }
/*!
* \brief Returns true if the gas phase is assumed to be compressible
@@ -147,7 +177,7 @@ public:
* \brief Returns true if the gas phase is assumed to be ideal
*/
static constexpr bool gasIsIdeal()
{ return true; }
{ return false; }
/*!
* \brief The pressure of gaseous \f$H_2\f$ in \f$\mathrm{[Pa]}\f$ at a given density and temperature.
@@ -158,8 +188,12 @@ public:
template <class Evaluation>
static Evaluation gasPressure(Evaluation temperature, Evaluation density)
{
// Assume an ideal gas
return IdealGas::pressure(temperature, density/molarMass());
// Eq. (56) in Span et al. (2000)
Scalar R = IdealGas::R;
Evaluation rho_red = density / (molarMass() * criticalDensity() * 1e6);
Evaluation T_red = H2::criticalTemperature() / temperature;
return rho_red * criticalDensity() * R * temperature
* (1 + rho_red * derivResHelmholtzWrtRedRho(T_red, rho_red));
}
/*!
@@ -169,10 +203,30 @@ public:
static Evaluation gasInternalEnergy(const Evaluation& temperature,
const Evaluation& pressure)
{
const Evaluation& h = gasEnthalpy(temperature, pressure);
const Evaluation& rho = gasDensity(temperature, pressure);
// Eq. (58) in Span et al. (2000)
Evaluation rho_red = reducedMolarDensity(temperature, pressure);
Evaluation T_red = criticalTemperature() / temperature;
Scalar R = IdealGas::R;
return h - (pressure / rho);
return R * criticalTemperature() * (derivIdealHelmholtzWrtRecipRedTemp(T_red)
+ derivResHelmholtzWrtRecipRedTemp(T_red, rho_red)) / molarMass();
}
/*!
* \brief Specific enthalpy \f$\mathrm{[J/kg]}\f$ of pure hydrogen gas.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
template <class Evaluation>
static const Evaluation gasEnthalpy(Evaluation temperature,
Evaluation pressure)
{
// Eq. (59) in Span et al. (2000)
Evaluation u = gasInternalEnergy(temperature, pressure);
Evaluation rho = gasDensity(temperature, pressure);
return u + (pressure / rho);
}
/*!
@@ -214,45 +268,452 @@ public:
}
/*!
* \brief Specific enthalpy \f$\mathrm{[J/kg]}\f$ of pure hydrogen gas.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
* \brief Specific isobaric heat capacity \f$\mathrm{[J/(kg*K)]}\f$ of pure hydrogen gas. This is equivalent to the
* partial derivative of the specific enthalpy to the temperature.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*
*/
template <class Evaluation>
static const Evaluation gasEnthalpy(Evaluation temperature,
Evaluation pressure)
static const Evaluation gasHeatCapacity(Evaluation temperature,
Evaluation pressure)
{
return gasHeatCapacity(temperature, pressure) * temperature;
// Reduced variables
Evaluation rho_red = reducedMolarDensity(temperature, pressure);
Evaluation T_red = criticalTemperature() / temperature;
// Need Eq. (62) in Span et al. (2000)
Evaluation cv = gasIsochoricHeatCapacity(temperature, pressure); // [J/(kg*K)]
// Some intermediate calculations
Evaluation numerator = pow(1 + rho_red * derivResHelmholtzWrtRedRho(T_red, rho_red)
- rho_red * T_red * secDerivResHelmholtzWrtRecipRedTempAndRedRho(T_red, rho_red), 2);
Evaluation denominator = 1 + 2 * rho_red * derivResHelmholtzWrtRedRho(T_red, rho_red)
+ pow(rho_red, 2) * secDerivResHelmholtzWrtRedRho(T_red, rho_red);
// Eq. (63) in Span et al. (2000).
Scalar R = IdealGas::R;
Evaluation cp = cv + R * (numerator / denominator) / molarMass(); // divide by M to get [J/(kg*K)]
// Return
return cp;
}
/*!
* \brief Specific isobaric heat capacity \f$\mathrm{[J/(kg*K)]}\f$ of pure
* hydrogen gas.
*
* This is equivalent to the partial derivative of the specific
* enthalpy to the temperature.
* \param T temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*
* See: R. Reid, et al. (1987, pp 154, 657, 665) \cite reid1987
*/
* \brief Specific isochoric heat capacity \f$\mathrm{[J/(kg*K)]}\f$ of pure hydrogen gas.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*
*/
template <class Evaluation>
static const Evaluation gasHeatCapacity(Evaluation T,
Evaluation pressure)
static const Evaluation gasIsochoricHeatCapacity(Evaluation temperature,
Evaluation pressure)
{
// method of Joback
const Scalar cpVapA = 27.14;
const Scalar cpVapB = 9.273e-3;
const Scalar cpVapC = -1.381e-5;
const Scalar cpVapD = 7.645e-9;
// Reduced variables
Evaluation rho_red = reducedMolarDensity(temperature, pressure);
Evaluation T_red = criticalTemperature() / temperature;
return
1/molarMass()* // conversion from [J/(mol*K)] to [J/(kg*K)]
(cpVapA + T*
(cpVapB/2 + T*
(cpVapC/3 + T*
(cpVapD/4))));
// Eq. (62) in Span et al. (2000)
Scalar R = IdealGas::R;
Evaluation cv = R * (-pow(T_red, 2) * (secDerivIdealHelmholtzWrtRecipRedTemp(T_red)
+ secDerivResHelmholtzWrtRecipRedTemp(T_red, rho_red))); // [J/(mol*K)]
return cv / molarMass();
}
/*!
* \brief Calculate reduced density (rho/rho_crit) from pressure and temperature. The conversion is done using the
* simplest root-finding algorithm, i.e. the bisection method.
*
* \param pg gas phase pressure [Pa]
* \param temperature temperature [K]
*/
template <class Evaluation>
static Evaluation reducedMolarDensity(const Evaluation& temperature, const Evaluation& pg)
{
// Interval for search. xmin = 0.0 is good since it guaranties negative fmin.
Evaluation rho_red_min = 0.0;
Evaluation fmin = -pg / 1.0e6; // at 0.0 we don't need to envoke function (see eq.(9) in Li et al. (2018))
// xmax must be high enough to give positive fmax. Check fmax for a starting xmax. If fmax is not positive we
// increase xmax until we have fmax >= 0 and abort if we don't achieve that.
Evaluation rho_red_max = 2.0; // start
Evaluation fmax = rootFindingObj_(rho_red_max, temperature, pg);
if (Opm::getValue(fmax) <= 0.0) {
// increase 10 times, and if that does not give fmax > 0, we abort
for (int i=0; i < 10; ++i) {
rho_red_max *= 1.25;
fmax = rootFindingObj_(rho_red_max, temperature, pg);
if (Opm::getValue(fmax) > 0.0) {
break;
}
else if (Opm::getValue(fmax) < 0.0 && i == 9) {
throw std::runtime_error("No max reduced density could be found for current pressure!");
}
else if (Opm::getValue(fmax) == 0.0) {
return rho_red_max;
}
}
}
// Bisection loop
for (int iteration=1; iteration<100; ++iteration) {
// New midpoint and its obj. value
Evaluation rho_red = (rho_red_min + rho_red_max) / 2;
Evaluation fmid = rootFindingObj_(rho_red, temperature, pg);
// Check if midpoint fulfills f=0 or interval is sufficiently small
if (Opm::abs(fmid) < 1e-6 || Opm::abs((rho_red_max - rho_red_min) / 2) < 1e-6) {
return rho_red;
}
// Else we repeat with midpoint being either xmin or xmax (depending on the signs)
else if ((Opm::getValue(fmid) > 0.0 && Opm::getValue(fmin) < 0.0) ||
(Opm::getValue(fmid) < 0.0 && Opm::getValue(fmin) > 0.0)) {
// fmid has same sign as fmax so we set xmid as the new xmax
rho_red_max = rho_red;
}
else {
// fmid has same sign as fmin so we set xmid as the new xmin
rho_red_min = rho_red;
fmin = fmid;
}
}
throw std::runtime_error("No reduced density could be found for current pressure using bisection!");
}
/*!
* \brief The ideal-gas part of Helmholtz energy.
*
* \param T_red reduced temperature [-]
* \param rho_red reduced density [-]
*/
template <class Evaluation>
static Evaluation idealGasPartHelmholtz(const Evaluation& T_red, const Evaluation& rho_red)
{
// Eq. (31), which can be compared with Eq. (53) in Span et al. (2000)
// Terms not in sum
Evaluation s1 = log(rho_red) + 1.5*log(T_red) + a_[0] + a_[1] * T_red;
// Sum term
Evaluation s2 = 0.0;
for (int i = 2; i < 7; ++i) {
s1 += a_[i] * log(1 - exp(b_[i-2] * T_red));
}
// Return total
Evaluation s = s1 + s2;
return s;
}
/*!
* \brief Derivative of the ideal-gas part of Helmholtz energy wrt to reciprocal reduced temperature.
*
* \param T_red reduced temperature [-]
*/
template <class Evaluation>
static Evaluation derivIdealHelmholtzWrtRecipRedTemp(const Evaluation& T_red)
{
// Derivative of Eq. (31) wrt. reciprocal reduced temperature, which can be compared with Eq. (79) in Span et
// al. (2000)
// Terms not in sum
Evaluation s1 = (1.5 / T_red) + a_[1];
// Sum term
Evaluation s2 = 0.0;
for (int i = 2; i < 7; ++i) {
s2 += (-a_[i] * b_[i-2] * exp(b_[i-2] * T_red)) / (1 - exp(b_[i-2] * T_red));
}
// Return total
Evaluation s = s1 + s2;
return s;
}
/*!
* \brief Second derivative of the ideal-gas part of Helmholtz energy wrt to reciprocal reduced temperature.
*
* \param T_red reduced temperature [-]
* \param rho_red reduced density [-]
*/
template <class Evaluation>
static Evaluation secDerivIdealHelmholtzWrtRecipRedTemp(const Evaluation& T_red)
{
// Second derivative of Eq. (31) wrt. reciprocal reduced temperature, which can be compared with Eq. (80) in
// Span et al. (2000)
// Sum term
Evaluation s1 = 0.0;
for (int i = 2; i < 7; ++i) {
s1 += (-a_[i] * pow(b_[i-2], 2) * exp(b_[i-2] * T_red)) / pow(1 - exp(b_[i-2] * T_red), 2);
}
// Return total
Evaluation s = (-1.5 / pow(T_red, 2)) + s1;
return s;
}
/*!
* \brief The residual part of Helmholtz energy.
*
* \param T_red reduced temperature [-]
* \param rho_red reduced density [-]
*/
template <class Evaluation>
static Evaluation residualPartHelmholtz(const Evaluation& T_red, const Evaluation& rho_red)
{
// Eq. (32), which can be compared with Eq. (55) in Span et al. (2000)
// First sum term
Evaluation s1 = 0.0;
for (int i = 0; i < 7; ++i) {
s1 += N_[i] * pow(rho_red, d_[i]) * pow(T_red, t_[i]);
}
// Second sum term
Evaluation s2 = 0.0;
for (int i = 7; i < 9; ++i) {
s2 += N_[i] * pow(T_red, t_[i]) * pow(rho_red, d_[i]) * exp(-pow(rho_red, p_[i-7]));
}
// Third, and last, sum term
Evaluation s3 = 0.0;
for (int i = 9; i < 14; ++i) {
s3 += N_[i] * pow(T_red, t_[i]) * pow(rho_red, d_[i]) *
exp(phi_[i-9] * pow(rho_red - D_[i-9], 2) + beta_[i-9] * pow(T_red - gamma_[i-9], 2));
}
// Return total sum
Evaluation s = s1 + s2 + s3;
return s;
}
/*!
* \brief Derivative of the residual part of Helmholtz energy wrt. reduced density.
*
* \param T_red reduced temperature [-]
* \param rho_red reduced density [-]
*/
template <class Evaluation>
static Evaluation derivResHelmholtzWrtRedRho(const Evaluation& T_red, const Evaluation& rho_red)
{
// Derivative of Eq. (32) wrt to reduced density, which can be compared with Eq. (81) in Span et al. (2000)
// First sum term
Evaluation s1 = 0.0;
for (int i = 0; i < 7; ++i) {
s1 += d_[i] * N_[i] * pow(rho_red, d_[i]-1) * pow(T_red, t_[i]);
}
// Second sum term
Evaluation s2 = 0.0;
for (int i = 7; i < 9; ++i) {
s2 += N_[i] * pow(T_red, t_[i]) * pow(rho_red, d_[i]-1) * exp(-pow(rho_red, p_[i-7])) *
(d_[i] - p_[i-7]*pow(rho_red, p_[i-7]));
}
// Third, and last, sum term
Evaluation s3 = 0.0;
for (int i = 9; i < 14; ++i) {
s3 += N_[i] * pow(T_red, t_[i]) * pow(rho_red, d_[i]-1) *
exp(phi_[i-9] * pow(rho_red - D_[i-9], 2) + beta_[i-9] * pow(T_red - gamma_[i-9], 2)) *
(d_[i] + 2 * phi_[i-9] * rho_red * (rho_red - D_[i-9]));
}
// Return total sum
Evaluation s = s1 + s2 + s3;
return s;
}
/*!
* \brief Second derivative of the residual part of Helmholtz energy wrt. reduced density.
*
* \param T_red reduced temperature [-]
* \param rho_red reduced density [-]
*/
template <class Evaluation>
static Evaluation secDerivResHelmholtzWrtRedRho(const Evaluation& T_red, const Evaluation& rho_red)
{
// Second derivative of Eq. (32) wrt to reduced density, which can be compared with Eq. (82) in Span et al.
// (2000)
// First sum term
Evaluation s1 = 0.0;
for (int i = 0; i < 7; ++i) {
s1 += d_[i] * (d_[i] - 1) * N_[i] * pow(rho_red, d_[i]-2) * pow(T_red, t_[i]);
}
// Second sum term
Evaluation s2 = 0.0;
for (int i = 7; i < 9; ++i) {
s2 += N_[i] * pow(T_red, t_[i]) * pow(rho_red, d_[i]-2) * exp(-pow(rho_red, p_[i-7])) *
((d_[i] - p_[i-7] * pow(rho_red, p_[i-7])) * (d_[i] - p_[i-7] * pow(rho_red, p_[i-7]) - 1.0)
- pow(p_[i-7], 2) * pow(rho_red, p_[i-7]));
}
// Third, and last, sum term
Evaluation s3 = 0.0;
for (int i = 9; i < 14; ++i) {
s3 += N_[i] * pow(T_red, t_[i]) * pow(rho_red, d_[i]-2) *
exp(phi_[i-9] * pow(rho_red - D_[i-9], 2) + beta_[i-9] * pow(T_red - gamma_[i-9], 2)) *
(pow(d_[i] + 2 * phi_[i-9] * rho_red * (rho_red - D_[i-9]), 2)
- d_[i] + 2 * phi_[i-9] * pow(rho_red, 2));
}
// Return total sum
Evaluation s = s1 + s2 + s3;
return s;
}
/*!
* \brief Derivative of the residual part of Helmholtz energy wrt. reciprocal reduced temperature.
*
* \param T_red reduced temperature [-]
* \param rho_red reduced density [-]
*/
template <class Evaluation>
static Evaluation derivResHelmholtzWrtRecipRedTemp(const Evaluation& T_red, const Evaluation& rho_red)
{
// Derivative of Eq. (32) wrt to reciprocal reduced temperature, which can be compared with Eq. (84) in Span et
// al. (2000).
// First sum term
Evaluation s1 = 0.0;
for (int i = 0; i < 7; ++i) {
s1 += t_[i] * N_[i] * pow(rho_red, d_[i]) * pow(T_red, t_[i]-1);
}
// Second sum term
Evaluation s2 = 0.0;
for (int i = 7; i < 9; ++i) {
s2 += t_[i] * N_[i] * pow(T_red, t_[i]-1) * pow(rho_red, d_[i]) * exp(-pow(rho_red, p_[i-7]));
}
// Third, and last, sum term
Evaluation s3 = 0.0;
for (int i = 9; i < 14; ++i) {
s3 += N_[i] * pow(T_red, t_[i]-1) * pow(rho_red, d_[i]) *
exp(phi_[i-9] * pow(rho_red - D_[i-9], 2) + beta_[i-9] * pow(T_red - gamma_[i-9], 2)) *
(t_[i] + 2 * beta_[i-9] * T_red * (T_red - gamma_[i-9]));
}
// Return total sum
Evaluation s = s1 + s2 + s3;
return s;
}
/*!
* \brief Second derivative of the residual part of Helmholtz energy wrt. reciprocal reduced temperature.
*
* \param T_red reduced temperature [-]
* \param rho_red reduced density [-]
*/
template <class Evaluation>
static Evaluation secDerivResHelmholtzWrtRecipRedTemp(const Evaluation& T_red, const Evaluation& rho_red)
{
// Second derivative of Eq. (32) wrt to reciprocal reduced temperature, which can be compared with Eq. (85) in
// Span et al. (2000).
// First sum term
Evaluation s1 = 0.0;
for (int i = 0; i < 7; ++i) {
s1 += t_[i] * (t_[i] - 1) * N_[i] * pow(rho_red, d_[i]) * pow(T_red, t_[i]-2);
}
// Second sum term
Evaluation s2 = 0.0;
for (int i = 7; i < 9; ++i) {
s2 += t_[i] * (t_[i] - 1) * N_[i] * pow(T_red, t_[i]-2) * pow(rho_red, d_[i]) * exp(-pow(rho_red, p_[i-7]));
}
// Third, and last, sum term
Evaluation s3 = 0.0;
for (int i = 9; i < 14; ++i) {
s3 += N_[i] * pow(T_red, t_[i]-2) * pow(rho_red, d_[i]) *
exp(phi_[i-9] * pow(rho_red - D_[i-9], 2) + beta_[i-9] * pow(T_red - gamma_[i-9], 2)) *
(pow(t_[i] + 2 * beta_[i-9] * T_red * (T_red - gamma_[i-9]), 2)
- t_[i] + 2 * beta_[i-9] * pow(T_red, 2));
}
// Return total sum
Evaluation s = s1 + s2 + s3;
return s;
}
/*!
* \brief Second derivative of the residual part of Helmholtz energy first wrt. reciprocal reduced temperature, and
* second wrt. reduced density (i.e. d^2 H / drho_red dT_red).
*
* \param T_red reduced temperature [-]
* \param rho_red reduced density [-]
*/
template <class Evaluation>
static Evaluation secDerivResHelmholtzWrtRecipRedTempAndRedRho(const Evaluation& T_red, const Evaluation& rho_red)
{
// Second derivative of Eq. (32) wrt to reciprocal reduced temperature and reduced density, which can be
// compared with Eq. (86) in Span et al. (2000).
// First sum term
Evaluation s1 = 0.0;
for (int i = 0; i < 7; ++i) {
s1 += t_[i] * d_[i] * N_[i] * pow(rho_red, d_[i]-1) * pow(T_red, t_[i]-1);
}
// Second sum term
Evaluation s2 = 0.0;
for (int i = 7; i < 9; ++i) {
s2 += t_[i] * N_[i] * pow(T_red, t_[i]-1) * pow(rho_red, d_[i]-1) * exp(-pow(rho_red, p_[i-7]))
* (d_[i] - p_[i-7] * pow(rho_red, p_[i-7]));
}
// Third, and last, sum term
Evaluation s3 = 0.0;
for (int i = 9; i < 14; ++i) {
s3 += N_[i] * pow(T_red, t_[i]-1) * pow(rho_red, d_[i]-1) *
exp(phi_[i-9] * pow(rho_red - D_[i-9], 2) + beta_[i-9] * pow(T_red - gamma_[i-9], 2)) *
(t_[i] + 2 * beta_[i-9] * T_red * (T_red - gamma_[i-9]))
* (d_[i] + 2 * phi_[i-9] * rho_red * (rho_red - D_[i-9]));
}
// Return total sum
Evaluation s = s1 + s2 + s3;
return s;
}
private:
// Parameter values need in the ideal-gas contribution to the reduced Helmholtz free energy given in Table 4
static constexpr Scalar a_[7] = {-1.4579856475, 1.888076782, 1.616, -0.4117, -0.792, 0.758, 1.217};
static constexpr Scalar b_[5] = {-16.0205159149, -22.6580178006, -60.0090511389, -74.9434303817, -206.9392065168};
// Parameter values needed in the residual contribution to the reduced Helmholtz free energy given in Table 5.
static constexpr Scalar N_[14] = {-6.93643, 0.01, 2.1101, 4.52059, 0.732564, -1.34086, 0.130985, -0.777414,
0.351944, -0.0211716, 0.0226312, 0.032187, -0.0231752, 0.0557346};
static constexpr Scalar t_[14] = {0.6844, 1.0, 0.989, 0.489, 0.803, 1.1444, 1.409, 1.754, 1.311, 4.187, 5.646,
0.791, 7.249, 2.986};
static constexpr Scalar d_[14] = {1, 4, 1, 1, 2, 2, 3, 1, 3, 2, 1, 3, 1, 1};
static constexpr Scalar p_[2] = {1, 1};
static constexpr Scalar phi_[5] = {-1.685, -0.489, -0.103, -2.506, -1.607};
static constexpr Scalar beta_[5] = {-0.1710, -0.2245, -0.1304, -0.2785, -0.3967};
static constexpr Scalar gamma_[5] = {0.7164, 1.3444, 1.4517, 0.7204, 1.5445};
static constexpr Scalar D_[5] = {1.506, 0.156, 1.736, 0.670, 1.662};
/*!
* \brief Objective function in root-finding done in convertPgToReducedRho.
*
* \param rho_red reduced density [-]
* \param pg gas phase pressure [Pa]
* \param temperature temperature [K]
*/
template <class Evaluation>
static Evaluation rootFindingObj_(const Evaluation& rho_red, const Evaluation& temperature, const Evaluation& pg)
{
// Temporary calculations
Evaluation T_red = criticalTemperature() / temperature; // reciprocal reduced temp.
Evaluation p_MPa = pg / 1.0e6; // Pa --> MPa
Scalar R = IdealGas::R;
Evaluation rho_cRT = criticalDensity() * R * temperature;
// Eq. (56) in Span et al. (2000)
Evaluation dResHelm_dRedRho = derivResHelmholtzWrtRedRho(T_red, rho_red);
Evaluation obj = rho_red * rho_cRT * (1 + rho_red * dResHelm_dRedRho) - p_MPa;
return obj;
}
};

265
opm/material/components/SimpleH2.hpp Normal file
View File

@@ -0,0 +1,265 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \ingroup Components
*
* \copydoc Opm::SimpleH2
*
*/
#ifndef OPM_SIMPLE_H2_HPP
#define OPM_SIMPLE_H2_HPP
#include <opm/material/IdealGas.hpp>
#include <opm/material/components/Component.hpp>
#include <opm/material/densead/Math.hpp>
#include <cmath>
namespace Opm {
/*!
* \ingroup Components
*
* \brief Properties of pure molecular hydrogen \f$H_2\f$. Use ideal gas equations for many properties.
*
* \tparam Scalar The type used for scalar values
*/
template <class Scalar>
class SimpleH2 : public Component<Scalar, SimpleH2<Scalar> >
{
using IdealGas = Opm::IdealGas<Scalar>;
public:
/*!
* \brief A human readable name for the \f$H_2\f$.
*/
static std::string name()
{ return "H2"; }
/*!
* \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of molecular hydrogen.
*/
static constexpr Scalar molarMass()
{ return 2.01588e-3; }
/*!
* \brief Returns the critical temperature \f$\mathrm{[K]}\f$ of molecular hydrogen.
*/
static Scalar criticalTemperature()
{ return 33.2; /* [K] */ }
/*!
* \brief Returns the critical pressure \f$\mathrm{[Pa]}\f$ of molecular hydrogen.
*/
static Scalar criticalPressure()
{ return 13.0e5; /* [N/m^2] */ }
/*!
* \brief Returns the critical density \f$\mathrm{[mol/cm^3]}\f$ of molecular hydrogen.
*/
static Scalar criticalDensity()
{ return 15.508e-3; /* [mol/cm^3] */ }
/*!
* \brief Returns the temperature \f$\mathrm{[K]}\f$ at molecular hydrogen's triple point.
*/
static Scalar tripleTemperature()
{ return 14.0; /* [K] */ }
/*!
* \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of pure molecular hydrogen
* at a given temperature.
*
*\param temperature temperature of component in \f$\mathrm{[K]}\f$
*
* Taken from:
*
* See: R. Reid, et al. (1987, pp 208-209, 669) \cite reid1987
*
* \todo implement the Gomez-Thodos approach...
*/
template <class Evaluation>
static Evaluation vaporPressure(Evaluation temperature)
{
if (temperature > criticalTemperature())
return criticalPressure();
if (temperature < tripleTemperature())
return 0; // H2 is solid: We don't take sublimation into
// account
// antoine equation
const Scalar A = -7.76451;
const Scalar B = 1.45838;
const Scalar C = -2.77580;
return 1e5 * exp(A - B/(temperature + C));
}
/*!
* \brief The density \f$\mathrm{[kg/m^3]}\f$ of \f$H_2\f$ at a given pressure and temperature.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
template <class Evaluation>
static Evaluation gasDensity(Evaluation temperature, Evaluation pressure)
{
// Assume an ideal gas
return IdealGas::density(Evaluation(molarMass()), temperature, pressure);
}
/*!
* \brief The molar density of \f$H_2\f$ in \f$\mathrm{[mol/m^3]}\f$,
* depending on pressure and temperature.
* \param temperature The temperature of the gas
* \param pressure The pressure of the gas
*/
template <class Evaluation>
static Evaluation gasMolarDensity(Evaluation temperature, Evaluation pressure)
{ return IdealGas::molarDensity(temperature, pressure); }
/*!
* \brief Returns true if the gas phase is assumed to be compressible
*/
static constexpr bool gasIsCompressible()
{ return true; }
/*!
* \brief Returns true if the gas phase is assumed to be ideal
*/
static constexpr bool gasIsIdeal()
{ return true; }
/*!
* \brief The pressure of gaseous \f$H_2\f$ in \f$\mathrm{[Pa]}\f$ at a given density and temperature.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param density density of component in \f$\mathrm{[kg/m^3]}\f$
*/
template <class Evaluation>
static Evaluation gasPressure(Evaluation temperature, Evaluation density)
{
// Assume an ideal gas
return IdealGas::pressure(temperature, density/molarMass());
}
/*!
* \brief Specific internal energy of H2 [J/kg].
*/
template <class Evaluation>
static Evaluation gasInternalEnergy(const Evaluation& temperature,
const Evaluation& pressure)
{
const Evaluation& h = gasEnthalpy(temperature, pressure);
const Evaluation& rho = gasDensity(temperature, pressure);
return h - (pressure / rho);
}
/*!
* \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of \f$H_2\f$ at a given pressure and temperature.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*
* See:
*
* See: R. Reid, et al.: The Properties of Gases and Liquids,
* 4th edition, McGraw-Hill, 1987, pp 396-397,
* 5th edition, McGraw-Hill, 2001 pp 9.7-9.8 (omega and V_c taken from p. A.19)
*
*/
template <class Evaluation>
static Evaluation gasViscosity(const Evaluation& temperature, const Evaluation& /*pressure*/)
{
const Scalar Tc = criticalTemperature();
const Scalar Vc = 64.2; // critical specific volume [cm^3/mol]
const Scalar omega = -0.217; // accentric factor
const Scalar M = molarMass() * 1e3; // molar mas [g/mol]
const Scalar dipole = 0.0; // dipole moment [debye]
Scalar mu_r4 = 131.3 * dipole / std::sqrt(Vc * Tc);
mu_r4 *= mu_r4;
mu_r4 *= mu_r4;
Scalar Fc = 1 - 0.2756*omega + 0.059035*mu_r4;
const Evaluation& Tstar = 1.2593 * temperature/Tc;
const Evaluation& Omega_v =
1.16145*pow(Tstar, -0.14874) +
0.52487*exp(- 0.77320*Tstar) +
2.16178*exp(- 2.43787*Tstar);
const Evaluation& mu = 40.785*Fc*sqrt(M*temperature)/(std::pow(Vc, 2./3)*Omega_v);
// convertion from micro poise to Pa s
return mu/1e6 / 10;
}
/*!
* \brief Specific enthalpy \f$\mathrm{[J/kg]}\f$ of pure hydrogen gas.
*
* \param temperature temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*/
template <class Evaluation>
static const Evaluation gasEnthalpy(Evaluation temperature,
Evaluation pressure)
{
return gasHeatCapacity(temperature, pressure) * temperature;
}
/*!
* \brief Specific isobaric heat capacity \f$\mathrm{[J/(kg*K)]}\f$ of pure
* hydrogen gas.
*
* This is equivalent to the partial derivative of the specific
* enthalpy to the temperature.
* \param T temperature of component in \f$\mathrm{[K]}\f$
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
*
* See: R. Reid, et al. (1987, pp 154, 657, 665) \cite reid1987
*/
template <class Evaluation>
static const Evaluation gasHeatCapacity(Evaluation T,
Evaluation pressure)
{
// method of Joback
const Scalar cpVapA = 27.14;
const Scalar cpVapB = 9.273e-3;
const Scalar cpVapC = -1.381e-5;
const Scalar cpVapD = 7.645e-9;
return
1/molarMass()* // conversion from [J/(mol*K)] to [J/(kg*K)]
(cpVapA + T*
(cpVapB/2 + T*
(cpVapC/3 + T*
(cpVapD/4))));
}
};
} // end namespace Opm
#endif