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

Merge pull request #3683 from svenn-t/h2store_updates

Browse Source 
H2STORE updates
This commit is contained in:
Tor Harald Sandve
2023-09-21 09:19:57 +02:00
committed by GitHub
parent e5bd4129a1 cf4a905660
commit 8ea35f2d0d
15 changed files with 47723 additions and 26700 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@@ -477,6 +477,7 @@ if(ENABLE_ECL_INPUT)
tests/test_PAvgDynamicSourceData.cpp
tests/test_Serialization.cpp
tests/material/test_co2brinepvt.cpp
tests/material/test_h2brinepvt.cpp
tests/material/test_eclblackoilfluidsystem.cpp
tests/material/test_eclblackoilpvt.cpp
tests/material/test_eclmateriallawmanager.cpp

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

@@ -49,27 +49,30 @@ public:
*
* \param temperature temperature [K]
* \param pg gas phase pressure [Pa]
* \param salinity salinity [mol NaCl / kg solution]
* \param salinity salinity [kg 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)
static Evaluation calculateMoleFractions(const Evaluation& temperature,
const Evaluation& pg,
const Evaluation& salinity,
bool extrapolate = false)
{
// Convert salinity to molality from mass fraction
Evaluation X_NaCl = salinityToMolality_(salinity);
// 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);
Evaluation lnGammaH2 = activityCoefficient_(temperature, X_NaCl);
// Eq. (4) to get mole fraction of H2 in brine
xH2 = exp(lnYH2 + lnPg + lnPhiH2 - lnKh - PF - lnGammaH2);
Evaluation xH2 = exp(lnYH2 + lnPg + lnPhiH2 - lnKh - PF - lnGammaH2);
return xH2;
}
/*!
@@ -100,7 +103,7 @@ public:
* \param salinity salinity [mol NaCl / kg solution]
*/
template <class Evaluation>
static Evaluation activityCoefficient_(const Evaluation& temperature, Scalar salinity)
static Evaluation activityCoefficient_(const Evaluation& temperature, const Evaluation& salinity)
{
// Linear approximation in temperature with following parameters (Table 5)
static const Scalar a[2] = {0.64485, 0.00142};
@@ -189,6 +192,22 @@ public:
return fullerMethod(M, SigmaNu, temperature, pressure);
}
private:
/*!
* \brief Convert mass fraction to molality NaCl [mol NaCl / kg water]
*
*/
template <class Evaluation>
static Evaluation salinityToMolality_(const Evaluation& salinity) {
// Molar mass NaCl
const Scalar M_NaCl = 58.44e-3;
// Convert
const Evaluation X_NaCl = salinity / ((1 - salinity) * M_NaCl);
return X_NaCl;
}
}; // end class Brine_H2
} // end namespace BinaryCoeff

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

