make the binary coefficients local-AD aware
This commit is contained in:
Andreas Lauser
parent
b5b7864ca5
commit
0149901e4f
@ -40,8 +40,8 @@ public:
|
||||
/*!
|
||||
*
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar henry(Scalar temperature)
|
||||
template <class Evaluation>
|
||||
static Evaluation henry(const Evaluation& temperature)
|
||||
{ OPM_THROW(std::runtime_error, "Not implemented: Henry coefficient of air in mesitylene"); }
|
||||
|
||||
/*!
|
||||
@ -51,36 +51,37 @@ public:
|
||||
* W.J. Lyman, W.F. Reehl, D.H. Rosenblatt
|
||||
*
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
template <class Evaluation>
|
||||
static Evaluation gasDiffCoeff(Evaluation temperature, Evaluation pressure)
|
||||
{
|
||||
typedef Opm::Air<Scalar> Air;
|
||||
typedef Opm::Mesitylene<Scalar> Mesitylene;
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
typedef Opm::Air<double> Air;
|
||||
typedef Opm::Mesitylene<double> Mesitylene;
|
||||
|
||||
temperature = std::max(temperature, 1e-9); // regularization
|
||||
temperature = std::min(temperature, 500.0); // regularization
|
||||
pressure = std::max(pressure, 0.0); // regularization
|
||||
pressure = std::min(pressure, 1e8); // regularization
|
||||
temperature = Toolbox::max(temperature, 1e-9); // regularization
|
||||
temperature = Toolbox::min(temperature, 500.0); // regularization
|
||||
pressure = Toolbox::max(pressure, 0.0); // regularization
|
||||
pressure = Toolbox::min(pressure, 1e8); // regularization
|
||||
|
||||
const Scalar M_m = 1e3*Mesitylene::molarMass(); // [g/mol] molecular weight of mesitylene
|
||||
const Scalar M_a = 1e3*Air::molarMass(); // [g/mol] molecular weight of air
|
||||
const Scalar Tb_m = 437.9; // [K] boiling temperature of mesitylene
|
||||
const Scalar sigma_a = 3.711; // charact. length of air
|
||||
const Scalar T_scal_a = 78.6; // [K] (molec. energy of attraction/Boltzmann constant)
|
||||
const Scalar V_B_m = 162.6; // [cm^3/mol] LeBas molal volume of mesitylene
|
||||
const Scalar sigma_m = 1.18*std::pow(V_B_m, 0.333); // charact. length of mesitylene
|
||||
const Scalar sigma_am = 0.5*(sigma_a + sigma_m);
|
||||
const Scalar T_scal_m = 1.15*Tb_m;
|
||||
const Scalar T_scal_am = std::sqrt(T_scal_a*T_scal_m);
|
||||
const double M_m = 1e3*Mesitylene::molarMass(); // [g/mol] molecular weight of mesitylene
|
||||
const double M_a = 1e3*Air::molarMass(); // [g/mol] molecular weight of air
|
||||
const double Tb_m = 437.9; // [K] boiling temperature of mesitylene
|
||||
const double sigma_a = 3.711; // charact. length of air
|
||||
const double T_scal_a = 78.6; // [K] (molec. energy of attraction/Boltzmann constant)
|
||||
const double V_B_m = 162.6; // [cm^3/mol] LeBas molal volume of mesitylene
|
||||
const double sigma_m = 1.18*std::pow(V_B_m, 0.333); // charact. length of mesitylene
|
||||
const double sigma_am = 0.5*(sigma_a + sigma_m);
|
||||
const double T_scal_m = 1.15*Tb_m;
|
||||
const double T_scal_am = std::sqrt(T_scal_a*T_scal_m);
|
||||
|
||||
Scalar T_star = temperature/T_scal_am;
|
||||
T_star = std::max(T_star, 1e-5); // regularization
|
||||
Evaluation T_star = temperature/T_scal_am;
|
||||
T_star = Toolbox::max(T_star, 1e-5); // regularization
|
||||
|
||||
const Scalar Omega = 1.06036/std::pow(T_star, 0.1561) + 0.193/std::exp(T_star*0.47635)
|
||||
+ 1.03587/std::exp(T_star*1.52996) + 1.76474/std::exp(T_star*3.89411);
|
||||
const Scalar B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_a + 1.0/M_m);
|
||||
const Scalar Mr = (M_a + M_m)/(M_a*M_m);
|
||||
const Scalar D_am = (B_*std::pow(temperature, 1.5) * std::sqrt(Mr))
|
||||
const Evaluation Omega = 1.06036/Toolbox::pow(T_star, 0.1561) + 0.193/Toolbox::exp(T_star*0.47635)
|
||||
+ 1.03587/Toolbox::exp(T_star*1.52996) + 1.76474/Toolbox::exp(T_star*3.89411);
|
||||
const double B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_a + 1.0/M_m);
|
||||
const double Mr = (M_a + M_m)/(M_a*M_m);
|
||||
const Evaluation D_am = (B_*Toolbox::pow(temperature, 1.5) * std::sqrt(Mr))
|
||||
/(1e-5*pressure*std::pow(sigma_am, 2.0) * Omega); // [cm^2/s]
|
||||
|
||||
return 1e-4*D_am; // [m^2/s]
|
||||
@ -89,10 +90,9 @@ public:
|
||||
/*!
|
||||
* \brief Diffusion coefficent [m^2/s] for molecular mesitylene in liquid water.
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
{ OPM_THROW(std::runtime_error, "Not implemented: Binary liquid diffusion coefficients of air and mesitylene");
|
||||
}
|
||||
template <class Evaluation>
|
||||
static Evaluation liquidDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
|
||||
{ OPM_THROW(std::runtime_error, "Not implemented: Binary liquid diffusion coefficients of air and mesitylene"); }
|
||||
};
|
||||
|
||||
} // namespace BinaryCoeff
|
||||
|
||||
@ -26,10 +26,8 @@
|
||||
#include <opm/material/components/Air.hpp>
|
||||
#include <opm/material/components/Xylene.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
namespace BinaryCoeff
|
||||
{
|
||||
namespace Opm {
|
||||
namespace BinaryCoeff {
|
||||
|
||||
/*!
|
||||
* \brief Binary coefficients for water and xylene.
|
||||
@ -40,8 +38,8 @@ public:
|
||||
/*!
|
||||
*
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar henry(Scalar temperature)
|
||||
template <class Evaluation>
|
||||
static Evaluation henry(const Evaluation& temperature)
|
||||
{ OPM_THROW(std::runtime_error, "Not implemented: Henry coefficient of air in xylene"); }
|
||||
|
||||
/*!
|
||||
@ -51,39 +49,38 @@ public:
|
||||
* W.J. Lyman, W.F. Reehl, D.H. Rosenblatt
|
||||
*
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
template <class Evaluation>
|
||||
static Evaluation gasDiffCoeff(Evaluation temperature, Evaluation pressure)
|
||||
{
|
||||
typedef Opm::Air<Scalar> Air;
|
||||
typedef Opm::Xylene<Scalar> Xylene;
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
typedef Opm::Air<double> Air;
|
||||
typedef Opm::Xylene<double> Xylene;
|
||||
|
||||
temperature = std::max(temperature, 1e-9); // regularization
|
||||
temperature = std::min(temperature, 500.0); // regularization
|
||||
pressure = std::max(pressure, 0.0); // regularization
|
||||
pressure = std::min(pressure, 1e8); // regularization
|
||||
temperature = Toolbox::max(temperature, 1e-9); // regularization
|
||||
temperature = Toolbox::min(temperature, 500.0); // regularization
|
||||
pressure = Toolbox::max(pressure, 0.0); // regularization
|
||||
pressure = Toolbox::min(pressure, 1e8); // regularization
|
||||
|
||||
const Scalar M_x = 1e3*Xylene::molarMass(); // [g/mol] molecular weight of xylene
|
||||
const Scalar M_a = 1e3*Air::molarMass(); // [g/mol] molecular weight of air
|
||||
const Scalar Tb_x = 412.0; // [K] boiling temperature of xylene
|
||||
const Scalar sigma_a = 3.711; // charact. length of air
|
||||
const Scalar T_scal_a = 78.6; // [K] (molec. energy of attraction/Boltzmann constant)
|
||||
const Scalar V_B_x = 140.4; // [cm^3/mol] LeBas molal volume of xylene
|
||||
const Scalar sigma_x = 1.18*std::pow(V_B_x, 0.333); // charact. length of xylene
|
||||
const Scalar sigma_ax = 0.5*(sigma_a + sigma_x);
|
||||
const Scalar T_scal_x = 1.15*Tb_x;
|
||||
const Scalar T_scal_ax = std::sqrt(T_scal_a*T_scal_x);
|
||||
const double M_x = 1e3*Xylene::molarMass(); // [g/mol] molecular weight of xylene
|
||||
const double M_a = 1e3*Air::molarMass(); // [g/mol] molecular weight of air
|
||||
const double Tb_x = 412.0; // [K] boiling temperature of xylene
|
||||
const double sigma_a = 3.711; // charact. length of air
|
||||
const double T_scal_a = 78.6; // [K] (molec. energy of attraction/Boltzmann constant)
|
||||
const double V_B_x = 140.4; // [cm^3/mol] LeBas molal volume of xylene
|
||||
const double sigma_x = 1.18*std::pow(V_B_x, 0.333); // charact. length of xylene
|
||||
const double sigma_ax = 0.5*(sigma_a + sigma_x);
|
||||
const double T_scal_x = 1.15*Tb_x;
|
||||
const double T_scal_ax = std::sqrt(T_scal_a*T_scal_x);
|
||||
|
||||
Scalar T_star = temperature/T_scal_ax;
|
||||
T_star = std::max(T_star, 1e-5); // regularization
|
||||
const Evaluation& T_star = Toolbox::max(temperature/T_scal_ax, 1e-5);
|
||||
|
||||
const Scalar Omega = 1.06036/std::pow(T_star, 0.1561) + 0.193/std::exp(T_star*0.47635)
|
||||
+ 1.03587/std::exp(T_star*1.52996) + 1.76474/std::exp(T_star*3.89411);
|
||||
const Scalar B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_a + 1.0/M_x);
|
||||
const Scalar Mr = (M_a + M_x)/(M_a*M_x);
|
||||
const Scalar D_ax = (B_*std::pow(temperature,1.5)*std::sqrt(Mr))
|
||||
/(1e-5*pressure*std::pow(sigma_ax, 2.0)*Omega); // [cm^2/s]
|
||||
|
||||
return D_ax*1e-4; // [m^2/s]
|
||||
const Evaluation& Omega = 1.06036/Toolbox::pow(T_star, 0.1561) + 0.193/Toolbox::exp(T_star*0.47635)
|
||||
+ 1.03587/Toolbox::exp(T_star*1.52996) + 1.76474/Toolbox::exp(T_star*3.89411);
|
||||
const double B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_a + 1.0/M_x);
|
||||
const double Mr = (M_a + M_x)/(M_a*M_x);
|
||||
return 1e-4
|
||||
*(B_*Toolbox::pow(temperature,1.5)*std::sqrt(Mr))
|
||||
/(1e-5*pressure*std::pow(sigma_ax, 2.0)*Omega);
|
||||
}
|
||||
|
||||
/*!
|
||||
@ -91,10 +88,9 @@ public:
|
||||
*
|
||||
* \todo
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
{ OPM_THROW(std::runtime_error, "Not implemented: Binary liquid diffusion coefficients of air and xylene");
|
||||
}
|
||||
template <class Evaluation>
|
||||
static Evaluation liquidDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
|
||||
{ OPM_THROW(std::runtime_error, "Not implemented: Binary liquid diffusion coefficients of air and xylene"); }
|
||||
};
|
||||
|
||||
} // namespace BinaryCoeff
|
||||
|
||||
@ -52,16 +52,15 @@ public:
|
||||
* \param temperature the temperature [K]
|
||||
* \param pressure the phase pressure [Pa]
|
||||
*/
|
||||
static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
template <class Evaluation>
|
||||
static Evaluation gasDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
|
||||
{
|
||||
//Diffusion coefficient of water in the CO2 phase
|
||||
Scalar const PI=3.141593;
|
||||
Scalar const k = 1.3806504e-23; // Boltzmann constant
|
||||
Scalar const c = 4; // slip parameter, can vary between 4 (slip condition) and 6 (stick condition)
|
||||
Scalar const R_h = 1.72e-10; // hydrodynamic radius of the solute
|
||||
Scalar mu = CO2::gasViscosity(temperature, pressure); // CO2 viscosity
|
||||
Scalar D = k / (c * PI * R_h) * (temperature / mu);
|
||||
return D;
|
||||
Scalar k = 1.3806504e-23; // Boltzmann constant
|
||||
Scalar c = 4; // slip parameter, can vary between 4 (slip condition) and 6 (stick condition)
|
||||
Scalar R_h = 1.72e-10; // hydrodynamic radius of the solute
|
||||
const Evaluation& mu = CO2::gasViscosity(temperature, pressure); // CO2 viscosity
|
||||
return k / (c * M_PI * R_h) * (temperature / mu);
|
||||
}
|
||||
|
||||
/*!
|
||||
@ -70,10 +69,13 @@ public:
|
||||
* \param temperature the temperature [K]
|
||||
* \param pressure the phase pressure [Pa]
|
||||
*/
|
||||
static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
template <class Evaluation>
|
||||
static Evaluation liquidDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
|
||||
{
|
||||
typedef MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
//Diffusion coefficient of CO2 in the brine phase
|
||||
return 2e-9;
|
||||
return Toolbox::createConstant(2e-9);
|
||||
}
|
||||
|
||||
/*!
|
||||
@ -93,25 +95,26 @@ public:
|
||||
* \param xlCO2 mole fraction of CO2 in brine [mol/mol]
|
||||
* \param ygH2O mole fraction of water in the gas phase [mol/mol]
|
||||
*/
|
||||
static void calculateMoleFractions(const Scalar temperature,
|
||||
const Scalar pg,
|
||||
const Scalar salinity,
|
||||
template <class Evaluation>
|
||||
static void calculateMoleFractions(const Evaluation& temperature,
|
||||
const Evaluation& pg,
|
||||
Scalar salinity,
|
||||
const int knownPhaseIdx,
|
||||
Scalar &xlCO2,
|
||||
Scalar &ygH2O)
|
||||
Evaluation& xlCO2,
|
||||
Evaluation& ygH2O)
|
||||
{
|
||||
Scalar A = computeA_(temperature, pg);
|
||||
Evaluation A = computeA_(temperature, pg);
|
||||
|
||||
/* salinity: conversion from mass fraction to mol fraction */
|
||||
const Scalar x_NaCl = salinityToMolFrac_(salinity);
|
||||
Scalar x_NaCl = salinityToMolFrac_(salinity);
|
||||
|
||||
// if both phases are present the mole fractions in each phase can be calculate
|
||||
// with the mutual solubility function
|
||||
if (knownPhaseIdx < 0) {
|
||||
Scalar molalityNaCl = molFracToMolality_(x_NaCl); // molality of NaCl //CHANGED
|
||||
Scalar m0_CO2 = molalityCO2inPureWater_(temperature, pg); // molality of CO2 in pure water
|
||||
Scalar gammaStar = activityCoefficient_(temperature, pg, molalityNaCl);// activity coefficient of CO2 in brine
|
||||
Scalar m_CO2 = m0_CO2 / gammaStar; // molality of CO2 in brine
|
||||
Scalar molalityNaCl = moleFracToMolality_(x_NaCl); // molality of NaCl //CHANGED
|
||||
Evaluation m0_CO2 = molalityCO2inPureWater_(temperature, pg); // molality of CO2 in pure water
|
||||
Evaluation gammaStar = activityCoefficient_(temperature, pg, molalityNaCl);// activity coefficient of CO2 in brine
|
||||
Evaluation m_CO2 = m0_CO2 / gammaStar; // molality of CO2 in brine
|
||||
xlCO2 = m_CO2 / (molalityNaCl + 55.508 + m_CO2); // mole fraction of CO2 in brine
|
||||
ygH2O = A * (1 - xlCO2 - x_NaCl); // mole fraction of water in the gas phase
|
||||
}
|
||||
@ -119,25 +122,22 @@ public:
|
||||
// if only liquid phase is present the mole fraction of CO2 in brine is given and
|
||||
// and the virtual equilibrium mole fraction of water in the non-existing gas phase can be estimated
|
||||
// with the mutual solubility function
|
||||
if (knownPhaseIdx == liquidPhaseIdx) {
|
||||
if (knownPhaseIdx == liquidPhaseIdx)
|
||||
ygH2O = A * (1 - xlCO2 - x_NaCl);
|
||||
|
||||
}
|
||||
|
||||
// if only gas phase is present the mole fraction of water in the gas phase is given and
|
||||
// and the virtual equilibrium mole fraction of CO2 in the non-existing liquid phase can be estimated
|
||||
// with the mutual solubility function
|
||||
if (knownPhaseIdx == gasPhaseIdx) {
|
||||
if (knownPhaseIdx == gasPhaseIdx)
|
||||
//y_H2o = fluidstate.
|
||||
xlCO2 = 1 - x_NaCl - ygH2O / A;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Henry coefficent \f$\mathrm{[N/m^2]}\f$ for CO2 in brine.
|
||||
*/
|
||||
static Scalar henry(Scalar temperature)
|
||||
template <class Evaluation>
|
||||
static Evaluation henry(const Evaluation& temperature)
|
||||
{ return fugacityCoefficientCO2(temperature, /*pressure=*/1e5)*1e5; }
|
||||
|
||||
/*!
|
||||
@ -148,34 +148,35 @@ public:
|
||||
* \param T the temperature [K]
|
||||
* \param pg the gas phase pressure [Pa]
|
||||
*/
|
||||
static Scalar fugacityCoefficientCO2(Scalar T, Scalar pg)
|
||||
template <class Evaluation>
|
||||
static Evaluation fugacityCoefficientCO2(const Evaluation& temperature, const Evaluation& pg)
|
||||
{
|
||||
Valgrind::CheckDefined(T);
|
||||
typedef MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
Valgrind::CheckDefined(temperature);
|
||||
Valgrind::CheckDefined(pg);
|
||||
|
||||
Scalar V = 1 / (CO2::gasDensity(T, pg) / CO2::molarMass()) * 1.e6; // molar volume in cm^3/mol
|
||||
Scalar pg_bar = pg / 1.e5; // gas phase pressure in bar
|
||||
Scalar a_CO2 = (7.54e7 - 4.13e4 * T); // mixture parameter of Redlich-Kwong equation
|
||||
static const Scalar b_CO2 = 27.8; // mixture parameter of Redlich-Kwong equation
|
||||
static const Scalar R = IdealGas::R * 10.; // ideal gas constant with unit bar cm^3 /(K mol)
|
||||
Scalar lnPhiCO2, phiCO2;
|
||||
Evaluation V = 1 / (CO2::gasDensity(temperature, pg) / CO2::molarMass()) * 1.e6; // molar volume in cm^3/mol
|
||||
Evaluation pg_bar = pg / 1.e5; // gas phase pressure in bar
|
||||
Evaluation a_CO2 = (7.54e7 - 4.13e4 * temperature); // mixture parameter of Redlich-Kwong equation
|
||||
Scalar b_CO2 = 27.8; // mixture parameter of Redlich-Kwong equation
|
||||
Scalar R = IdealGas::R * 10.; // ideal gas constant with unit bar cm^3 /(K mol)
|
||||
Evaluation lnPhiCO2;
|
||||
|
||||
lnPhiCO2 = std::log(V / (V - b_CO2));
|
||||
lnPhiCO2 = Toolbox::log(V / (V - b_CO2));
|
||||
lnPhiCO2 += b_CO2 / (V - b_CO2);
|
||||
lnPhiCO2 -= 2 * a_CO2 / (R * std::pow(T, 1.5) * b_CO2) * log((V + b_CO2) / V);
|
||||
lnPhiCO2 -= 2 * a_CO2 / (R * Toolbox::pow(temperature, 1.5) * b_CO2) * log((V + b_CO2) / V);
|
||||
lnPhiCO2 +=
|
||||
a_CO2 * b_CO2
|
||||
/ (R
|
||||
* std::pow(T, 1.5)
|
||||
* Toolbox::pow(temperature, 1.5)
|
||||
* b_CO2
|
||||
* b_CO2)
|
||||
* (std::log((V + b_CO2) / V)
|
||||
* (Toolbox::log((V + b_CO2) / V)
|
||||
- b_CO2 / (V + b_CO2));
|
||||
lnPhiCO2 -= std::log(pg_bar * V / (R * T));
|
||||
|
||||
phiCO2 = exp(lnPhiCO2); // fugacity coefficient of CO2
|
||||
return phiCO2;
|
||||
lnPhiCO2 -= Toolbox::log(pg_bar * V / (R * temperature));
|
||||
|
||||
return Toolbox::exp(lnPhiCO2); // fugacity coefficient of CO2
|
||||
}
|
||||
|
||||
/*!
|
||||
@ -183,26 +184,31 @@ public:
|
||||
*
|
||||
* (given in Spycher, Pruess and Ennis-King (2003))
|
||||
*
|
||||
* \param T the temperature [K]
|
||||
* \param temperature the temperature [K]
|
||||
* \param pg the gas phase pressure [Pa]
|
||||
*/
|
||||
static Scalar fugacityCoefficientH2O(Scalar T, Scalar pg)
|
||||
template <class Evaluation>
|
||||
static Evaluation fugacityCoefficientH2O(const Evaluation& temperature, const Evaluation& pg)
|
||||
{
|
||||
Scalar V = 1 / (CO2::gasDensity(T, pg) / CO2::molarMass()) * 1.e6; // molar volume in cm^3/mol
|
||||
Scalar pg_bar = pg / 1.e5; // gas phase pressure in bar
|
||||
Scalar a_CO2 = (7.54e7 - 4.13e4 * T);// mixture parameter of Redlich-Kwong equation
|
||||
static const Scalar a_CO2_H2O = 7.89e7;// mixture parameter of Redlich-Kwong equation
|
||||
static const Scalar b_CO2 = 27.8;// mixture parameter of Redlich-Kwong equation
|
||||
static const Scalar b_H2O = 18.18;// mixture parameter of Redlich-Kwong equation
|
||||
static const Scalar R = IdealGas::R * 10.; // ideal gas constant with unit bar cm^3 /(K mol)
|
||||
Scalar lnPhiH2O, phiH2O;
|
||||
typedef MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
lnPhiH2O = log(V / (V - b_CO2)) + b_H2O / (V - b_CO2) - 2 * a_CO2_H2O
|
||||
/ (R * pow(T, 1.5) * b_CO2) * log((V + b_CO2) / V) + a_CO2
|
||||
* b_H2O / (R * pow(T, 1.5) * b_CO2 * b_CO2) * (log((V + b_CO2)
|
||||
/ V) - b_CO2 / (V + b_CO2)) - log(pg_bar * V / (R * T));
|
||||
phiH2O = exp(lnPhiH2O); // fugacity coefficient of H2O
|
||||
return phiH2O;
|
||||
const Evaluation& V = 1 / (CO2::gasDensity(temperature, pg) / CO2::molarMass()) * 1.e6; // molar volume in cm^3/mol
|
||||
const Evaluation& pg_bar = pg / 1.e5; // gas phase pressure in bar
|
||||
const Evaluation& a_CO2 = (7.54e7 - 4.13e4 * temperature);// mixture parameter of Redlich-Kwong equation
|
||||
Scalar a_CO2_H2O = 7.89e7;// mixture parameter of Redlich-Kwong equation
|
||||
Scalar b_CO2 = 27.8;// mixture parameter of Redlich-Kwong equation
|
||||
Scalar b_H2O = 18.18;// mixture parameter of Redlich-Kwong equation
|
||||
Scalar R = IdealGas::R * 10.; // ideal gas constant with unit bar cm^3 /(K mol)
|
||||
Evaluation lnPhiH2O;
|
||||
|
||||
lnPhiH2O =
|
||||
Toolbox::log(V/(V - b_CO2))
|
||||
+ b_H2O/(V - b_CO2) - 2*a_CO2_H2O
|
||||
/ (R*Toolbox::pow(temperature, 1.5)*b_CO2)*Toolbox::log((V + b_CO2)/V)
|
||||
+ a_CO2*b_H2O/(R*Toolbox::pow(temperature, 1.5)*b_CO2*b_CO2)
|
||||
*(Toolbox::log((V + b_CO2)/V) - b_CO2/(V + b_CO2))
|
||||
- Toolbox::log(pg_bar*V/(R*temperature));
|
||||
return Toolbox::exp(lnPhiH2O); // fugacity coefficient of H2O
|
||||
}
|
||||
|
||||
private:
|
||||
@ -227,12 +233,10 @@ private:
|
||||
*
|
||||
* \param x_NaCl mole fraction of NaCL in brine [mol/mol]
|
||||
*/
|
||||
static Scalar molFracToMolality_(Scalar x_NaCl) {
|
||||
|
||||
static Scalar moleFracToMolality_(Scalar x_NaCl)
|
||||
{
|
||||
// conversion from mol fraction to molality (dissolved CO2 neglected)
|
||||
const Scalar mol_NaCl = 55.508 * x_NaCl / (1 - x_NaCl);
|
||||
|
||||
return mol_NaCl;
|
||||
return 55.508 * x_NaCl / (1 - x_NaCl);
|
||||
}
|
||||
|
||||
/*!
|
||||
@ -242,13 +246,14 @@ private:
|
||||
* \param temperature The temperature [K]
|
||||
* \param pg The gas phase pressure [Pa]
|
||||
*/
|
||||
static Scalar molalityCO2inPureWater_(Scalar temperature, Scalar pg) {
|
||||
Scalar A = computeA_(temperature, pg); // according to Spycher, Pruess and Ennis-King (2003)
|
||||
Scalar B = computeB_(temperature, pg); // according to Spycher, Pruess and Ennis-King (2003)
|
||||
Scalar yH2OinGas = (1 - B) / (1. / A - B); // equilibrium mol fraction of H2O in the gas phase
|
||||
Scalar xCO2inWater = B * (1 - yH2OinGas); // equilibrium mol fraction of CO2 in the water phase
|
||||
Scalar molalityCO2 = (xCO2inWater * 55.508) / (1 - xCO2inWater); // CO2 molality
|
||||
return molalityCO2;
|
||||
template <class Evaluation>
|
||||
static Evaluation molalityCO2inPureWater_(const Evaluation& temperature, const Evaluation& pg)
|
||||
{
|
||||
const Evaluation& A = computeA_(temperature, pg); // according to Spycher, Pruess and Ennis-King (2003)
|
||||
const Evaluation& B = computeB_(temperature, pg); // according to Spycher, Pruess and Ennis-King (2003)
|
||||
const Evaluation& yH2OinGas = (1 - B) / (1. / A - B); // equilibrium mol fraction of H2O in the gas phase
|
||||
const Evaluation& xCO2inWater = B * (1 - yH2OinGas); // equilibrium mol fraction of CO2 in the water phase
|
||||
return (xCO2inWater * 55.508) / (1 - xCO2inWater); // CO2 molality
|
||||
}
|
||||
|
||||
/*!
|
||||
@ -260,16 +265,18 @@ private:
|
||||
* \param pg the gas phase pressure [Pa]
|
||||
* \param molalityNaCl molality of NaCl (mol NaCl / kg water)
|
||||
*/
|
||||
static Scalar activityCoefficient_(Scalar temperature,
|
||||
Scalar pg,
|
||||
Scalar molalityNaCl)
|
||||
template <class Evaluation>
|
||||
static Evaluation activityCoefficient_(const Evaluation& temperature,
|
||||
const Evaluation& pg,
|
||||
Scalar molalityNaCl)
|
||||
{
|
||||
Scalar lambda = computeLambda_(temperature, pg); // lambda_{CO2-Na+}
|
||||
Scalar xi = computeXi_(temperature, pg); // Xi_{CO2-Na+-Cl-}
|
||||
Scalar lnGammaStar = 2 * lambda * molalityNaCl + xi * molalityNaCl
|
||||
* molalityNaCl;
|
||||
Scalar gammaStar = exp(lnGammaStar);
|
||||
return gammaStar; // molal activity coefficient of CO2 in brine
|
||||
typedef MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
const Evaluation& lambda = computeLambda_(temperature, pg); // lambda_{CO2-Na+}
|
||||
const Evaluation& xi = computeXi_(temperature, pg); // Xi_{CO2-Na+-Cl-}
|
||||
const Evaluation& lnGammaStar =
|
||||
2*molalityNaCl*lambda + xi*molalityNaCl*molalityNaCl;
|
||||
return Toolbox::exp(lnGammaStar);
|
||||
}
|
||||
|
||||
/*!
|
||||
@ -280,17 +287,18 @@ private:
|
||||
* \param T the temperature [K]
|
||||
* \param pg the gas phase pressure [Pa]
|
||||
*/
|
||||
static Scalar computeA_(Scalar T, Scalar pg)
|
||||
template <class Evaluation>
|
||||
static Evaluation computeA_(const Evaluation& temperature, const Evaluation& pg)
|
||||
{
|
||||
Scalar deltaP = pg / 1e5 - 1; // pressure range [bar] from p0 = 1bar to pg[bar]
|
||||
const Scalar v_av_H2O = 18.1; // average partial molar volume of H2O [cm^3/mol]
|
||||
const Scalar R = IdealGas::R * 10;
|
||||
Scalar k0_H2O = equilibriumConstantH2O_(T); // equilibrium constant for H2O at 1 bar
|
||||
Scalar phi_H2O = fugacityCoefficientH2O(T, pg); // fugacity coefficient of H2O for the water-CO2 system
|
||||
Scalar pg_bar = pg / 1.e5;
|
||||
Scalar A = k0_H2O / (phi_H2O * pg_bar) * exp(deltaP * v_av_H2O / (R * T));
|
||||
return A;
|
||||
typedef MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
const Evaluation& deltaP = pg / 1e5 - 1; // pressure range [bar] from p0 = 1bar to pg[bar]
|
||||
Scalar v_av_H2O = 18.1; // average partial molar volume of H2O [cm^3/mol]
|
||||
Scalar R = IdealGas::R * 10;
|
||||
const Evaluation& k0_H2O = equilibriumConstantH2O_(temperature); // equilibrium constant for H2O at 1 bar
|
||||
const Evaluation& phi_H2O = fugacityCoefficientH2O(temperature, pg); // fugacity coefficient of H2O for the water-CO2 system
|
||||
const Evaluation& pg_bar = pg / 1.e5;
|
||||
return k0_H2O/(phi_H2O*pg_bar)*Toolbox::exp(deltaP*v_av_H2O/(R*temperature));
|
||||
}
|
||||
|
||||
/*!
|
||||
@ -298,72 +306,80 @@ private:
|
||||
* the mutual solubility in the water-CO2 system.
|
||||
* Given in Spycher, Pruess and Ennis-King (2003)
|
||||
*
|
||||
* \param T the temperature [K]
|
||||
* \param temperature the temperature [K]
|
||||
* \param pg the gas phase pressure [Pa]
|
||||
*/
|
||||
static Scalar computeB_(Scalar T, Scalar pg) {
|
||||
Scalar deltaP = pg / 1e5 - 1; // pressure range [bar] from p0 = 1bar to pg[bar]
|
||||
template <class Evaluation>
|
||||
static Evaluation computeB_(const Evaluation& temperature, const Evaluation& pg)
|
||||
{
|
||||
typedef MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
const Evaluation& deltaP = pg / 1e5 - 1; // pressure range [bar] from p0 = 1bar to pg[bar]
|
||||
const Scalar v_av_CO2 = 32.6; // average partial molar volume of CO2 [cm^3/mol]
|
||||
const Scalar R = IdealGas::R * 10;
|
||||
Scalar k0_CO2 = equilibriumConstantCO2_(T); // equilibrium constant for CO2 at 1 bar
|
||||
Scalar phi_CO2 = fugacityCoefficientCO2(T, pg); // fugacity coefficient of CO2 for the water-CO2 system
|
||||
Scalar pg_bar = pg / 1.e5;
|
||||
Scalar B = phi_CO2 * pg_bar / (55.508 * k0_CO2) * exp(-(deltaP
|
||||
* v_av_CO2) / (R * T));
|
||||
return B;
|
||||
const Evaluation& k0_CO2 = equilibriumConstantCO2_(temperature); // equilibrium constant for CO2 at 1 bar
|
||||
const Evaluation& phi_CO2 = fugacityCoefficientCO2(temperature, pg); // fugacity coefficient of CO2 for the water-CO2 system
|
||||
const Evaluation& pg_bar = pg / 1.e5;
|
||||
return phi_CO2*pg_bar/(55.508*k0_CO2)*Toolbox::exp(-(deltaP*v_av_CO2)/(R*temperature));
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the parameter lambda, which is needed for the
|
||||
* calculation of the CO2 activity coefficient in the brine-CO2 system.
|
||||
* Given in Spycher and Pruess (2005)
|
||||
* \param T the temperature [K]
|
||||
* \param temperature the temperature [K]
|
||||
* \param pg the gas phase pressure [Pa]
|
||||
*/
|
||||
static Scalar computeLambda_(Scalar T, Scalar pg)
|
||||
template <class Evaluation>
|
||||
static Evaluation computeLambda_(const Evaluation& temperature, const Evaluation& pg)
|
||||
{
|
||||
Scalar lambda;
|
||||
static const Scalar c[6] = { -0.411370585, 6.07632013E-4, 97.5347708,
|
||||
-0.0237622469, 0.0170656236, 1.41335834E-5 };
|
||||
typedef MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
Scalar pg_bar = pg / 1.0E5; /* conversion from Pa to bar */
|
||||
lambda = c[0] + c[1] * T + c[2] / T + c[3] * pg_bar / T + c[4] * pg_bar
|
||||
/ (630.0 - T) + c[5] * T * std::log(pg_bar);
|
||||
static const Scalar c[6] =
|
||||
{ -0.411370585, 6.07632013E-4, 97.5347708, -0.0237622469, 0.0170656236, 1.41335834E-5 };
|
||||
|
||||
return lambda;
|
||||
Evaluation pg_bar = pg / 1.0E5; /* conversion from Pa to bar */
|
||||
return
|
||||
c[0]
|
||||
+ c[1]*temperature
|
||||
+ c[2]/temperature
|
||||
+ c[3]*pg_bar/temperature
|
||||
+ c[4]*pg_bar/(630.0 - temperature)
|
||||
+ c[5]*temperature*Toolbox::log(pg_bar);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the parameter xi, which is needed for the
|
||||
* calculation of the CO2 activity coefficient in the brine-CO2 system.
|
||||
* Given in Spycher and Pruess (2005)
|
||||
* \param T the temperature [K]
|
||||
* \param temperature the temperature [K]
|
||||
* \param pg the gas phase pressure [Pa]
|
||||
*/
|
||||
static Scalar computeXi_(Scalar T, Scalar pg)
|
||||
template <class Evaluation>
|
||||
static Evaluation computeXi_(const Evaluation& temperature, const Evaluation& pg)
|
||||
{
|
||||
Scalar xi;
|
||||
static const Scalar c[4] = { 3.36389723E-4, -1.98298980E-5,
|
||||
2.12220830E-3, -5.24873303E-3 };
|
||||
static const Scalar c[4] =
|
||||
{ 3.36389723E-4, -1.98298980E-5, 2.12220830E-3, -5.24873303E-3 };
|
||||
|
||||
Scalar pg_bar = pg / 1.0E5; /* conversion from Pa to bar */
|
||||
xi = c[0] + c[1] * T + c[2] * pg_bar / T + c[3] * pg_bar / (630.0 - T);
|
||||
|
||||
return xi;
|
||||
Evaluation pg_bar = pg / 1.0E5; /* conversion from Pa to bar */
|
||||
return c[0] + c[1]*temperature + c[2]*pg_bar/temperature + c[3]*pg_bar/(630.0 - temperature);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the equilibrium constant for CO2, which is needed for the
|
||||
* calculation of the mutual solubility in the water-CO2 system
|
||||
* Given in Spycher, Pruess and Ennis-King (2003)
|
||||
* \param T the temperature [K]
|
||||
* \param temperature the temperature [K]
|
||||
*/
|
||||
static Scalar equilibriumConstantCO2_(Scalar T)
|
||||
template <class Evaluation>
|
||||
static Evaluation equilibriumConstantCO2_(const Evaluation& temperature)
|
||||
{
|
||||
Scalar TinC = T - 273.15; //temperature in °C
|
||||
typedef MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
Evaluation temperatureCelcius = temperature - 273.15;
|
||||
static const Scalar c[3] = { 1.189, 1.304e-2, -5.446e-5 };
|
||||
Scalar logk0_CO2 = c[0] + c[1] * TinC + c[2] * TinC * TinC;
|
||||
Scalar k0_CO2 = pow(10, logk0_CO2);
|
||||
Evaluation logk0_CO2 = c[0] + temperatureCelcius*(c[1] + temperatureCelcius*c[2]);
|
||||
Evaluation k0_CO2 = Toolbox::pow(10.0, logk0_CO2);
|
||||
return k0_CO2;
|
||||
}
|
||||
|
||||
@ -371,16 +387,18 @@ private:
|
||||
* \brief Returns the equilibrium constant for H2O, which is needed for the
|
||||
* calculation of the mutual solubility in the water-CO2 system
|
||||
* Given in Spycher, Pruess and Ennis-King (2003)
|
||||
* \param T the temperature [K]
|
||||
* \param temperature the temperature [K]
|
||||
*/
|
||||
static Scalar equilibriumConstantH2O_(Scalar T)
|
||||
template <class Evaluation>
|
||||
static Evaluation equilibriumConstantH2O_(const Evaluation& temperature)
|
||||
{
|
||||
Scalar TinC = T - 273.15; //temperature in °C
|
||||
typedef MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
Evaluation temperatureCelcius = temperature - 273.15;
|
||||
static const Scalar c[4] = { -2.209, 3.097e-2, -1.098e-4, 2.048e-7 };
|
||||
Scalar logk0_H2O = c[0] + c[1] * TinC + c[2] * TinC * TinC + c[3]
|
||||
* TinC * TinC * TinC;
|
||||
Scalar k0_H2O = pow(10, logk0_H2O);
|
||||
return k0_H2O;
|
||||
Evaluation logk0_H2O =
|
||||
c[0] + temperatureCelcius*(c[1] + temperatureCelcius*(c[2] + temperatureCelcius*c[3]));
|
||||
return Toolbox::pow(10.0, logk0_H2O);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
@ -25,6 +25,7 @@
|
||||
#define OPM_FULLERMETHOD_HPP
|
||||
|
||||
#include <opm/material/common/Means.hpp>
|
||||
#include <opm/material/common/MathToolbox.hpp>
|
||||
|
||||
#include <cmath>
|
||||
|
||||
@ -48,18 +49,20 @@ namespace BinaryCoeff {
|
||||
* See: R. Reid, et al.: The Properties of Gases and Liquids, 4th
|
||||
* edition, McGraw-Hill, 1987, pp. 587-588
|
||||
*/
|
||||
template <class Scalar>
|
||||
inline Scalar fullerMethod(const Scalar *M, // molar masses [g/mol]
|
||||
const Scalar *SigmaNu, // atomic diffusion volume
|
||||
const Scalar temperature, // [K]
|
||||
const Scalar pressure) // [Pa]
|
||||
template <class Scalar, class Evaluation = Scalar>
|
||||
inline Evaluation fullerMethod(const Scalar *M, // molar masses [g/mol]
|
||||
const Scalar *SigmaNu, // atomic diffusion volume
|
||||
const Evaluation& temperature, // [K]
|
||||
const Evaluation& pressure) // [Pa]
|
||||
{
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
// "effective" molar mass in [g/m^3]
|
||||
Scalar Mab = Opm::harmonicMean(M[0], M[1]);
|
||||
|
||||
// Fuller's method
|
||||
Scalar tmp = std::pow(SigmaNu[0], 1./3) + std::pow(SigmaNu[1], 1./3);
|
||||
return 1e-4 * (143.0*std::pow(temperature, 1.75))/(pressure*std::sqrt(Mab)*tmp*tmp);
|
||||
const Evaluation& tmp = std::pow(SigmaNu[0], 1./3) + std::pow(SigmaNu[1], 1./3);
|
||||
return 1e-4 * (143.0*Toolbox::pow(temperature, 1.75))/(pressure*std::sqrt(Mab)*tmp*tmp);
|
||||
}
|
||||
|
||||
} // namespace BinaryCoeff
|
||||
|
||||
@ -25,12 +25,12 @@
|
||||
#ifndef OPM_BINARY_COEFF_H2O_AIR_HPP
|
||||
#define OPM_BINARY_COEFF_H2O_AIR_HPP
|
||||
|
||||
#include <opm/material/common/MathToolbox.hpp>
|
||||
|
||||
#include <cmath>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
namespace BinaryCoeff
|
||||
{
|
||||
namespace Opm {
|
||||
namespace BinaryCoeff {
|
||||
|
||||
/*!
|
||||
* \ingroup Binarycoefficients
|
||||
@ -49,12 +49,12 @@ public:
|
||||
* page 29 Formula (2.9) (nach Tchobanoglous & Schroeder, 1985)
|
||||
*
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar henry(Scalar temperature)
|
||||
template <class Evaluation>
|
||||
static Evaluation henry(const Evaluation& temperature)
|
||||
{
|
||||
Scalar r = (0.8942+1.47*std::exp(-0.04394*(temperature-273.15)))*1.E-10;
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
return 1./r;
|
||||
return 1.0/((0.8942+1.47*Toolbox::exp(-0.04394*(temperature-273.15)))*1.E-10);
|
||||
}
|
||||
|
||||
/*!
|
||||
@ -68,18 +68,17 @@ public:
|
||||
* Dep. of Agricultural and Chemical Engineering, Colorado State University,
|
||||
* Fort Collins, 1981.
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
template <class Evaluation>
|
||||
static Evaluation gasDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
|
||||
{
|
||||
const Scalar Theta=1.8;
|
||||
const Scalar Daw=2.13e-5; /* reference value */
|
||||
const Scalar pg0=1.e5; /* reference pressure */
|
||||
const Scalar T0=273.15; /* reference temperature */
|
||||
Scalar Dgaw;
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
Dgaw=Daw*(pg0/pressure)*std::pow((temperature/T0),Theta);
|
||||
double Theta=1.8;
|
||||
double Daw=2.13e-5; /* reference value */
|
||||
double pg0=1.e5; /* reference pressure */
|
||||
double T0=273.15; /* reference temperature */
|
||||
|
||||
return Dgaw;
|
||||
return Daw*(pg0/pressure)*Toolbox::pow((temperature/T0),Theta);
|
||||
}
|
||||
|
||||
/*!
|
||||
@ -102,12 +101,12 @@ public:
|
||||
* Oxygen in Water", Journal of Chemical Engineering and Data,
|
||||
* Vol. 12, No. 1, pp. 111-115, 1967
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
template <class Evaluation>
|
||||
static Evaluation liquidDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
|
||||
{
|
||||
const Scalar Texp = 273.15 + 25; // [K]
|
||||
const Scalar Dexp = 2.01e-9; // [m^2/s]
|
||||
return Dexp * temperature/Texp;
|
||||
const double Texp = 273.15 + 25; // [K]
|
||||
const double Dexp = 2.01e-9; // [m^2/s]
|
||||
return Dexp/Texp*temperature;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
@ -50,8 +50,8 @@ public:
|
||||
* Temperatures"
|
||||
* http://www.iapws.org/relguide/HenGuide.pdf
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar henry(Scalar temperature)
|
||||
template <class Scalar, class Evaluation = Scalar>
|
||||
static Evaluation henry(const Evaluation& temperature)
|
||||
{
|
||||
const Scalar E = 1672.9376;
|
||||
const Scalar F = 28.1751;
|
||||
@ -66,8 +66,8 @@ public:
|
||||
*
|
||||
* To calculate the values, the \ref fullerMethod is used.
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
template <class Scalar, class Evaluation = Scalar>
|
||||
static Evaluation gasDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
|
||||
{
|
||||
typedef Opm::H2O<Scalar> H2O;
|
||||
typedef Opm::SimpleCO2<Scalar> CO2;
|
||||
@ -83,8 +83,8 @@ public:
|
||||
/*!
|
||||
* \brief Diffusion coefficent [m^2/s] for molecular CO2 in liquid water.
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
template <class Scalar, class Evaluation = Scalar>
|
||||
static Evaluation liquidDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
|
||||
{ OPM_THROW(std::runtime_error, "Not implemented: Binary liquid diffusion coefficients of CO2 and CH4"); }
|
||||
};
|
||||
|
||||
|
||||
@ -44,13 +44,13 @@ public:
|
||||
*
|
||||
* Sanders1999 Henry collection
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar henry(Scalar temperature)
|
||||
template <class Evaluation>
|
||||
static Evaluation henry(const Evaluation& temperature)
|
||||
{
|
||||
// after Sanders
|
||||
Scalar sanderH = 1.7e-1; // [M/atm]
|
||||
double sanderH = 1.7e-1; // [M/atm]
|
||||
//conversion to our Henry definition
|
||||
Scalar opmH = sanderH / 101.325; // has now [(mol/m^3)/Pa]
|
||||
double opmH = sanderH / 101.325; // has now [(mol/m^3)/Pa]
|
||||
opmH *= 18.02e-6; // multiplied by molar volume of reference phase = water
|
||||
return 1.0/opmH; // [Pa]
|
||||
}
|
||||
@ -59,38 +59,39 @@ public:
|
||||
* \brief Binary diffusion coefficent [m^2/s] for molecular water and mesitylene.
|
||||
*
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
template <class Evaluation>
|
||||
static Evaluation gasDiffCoeff(Evaluation temperature, Evaluation pressure)
|
||||
{
|
||||
typedef Opm::H2O<Scalar> H2O;
|
||||
typedef Opm::Mesitylene<Scalar> Mesitylene;
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
temperature = std::max(temperature, 1e-9); // regularization
|
||||
temperature = std::min(temperature, 500.0); // regularization
|
||||
pressure = std::max(pressure, 0.0); // regularization
|
||||
pressure = std::min(pressure, 1e8); // regularization
|
||||
typedef Opm::H2O<double> H2O;
|
||||
typedef Opm::Mesitylene<double> Mesitylene;
|
||||
|
||||
const Scalar M_m = 1e3*Mesitylene::molarMass(); // [g/mol] molecular weight of mesitylene
|
||||
const Scalar M_w = 1e3*H2O::molarMass(); // [g/mol] molecular weight of water
|
||||
const Scalar Tb_m = 437.9; // [K] boiling temperature of mesitylen
|
||||
const Scalar Tb_w = 373.15; // [K] boiling temperature of water (at p_atm)
|
||||
const Scalar V_B_w = 18.0; // [cm^3/mol] LeBas molal volume of water
|
||||
const Scalar sigma_w = 1.18*std::pow(V_B_w, 0.333); // charact. length of air
|
||||
const Scalar T_scal_w = 1.15*Tb_w; // [K] (molec. energy of attraction/Boltzmann constant)
|
||||
const Scalar V_B_m = 162.6; // [cm^3/mol] LeBas molal volume of mesitylen
|
||||
const Scalar sigma_m = 1.18*std::pow(V_B_m, 0.333); // charact. length of mesitylen
|
||||
const Scalar sigma_wm = 0.5*(sigma_w + sigma_m);
|
||||
const Scalar T_scal_m = 1.15*Tb_m;
|
||||
const Scalar T_scal_wm = std::sqrt(T_scal_w*T_scal_m);
|
||||
temperature = Toolbox::max(temperature, 1e-9); // regularization
|
||||
temperature = Toolbox::min(temperature, 500.0); // regularization
|
||||
pressure = Toolbox::max(pressure, 0.0); // regularization
|
||||
pressure = Toolbox::min(pressure, 1e8); // regularization
|
||||
|
||||
Scalar T_star = temperature/T_scal_wm;
|
||||
T_star = std::max(T_star, 1e-5); // regularization
|
||||
const double M_m = 1e3*Mesitylene::molarMass(); // [g/mol] molecular weight of mesitylene
|
||||
const double M_w = 1e3*H2O::molarMass(); // [g/mol] molecular weight of water
|
||||
const double Tb_m = 437.9; // [K] boiling temperature of mesitylen
|
||||
const double Tb_w = 373.15; // [K] boiling temperature of water (at p_atm)
|
||||
const double V_B_w = 18.0; // [cm^3/mol] LeBas molal volume of water
|
||||
const double sigma_w = 1.18*std::pow(V_B_w, 0.333); // charact. length of air
|
||||
const double T_scal_w = 1.15*Tb_w; // [K] (molec. energy of attraction/Boltzmann constant)
|
||||
const double V_B_m = 162.6; // [cm^3/mol] LeBas molal volume of mesitylen
|
||||
const double sigma_m = 1.18*std::pow(V_B_m, 0.333); // charact. length of mesitylen
|
||||
const double sigma_wm = 0.5*(sigma_w + sigma_m);
|
||||
const double T_scal_m = 1.15*Tb_m;
|
||||
const double T_scal_wm = std::sqrt(T_scal_w*T_scal_m);
|
||||
|
||||
const Scalar Omega = 1.06036/std::pow(T_star,0.1561) + 0.193/std::exp(T_star*0.47635)
|
||||
+ 1.03587/std::exp(T_star*1.52996) + 1.76474/std::exp(T_star*3.89411);
|
||||
const Scalar B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_w + 1.0/M_m);
|
||||
const Scalar Mr = (M_w + M_m)/(M_w*M_m);
|
||||
const Scalar D_wm = (B_*std::pow(temperature,1.6)*std::sqrt(Mr))
|
||||
const Evaluation T_star = Toolbox::max(temperature/T_scal_wm, 1e-5);
|
||||
|
||||
const Evaluation& Omega = 1.06036/Toolbox::pow(T_star,0.1561) + 0.193/Toolbox::exp(T_star*0.47635)
|
||||
+ 1.03587/Toolbox::exp(T_star*1.52996) + 1.76474/Toolbox::exp(T_star*3.89411);
|
||||
const double B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_w + 1.0/M_m);
|
||||
const double Mr = (M_w + M_m)/(M_w*M_m);
|
||||
const Evaluation D_wm = (B_*Toolbox::pow(temperature,1.6)*std::sqrt(Mr))
|
||||
/(1e-5*pressure*std::pow(sigma_wm, 2.0)*Omega); // [cm^2/s]
|
||||
|
||||
return D_wm*1e-4; // [m^2/s]
|
||||
@ -99,12 +100,11 @@ public:
|
||||
/*!
|
||||
* \brief Diffusion coefficent [m^2/s] for mesitylene in liquid water.
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
template <class Evaluation>
|
||||
static Evaluation liquidDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
|
||||
{
|
||||
// This is just an order of magnitude estimate. Please
|
||||
// improve!
|
||||
return 1.e-9;
|
||||
// This is just an order of magnitude estimate. Please improve!
|
||||
return 1e-9;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
@ -46,13 +46,13 @@ public:
|
||||
*
|
||||
* \copydetails Opm::henryIAPWS
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar henry(Scalar temperature)
|
||||
template <class Evaluation>
|
||||
static Evaluation henry(const Evaluation& temperature)
|
||||
{
|
||||
const Scalar E = 2388.8777;
|
||||
const Scalar F = -14.9593;
|
||||
const Scalar G = 42.0179;
|
||||
const Scalar H = -29.4396;
|
||||
const double E = 2388.8777;
|
||||
const double F = -14.9593;
|
||||
const double G = 42.0179;
|
||||
const double H = -29.4396;
|
||||
|
||||
return henryIAPWS(E, F, G, H, temperature);
|
||||
}
|
||||
@ -64,16 +64,16 @@ public:
|
||||
* \param temperature the temperature \f$\mathrm{[K]}\f$
|
||||
* \param pressure the phase pressure \f$\mathrm{[Pa]}\f$
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
template <class Evaluation>
|
||||
static Evaluation gasDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
|
||||
{
|
||||
typedef Opm::H2O<Scalar> H2O;
|
||||
typedef Opm::N2<Scalar> N2;
|
||||
typedef Opm::H2O<double> H2O;
|
||||
typedef Opm::N2<double> N2;
|
||||
|
||||
// atomic diffusion volumes
|
||||
const Scalar SigmaNu[2] = { 13.1 /* H2O */, 18.5 /* N2 */ };
|
||||
const double SigmaNu[2] = { 13.1 /* H2O */, 18.5 /* N2 */ };
|
||||
// molar masses [g/mol]
|
||||
const Scalar M[2] = { H2O::molarMass()*Scalar(1e3), N2::molarMass()*Scalar(1e3) };
|
||||
const double M[2] = { H2O::molarMass()*1e3, N2::molarMass()*1e3 };
|
||||
|
||||
return fullerMethod(M, SigmaNu, temperature, pressure);
|
||||
}
|
||||
@ -96,11 +96,11 @@ public:
|
||||
* Oxygen in Water", Journal of Chemical Engineering and Data,
|
||||
* Vol. 12, No. 1, pp. 111-115, 1967
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
template <class Evaluation>
|
||||
static Evaluation liquidDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
|
||||
{
|
||||
const Scalar Texp = 273.15 + 25; // [K]
|
||||
const Scalar Dexp = 2.01e-9; // [m^2/s]
|
||||
const double Texp = 273.15 + 25; // [K]
|
||||
const double Dexp = 2.01e-9; // [m^2/s]
|
||||
|
||||
return Dexp * temperature/Texp;
|
||||
}
|
||||
|
||||
@ -26,10 +26,8 @@
|
||||
#include <opm/material/components/H2O.hpp>
|
||||
#include <opm/material/components/Xylene.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
namespace BinaryCoeff
|
||||
{
|
||||
namespace Opm {
|
||||
namespace BinaryCoeff {
|
||||
|
||||
/*!
|
||||
* \brief Binary coefficients for water and xylene.
|
||||
@ -45,13 +43,13 @@ public:
|
||||
* Sanders1999 Henry collection
|
||||
*/
|
||||
|
||||
template <class Scalar>
|
||||
static Scalar henry(Scalar temperature)
|
||||
template <class Evaluation>
|
||||
static Evaluation henry(const Evaluation& temperature)
|
||||
{
|
||||
// after Sanders
|
||||
Scalar sanderH = 1.5e-1; //[M/atm]
|
||||
double sanderH = 1.5e-1; //[M/atm]
|
||||
//conversion to our Henry definition
|
||||
Scalar opmH = sanderH / 101.325; // has now [(mol/m^3)/Pa]
|
||||
double opmH = sanderH / 101.325; // has now [(mol/m^3)/Pa]
|
||||
opmH *= 18.02e-6; //multiplied by molar volume of reference phase = water
|
||||
return 1.0/opmH; // [Pa]
|
||||
}
|
||||
@ -60,41 +58,40 @@ public:
|
||||
* \brief Binary diffusion coefficent [m^2/s] for molecular water and xylene.
|
||||
*
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
template <class Evaluation>
|
||||
static Evaluation gasDiffCoeff(Evaluation temperature, Evaluation pressure)
|
||||
{
|
||||
typedef Opm::H2O<Scalar> H2O;
|
||||
typedef Opm::Xylene<Scalar> Xylene;
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
typedef Opm::H2O<double> H2O;
|
||||
typedef Opm::Xylene<double> Xylene;
|
||||
|
||||
temperature = std::max(temperature, 1e-9); // regularization
|
||||
temperature = std::min(temperature, 500.0); // regularization
|
||||
pressure = std::max(pressure, 0.0); // regularization
|
||||
pressure = std::min(pressure, 1e8); // regularization
|
||||
temperature = Toolbox::max(temperature, 1e-9); // regularization
|
||||
temperature = Toolbox::min(temperature, 500.0); // regularization
|
||||
pressure = Toolbox::max(pressure, 0.0); // regularization
|
||||
pressure = Toolbox::min(pressure, 1e8); // regularization
|
||||
|
||||
const Scalar M_x = 1e3*Xylene::molarMass(); // [g/mol] molecular weight of xylene
|
||||
const Scalar M_w = 1e3*H2O::molarMass(); // [g/mol] molecular weight of water
|
||||
const Scalar Tb_x = 412.9; // [K] boiling temperature of xylene
|
||||
const Scalar Tb_w = 373.15; // [K] boiling temperature of water (at p_atm)
|
||||
const Scalar V_B_w = 18.0; // [cm^3/mol] LeBas molal volume of water
|
||||
const Scalar sigma_w = 1.18*std::pow(V_B_w, 0.333); // charact. length of air
|
||||
const Scalar T_scal_w = 1.15*Tb_w; // [K] (molec. energy of attraction/Boltzmann constant)
|
||||
const Scalar V_B_x = 140.4; // [cm^3/mol] LeBas molal volume of xylene
|
||||
const Scalar sigma_x = 1.18*std::pow(V_B_x, 0.333); // charact. length of xylene
|
||||
const Scalar sigma_wx = 0.5*(sigma_w + sigma_x);
|
||||
const Scalar T_scal_x = 1.15*Tb_x;
|
||||
const Scalar T_scal_wx = std::sqrt(T_scal_w*T_scal_x);
|
||||
const double M_x = 1e3*Xylene::molarMass(); // [g/mol] molecular weight of xylene
|
||||
const double M_w = 1e3*H2O::molarMass(); // [g/mol] molecular weight of water
|
||||
const double Tb_x = 412.9; // [K] boiling temperature of xylene
|
||||
const double Tb_w = 373.15; // [K] boiling temperature of water (at p_atm)
|
||||
const double V_B_w = 18.0; // [cm^3/mol] LeBas molal volume of water
|
||||
const double sigma_w = 1.18*std::pow(V_B_w, 0.333); // charact. length of air
|
||||
const double T_scal_w = 1.15*Tb_w; // [K] (molec. energy of attraction/Boltzmann constant)
|
||||
const double V_B_x = 140.4; // [cm^3/mol] LeBas molal volume of xylene
|
||||
const double sigma_x = 1.18*std::pow(V_B_x, 0.333); // charact. length of xylene
|
||||
const double sigma_wx = 0.5*(sigma_w + sigma_x);
|
||||
const double T_scal_x = 1.15*Tb_x;
|
||||
const double T_scal_wx = std::sqrt(T_scal_w*T_scal_x);
|
||||
|
||||
Scalar T_star = temperature/T_scal_wx;
|
||||
T_star = std::max(T_star, 1e-5); // regularization
|
||||
const Evaluation& T_star = Toolbox::max(temperature/T_scal_wx, 1e-5);
|
||||
|
||||
const Scalar Omega = 1.06036/std::pow(T_star, 0.1561) + 0.193/std::exp(T_star*0.47635)
|
||||
+ 1.03587/std::exp(T_star*1.52996) + 1.76474/std::exp(T_star*3.89411);
|
||||
const Scalar B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_w + 1.0/M_x);
|
||||
const Scalar Mr = (M_w + M_x)/(M_w*M_x);
|
||||
const Scalar D_wx = (B_*std::pow(temperature,1.6)*std::sqrt(Mr))
|
||||
/(1e-5*pressure*std::pow(sigma_wx, 2.0)*Omega); // [cm^2/s]
|
||||
|
||||
return D_wx*1e-4; // [m^2/s]
|
||||
const Evaluation& Omega = 1.06036/Toolbox::pow(T_star, 0.1561) + 0.193/Toolbox::exp(T_star*0.47635)
|
||||
+ 1.03587/Toolbox::exp(T_star*1.52996) + 1.76474/Toolbox::exp(T_star*3.89411);
|
||||
const double B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_w + 1.0/M_x);
|
||||
const double Mr = (M_w + M_x)/(M_w*M_x);
|
||||
return 1e-4
|
||||
*(B_*Toolbox::pow(temperature,1.6)*std::sqrt(Mr))
|
||||
/(1e-5*pressure*std::pow(sigma_wx, 2.0)*Omega);
|
||||
}
|
||||
|
||||
/*!
|
||||
@ -102,8 +99,8 @@ public:
|
||||
*
|
||||
* \todo
|
||||
*/
|
||||
template <class Scalar>
|
||||
static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
|
||||
template <class Evaluation>
|
||||
static Evaluation liquidDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
|
||||
{
|
||||
return 1.e-9; // This is just an order of magnitude. Please improve it!
|
||||
}
|
||||
|
||||
@ -31,27 +31,25 @@ namespace Opm
|
||||
{
|
||||
/*!
|
||||
* \ingroup Binarycoefficients
|
||||
* \brief The Henry constants in liquid water using the IAPWS 2004
|
||||
* formulation.
|
||||
* \brief The Henry constants in liquid water using the IAPWS 2004 formulation.
|
||||
*
|
||||
* This function calculates \f$K_D\f$, see:
|
||||
*
|
||||
* IAPWS: "Guideline on the Henry's Constant and Vapor-Liquid
|
||||
* Distribution Constant for Gases in H2O and D2O at High
|
||||
* Temperatures"
|
||||
* http://www.iapws.org/relguide/HenGuide.pdf
|
||||
* IAPWS: "Guideline on the Henry's Constant and Vapor-Liquid Distribution Constant for
|
||||
* Gases in H2O and D2O at High Temperatures" http://www.iapws.org/relguide/HenGuide.pdf
|
||||
*/
|
||||
template <class Scalar>
|
||||
inline Scalar henryIAPWS(Scalar E,
|
||||
Scalar F,
|
||||
Scalar G,
|
||||
Scalar H,
|
||||
Scalar temperature)
|
||||
template <class Scalar, class Evaluation>
|
||||
inline Evaluation henryIAPWS(Scalar E,
|
||||
Scalar F,
|
||||
Scalar G,
|
||||
Scalar H,
|
||||
const Evaluation& temperature)
|
||||
{
|
||||
typedef Opm::H2O<Scalar> H2O;
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
typedef Opm::H2O<Evaluation> H2O;
|
||||
|
||||
Scalar Tr = temperature/H2O::criticalTemperature();
|
||||
Scalar tau = 1 - Tr;
|
||||
Evaluation Tr = temperature/H2O::criticalTemperature();
|
||||
Evaluation tau = 1 - Tr;
|
||||
|
||||
static const Scalar c[6] = {
|
||||
1.99274064, 1.09965342, -0.510839303,
|
||||
@ -63,22 +61,22 @@ inline Scalar henryIAPWS(Scalar E,
|
||||
};
|
||||
static const Scalar q = -0.023767;
|
||||
|
||||
Scalar f = 0;
|
||||
Evaluation f = 0;
|
||||
for (int i = 0; i < 6; ++i) {
|
||||
f += c[i]*std::pow(tau, d[i]);
|
||||
f += c[i]*Toolbox::pow(tau, d[i]);
|
||||
}
|
||||
|
||||
Scalar exponent =
|
||||
const Evaluation& exponent =
|
||||
q*F +
|
||||
E/temperature*f +
|
||||
(F +
|
||||
G*std::pow(tau, 2.0/3) +
|
||||
G*Toolbox::pow(tau, 2.0/3) +
|
||||
H*tau)*
|
||||
std::exp((H2O::tripleTemperature() - temperature)/100);
|
||||
Toolbox::exp((H2O::tripleTemperature() - temperature)/100);
|
||||
// CAUTION: K_D is formulated in mole fractions. We have to
|
||||
// multiply it with the vapor pressure of water in order to get
|
||||
// derivative of the partial pressure.
|
||||
return std::exp(exponent)*H2O::vaporPressure(temperature);
|
||||
return Toolbox::exp(exponent)*H2O::vaporPressure(temperature);
|
||||
}
|
||||
} // namespace Opm
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user