28 #ifndef EWOMS_BLACK_OIL_POLYMER_MODULE_HH 29 #define EWOMS_BLACK_OIL_POLYMER_MODULE_HH 35 #include <opm/material/common/Tabulated1DFunction.hpp> 38 #include <opm/parser/eclipse/Deck/Deck.hpp> 39 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp> 40 #include <opm/parser/eclipse/EclipseState/Tables/PlyadsTable.hpp> 41 #include <opm/parser/eclipse/EclipseState/Tables/PlymaxTable.hpp> 42 #include <opm/parser/eclipse/EclipseState/Tables/PlyrockTable.hpp> 43 #include <opm/parser/eclipse/EclipseState/Tables/PlyshlogTable.hpp> 44 #include <opm/parser/eclipse/EclipseState/Tables/PlyviscTable.hpp> 45 #include <opm/parser/eclipse/EclipseState/Tables/PlyshlogTable.hpp> 48 #include <opm/common/Valgrind.hpp> 49 #include <opm/common/Unused.hpp> 50 #include <opm/common/ErrorMacros.hpp> 51 #include <opm/common/Exceptions.hpp> 53 #include <dune/common/fvector.hh> 63 template <
class TypeTag,
bool enablePolymerV = GET_PROP_VALUE(TypeTag, EnablePolymer)>
66 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
67 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
68 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
69 typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
70 typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
71 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
72 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
73 typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
75 typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
76 typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
77 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
79 typedef Opm::MathToolbox<Evaluation> Toolbox;
81 typedef typename Opm::Tabulated1DFunction<Scalar> TabulatedFunction;
83 static constexpr
unsigned polymerConcentrationIdx = Indices::polymerConcentrationIdx;
84 static constexpr
unsigned contiPolymerEqIdx = Indices::contiPolymerEqIdx;
85 static constexpr
unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
88 static constexpr
unsigned enablePolymer = enablePolymerV;
90 static constexpr
unsigned numPhases = FluidSystem::numPhases;
93 enum AdsorptionBehaviour { Desorption = 1, NoDesorption = 2 };
99 static void initFromDeck(
const Opm::Deck& deck,
const Opm::EclipseState& eclState)
103 if (enablePolymer && !deck.hasKeyword(
"POLYMER")) {
104 OPM_THROW(std::runtime_error,
105 "Non-trivial polymer treatment requested at compile time, but " 106 "the deck does not contain the POLYMER keyword");
108 else if (!enablePolymer && deck.hasKeyword(
"POLYMER")) {
109 OPM_THROW(std::runtime_error,
110 "Polymer treatment disabled at compile time, but the deck " 111 "contains the POLYMER keyword");
114 if (!deck.hasKeyword(
"POLYMER"))
117 const auto& tableManager = eclState.getTableManager();
119 unsigned numSatRegions = tableManager.getTabdims().getNumSatTables();
123 const auto& plyrockTables = tableManager.getPlyrockTables();
124 if (!plyrockTables.empty()) {
125 assert(numSatRegions == plyrockTables.size());
126 for (
unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++ satRegionIdx) {
127 const auto& plyrockTable = plyrockTables.template getTable<Opm::PlyrockTable>(satRegionIdx);
129 plyrockTable.getDeadPoreVolumeColumn()[satRegionIdx],
130 plyrockTable.getResidualResistanceFactorColumn()[satRegionIdx],
131 plyrockTable.getRockDensityFactorColumn()[satRegionIdx],
132 static_cast<AdsorptionBehaviour
>(plyrockTable.getAdsorbtionIndexColumn()[satRegionIdx]),
133 plyrockTable.getMaxAdsorbtionColumn()[satRegionIdx]);
136 OPM_THROW(std::runtime_error,
"PLYROCK must be specified in POLYMER runs\n");
140 const auto& plyadsTables = tableManager.getPlyadsTables();
141 if (!plyadsTables.empty()) {
142 assert(numSatRegions == plyadsTables.size());
143 for (
unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++ satRegionIdx) {
144 const auto& plyadsTable = plyadsTables.template getTable<Opm::PlyadsTable>(satRegionIdx);
146 const auto& c = plyadsTable.getPolymerConcentrationColumn();
147 const auto& ads = plyadsTable.getAdsorbedPolymerColumn();
148 plyadsAdsorbedPolymer_[satRegionIdx].setXYContainers(c, ads);
151 OPM_THROW(std::runtime_error,
"PLYADS must be specified in POLYMER runs\n");
155 unsigned numPvtRegions = tableManager.getTabdims().getNumPVTTables();
159 const auto& plyviscTables = tableManager.getPlyviscTables();
160 if (!plyviscTables.empty()) {
162 assert(numPvtRegions == plyviscTables.size());
163 for (
unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) {
164 const auto& plyadsTable = plyviscTables.template getTable<Opm::PlyviscTable>(pvtRegionIdx);
166 const auto& c = plyadsTable.getPolymerConcentrationColumn();
167 const auto& visc = plyadsTable.getViscosityMultiplierColumn();
168 plyviscViscosityMultiplierTable_[pvtRegionIdx].setXYContainers(c, visc);
171 OPM_THROW(std::runtime_error,
"PLYVISC must be specified in POLYMER runs\n");
175 const auto& plymaxTables = tableManager.getPlymaxTables();
176 unsigned numMixRegions = plymaxTables.size();
178 if (!plymaxTables.empty()) {
179 for (
unsigned mixRegionIdx = 0; mixRegionIdx < numMixRegions; ++ mixRegionIdx) {
180 const auto& plymaxTable = plymaxTables.template getTable<Opm::PlymaxTable>(mixRegionIdx);
181 setPlymax(mixRegionIdx, plymaxTable.getPolymerConcentrationColumn()[mixRegionIdx]);
184 OPM_THROW(std::runtime_error,
"PLYMAX must be specified in POLYMER runs\n");
187 if (deck.hasKeyword(
"PLMIXPAR")) {
189 for (
unsigned mixRegionIdx = 0; mixRegionIdx < numMixRegions; ++ mixRegionIdx) {
190 const auto& plmixparRecord = deck.getKeyword(
"PLMIXPAR").getRecord(mixRegionIdx);
191 setPlmixpar(mixRegionIdx, plmixparRecord.getItem(
"TODD_LONGSTAFF").getSIDouble(0));
194 OPM_THROW(std::runtime_error,
"PLMIXPAR must be specified in POLYMER runs\n");
197 hasPlyshlog_ = deck.hasKeyword(
"PLYSHLOG");
198 hasShrate_ = deck.hasKeyword(
"SHRATE");
201 const auto& plyshlogTables = tableManager.getPlyshlogTables();
202 assert(numPvtRegions == plyshlogTables.size());
203 plyshlogShearEffectRefMultiplier_.resize(numPvtRegions);
204 plyshlogShearEffectRefLogVelocity_.resize(numPvtRegions);
205 for (
unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) {
206 const auto& plyshlogTable = plyshlogTables.template getTable<Opm::PlyshlogTable>(pvtRegionIdx);
208 Scalar plyshlogRefPolymerConcentration = plyshlogTable.getRefPolymerConcentration();
209 std::vector<Scalar> waterVelocity = plyshlogTable.getWaterVelocityColumn().vectorCopy();
210 std::vector<Scalar> shearMultiplier = plyshlogTable.getShearMultiplierColumn().vectorCopy();
213 Opm::UnitSystem unitSystem = deck.getActiveUnitSystem();
214 double siFactor = hasShrate_? unitSystem.parse(
"1/Time").getSIScaling() : unitSystem.parse(
"Length/Time").getSIScaling();
215 for (
size_t i = 0; i < waterVelocity.size(); ++i ) {
216 waterVelocity[i] *= siFactor;
219 waterVelocity[i] = std::log(waterVelocity[i]);
222 Scalar refViscMult = plyviscViscosityMultiplierTable_[pvtRegionIdx].eval(plyshlogRefPolymerConcentration,
true);
224 for (
size_t i = 0; i < waterVelocity.size(); ++i ) {
225 shearMultiplier[i] *= refViscMult;
226 shearMultiplier[i] -= 1;
227 shearMultiplier[i] /= (refViscMult - 1);
228 shearMultiplier[i] = shearMultiplier[i];
230 plyshlogShearEffectRefMultiplier_[pvtRegionIdx].resize(waterVelocity.size());
231 plyshlogShearEffectRefLogVelocity_[pvtRegionIdx].resize(waterVelocity.size());
233 for (
size_t i = 0; i < waterVelocity.size(); ++i ) {
234 plyshlogShearEffectRefMultiplier_[pvtRegionIdx][i] = shearMultiplier[i];
235 plyshlogShearEffectRefLogVelocity_[pvtRegionIdx][i] = waterVelocity[i];
242 OPM_THROW(std::runtime_error,
"PLYSHLOG must be specified if SHRATE is used in POLYMER runs\n");
244 const auto& shrateKeyword = deck.getKeyword(
"SHRATE");
245 const std::vector<double>& shrateFromDeck = shrateKeyword.getSIDoubleData();
246 shrate_.resize(numPvtRegions);
247 for (
unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) {
248 if (shrateFromDeck.empty()) {
249 shrate_[pvtRegionIdx] = 4.8;
250 }
else if (shrateFromDeck.size() == numPvtRegions) {
251 shrate_[pvtRegionIdx] = shrateKeyword.getSIDoubleData()[pvtRegionIdx];
253 OPM_THROW(std::runtime_error,
"SHRATE must either have 0 or number of NUMPVT entries\n");
268 plyrockDeadPoreVolume_.resize(numRegions);
269 plyrockResidualResistanceFactor_.resize(numRegions);
270 plyrockRockDensityFactor_.resize(numRegions);
271 plyrockAdsorbtionIndex_.resize(numRegions);
272 plyrockMaxAdsorbtion_.resize(numRegions);
273 plyadsAdsorbedPolymer_.resize(numRegions);
282 const Scalar& plyrockDeadPoreVolume,
283 const Scalar& plyrockResidualResistanceFactor,
284 const Scalar& plyrockRockDensityFactor,
285 const Scalar& plyrockAdsorbtionIndex,
286 const Scalar& plyrockMaxAdsorbtion)
288 plyrockDeadPoreVolume_[satRegionIdx] = plyrockDeadPoreVolume;
289 plyrockResidualResistanceFactor_[satRegionIdx] = plyrockResidualResistanceFactor;
290 plyrockRockDensityFactor_[satRegionIdx] = plyrockRockDensityFactor;
291 plyrockAdsorbtionIndex_[satRegionIdx] = plyrockAdsorbtionIndex;
292 plyrockMaxAdsorbtion_[satRegionIdx] = plyrockMaxAdsorbtion;
301 const TabulatedFunction& plyadsAdsorbedPolymer)
303 plyadsAdsorbedPolymer_[satRegionIdx] = plyadsAdsorbedPolymer;
313 plyviscViscosityMultiplierTable_.resize(numRegions);
322 const TabulatedFunction& plyviscViscosityMultiplierTable)
324 plyviscViscosityMultiplierTable_[satRegionIdx] = plyviscViscosityMultiplierTable;
334 plymaxMaxConcentration_.resize(numRegions);
335 plymixparToddLongstaff_.resize(numRegions);
344 const Scalar& plymaxMaxConcentration)
346 plymaxMaxConcentration_[mixRegionIdx] = plymaxMaxConcentration;
355 const Scalar& plymixparToddLongstaff)
357 plymixparToddLongstaff_[mixRegionIdx] = plymixparToddLongstaff;
387 static bool primaryVarApplies(
unsigned pvIdx)
393 return pvIdx == polymerConcentrationIdx;
396 static std::string primaryVarName(
unsigned pvIdx OPM_OPTIM_UNUSED)
398 assert(primaryVarApplies(pvIdx));
400 return "polymer_waterconcentration";
403 static Scalar primaryVarWeight(
unsigned pvIdx OPM_OPTIM_UNUSED)
405 assert(primaryVarApplies(pvIdx));
408 return static_cast<Scalar
>(1.0);
411 static bool eqApplies(
unsigned eqIdx)
416 return eqIdx == contiPolymerEqIdx;
419 static std::string eqName(
unsigned eqIdx OPM_OPTIM_UNUSED)
421 assert(eqApplies(eqIdx));
423 return "conti^polymer";
426 static Scalar eqWeight(
unsigned eqIdx OPM_OPTIM_UNUSED)
428 assert(eqApplies(eqIdx));
431 return static_cast<Scalar
>(1.0);
435 template <
class LhsEval>
436 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
437 const IntensiveQuantities& intQuants)
442 const auto& fs = intQuants.fluidState();
444 LhsEval surfaceVolumeWater =
445 Toolbox::template decay<LhsEval>(fs.saturation(waterPhaseIdx))
446 * Toolbox::template decay<LhsEval>(fs.invB(waterPhaseIdx))
447 * Toolbox::template decay<LhsEval>(intQuants.porosity());
450 surfaceVolumeWater = Opm::max(surfaceVolumeWater, 1e-10);
453 storage[contiPolymerEqIdx] += surfaceVolumeWater
454 * Toolbox::template decay<LhsEval>(intQuants.polymerConcentration())
455 * (1.0 - Toolbox::template decay<LhsEval>(intQuants.polymerDeadPoreVolume()));
458 storage[contiPolymerEqIdx] +=
459 Toolbox::template decay<LhsEval>(1.0 - intQuants.porosity())
460 * Toolbox::template decay<LhsEval>(intQuants.polymerRockDensity())
461 * Toolbox::template decay<LhsEval>(intQuants.polymerAdsorption());
465 static void computeFlux(RateVector& flux,
466 const ElementContext& elemCtx,
474 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
476 unsigned upIdx = extQuants.upstreamIndex(FluidSystem::waterPhaseIdx);
477 unsigned inIdx = extQuants.interiorIndex();
478 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
479 const unsigned contiWaterEqIdx = Indices::conti0EqIdx + Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
482 if (upIdx == inIdx) {
483 flux[contiPolymerEqIdx] =
484 extQuants.volumeFlux(waterPhaseIdx)
485 *up.fluidState().invB(waterPhaseIdx)
486 *up.polymerViscosityCorrection()
487 /extQuants.polymerShearFactor()
488 *up.polymerConcentration();
491 flux[contiWaterEqIdx] /=
492 extQuants.waterShearFactor();
494 flux[contiPolymerEqIdx] =
495 extQuants.volumeFlux(waterPhaseIdx)
496 *Opm::decay<Scalar>(up.fluidState().invB(waterPhaseIdx))
497 *Opm::decay<Scalar>(up.polymerViscosityCorrection())
498 /Opm::decay<Scalar>(extQuants.polymerShearFactor())
499 *Opm::decay<Scalar>(up.polymerConcentration());
502 flux[contiWaterEqIdx] /=
503 Opm::decay<Scalar>(extQuants.waterShearFactor());
512 Scalar polymerConcentration)
517 priVars[polymerConcentrationIdx] = polymerConcentration;
524 const PrimaryVariables& oldPv,
525 const EqVector& delta)
531 newPv[polymerConcentrationIdx] = oldPv[polymerConcentrationIdx] - delta[polymerConcentrationIdx];
538 const EqVector& delta OPM_UNUSED)
543 return static_cast<Scalar
>(0.0);
552 return std::abs(Toolbox::scalarValue(resid[contiPolymerEqIdx]));
555 template <
class DofEntity>
556 static void serializeEntity(
const Model& model, std::ostream& outstream,
const DofEntity& dof)
561 unsigned dofIdx = model.dofMapper().index(dof);
562 const PrimaryVariables& priVars = model.solution(0)[dofIdx];
563 outstream << priVars[polymerConcentrationIdx];
566 template <
class DofEntity>
567 static void deserializeEntity(Model& model, std::istream& instream,
const DofEntity& dof)
572 unsigned dofIdx = model.dofMapper().index(dof);
573 PrimaryVariables& priVars0 = model.solution(0)[dofIdx];
574 PrimaryVariables& priVars1 = model.solution(1)[dofIdx];
576 instream >> priVars0[polymerConcentrationIdx];
579 priVars1 = priVars0[polymerConcentrationIdx];
582 static const Scalar plyrockDeadPoreVolume(
const ElementContext& elemCtx,
586 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
587 return plyrockDeadPoreVolume_[satnumRegionIdx];
590 static const Scalar plyrockResidualResistanceFactor(
const ElementContext& elemCtx,
594 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
595 return plyrockResidualResistanceFactor_[satnumRegionIdx];
598 static const Scalar plyrockRockDensityFactor(
const ElementContext& elemCtx,
602 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
603 return plyrockRockDensityFactor_[satnumRegionIdx];
606 static const Scalar plyrockAdsorbtionIndex(
const ElementContext& elemCtx,
610 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
611 return plyrockAdsorbtionIndex_[satnumRegionIdx];
614 static const Scalar plyrockMaxAdsorbtion(
const ElementContext& elemCtx,
618 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
619 return plyrockMaxAdsorbtion_[satnumRegionIdx];
622 static const TabulatedFunction& plyadsAdsorbedPolymer(
const ElementContext& elemCtx,
626 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
627 return plyadsAdsorbedPolymer_[satnumRegionIdx];
630 static const TabulatedFunction& plyviscViscosityMultiplierTable(
const ElementContext& elemCtx,
634 unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
635 return plyviscViscosityMultiplierTable_[pvtnumRegionIdx];
638 static const TabulatedFunction& plyviscViscosityMultiplierTable(
unsigned pvtnumRegionIdx)
640 return plyviscViscosityMultiplierTable_[pvtnumRegionIdx];
643 static const Scalar plymaxMaxConcentration(
const ElementContext& elemCtx,
647 unsigned polymerMixRegionIdx = elemCtx.problem().plmixnumRegionIndex(elemCtx, scvIdx, timeIdx);
648 return plymaxMaxConcentration_[polymerMixRegionIdx];
651 static const Scalar plymixparToddLongstaff(
const ElementContext& elemCtx,
655 unsigned polymerMixRegionIdx = elemCtx.problem().plmixnumRegionIndex(elemCtx, scvIdx, timeIdx);
656 return plymixparToddLongstaff_[polymerMixRegionIdx];
659 static bool hasPlyshlog()
664 static bool hasShrate()
669 static const Scalar shrate(
unsigned pvtnumRegionIdx)
671 return shrate_[pvtnumRegionIdx];
680 template <
class Evaluation>
682 unsigned pvtnumRegionIdx,
683 const Evaluation& v0) {
685 const auto& viscosityMultiplierTable = plyviscViscosityMultiplierTable_[pvtnumRegionIdx];
686 Scalar viscosityMultiplier = viscosityMultiplierTable.eval( Opm::scalarValue(polymerConcentration),
true);
688 const Scalar eps = 1e-14;
690 if ( std::abs( (viscosityMultiplier - 1.0) ) < eps){
694 const std::vector<Scalar>& shearEffectRefLogVelocity = plyshlogShearEffectRefLogVelocity_[pvtnumRegionIdx];
695 auto v0AbsLog = Opm::log(Opm::abs(v0));
697 if (v0AbsLog < shearEffectRefLogVelocity[0])
704 const std::vector<Scalar>& shearEffectRefMultiplier = plyshlogShearEffectRefMultiplier_[pvtnumRegionIdx];
705 size_t numTableEntries = shearEffectRefLogVelocity.size();
706 assert(shearEffectRefMultiplier.size() == numTableEntries);
708 std::vector<Scalar> shearEffectMultiplier(numTableEntries, 1.0);
709 for (
size_t i = 0; i < numTableEntries; ++i) {
710 shearEffectMultiplier[i] = (1.0 + (viscosityMultiplier - 1.0)*shearEffectRefMultiplier[i]) / viscosityMultiplier;
711 shearEffectMultiplier[i] = Opm::log(shearEffectMultiplier[i]);
715 TabulatedFunction logShearEffectMultiplier = TabulatedFunction(numTableEntries, shearEffectRefLogVelocity, shearEffectMultiplier,
false );
722 auto F = [&logShearEffectMultiplier, &v0AbsLog](
const Evaluation& u) {
723 return u + logShearEffectMultiplier.eval(u,
true) - v0AbsLog;
726 auto dF = [&logShearEffectMultiplier](
const Evaluation& u) {
727 return 1 + logShearEffectMultiplier.evalDerivative(u,
true);
733 bool converged =
false;
734 for (
int i = 0; i < 20; ++i ) {
738 if (std::abs(Opm::scalarValue(f)) < 1e-12) {
744 OPM_THROW(std::runtime_error,
"Not able to compute shear velocity. \n");
748 return Opm::exp( logShearEffectMultiplier.eval(u,
true) );
756 static std::vector<Scalar> plyrockDeadPoreVolume_;
757 static std::vector<Scalar> plyrockResidualResistanceFactor_;
758 static std::vector<Scalar> plyrockRockDensityFactor_;
759 static std::vector<Scalar> plyrockAdsorbtionIndex_;
760 static std::vector<Scalar> plyrockMaxAdsorbtion_;
761 static std::vector<TabulatedFunction> plyadsAdsorbedPolymer_;
762 static std::vector<TabulatedFunction> plyviscViscosityMultiplierTable_;
763 static std::vector<Scalar> plymaxMaxConcentration_;
764 static std::vector<Scalar> plymixparToddLongstaff_;
765 static std::vector<std::vector<Scalar>> plyshlogShearEffectRefMultiplier_;
766 static std::vector<std::vector<Scalar>> plyshlogShearEffectRefLogVelocity_;
767 static std::vector<Scalar> shrate_;
768 static bool hasShrate_;
769 static bool hasPlyshlog_;
775 template <
class TypeTag,
bool enablePolymerV>
776 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
777 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyrockDeadPoreVolume_;
779 template <
class TypeTag,
bool enablePolymerV>
780 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
781 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyrockResidualResistanceFactor_;
783 template <
class TypeTag,
bool enablePolymerV>
784 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
785 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyrockRockDensityFactor_;
787 template <
class TypeTag,
bool enablePolymerV>
788 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
789 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyrockAdsorbtionIndex_;
791 template <
class TypeTag,
bool enablePolymerV>
792 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
793 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyrockMaxAdsorbtion_;
795 template <
class TypeTag,
bool enablePolymerV>
796 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::TabulatedFunction>
797 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyadsAdsorbedPolymer_;
799 template <
class TypeTag,
bool enablePolymerV>
800 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::TabulatedFunction>
801 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyviscViscosityMultiplierTable_;
803 template <
class TypeTag,
bool enablePolymerV>
804 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
805 BlackOilPolymerModule<TypeTag, enablePolymerV>::plymaxMaxConcentration_;
807 template <
class TypeTag,
bool enablePolymerV>
808 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
809 BlackOilPolymerModule<TypeTag, enablePolymerV>::plymixparToddLongstaff_;
811 template <
class TypeTag,
bool enablePolymerV>
812 std::vector<std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>>
813 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyshlogShearEffectRefMultiplier_;
815 template <
class TypeTag,
bool enablePolymerV>
816 std::vector<std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>>
817 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyshlogShearEffectRefLogVelocity_;
819 template <
class TypeTag,
bool enablePolymerV>
820 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
821 BlackOilPolymerModule<TypeTag, enablePolymerV>::shrate_;
823 template <
class TypeTag,
bool enablePolymerV>
825 BlackOilPolymerModule<TypeTag, enablePolymerV>::hasShrate_;
827 template <
class TypeTag,
bool enablePolymerV>
829 BlackOilPolymerModule<TypeTag, enablePolymerV>::hasPlyshlog_;
838 template <
class TypeTag,
bool enablePolymerV = GET_PROP_VALUE(TypeTag, EnablePolymer)>
841 typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) Implementation;
843 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
844 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
845 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
846 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
847 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
848 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
849 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
854 static constexpr
int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
855 static constexpr
int waterPhaseIdx = FluidSystem::waterPhaseIdx;
869 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
870 polymerConcentration_ = priVars.makeEvaluation(polymerConcentrationIdx, timeIdx);
871 const Scalar cmax = PolymerModule::plymaxMaxConcentration(elemCtx, dofIdx, timeIdx);
874 const Scalar& maxAdsorbtion = PolymerModule::plyrockMaxAdsorbtion(elemCtx, dofIdx, timeIdx);
875 const auto& plyadsAdsorbedPolymer = PolymerModule::plyadsAdsorbedPolymer(elemCtx, dofIdx, timeIdx);
876 polymerAdsorption_ = plyadsAdsorbedPolymer.eval(polymerConcentration_,
true);
877 if (PolymerModule::plyrockAdsorbtionIndex(elemCtx, dofIdx, timeIdx) == PolymerModule::NoDesorption ) {
878 const Scalar& maxPolymerAdsorption = elemCtx.problem().maxPolymerAdsorption(elemCtx, dofIdx, timeIdx);
879 polymerAdsorption_ = std::max(Evaluation(maxPolymerAdsorption) , polymerAdsorption_);
883 const Scalar& residualResistanceFactor = PolymerModule::plyrockResidualResistanceFactor(elemCtx, dofIdx, timeIdx);
884 const Evaluation resistanceFactor = 1.0 + (residualResistanceFactor - 1.0) * polymerAdsorption_ / maxAdsorbtion;
887 auto& fs = asImp_().fluidState_;
888 const Evaluation& muWater = fs.viscosity(waterPhaseIdx);
889 const auto& viscosityMultiplier = PolymerModule::plyviscViscosityMultiplierTable(elemCtx, dofIdx, timeIdx);
890 Evaluation viscosityMixture = viscosityMultiplier.eval(polymerConcentration_,
true) * muWater;
893 const Scalar plymixparToddLongstaff = PolymerModule::plymixparToddLongstaff(elemCtx, dofIdx, timeIdx);
894 Evaluation viscosityPolymer = viscosityMultiplier.eval(cmax,
true) * muWater;
895 Evaluation viscosityPolymerEffective = pow(viscosityMixture, plymixparToddLongstaff) * pow(viscosityPolymer, 1.0 - plymixparToddLongstaff);
896 Evaluation viscosityWaterEffective = pow(viscosityMixture, plymixparToddLongstaff) * pow(muWater, 1.0 - plymixparToddLongstaff);
898 Evaluation cbar = polymerConcentration_ / cmax;
900 waterViscosityCorrection_ = muWater * ( (1.0 - cbar) / viscosityWaterEffective + cbar / viscosityPolymerEffective );
903 asImp_().mobility_[waterPhaseIdx] *= waterViscosityCorrection_ / resistanceFactor;
906 polymerViscosityCorrection_ = (muWater / waterViscosityCorrection_) / viscosityPolymerEffective;
909 polymerDeadPoreVolume_ = PolymerModule::plyrockDeadPoreVolume(elemCtx, dofIdx, timeIdx);
910 polymerRockDensity_ = PolymerModule::plyrockRockDensityFactor(elemCtx, dofIdx, timeIdx);
913 const Evaluation& polymerConcentration()
const 914 {
return polymerConcentration_; }
916 const Scalar& polymerDeadPoreVolume()
const 917 {
return polymerDeadPoreVolume_; }
919 const Evaluation& polymerAdsorption()
const 920 {
return polymerAdsorption_; }
922 const Scalar& polymerRockDensity()
const 923 {
return polymerRockDensity_; }
926 const Evaluation& polymerViscosityCorrection()
const 927 {
return polymerViscosityCorrection_; }
930 const Evaluation& waterViscosityCorrection()
const 931 {
return waterViscosityCorrection_; }
935 Implementation& asImp_()
936 {
return *
static_cast<Implementation*
>(
this); }
938 Evaluation polymerConcentration_;
939 Scalar polymerDeadPoreVolume_;
940 Scalar polymerRockDensity_;
941 Evaluation polymerAdsorption_;
942 Evaluation polymerViscosityCorrection_;
943 Evaluation waterViscosityCorrection_;
948 template <
class TypeTag>
951 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
952 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
953 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
957 unsigned scvIdx OPM_UNUSED,
958 unsigned timeIdx OPM_UNUSED)
961 const Evaluation& polymerConcentration()
const 962 { OPM_THROW(std::runtime_error,
"polymerConcentration() called but polymers are disabled"); }
964 const Evaluation& polymerDeadPoreVolume()
const 965 { OPM_THROW(std::runtime_error,
"polymerDeadPoreVolume() called but polymers are disabled"); }
967 const Evaluation& polymerAdsorption()
const 968 { OPM_THROW(std::runtime_error,
"polymerAdsorption() called but polymers are disabled"); }
970 const Evaluation& polymerRockDensity()
const 971 { OPM_THROW(std::runtime_error,
"polymerRockDensity() called but polymers are disabled"); }
973 const Evaluation& polymerViscosityCorrection()
const 974 { OPM_THROW(std::runtime_error,
"polymerViscosityCorrection() called but polymers are disabled"); }
976 const Evaluation& waterViscosityCorrection()
const 977 { OPM_THROW(std::runtime_error,
"waterViscosityCorrection() called but polymers are disabled"); }
990 template <
class TypeTag,
bool enablePolymerV = GET_PROP_VALUE(TypeTag, EnablePolymer)>
993 typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
995 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
996 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
997 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
998 typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
999 typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
1000 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
1001 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
1003 typedef Opm::MathToolbox<Evaluation> Toolbox;
1008 static constexpr
unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
1009 static constexpr
int dimWorld = GridView::dimensionworld;
1010 static constexpr
unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
1012 typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
1013 typedef Dune::FieldVector<Evaluation, dimWorld> DimEvalVector;
1022 template <
class Dummy =
bool>
1025 unsigned scvfIdx OPM_UNUSED,
1026 unsigned timeIdx OPM_UNUSED)
1028 OPM_THROW(std::runtime_error,
1029 "The extension of the blackoil model for polymers is not yet " 1030 "implemented for problems specified using permeabilities.");
1039 template <
class Dummy =
bool>
1046 waterShearFactor_ = 1.0;
1047 polymerShearFactor_ = 1.0;
1049 if (!PolymerModule::hasPlyshlog())
1052 const ExtensiveQuantities& extQuants = asImp_();
1053 unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
1054 unsigned interiorDofIdx = extQuants.interiorIndex();
1055 unsigned exteriorDofIdx = extQuants.exteriorIndex();
1056 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
1057 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
1058 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
1061 Evaluation poroAvg = intQuantsIn.porosity()*0.5 + intQuantsEx.porosity()*0.5;
1062 unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvfIdx, timeIdx);
1063 const Evaluation& Sw = up.fluidState().saturation(waterPhaseIdx);
1064 unsigned cellIdx = elemCtx.globalSpaceIndex(scvfIdx, timeIdx);
1065 const auto& materialLawManager = elemCtx.problem().materialLawManager();
1066 const auto& scaledDrainageInfo =
1067 materialLawManager->oilWaterScaledEpsInfoDrainage(cellIdx);
1068 const Scalar& Swcr = scaledDrainageInfo.Swcr;
1071 Evaluation denom = Opm::max(poroAvg * (Sw - Swcr), 1e-12);
1072 Evaluation waterVolumeVelocity = extQuants.volumeFlux(waterPhaseIdx) / denom;
1075 if (PolymerModule::hasShrate()) {
1076 const Evaluation& relWater = up.relativePermeability(waterPhaseIdx);
1077 Scalar trans = elemCtx.problem().transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
1079 Scalar faceArea = elemCtx.stencil(timeIdx).interiorFace(scvfIdx).area();
1080 auto dist = elemCtx.pos(interiorDofIdx, timeIdx ) - elemCtx.pos(exteriorDofIdx, timeIdx );
1082 Scalar absPerm = trans / faceArea * dist.two_norm();
1083 waterVolumeVelocity *= PolymerModule::shrate(pvtnumRegionIdx) * Opm::sqrt( poroAvg * Sw / (relWater * absPerm ) ) ;
1084 assert(Opm::isfinite(waterVolumeVelocity));
1094 const Evaluation& polymerShearFactor()
const 1095 {
return polymerShearFactor_; }
1097 const Evaluation& waterShearFactor()
const 1098 {
return waterShearFactor_; }
1102 Implementation& asImp_()
1103 {
return *
static_cast<Implementation*
>(
this); }
1105 Evaluation polymerShearFactor_;
1106 Evaluation waterShearFactor_;
1110 template <
class TypeTag>
1113 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
1114 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
1118 unsigned scvfIdx OPM_UNUSED,
1119 unsigned timeIdx OPM_UNUSED)
1123 unsigned scvfIdx OPM_UNUSED,
1124 unsigned timeIdx OPM_UNUSED)
1127 const Evaluation& polymerShearFactor()
const 1128 { OPM_THROW(std::runtime_error,
"polymerShearFactor() called but polymers are disabled"); }
1130 const Evaluation& waterShearFactor()
const 1131 { OPM_THROW(std::runtime_error,
"waterShearFactor() called but polymers are disabled"); }
Provides the volumetric quantities required for the equations needed by the polymers extension of the...
Definition: blackoilpolymermodules.hh:839
static Scalar computeResidualError(const EqVector &resid)
Return how much a residual is considered an error.
Definition: blackoilpolymermodules.hh:549
static void setAdsrock(unsigned satRegionIdx, const TabulatedFunction &plyadsAdsorbedPolymer)
Specify the polymer rock properties a single region.
Definition: blackoilpolymermodules.hh:300
Definition: baseauxiliarymodule.hh:37
VTK output module for the black oil model's polymer related quantities.
static Scalar computeUpdateError(const PrimaryVariables &oldPv OPM_UNUSED, const EqVector &delta OPM_UNUSED)
Return how much a Newton-Raphson update is considered an error.
Definition: blackoilpolymermodules.hh:537
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
void updateShearMultipliersPerm(const ElementContext &elemCtx OPM_UNUSED, unsigned scvfIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
Method which calculates the shear factor based on flow velocity.
Definition: blackoilpolymermodules.hh:1024
static void assignPrimaryVars(PrimaryVariables &priVars, Scalar polymerConcentration)
Assign the polymer specific primary variables to a PrimaryVariables object.
Definition: blackoilpolymermodules.hh:511
Provides the polymer specific extensive quantities to the generic black-oil module's extensive quanti...
Definition: blackoilpolymermodules.hh:991
static void setPlmixpar(unsigned mixRegionIdx, const Scalar &plymixparToddLongstaff)
Specify the maximum polymer concentration a single region.
Definition: blackoilpolymermodules.hh:354
void polymerPropertiesUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the intensive properties needed to handle polymers from the primary variables.
Definition: blackoilpolymermodules.hh:865
static void registerOutputModules(Model &model, Simulator &simulator)
Register all polymer specific VTK and ECL output modules.
Definition: blackoilpolymermodules.hh:377
Declares the properties required by the black oil model.
static void setNumSatRegions(unsigned numRegions)
Specify the number of satuation regions.
Definition: blackoilpolymermodules.hh:266
static Evaluation computeShearFactor(const Evaluation &polymerConcentration, unsigned pvtnumRegionIdx, const Evaluation &v0)
Computes the shear factor.
Definition: blackoilpolymermodules.hh:681
void updateShearMultipliers(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Method which calculates the shear factor based on flow velocity.
Definition: blackoilpolymermodules.hh:1041
static void updatePrimaryVars(PrimaryVariables &newPv, const PrimaryVariables &oldPv, const EqVector &delta)
Do a Newton-Raphson update the primary variables of the polymers.
Definition: blackoilpolymermodules.hh:523
This method contains all callback classes for quantities that are required by some extensive quantiti...
static void setNumPvtRegions(unsigned numRegions)
Specify the number of pvt regions.
Definition: blackoilpolymermodules.hh:311
static void setPlyvisc(unsigned satRegionIdx, const TabulatedFunction &plyviscViscosityMultiplierTable)
Specify the polymer viscosity a single region.
Definition: blackoilpolymermodules.hh:321
static void registerParameters()
Register all run-time parameters for the black-oil polymer module.
Definition: blackoilpolymermodules.hh:365
static void setPlyrock(unsigned satRegionIdx, const Scalar &plyrockDeadPoreVolume, const Scalar &plyrockResidualResistanceFactor, const Scalar &plyrockRockDensityFactor, const Scalar &plyrockAdsorbtionIndex, const Scalar &plyrockMaxAdsorbtion)
Specify the polymer rock properties a single region.
Definition: blackoilpolymermodules.hh:281
static void setNumMixRegions(unsigned numRegions)
Specify the number of mix regions.
Definition: blackoilpolymermodules.hh:332
Contains the high level supplements required to extend the black oil model by polymer.
Definition: blackoilpolymermodules.hh:64
static void registerParameters()
Register all run-time parameters for the multi-phase VTK output module.
Definition: vtkblackoilpolymermodule.hh:101
static void setPlymax(unsigned mixRegionIdx, const Scalar &plymaxMaxConcentration)
Specify the maximum polymer concentration a single region.
Definition: blackoilpolymermodules.hh:343
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
VTK output module for the black oil model's polymer related quantities.
Definition: vtkblackoilpolymermodule.hh:75