Merge pull request #5491 from totto82/convmix

Convmix

opm/simulators/flow/BlackoilModel.hpp

@@ -130,6 +130,10 @@ template<class TypeTag>
 struct EnableDispersion<TypeTag, TTag::FlowProblem>
 { static constexpr bool value = false; };
 
+template<class TypeTag>
+struct EnableConvectiveMixing<TypeTag, TTag::FlowProblem>
+{ static constexpr bool value = true; };
+
 template<class TypeTag>
 struct WellModel<TypeTag, TTag::FlowProblem>
 { using type = BlackoilWellModel<TypeTag>; };

opm/simulators/flow/FlowProblem.hpp

@@ -57,7 +57,9 @@
 #include <opm/models/common/directionalmobility.hh>
 #include <opm/models/utils/pffgridvector.hh>
 #include <opm/models/blackoil/blackoilmodel.hh>
+#include <opm/models/blackoil/blackoilconvectivemixingmodule.hh>
 #include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
+#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
 
 #include <opm/output/eclipse/EclipseIO.hpp>
 
@@ -143,6 +145,7 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
     enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
+    enum { enableConvectiveMixing = getPropValue<TypeTag, Properties::EnableConvectiveMixing>() };
     enum { enableThermalFluxBoundaries = getPropValue<TypeTag, Properties::EnableThermalFluxBoundaries>() };
     enum { enableApiTracking = getPropValue<TypeTag, Properties::EnableApiTracking>() };
     enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
@@ -180,6 +183,8 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
     using MICPModule = BlackOilMICPModule<TypeTag>;
     using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
     using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
+    using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag, enableConvectiveMixing>;
+    using ModuleParams = typename BlackOilLocalResidualTPFA<TypeTag>::ModuleParams;
 
     using InitialFluidState = typename EquilInitializer<TypeTag>::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<class Serializer>
     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<int> bcindex_;
     bool nonTrivialBoundaryConditions_ = false;
     bool explicitRockCompaction_ = false;
+
+    ModuleParams moduleParams_;
+
+
 };
 
 } // namespace Opm

opm/simulators/flow/FlowProblemProperties.hpp

@@ -198,6 +198,11 @@ template<class TypeTag>
 struct EnableDispersion<TypeTag, TTag::FlowBaseProblem>
 { static constexpr bool value = false; };
 
+// Enable Convective Mixing
+template<class TypeTag>
+struct EnableConvectiveMixing<TypeTag, TTag::FlowBaseProblem>
+{ static constexpr bool value = true; };
+
 // disable API tracking
 template<class TypeTag>
 struct EnableApiTracking<TypeTag, TTag::FlowBaseProblem>

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<class FluidSystem>
@@ -121,9 +118,7 @@ bool MixingRateControls<FluidSystem>::
 drsdtActive(int episodeIdx) const
 {
     const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
-    const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
+    return (oilVaporizationControl.drsdtActive());
-                                  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
-    return (oilVaporizationControl.drsdtActive() && bothOilGasActive);
 }
 
 template<class FluidSystem>
@@ -131,9 +126,7 @@ bool MixingRateControls<FluidSystem>::
 drvdtActive(int episodeIdx) const
 {
     const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
-    const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
+    return (oilVaporizationControl.drvdtActive());
-                                  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
-    return (oilVaporizationControl.drvdtActive() && bothOilGasActive);
 }
 
 template<class FluidSystem>
@@ -141,9 +134,7 @@ bool MixingRateControls<FluidSystem>::
 drsdtConvective(int episodeIdx) const
 {
     const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
-    const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
+    return (oilVaporizationControl.drsdtConvective());
-                                  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
-    return (oilVaporizationControl.drsdtConvective() && bothOilGasActive);
 }
 
 template<class FluidSystem>
@@ -271,31 +262,64 @@ void MixingRateControls<FluidSystem>::
 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 rssat = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
-    const Scalar saturatedInvB
+        FluidSystem::waterPvt().saturatedGasDissolutionFactor(pvtRegionIndex, t, p, salt) :
-        = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIndex, t, p);
+        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
+    const Scalar sg_max = 1.0;
-        = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIndex, t, p, rsZero)
+    const Scalar pureDensity = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
-        * FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex);
+        (FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIndex, t, p, rsZero, salt)
-    const Scalar saturatedDensity = saturatedInvB
+        * FluidSystem::waterPvt().waterReferenceDensity(pvtRegionIndex)) :
-        * (FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex)
+         (FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIndex, t, p, rsZero)
-           + rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIndex));
+        * FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex));
-    const Scalar deltaDensity = saturatedDensity - pureDensity;
+    const Scalar saturatedDensity = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
-    const Scalar visc = FluidSystem::oilPvt().viscosity(pvtRegionIndex, t, p, rs);
+        (saturatedInvB * (FluidSystem::waterPvt().waterReferenceDensity(pvtRegionIndex)
-    // Note that for so = 0 this gives no limits (inf) for the dissolution rate
+        + 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<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>>;

opm/simulators/flow/MixingRateControls.hpp

@@ -117,20 +117,41 @@ public:
                 const int pvtRegionIdx,
                 const std::array<bool,3>& 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)),
+                                         temperature,
-                                         getValue(fs.pressure(FluidSystem::oilPhaseIdx)),
+                                         pressure,
-                                         getValue(fs.Rs()),
+                                         pressuregas,
-                                         getValue(fs.saturation(FluidSystem::oilPhaseIdx)),
+                                         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<decltype(fs)>::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_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex());
+                    = ((FluidSystem::enableDissolvedGasInWater())) ?
+                    BlackOil::template getRsw_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex()) :
+                    BlackOil::template getRs_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex());
             }
             else
                 lastRs_[compressedDofIdx] = std::numeric_limits<Scalar>::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<Scalar> lastRv_;

opm/simulators/flow/NewTranFluxModule.hpp

@@ -44,6 +44,7 @@
 #include <opm/models/discretization/common/fvbaseproperties.hh>
 #include <opm/models/blackoil/blackoilproperties.hh>
 #include <opm/models/utils/signum.hh>
+#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
 
 #include <array>
 
@@ -122,11 +123,16 @@ class NewTranExtensiveQuantities
     enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
 
+    static constexpr bool enableConvectiveMixing = getPropValue<TypeTag, Properties::EnableConvectiveMixing>();
+
+
     using Toolbox = MathToolbox<Evaluation>;
     using DimVector = Dune::FieldVector<Scalar, dimWorld>;
     using EvalDimVector = Dune::FieldVector<Evaluation, dimWorld>;
     using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
 
+    using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag, enableConvectiveMixing>;
+    using ModuleParams = typename BlackOilLocalResidualTPFA<TypeTag>::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