OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4034 from plgbrts/equil2

Browse Source 
Equilibration with water-gas ratio vs depth (RVWVD) table
This commit is contained in:
Bård Skaflestad 2022-09-14 14:30:04 +02:00 committed by GitHub
commit 67ec01a730
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 321 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
ebos/equil/equilibrationhelpers.cc
View File

@ -195,6 +195,42 @@ satRv(const double press, const double temp) const
}
template<class FluidSystem>
RvwVD<FluidSystem>::RvwVD(const int pvtRegionIdx,
const std::vector<double>& depth,
const std::vector<double>& rvw)
: pvtRegionIdx_(pvtRegionIdx)
, rvwVsDepth_(depth, rvw)
{
}
template<class FluidSystem>
double RvwVD<FluidSystem>::
operator()(const double depth,
const double press,
const double temp,
const double satWat) const
{
if (std::abs(satWat) > 1e-16) {
return satRvw(press, temp); //saturated Rvw
}
else {
if (rvwVsDepth_.xMin() > depth)
return rvwVsDepth_.valueAt(0);
else if (rvwVsDepth_.xMax() < depth)
return rvwVsDepth_.valueAt(rvwVsDepth_.numSamples() - 1);
return std::min(satRvw(press, temp), rvwVsDepth_.eval(depth, /*extrapolate=*/false));
}
}
template<class FluidSystem>
double RvwVD<FluidSystem>::
satRvw(const double press, const double temp) const
{
return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);
}
template<class FluidSystem>
RsSatAtContact<FluidSystem>::
RsSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact)
@ -255,16 +291,48 @@ satRv(const double press, const double temp) const
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);;
}
template<class FluidSystem>
RvwSatAtContact<FluidSystem>::
RvwSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact)
: pvtRegionIdx_(pvtRegionIdx)
{
rvwSatContact_ = satRvw(pContact, T_contact);
}
template<class FluidSystem>
double RvwSatAtContact<FluidSystem>::
operator()(const double /*depth*/,
const double press,
const double temp,
const double satWat) const
{
if (satWat > 0.0) {
return satRvw(press, temp);
}
else {
return std::min(satRvw(press, temp), rvwSatContact_);
}
}
template<class FluidSystem>
double RvwSatAtContact<FluidSystem>::
satRvw(const double press, const double temp) const
{
return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);;
}
} // namespace Miscibility
EquilReg::EquilReg(const EquilRecord& rec,
std::shared_ptr<Miscibility::RsFunction> rs,
std::shared_ptr<Miscibility::RsFunction> rv,
std::shared_ptr<Miscibility::RsFunction> rvw,
const TabulatedFunction& saltVdTable,
const int pvtIdx)
: rec_ (rec)
, rs_ (rs)
, rv_ (rv)
, rvw_ (rvw)
, saltVdTable_ (saltVdTable)
, pvtIdx_ (pvtIdx)
{
@ -317,6 +385,12 @@ EquilReg::evaporationCalculator() const
return *this->rv_;
}
const EquilReg::CalcWaterEvaporation&
EquilReg::waterEvaporationCalculator() const
{
return *this->rvw_;
}
const EquilReg::TabulatedFunction&
EquilReg::saltVdTable() const
{
@ -556,7 +630,9 @@ namespace Miscibility {
template class RsVD<FS>;
template class RsSatAtContact<FS>;
template class RvSatAtContact<FS>;
template class RvwSatAtContact<FS>;
template class RvVD<FS>;
template class RvwVD<FS>;
}
} // namespace Equil

119
ebos/equil/equilibrationhelpers.hh
View File

