diff --git a/ebos/equil/equilibrationhelpers.cc b/ebos/equil/equilibrationhelpers.cc index cf01c3005..0a4325efa 100644 --- a/ebos/equil/equilibrationhelpers.cc +++ b/ebos/equil/equilibrationhelpers.cc @@ -195,6 +195,42 @@ satRv(const double press, const double temp) const } +template +RvwVD::RvwVD(const int pvtRegionIdx, + const std::vector& depth, + const std::vector& rvw) + : pvtRegionIdx_(pvtRegionIdx) + , rvwVsDepth_(depth, rvw) +{ +} + +template +double RvwVD:: +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 +double RvwVD:: +satRvw(const double press, const double temp) const +{ + return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press); +} + + template RsSatAtContact:: 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 +RvwSatAtContact:: +RvwSatAtContact(const int pvtRegionIdx, const double pContact, const double T_contact) + : pvtRegionIdx_(pvtRegionIdx) +{ + rvwSatContact_ = satRvw(pContact, T_contact); +} + +template +double RvwSatAtContact:: +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 +double RvwSatAtContact:: +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 rs, std::shared_ptr rv, + std::shared_ptr 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; template class RsSatAtContact; template class RvSatAtContact; + template class RvwSatAtContact; template class RvVD; + template class RvwVD; } } // namespace Equil diff --git a/ebos/equil/equilibrationhelpers.hh b/ebos/equil/equilibrationhelpers.hh index c22098ef8..69dde260f 100644 --- a/ebos/equil/equilibrationhelpers.hh +++ b/ebos/equil/equilibrationhelpers.hh @@ -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 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& depth, + const std::vector& 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; + + 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 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 rs, std::shared_ptr rv, + std::shared_ptr 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 rs_; /**< RS calculator */ std::shared_ptr rv_; /**< RV calculator */ + std::shared_ptr rvw_; /**< RVW calculator */ const TabulatedFunction& saltVdTable_; const int pvtIdx_; }; diff --git a/ebos/equil/initstateequil.cc b/ebos/equil/initstateequil.cc index f2b880a76..6d2701f26 100644 --- a/ebos/equil/initstateequil.cc +++ b/ebos/equil/initstateequil.cc @@ -34,6 +34,7 @@ #include #include #include +#include #include #include #include @@ -324,29 +325,31 @@ density(const double depth, return rho; } -template -Gas:: +template +Gas:: 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 -double Gas:: +template +double Gas:: operator()(const double depth, const double press) const { return this->density(depth, press) * g_; } -template -double Gas:: +template +double Gas:: 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(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>()); + continue; + } + const int pvtIdx = regionPvtIdx_[i]; + if (!rec[i].humidGasInitConstantRvw()) { + const TableContainer& rvwvdTables = tables.getRvwvdTables(); + + if (rvwvdTables.size() > 0) { + const RvwvdTable& rvwvdTable = rvwvdTables.getTable(i); + std::vector depthColumn = rvwvdTable.getColumn("DEPTH").vectorCopy(); + std::vector rvwvdColumn = rvwvdTable.getColumn("RVWVD").vectorCopy(); + rvwFunc_.push_back(std::make_shared>(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>(pvtIdx,pContact, TContact)); + } + } + } + else { + for (size_t i = 0; i < rec.size(); ++i) { + rvwFunc_.push_back(std::make_shared()); + } + } + + // 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); }); } diff --git a/ebos/equil/initstateequil.hh b/ebos/equil/initstateequil.hh index da574820a..af995160c 100644 --- a/ebos/equil/initstateequil.hh +++ b/ebos/equil/initstateequil.hh @@ -128,12 +128,13 @@ private: const double press) const; }; -template +template 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; @@ -741,6 +743,7 @@ private: std::vector< std::shared_ptr > rsFunc_; std::vector< std::shared_ptr > rvFunc_; + std::vector< std::shared_ptr > rvwFunc_; using TabulatedFunction = Tabulated1DFunction; std::vector saltVdTable_; std::vector saltpVdTable_; diff --git a/tests/test_equil.cc b/tests/test_equil.cc index 16213a2ac..7222d2a9e 100644 --- a/tests/test_equil.cc +++ b/tests/test_equil.cc @@ -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 @@ -264,6 +264,7 @@ BOOST_AUTO_TEST_CASE(PhasePressure) record, std::make_shared(), std::make_shared(), + std::make_shared(), trivialSaltVdTable, 0 }; @@ -320,24 +321,28 @@ BOOST_AUTO_TEST_CASE(CellSubset) Opm::EQUIL::EquilReg(record[0], std::make_shared(), std::make_shared(), + std::make_shared(), trivialSaltVdTable, 0) , Opm::EQUIL::EquilReg(record[0], std::make_shared(), std::make_shared(), + std::make_shared(), trivialSaltVdTable, 0) , Opm::EQUIL::EquilReg(record[1], std::make_shared(), std::make_shared(), + std::make_shared(), trivialSaltVdTable, 0) , Opm::EQUIL::EquilReg(record[1], std::make_shared(), std::make_shared(), + std::make_shared(), trivialSaltVdTable, 0) }; @@ -423,24 +428,28 @@ BOOST_AUTO_TEST_CASE(RegMapping) Opm::EQUIL::EquilReg(record[0], std::make_shared(), std::make_shared(), + std::make_shared(), trivialSaltVdTable, 0) , Opm::EQUIL::EquilReg(record[0], std::make_shared(), std::make_shared(), + std::make_shared(), trivialSaltVdTable, 0) , Opm::EQUIL::EquilReg(record[1], std::make_shared(), std::make_shared(), + std::make_shared(), trivialSaltVdTable, 0) , Opm::EQUIL::EquilReg(record[1], std::make_shared(), std::make_shared(), + std::make_shared(), trivialSaltVdTable, 0) };