Merge pull request #684 from plgbrts/vapwat

Enabling water evaporation into the gas phase

opm/models/blackoil/blackoilbrinemodules.hh

@@ -266,12 +266,7 @@ public:
         if (!enableBrine)
             return;
 
-        static unsigned contiGasEqIdx, contiOilEqIdx;
-        if(gasEnabled) { contiGasEqIdx = Indices::conti0EqIdx + Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); }
-        if(oilEnabled) { contiOilEqIdx = Indices::conti0EqIdx + Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); }
-
         const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
-        const auto& intQuants = elemCtx.intensiveQuantities(scvfIdx, timeIdx);
 
         const unsigned upIdx = extQuants.upstreamIndex(FluidSystem::waterPhaseIdx);
         const unsigned inIdx = extQuants.interiorIndex();
@@ -282,24 +277,12 @@ public:
                     extQuants.volumeFlux(waterPhaseIdx)
                     *up.fluidState().invB(waterPhaseIdx)
                     *up.fluidState().saltConcentration();
-
-            if (enableSaltPrecipitation) {
-                // modify gas and oil flux for mobility change
-                if(gasEnabled) { flux[contiGasEqIdx] *= intQuants.permFactor(); }
-                if(oilEnabled) { flux[contiOilEqIdx] *= intQuants.permFactor(); }
-            }
         }
         else {
             flux[contiBrineEqIdx] =
                     extQuants.volumeFlux(waterPhaseIdx)
                     *decay<Scalar>(up.fluidState().invB(waterPhaseIdx))
                     *decay<Scalar>(up.fluidState().saltConcentration());
-
-            if (enableSaltPrecipitation) {
-                // modify gas and oil flux for mobility change
-                if(gasEnabled) { flux[contiGasEqIdx] *= decay<Scalar>(intQuants.permFactor()); }
-                if(oilEnabled) { flux[contiOilEqIdx] *= decay<Scalar>(intQuants.permFactor()); }
-            }
         }
     }
 
@@ -489,11 +472,13 @@ public:
             const auto& permfactTable = BrineModule::permfactTable(elemCtx, dofIdx, timeIdx);
 
             permFactor_ = permfactTable.eval(scalarValue(porosityFactor));
-            if (permFactor_ < 1 ) {
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-                // adjust mobility for changing permeability
+                if (!FluidSystem::phaseIsActive(phaseIdx))
-                asImp_().mobility_[waterPhaseIdx] *= permFactor_;
+                    continue;
+
+                asImp_().mobility_[phaseIdx] *= permFactor_;
             }
-         }
+        }
     }
 
     const Evaluation& saltConcentration() const

opm/models/blackoil/blackoilintensivequantities.hh

@@ -149,7 +149,7 @@ public:
         if (waterEnabled) {
             if (priVars.primaryVarsMeaning() == PrimaryVariables::OnePhase_p) {
                 Sw = 1.0;
-            } else {
+            } else if (priVars.primaryVarsMeaning() != PrimaryVariables::Rvw_po_Sg) {
                 Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx);
             }
         }
@@ -236,7 +236,7 @@ public:
                         elemCtx.problem().maxOilSaturation(globalSpaceIdx));
         }
 
-        // take the meaning of the switiching primary variable into account for the gas
+        // take the meaning of the switching primary variable into account for the gas
         // and oil phase compositions
         if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
             // in the threephase case, gas and oil phases are potentially present, i.e.,
@@ -264,6 +264,18 @@ public:
             }
             else if (compositionSwitchEnabled)
                 fluidState_.setRv(0.0);
+
+            if (FluidSystem::enableVaporizedWater()) {
+                const Evaluation& RvwSat = FluidSystem::saturatedVaporizationFactor(fluidState_,
+                                                            gasPhaseIdx,
+                                                            pvtRegionIdx);
+                fluidState_.setRvw(RvwSat);
+            }
+        }
+        else if (priVars.primaryVarsMeaning() == PrimaryVariables::Rvw_po_Sg) {
+            // The switching variable is the water-gas ratio Rvw
+            const auto& Rvw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx);
+            fluidState_.setRvw(Rvw);
         }
         else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
             // if the switching variable is the mole fraction of the gas component in the

opm/models/blackoil/blackoillocalresidual.hh

@@ -142,6 +142,14 @@ public:
                     Toolbox::template decay<LhsEval>(intQuants.fluidState().Rv())
                     * surfaceVolume;
             }
+
+            // account for vaporized water
+            if (phaseIdx == gasPhaseIdx && FluidSystem::enableVaporizedWater()) {
+                unsigned activeWaterCompIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx);
+                storage[conti0EqIdx + activeWaterCompIdx] +=
+                    Toolbox::template decay<LhsEval>(intQuants.fluidState().Rvw())
+                    * surfaceVolume;
+            }
         }
 
         adaptMassConservationQuantities_(storage, intQuants.pvtRegionIndex());