@@ -58,9 +58,13 @@ class H2 : public Component<Scalar, H2<Scalar> >
{
using IdealGas = Opm::IdealGas<Scalar>;
static const UniformTabulated2DFunction<double>& tabulatedEnthalpy;
static const UniformTabulated2DFunction<double>& tabulatedDensity;
public:
// For H2Tables class
static const Scalar brineSalinity;
/*!
* \brief A human readable name for the \f$H_2\f$.
*/
@@ -212,13 +216,10 @@ public:
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;
const Evaluation h = gasEnthalpy(temperature, pressure, extrapolate);
const Evaluation rho = gasDensity(temperature, pressure, extrapolate);
return R * criticalTemperature() * (derivIdealHelmholtzWrtRecipRedTemp(T_red)
+ derivResHelmholtzWrtRecipRedTemp(T_red, rho_red)) / molarMass();
return h - (pressure / rho);
}
/*!
@@ -232,11 +233,7 @@ public:
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);
return tabulatedEnthalpy.eval(temperature, pressure, extrapolate);
}
/*!

225
opm/material/fluidsystems/blackoilpvt/BrineH2Pvt.hpp
View File

@@ -31,7 +31,7 @@ copyright holders.
#include <opm/material/binarycoefficients/Brine_H2.hpp>
#include <opm/material/components/SimpleHuDuanH2O.hpp>
#include <opm/material/components/Brine.hpp>
#include <opm/material/components/BrineDynamic.hpp>
#include <opm/material/components/H2.hpp>
#include <opm/material/common/UniformTabulated2DFunction.hpp>
#include <opm/material/common/Valgrind.hpp>
@@ -54,7 +54,7 @@ class BrineH2Pvt
static const bool extrapolate = true;
public:
using H2O = SimpleHuDuanH2O<Scalar>;
using Brine = ::Opm::Brine<Scalar, H2O>;
using Brine = ::Opm::BrineDynamic<Scalar, H2O>;
using H2 = ::Opm::H2<Scalar>;
// The binary coefficients for brine and H2 used by this fluid system
@@ -70,10 +70,9 @@ public:
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);
brineReferenceDensity_[i] = Brine::liquidDensity(T_ref, P_ref, salinity_[i], true);
}
}
@@ -116,14 +115,23 @@ public:
}
/*!
* \brief Specify whether the PVT model should consider that the H2 component can
* dissolve in the brine phase
*
* By default, dissolved H2 is considered.
*/
* \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 Specify whether the PVT model should consider salt concentration from
* the fluidstate or a fixed salinty
*
* By default, fixed salinity is considered
*/
void setEnableSaltConcentration(bool yesno)
{ enableSaltConcentration_ = yesno; }
/*!
* \brief Return the number of PVT regions which are considered by this PVT-object.
*/
@@ -134,25 +142,36 @@ public:
* \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
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.");
const Evaluation salinity = salinityFromConcentration(regionIdx, temperature, pressure, saltConcentration);
const Evaluation xlH2 = convertRsToXoG_(Rs,regionIdx);
return (liquidEnthalpyBrineH2_(temperature,
pressure,
salinity,
xlH2)
- pressure / density_(regionIdx, temperature, pressure, Rs, salinity ));
}
/*!
* \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 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.");
const Evaluation xlH2 = convertRsToXoG_(Rs,regionIdx);
return (liquidEnthalpyBrineH2_(temperature,
pressure,
Evaluation(salinity_[regionIdx]),
xlH2)
- pressure / density_(regionIdx, temperature, pressure, Rs, Evaluation(salinity_[regionIdx])));
}
/*!
@@ -175,9 +194,10 @@ public:
Evaluation saturatedViscosity(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& /*saltConcentration*/) const
const Evaluation& saltConcentration) const
{
return saturatedViscosity(regionIdx, temperature, pressure);
const Evaluation salinity = salinityFromConcentration(regionIdx, temperature, pressure, saltConcentration);
return Brine::liquidViscosity(temperature, pressure, salinity);
}
/*!
@@ -188,21 +208,21 @@ public:
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& /*Rsw*/,
const Evaluation& /*saltConcentration*/) const
const Evaluation& saltConcentration) const
{
//TODO: The viscosity does not yet depend on the composition
return saturatedViscosity(regionIdx, temperature, pressure);
return saturatedViscosity(regionIdx, temperature, pressure, saltConcentration);
}
/*!
* \brief Returns the dynamic viscosity [Pa s] of oil saturated gas at given pressure.
*/
template <class Evaluation>
Evaluation saturatedViscosity(unsigned /*regionIdx*/,
Evaluation saturatedViscosity(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure) const
{
return Brine::liquidViscosity(temperature, pressure);
return Brine::liquidViscosity(temperature, pressure, Evaluation(salinity_[regionIdx]));
}
/*!
@@ -212,9 +232,12 @@ public:
Evaluation saturatedInverseFormationVolumeFactor(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& /*saltconcentration*/) const
const Evaluation& saltconcentration) const
{
return saturatedInverseFormationVolumeFactor(regionIdx, temperature, pressure);
const Evaluation salinity = salinityFromConcentration(regionIdx, temperature, pressure, saltconcentration);
Evaluation rsSat = rsSat_(regionIdx, temperature, pressure, salinity);
return (1.0 - convertRsToXoG_(rsSat,regionIdx)) *
density_(regionIdx, temperature, pressure, rsSat, salinity) / brineReferenceDensity_[regionIdx];
}
/*!
@@ -225,9 +248,11 @@ public:
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& Rs,
const Evaluation& /*saltConcentration*/) const
const Evaluation& saltConcentration) const
{
return inverseFormationVolumeFactor(regionIdx, temperature, pressure, Rs);
const Evaluation salinity = salinityFromConcentration(regionIdx, temperature, pressure, saltConcentration);
return (1.0 - convertRsToXoG_(Rs,regionIdx)) *
density_(regionIdx, temperature, pressure, Rs, salinity) / brineReferenceDensity_[regionIdx];
}
/*!
@@ -239,8 +264,9 @@ public:
const Evaluation& pressure,
const Evaluation& Rs) const
{
return (1.0 - convertRsToXoG_(Rs, regionIdx)) * density_(regionIdx, temperature, pressure, Rs) /
brineReferenceDensity_[regionIdx];
return (1.0 - convertRsToXoG_(Rs, regionIdx)) *
density_(regionIdx, temperature, pressure, Rs, Evaluation(salinity_[regionIdx])) /
brineReferenceDensity_[regionIdx];
}
/*!
@@ -251,9 +277,10 @@ public:
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];
Evaluation rsSat = rsSat_(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
return (1.0 - convertRsToXoG_(rsSat, regionIdx)) *
density_(regionIdx, temperature, pressure, rsSat, Evaluation(salinity_[regionIdx])) /
brineReferenceDensity_[regionIdx];
}
/*!
@@ -294,7 +321,7 @@ public:
const Evaluation& /*maxOilSaturation*/) const
{
//TODO support VAPPARS
return rsSat_(regionIdx, temperature, pressure);
return rsSat_(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
}
/*!
@@ -304,9 +331,10 @@ public:
Evaluation saturatedGasDissolutionFactor(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& /*saltConcentration*/) const
const Evaluation& saltConcentration) const
{
return rsSat_(regionIdx, temperature, pressure);
const Evaluation salinity = salinityFromConcentration(regionIdx, temperature, pressure, saltConcentration);
return rsSat_(regionIdx, temperature, pressure, salinity);
}
/*!
@@ -317,7 +345,7 @@ public:
const Evaluation& temperature,
const Evaluation& pressure) const
{
return rsSat_(regionIdx, temperature, pressure);
return rsSat_(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
}
const Scalar oilReferenceDensity(unsigned regionIdx) const
@@ -348,14 +376,14 @@ public:
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
const Evaluation& mu_pure = H2O::liquidViscosity(temperature, pressure, extrapolate) * 1e3; // [cP]
const Evaluation& mu_brine = Brine::liquidViscosity(temperature, pressure, Evaluation(salinity_[0])) * 1e3;
// 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
@@ -366,6 +394,7 @@ private:
std::vector<Scalar> h2ReferenceDensity_;
std::vector<Scalar> salinity_;
bool enableDissolution_ = true;
bool enableSaltConcentration_ = false;
/*!
* \brief Calculate density of aqueous solution (H2O-NaCl/brine and H2).
@@ -378,15 +407,17 @@ private:
LhsEval density_(unsigned regionIdx,
const LhsEval& temperature,
const LhsEval& pressure,
const LhsEval& Rs) const
const LhsEval& Rs,
const LhsEval& salinity) const
{
// convert Rs to mole fraction (via mass fraction)
LhsEval xlH2 = convertXoGToxoG_(convertRsToXoG_(Rs,regionIdx));
LhsEval xlH2 = convertXoGToxoG_(convertRsToXoG_(Rs,regionIdx), salinity);
// calculate the density of solution
LhsEval result = liquidDensity_(temperature,
pressure,
xlH2);
xlH2,
salinity);
Valgrind::CheckDefined(result);
return result;
@@ -403,7 +434,8 @@ private:
template <class LhsEval>
LhsEval liquidDensity_(const LhsEval& T,
const LhsEval& pl,
const LhsEval& xlH2) const
const LhsEval& xlH2,
const LhsEval& salinity) const
{
// check input variables
Valgrind::CheckDefined(T);
@@ -427,7 +459,7 @@ private:
}
// calculate individual contribution to density
const LhsEval& rho_brine = Brine::liquidDensity(T, pl, extrapolate);
const LhsEval& rho_brine = Brine::liquidDensity(T, pl, salinity, 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;
@@ -491,10 +523,10 @@ private:
* \param XoG mass fraction [-]
*/
template <class LhsEval>
LhsEval convertXoGToxoG_(const LhsEval& XoG) const
LhsEval convertXoGToxoG_(const LhsEval& XoG, const LhsEval& salinity) const
{
Scalar M_H2 = H2::molarMass();
Scalar M_Brine = Brine::molarMass();
LhsEval M_Brine = Brine::molarMass(salinity);
return XoG*M_Brine / (M_H2*(1 - XoG) + XoG*M_Brine);
}
@@ -504,10 +536,10 @@ private:
* \param xoG mole fraction [-]
*/
template <class LhsEval>
LhsEval convertxoGToXoG(const LhsEval& xoG) const
LhsEval convertxoGToXoG(const LhsEval& xoG, const LhsEval& salinity) const
{
Scalar M_H2 = H2::molarMass();
Scalar M_Brine = Brine::molarMass();
LhsEval M_Brine = Brine::molarMass(salinity);
return xoG*M_H2 / (xoG*(M_H2 - M_Brine) + M_Brine);
}
@@ -538,24 +570,101 @@ private:
template <class LhsEval>
LhsEval rsSat_(unsigned regionIdx,
const LhsEval& temperature,
const LhsEval& pressure) const
const LhsEval& pressure,
const LhsEval& salinity) 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);
LhsEval xlH2 = BinaryCoeffBrineH2::calculateMoleFractions(temperature, pressure, salinity, extrapolate);
// normalize the phase compositions
xlH2 = max(0.0, min(1.0, xlH2));
return convertXoGToRs(convertxoGToXoG(xlH2), regionIdx);
return convertXoGToRs(convertxoGToXoG(xlH2, salinity), regionIdx);
}
template <class LhsEval>
static LhsEval liquidEnthalpyBrineH2_(const LhsEval& T,
const LhsEval& p,
const LhsEval& salinity,
const LhsEval& X_H2_w)
{
/* X_H2_w : mass fraction of H2 in brine */
/*NOTE: The heat of dissolution of H2 in brine is not included at the moment!*/
/*Numerical coefficents from PALLISER*/
static constexpr Scalar f[] = {
2.63500E-1, 7.48368E-6, 1.44611E-6, -3.80860E-10
};
/*Numerical coefficents from MICHAELIDES for the enthalpy of brine*/
static constexpr Scalar a[4][3] = {
{ 9633.6, -4080.0, +286.49 },
{ +166.58, +68.577, -4.6856 },
{ -0.90963, -0.36524, +0.249667E-1 },
{ +0.17965E-2, +0.71924E-3, -0.4900E-4 }
};
LhsEval theta, h_NaCl;
LhsEval h_ls1, d_h;
LhsEval delta_h;
LhsEval hg, hw;
// Temperature in Celsius
theta = T - 273.15;
// Regularization
Scalar scalarTheta = scalarValue(theta);
Scalar S_lSAT = f[0] + scalarTheta*(f[1] + scalarTheta*(f[2] + scalarTheta*f[3]));
LhsEval S = salinity;
if (S > S_lSAT)
S = S_lSAT;
hw = H2O::liquidEnthalpy(T, p) /1E3; /* kJ/kg */
/*DAUBERT and DANNER*/
/*U=*/h_NaCl = (3.6710E4*T + 0.5*(6.2770E1)*T*T - ((6.6670E-2)/3)*T*T*T
+((2.8000E-5)/4)*(T*T*T*T))/(58.44E3)- 2.045698e+02; /* kJ/kg */
LhsEval m = 1E3/58.44 * S/(1-S);
int i = 0;
int j = 0;
d_h = 0;
for (i = 0; i<=3; i++) {
for (j=0; j<=2; j++) {
d_h = d_h + a[i][j] * pow(theta, static_cast<Scalar>(i)) * pow(m, j);
}
}
/* heat of dissolution for halite according to Michaelides 1971 */
delta_h = (4.184/(1E3 + (58.44 * m)))*d_h;
/* Enthalpy of brine without H2 */
h_ls1 =(1-S)*hw + S*h_NaCl + S*delta_h; /* kJ/kg */
/* enthalpy contribution of H2 gas (kJ/kg) */
hg = H2::gasEnthalpy(T, p, extrapolate) / 1E3;
/* Enthalpy of brine with dissolved H2 */
return (h_ls1 - X_H2_w*hw + hg*X_H2_w)*1E3; /*J/kg*/
}
template <class LhsEval>
const LhsEval salinityFromConcentration(unsigned regionIdx,
const LhsEval&T,
const LhsEval& P,
const LhsEval& saltConcentration) const
{
if (enableSaltConcentration_)
return saltConcentration/H2O::liquidDensity(T, P, true);
return salinity(regionIdx);
}
}; // end class BrineH2Pvt
} // end namespace Opm
#endif