@ -363,6 +363,56 @@ private:
};
/**
* Type that implements "vaporized water-gas ratio"
* tabulated as a function of depth policy. Data
* typically taken from keyword 'RVWVD'.
*/
template <class FluidSystem>
class RvwVD : public RsFunction
{
public:
/**
* Constructor.
*
* \param[in] pvtRegionIdx The pvt region index
* \param[in] depth Depth nodes.
* \param[in] rvw Evaporized water-gasl ratio at @c depth.
*/
RvwVD(const int pvtRegionIdx,
const std::vector<double>& depth,
const std::vector<double>& rvw);
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RVW
* value.
*
* \param[in] press Pressure at which to calculate RVW
* value.
*
* \param[in] temp Temperature at which to calculate RVW
* value.
*
* \return Vaporized water-gas ratio (RVW) at depth @c
* depth and pressure @c press.
*/
double operator()(const double depth,
const double press,
const double temp,
const double satWat = 0.0) const;
private:
using RvwVsDepthFunc = Tabulated1DFunction<double>;
const int pvtRegionIdx_;
RvwVsDepthFunc rvwVsDepth_;
double satRvw(const double press, const double temp) const;
};
/**
* Class that implements "dissolved gas-oil ratio" (Rs)
* as function of depth and pressure as follows:
@ -472,6 +522,60 @@ private:
double satRv(const double press, const double temp) const;
};
/**
* Class that implements "vaporized water-gas ratio" (Rvw)
* as function of depth and pressure as follows:
*
* 1. The Rvw at the gas-water contact is equal to the
* saturated Rv value, RvwSatContact.
*
* 2. The Rvw elsewhere is equal to RvwSatContact, but
* constrained to the saturated value as given by the
* local pressure.
*
* This should yield Rvw-values that are constant below the
* contact, and decreasing above the contact.
*/
template <class FluidSystem>
class RvwSatAtContact : public RsFunction
{
public:
/**
* Constructor.
*
* \param[in] pvtRegionIdx The pvt region index
* \param[in] pContact oil pressure at the contact
* \param[in] T_contact temperature at the contact
*/
RvwSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact);
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RVW
* value.
*
* \param[in] press Pressure at which to calculate RVW
* value.
*
* \param[in] temp Temperature at which to calculate RVW
* value.
*
* \return Dissolved water-gas ratio (RVW) at depth @c
* depth and pressure @c press.
*/
double operator()(const double /*depth*/,
const double press,
const double temp,
const double satWat = 0.0) const;
private:
const int pvtRegionIdx_;
double rvwSatContact_;
double satRvw(const double press, const double temp) const;
};
} // namespace Miscibility
/**
@ -504,11 +608,13 @@ public:
* \param[in] rec Equilibration data of current region.
* \param[in] rs Calculator of dissolved gas-oil ratio.
* \param[in] rv Calculator of vapourised oil-gas ratio.
* \param[in] rvw Calculator of vapourised water-gas ratio.
* \param[in] pvtRegionIdx The pvt region index
*/
EquilReg(const EquilRecord& rec,
std::shared_ptr<Miscibility::RsFunction> rs,
std::shared_ptr<Miscibility::RsFunction> rv,
std::shared_ptr<Miscibility::RsFunction> rvw,
const TabulatedFunction& saltVdTable,
const int pvtIdx);
@ -522,6 +628,12 @@ public:
*/
using CalcEvaporation = Miscibility::RsFunction;
/**
* Type of vapourised water-gas ratio calculator.
*/
using CalcWaterEvaporation = Miscibility::RsFunction;
/**
* Datum depth in current region
*/
@ -577,6 +689,12 @@ public:
*/
const CalcEvaporation& evaporationCalculator() const;
/**
* Retrieve vapourised water-gas ratio calculator of current
* region.
*/
const CalcWaterEvaporation& waterEvaporationCalculator() const;
const TabulatedFunction& saltVdTable() const;
/**
* Retrieve pvtIdx of the region.
@ -587,6 +705,7 @@ private:
EquilRecord rec_; /**< Equilibration data */
std::shared_ptr<Miscibility::RsFunction> rs_; /**< RS calculator */
std::shared_ptr<Miscibility::RsFunction> rv_; /**< RV calculator */
std::shared_ptr<Miscibility::RsFunction> rvw_; /**< RVW calculator */
const TabulatedFunction& saltVdTable_;
const int pvtIdx_;
};

129
ebos/equil/initstateequil.cc
View File

