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

Merge pull request #443 from totto82/blackoilDiff

Browse Source 
Add support for diffusion in the blackoil fluid system
This commit is contained in:
Bård Skaflestad 2021-02-18 18:05:15 +01:00 committed by GitHub
commit c7317befa7
13 changed files with 223 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

101
opm/material/fluidsystems/BlackOilFluidSystem.hpp
View File

@ -243,7 +243,43 @@ public:
regionIdx);
}
// set default molarMass and mappings
initEnd();
// use molarMass of CO2 and Brine as default
// when we are using the the CO2STORE option
// NB the oil component is used internally for
// brine
if (eclState.runspec().co2Storage()) {
for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
molarMass_[regionIdx][oilCompIdx] = BrineCo2Pvt<Scalar>::Brine::molarMass();
molarMass_[regionIdx][gasCompIdx] = BrineCo2Pvt<Scalar>::CO2::molarMass();
}
}
setEnableDiffusion(eclState.getSimulationConfig().isDiffusive());
if(enableDiffusion()) {
const auto& diffCoeffTables = eclState.getTableManager().getDiffusionCoefficientTable();
if(!diffCoeffTables.empty()) {
// if diffusion coefficient table is empty we relay on the PVT model to
// to give us the coefficients.
diffusionCoefficients_.resize(numRegions,{0,0,0,0,0,0,0,0,0});
assert(diffCoeffTables.size() == numRegions);
for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const auto& diffCoeffTable = diffCoeffTables[regionIdx];
molarMass_[regionIdx][oilCompIdx] = diffCoeffTable.oil_mw;
molarMass_[regionIdx][gasCompIdx] = diffCoeffTable.gas_mw;
setDiffusionCoefficient(diffCoeffTable.gas_in_gas, gasCompIdx, gasPhaseIdx, regionIdx);
setDiffusionCoefficient(diffCoeffTable.oil_in_gas, oilCompIdx, gasPhaseIdx, regionIdx);
setDiffusionCoefficient(diffCoeffTable.gas_in_oil, gasCompIdx, oilPhaseIdx, regionIdx);
setDiffusionCoefficient(diffCoeffTable.oil_in_oil, oilCompIdx, oilPhaseIdx, regionIdx);
if(diffCoeffTable.gas_in_oil_cross_phase > 0 || diffCoeffTable.oil_in_oil_cross_phase > 0) {
throw std::runtime_error("Cross phase diffusion is set in the deck, but not implemented in Flow. "
"Please default DIFFC item 7 and item 8 or set it to zero.");
}
}
}
}
}
#endif // HAVE_ECL_INPUT
@ -261,6 +297,7 @@ public:
enableDissolvedGas_ = true;
enableVaporizedOil_ = false;
enableDiffusion_ = false;
oilPvt_ = nullptr;
gasPvt_ = nullptr;
@ -294,6 +331,15 @@ public:
static void setEnableVaporizedOil(bool yesno)
{ enableVaporizedOil_ = yesno; }
/*!
* \brief Specify whether the fluid system should consider diffusion
*
* By default, diffusion is not considered.
*/
static void setEnableDiffusion(bool yesno)
{ enableDiffusion_ = yesno; }
/*!
* \brief Set the pressure-volume-saturation (PVT) relations for the gas phase.
*/
@ -329,6 +375,7 @@ public:
referenceDensity_[regionIdx][gasPhaseIdx] = rhoGas;
}
/*!
* \brief Finish initializing the black oil fluid system.
*/
@ -543,6 +590,14 @@ public:
static bool enableVaporizedOil()
{ return enableVaporizedOil_; }
/*!
* \brief Returns whether the fluid system should consider diffusion
*
* By default, diffusion is not considered.
*/
static bool enableDiffusion()
{ return enableDiffusion_; }
/*!
* \brief Returns the density of a fluid phase at surface pressure [kg/m^3]
*
@ -1274,6 +1329,43 @@ public:
return canonicalToActivePhaseIdx_[phaseIdx];
}
//! \copydoc BaseFluidSystem::diffusionCoefficient
static Scalar diffusionCoefficient(unsigned compIdx, unsigned phaseIdx, unsigned regionIdx = 0)
{ return diffusionCoefficients_[regionIdx][numPhases*compIdx + phaseIdx]; }
//! \copydoc BaseFluidSystem::setDiffusionCoefficient
static void setDiffusionCoefficient(Scalar coefficient, unsigned compIdx, unsigned phaseIdx, unsigned regionIdx = 0)
{ diffusionCoefficients_[regionIdx][numPhases*compIdx + phaseIdx] = coefficient ; }
/*!
* \copydoc BaseFluidSystem::diffusionCoefficient
*/
template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
static LhsEval diffusionCoefficient(const FluidState& fluidState,
const ParameterCache<ParamCacheEval>& paramCache,
unsigned phaseIdx,
unsigned compIdx)
{
// diffusion is disabled by the user
if(!enableDiffusion())
return 0.0;
// diffusion coefficients are set, and we use them
if(!diffusionCoefficients_.empty()) {
return diffusionCoefficient(compIdx, phaseIdx, paramCache.regionIndex());
}
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
switch (phaseIdx) {
case oilPhaseIdx: return oilPvt().diffusionCoefficient(T, p, compIdx);
case gasPhaseIdx: return gasPvt().diffusionCoefficient(T, p, compIdx);
case waterPhaseIdx: return 0.0;
default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
}
}
private:
static void resizeArrays_(size_t numRegions)
{
@ -1289,12 +1381,14 @@ private:
static bool enableDissolvedGas_;
static bool enableVaporizedOil_;
static bool enableDiffusion_;
// HACK for GCC 4.4: the array size has to be specified using the literal value '3'
// here, because GCC 4.4 seems to be unable to determine the number of phases from
// the BlackOil fluid system in the attribute declaration below...
static std::vector<std::array<Scalar, /*numPhases=*/3> > referenceDensity_;
static std::vector<std::array<Scalar, /*numComponents=*/3> > molarMass_;
static std::vector<std::array<Scalar, /*numComponents=*/3 * /*numPhases=*/3> > diffusionCoefficients_;
static std::array<short, numPhases> activeToCanonicalPhaseIdx_;
static std::array<short, numPhases> canonicalToActivePhaseIdx_;
@ -1332,6 +1426,9 @@ bool BlackOilFluidSystem<Scalar, IndexTraits>::enableDissolvedGas_;
template <class Scalar, class IndexTraits>
bool BlackOilFluidSystem<Scalar, IndexTraits>::enableVaporizedOil_;
template <class Scalar, class IndexTraits>
bool BlackOilFluidSystem<Scalar, IndexTraits>::enableDiffusion_;
template <class Scalar, class IndexTraits>
std::shared_ptr<OilPvtMultiplexer<Scalar> >
BlackOilFluidSystem<Scalar, IndexTraits>::oilPvt_;
@ -1352,6 +1449,10 @@ template <class Scalar, class IndexTraits>
std::vector<std::array<Scalar, 3> >
BlackOilFluidSystem<Scalar, IndexTraits>::molarMass_;
template <class Scalar, class IndexTraits>
std::vector<std::array<Scalar, 9> >
BlackOilFluidSystem<Scalar, IndexTraits>::diffusionCoefficients_;
template <class Scalar, class IndexTraits>
bool BlackOilFluidSystem<Scalar, IndexTraits>::isInitialized_ = false;

31
opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.hpp
View File

@ -57,8 +57,6 @@ template <class Scalar>
class BrineCo2Pvt
{
typedef std::vector<std::pair<Scalar, Scalar> > SamplingPoints;
typedef Opm::SimpleHuDuanH2O<Scalar> H2O;
typedef Opm::Brine<Scalar, H2O> Brine;
//typedef Opm::H2O<Scalar> H2O_IAPWS;
//typedef Opm::Brine<Scalar, H2O_IAPWS> Brine_IAPWS;
@ -68,9 +66,12 @@ class BrineCo2Pvt
//typedef H2O_Tabulated H2O;
//typedef Brine_Tabulated Brine;
typedef Opm::CO2<Scalar, CO2Tables> CO2;
public:
typedef Opm::SimpleHuDuanH2O<Scalar> H2O;
typedef Opm::Brine<Scalar, H2O> Brine;
typedef Opm::CO2<Scalar, CO2Tables> CO2;
typedef Opm::Tabulated1DFunction<Scalar> TabulatedOneDFunction;
//! The binary coefficients for brine and CO2 used by this fluid system
@ -276,6 +277,22 @@ public:
brineReferenceDensity_ == data.brineReferenceDensity_;
}
template <class Evaluation>
Evaluation diffusionCoefficient(const Evaluation& temperature,
const Evaluation& pressure,
unsigned /*compIdx*/) const
{
//Diffusion coefficient of CO2 in pure water according to (McLachlan and Danckwerts, 1972)
const Evaluation log_D_H20 = -4.1764 + 712.52 / temperature - 2.5907e5 / (temperature*temperature);
//Diffusion coefficient of CO2 in the brine phase modified following (Ratcliff and Holdcroft,1963 and Al-Rawajfeh, 2004)
const Evaluation& mu_H20 = H2O::liquidViscosity(temperature, pressure); // Water viscosity
const Evaluation& mu_Brine = Brine::liquidViscosity(temperature, pressure); // Brine viscosity
const Evaluation log_D_Brine = log_D_H20 - 0.87*Opm::log10(mu_Brine / mu_H20);
return Opm::pow(Evaluation(10), log_D_Brine) * 1e-4; // convert from cm2/s to m2/s
}
private:
std::vector<Scalar> brineReferenceDensity_;
std::vector<Scalar> co2ReferenceDensity_;
@ -372,8 +389,8 @@ private:
LhsEval convertXoGToxoG_(const LhsEval& XoG) const
{
Scalar M_CO2 = CO2::molarMass();
Scalar M_H2O = H2O::molarMass();
return XoG*M_H2O / (M_CO2*(1 - XoG) + XoG*M_H2O);
Scalar M_Brine = Brine::molarMass();
return XoG*M_Brine / (M_CO2*(1 - XoG) + XoG*M_Brine);
}
@ -384,9 +401,9 @@ private:
LhsEval convertxoGToXoG(const LhsEval& xoG) const
{
Scalar M_CO2 = CO2::molarMass();
Scalar M_H2O = H2O::molarMass();
Scalar M_Brine = Brine::molarMass();
return xoG*M_CO2 / (xoG*(M_CO2 - M_H2O) + M_H2O);
return xoG*M_CO2 / (xoG*(M_CO2 - M_Brine) + M_Brine);
}

14
opm/material/fluidsystems/blackoilpvt/Co2GasPvt.hpp
View File

@ -31,7 +31,9 @@
#include <opm/material/common/Tabulated1DFunction.hpp>
#include <opm/material/components/CO2.hpp>
#include <opm/material/components/SimpleHuDuanH2O.hpp>
#include <opm/material/common/UniformTabulated2DFunction.hpp>
#include <opm/material/binarycoefficients/Brine_CO2.hpp>
#include <opm/material/components/co2tables.inc>
#if HAVE_ECL_INPUT
@ -52,10 +54,14 @@ class Co2GasPvt
{
typedef std::vector<std::pair<Scalar, Scalar> > SamplingPoints;
typedef Opm::CO2<Scalar, CO2Tables> CO2;
typedef Opm::SimpleHuDuanH2O<Scalar> H2O;
public:
typedef Opm::Tabulated1DFunction<Scalar> TabulatedOneDFunction;
//! The binary coefficients for brine and CO2 used by this fluid system
typedef Opm::BinaryCoeff::Brine_CO2<Scalar, H2O, CO2> BinaryCoeffBrineCO2;
explicit Co2GasPvt() = default;
Co2GasPvt(const std::vector<Scalar>& gasReferenceDensity)
: gasReferenceDensity_(gasReferenceDensity)
@ -205,6 +211,14 @@ public:
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 BinaryCoeffBrineCO2::gasDiffCoeff(temperature, pressure);
}
const Scalar gasReferenceDensity(unsigned regionIdx) const
{ return gasReferenceDensity_[regionIdx]; }

8
opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp
View File

@ -272,6 +272,14 @@ public:
const Evaluation& /*Rs*/) const
{ return 0.0; /* this is dead oil, so there isn't any meaningful saturation pressure! */ }
template <class Evaluation>
Evaluation diffusionCoefficient(const Evaluation& /*temperature*/,
const Evaluation& /*pressure*/,
unsigned /*compIdx*/) const
{
throw std::runtime_error("Not implemented: The PVT model does not provide a diffusionCoefficient()");
}
const Scalar oilReferenceDensity(unsigned regionIdx) const
{ return oilReferenceDensity_[regionIdx]; }

