Extend foam model with solvent or water

opm/models/blackoil/blackoilfoammodules.hh

@@ -75,6 +75,7 @@ class BlackOilFoamModule
     static constexpr unsigned foamConcentrationIdx = Indices::foamConcentrationIdx;
     static constexpr unsigned contiFoamEqIdx = Indices::contiFoamEqIdx;
     static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
+    static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
 
     static constexpr unsigned enableFoam = enableFoamV;
     static constexpr bool enableVtkOutput = getPropValue<TypeTag, Properties::EnableVtkOutput>();
@@ -82,6 +83,8 @@ class BlackOilFoamModule
     static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
     static constexpr unsigned numPhases = FluidSystem::numPhases;
 
+    enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
+
 public:
 #if HAVE_ECL_INPUT
     /*!
@@ -104,11 +107,8 @@ public:
             return; // foam treatment is supposed to be disabled
         }
 
-        // Check that only implemented options are used.
+        params_.transport_phase_ = eclState.getInitConfig().getFoamConfig().getTransportPhase();
-        // We only support the default values of FOAMOPTS (GAS, TAB).
+
-        if (eclState.getInitConfig().getFoamConfig().getTransportPhase() != Phase::GAS) {
-            throw std::runtime_error("In FOAMOPTS, only GAS is allowed for the transport phase.");
-        }
         if (eclState.getInitConfig().getFoamConfig().getMobilityModel() != FoamConfig::MobilityModel::TAB) {
             throw std::runtime_error("In FOAMOPTS, only TAB is allowed for the gas mobility factor reduction model.");
         }
@@ -249,16 +249,27 @@ public:
         if constexpr (enableFoam) {
             const auto& fs = intQuants.fluidState();
 
-            LhsEval surfaceVolumeFreeGas =
+            LhsEval surfaceVolume = Toolbox::template decay<LhsEval>(intQuants.porosity());
-                Toolbox::template decay<LhsEval>(fs.saturation(gasPhaseIdx))
+            if (params_.transport_phase_ == Phase::WATER) {
-                * Toolbox::template decay<LhsEval>(fs.invB(gasPhaseIdx))
+                surfaceVolume *= (Toolbox::template decay<LhsEval>(fs.saturation(waterPhaseIdx))
-                * Toolbox::template decay<LhsEval>(intQuants.porosity());
+                * Toolbox::template decay<LhsEval>(fs.invB(waterPhaseIdx)));
+            } else if (params_.transport_phase_ == Phase::GAS) {
+                surfaceVolume *= (Toolbox::template decay<LhsEval>(fs.saturation(gasPhaseIdx))
+                * Toolbox::template decay<LhsEval>(fs.invB(gasPhaseIdx)));
+            } else if (params_.transport_phase_ == Phase::SOLVENT) {
+                if constexpr (enableSolvent) {
+                    surfaceVolume *= (Toolbox::template decay<LhsEval>( intQuants.solventSaturation())
+                    *  Toolbox::template decay<LhsEval>(intQuants.solventInverseFormationVolumeFactor()));
+                }
+            } else {
+                throw std::runtime_error("Transport phase is GAS/WATER/SOLVENT");
+            }
 
             // Avoid singular matrix if no gas is present.
-            surfaceVolumeFreeGas = max(surfaceVolumeFreeGas, 1e-10);
+            surfaceVolume = max(surfaceVolume, 1e-10);
 
-            // Foam/surfactant in gas phase.
+            // Foam/surfactant in free phase.
-            const LhsEval gasFoam = surfaceVolumeFreeGas
+            const LhsEval freeFoam = surfaceVolume
                 * Toolbox::template decay<LhsEval>(intQuants.foamConcentration());
 
             // Adsorbed foam/surfactant.
@@ -267,7 +278,7 @@ public:
                 * Toolbox::template decay<LhsEval>(intQuants.foamRockDensity())
                 * Toolbox::template decay<LhsEval>(intQuants.foamAdsorbed());
 
-            LhsEval accumulationFoam = gasFoam + adsorbedFoam;
+            LhsEval accumulationFoam = freeFoam + adsorbedFoam;
             storage[contiFoamEqIdx] += accumulationFoam;
         }
     }
@@ -280,23 +291,67 @@ public:
     {
         if constexpr (enableFoam) {
             const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
-            const unsigned upIdx = extQuants.upstreamIndex(FluidSystem::gasPhaseIdx);
             const unsigned inIdx = extQuants.interiorIndex();
-            const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
 
-            // The effect of the gas mobility reduction factor is
+            // The effect of the mobility reduction factor is
-            // incorporated in the mobility, so the oil (if vaporized oil
+            // incorporated in the mobility for the relevant phase,
-            // is active) and gas fluxes do not need modification here.
+            // so fluxes do not need modification here.
-            if (upIdx == inIdx) {
+            switch (transportPhase()) {
-                flux[contiFoamEqIdx] =
+                case Phase::WATER: {
-                    extQuants.volumeFlux(gasPhaseIdx)
+                    const unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
-                    *up.fluidState().invB(gasPhaseIdx)
+                    const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
-                    *up.foamConcentration();
+                    if (upIdx == inIdx) {
-            } else {
+                        flux[contiFoamEqIdx] =
-                flux[contiFoamEqIdx] =
+                            extQuants.volumeFlux(waterPhaseIdx)
-                    extQuants.volumeFlux(gasPhaseIdx)
+                            *up.fluidState().invB(waterPhaseIdx)
-                    *decay<Scalar>(up.fluidState().invB(gasPhaseIdx))
+                            *up.foamConcentration();
-                    *decay<Scalar>(up.foamConcentration());
+                    } else {
+                        flux[contiFoamEqIdx] =
+                            extQuants.volumeFlux(waterPhaseIdx)
+                            *decay<Scalar>(up.fluidState().invB(waterPhaseIdx))
+                            *decay<Scalar>(up.foamConcentration());
+                    }
+                    break;
+                }
+                case Phase::GAS: {
+                    const unsigned upIdx = extQuants.upstreamIndex(gasPhaseIdx);
+                    const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
+                    if (upIdx == inIdx) {
+                        flux[contiFoamEqIdx] =
+                            extQuants.volumeFlux(gasPhaseIdx)
+                            *up.fluidState().invB(gasPhaseIdx)
+                            *up.foamConcentration();
+                    } else {
+                        flux[contiFoamEqIdx] =
+                            extQuants.volumeFlux(gasPhaseIdx)
+                            *decay<Scalar>(up.fluidState().invB(gasPhaseIdx))
+                            *decay<Scalar>(up.foamConcentration());
+                    }
+                    break;
+                }
+                case Phase::SOLVENT: {
+                    if constexpr (enableSolvent) {
+                        const unsigned upIdx = extQuants.solventUpstreamIndex();
+                        const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
+                        if (upIdx == inIdx) {
+                            flux[contiFoamEqIdx] =
+                                extQuants.solventVolumeFlux()
+                                *up.solventInverseFormationVolumeFactor()
+                                *up.foamConcentration();
+                        } else {
+                            flux[contiFoamEqIdx] =
+                                extQuants.solventVolumeFlux()
+                                *decay<Scalar>(up.solventInverseFormationVolumeFactor())
+                                *decay<Scalar>(up.foamConcentration());
+                        }
+                    } else {
+                        throw std::runtime_error("Foam transport phase is SOLVENT but SOLVENT is not activated.");
+                    }
+                    break;
+                }
+                default: {
+                    throw std::runtime_error("Foam transport phase must be GAS/WATER/SOLVENT.");
+                }
             }
         }
     }
@@ -382,6 +437,10 @@ public:
         return params_.foamCoefficients_[satnumRegionIdx];
     }
 
+    static Phase transportPhase() {
+        return params_.transport_phase_;
+    }
+
 private:
     static BlackOilFoamParams<Scalar> params_;
 };
@@ -413,6 +472,8 @@ class BlackOilFoamIntensiveQuantities
     using FoamModule = BlackOilFoamModule<TypeTag>;
 
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
+    enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
+
     static constexpr int foamConcentrationIdx = Indices::foamConcentrationIdx;
     static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
     static constexpr unsigned oilPhaseIdx = FluidSystem::oilPhaseIdx;
@@ -475,8 +536,28 @@ public:
             mobilityReductionFactor = gasMobilityMultiplier.eval(foamConcentration_, /* extrapolate = */ true);
         }
 