196
opm/material/fluidsystems/blackoilpvt/H2GasPvt.hpp
View File

@@ -28,6 +28,7 @@
#define OPM_H2_GAS_PVT_HPP
#include <opm/material/components/SimpleHuDuanH2O.hpp>
#include <opm/material/components/BrineDynamic.hpp>
#include <opm/material/components/H2.hpp>
#include <opm/material/binarycoefficients/Brine_H2.hpp>
#include <opm/material/common/UniformTabulated2DFunction.hpp>
@@ -48,6 +49,7 @@ template <class Scalar>
class H2GasPvt
{
using H2O = SimpleHuDuanH2O<Scalar>;
using Brine = ::Opm::BrineDynamic<Scalar, H2O>;
using H2 = ::Opm::H2<Scalar>;
static const bool extrapolate = true;
@@ -57,13 +59,16 @@ public:
explicit H2GasPvt() = default;
H2GasPvt(size_t numRegions,
H2GasPvt(const std::vector<Scalar>& salinity,
Scalar T_ref = 288.71, //(273.15 + 15.56)
Scalar P_ref = 101325)
: salinity_(salinity)
{
int numRegions = salinity_.size();
setNumRegions(numRegions);
for (size_t i = 0; i < numRegions; ++i) {
for (int i = 0; i < numRegions; ++i) {
gasReferenceDensity_[i] = H2::gasDensity(T_ref, P_ref, extrapolate);
brineReferenceDensity_[i] = Brine::liquidDensity(T_ref, P_ref, salinity_[i], extrapolate);
}
}
@@ -77,6 +82,8 @@ public:
void setNumRegions(size_t numRegions)
{
gasReferenceDensity_.resize(numRegions);
brineReferenceDensity_.resize(numRegions);
salinity_.resize(numRegions);
}
void setVapPars(const Scalar, const Scalar)
@@ -88,13 +95,23 @@ public:
* \brief Initialize the reference densities of all fluids for a given PVT region
*/
void setReferenceDensities(unsigned regionIdx,
Scalar /*rhoRefOil*/,
Scalar rhoRefBrine,
Scalar rhoRefGas,
Scalar /*rhoRefWater*/)
{
gasReferenceDensity_[regionIdx] = rhoRefGas;
brineReferenceDensity_[regionIdx] = rhoRefBrine;
}
/*!
* \brief Specify whether the PVT model should consider that the water component can
* vaporize in the gas phase
*
* By default, vaporized water is considered.
*/
void setEnableVaporizationWater(bool yesno)
{ enableVaporization_ = yesno; }
/*!
* \brief Finish initializing the oil phase PVT properties.
*/
@@ -110,15 +127,28 @@ public:
/*!
* \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters.
*
*/
template <class Evaluation>
Evaluation internalEnergy(unsigned,
Evaluation internalEnergy(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& /*rv*/,
const Evaluation& /*rvw*/) const
const Evaluation& rv,
const Evaluation& rvw) const
{
return H2::gasInternalEnergy(temperature, pressure, extrapolate);
/* NOTE: Assume ideal mixing */
// Init output
Evaluation result = 0;
// We have to check that one of RV and RVW is zero since H2STORE works with either GAS/WATER or GAS/OIL system
assert(rv == 0.0 || rvw == 0.0);
// Calculate each component contribution and return weighted sum
const Evaluation xBrine = convertRvwToXgW_(max(rvw, rv), regionIdx);
result += xBrine * H2O::gasInternalEnergy(temperature, pressure);
result += (1 - xBrine) * H2::gasInternalEnergy(temperature, pressure, extrapolate);
return result;
}
/*!
@@ -130,7 +160,9 @@ public:
const Evaluation& pressure,
const Evaluation& /*Rv*/,
const Evaluation& /*Rvw*/) const
{ return saturatedViscosity(regionIdx, temperature, pressure); }
{
return saturatedViscosity(regionIdx, temperature, pressure);
}
/*!
* \brief Returns the dynamic viscosity [Pa s] of oil saturated gas at given pressure.
@@ -140,7 +172,7 @@ public:
const Evaluation& temperature,
const Evaluation& pressure) const
{
return H2::gasViscosity(temperature, pressure);
return H2::gasViscosity(temperature, pressure, extrapolate);
}
/*!
@@ -150,9 +182,24 @@ public:
Evaluation inverseFormationVolumeFactor(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& /*Rv*/,
const Evaluation& /*Rvw*/) const
{ return saturatedInverseFormationVolumeFactor(regionIdx, temperature, pressure); }
const Evaluation& rv,
const Evaluation& rvw) const
{
// If vaporization is disabled, return H2 gas volume factor
if (!enableVaporization_)
return H2::gasDensity(temperature, pressure, extrapolate)/gasReferenceDensity_[regionIdx];
/* NOTE: Assume ideal mixing! */
// We have to check that one of RV and RVW is zero since H2STORE works with either GAS/WATER or GAS/OIL system
assert(rv == 0.0 || rvw == 0.0);
// Calculate each component contribution and return weighted sum
const Evaluation xBrine = convertRvwToXgW_(max(rvw, rv),regionIdx);
const auto& rhoH2 = H2::gasDensity(temperature, pressure, extrapolate);
const auto& rhoH2O = H2O::gasDensity(temperature, pressure);
return 1.0 / ((xBrine / rhoH2O + (1.0 - xBrine) / rhoH2) * gasReferenceDensity_[regionIdx]);
}
/*!
* \brief Returns the formation volume factor [-] of oil saturated gas at given pressure.
@@ -162,7 +209,8 @@ public:
const Evaluation& temperature,
const Evaluation& pressure) const
{
return H2::gasDensity(temperature, pressure, extrapolate)/gasReferenceDensity_[regionIdx];
const Evaluation rvw = rvwSat_(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
return inverseFormationVolumeFactor(regionIdx, temperature, pressure, Evaluation(0.0), rvw);
}
/*!
@@ -174,46 +222,55 @@ public:
Evaluation saturationPressure(unsigned /*regionIdx*/,
const Evaluation& /*temperature*/,
const Evaluation& /*Rv*/) const
{ return 0.0; /* this is dry gas! */ }
{ return 0.0; /* Not implemented! */ }
/*!
* \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! */ }
Evaluation saturatedWaterVaporizationFactor(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure) const
{
return rvwSat_(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
}
/*!
* \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; }
Evaluation saturatedWaterVaporizationFactor(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& saltConcentration) const
{
const Evaluation salinity = salinityFromConcentration(temperature, pressure, saltConcentration);
return rvwSat_(regionIdx, temperature, pressure, salinity);
}
/*!
* \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*/,
Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure,
const Evaluation& /*oilSaturation*/,
const Evaluation& /*maxOilSaturation*/) const
{ return 0.0; /* this is dry gas! */ }
{
return rvwSat_(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
}
/*!
* \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! */ }
Evaluation saturatedOilVaporizationFactor(unsigned regionIdx,
const Evaluation& temperature,
const Evaluation& pressure) const
{
return rvwSat_(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
}
template <class Evaluation>
Evaluation diffusionCoefficient(const Evaluation& temperature,
@@ -226,8 +283,85 @@ public:
const Scalar gasReferenceDensity(unsigned regionIdx) const
{ return gasReferenceDensity_[regionIdx]; }
Scalar oilReferenceDensity(unsigned regionIdx) const
{ return brineReferenceDensity_[regionIdx]; }
Scalar waterReferenceDensity(unsigned regionIdx) const
{ return brineReferenceDensity_[regionIdx]; }
Scalar salinity(unsigned regionIdx) const
{ return salinity_[regionIdx]; }
private:
std::vector<Scalar> gasReferenceDensity_;
std::vector<Scalar> brineReferenceDensity_;
std::vector<Scalar> salinity_;
bool enableVaporization_ = true;
template <class LhsEval>
LhsEval rvwSat_(unsigned regionIdx,
const LhsEval& temperature,
const LhsEval& pressure,
const LhsEval& salinity) const
{
// If water vaporization is disabled, we return zero
if (!enableVaporization_)
return 0.0;
// From Li et al., Int. J. Hydrogen Energ., 2018, water mole fraction is calculated assuming ideal mixing
LhsEval pw_sat = H2O::vaporPressure(temperature);
LhsEval yH2O = pw_sat / pressure;
// normalize the phase compositions
yH2O = max(0.0, min(1.0, yH2O));
return convertXgWToRvw(convertxgWToXgW(yH2O, salinity), regionIdx);
}
/*!
* \brief Convert the mass fraction of the water component in the gas phase to the
* corresponding water vaporization factor.
*/
template <class LhsEval>
LhsEval convertXgWToRvw(const LhsEval& XgW, unsigned regionIdx) const
{
Scalar rho_wRef = brineReferenceDensity_[regionIdx];
Scalar rho_gRef = gasReferenceDensity_[regionIdx];
return XgW/(1.0 - XgW)*(rho_gRef/rho_wRef);
}
/*!
* \brief Convert a water vaporization factor to the the corresponding mass fraction
* of the water component in the gas phase.
*/
template <class LhsEval>
LhsEval convertRvwToXgW_(const LhsEval& Rvw, unsigned regionIdx) const
{
Scalar rho_wRef = brineReferenceDensity_[regionIdx];
Scalar rho_gRef = gasReferenceDensity_[regionIdx];
const LhsEval& rho_wG = Rvw*rho_wRef;
return rho_wG/(rho_gRef + rho_wG);
}
/*!
* \brief Convert a water mole fraction in the gas phase the corresponding mass fraction.
*/
template <class LhsEval>
LhsEval convertxgWToXgW(const LhsEval& xgW, const LhsEval& salinity) const
{
Scalar M_H2 = H2::molarMass();
LhsEval M_Brine = Brine::molarMass(salinity);
return xgW*M_Brine / (xgW*(M_Brine - M_H2) + M_H2);
}
template <class LhsEval>
const LhsEval salinityFromConcentration(const LhsEval&T, const LhsEval& P, const LhsEval& saltConcentration) const
{
return saltConcentration / H2O::liquidDensity(T, P, true);
}
}; // end class H2GasPvt
} // end namspace Opm