@@ -281,6 +289,16 @@ public:
                 else
                     flux[conti0EqIdx + activeOilCompIdx] += Rv*surfaceVolumeFlux*FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
             }
+             // vaporized water (in the gas phase).
+            if (FluidSystem::enableVaporizedWater()) {
+                const auto& Rvw = BlackOil::getRvw_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
+
+                unsigned activeWaterCompIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx);
+                if (blackoilConserveSurfaceVolume)
+                    flux[conti0EqIdx + activeWaterCompIdx] += Rvw*surfaceVolumeFlux;
+                else
+                    flux[conti0EqIdx + activeWaterCompIdx] += Rvw*surfaceVolumeFlux*FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
+            }
         }
     }
 

opm/models/blackoil/blackoilnewtonmethod.hh

@@ -264,7 +264,7 @@ protected:
         Scalar deltaSg = 0.0;
         Scalar deltaSs = 0.0;
 
-        if (Indices::waterEnabled && FluidSystem::numActivePhases() > 1) {
+        if (Indices::waterEnabled && FluidSystem::numActivePhases() > 1 && currentValue.primaryVarsMeaning() != PrimaryVariables::Rvw_po_Sg) {
             deltaSw = update[Indices::waterSaturationIdx];
             deltaSo = -deltaSw;
         }
@@ -306,7 +306,13 @@ protected:
             }
             // water saturation delta
             else if (pvIdx == Indices::waterSaturationIdx)
-                delta *= satAlpha;
+                if (currentValue.primaryVarsMeaning() != PrimaryVariables::Rvw_po_Sg)
+                    delta *= satAlpha;
+                else{
+                    //Ensure Rvw factor does not become negative
+                    if (delta > currentValue[ Indices::waterSaturationIdx]) 
+                        delta = currentValue[ Indices::waterSaturationIdx];
+                }
             else if (pvIdx == Indices::compositionSwitchIdx) {
                 // the switching primary variable for composition is tricky because the
                 // "reasonable" value ranges it exhibits vary widely depending on its
@@ -342,6 +348,11 @@ protected:
                 const double sign = delta >= 0. ? 1. : -1.;
                 delta = sign * std::min(std::abs(delta), maxTempChange_);
             }
+            else if (enableBrine && pvIdx == Indices::saltConcentrationIdx && currentValue.primaryVarsMeaningBrine() == PrimaryVariables::Sp) { 
+                const double maxSaltSaturationChange = 0.1;
+                const double sign = delta >= 0. ? 1. : -1.;
+                delta = sign * std::min(std::abs(delta), maxSaltSaturationChange);
+            }
 
             // do the actual update
             nextValue[pvIdx] = currentValue[pvIdx] - delta;
@@ -369,10 +380,13 @@ protected:
             if (enableFoam && pvIdx == Indices::foamConcentrationIdx)
                 nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);
 
-            // keep the salt concentration above 0
             if (enableBrine && pvIdx == Indices::saltConcentrationIdx) { 
+               // keep the salt concentration above 0
                if (!enableSaltPrecipitation || (enableSaltPrecipitation && currentValue.primaryVarsMeaningBrine() == PrimaryVariables::Cs))
                    nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0); 
+               // keep the salt saturation below upperlimit
+               if ((enableSaltPrecipitation && currentValue.primaryVarsMeaningBrine() == PrimaryVariables::Sp))
+                   nextValue[pvIdx] = std::min(nextValue[pvIdx], 0.9); 
             }
 
             // keep the temperature within given values

opm/models/blackoil/blackoilprimaryvariables.hh

@@ -103,6 +103,7 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
     enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
     enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
     enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
+    enum { enableEvaporation = getPropValue<TypeTag, Properties::EnableEvaporation>() };
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
     enum { gasCompIdx = FluidSystem::gasCompIdx };
@@ -127,6 +128,8 @@ public:
         Sw_po_Sg, // threephase case
         Sw_po_Rs, // water + oil case
         Sw_pg_Rv, // water + gas case
+        Rvw_po_Sg, // gas + oil  case
+
         OnePhase_p, // onephase case
     };
 
@@ -289,6 +292,7 @@ public:
 
         bool gasPresent = FluidSystem::phaseIsActive(gasPhaseIdx)?(fluidState.saturation(gasPhaseIdx) > 0.0):false;
         bool oilPresent = FluidSystem::phaseIsActive(oilPhaseIdx)?(fluidState.saturation(oilPhaseIdx) > 0.0):false;