8
opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp
View File

@ -278,6 +278,14 @@ public:
const Evaluation& /*pressure*/) const
{ return 0.0; /* this is dead oil! */ }
template <class Evaluation>
Evaluation diffusionCoefficient(const Evaluation& /*temperature*/,
const Evaluation& /*pressure*/,
unsigned /*compIdx*/) const
{
throw std::runtime_error("Not implemented: The PVT model does not provide a diffusionCoefficient()");
}
const Scalar oilReferenceDensity(unsigned regionIdx) const
{ return oilReferenceDensity_[regionIdx]; }

8
opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp
View File

@ -288,6 +288,14 @@ public:
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
{
throw std::runtime_error("Not implemented: The PVT model does not provide a diffusionCoefficient()");
}
const Scalar gasReferenceDensity(unsigned regionIdx) const
{ return gasReferenceDensity_[regionIdx]; }

11
opm/material/fluidsystems/blackoilpvt/GasPvtMultiplexer.hpp
View File

@ -272,6 +272,17 @@ public:
const Evaluation& Rv) const
{ OPM_GAS_PVT_MULTIPLEXER_CALL(return pvtImpl.saturationPressure(regionIdx, temperature, Rv)); return 0; }
/*!
* \copydoc BaseFluidSystem::diffusionCoefficient
*/
template <class Evaluation>
Evaluation diffusionCoefficient(const Evaluation& temperature,
const Evaluation& pressure,
unsigned compIdx) const
{
OPM_GAS_PVT_MULTIPLEXER_CALL(return pvtImpl.diffusionCoefficient(temperature, pressure, compIdx)); return 0;
}
/*!
* \brief Returns the concrete approach for calculating the PVT relations.
*

7
opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp
View File

@ -367,6 +367,13 @@ public:
const Evaluation& pressure) const
{ return isothermalPvt_->saturationPressure(regionIdx, temperature, pressure); }
template <class Evaluation>
Evaluation diffusionCoefficient(const Evaluation& temperature,
const Evaluation& pressure,
unsigned compIdx) const
{
return isothermalPvt_->diffusionCoefficient(temperature, pressure, compIdx);
}
const IsothermalPvt* isoThermalPvt() const
{ return isothermalPvt_; }

7
opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp
View File

@ -602,6 +602,13 @@ public:
throw NumericalIssue(errlog.str());
}
template <class Evaluation>
Evaluation diffusionCoefficient(const Evaluation& /*temperature*/,
const Evaluation& /*pressure*/,
unsigned /*compIdx*/) const
{
throw std::runtime_error("Not implemented: The PVT model does not provide a diffusionCoefficient()");
}
const Scalar gasReferenceDensity(unsigned regionIdx) const
{ return gasReferenceDensity_[regionIdx]; }