@ -34,6 +34,7 @@
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Tables/RsvdTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/RvvdTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/RvwvdTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PbvdTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/PdvdTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/SaltvdTable.hpp>
@ -324,29 +325,31 @@ density(const double depth,
return rho;
}
template<class FluidSystem, class RV>
Gas<FluidSystem,RV>::
template<class FluidSystem, class RV, class RVW>
Gas<FluidSystem,RV,RVW>::
Gas(const double temp,
const RV& rv,
const RVW& rvw,
const int pvtRegionIdx,
const double normGrav)
: temp_(temp)
, rv_(rv)
, rvw_(rvw)
, pvtRegionIdx_(pvtRegionIdx)
, g_(normGrav)
{
}
template<class FluidSystem, class RV>
double Gas<FluidSystem,RV>::
template<class FluidSystem, class RV, class RVW>
double Gas<FluidSystem,RV,RVW>::
operator()(const double depth,
const double press) const
{
return this->density(depth, press) * g_;
}
template<class FluidSystem, class RV>
double Gas<FluidSystem,RV>::
template<class FluidSystem, class RV, class RVW>
double Gas<FluidSystem,RV,RVW>::
density(const double depth,
const double press) const
{
@ -354,18 +357,53 @@ density(const double depth,
if (FluidSystem::enableVaporizedOil())
rv = rv_(depth, press, temp_);
double rvw = 0.0;
if (FluidSystem::enableVaporizedWater())
rvw = rvw_(depth, press, temp_);
double bGas = 0.0;
if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press)) {
bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
if (FluidSystem::enableVaporizedOil() && FluidSystem::enableVaporizedWater()) {
if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press)
&& rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp_, press))
{
bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
} else {
bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv, rvw);
}
double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_)
+ rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
return rho;
}
else {
bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv, 0.0/*=Rvw*/);
}
double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
if (FluidSystem::enableVaporizedOil()) {
if (FluidSystem::enableVaporizedOil()){
if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press)) {
bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
} else {
bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv, 0.0/*=rvw*/);
}
double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_);
return rho;
}
if (FluidSystem::enableVaporizedWater()){
if (rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp_, press)) {
bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
}
else {
bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, 0.0/*=rv*/, rvw);
}
double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
rho += rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
return rho;
}
// immiscible gas
bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, 0.0/*=rv*/, 0.0/*=rvw*/);
double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
return rho;
}
@ -1124,7 +1162,7 @@ makeGasPressure(const typename GPress::InitCond& ic,
const VSpan& span)
{
const auto drho = GasPressODE {
this->temperature_, reg.evaporationCalculator(),
this->temperature_, reg.evaporationCalculator(), reg.waterEvaporationCalculator(),
reg.pvtIdx(), this->gravity_
};
@ -1203,6 +1241,7 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
std::vector<double>(grid.size(/*codim=*/0))),
rs_(grid.size(/*codim=*/0)),
rv_(grid.size(/*codim=*/0)),
rvw_(grid.size(/*codim=*/0)),
cartesianIndexMapper_(cartMapper)
{
//Check for presence of kw SWATINIT
@ -1322,6 +1361,47 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
}
}
rvwFunc_.reserve(rec.size());
if (FluidSystem::enableVaporizedWater()) {
for (size_t i = 0; i < rec.size(); ++i) {
if (eqlmap.cells(i).empty()) {
rvwFunc_.push_back(std::shared_ptr<Miscibility::RvwVD<FluidSystem>>());
continue;
}
const int pvtIdx = regionPvtIdx_[i];
if (!rec[i].humidGasInitConstantRvw()) {
const TableContainer& rvwvdTables = tables.getRvwvdTables();
if (rvwvdTables.size() > 0) {
const RvwvdTable& rvwvdTable = rvwvdTables.getTable<RvwvdTable>(i);
std::vector<double> depthColumn = rvwvdTable.getColumn("DEPTH").vectorCopy();
std::vector<double> rvwvdColumn = rvwvdTable.getColumn("RVWVD").vectorCopy();
rvwFunc_.push_back(std::make_shared<Miscibility::RvwVD<FluidSystem>>(pvtIdx,
depthColumn, rvwvdColumn));
} else {
throw std::runtime_error("Cannot initialise: RVWVD table not available.");
}
}
else {
if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
throw std::runtime_error(
"Cannot initialise: when no explicit RVWVD table is given, \n"
"datum depth must be at the gas-oil-contact. "
"In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold.");
}
const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure();
const double TContact = 273.15 + 20; // standard temperature for now
rvwFunc_.push_back(std::make_shared<Miscibility::RvwSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact));
}
}
}
else {
for (size_t i = 0; i < rec.size(); ++i) {
rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
}
}
// EXTRACT the initial temperature
updateInitialTemperature_(eclipseState);
@ -1603,7 +1683,7 @@ calcPressSatRsRv(const RMap& reg,
}
const auto eqreg = EquilReg {
rec[r], this->rsFunc_[r], this->rvFunc_[r], this->saltVdTable_[r], this->regionPvtIdx_[r]
rec[r], this->rsFunc_[r], this->rvFunc_[r], this->rvwFunc_[r], this->saltVdTable_[r], this->regionPvtIdx_[r]
};
// Ensure gas/oil and oil/water contacts are within the span for the
@ -1664,9 +1744,10 @@ cellLoop(const CellRange& cells,
auto saturations = Details::PhaseQuantityValue{};
auto Rs = 0.0;
auto Rv = 0.0;
auto Rvw = 0.0;
for (const auto& cell : cells) {
eqmethod(cell, pressures, saturations, Rs, Rv);
eqmethod(cell, pressures, saturations, Rs, Rv, Rvw);
if (oilActive) {
this->pp_ [oilPos][cell] = pressures.oil;
@ -1687,6 +1768,10 @@ cellLoop(const CellRange& cells,
this->rs_[cell] = Rs;
this->rv_[cell] = Rv;
}
if (watActive && gasActive) {
this->rvw_[cell] = Rvw;
}
}
}
@ -1714,7 +1799,8 @@ equilibrateCellCentres(const CellRange& cells,
Details::PhaseQuantityValue& pressures,
Details::PhaseQuantityValue& saturations,
double& Rs,
double& Rv) -> void
double& Rv,
double& Rvw) -> void
{
const auto pos = CellPos {
cell, cellCenterDepth_[cell]
@ -1730,6 +1816,9 @@ equilibrateCellCentres(const CellRange& cells,
Rv = eqreg.evaporationCalculator()
(pos.depth, pressures.gas, temp, saturations.oil);
Rvw = eqreg.waterEvaporationCalculator()
(pos.depth, pressures.gas, temp, saturations.water);
});
}
@ -1759,7 +1848,8 @@ equilibrateHorizontal(const CellRange& cells,
Details::PhaseQuantityValue& pressures,
Details::PhaseQuantityValue& saturations,
double& Rs,
double& Rv) -> void
double& Rv,
double& Rvw) -> void
{
pressures .reset();
saturations.reset();
@ -1785,6 +1875,9 @@ equilibrateHorizontal(const CellRange& cells,
Rv = eqreg.evaporationCalculator()
(cz, pressures.gas, temp, saturations.oil);
Rvw = eqreg.waterEvaporationCalculator()
(cz, pressures.gas, temp, saturations.water);
});
}

