diff --git a/opm/simulators/flow/BlackoilModel.hpp b/opm/simulators/flow/BlackoilModel.hpp index 7bdaee297..1391e4359 100644 --- a/opm/simulators/flow/BlackoilModel.hpp +++ b/opm/simulators/flow/BlackoilModel.hpp @@ -130,6 +130,10 @@ template struct EnableDispersion { static constexpr bool value = false; }; +template +struct EnableConvectiveMixing +{ static constexpr bool value = true; }; + template struct WellModel { using type = BlackoilWellModel; }; diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp index f444b65ef..3116c182f 100644 --- a/opm/simulators/flow/FlowProblem.hpp +++ b/opm/simulators/flow/FlowProblem.hpp @@ -57,7 +57,9 @@ #include #include #include +#include #include +#include #include @@ -143,6 +145,7 @@ class FlowProblem : public GetPropType enum { enableEnergy = getPropValue() }; enum { enableDiffusion = getPropValue() }; enum { enableDispersion = getPropValue() }; + enum { enableConvectiveMixing = getPropValue() }; enum { enableThermalFluxBoundaries = getPropValue() }; enum { enableApiTracking = getPropValue() }; enum { enableMICP = getPropValue() }; @@ -180,6 +183,8 @@ class FlowProblem : public GetPropType using MICPModule = BlackOilMICPModule; using DispersionModule = BlackOilDispersionModule; using DiffusionModule = BlackOilDiffusionModule; + using ConvectiveMixingModule = BlackOilConvectiveMixingModule; + using ModuleParams = typename BlackOilLocalResidualTPFA::ModuleParams; using InitialFluidState = typename EquilInitializer::ScalarFluidState; @@ -589,6 +594,8 @@ public: FluidSystem::setVapPars(0.0, 0.0); } } + + ConvectiveMixingModule::beginEpisode(simulator.vanguard().eclState(), simulator.vanguard().schedule(), episodeIdx, moduleParams_.convectiveMixingModuleParam); } /*! @@ -1855,6 +1862,10 @@ public: eclWriter_->mutableOutputModule().setCnvData(data); } + const ModuleParams& moduleParams() const { + return moduleParams_; + } + template void serializeOp(Serializer& serializer) { @@ -1953,6 +1964,7 @@ protected: return true; } + bool updateMaxOilSaturation_() { OPM_TIMEBLOCK(updateMaxOilSaturation); @@ -2597,7 +2609,6 @@ protected: return this->rockCompTransMultVal_[dofIdx]; } - private: struct PffDofData_ { @@ -2883,6 +2894,10 @@ private: BCData bcindex_; bool nonTrivialBoundaryConditions_ = false; bool explicitRockCompaction_ = false; + + ModuleParams moduleParams_; + + }; } // namespace Opm diff --git a/opm/simulators/flow/FlowProblemProperties.hpp b/opm/simulators/flow/FlowProblemProperties.hpp index 803e1e404..94b998530 100644 --- a/opm/simulators/flow/FlowProblemProperties.hpp +++ b/opm/simulators/flow/FlowProblemProperties.hpp @@ -198,6 +198,11 @@ template struct EnableDispersion { static constexpr bool value = false; }; +// Enable Convective Mixing +template +struct EnableConvectiveMixing +{ static constexpr bool value = true; }; + // disable API tracking template struct EnableApiTracking diff --git a/opm/simulators/flow/MixingRateControls.cpp b/opm/simulators/flow/MixingRateControls.cpp index 9c4aecca3..577c90935 100644 --- a/opm/simulators/flow/MixingRateControls.cpp +++ b/opm/simulators/flow/MixingRateControls.cpp @@ -101,19 +101,16 @@ init(std::size_t numDof, int episodeIdx, const unsigned ntpvt) //TODO We may want to only allocate these properties only if active. //But since they may be activated at later time we need some more //intrastructure to handle it - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && - FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + maxDRv_.resize(ntpvt, 1e30); lastRv_.resize(numDof, 0.0); maxDRs_.resize(ntpvt, 1e30); dRsDtOnlyFreeGas_.resize(ntpvt, false); lastRs_.resize(numDof, 0.0); maxDRv_.resize(ntpvt, 1e30); - lastRv_.resize(numDof, 0.0); if (this->drsdtConvective(episodeIdx)) { convectiveDrs_.resize(numDof, 1.0); } - } } template @@ -121,9 +118,7 @@ bool MixingRateControls:: drsdtActive(int episodeIdx) const { const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); - const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && - FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - return (oilVaporizationControl.drsdtActive() && bothOilGasActive); + return (oilVaporizationControl.drsdtActive()); } template @@ -131,9 +126,7 @@ bool MixingRateControls:: drvdtActive(int episodeIdx) const { const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); - const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && - FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - return (oilVaporizationControl.drvdtActive() && bothOilGasActive); + return (oilVaporizationControl.drvdtActive()); } template @@ -141,9 +134,7 @@ bool MixingRateControls:: drsdtConvective(int episodeIdx) const { const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); - const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && - FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - return (oilVaporizationControl.drsdtConvective() && bothOilGasActive); + return (oilVaporizationControl.drsdtConvective()); } template @@ -271,31 +262,64 @@ void MixingRateControls:: updateConvectiveDRsDt_(const unsigned compressedDofIdx, const Scalar t, const Scalar p, + const Scalar pg, const Scalar rs, - const Scalar so, + const Scalar sg, const Scalar poro, const Scalar permz, const Scalar distZ, const Scalar gravity, + const Scalar salt, + const Scalar Xhi, + const Scalar Psi, + const Scalar omegainn, const int pvtRegionIndex) { - const Scalar rssat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIndex, t, p); - const Scalar saturatedInvB - = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIndex, t, p); + const Scalar rssat = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ? + FluidSystem::waterPvt().saturatedGasDissolutionFactor(pvtRegionIndex, t, p, salt) : + FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIndex, t, p); + const Scalar saturatedInvB = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ? + FluidSystem::waterPvt().saturatedInverseFormationVolumeFactor(pvtRegionIndex, t, p, salt) : + FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIndex, t, p); const Scalar rsZero = 0.0; - const Scalar pureDensity - = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIndex, t, p, rsZero) - * FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex); - const Scalar saturatedDensity = saturatedInvB - * (FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex) - + rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIndex)); - const Scalar deltaDensity = saturatedDensity - pureDensity; - const Scalar visc = FluidSystem::oilPvt().viscosity(pvtRegionIndex, t, p, rs); - // Note that for so = 0 this gives no limits (inf) for the dissolution rate + const Scalar sg_max = 1.0; + const Scalar pureDensity = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ? + (FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIndex, t, p, rsZero, salt) + * FluidSystem::waterPvt().waterReferenceDensity(pvtRegionIndex)) : + (FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIndex, t, p, rsZero) + * FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex)); + const Scalar saturatedDensity = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ? + (saturatedInvB * (FluidSystem::waterPvt().waterReferenceDensity(pvtRegionIndex) + + rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIndex))) : + (saturatedInvB * (FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex) + + rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIndex))); + Scalar deltaDensity = saturatedDensity - pureDensity; + const Scalar visc = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ? + FluidSystem::waterPvt().viscosity(pvtRegionIndex, t, p, rs, salt) : + FluidSystem::oilPvt().viscosity(pvtRegionIndex, t, p, rs); + + // Note that for sLiquid = 0 this gives no limits (inf) for the dissolution rate // Also we restrict the effect of convective mixing to positive density differences // i.e. we only allow for fingers moving downward + + Scalar co2Density = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIndex, + t,p,0.0 /*=Rv*/, 0.0 /*=Rvw*/) * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIndex); + Scalar factor = 1.0; + Scalar X = (rs - rssat * sg) / (rssat * ( 1.0 - sg)); + Scalar omega = 0.0; + const Scalar pCap = Opm::abs(pg - p); + if ((rs >= (rssat * sg)) || (pCap < 1e-12)){ + if(X > Psi){ + factor = 0.0; + omega = omegainn; + } + } else { + factor /= Xhi; + deltaDensity = (saturatedDensity - co2Density); + } + convectiveDrs_[compressedDofIdx] - = permz * rssat * max(0.0, deltaDensity) * gravity / (so * visc * distZ * poro); + = factor * permz * rssat * max(0.0, deltaDensity) * gravity / ( std::max(sg_max - sg, 0.0) * visc * distZ * poro) + (omega/Xhi); } template class MixingRateControls>; diff --git a/opm/simulators/flow/MixingRateControls.hpp b/opm/simulators/flow/MixingRateControls.hpp index 401233ed6..f5f5dfe90 100644 --- a/opm/simulators/flow/MixingRateControls.hpp +++ b/opm/simulators/flow/MixingRateControls.hpp @@ -117,20 +117,41 @@ public: const int pvtRegionIdx, const std::array& active) { + const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); if (active[0]) { // This implements the convective DRSDT as described in // Sandve et al. "Convective dissolution in field scale CO2 storage simulations using the OPM Flow // simulator" Submitted to TCCS 11, 2021 + // modification and introduction of regimes following Mykkeltvedt et al. Submitted to TIMP 2024 const auto& fs = iq.fluidState(); + + const auto& temperature = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ? + getValue(fs.temperature(FluidSystem::waterPhaseIdx)) : + getValue(fs.temperature(FluidSystem::oilPhaseIdx)); + const auto& pressure = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ? + getValue(fs.pressure(FluidSystem::waterPhaseIdx)) : + getValue(fs.pressure(FluidSystem::oilPhaseIdx)); + const auto& pressuregas = getValue(fs.pressure(FluidSystem::gasPhaseIdx)); + const auto& rs = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ? + getValue(fs.Rsw()) : + getValue(fs.Rs()); + + const auto& salt = getValue(fs.saltSaturation()); + this->updateConvectiveDRsDt_(compressedDofIdx, - getValue(fs.temperature(FluidSystem::oilPhaseIdx)), - getValue(fs.pressure(FluidSystem::oilPhaseIdx)), - getValue(fs.Rs()), - getValue(fs.saturation(FluidSystem::oilPhaseIdx)), + temperature, + pressure, + pressuregas, + rs, + getValue(fs.saturation(FluidSystem::gasPhaseIdx)), getValue(iq.porosity()), permZ, distZ, gravity, + salt, + oilVaporizationControl.getMaxDRSDT(fs.pvtRegionIndex()), + oilVaporizationControl.getPsi(fs.pvtRegionIndex()), + oilVaporizationControl.getOmega(fs.pvtRegionIndex()), fs.pvtRegionIndex()); } @@ -138,13 +159,13 @@ public: const auto& fs = iq.fluidState(); using FluidState = typename std::decay::type; - - const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap(); constexpr Scalar freeGasMinSaturation_ = 1e-7; if (oilVaporizationControl.getOption(pvtRegionIdx) || fs.saturation(FluidSystem::gasPhaseIdx) > freeGasMinSaturation_) { lastRs_[compressedDofIdx] - = BlackOil::template getRs_(fs, iq.pvtRegionIndex()); + = ((FluidSystem::enableDissolvedGasInWater())) ? + BlackOil::template getRsw_(fs, iq.pvtRegionIndex()) : + BlackOil::template getRs_(fs, iq.pvtRegionIndex()); } else lastRs_[compressedDofIdx] = std::numeric_limits::infinity(); @@ -162,12 +183,17 @@ private: void updateConvectiveDRsDt_(const unsigned compressedDofIdx, const Scalar t, const Scalar p, + const Scalar pg, const Scalar rs, - const Scalar so, + const Scalar sg, const Scalar poro, const Scalar permz, const Scalar distZ, const Scalar gravity, + const Scalar salt, + const Scalar Xhi, + const Scalar Psi, + const Scalar omegainn, const int pvtRegionIndex); std::vector lastRv_; diff --git a/opm/simulators/flow/NewTranFluxModule.hpp b/opm/simulators/flow/NewTranFluxModule.hpp index 7b541d001..39f74c58a 100644 --- a/opm/simulators/flow/NewTranFluxModule.hpp +++ b/opm/simulators/flow/NewTranFluxModule.hpp @@ -44,6 +44,7 @@ #include #include #include +#include #include @@ -122,11 +123,16 @@ class NewTranExtensiveQuantities enum { enableExtbo = getPropValue() }; enum { enableEnergy = getPropValue() }; + static constexpr bool enableConvectiveMixing = getPropValue(); + + using Toolbox = MathToolbox; using DimVector = Dune::FieldVector; using EvalDimVector = Dune::FieldVector; using DimMatrix = Dune::FieldMatrix; + using ConvectiveMixingModule = BlackOilConvectiveMixingModule; + using ModuleParams = typename BlackOilLocalResidualTPFA::ModuleParams; public: /*! * \brief Return the intrinsic permeability tensor at a face [m^2] @@ -277,7 +283,8 @@ public: I, J, distZ*g, - thpres); + thpres, + problem.moduleParams()); if (pressureDifferences[phaseIdx] == 0) { volumeFlux[phaseIdx] = 0.0; continue; @@ -317,8 +324,8 @@ public: const unsigned globalIndexIn, const unsigned globalIndexEx, const Scalar distZg, - const Scalar thpres - ) + const Scalar thpres, + const ModuleParams& moduleParams) { // check shortcut: if the mobility of the phase is zero in the interior as @@ -338,6 +345,10 @@ public: Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx)); Evaluation rhoAvg = (rhoIn + rhoEx)/2; + if constexpr(enableConvectiveMixing) { + ConvectiveMixingModule::modifyAvgDensity(rhoAvg, intQuantsIn, intQuantsEx, phaseIdx, moduleParams.convectiveMixingModuleParam); + } + const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx); Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(phaseIdx)); if (enableExtbo) // added stability; particulary useful for solvent migrating in pure water