Merge pull request #3445 from svenn-t/h2store
Hydrogen-brine simulations - H2STORE
This commit is contained in:
Tor Harald Sandve
GitHub
@@ -51,6 +51,7 @@ list (APPEND MAIN_SOURCE_FILES
|
||||
src/opm/material/common/TridiagonalMatrix.cpp
|
||||
src/opm/material/common/UniformXTabulated2DFunction.cpp
|
||||
src/opm/material/components/CO2.cpp
|
||||
src/opm/material/components/H2.cpp
|
||||
src/opm/material/densead/Evaluation.cpp
|
||||
src/opm/material/fluidmatrixinteractions/EclEpsScalingPoints.cpp
|
||||
src/opm/material/fluidsystems/BlackOilFluidSystem.cpp
|
||||
@@ -291,6 +292,8 @@ if(ENABLE_ECL_INPUT)
|
||||
src/opm/material/fluidmatrixinteractions/EclMaterialLawManagerHystParams.cpp
|
||||
src/opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.cpp
|
||||
src/opm/material/fluidsystems/blackoilpvt/Co2GasPvt.cpp
|
||||
src/opm/material/fluidsystems/blackoilpvt/BrineH2Pvt.cpp
|
||||
src/opm/material/fluidsystems/blackoilpvt/H2GasPvt.cpp
|
||||
src/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityBrinePvt.cpp
|
||||
src/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.cpp
|
||||
src/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.cpp
|
||||
@@ -816,6 +819,7 @@ list( APPEND PUBLIC_HEADER_FILES
|
||||
opm/material/binarycoefficients/H2O_CO2.hpp
|
||||
opm/material/binarycoefficients/Air_Xylene.hpp
|
||||
opm/material/binarycoefficients/Brine_CO2.hpp
|
||||
opm/material/binarycoefficients/Brine_H2.hpp
|
||||
opm/material/binarycoefficients/HenryIapws.hpp
|
||||
opm/material/Constants.hpp
|
||||
opm/material/fluidsystems/NullParameterCache.hpp
|
||||
@@ -839,6 +843,7 @@ list( APPEND PUBLIC_HEADER_FILES
|
||||
opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp
|
||||
opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.hpp
|
||||
opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.hpp
|
||||
opm/material/fluidsystems/blackoilpvt/BrineH2Pvt.hpp
|
||||
opm/material/fluidsystems/blackoilpvt/OilPvtMultiplexer.hpp
|
||||
opm/material/fluidsystems/blackoilpvt/GasPvtMultiplexer.hpp
|
||||
opm/material/fluidsystems/blackoilpvt/DryHumidGasPvt.hpp
|
||||
@@ -851,6 +856,7 @@ list( APPEND PUBLIC_HEADER_FILES
|
||||
opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityBrinePvt.hpp
|
||||
opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp
|
||||
opm/material/fluidsystems/blackoilpvt/Co2GasPvt.hpp
|
||||
opm/material/fluidsystems/blackoilpvt/H2GasPvt.hpp
|
||||
opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp
|
||||
opm/material/fluidsystems/H2OAirFluidSystem.hpp
|
||||
opm/material/fluidsystems/H2ON2FluidSystem.hpp
|
||||
|
||||
@@ -455,6 +455,7 @@ public:
|
||||
const Nupcol& nupcol() const noexcept;
|
||||
const Tracers& tracers() const;
|
||||
bool co2Storage() const noexcept;
|
||||
bool h2Storage() const noexcept;
|
||||
bool micp() const noexcept;
|
||||
|
||||
bool operator==(const Runspec& data) const;
|
||||
@@ -478,6 +479,7 @@ public:
|
||||
serializer(m_sfuncctrl);
|
||||
serializer(m_nupcol);
|
||||
serializer(m_co2storage);
|
||||
serializer(m_h2storage);
|
||||
serializer(m_micp);
|
||||
}
|
||||
|
||||
@@ -498,6 +500,7 @@ private:
|
||||
Nupcol m_nupcol;
|
||||
Tracers m_tracers;
|
||||
bool m_co2storage;
|
||||
bool m_h2storage;
|
||||
bool m_micp;
|
||||
};
|
||||
|
||||
|
||||
196
opm/material/binarycoefficients/Brine_H2.hpp
Normal file
196
opm/material/binarycoefficients/Brine_H2.hpp
Normal file
@@ -0,0 +1,196 @@
|
||||
// -*- 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
|
||||
*
|
||||
* \copydoc Opm::BinaryCoeff::Brine_H2
|
||||
*/
|
||||
#ifndef OPM_BINARY_COEFF_BRINE_H2_HPP
|
||||
#define OPM_BINARY_COEFF_BRINE_H2_HPP
|
||||
|
||||
#include <opm/material/binarycoefficients/FullerMethod.hpp>
|
||||
|
||||
namespace Opm {
|
||||
namespace BinaryCoeff {
|
||||
|
||||
/*!
|
||||
* \ingroup Binarycoefficients
|
||||
* \brief Binary coefficients for brine and H2.
|
||||
*/
|
||||
template<class Scalar, class H2O, class H2, bool verbose = true>
|
||||
class Brine_H2 {
|
||||
static const int liquidPhaseIdx = 0; // index of the liquid phase
|
||||
static const int gasPhaseIdx = 1; // index of the gas phase
|
||||
|
||||
public:
|
||||
/*!
|
||||
* \brief Returns the _mol_ (!) fraction of H2 in the liquid phase for a given temperature, pressure, and brine
|
||||
* salinity. Implemented according to Li et al., Int. J. Hydrogen Energ., 2018.
|
||||
*
|
||||
* \param temperature temperature [K]
|
||||
* \param pg gas phase pressure [Pa]
|
||||
* \param salinity salinity [mol NaCl / kg solution]
|
||||
* \param knownPhaseIdx indicates which phases are present
|
||||
* \param xlH2 mole fraction of H2 in brine [mol/mol]
|
||||
*/
|
||||
template <class Evaluation>
|
||||
static void calculateMoleFractions(const Evaluation& temperature,
|
||||
const Evaluation& pg,
|
||||
Scalar salinity,
|
||||
Evaluation& xH2,
|
||||
bool extrapolate = false)
|
||||
{
|
||||
// All intermediate calculations
|
||||
Evaluation lnYH2 = moleFractionGasH2_(temperature, pg);
|
||||
Evaluation lnPg = log(pg / 1e6); // Pa --> MPa before ln
|
||||
Evaluation lnPhiH2 = fugacityCoefficientH2(temperature, pg, extrapolate);
|
||||
Evaluation lnKh = henrysConstant_(temperature);
|
||||
Evaluation PF = computePoyntingFactor_(temperature, pg);
|
||||
Evaluation lnGammaH2 = activityCoefficient_(temperature, salinity);
|
||||
|
||||
// Eq. (4) to get mole fraction of H2 in brine
|
||||
xH2 = exp(lnYH2 + lnPg + lnPhiH2 - lnKh - PF - lnGammaH2);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the Poynting Factor (PF) which is needed in calculation of H2 solubility in Li et al (2018).
|
||||
*
|
||||
* \param temperature temperature [K]
|
||||
* \param pg gas phase pressure [Pa]
|
||||
*/
|
||||
template <class Evaluation>
|
||||
static Evaluation computePoyntingFactor_(const Evaluation& temperature, const Evaluation& pg)
|
||||
{
|
||||
// PF is approximated as a polynomial expansion in terms of temperature and pressure with the following
|
||||
// parameters (Table 4)
|
||||
static const Scalar a[4] = {6.156755, -2.502396e-2, 4.140593e-5, -1.322988e-3};
|
||||
|
||||
// Eq. (16)
|
||||
Evaluation pg_mpa = pg / 1.0e6; // convert from Pa to MPa
|
||||
Evaluation PF = a[0]*pg_mpa/temperature + a[1]*pg_mpa + a[2]*temperature*pg_mpa + a[3]*pg_mpa*pg_mpa/temperature;
|
||||
return PF;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the activity coefficient of H2 in brine which is needed in calculation of H2 solubility in Li et
|
||||
* al (2018). Note that we only include NaCl effects. Could be extended with other salts, e.g. from Duan & Sun,
|
||||
* Chem. Geol., 2003.
|
||||
*
|
||||
* \param temperature temperature [K]
|
||||
* \param salinity salinity [mol NaCl / kg solution]
|
||||
*/
|
||||
template <class Evaluation>
|
||||
static Evaluation activityCoefficient_(const Evaluation& temperature, Scalar salinity)
|
||||
{
|
||||
// Linear approximation in temperature with following parameters (Table 5)
|
||||
static const Scalar a[2] = {0.64485, 0.00142};
|
||||
|
||||
// Eq. (17)
|
||||
Evaluation lnGamma = (a[0] - a[1]*temperature)*salinity;
|
||||
return lnGamma;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns Henry's constant of H2 in brine which is needed in calculation of H2 solubility in Li et al (2018).
|
||||
*
|
||||
* \param temperature temperature [K]
|
||||
*/
|
||||
template <class Evaluation>
|
||||
static Evaluation henrysConstant_(const Evaluation& temperature)
|
||||
{
|
||||
// Polynomic approximation in temperature with following parameters (Table 2)
|
||||
static const Scalar a[5] = {2.68721e-5, -0.05121, 33.55196, -3411.0432, -31258.74683};
|
||||
|
||||
// Eq. (13)
|
||||
Evaluation lnKh = a[0]*temperature*temperature + a[1]*temperature + a[2] + a[3]/temperature
|
||||
+ a[4]/(temperature*temperature);
|
||||
return lnKh;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns mole fraction of H2 in gasous phase which is needed in calculation of H2 solubility in Li et al
|
||||
* (2018).
|
||||
*
|
||||
* \param temperature temperature [K]
|
||||
* \param pg gas phase pressure [Pa]
|
||||
*/
|
||||
template <class Evaluation>
|
||||
static Evaluation moleFractionGasH2_(const Evaluation& temperature, const Evaluation& pg)
|
||||
{
|
||||
// Need saturaturated vapor pressure of pure water
|
||||
Evaluation pw_sat = H2O::vaporPressure(temperature);
|
||||
|
||||
// Eq. (12)
|
||||
Evaluation lnyH2 = log(1 - (pw_sat / pg));
|
||||
return lnyH2;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \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 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]
|
||||
*/
|
||||
template <class Evaluation>
|
||||
static Evaluation fugacityCoefficientH2(const Evaluation& temperature,
|
||||
const Evaluation& pg,
|
||||
bool extrapolate = false)
|
||||
{
|
||||
// Convert pressure to reduced density and temperature to reduced temperature
|
||||
Evaluation rho_red = H2::reducedMolarDensity(temperature, pg, extrapolate);
|
||||
Evaluation T_red = H2::criticalTemperature() / temperature;
|
||||
|
||||
// Residual Helmholtz energy, Eq. (7) in Li et al. (2018)
|
||||
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 dResHelm_dRedRho = H2::derivResHelmholtzWrtRedRho(T_red, rho_red);
|
||||
|
||||
// Fugacity coefficient, Eq. (8) in Li et al. (2018)
|
||||
Evaluation lnPhiH2 = resHelm + rho_red * dResHelm_dRedRho - log(rho_red * dResHelm_dRedRho + 1);
|
||||
|
||||
return lnPhiH2;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Binary diffusion coefficent [m^2/s] for molecular water and H2 as an approximation for brine-H2 diffusion.
|
||||
*
|
||||
* To calculate the values, the \ref fullerMethod is used.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
static Evaluation gasDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
|
||||
{
|
||||
// atomic diffusion volumes
|
||||
const Scalar SigmaNu[2] = { 13.1 /* H2O */, 7.07 /* H2 */ };
|
||||
// molar masses [g/mol]
|
||||
const Scalar M[2] = { H2O::molarMass()*Scalar(1e3), H2::molarMass()*Scalar(1e3) };
|
||||
|
||||
return fullerMethod(M, SigmaNu, temperature, pressure);
|
||||
}
|
||||
}; // end class Brine_H2
|
||||
|
||||
} // end namespace BinaryCoeff
|
||||
} // end namespace Opm
|
||||
#endif
|
||||
@@ -21,58 +21,93 @@
|
||||
copyright holders.
|
||||
*/
|
||||
/*!
|
||||
* \file
|
||||
* \copydoc Opm::H2
|
||||
*/
|
||||
* \file
|
||||
*
|
||||
* \ingroup Components
|
||||
*
|
||||
* \copydoc Opm::H2
|
||||
*
|
||||
*/
|
||||
#ifndef OPM_H2_HPP
|
||||
#define OPM_H2_HPP
|
||||
|
||||
#include "Component.hpp"
|
||||
|
||||
#include <opm/material/IdealGas.hpp>
|
||||
#include <opm/material/common/MathToolbox.hpp>
|
||||
#include <opm/material/components/Component.hpp>
|
||||
#include <opm/material/common/UniformTabulated2DFunction.hpp>
|
||||
#include <opm/material/densead/Math.hpp>
|
||||
|
||||
#include <cmath>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
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>
|
||||
class H2 : public Component<Scalar, H2<Scalar> >
|
||||
{
|
||||
typedef ::Opm::IdealGas<Scalar> IdealGas;
|
||||
using IdealGas = Opm::IdealGas<Scalar>;
|
||||
|
||||
static const UniformTabulated2DFunction<double>& tabulatedDensity;
|
||||
|
||||
public:
|
||||
/*!
|
||||
* \brief A human readable name for hydrogen.
|
||||
*/
|
||||
static const char* name()
|
||||
* \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 Scalar molarMass()
|
||||
{ return 0.0020156;}
|
||||
* \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 1.297e6; /* [N/m^2] */ }
|
||||
{ return 1.2964e6; /* [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 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 Critical volume of \f$H_2\f$ [m2/kmol].
|
||||
@@ -84,9 +119,567 @@ public:
|
||||
*/
|
||||
static Scalar acentricFactor() { return -0.22; }
|
||||
|
||||
/*!
|
||||
* \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$
|
||||
*
|
||||
*/
|
||||
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
|
||||
|
||||
// Intermediate calculations involving temperature
|
||||
Evaluation sigma = 1 - temperature/criticalTemperature();
|
||||
Evaluation T_recp = criticalTemperature() / temperature;
|
||||
|
||||
// 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();
|
||||
}
|
||||
|
||||
/*!
|
||||
* \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, bool extrapolate = false)
|
||||
{
|
||||
return tabulatedDensity.eval(temperature, pressure, extrapolate);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \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, bool extrapolate = false)
|
||||
{ return gasDensity(temperature, pressure, extrapolate) / molarMass(); }
|
||||
|
||||
/*!
|
||||
* \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 false; }
|
||||
|
||||
/*!
|
||||
* \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)
|
||||
{
|
||||
// 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));
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Specific internal energy of H2 [J/kg].
|
||||
*/
|
||||
template <class Evaluation>
|
||||
static Evaluation gasInternalEnergy(const Evaluation& temperature,
|
||||
const Evaluation& pressure,
|
||||
bool extrapolate = false)
|
||||
{
|
||||
// Eq. (58) in Span et al. (2000)
|
||||
Evaluation rho_red = reducedMolarDensity(temperature, pressure, extrapolate);
|
||||
Evaluation T_red = criticalTemperature() / temperature;
|
||||
Scalar R = IdealGas::R;
|
||||
|
||||
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,
|
||||
bool extrapolate = false)
|
||||
{
|
||||
// Eq. (59) in Span et al. (2000)
|
||||
Evaluation u = gasInternalEnergy(temperature, pressure);
|
||||
Evaluation rho = gasDensity(temperature, pressure, extrapolate);
|
||||
|
||||
return u + (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 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 gasHeatCapacity(Evaluation temperature,
|
||||
Evaluation pressure)
|
||||
{
|
||||
// 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 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 gasIsochoricHeatCapacity(Evaluation temperature,
|
||||
Evaluation pressure)
|
||||
{
|
||||
// Reduced variables
|
||||
Evaluation rho_red = reducedMolarDensity(temperature, pressure);
|
||||
Evaluation T_red = criticalTemperature() / temperature;
|
||||
|
||||
// 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,
|
||||
bool extrapolate = false)
|
||||
{
|
||||
return gasDensity(temperature, pg, extrapolate) / (molarMass() * criticalDensity() * 1e6);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \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;
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
} // end namespace Opm
|
||||
|
||||
#endif
|
||||
|
||||
275
opm/material/components/SimpleH2.hpp
Normal file
275
opm/material/components/SimpleH2.hpp
Normal file
@@ -0,0 +1,275 @@
|
||||
// -*- 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$. Uses 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 Critical volume of \f$H_2\f$ [m2/kmol].
|
||||
*/
|
||||
static Scalar criticalVolume() {return 6.45e-2; }
|
||||
|
||||
/*!
|
||||
* \brief Acentric factor of \f$H_2\f$.
|
||||
*/
|
||||
static Scalar acentricFactor() { return -0.22; }
|
||||
|
||||
/*!
|
||||
* \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
|
||||
@@ -32,6 +32,7 @@
|
||||
#include "blackoilpvt/GasPvtMultiplexer.hpp"
|
||||
#include "blackoilpvt/WaterPvtMultiplexer.hpp"
|
||||
#include "blackoilpvt/BrineCo2Pvt.hpp"
|
||||
#include "blackoilpvt/BrineH2Pvt.hpp"
|
||||
|
||||
#include <opm/common/TimingMacros.hpp>
|
||||
|
||||
|
||||
556
opm/material/fluidsystems/blackoilpvt/BrineH2Pvt.hpp
Normal file
556
opm/material/fluidsystems/blackoilpvt/BrineH2Pvt.hpp
Normal file
@@ -0,0 +1,556 @@
|
||||
// -*- 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
|
||||
* \copydoc Opm::BrineH2Pvt
|
||||
*/
|
||||
#ifndef OPM_BRINE_H2_PVT_HPP
|
||||
#define OPM_BRINE_H2_PVT_HPP
|
||||
|
||||
#include <opm/common/Exceptions.hpp>
|
||||
|
||||
#include <opm/material/binarycoefficients/Brine_H2.hpp>
|
||||
#include <opm/material/components/SimpleHuDuanH2O.hpp>
|
||||
#include <opm/material/components/Brine.hpp>
|
||||
#include <opm/material/components/H2.hpp>
|
||||
#include <opm/material/common/UniformTabulated2DFunction.hpp>
|
||||
|
||||
#include <vector>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
#if HAVE_ECL_INPUT
|
||||
class EclipseState;
|
||||
class Schedule;
|
||||
#endif
|
||||
|
||||
/*!
|
||||
* \brief This class represents the Pressure-Volume-Temperature relations of the liquid phase for a H2-Brine system
|
||||
*/
|
||||
template <class Scalar>
|
||||
class BrineH2Pvt
|
||||
{
|
||||
static const bool extrapolate = true;
|
||||
public:
|
||||
using H2O = SimpleHuDuanH2O<Scalar>;
|
||||
using Brine = ::Opm::Brine<Scalar, H2O>;
|
||||
using H2 = ::Opm::H2<Scalar>;
|
||||
|
||||
// The binary coefficients for brine and H2 used by this fluid system
|
||||
using BinaryCoeffBrineH2 = BinaryCoeff::Brine_H2<Scalar, H2O, H2>;
|
||||
|
||||
explicit BrineH2Pvt() = default;
|
||||
|
||||
BrineH2Pvt(const std::vector<Scalar>& salinity,
|
||||
Scalar T_ref = 288.71, //(273.15 + 15.56)
|
||||
Scalar P_ref = 101325)
|
||||
: salinity_(salinity)
|
||||
{
|
||||
int num_regions = salinity_.size();
|
||||
h2ReferenceDensity_.resize(num_regions);
|
||||
brineReferenceDensity_.resize(num_regions);
|
||||
Brine::salinity = salinity[0];
|
||||
for (int i = 0; i < num_regions; ++i) {
|
||||
h2ReferenceDensity_[i] = H2::gasDensity(T_ref, P_ref, true);
|
||||
brineReferenceDensity_[i] = Brine::liquidDensity(T_ref, P_ref, true);
|
||||
}
|
||||
}
|
||||
|
||||
#if HAVE_ECL_INPUT
|
||||
/*!
|
||||
* \brief Initialize the parameters for Brine-H2 system using an ECL deck.
|
||||
*
|
||||
*/
|
||||
void initFromState(const EclipseState& eclState, const Schedule&);
|
||||
#endif
|
||||
|
||||
void setNumRegions(size_t numRegions)
|
||||
{
|
||||
brineReferenceDensity_.resize(numRegions);
|
||||
h2ReferenceDensity_.resize(numRegions);
|
||||
salinity_.resize(numRegions);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Initialize the reference densities of all fluids for a given PVT region
|
||||
*/
|
||||
void setReferenceDensities(unsigned regionIdx,
|
||||
Scalar rhoRefBrine,
|
||||
Scalar rhoRefH2,
|
||||
Scalar /*rhoRefWater*/)
|
||||
{
|
||||
brineReferenceDensity_[regionIdx] = rhoRefBrine;
|
||||
h2ReferenceDensity_[regionIdx] = rhoRefH2;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Finish initializing the oil phase PVT properties.
|
||||
*/
|
||||
void initEnd()
|
||||
{
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Specify whether the PVT model should consider that the H2 component can
|
||||
* dissolve in the brine phase
|
||||
*
|
||||
* By default, dissolved H2 is considered.
|
||||
*/
|
||||
void setEnableDissolvedGas(bool yesno)
|
||||
{ enableDissolution_ = yesno; }
|
||||
|
||||
/*!
|
||||
* \brief Return the number of PVT regions which are considered by this PVT-object.
|
||||
*/
|
||||
unsigned numRegions() const
|
||||
{ return brineReferenceDensity_.size(); }
|
||||
|
||||
/*!
|
||||
* \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation internalEnergy(unsigned /*regionIdx*/,
|
||||
const Evaluation& /*temperature*/,
|
||||
const Evaluation& /*pressure*/,
|
||||
const Evaluation& /*Rs*/,
|
||||
const Evaluation& /*saltConcentration*/) const
|
||||
{
|
||||
throw std::runtime_error("Requested internal energy for the H2-brine PVT module. Not yet implemented.");
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation internalEnergy(unsigned /*regionIdx*/,
|
||||
const Evaluation& /*temperature*/,
|
||||
const Evaluation& /*pressure*/,
|
||||
const Evaluation& /*Rs*/) const
|
||||
{
|
||||
throw std::runtime_error("Requested internal energy for the H2-brine PVT module. Not yet implemented.");
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation viscosity(unsigned regionIdx,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure,
|
||||
const Evaluation& /*Rs*/) const
|
||||
{
|
||||
//TODO: The viscosity does not yet depend on the composition
|
||||
return saturatedViscosity(regionIdx, temperature, pressure);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturatedViscosity(unsigned regionIdx,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure,
|
||||
const Evaluation& /*saltConcentration*/) const
|
||||
{
|
||||
return saturatedViscosity(regionIdx, temperature, pressure);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation viscosity(unsigned regionIdx,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure,
|
||||
const Evaluation& /*Rsw*/,
|
||||
const Evaluation& /*saltConcentration*/) const
|
||||
{
|
||||
//TODO: The viscosity does not yet depend on the composition
|
||||
return saturatedViscosity(regionIdx, temperature, pressure);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the dynamic viscosity [Pa s] of oil saturated gas at given pressure.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturatedViscosity(unsigned /*regionIdx*/,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure) const
|
||||
{
|
||||
return Brine::liquidViscosity(temperature, pressure);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the formation volume factor [-] of the fluid phase.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure,
|
||||
const Evaluation& /*saltconcentration*/) const
|
||||
{
|
||||
return saturatedInverseFormationVolumeFactor(regionIdx, temperature, pressure);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the formation volume factor [-] of the fluid phase.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure,
|
||||
const Evaluation& Rs,
|
||||
const Evaluation& /*saltConcentration*/) const
|
||||
{
|
||||
return inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rs);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the formation volume factor [-] of the fluid phase.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure,
|
||||
const Evaluation& Rs) const
|
||||
{
|
||||
return (1.0 - convertRsToXoG_(Rs, regionIdx)) * density_(regionIdx, temperature, pressure, Rs) /
|
||||
brineReferenceDensity_[regionIdx];
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the formation volume factor [-] of brine saturated with H2 at a given pressure.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure) const
|
||||
{
|
||||
Evaluation rsSat = rsSat_(regionIdx, temperature, pressure);
|
||||
return (1.0 - convertRsToXoG_(rsSat, regionIdx)) * density_(regionIdx, temperature, pressure, rsSat) /
|
||||
brineReferenceDensity_[regionIdx];
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the saturation pressure of the brine phase [Pa] depending on its mass fraction of the gas component
|
||||
*
|
||||
* \param Rs
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturationPressure(unsigned /*regionIdx*/,
|
||||
const Evaluation& /*temperature*/,
|
||||
const Evaluation& /*Rs*/) const
|
||||
{
|
||||
throw std::runtime_error("Saturation pressure for the Brine-H2 PVT module has not been implemented yet!");
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the saturation pressure of the brine phase [Pa] depending on its mass fraction of the gas component
|
||||
*
|
||||
* \param Rs
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturationPressure(unsigned /*regionIdx*/,
|
||||
const Evaluation& /*temperature*/,
|
||||
const Evaluation& /*Rs*/,
|
||||
const Evaluation& /*saltConcentration*/) const
|
||||
{
|
||||
throw std::runtime_error("Saturation pressure for the Brine-H2 PVT module has not been implemented yet!");
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the gas dissoluiton factor \f$R_s\f$ [m^3/m^3] of the liquid phase.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturatedGasDissolutionFactor(unsigned regionIdx,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure,
|
||||
const Evaluation& /*oilSaturation*/,
|
||||
const Evaluation& /*maxOilSaturation*/) const
|
||||
{
|
||||
//TODO support VAPPARS
|
||||
return rsSat_(regionIdx, temperature, pressure);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the gas dissoluiton factor \f$R_s\f$ [m^3/m^3] of the liquid phase.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturatedGasDissolutionFactor(unsigned regionIdx,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure,
|
||||
const Evaluation& /*saltConcentration*/) const
|
||||
{
|
||||
return rsSat_(regionIdx, temperature, pressure);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns thegas dissoluiton factor \f$R_s\f$ [m^3/m^3] of the liquid phase.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturatedGasDissolutionFactor(unsigned regionIdx,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure) const
|
||||
{
|
||||
return rsSat_(regionIdx, temperature, pressure);
|
||||
}
|
||||
|
||||
const Scalar oilReferenceDensity(unsigned regionIdx) const
|
||||
{ return brineReferenceDensity_[regionIdx]; }
|
||||
|
||||
const Scalar waterReferenceDensity(unsigned regionIdx) const
|
||||
{ return brineReferenceDensity_[regionIdx]; }
|
||||
|
||||
const Scalar gasReferenceDensity(unsigned regionIdx) const
|
||||
{ return h2ReferenceDensity_[regionIdx]; }
|
||||
|
||||
const Scalar salinity(unsigned regionIdx) const
|
||||
{ return salinity_[regionIdx]; }
|
||||
|
||||
/*!
|
||||
* \brief Diffusion coefficient of H2 in water
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation diffusionCoefficient(const Evaluation& temperature,
|
||||
const Evaluation& pressure,
|
||||
unsigned /*compIdx*/) const
|
||||
{
|
||||
// Diffusion coefficient of H2 in pure water according to Ferrell & Himmelbau, AIChE Journal, 13(4), 1967 (Eq.
|
||||
// 23)
|
||||
// Some intermediate calculations and definitions
|
||||
const Scalar vm = 28.45; // molar volume at normal boiling point (20.271 K and 1 atm) in cm2/mol
|
||||
const Scalar sigma = 2.96 * 1e-8; // Lennard-Jones 6-12 potential in cm (1 Å = 1e-8 cm)
|
||||
const Scalar avogadro = 6.022e23; // Avogrado's number in mol^-1
|
||||
const Scalar alpha = sigma / pow((vm / avogadro), 1 / 3); // Eq. (19)
|
||||
const Scalar lambda = 1.729; // quantum parameter [-]
|
||||
const Evaluation& mu_pure = H2O::liquidViscosity(temperature, pressure, extrapolate) * 1e3; // water viscosity in cP
|
||||
|
||||
// Diffusion coeff in pure water in cm2/s
|
||||
const Evaluation D_pure = ((4.8e-7 * temperature) / pow(mu_pure, alpha)) * pow((1 + pow(lambda, 2)) / vm, 0.6);
|
||||
|
||||
// Diffusion coefficient in brine using Ratcliff and Holdcroft, J. G. Trans. Inst. Chem. Eng, 1963. OBS: Value
|
||||
// of n is noted as the recommended single value according to Akita, Ind. Eng. Chem. Fundam., 1981.
|
||||
const Evaluation& mu_brine = Brine::liquidViscosity(temperature, pressure) * 1e3; // Brine viscosity in cP
|
||||
const Evaluation log_D_brine = log10(D_pure) - 0.637 * log10(mu_brine / mu_pure);
|
||||
|
||||
return pow(Evaluation(10), log_D_brine) * 1e-4; // convert from cm2/s to m2/s
|
||||
}
|
||||
|
||||
private:
|
||||
std::vector<Scalar> brineReferenceDensity_;
|
||||
std::vector<Scalar> h2ReferenceDensity_;
|
||||
std::vector<Scalar> salinity_;
|
||||
bool enableDissolution_ = true;
|
||||
|
||||
/*!
|
||||
* \brief Calculate density of aqueous solution (H2O-NaCl/brine and H2).
|
||||
*
|
||||
* \param temperature temperature [K]
|
||||
* \param pressure pressure [Pa]
|
||||
* \param Rs gas dissolution factor [-]
|
||||
*/
|
||||
template <class LhsEval>
|
||||
LhsEval density_(unsigned regionIdx,
|
||||
const LhsEval& temperature,
|
||||
const LhsEval& pressure,
|
||||
const LhsEval& Rs) const
|
||||
{
|
||||
// convert Rs to mole fraction (via mass fraction)
|
||||
LhsEval xlH2 = convertXoGToxoG_(convertRsToXoG_(Rs,regionIdx));
|
||||
|
||||
// calculate the density of solution
|
||||
LhsEval result = liquidDensity_(temperature,
|
||||
pressure,
|
||||
xlH2);
|
||||
|
||||
Valgrind::CheckDefined(result);
|
||||
return result;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Calculated the density of the aqueous solution where contributions of salinity and dissolved H2 is taken
|
||||
* into account.
|
||||
*
|
||||
* \param T temperature [K]
|
||||
* \param pl liquid pressure [Pa]
|
||||
* \param xlH2 mole fraction H2 [-]
|
||||
*/
|
||||
template <class LhsEval>
|
||||
LhsEval liquidDensity_(const LhsEval& T,
|
||||
const LhsEval& pl,
|
||||
const LhsEval& xlH2) const
|
||||
{
|
||||
// check input variables
|
||||
Valgrind::CheckDefined(T);
|
||||
Valgrind::CheckDefined(pl);
|
||||
Valgrind::CheckDefined(xlH2);
|
||||
|
||||
// check if pressure and temperature is valid
|
||||
if(!extrapolate && T < 273.15) {
|
||||
const std::string msg =
|
||||
"Liquid density for Brine and H2 is only "
|
||||
"defined above 273.15K (is " +
|
||||
std::to_string(getValue(T)) + "K)";
|
||||
throw NumericalProblem(msg);
|
||||
}
|
||||
if(!extrapolate && pl >= 2.5e8) {
|
||||
const std::string msg =
|
||||
"Liquid density for Brine and H2 is only "
|
||||
"defined below 250MPa (is " +
|
||||
std::to_string(getValue(pl)) + "Pa)";
|
||||
throw NumericalProblem(msg);
|
||||
}
|
||||
|
||||
// calculate individual contribution to density
|
||||
const LhsEval& rho_brine = Brine::liquidDensity(T, pl, extrapolate);
|
||||
const LhsEval& rho_pure = H2O::liquidDensity(T, pl, extrapolate);
|
||||
const LhsEval& rho_lH2 = liquidDensityWaterH2_(T, pl, xlH2);
|
||||
const LhsEval& contribH2 = rho_lH2 - rho_pure;
|
||||
|
||||
return rho_brine + contribH2;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Density of aqueous solution with dissolved H2. Formula from Li et al. (2018) and Garica, Lawrence Berkeley
|
||||
* National Laboratory, 2001.
|
||||
*
|
||||
* \param temperature [K]
|
||||
* \param pl liquid pressure [Pa]
|
||||
* \param xlH2 mole fraction [-]
|
||||
*/
|
||||
template <class LhsEval>
|
||||
LhsEval liquidDensityWaterH2_(const LhsEval& temperature,
|
||||
const LhsEval& pl,
|
||||
const LhsEval& xlH2) const
|
||||
{
|
||||
// molar masses
|
||||
Scalar M_H2 = H2::molarMass();
|
||||
Scalar M_H2O = H2O::molarMass();
|
||||
|
||||
// density of pure water
|
||||
const LhsEval& rho_pure = H2O::liquidDensity(temperature, pl, extrapolate);
|
||||
|
||||
// (apparent) molar volume of H2, Eq. (14) in Li et al. (2018)
|
||||
const LhsEval& A1 = 51.1904 - 0.208062*temperature + 3.4427e-4*(temperature*temperature);
|
||||
const LhsEval& A2 = -0.022;
|
||||
const LhsEval& V_phi = (A1 + A2 * (pl / 1e6)) / 1e6; // pressure in [MPa] and Vphi in [m3/mol] (from [cm3/mol])
|
||||
|
||||
// density of solution, Eq. (19) in Garcia (2001)
|
||||
const LhsEval xlH2O = 1.0 - xlH2;
|
||||
const LhsEval& M_T = M_H2O * xlH2O + M_H2 * xlH2;
|
||||
const LhsEval& rho_aq = 1 / (xlH2 * V_phi/M_T + M_H2O * xlH2O / (rho_pure * M_T));
|
||||
|
||||
return rho_aq;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Convert a gas dissolution factor to the the corresponding mass fraction of the gas component in the oil
|
||||
* phase.
|
||||
*
|
||||
* \param Rs gass dissolution factor [-]
|
||||
* \param regionIdx region index
|
||||
*/
|
||||
template <class LhsEval>
|
||||
LhsEval convertRsToXoG_(const LhsEval& Rs, unsigned regionIdx) const
|
||||
{
|
||||
Scalar rho_oRef = brineReferenceDensity_[regionIdx];
|
||||
Scalar rho_gRef = h2ReferenceDensity_[regionIdx];
|
||||
|
||||
const LhsEval& rho_oG = Rs*rho_gRef;
|
||||
return rho_oG/(rho_oRef + rho_oG);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Convert a gas mass fraction in the oil phase the corresponding mole fraction.
|
||||
*
|
||||
* \param XoG mass fraction [-]
|
||||
*/
|
||||
template <class LhsEval>
|
||||
LhsEval convertXoGToxoG_(const LhsEval& XoG) const
|
||||
{
|
||||
Scalar M_H2 = H2::molarMass();
|
||||
Scalar M_Brine = Brine::molarMass();
|
||||
return XoG*M_Brine / (M_H2*(1 - XoG) + XoG*M_Brine);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Convert a gas mole fraction in the oil phase the corresponding mass fraction.
|
||||
*
|
||||
* \param xoG mole fraction [-]
|
||||
*/
|
||||
template <class LhsEval>
|
||||
LhsEval convertxoGToXoG(const LhsEval& xoG) const
|
||||
{
|
||||
Scalar M_H2 = H2::molarMass();
|
||||
Scalar M_Brine = Brine::molarMass();
|
||||
|
||||
return xoG*M_H2 / (xoG*(M_H2 - M_Brine) + M_Brine);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Convert the mass fraction of the gas component in the oil phase to the corresponding gas dissolution
|
||||
* factor.
|
||||
*
|
||||
* \param XoG mass fraction [-]
|
||||
* \param regionIdx region index
|
||||
*/
|
||||
template <class LhsEval>
|
||||
LhsEval convertXoGToRs(const LhsEval& XoG, unsigned regionIdx) const
|
||||
{
|
||||
Scalar rho_oRef = brineReferenceDensity_[regionIdx];
|
||||
Scalar rho_gRef = h2ReferenceDensity_[regionIdx];
|
||||
|
||||
return XoG/(1.0 - XoG)*(rho_oRef/rho_gRef);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Saturated gas dissolution factor, Rs.
|
||||
*
|
||||
* \param regionIdx region index
|
||||
* \param temperature [K]
|
||||
* \param pressure pressure [Pa]
|
||||
*/
|
||||
template <class LhsEval>
|
||||
LhsEval rsSat_(unsigned regionIdx,
|
||||
const LhsEval& temperature,
|
||||
const LhsEval& pressure) const
|
||||
{
|
||||
// Return Rs=0.0 if dissolution is disabled
|
||||
if (!enableDissolution_)
|
||||
return 0.0;
|
||||
// calulate the equilibrium composition for the given temperature and pressure
|
||||
LhsEval xlH2;
|
||||
BinaryCoeffBrineH2::calculateMoleFractions(temperature,
|
||||
pressure,
|
||||
salinity_[regionIdx],
|
||||
xlH2,
|
||||
extrapolate);
|
||||
|
||||
// normalize the phase compositions
|
||||
xlH2 = max(0.0, min(1.0, xlH2));
|
||||
|
||||
return convertXoGToRs(convertxoGToXoG(xlH2), regionIdx);
|
||||
}
|
||||
}; // end class BrineH2Pvt
|
||||
} // end namespace Opm
|
||||
#endif
|
||||
@@ -33,6 +33,7 @@
|
||||
#include "WetGasPvt.hpp"
|
||||
#include "GasPvtThermal.hpp"
|
||||
#include "Co2GasPvt.hpp"
|
||||
#include "H2GasPvt.hpp"
|
||||
|
||||
namespace Opm {
|
||||
|
||||
@@ -73,6 +74,11 @@ class Schedule;
|
||||
codeToCall; \
|
||||
break; \
|
||||
} \
|
||||
case GasPvtApproach::H2Gas: { \
|
||||
auto& pvtImpl = getRealPvt<GasPvtApproach::H2Gas>(); \
|
||||
codeToCall; \
|
||||
break; \
|
||||
} \
|
||||
case GasPvtApproach::NoGas: \
|
||||
throw std::logic_error("Not implemented: Gas PVT of this deck!"); \
|
||||
}
|
||||
@@ -84,7 +90,8 @@ enum class GasPvtApproach {
|
||||
WetHumidGas,
|
||||
WetGas,
|
||||
ThermalGas,
|
||||
Co2Gas
|
||||
Co2Gas,
|
||||
H2Gas
|
||||
};
|
||||
|
||||
/*!
|
||||
@@ -144,6 +151,10 @@ public:
|
||||
delete &getRealPvt<GasPvtApproach::Co2Gas>();
|
||||
break;
|
||||
}
|
||||
case GasPvtApproach::H2Gas: {
|
||||
delete &getRealPvt<GasPvtApproach::H2Gas>();
|
||||
break;
|
||||
}
|
||||
case GasPvtApproach::NoGas:
|
||||
break;
|
||||
}
|
||||
@@ -185,6 +196,10 @@ public:
|
||||
realGasPvt_ = new Co2GasPvt<Scalar>;
|
||||
break;
|
||||
|
||||
case GasPvtApproach::H2Gas:
|
||||
realGasPvt_ = new H2GasPvt<Scalar>;
|
||||
break;
|
||||
|
||||
case GasPvtApproach::NoGas:
|
||||
throw std::logic_error("Not implemented: Gas PVT of this deck!");
|
||||
}
|
||||
@@ -416,6 +431,20 @@ public:
|
||||
return *static_cast<const Co2GasPvt<Scalar>* >(realGasPvt_);
|
||||
}
|
||||
|
||||
template <GasPvtApproach approachV>
|
||||
typename std::enable_if<approachV == GasPvtApproach::H2Gas, H2GasPvt<Scalar> >::type& getRealPvt()
|
||||
{
|
||||
assert(gasPvtApproach() == approachV);
|
||||
return *static_cast<H2GasPvt<Scalar>* >(realGasPvt_);
|
||||
}
|
||||
|
||||
template <GasPvtApproach approachV>
|
||||
typename std::enable_if<approachV == GasPvtApproach::H2Gas, const H2GasPvt<Scalar> >::type& getRealPvt() const
|
||||
{
|
||||
assert(gasPvtApproach() == approachV);
|
||||
return *static_cast<const H2GasPvt<Scalar>* >(realGasPvt_);
|
||||
}
|
||||
|
||||
const void* realGasPvt() const { return realGasPvt_; }
|
||||
|
||||
GasPvtMultiplexer<Scalar,enableThermal>& operator=(const GasPvtMultiplexer<Scalar,enableThermal>& data)
|
||||
@@ -440,6 +469,9 @@ public:
|
||||
case GasPvtApproach::Co2Gas:
|
||||
realGasPvt_ = new Co2GasPvt<Scalar>(*static_cast<const Co2GasPvt<Scalar>*>(data.realGasPvt_));
|
||||
break;
|
||||
case GasPvtApproach::H2Gas:
|
||||
realGasPvt_ = new H2GasPvt<Scalar>(*static_cast<const H2GasPvt<Scalar>*>(data.realGasPvt_));
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
|
||||
231
opm/material/fluidsystems/blackoilpvt/H2GasPvt.hpp
Normal file
231
opm/material/fluidsystems/blackoilpvt/H2GasPvt.hpp
Normal file
@@ -0,0 +1,231 @@
|
||||
// -*- 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
|
||||
* \copydoc Opm::H2GasPvt
|
||||
*/
|
||||
#ifndef OPM_H2_GAS_PVT_HPP
|
||||
#define OPM_H2_GAS_PVT_HPP
|
||||
|
||||
#include <opm/material/components/SimpleHuDuanH2O.hpp>
|
||||
#include <opm/material/components/H2.hpp>
|
||||
#include <opm/material/binarycoefficients/Brine_H2.hpp>
|
||||
#include <opm/material/common/UniformTabulated2DFunction.hpp>
|
||||
|
||||
#include <vector>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
#if HAVE_ECL_INPUT
|
||||
class EclipseState;
|
||||
class Schedule;
|
||||
#endif
|
||||
|
||||
/*!
|
||||
* \brief This class represents the Pressure-Volume-Temperature relations of the gas phase for H2
|
||||
*/
|
||||
template <class Scalar>
|
||||
class H2GasPvt
|
||||
{
|
||||
using H2O = SimpleHuDuanH2O<Scalar>;
|
||||
using H2 = ::Opm::H2<Scalar>;
|
||||
static const bool extrapolate = true;
|
||||
|
||||
public:
|
||||
// The binary coefficients for brine and H2 used by this fluid system
|
||||
using BinaryCoeffBrineH2 = BinaryCoeff::Brine_H2<Scalar, H2O, H2>;
|
||||
|
||||
explicit H2GasPvt() = default;
|
||||
|
||||
H2GasPvt(size_t numRegions,
|
||||
Scalar T_ref = 288.71, //(273.15 + 15.56)
|
||||
Scalar P_ref = 101325)
|
||||
{
|
||||
setNumRegions(numRegions);
|
||||
for (size_t i = 0; i < numRegions; ++i) {
|
||||
gasReferenceDensity_[i] = H2::gasDensity(T_ref, P_ref, extrapolate);
|
||||
}
|
||||
}
|
||||
|
||||
#if HAVE_ECL_INPUT
|
||||
/*!
|
||||
* \brief Initialize the parameters for H2 gas using an ECL deck.
|
||||
*/
|
||||
void initFromState(const EclipseState& eclState, const Schedule&);
|
||||
#endif
|
||||
|
||||
void setNumRegions(size_t numRegions)
|
||||
{
|
||||
gasReferenceDensity_.resize(numRegions);
|
||||
}
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Initialize the reference densities of all fluids for a given PVT region
|
||||
*/
|
||||
void setReferenceDensities(unsigned regionIdx,
|
||||
Scalar /*rhoRefOil*/,
|
||||
Scalar rhoRefGas,
|
||||
Scalar /*rhoRefWater*/)
|
||||
{
|
||||
gasReferenceDensity_[regionIdx] = rhoRefGas;
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Finish initializing the oil phase PVT properties.
|
||||
*/
|
||||
void initEnd()
|
||||
{
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Return the number of PVT regions which are considered by this PVT-object.
|
||||
*/
|
||||
unsigned numRegions() const
|
||||
{ return gasReferenceDensity_.size(); }
|
||||
|
||||
/*!
|
||||
* \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation internalEnergy(unsigned,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure,
|
||||
const Evaluation& /*rv*/,
|
||||
const Evaluation& /*rvw*/) const
|
||||
{
|
||||
return H2::gasInternalEnergy(temperature, pressure, extrapolate);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation viscosity(unsigned regionIdx,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure,
|
||||
const Evaluation& /*Rv*/,
|
||||
const Evaluation& /*Rvw*/) const
|
||||
{ return saturatedViscosity(regionIdx, temperature, pressure); }
|
||||
|
||||
/*!
|
||||
* \brief Returns the dynamic viscosity [Pa s] of oil saturated gas at given pressure.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturatedViscosity(unsigned /*regionIdx*/,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure) const
|
||||
{
|
||||
return H2::gasViscosity(temperature, pressure);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the formation volume factor [-] of the fluid phase.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure,
|
||||
const Evaluation& /*Rv*/,
|
||||
const Evaluation& /*Rvw*/) const
|
||||
{ return saturatedInverseFormationVolumeFactor(regionIdx, temperature, pressure); }
|
||||
|
||||
/*!
|
||||
* \brief Returns the formation volume factor [-] of oil saturated gas at given pressure.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
|
||||
const Evaluation& temperature,
|
||||
const Evaluation& pressure) const
|
||||
{
|
||||
return H2::gasDensity(temperature, pressure, extrapolate)/gasReferenceDensity_[regionIdx];
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the saturation pressure of the gas phase [Pa] depending on its mass fraction of the oil component
|
||||
*
|
||||
* \param Rv The surface volume of oil component dissolved in what will yield one cubic meter of gas at the surface [-]
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturationPressure(unsigned /*regionIdx*/,
|
||||
const Evaluation& /*temperature*/,
|
||||
const Evaluation& /*Rv*/) const
|
||||
{ return 0.0; /* this is dry gas! */ }
|
||||
|
||||
/*!
|
||||
* \brief Returns the water vaporization factor \f$R_vw\f$ [m^3/m^3] of the water phase.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturatedWaterVaporizationFactor(unsigned /*regionIdx*/,
|
||||
const Evaluation& /*temperature*/,
|
||||
const Evaluation& /*pressure*/) const
|
||||
{ return 0.0; /* this is non-humid gas! */ }
|
||||
|
||||
/*!
|
||||
* \brief Returns the water vaporization factor \f$R_vw\f$ [m^3/m^3] of water saturated gas.
|
||||
*/
|
||||
template <class Evaluation = Scalar>
|
||||
Evaluation saturatedWaterVaporizationFactor(unsigned /*regionIdx*/,
|
||||
const Evaluation& /*temperature*/,
|
||||
const Evaluation& /*pressure*/,
|
||||
const Evaluation& /*saltConcentration*/) const
|
||||
{ return 0.0; }
|
||||
|
||||
/*!
|
||||
* \brief Returns the oil vaporization factor \f$R_v\f$ [m^3/m^3] of the oil phase.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturatedOilVaporizationFactor(unsigned /*regionIdx*/,
|
||||
const Evaluation& /*temperature*/,
|
||||
const Evaluation& /*pressure*/,
|
||||
const Evaluation& /*oilSaturation*/,
|
||||
const Evaluation& /*maxOilSaturation*/) const
|
||||
{ return 0.0; /* this is dry gas! */ }
|
||||
|
||||
/*!
|
||||
* \brief Returns the oil vaporization factor \f$R_v\f$ [m^3/m^3] of the oil phase.
|
||||
*/
|
||||
template <class Evaluation>
|
||||
Evaluation saturatedOilVaporizationFactor(unsigned /*regionIdx*/,
|
||||
const Evaluation& /*temperature*/,
|
||||
const Evaluation& /*pressure*/) const
|
||||
{ return 0.0; /* this is dry gas! */ }
|
||||
|
||||
template <class Evaluation>
|
||||
Evaluation diffusionCoefficient(const Evaluation& temperature,
|
||||
const Evaluation& pressure,
|
||||
unsigned /*compIdx*/) const
|
||||
{
|
||||
return BinaryCoeffBrineH2::gasDiffCoeff(temperature, pressure);
|
||||
}
|
||||
|
||||
const Scalar gasReferenceDensity(unsigned regionIdx) const
|
||||
{ return gasReferenceDensity_[regionIdx]; }
|
||||
|
||||
private:
|
||||
std::vector<Scalar> gasReferenceDensity_;
|
||||
}; // end class H2GasPvt
|
||||
|
||||
} // end namspace Opm
|
||||
|
||||
#endif
|
||||
@@ -32,6 +32,7 @@
|
||||
#include "LiveOilPvt.hpp"
|
||||
#include "OilPvtThermal.hpp"
|
||||
#include "BrineCo2Pvt.hpp"
|
||||
#include "BrineH2Pvt.hpp"
|
||||
|
||||
namespace Opm {
|
||||
|
||||
@@ -67,6 +68,11 @@ class Schedule;
|
||||
codeToCall; \
|
||||
break; \
|
||||
} \
|
||||
case OilPvtApproach::BrineH2: { \
|
||||
auto& pvtImpl = getRealPvt<OilPvtApproach::BrineH2>(); \
|
||||
codeToCall; \
|
||||
break; \
|
||||
} \
|
||||
case OilPvtApproach::NoOil: \
|
||||
throw std::logic_error("Not implemented: Oil PVT of this deck!"); \
|
||||
}
|
||||
@@ -77,7 +83,8 @@ enum class OilPvtApproach {
|
||||
DeadOil,
|
||||
ConstantCompressibilityOil,
|
||||
ThermalOil,
|
||||
BrineCo2
|
||||
BrineCo2,
|
||||
BrineH2
|
||||
};
|
||||
|
||||
/*!
|
||||
@@ -135,7 +142,10 @@ public:
|
||||
delete &getRealPvt<OilPvtApproach::BrineCo2>();
|
||||
break;
|
||||
}
|
||||
|
||||
case OilPvtApproach::BrineH2: {
|
||||
delete &getRealPvt<OilPvtApproach::BrineH2>();
|
||||
break;
|
||||
}
|
||||
case OilPvtApproach::NoOil:
|
||||
break;
|
||||
}
|
||||
@@ -281,6 +291,10 @@ public:
|
||||
realOilPvt_ = new BrineCo2Pvt<Scalar>;
|
||||
break;
|
||||
|
||||
case OilPvtApproach::BrineH2:
|
||||
realOilPvt_ = new BrineH2Pvt<Scalar>;
|
||||
break;
|
||||
|
||||
case OilPvtApproach::NoOil:
|
||||
throw std::logic_error("Not implemented: Oil PVT of this deck!");
|
||||
}
|
||||
@@ -369,6 +383,20 @@ public:
|
||||
|
||||
const void* realOilPvt() const { return realOilPvt_; }
|
||||
|
||||
template <OilPvtApproach approachV>
|
||||
typename std::enable_if<approachV == OilPvtApproach::BrineH2, BrineH2Pvt<Scalar> >::type& getRealPvt()
|
||||
{
|
||||
assert(approach() == approachV);
|
||||
return *static_cast<BrineH2Pvt<Scalar>* >(realOilPvt_);
|
||||
}
|
||||
|
||||
template <OilPvtApproach approachV>
|
||||
typename std::enable_if<approachV == OilPvtApproach::BrineH2, const BrineH2Pvt<Scalar> >::type& getRealPvt() const
|
||||
{
|
||||
assert(approach() == approachV);
|
||||
return *static_cast<const BrineH2Pvt<Scalar>* >(realOilPvt_);
|
||||
}
|
||||
|
||||
OilPvtMultiplexer<Scalar,enableThermal>& operator=(const OilPvtMultiplexer<Scalar,enableThermal>& data)
|
||||
{
|
||||
approach_ = data.approach_;
|
||||
@@ -388,6 +416,9 @@ public:
|
||||
case OilPvtApproach::BrineCo2:
|
||||
realOilPvt_ = new BrineCo2Pvt<Scalar>(*static_cast<const BrineCo2Pvt<Scalar>*>(data.realOilPvt_));
|
||||
break;
|
||||
case OilPvtApproach::BrineH2:
|
||||
realOilPvt_ = new BrineH2Pvt<Scalar>(*static_cast<const BrineH2Pvt<Scalar>*>(data.realOilPvt_));
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
|
||||
@@ -31,6 +31,7 @@
|
||||
#include "ConstantCompressibilityBrinePvt.hpp"
|
||||
#include "WaterPvtThermal.hpp"
|
||||
#include "BrineCo2Pvt.hpp"
|
||||
#include "BrineH2Pvt.hpp"
|
||||
|
||||
#define OPM_WATER_PVT_MULTIPLEXER_CALL(codeToCall) \
|
||||
switch (approach_) { \
|
||||
@@ -54,6 +55,11 @@
|
||||
codeToCall; \
|
||||
break; \
|
||||
} \
|
||||
case WaterPvtApproach::BrineH2: { \
|
||||
auto& pvtImpl = getRealPvt<WaterPvtApproach::BrineH2>(); \
|
||||
codeToCall; \
|
||||
break; \
|
||||
} \
|
||||
case WaterPvtApproach::NoWater: \
|
||||
throw std::logic_error("Not implemented: Water PVT of this deck!"); \
|
||||
}
|
||||
@@ -65,7 +71,8 @@ enum class WaterPvtApproach {
|
||||
ConstantCompressibilityBrine,
|
||||
ConstantCompressibilityWater,
|
||||
ThermalWater,
|
||||
BrineCo2
|
||||
BrineCo2,
|
||||
BrineH2
|
||||
};
|
||||
|
||||
#if HAVE_ECL_INPUT
|
||||
@@ -116,6 +123,10 @@ public:
|
||||
delete &getRealPvt<WaterPvtApproach::BrineCo2>();
|
||||
break;
|
||||
}
|
||||
case WaterPvtApproach::BrineH2: {
|
||||
delete &getRealPvt<WaterPvtApproach::BrineH2>();
|
||||
break;
|
||||
}
|
||||
case WaterPvtApproach::NoWater:
|
||||
break;
|
||||
}
|
||||
@@ -269,6 +280,10 @@ public:
|
||||
realWaterPvt_ = new BrineCo2Pvt<Scalar>;
|
||||
break;
|
||||
|
||||
case WaterPvtApproach::BrineH2:
|
||||
realWaterPvt_ = new BrineH2Pvt<Scalar>;
|
||||
break;
|
||||
|
||||
case WaterPvtApproach::NoWater:
|
||||
throw std::logic_error("Not implemented: Water PVT of this deck!");
|
||||
}
|
||||
@@ -341,6 +356,20 @@ public:
|
||||
return *static_cast<const BrineCo2Pvt<Scalar>* >(realWaterPvt_);
|
||||
}
|
||||
|
||||
template <WaterPvtApproach approachV>
|
||||
typename std::enable_if<approachV == WaterPvtApproach::BrineH2, BrineH2Pvt<Scalar> >::type& getRealPvt()
|
||||
{
|
||||
assert(approach() == approachV);
|
||||
return *static_cast<BrineH2Pvt<Scalar>* >(realWaterPvt_);
|
||||
}
|
||||
|
||||
template <WaterPvtApproach approachV>
|
||||
typename std::enable_if<approachV == WaterPvtApproach::BrineH2, const BrineH2Pvt<Scalar> >::type& getRealPvt() const
|
||||
{
|
||||
assert(approach() == approachV);
|
||||
return *static_cast<const BrineH2Pvt<Scalar>* >(realWaterPvt_);
|
||||
}
|
||||
|
||||
const void* realWaterPvt() const { return realWaterPvt_; }
|
||||
|
||||
WaterPvtMultiplexer<Scalar,enableThermal,enableBrine>& operator=(const WaterPvtMultiplexer<Scalar,enableThermal,enableBrine>& data)
|
||||
@@ -359,6 +388,9 @@ public:
|
||||
case WaterPvtApproach::BrineCo2:
|
||||
realWaterPvt_ = new BrineCo2Pvt<Scalar>(*static_cast<const BrineCo2Pvt<Scalar>*>(data.realWaterPvt_));
|
||||
break;
|
||||
case WaterPvtApproach::BrineH2:
|
||||
realWaterPvt_ = new BrineH2Pvt<Scalar>(*static_cast<const BrineH2Pvt<Scalar>*>(data.realWaterPvt_));
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
|
||||
@@ -26,6 +26,7 @@
|
||||
#include <opm/input/eclipse/Parser/ParserKeywords/C.hpp>
|
||||
#include <opm/input/eclipse/Parser/ParserKeywords/F.hpp>
|
||||
#include <opm/input/eclipse/Parser/ParserKeywords/G.hpp>
|
||||
#include <opm/input/eclipse/Parser/ParserKeywords/H.hpp>
|
||||
#include <opm/input/eclipse/Parser/ParserKeywords/M.hpp>
|
||||
#include <opm/input/eclipse/Parser/ParserKeywords/N.hpp>
|
||||
#include <opm/input/eclipse/Parser/ParserKeywords/O.hpp>
|
||||
@@ -585,6 +586,7 @@ Runspec::Runspec( const Deck& deck )
|
||||
, m_nupcol( )
|
||||
, m_tracers( deck )
|
||||
, m_co2storage (false)
|
||||
, m_h2storage (false)
|
||||
, m_micp (false)
|
||||
{
|
||||
if (DeckSection::hasRUNSPEC(deck)) {
|
||||
@@ -619,6 +621,13 @@ Runspec::Runspec( const Deck& deck )
|
||||
|
||||
}
|
||||
|
||||
if (runspecSection.hasKeyword<ParserKeywords::H2STORE>()) {
|
||||
m_h2storage = true;
|
||||
std::string msg = "The H2 storage option is given. PVT properties from the Brine-H2 system is used \n"
|
||||
"See the OPM manual for details on the used models.";
|
||||
OpmLog::note(msg);
|
||||
}
|
||||
|
||||
if (runspecSection.hasKeyword<ParserKeywords::MICP>()) {
|
||||
m_micp = true;
|
||||
std::string msg = "The MICP option is given. Single phase (WATER) + 3 transported components + \n"
|
||||
@@ -645,6 +654,7 @@ Runspec Runspec::serializationTestObject()
|
||||
result.m_sfuncctrl = SatFuncControls::serializationTestObject();
|
||||
result.m_nupcol = Nupcol::serializationTestObject();
|
||||
result.m_co2storage = true;
|
||||
result.m_h2storage = true;
|
||||
result.m_micp = true;
|
||||
|
||||
return result;
|
||||
@@ -710,6 +720,11 @@ bool Runspec::co2Storage() const noexcept
|
||||
return this->m_co2storage;
|
||||
}
|
||||
|
||||
bool Runspec::h2Storage() const noexcept
|
||||
{
|
||||
return this->m_h2storage;
|
||||
}
|
||||
|
||||
bool Runspec::micp() const noexcept
|
||||
{
|
||||
return this->m_micp;
|
||||
@@ -752,6 +767,7 @@ bool Runspec::rst_cmp(const Runspec& full_spec, const Runspec& rst_spec) {
|
||||
full_spec.saturationFunctionControls() == rst_spec.saturationFunctionControls() &&
|
||||
full_spec.m_nupcol == rst_spec.m_nupcol &&
|
||||
full_spec.m_co2storage == rst_spec.m_co2storage &&
|
||||
full_spec.m_h2storage == rst_spec.m_h2storage &&
|
||||
full_spec.m_micp == rst_spec.m_micp &&
|
||||
Welldims::rst_cmp(full_spec.wellDimensions(), rst_spec.wellDimensions());
|
||||
}
|
||||
@@ -769,6 +785,7 @@ bool Runspec::operator==(const Runspec& data) const {
|
||||
this->saturationFunctionControls() == data.saturationFunctionControls() &&
|
||||
this->m_nupcol == data.m_nupcol &&
|
||||
this->m_co2storage == data.m_co2storage &&
|
||||
this->m_h2storage == data.m_h2storage &&
|
||||
this->m_micp == data.m_micp;
|
||||
}
|
||||
|
||||
|
||||
@@ -0,0 +1,6 @@
|
||||
{
|
||||
"name": "H2STORE",
|
||||
"sections": [
|
||||
"RUNSPEC"
|
||||
]
|
||||
}
|
||||
@@ -1062,6 +1062,7 @@ set( keywords
|
||||
001_Eclipse300/G/GASWAT
|
||||
001_Eclipse300/G/GCONPROD
|
||||
001_Eclipse300/G/GSF
|
||||
001_Eclipse300/H/H2STORE
|
||||
001_Eclipse300/H/HEATCR
|
||||
001_Eclipse300/H/HEATCRT
|
||||
001_Eclipse300/H/HWELLS
|
||||
|
||||
39
src/opm/material/components/H2.cpp
Normal file
39
src/opm/material/components/H2.cpp
Normal file
@@ -0,0 +1,39 @@
|
||||
// -*- 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.
|
||||
*/
|
||||
|
||||
#include <config.h>
|
||||
#include <opm/material/components/H2.hpp>
|
||||
|
||||
#include "h2tables.inc"
|
||||
|
||||
namespace Opm {
|
||||
|
||||
template<>
|
||||
const UniformTabulated2DFunction<double>&
|
||||
H2<double>::tabulatedDensity = H2Tables::tabulatedDensity;
|
||||
|
||||
template<>
|
||||
const UniformTabulated2DFunction<double>&
|
||||
H2<float>::tabulatedDensity = H2Tables::tabulatedDensity;
|
||||
|
||||
} // namespace Opm
|
||||
20452
src/opm/material/components/h2tables.inc
Normal file
20452
src/opm/material/components/h2tables.inc
Normal file
File diff suppressed because it is too large
Load Diff
@@ -82,11 +82,11 @@ initFromState(const EclipseState& eclState, const Schedule& schedule)
|
||||
setEnableVaporizedWater(eclState.getSimulationConfig().hasVAPWAT());
|
||||
|
||||
if (eclState.getSimulationConfig().hasDISGASW()) {
|
||||
if (eclState.runspec().co2Storage())
|
||||
if (eclState.runspec().co2Storage() || eclState.runspec().h2Storage())
|
||||
setEnableDissolvedGasInWater(eclState.getSimulationConfig().hasDISGASW());
|
||||
else
|
||||
OPM_THROW(std::runtime_error,
|
||||
"DISGASW only supported in combination with CO2STORE");
|
||||
"DISGASW only supported in combination with CO2STORE or H2STORE");
|
||||
}
|
||||
|
||||
if (phaseIsActive(gasPhaseIdx)) {
|
||||
@@ -135,6 +135,19 @@ initFromState(const EclipseState& eclState, const Schedule& schedule)
|
||||
}
|
||||
}
|
||||
|
||||
// Use molar mass of H2 and Brine as default in H2STORE keyword
|
||||
if (eclState.runspec().h2Storage()) {
|
||||
for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
|
||||
if (phaseIsActive(oilPhaseIdx)) // The oil component is used for the brine if OIL is active
|
||||
molarMass_[regionIdx][oilCompIdx] = BrineH2Pvt<Scalar>::Brine::molarMass();
|
||||
if (!phaseIsActive(gasPhaseIdx)) {
|
||||
OPM_THROW(std::runtime_error,
|
||||
"H2STORE requires gas phase\n");
|
||||
}
|
||||
molarMass_[regionIdx][gasCompIdx] = BrineH2Pvt<Scalar>::H2::molarMass();
|
||||
}
|
||||
}
|
||||
|
||||
setEnableDiffusion(eclState.getSimulationConfig().isDiffusive());
|
||||
if (enableDiffusion()) {
|
||||
const auto& diffCoeffTables = eclState.getTableManager().getDiffusionCoefficientTable();
|
||||
|
||||
75
src/opm/material/fluidsystems/blackoilpvt/BrineH2Pvt.cpp
Normal file
75
src/opm/material/fluidsystems/blackoilpvt/BrineH2Pvt.cpp
Normal file
@@ -0,0 +1,75 @@
|
||||
// -*- 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.
|
||||
*/
|
||||
|
||||
#include <config.h>
|
||||
#include <opm/material/fluidsystems/blackoilpvt/BrineH2Pvt.hpp>
|
||||
|
||||
#include <opm/common/OpmLog/OpmLog.hpp>
|
||||
|
||||
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
|
||||
#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
template<class Scalar>
|
||||
void BrineH2Pvt<Scalar>::
|
||||
initFromState(const EclipseState& eclState, const Schedule&)
|
||||
{
|
||||
if( !eclState.getTableManager().getDensityTable().empty()) {
|
||||
OpmLog::warning("H2STORE is enabled but DENSITY is in the deck. \n"
|
||||
"The surface density is computed based on H2-BRINE PVT "
|
||||
"at standard conditions (STCOND) and DENSITY is ignored ");
|
||||
}
|
||||
|
||||
if(eclState.getTableManager().hasTables("PVDO") ||
|
||||
!eclState.getTableManager().getPvtgTables().empty()) {
|
||||
OpmLog::warning("H2STORE is enabled but PVDO or PVTO is in the deck. \n"
|
||||
"H2 PVT properties are calculated internally, "
|
||||
"and PVDO/PVTO input is ignored.");
|
||||
}
|
||||
// Check if DISGAS has been activated (enables H2 dissolved in brine)
|
||||
setEnableDissolvedGas(eclState.getSimulationConfig().hasDISGASW() || eclState.getSimulationConfig().hasDISGAS());
|
||||
|
||||
// We only supported single pvt region for the H2-brine module
|
||||
size_t numRegions = 1;
|
||||
setNumRegions(numRegions);
|
||||
size_t regionIdx = 0;
|
||||
|
||||
// Currently we only support constant salinity
|
||||
const Scalar molality = eclState.getTableManager().salinity(); // mol/kg
|
||||
const Scalar MmNaCl = 58e-3; // molar mass of NaCl [kg/mol]
|
||||
Brine::salinity = 1 / ( 1 + 1 / (molality*MmNaCl)); // convert to mass fraction
|
||||
salinity_[regionIdx] = molality; // molality used in BrineH2Pvt functions
|
||||
|
||||
// set the surface conditions using the STCOND keyword
|
||||
Scalar T_ref = eclState.getTableManager().stCond().temperature;
|
||||
Scalar P_ref = eclState.getTableManager().stCond().pressure;
|
||||
|
||||
brineReferenceDensity_[regionIdx] = Brine::liquidDensity(T_ref, P_ref, extrapolate);
|
||||
h2ReferenceDensity_[regionIdx] = H2::gasDensity(T_ref, P_ref, extrapolate);
|
||||
}
|
||||
|
||||
template class BrineH2Pvt<double>;
|
||||
template class BrineH2Pvt<float>;
|
||||
|
||||
} // namespace Opm
|
||||
@@ -36,6 +36,8 @@ initFromState(const EclipseState& eclState, const Schedule& schedule)
|
||||
|
||||
if (eclState.runspec().co2Storage())
|
||||
setApproach(GasPvtApproach::Co2Gas);
|
||||
else if (eclState.runspec().h2Storage())
|
||||
setApproach(GasPvtApproach::H2Gas);
|
||||
else if (enableThermal && eclState.getSimulationConfig().isThermal())
|
||||
setApproach(GasPvtApproach::ThermalGas);
|
||||
else if (!eclState.getTableManager().getPvtgwTables().empty() &&
|
||||
|
||||
63
src/opm/material/fluidsystems/blackoilpvt/H2GasPvt.cpp
Normal file
63
src/opm/material/fluidsystems/blackoilpvt/H2GasPvt.cpp
Normal file
@@ -0,0 +1,63 @@
|
||||
// -*- 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.
|
||||
*/
|
||||
|
||||
#include <config.h>
|
||||
#include <opm/material/fluidsystems/blackoilpvt/H2GasPvt.hpp>
|
||||
|
||||
#include <opm/common/OpmLog/OpmLog.hpp>
|
||||
|
||||
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
|
||||
#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
template<class Scalar>
|
||||
void H2GasPvt<Scalar>::
|
||||
initFromState(const EclipseState& eclState, const Schedule&)
|
||||
{
|
||||
if( !eclState.getTableManager().getDensityTable().empty()) {
|
||||
OpmLog::warning("H2STORE is enabled but DENSITY is in the deck. \n"
|
||||
"The surface density is computed based on H2-BRINE PVT "
|
||||
"at standard conditions (STCOND) and DENSITY is ignored ");
|
||||
}
|
||||
|
||||
if( eclState.getTableManager().hasTables("PVDG") || !eclState.getTableManager().getPvtgTables().empty()) {
|
||||
OpmLog::warning("H2STORE is enabled but PVDG or PVTG is in the deck. \n"
|
||||
"H2 pvt properties are calculated based on ideal gas relations, "
|
||||
"and PVDG/PVTG input is ignored.");
|
||||
}
|
||||
|
||||
// We only supported single pvt region for the H2-brine module
|
||||
size_t numRegions = 1;
|
||||
setNumRegions(numRegions);
|
||||
size_t regionIdx = 0;
|
||||
Scalar T_ref = eclState.getTableManager().stCond().temperature;
|
||||
Scalar P_ref = eclState.getTableManager().stCond().pressure;
|
||||
gasReferenceDensity_[regionIdx] = H2::gasDensity(T_ref, P_ref, extrapolate);
|
||||
initEnd();
|
||||
}
|
||||
|
||||
template class H2GasPvt<double>;
|
||||
template class H2GasPvt<float>;
|
||||
|
||||
}
|
||||
@@ -40,6 +40,8 @@ initFromState(const EclipseState& eclState, const Schedule& schedule)
|
||||
// and water/brine + gas
|
||||
if (eclState.runspec().co2Storage())
|
||||
setApproach(OilPvtApproach::BrineCo2);
|
||||
else if (eclState.runspec().h2Storage())
|
||||
setApproach(OilPvtApproach::BrineH2);
|
||||
else if (enableThermal && eclState.getSimulationConfig().isThermal())
|
||||
setApproach(OilPvtApproach::ThermalOil);
|
||||
else if (!eclState.getTableManager().getPvcdoTable().empty())
|
||||
|
||||
@@ -40,6 +40,8 @@ initFromState(const EclipseState& eclState, const Schedule& schedule)
|
||||
// and water/brine + gas
|
||||
if (eclState.runspec().co2Storage())
|
||||
setApproach(WaterPvtApproach::BrineCo2);
|
||||
else if (eclState.runspec().h2Storage())
|
||||
setApproach(WaterPvtApproach::BrineH2);
|
||||
else if (enableThermal && eclState.getSimulationConfig().isThermal())
|
||||
setApproach(WaterPvtApproach::ThermalWater);
|
||||
else if (!eclState.getTableManager().getPvtwTable().empty())
|
||||
|
||||
@@ -1065,6 +1065,28 @@ BOOST_AUTO_TEST_CASE(Co2Storage_oilwater) {
|
||||
BOOST_CHECK_THROW( Runspec{deck}, std::runtime_error );
|
||||
}
|
||||
|
||||
BOOST_AUTO_TEST_CASE(H2Storage) {
|
||||
const std::string input = R"(
|
||||
RUNSPEC
|
||||
OIL
|
||||
GAS
|
||||
H2STORE
|
||||
)";
|
||||
|
||||
Parser parser;
|
||||
|
||||
auto deck = parser.parseString(input);
|
||||
|
||||
Runspec runspec( deck );
|
||||
const auto& phases = runspec.phases();
|
||||
BOOST_CHECK_EQUAL( 2U, phases.size() );
|
||||
BOOST_CHECK( phases.active( Phase::OIL ) );
|
||||
BOOST_CHECK( phases.active( Phase::GAS ) );
|
||||
BOOST_CHECK( runspec.h2Storage() );
|
||||
|
||||
|
||||
}
|
||||
|
||||
BOOST_AUTO_TEST_CASE(NUPCOL_DEFAULT) {
|
||||
Nupcol np;
|
||||
auto default_value = ParserKeywords::NUPCOL::NUM_ITER::defaultValue;
|
||||
|
||||
Reference in New Issue
Block a user