-        // adjust gas mobility
+        // adjust mobility
-        asImp_().mobility_[gasPhaseIdx] *= mobilityReductionFactor;
+        switch (FoamModule::transportPhase()) {
+            case Phase::WATER: {
+                asImp_().mobility_[waterPhaseIdx] *= mobilityReductionFactor;
+                break;
+            }
+            case Phase::GAS: {
+                asImp_().mobility_[gasPhaseIdx] *= mobilityReductionFactor;
+                break;
+            }
+            case Phase::SOLVENT: {
+                if constexpr (enableSolvent) {
+                    asImp_().solventMobility_ *= mobilityReductionFactor;
+                } else {
+                    throw std::runtime_error("Foam transport phase is SOLVENT but SOLVENT is not activated.");
+                }
+                break;
+            }
+            default: {
+                throw std::runtime_error("Foam transport phase must be GAS/WATER/SOLVENT.");
+            }
+        }
 
         foamRockDensity_ = FoamModule::foamRockDensity(elemCtx, dofIdx, timeIdx);
 

opm/models/blackoil/blackoilfoamparams.hh

@@ -78,6 +78,7 @@ struct BlackOilFoamParams {
     std::vector<FoamCoefficients> foamCoefficients_;
     std::vector<TabulatedFunction> adsorbedFoamTable_;
     std::vector<TabulatedFunction> gasMobilityMultiplierTable_;
+    Opm::Phase transport_phase_;
 };
 
 } // namespace Opm