11
opm/material/fluidsystems/blackoilpvt/OilPvtMultiplexer.hpp
View File

@ -266,6 +266,17 @@ public:
const Evaluation& Rs) const
{ OPM_OIL_PVT_MULTIPLEXER_CALL(return pvtImpl.saturationPressure(regionIdx, temperature, Rs)); return 0; }
/*!
* \copydoc BaseFluidSystem::diffusionCoefficient
*/
template <class Evaluation>
Evaluation diffusionCoefficient(const Evaluation& temperature,
const Evaluation& pressure,
unsigned compIdx) const
{
OPM_OIL_PVT_MULTIPLEXER_CALL(return pvtImpl.diffusionCoefficient(temperature, pressure, compIdx)); return 0;
}
void setApproach(OilPvtApproach appr)
{
switch (appr) {

8
opm/material/fluidsystems/blackoilpvt/OilPvtThermal.hpp
View File

@ -377,6 +377,14 @@ public:
const Evaluation& pressure) const
{ return isothermalPvt_->saturationPressure(regionIdx, temperature, pressure); }
template <class Evaluation>
Evaluation diffusionCoefficient(const Evaluation& temperature,
const Evaluation& pressure,
unsigned compIdx) const
{
return isothermalPvt_->diffusionCoefficient(temperature, pressure, compIdx);
}
const IsothermalPvt* isoThermalPvt() const
{ return isothermalPvt_; }

8
opm/material/fluidsystems/blackoilpvt/SolventPvt.hpp
View File

@ -189,6 +189,14 @@ public:
return invBg/invMugBg;
}
template <class Evaluation>
Evaluation diffusionCoefficient(const Evaluation& /*temperature*/,
const Evaluation& /*pressure*/,
unsigned /*compIdx*/) const
{
throw std::runtime_error("Not implemented: The PVT model does not provide a diffusionCoefficient()");
}
/*!
* \brief Return the reference density the solvent phase for a given PVT region
*/

8
opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp
View File

@ -628,6 +628,14 @@ public:
throw NumericalIssue(errlog.str());
}
template <class Evaluation>
Evaluation diffusionCoefficient(const Evaluation& /*temperature*/,
const Evaluation& /*pressure*/,
unsigned /*compIdx*/) const
{
throw std::runtime_error("Not implemented: The PVT model does not provide a diffusionCoefficient()");
}
const Scalar gasReferenceDensity(unsigned regionIdx) const
{ return gasReferenceDensity_[regionIdx]; }