+        bool waterPresent = FluidSystem::phaseIsActive(waterPhaseIdx)?(fluidState.saturation(waterPhaseIdx) > 0.0):false;
         static const Scalar thresholdWaterFilledCell = 1.0 - 1e-6;
         bool onlyWater = FluidSystem::phaseIsActive(waterPhaseIdx)?(fluidState.saturation(waterPhaseIdx) > thresholdWaterFilledCell):false;
         const auto& saltSaturation = BlackOil::getSaltSaturation_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
@@ -320,6 +324,8 @@ public:
             // composition, if it is disabled, the oil component must stick to its phase
             if (FluidSystem::enableVaporizedOil())
                 primaryVarsMeaning_ = Sw_pg_Rv;
+            else if (FluidSystem::enableVaporizedWater() && !waterPresent)
+                primaryVarsMeaning_ = Rvw_po_Sg; 
             else
                 primaryVarsMeaning_ = Sw_po_Sg;
         }
@@ -361,6 +367,20 @@ public:
             if( compositionSwitchEnabled )
                 (*this)[compositionSwitchIdx] = Rs;
         }
+        else if (primaryVarsMeaning() == Rvw_po_Sg && FluidSystem::enableVaporizedWater()) {
+            const auto& Rvw = BlackOil::getRvw_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
+            if (waterEnabled)
+                (*this)[waterSaturationIdx] = Rvw; //waterSaturationIdx becomes a switching idx
+            if (gasEnabled && waterEnabled && !oilEnabled) {
+                //-> water-gas system
+                (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(gasPhaseIdx));
+            }
+            else if (oilEnabled) {
+                (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(oilPhaseIdx));
+            }
+            if( compositionSwitchEnabled )
+                (*this)[compositionSwitchIdx] = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
+        }
         else {
             assert(primaryVarsMeaning() == Sw_pg_Rv);
 
@@ -412,7 +432,7 @@ public:
             return false;
         }
         Scalar Sw = 0.0;
-        if (waterEnabled)
+        if (waterEnabled && primaryVarsMeaning() != Rvw_po_Sg)
             Sw = (*this)[Indices::waterSaturationIdx];
 
         if (enableSaltPrecipitation) {
@@ -455,6 +475,23 @@ public:
 
             Scalar So = 1.0 - Sw - Sg - solventSaturation_();
 
+            //water disappears 
+            if( Sw < -eps && FluidSystem::enableVaporizedWater()) {
+                Scalar po = (*this)[Indices::pressureSwitchIdx]; 
+                Scalar T = asImp_().temperature_();
+                Scalar pC[numPhases] = { 0.0 };
+                const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
+                computeCapillaryPressures_(pC, So, Sg + solventSaturation_(), Sw, matParams);
+                Scalar pg = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
+                Scalar RvwSat = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_,
+                                                                                      T,
+                                                                                      pg);
+                setPrimaryVarsMeaning(Rvw_po_Sg);
+                (*this)[Indices::waterSaturationIdx] = RvwSat; //primary variable becomes Rvw 
+
+                return true;
+            }
+
             Scalar So2 = 1.0 - Sw - solventSaturation_();
             if (Sg < -eps && So2 > 0.0 && FluidSystem::enableDissolvedGas()) {
                 // the hydrocarbon gas phase disappeared and some oil phase is left,
@@ -566,6 +603,32 @@ public:
 
             return false;
         }
+        else if (primaryVarsMeaning() == Rvw_po_Sg) {
+            
+            Scalar po = (*this)[Indices::pressureSwitchIdx]; 
+            Scalar T = asImp_().temperature_();
+            Scalar Sg = 0.0;
+            if (compositionSwitchEnabled)
+                Sg = (*this)[Indices::compositionSwitchIdx];
+
+            Scalar So = 1.0 - Sg - solventSaturation_();
+            Scalar pC[numPhases] = { 0.0 };
+            const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
+            computeCapillaryPressures_(pC, So, Sg + solventSaturation_(), /*Sw=*/ 0.0, matParams);
+            Scalar pg = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
+            Scalar RvwSat = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_,
+                                                                                   T,
+                                                                                   pg); 
+            Scalar Rvw = (*this)[Indices::waterSaturationIdx];
+            if (Rvw >  RvwSat*(1.0 + eps)) {
+                // water phase appears
+                setPrimaryVarsMeaning(Sw_po_Sg);
+                (*this)[Indices::waterSaturationIdx] = 0.0; // water saturation
+
+                return true;
+            }
+            return false;
+        }
         else {
             assert(primaryVarsMeaning() == Sw_pg_Rv);
             assert(compositionSwitchEnabled);
@@ -646,7 +709,7 @@ public:
             return false;
         }
         Scalar Sw = 0.0;
-        if (waterEnabled)
+        if (waterEnabled && primaryVarsMeaning() != Rvw_po_Sg)
             Sw = (*this)[Indices::waterSaturationIdx];
 
         if (primaryVarsMeaning() == Sw_po_Sg) {