0
src/opm/input/eclipse/share/keywords/001_Eclipse300/H/H2STORE → src/opm/input/eclipse/share/keywords/900_OPM/H/H2STORE
View File

2
src/opm/input/eclipse/share/keywords/keyword_list.cmake
View File

@@ -1061,7 +1061,6 @@ set( keywords
001_Eclipse300/G/GASVISCT
001_Eclipse300/G/GASWAT
001_Eclipse300/G/GSF
001_Eclipse300/H/H2STORE
001_Eclipse300/H/HEATCR
001_Eclipse300/H/HEATCRT
001_Eclipse300/H/HWELLS
@@ -1115,6 +1114,7 @@ set( keywords
900_OPM/G/GCOMPIDX
900_OPM/G/GASDENT
900_OPM/G/GASJT
900_OPM/H/H2STORE
900_OPM/M/MECH
900_OPM/M/MICP
900_OPM/M/MICPPARA

27
src/opm/material/components/H2.cpp
View File

@@ -23,17 +23,40 @@
#include <config.h>
#include <opm/material/components/H2.hpp>
#if HAVE_QUAD
#include <opm/material/common/quad.hpp>
#endif
#include "h2tables.inc"
namespace Opm {
template<>
const UniformTabulated2DFunction<double>&
H2<double>::tabulatedEnthalpy = H2Tables::tabulatedEnthalpy;
template<>
const UniformTabulated2DFunction<double>&
H2<double>::tabulatedDensity = H2Tables::tabulatedDensity;
template<>
const double H2<double>::brineSalinity = H2Tables::brineSalinity;
template<>
const UniformTabulated2DFunction<double>&
H2<float>::tabulatedEnthalpy = H2Tables::tabulatedEnthalpy;
template<>
const UniformTabulated2DFunction<double>&
H2<float>::tabulatedDensity = H2Tables::tabulatedDensity;
template<>
const float H2<float>::brineSalinity = H2Tables::brineSalinity;
} // namespace Opm
#if HAVE_QUAD
template<>
const UniformTabulated2DFunction<double>&
H2<quad>::tabulatedEnthalpy = H2Tables::tabulatedEnthalpy;
template<>
const UniformTabulated2DFunction<double>&
H2<quad>::tabulatedDensity = H2Tables::tabulatedDensity;
template<>
const quad H2<quad>::brineSalinity = H2Tables::brineSalinity;
#endif
} // namespace Opm

61280
src/opm/material/components/h2tables.inc
View File

File diff suppressed because it is too large Load Diff

12
src/opm/material/fluidsystems/BlackOilFluidSystem.cpp
View File

@@ -137,9 +137,15 @@ initFromState(const EclipseState& eclState, const Schedule& schedule)
// Use molar mass of H2 and Brine as default in H2STORE keyword
if (eclState.runspec().h2Storage()) {
// Salinity in mass fraction
const Scalar molality = eclState.getTableManager().salinity(); // mol/kg
const Scalar MmNaCl = 58.44e-3; // molar mass of NaCl [kg/mol]
const Scalar salinity = 1 / ( 1 + 1 / (molality*MmNaCl));
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();
molarMass_[regionIdx][oilCompIdx] = BrineH2Pvt<Scalar>::Brine::molarMass(salinity);
if (phaseIsActive(waterPhaseIdx))
molarMass_[regionIdx][waterCompIdx] = BrineH2Pvt<Scalar>::Brine::molarMass(salinity);
if (!phaseIsActive(gasPhaseIdx)) {
OPM_THROW(std::runtime_error,
"H2STORE requires gas phase\n");
@@ -176,13 +182,13 @@ initFromState(const EclipseState& eclState, const Schedule& schedule)
"or set it to zero.");
}
}
} else if ( eclState.runspec().co2Storage()
} else if ( (eclState.runspec().co2Storage() || eclState.runspec().h2Storage())
&& eclState.runspec().phases().active(Phase::GAS)
&& eclState.runspec().phases().active(Phase::WATER))
{
diffusionCoefficients_.resize(numRegions,{0,0,0,0,0,0,0,0,0});
// diffusion coefficients can be set using DIFFCGAS and DIFFCWAT
// for CO2STORE cases with gas + water
// for CO2STORE and H2STORE cases with gas + water
const auto& diffCoeffWatTables = eclState.getTableManager().getDiffusionCoefficientWaterTable();
if (!diffCoeffWatTables.empty()) {
for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {

10
src/opm/material/fluidsystems/blackoilpvt/BrineH2Pvt.cpp
View File

@@ -50,6 +50,9 @@ initFromState(const EclipseState& eclState, const Schedule&)
// Check if DISGAS has been activated (enables H2 dissolved in brine)
setEnableDissolvedGas(eclState.getSimulationConfig().hasDISGASW() || eclState.getSimulationConfig().hasDISGAS());
// Check if BRINE has been activated (varying salt concentration in brine)
setEnableSaltConcentration(eclState.runspec().phases().active(Phase::BRINE));
// We only supported single pvt region for the H2-brine module
size_t numRegions = 1;
setNumRegions(numRegions);
@@ -57,15 +60,14 @@ initFromState(const EclipseState& eclState, const Schedule&)
// 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
const Scalar MmNaCl = 58.44e-3; // molar mass of NaCl [kg/mol]
salinity_[regionIdx] = 1 / ( 1 + 1 / (molality*MmNaCl)); // convert to mass fraction
// 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);
brineReferenceDensity_[regionIdx] = Brine::liquidDensity(T_ref, P_ref, salinity_[regionIdx], extrapolate);
h2ReferenceDensity_[regionIdx] = H2::gasDensity(T_ref, P_ref, extrapolate);
}

4
src/opm/material/fluidsystems/blackoilpvt/H2GasPvt.cpp
View File

@@ -46,6 +46,9 @@ initFromState(const EclipseState& eclState, const Schedule&)
"H2 pvt properties are calculated based on ideal gas relations, "
"and PVDG/PVTG input is ignored.");
}
// Enable vaporization of water if needed
setEnableVaporizationWater(eclState.getSimulationConfig().hasVAPOIL() ||
eclState.getSimulationConfig().hasVAPWAT());
// We only supported single pvt region for the H2-brine module
size_t numRegions = 1;
@@ -54,6 +57,7 @@ initFromState(const EclipseState& eclState, const Schedule&)
Scalar T_ref = eclState.getTableManager().stCond().temperature;
Scalar P_ref = eclState.getTableManager().stCond().pressure;
gasReferenceDensity_[regionIdx] = H2::gasDensity(T_ref, P_ref, extrapolate);
brineReferenceDensity_[regionIdx] = Brine::liquidDensity(T_ref, P_ref, salinity_[regionIdx], extrapolate);
initEnd();
}

