diff --git a/opm/models/blackoil/blackoilbrinemodules.hh b/opm/models/blackoil/blackoilbrinemodules.hh index 09aeec418..2dd5a5f06 100644 --- a/opm/models/blackoil/blackoilbrinemodules.hh +++ b/opm/models/blackoil/blackoilbrinemodules.hh @@ -38,6 +38,10 @@ #include #include #include +#include +#include +#include +#include #endif #include @@ -78,8 +82,10 @@ class BlackOilBrineModule static constexpr unsigned saltConcentrationIdx = Indices::saltConcentrationIdx; static constexpr unsigned contiBrineEqIdx = Indices::contiBrineEqIdx; static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx; - + static const bool gasEnabled = Indices::gasEnabled; + static const bool oilEnabled = Indices::oilEnabled; static constexpr unsigned enableBrine = enableBrineV; + static constexpr unsigned enableSaltPrecipitation = getPropValue(); static constexpr unsigned numEq = getPropValue(); static constexpr unsigned numPhases = FluidSystem::numPhases; @@ -103,7 +109,6 @@ public: "contains the BRINE keyword"); } - if (!eclState.runspec().phases().active(Phase::BRINE)) return; // brine treatment is supposed to be disabled @@ -126,6 +131,25 @@ public: bdensityTable_[pvtRegionIdx].setXYContainers(c, bdensityTable); } } + + if (enableSaltPrecipitation) { + const TableContainer& permfactTables = tableManager.getPermfactTables(); + permfactTable_.resize(numPvtRegions); + for (size_t i = 0; i < permfactTables.size(); ++i) { + const PermfactTable& permfactTable = permfactTables.getTable(i); + permfactTable_[i].setXYContainers(permfactTable.getPorosityChangeColumn(), permfactTable.getPermeabilityMultiplierColumn()); + } + + const TableContainer& saltsolTables = tableManager.getSaltsolTables(); + if (!saltsolTables.empty()) { + saltsolTable_.resize(numPvtRegions); + assert(numPvtRegions == saltsolTables.size()); + for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) { + const SaltsolTable& saltsolTable = saltsolTables.getTable(pvtRegionIdx ); + saltsolTable_[pvtRegionIdx] = saltsolTable.getSaltsolColumn().front(); + } + } + } } #endif @@ -139,7 +163,6 @@ public: return; } - static bool primaryVarApplies(unsigned pvIdx) { if (!enableBrine) @@ -222,7 +245,16 @@ public: const LhsEval massBrine = surfaceVolumeWater * Toolbox::template decay(fs.saltConcentration()); - storage[contiBrineEqIdx] += massBrine; + if (enableSaltPrecipitation){ + double saltDensity = 2170; // Solid salt density kg/m3 + const LhsEval solidSalt = + Toolbox::template decay(intQuants.porosity()) + * saltDensity + * Toolbox::template decay(intQuants.saltSaturation()); + + storage[contiBrineEqIdx] += massBrine + solidSalt; + } + else { storage[contiBrineEqIdx] += massBrine;} } static void computeFlux(RateVector& flux, @@ -234,7 +266,12 @@ 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(); @@ -245,12 +282,24 @@ 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(up.fluidState().invB(waterPhaseIdx)) *decay(up.fluidState().saltConcentration()); + + if (enableSaltPrecipitation) { + // modify gas and oil flux for mobility change + if(gasEnabled) { flux[contiGasEqIdx] *= decay(intQuants.permFactor()); } + if(oilEnabled) { flux[contiOilEqIdx] *= decay(intQuants.permFactor()); } + } } } @@ -310,13 +359,40 @@ public: return bdensityTable_[pvtnumRegionIdx]; } + static const TabulatedFunction& permfactTable(const ElementContext& elemCtx, + unsigned scvIdx, + unsigned timeIdx) + { + unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx); + return permfactTable_[pvtnumRegionIdx]; + } + + static const Scalar saltsolTable(const ElementContext& elemCtx, + unsigned scvIdx, + unsigned timeIdx) + { + unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx); + return saltsolTable_[pvtnumRegionIdx]; + } + static bool hasBDensityTables() { return !bdensityTable_.empty(); } + static bool hasSaltsolTables() + { + return !saltsolTable_.empty(); + } + + static Scalar saltSol(unsigned regionIdx) { + return saltsolTable_[regionIdx]; + } + private: static std::vector bdensityTable_; + static std::vector permfactTable_; + static std::vector saltsolTable_; static std::vector referencePressure_; }; @@ -329,6 +405,14 @@ template std::vector::Scalar> BlackOilBrineModule::referencePressure_; +template +std::vector::Scalar> +BlackOilBrineModule::saltsolTable_; + +template +std::vector::TabulatedFunction> +BlackOilBrineModule::permfactTable_; + /*! * \ingroup BlackOil * \class Ewoms::BlackOilBrineIntensiveQuantities @@ -354,8 +438,10 @@ class BlackOilBrineIntensiveQuantities enum { numPhases = getPropValue() }; static constexpr int saltConcentrationIdx = Indices::saltConcentrationIdx; static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx; + static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx; static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx; static constexpr unsigned enableBrine = enableBrineV; + static constexpr unsigned enableSaltPrecipitation = getPropValue(); static constexpr int contiBrineEqIdx = Indices::contiBrineEqIdx; public: @@ -372,9 +458,42 @@ public: const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx); auto& fs = asImp_().fluidState_; - // set saltconcentration - fs.setSaltConcentration(priVars.makeEvaluation(saltConcentrationIdx, timeIdx)); + if (enableSaltPrecipitation) { + const auto& saltsolTable = BrineModule::saltsolTable(elemCtx, dofIdx, timeIdx); + saltSolubility_ = saltsolTable; + + if (priVars.primaryVarsMeaningBrine() == PrimaryVariables::Sp) { + saltSaturation_ = priVars.makeEvaluation(saltConcentrationIdx, timeIdx); + fs.setSaltConcentration(saltSolubility_); + } + else { + saltConcentration_ = priVars.makeEvaluation(saltConcentrationIdx, timeIdx); + fs.setSaltConcentration(saltConcentration_); + saltSaturation_ = 0.0; + } + } + else { + saltConcentration_ = priVars.makeEvaluation(saltConcentrationIdx, timeIdx); + fs.setSaltConcentration(saltConcentration_); + } + } + + void saltPropertiesUpdate_(const ElementContext& elemCtx, + unsigned dofIdx, + unsigned timeIdx) + { + if (enableSaltPrecipitation) { + const Evaluation porosityFactor = 1.0 - saltSaturation(); //phi/phi_0 + + const auto& permfactTable = BrineModule::permfactTable(elemCtx, dofIdx, timeIdx); + + permFactor_ = permfactTable.eval(scalarValue(porosityFactor)); + if (permFactor_ < 1 ) { + // adjust mobility for changing permeability + asImp_().mobility_[waterPhaseIdx] *= permFactor_; + } + } } const Evaluation& saltConcentration() const @@ -383,12 +502,24 @@ public: const Evaluation& brineRefDensity() const { return refDensity_; } + const Evaluation& saltSaturation() const + { return saltSaturation_; } + + Scalar saltSolubility() const + { return saltSolubility_; } + + const Evaluation& permFactor() const + { return permFactor_; } + protected: Implementation& asImp_() { return *static_cast(this); } Evaluation saltConcentration_; Evaluation refDensity_; + Evaluation saltSaturation_; + Evaluation permFactor_; + Scalar saltSolubility_; }; @@ -405,18 +536,28 @@ public: unsigned) { } + void saltPropertiesUpdate_(const ElementContext&, + unsigned, + unsigned) + { } + const Evaluation& saltConcentration() const { throw std::runtime_error("saltConcentration() called but brine are disabled"); } const Evaluation& brineRefDensity() const { throw std::runtime_error("brineRefDensity() called but brine are disabled"); } + const Evaluation& saltSaturation() const + { throw std::logic_error("saltSaturation() called but salt precipitation is disabled"); } + + const Scalar saltSolubility() const + { throw std::logic_error("saltSolubility() called but salt precipitation is disabled"); } + + const Evaluation& permFactor() const + { throw std::logic_error("permFactor() called but salt precipitation is disabled"); } + }; - - - - } // namespace Ewoms #endif diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 7e5e119a9..ff4f602ba 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -403,6 +403,7 @@ public: asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache); asImp_().foamPropertiesUpdate_(elemCtx, dofIdx, timeIdx); asImp_().MICPPropertiesUpdate_(elemCtx, dofIdx, timeIdx); + asImp_().saltPropertiesUpdate_(elemCtx, dofIdx, timeIdx); // update the quantities which are required by the chosen // velocity model diff --git a/opm/models/blackoil/blackoilnewtonmethod.hh b/opm/models/blackoil/blackoilnewtonmethod.hh index 8041796be..6e2d8c496 100644 --- a/opm/models/blackoil/blackoilnewtonmethod.hh +++ b/opm/models/blackoil/blackoilnewtonmethod.hh @@ -118,6 +118,7 @@ class BlackOilNewtonMethod : public GetPropType; static const unsigned numEq = getPropValue(); + static constexpr bool enableSaltPrecipitation = getPropValue(); public: BlackOilNewtonMethod(Simulator& simulator) : ParentType(simulator) @@ -369,8 +370,10 @@ protected: nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0); // keep the salt concentration above 0 - if (enableBrine && pvIdx == Indices::saltConcentrationIdx) - nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0); + if (enableBrine && pvIdx == Indices::saltConcentrationIdx) { + if (!enableSaltPrecipitation || (enableSaltPrecipitation && currentValue.primaryVarsMeaningBrine() == PrimaryVariables::Cs)) + nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0); + } // keep the temperature within given values if (enableEnergy && pvIdx == Indices::temperatureIdx) diff --git a/opm/models/blackoil/blackoilprimaryvariables.hh b/opm/models/blackoil/blackoilprimaryvariables.hh index 12971b8f2..39b648789 100644 --- a/opm/models/blackoil/blackoilprimaryvariables.hh +++ b/opm/models/blackoil/blackoilprimaryvariables.hh @@ -82,6 +82,7 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables enum { waterSaturationIdx = Indices::waterSaturationIdx }; enum { pressureSwitchIdx = Indices::pressureSwitchIdx }; enum { compositionSwitchIdx = Indices::compositionSwitchIdx }; + enum { saltConcentrationIdx = Indices::saltConcentrationIdx }; static const bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0; static const bool waterEnabled = Indices::waterEnabled; @@ -101,6 +102,7 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables enum { enablePolymer = getPropValue() }; enum { enableFoam = getPropValue() }; enum { enableBrine = getPropValue() }; + enum { enableSaltPrecipitation = getPropValue() }; enum { enableEnergy = getPropValue() }; enum { enableMICP = getPropValue() }; enum { gasCompIdx = FluidSystem::gasCompIdx }; @@ -128,6 +130,11 @@ public: OnePhase_p, // onephase case }; + enum PrimaryVarsMeaningBrine { + Cs, // salt concentration + Sp, // (precipitated) salt saturation + }; + BlackOilPrimaryVariables() : ParentType() { @@ -181,6 +188,17 @@ public: void setPrimaryVarsMeaning(PrimaryVarsMeaning newMeaning) { primaryVarsMeaning_ = newMeaning; } + PrimaryVarsMeaningBrine primaryVarsMeaningBrine() const + { return primaryVarsMeaningBrine_; } + + /*! + * \brief Set the interpretation which should be applied to the switching primary + * variables. + */ + + void setPrimaryVarsMeaningBrine(PrimaryVarsMeaningBrine newMeaning) + { primaryVarsMeaningBrine_ = newMeaning; } + /*! * \copydoc ImmisciblePrimaryVariables::assignMassConservative */ @@ -273,6 +291,8 @@ public: bool oilPresent = FluidSystem::phaseIsActive(oilPhaseIdx)?(fluidState.saturation(oilPhaseIdx) > 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_(fluidState, pvtRegionIdx_); + bool precipitatedSaltPresent = enableSaltPrecipitation?(saltSaturation > 0.0):false; // deal with the primary variables for the energy extension EnergyModule::assignPrimaryVars(*this, fluidState); @@ -304,6 +324,13 @@ public: primaryVarsMeaning_ = Sw_po_Sg; } + if (enableSaltPrecipitation){ + if (precipitatedSaltPresent) + primaryVarsMeaningBrine_ = Sp; + else + primaryVarsMeaningBrine_ = Cs; + } + // assign the actual primary variables if (primaryVarsMeaning() == OnePhase_p) { if (waterEnabled) { @@ -311,7 +338,6 @@ public: } else { throw std::logic_error("For single-phase runs, only pure water is presently allowed."); } - } else if (primaryVarsMeaning() == Sw_po_Sg) { if (waterEnabled) @@ -347,6 +373,16 @@ public: (*this)[compositionSwitchIdx] = Rv; } + if (enableSaltPrecipitation) { + if (primaryVarsMeaningBrine() == Sp) { + (*this)[saltConcentrationIdx] = FsToolbox::value(saltSaturation); + } + else { + const auto& saltConcentration = BlackOil::getSaltConcentration_(fluidState, pvtRegionIdx_); + (*this)[saltConcentrationIdx] = FsToolbox::value(saltConcentration); + } + } + checkDefined(); } @@ -379,6 +415,24 @@ public: if (waterEnabled) Sw = (*this)[Indices::waterSaturationIdx]; + if (enableSaltPrecipitation) { + Scalar saltSolubility = BrineModule::saltSol(pvtRegionIndex()); + if (primaryVarsMeaningBrine() == Sp) { + Scalar saltSat = (*this)[saltConcentrationIdx]; + if (saltSat < -eps){ //precipitated salt dissappears + setPrimaryVarsMeaningBrine(Cs); + (*this)[saltConcentrationIdx] = saltSolubility; //set salt concentration to solubility limit + } + } + else if (primaryVarsMeaningBrine() == Cs) { + Scalar saltConc = (*this)[saltConcentrationIdx]; + if (saltConc > saltSolubility + eps){ //salt concentration exceeds solubility limit + setPrimaryVarsMeaningBrine(Sp); + (*this)[saltConcentrationIdx] = 0.0; + } + } + } + if (primaryVarsMeaning() == Sw_po_Sg) { // special case for cells with almost only water @@ -816,6 +870,7 @@ private: } PrimaryVarsMeaning primaryVarsMeaning_; + PrimaryVarsMeaningBrine primaryVarsMeaningBrine_; unsigned short pvtRegionIdx_; };