7
ebos/equil/initstateequil.hh
View File

@ -128,12 +128,13 @@ private:
const double press) const;
};
template <class FluidSystem, class RV>
template <class FluidSystem, class RV, class RVW>
class Gas
{
public:
Gas(const double temp,
const RV& rv,
const RVW& rvw,
const int pvtRegionIdx,
const double normGrav);
@ -143,6 +144,7 @@ public:
private:
const double temp_;
const RV& rv_;
const RVW& rvw_;
const int pvtRegionIdx_;
const double g_;
@ -276,7 +278,7 @@ private:
>;
using GasPressODE = PhasePressODE::Gas<
FluidSystem, typename Region::CalcEvaporation
FluidSystem, typename Region::CalcEvaporation, typename Region::CalcWaterEvaporation
>;
using WatPressODE = PhasePressODE::Water<FluidSystem>;
@ -741,6 +743,7 @@ private:
std::vector< std::shared_ptr<Miscibility::RsFunction> > rsFunc_;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rvFunc_;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rvwFunc_;
using TabulatedFunction = Tabulated1DFunction<double>;
std::vector<TabulatedFunction> saltVdTable_;
std::vector<TabulatedFunction> saltpVdTable_;

11
tests/test_equil.cc
View File

@ -200,7 +200,7 @@ static Opm::EquilRecord mkEquilRecord( double datd, double datp,
double zwoc, double pcow_woc,
double zgoc, double pcgo_goc )
{
return Opm::EquilRecord( datd, datp, zwoc, pcow_woc, zgoc, pcgo_goc, true, true, 0);
return Opm::EquilRecord( datd, datp, zwoc, pcow_woc, zgoc, pcgo_goc, true, true, 0, true);
}
template <typename Simulator>
@ -264,6 +264,7 @@ BOOST_AUTO_TEST_CASE(PhasePressure)
record,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0
};
@ -320,24 +321,28 @@ BOOST_AUTO_TEST_CASE(CellSubset)
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
};
@ -423,24 +428,28 @@ BOOST_AUTO_TEST_CASE(RegMapping)
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
trivialSaltVdTable,
0)
};