12322
tests/material/h2_unittest.json
View File

File diff suppressed because it is too large Load Diff

5
tests/material/test_components.cpp
View File

@@ -703,6 +703,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(H2Class, Scalar, Types)
// Rel. diff. tolerance
Scalar tol = 1e-2;
Scalar tol_visc = 2.6e-2;
Scalar tol_enth = 3.8e-2;
// Loop over temperature and pressure, and compare to reference values in JSON file
for (int iT = 0; iT < numT; ++iT) {
@@ -736,9 +737,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(H2Class, Scalar, Types)
Json::JsonObject enth_ref_row = enthalpy_ref.get_array_item(iT);
Scalar enth_ref = Scalar(enth_ref_row.get_array_item(iP).as_double());
BOOST_CHECK_MESSAGE(close_at_tolerance(enthalpy, enth_ref, tol),
BOOST_CHECK_MESSAGE(close_at_tolerance(enthalpy, enth_ref, tol_enth),
"relative difference between enthalpy {"<<enthalpy<<"} and reference {"<<enth_ref<<
"} exceeds tolerance {"<<tol<<"} at (T, p) = ("<<T.value()<<", "<<p.value()<<")");
"} exceeds tolerance {"<<tol_enth<<"} at (T, p) = ("<<T.value()<<", "<<p.value()<<")");
}
}
}

283
tests/material/test_h2brinepvt.cpp Normal file
View File

@@ -0,0 +1,283 @@
// -*- 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
*
* \brief This is the unit test for the H2-brine PVT model
*
* This test requires the presence of opm-common.
*/
#include "config.h"
#define BOOST_TEST_MODULE H2BrinePvt
#include <boost/test/unit_test.hpp>
#if !HAVE_ECL_INPUT
#error "The test for the H2-brine PVT classes requires eclipse input support in opm-common"
#endif
#include <opm/material/fluidsystems/blackoilpvt/GasPvtMultiplexer.hpp>
#include <opm/material/fluidsystems/blackoilpvt/OilPvtMultiplexer.hpp>
#include <opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/input/eclipse/Parser/Parser.hpp>
#include <opm/input/eclipse/Deck/Deck.hpp>
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/Python/Python.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <iostream>
// values of strings based on the first SPE1 test case of opm-data. note that in the
// real world it does not make much sense to specify a fluid phase using more than a
// single keyword, but for a unit test, this saves a lot of boiler-plate code.
constexpr const char* deckString1 =
"RUNSPEC\n"
"\n"
"DIMENS\n"
" 10 10 3 /\n"
"\n"
"TABDIMS\n"
" * 1 /\n"
"\n"
"OIL\n"
"GAS\n"
"H2STORE\n"
"\n"
"DISGAS\n"
"\n"
"METRIC\n"
"\n"
"GRID\n"
"\n"
"DX\n"
" 300*1000 /\n"
"DY\n"
" 300*1000 /\n"
"DZ\n"
" 100*20 100*30 100*50 /\n"
"\n"
"TOPS\n"
" 100*1234 /\n"
"\n"
"PORO\n"
" 300*0.15 /\n"
"PROPS\n"
"\n";
constexpr const char* deckString2 =
"RUNSPEC\n"
"\n"
"DIMENS\n"
" 10 10 3 /\n"
"\n"
"TABDIMS\n"
" * 1 /\n"
"\n"
"WATER\n"
"GAS\n"
"H2STORE\n"
"\n"
"DISGASW\n"
"\n"
"METRIC\n"
"\n"
"GRID\n"
"\n"
"DX\n"
" 300*1000 /\n"
"DY\n"
" 300*1000 /\n"
"DZ\n"
" 100*20 100*30 100*50 /\n"
"\n"
"TOPS\n"
" 100*1234 /\n"
"\n"
"PORO\n"
" 300*0.15 /\n"
"PROPS\n"
"\n";
template <class Evaluation, class BrinePvt>
void ensurePvtApiBrine(const BrinePvt& brinePvt)
{
// we don't want to run this, we just want to make sure that it compiles
while (0) {
Evaluation temperature = 273.15 + 20.0;
Evaluation pressure = 1e5;
Evaluation saltconcentration = 0.0;
Evaluation rs = 0.0;
////
// Water PVT API
/////
std::cout << brinePvt.viscosity(/*regionIdx=*/0,
temperature,
pressure,
rs,
saltconcentration);
std::cout << brinePvt.inverseFormationVolumeFactor(/*regionIdx=*/0,
temperature,
pressure,
rs,
saltconcentration);
std::cout << brinePvt.internalEnergy(/*regionIdx=*/0,
temperature,
pressure,
rs,
saltconcentration);
}
}
template <class Evaluation, class H2Pvt>
void ensurePvtApiGas(const H2Pvt& h2Pvt)
{
// we don't want to run this, we just want to make sure that it compiles
while (0) {
Evaluation temperature = 273.15 + 20.0;
Evaluation pressure = 1e5;
Evaluation Rv = 0.0;
Evaluation Rvw = 0.0;
Evaluation So = 0.5;
Evaluation maxSo = 1.0;
/////
// H2 PVT API
/////
std::cout << h2Pvt.viscosity(/*regionIdx=*/0,
temperature,
pressure,
Rv,
Rvw);
std::cout << h2Pvt.inverseFormationVolumeFactor(/*regionIdx=*/0,
temperature,
pressure,
Rv,
Rvw);
std::cout << h2Pvt.saturatedViscosity(/*regionIdx=*/0,
temperature,
pressure);
std::cout << h2Pvt.saturatedInverseFormationVolumeFactor(/*regionIdx=*/0,
temperature,
pressure);
std::cout << h2Pvt.saturationPressure(/*regionIdx=*/0,
temperature,
Rv);
std::cout << h2Pvt.saturatedOilVaporizationFactor(/*regionIdx=*/0,
temperature,
pressure);
std::cout << h2Pvt.saturatedOilVaporizationFactor(/*regionIdx=*/0,
temperature,
pressure,
So,
maxSo);
}
}
template <class Evaluation, class BrinePvt>
void ensurePvtApiBrineOil(const BrinePvt& brinePvt)
{
// we don't want to run this, we just want to make sure that it compiles
while (0) {
Evaluation temperature = 273.15 + 20.0;
Evaluation pressure = 1e5;
Evaluation Rs = 0.0;
Evaluation So = 0.5;
Evaluation maxSo = 1.0;
/////
// brine PVT API
/////
std::cout << brinePvt.viscosity(/*regionIdx=*/0,
temperature,
pressure,
Rs);
std::cout << brinePvt.inverseFormationVolumeFactor(/*regionIdx=*/0,
temperature,
pressure,
Rs);
std::cout << brinePvt.saturatedViscosity(/*regionIdx=*/0,
temperature,
pressure);
std::cout << brinePvt.saturatedInverseFormationVolumeFactor(/*regionIdx=*/0,
temperature,
pressure);
std::cout << brinePvt.saturationPressure(/*regionIdx=*/0,
temperature,
Rs);
std::cout << brinePvt.saturatedGasDissolutionFactor(/*regionIdx=*/0,
temperature,
pressure);
std::cout << brinePvt.saturatedGasDissolutionFactor(/*regionIdx=*/0,
temperature,
pressure,
So,
maxSo);
}
}
using Types = std::tuple<float,double>;
BOOST_AUTO_TEST_CASE_TEMPLATE(Oil, Scalar, Types)
{
Opm::Parser parser;
auto python = std::make_shared<Opm::Python>();
auto deck1 = parser.parseString(deckString1);
Opm::EclipseState eclState1(deck1);
Opm::Schedule schedule1(deck1, eclState1, python);
Opm::GasPvtMultiplexer<Scalar> h2Pvt_oil;
Opm::OilPvtMultiplexer<Scalar> brinePvt_oil;
BOOST_CHECK_NO_THROW(h2Pvt_oil.initFromState(eclState1, schedule1));
BOOST_CHECK_NO_THROW(brinePvt_oil.initFromState(eclState1, schedule1));
using Eval = Opm::DenseAd::Evaluation<Scalar, 1>;
ensurePvtApiGas<Scalar>(h2Pvt_oil);
ensurePvtApiBrineOil<Eval>(brinePvt_oil);
}
BOOST_AUTO_TEST_CASE_TEMPLATE(Water, Scalar, Types) {
Opm::Parser parser;
auto python = std::make_shared<Opm::Python>();
auto deck2 = parser.parseString(deckString2);
Opm::EclipseState eclState2(deck2);
Opm::Schedule schedule2(deck2, eclState2, python);
Opm::GasPvtMultiplexer<Scalar> h2Pvt;
Opm::WaterPvtMultiplexer<Scalar> brinePvt;
BOOST_CHECK_NO_THROW(h2Pvt.initFromState(eclState2, schedule2));
BOOST_CHECK_NO_THROW(brinePvt.initFromState(eclState2, schedule2));
using Eval = Opm::DenseAd::Evaluation<Scalar, 1>;
ensurePvtApiGas<Scalar>(h2Pvt);
ensurePvtApiBrine<Eval>(brinePvt);
}