28#ifndef EWOMS_BLACK_OIL_POLYMER_MODULE_HH
29#define EWOMS_BLACK_OIL_POLYMER_MODULE_HH
36#include <opm/common/OpmLog/OpmLog.hpp>
39#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
40#include <opm/input/eclipse/EclipseState/Tables/PlyadsTable.hpp>
41#include <opm/input/eclipse/EclipseState/Tables/PlymaxTable.hpp>
42#include <opm/input/eclipse/EclipseState/Tables/PlyrockTable.hpp>
43#include <opm/input/eclipse/EclipseState/Tables/PlyshlogTable.hpp>
44#include <opm/input/eclipse/EclipseState/Tables/PlyviscTable.hpp>
47#include <dune/common/fvector.hh>
59template <class TypeTag, bool enablePolymerV = getPropValue<TypeTag, Properties::EnablePolymer>()>
75 using Toolbox = MathToolbox<Evaluation>;
80 static constexpr unsigned polymerConcentrationIdx = Indices::polymerConcentrationIdx;
81 static constexpr unsigned polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
82 static constexpr unsigned contiPolymerEqIdx = Indices::contiPolymerEqIdx;
83 static constexpr unsigned contiPolymerMolarWeightEqIdx = Indices::contiPolymerMWEqIdx;
84 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
87 static constexpr unsigned enablePolymer = enablePolymerV;
88 static constexpr bool enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>();
90 static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
91 static constexpr unsigned numPhases = FluidSystem::numPhases;
98 static void initFromState(
const EclipseState& eclState)
102 if (enablePolymer && !eclState.runspec().phases().active(Phase::POLYMER)) {
103 throw std::runtime_error(
"Non-trivial polymer treatment requested at compile time, but "
104 "the deck does not contain the POLYMER keyword");
106 else if (!enablePolymer && eclState.runspec().phases().active(Phase::POLYMER)) {
107 throw std::runtime_error(
"Polymer treatment disabled at compile time, but the deck "
108 "contains the POLYMER keyword");
111 if (enablePolymerMolarWeight && !eclState.runspec().phases().active(Phase::POLYMW)) {
112 throw std::runtime_error(
"Polymer molecular weight tracking is enabled at compile time, but "
113 "the deck does not contain the POLYMW keyword");
115 else if (!enablePolymerMolarWeight && eclState.runspec().phases().active(Phase::POLYMW)) {
116 throw std::runtime_error(
"Polymer molecular weight tracking is disabled at compile time, but the deck "
117 "contains the POLYMW keyword");
120 if (enablePolymerMolarWeight && !enablePolymer) {
121 throw std::runtime_error(
"Polymer molecular weight tracking is enabled while polymer treatment "
122 "is disabled at compile time");
125 if (!eclState.runspec().phases().active(Phase::POLYMER))
128 const auto& tableManager = eclState.getTableManager();
130 unsigned numSatRegions = tableManager.getTabdims().getNumSatTables();
131 params_.setNumSatRegions(numSatRegions);
134 const auto& plyrockTables = tableManager.getPlyrockTables();
135 if (!plyrockTables.empty()) {
136 assert(numSatRegions == plyrockTables.size());
137 for (
unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++ satRegionIdx) {
138 const auto& plyrockTable = plyrockTables.template getTable<PlyrockTable>(satRegionIdx);
139 params_.setPlyrock(satRegionIdx,
140 plyrockTable.getDeadPoreVolumeColumn()[0],
141 plyrockTable.getResidualResistanceFactorColumn()[0],
142 plyrockTable.getRockDensityFactorColumn()[0],
144 plyrockTable.getMaxAdsorbtionColumn()[0]);
148 throw std::runtime_error(
"PLYROCK must be specified in POLYMER runs\n");
152 const auto& plyadsTables = tableManager.getPlyadsTables();
153 if (!plyadsTables.empty()) {
154 assert(numSatRegions == plyadsTables.size());
155 for (
unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++ satRegionIdx) {
156 const auto& plyadsTable = plyadsTables.template getTable<PlyadsTable>(satRegionIdx);
158 const auto& c = plyadsTable.getPolymerConcentrationColumn();
159 const auto& ads = plyadsTable.getAdsorbedPolymerColumn();
160 params_.plyadsAdsorbedPolymer_[satRegionIdx].setXYContainers(c, ads);
164 throw std::runtime_error(
"PLYADS must be specified in POLYMER runs\n");
168 unsigned numPvtRegions = tableManager.getTabdims().getNumPVTTables();
169 params_.plyviscViscosityMultiplierTable_.resize(numPvtRegions);
172 const auto& plyviscTables = tableManager.getPlyviscTables();
173 if (!plyviscTables.empty()) {
175 if (enablePolymerMolarWeight) {
176 OpmLog::warning(
"PLYVISC should not be used in POLYMW runs, "
177 "it will have no effect. A viscosity model based on PLYVMH is used instead.\n");
180 assert(numPvtRegions == plyviscTables.size());
181 for (
unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) {
182 const auto& plyadsTable = plyviscTables.template getTable<PlyviscTable>(pvtRegionIdx);
184 const auto& c = plyadsTable.getPolymerConcentrationColumn();
185 const auto& visc = plyadsTable.getViscosityMultiplierColumn();
186 params_.plyviscViscosityMultiplierTable_[pvtRegionIdx].setXYContainers(c, visc);
190 else if (!enablePolymerMolarWeight) {
191 throw std::runtime_error(
"PLYVISC must be specified in POLYMER runs\n");
195 const auto& plymaxTables = tableManager.getPlymaxTables();
196 const unsigned numMixRegions = plymaxTables.size();
197 params_.setNumMixRegions(numMixRegions, enablePolymerMolarWeight);
198 if (!plymaxTables.empty()) {
199 for (
unsigned mixRegionIdx = 0; mixRegionIdx < numMixRegions; ++ mixRegionIdx) {
200 const auto& plymaxTable = plymaxTables.template getTable<PlymaxTable>(mixRegionIdx);
201 params_.plymaxMaxConcentration_[mixRegionIdx] = plymaxTable.getPolymerConcentrationColumn()[0];
205 throw std::runtime_error(
"PLYMAX must be specified in POLYMER runs\n");
208 if (!eclState.getTableManager().getPlmixparTable().empty()) {
209 if (enablePolymerMolarWeight) {
210 OpmLog::warning(
"PLMIXPAR should not be used in POLYMW runs, it will have no effect.\n");
213 const auto& plmixparTable = eclState.getTableManager().getPlmixparTable();
215 for (
unsigned mixRegionIdx = 0; mixRegionIdx < numMixRegions; ++ mixRegionIdx) {
216 params_.plymixparToddLongstaff_[mixRegionIdx] = plmixparTable[mixRegionIdx].todd_langstaff;
220 else if (!enablePolymerMolarWeight) {
221 throw std::runtime_error(
"PLMIXPAR must be specified in POLYMER runs\n");
224 params_.hasPlyshlog_ = eclState.getTableManager().hasTables(
"PLYSHLOG");
225 params_.hasShrate_ = eclState.getTableManager().useShrate();
227 if ((params_.hasPlyshlog_ || params_.hasShrate_) && enablePolymerMolarWeight) {
228 OpmLog::warning(
"PLYSHLOG and SHRATE should not be used in POLYMW runs, they will have no effect.\n");
231 if (params_.hasPlyshlog_ && !enablePolymerMolarWeight) {
232 const auto& plyshlogTables = tableManager.getPlyshlogTables();
233 assert(numPvtRegions == plyshlogTables.size());
234 params_.plyshlogShearEffectRefMultiplier_.resize(numPvtRegions);
235 params_.plyshlogShearEffectRefLogVelocity_.resize(numPvtRegions);
236 for (
unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) {
237 const auto& plyshlogTable = plyshlogTables.template getTable<PlyshlogTable>(pvtRegionIdx);
239 Scalar plyshlogRefPolymerConcentration = plyshlogTable.getRefPolymerConcentration();
240 auto waterVelocity = plyshlogTable.getWaterVelocityColumn().vectorCopy();
241 auto shearMultiplier = plyshlogTable.getShearMultiplierColumn().vectorCopy();
244 UnitSystem unitSystem = eclState.getDeckUnitSystem();
245 double siFactor = params_.hasShrate_? unitSystem.parse(
"1/Time").getSIScaling() : unitSystem.parse(
"Length/Time").getSIScaling();
246 for (
size_t i = 0; i < waterVelocity.size(); ++i) {
247 waterVelocity[i] *= siFactor;
250 waterVelocity[i] = std::log(waterVelocity[i]);
253 Scalar refViscMult = params_.plyviscViscosityMultiplierTable_[pvtRegionIdx].eval(plyshlogRefPolymerConcentration,
true);
255 for (
size_t i = 0; i < waterVelocity.size(); ++i) {
256 shearMultiplier[i] *= refViscMult;
257 shearMultiplier[i] -= 1;
258 shearMultiplier[i] /= (refViscMult - 1);
259 shearMultiplier[i] = shearMultiplier[i];
261 params_.plyshlogShearEffectRefMultiplier_[pvtRegionIdx].resize(waterVelocity.size());
262 params_.plyshlogShearEffectRefLogVelocity_[pvtRegionIdx].resize(waterVelocity.size());
264 for (
size_t i = 0; i < waterVelocity.size(); ++i) {
265 params_.plyshlogShearEffectRefMultiplier_[pvtRegionIdx][i] = shearMultiplier[i];
266 params_.plyshlogShearEffectRefLogVelocity_[pvtRegionIdx][i] = waterVelocity[i];
271 if (params_.hasShrate_ && !enablePolymerMolarWeight) {
272 if (!params_.hasPlyshlog_) {
273 throw std::runtime_error(
"PLYSHLOG must be specified if SHRATE is used in POLYMER runs\n");
275 const auto& shrateTable = eclState.getTableManager().getShrateTable();
276 params_.shrate_.resize(numPvtRegions);
277 for (
unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) {
278 if (shrateTable.empty()) {
279 params_.shrate_[pvtRegionIdx] = 4.8;
281 else if (shrateTable.size() == numPvtRegions) {
282 params_.shrate_[pvtRegionIdx] = shrateTable[pvtRegionIdx].rate;
285 throw std::runtime_error(
"SHRATE must either have 0 or number of NUMPVT entries\n");
290 if constexpr (enablePolymerMolarWeight) {
291 const auto& plyvmhTable = eclState.getTableManager().getPlyvmhTable();
292 if (!plyvmhTable.empty()) {
293 assert(plyvmhTable.size() == numMixRegions);
294 for (
size_t regionIdx = 0; regionIdx < numMixRegions; ++regionIdx) {
295 params_.plyvmhCoefficients_[regionIdx].k_mh = plyvmhTable[regionIdx].k_mh;
296 params_.plyvmhCoefficients_[regionIdx].a_mh = plyvmhTable[regionIdx].a_mh;
297 params_.plyvmhCoefficients_[regionIdx].gamma = plyvmhTable[regionIdx].gamma;
298 params_.plyvmhCoefficients_[regionIdx].kappa = plyvmhTable[regionIdx].kappa;
302 throw std::runtime_error(
"PLYVMH keyword must be specified in POLYMW rus \n");
306 const auto& plymwinjTables = tableManager.getPlymwinjTables();
307 for (
const auto& table : plymwinjTables) {
308 const int tableNumber = table.first;
309 const auto& plymwinjtable = table.second;
310 const std::vector<double>& throughput = plymwinjtable.getThroughputs();
311 const std::vector<double>& watervelocity = plymwinjtable.getVelocities();
312 const std::vector<std::vector<double>>& molecularweight = plymwinjtable.getMoleWeights();
313 TabulatedTwoDFunction tablefunc(throughput, watervelocity, molecularweight,
true,
false);
314 params_.plymwinjTables_[tableNumber] = std::move(tablefunc);
318 const auto& skprwatTables = tableManager.getSkprwatTables();
319 for (
const auto& table : skprwatTables) {
320 const int tableNumber = table.first;
321 const auto& skprwattable = table.second;
322 const std::vector<double>& throughput = skprwattable.getThroughputs();
323 const std::vector<double>& watervelocity = skprwattable.getVelocities();
324 const std::vector<std::vector<double>>& skinpressure = skprwattable.getSkinPressures();
325 TabulatedTwoDFunction tablefunc(throughput, watervelocity, skinpressure,
true,
false);
326 params_.skprwatTables_[tableNumber] = std::move(tablefunc);
330 const auto& skprpolyTables = tableManager.getSkprpolyTables();
331 for (
const auto& table : skprpolyTables) {
332 const int tableNumber = table.first;
333 const auto& skprpolytable = table.second;
334 const std::vector<double>& throughput = skprpolytable.getThroughputs();
335 const std::vector<double>& watervelocity = skprpolytable.getVelocities();
336 const std::vector<std::vector<double>>& skinpressure = skprpolytable.getSkinPressures();
337 const double refPolymerConcentration = skprpolytable.referenceConcentration();
339 {refPolymerConcentration,
340 TabulatedTwoDFunction(throughput, watervelocity, skinpressure,
true,
false)};
341 params_.skprpolyTables_[tableNumber] = std::move(tablefunc);
352 const auto iterTable = params_.plymwinjTables_.find(tableNumber);
353 if (iterTable != params_.plymwinjTables_.end()) {
354 return iterTable->second;
357 throw std::runtime_error(
" the PLYMWINJ table " + std::to_string(tableNumber) +
" does not exist\n");
366 const auto iterTable = params_.skprwatTables_.find(tableNumber);
367 if (iterTable != params_.skprwatTables_.end()) {
368 return iterTable->second;
371 throw std::runtime_error(
" the SKPRWAT table " + std::to_string(tableNumber) +
" does not exist\n");
381 const auto iterTable = params_.skprpolyTables_.find(tableNumber);
382 if (iterTable != params_.skprpolyTables_.end()) {
383 return iterTable->second;
386 throw std::runtime_error(
" the SKPRPOLY table " + std::to_string(tableNumber) +
" does not exist\n");
395 if constexpr (enablePolymer)
403 Simulator& simulator)
405 if constexpr (enablePolymer)
411 if constexpr (enablePolymer) {
412 if constexpr (enablePolymerMolarWeight)
413 return pvIdx == polymerConcentrationIdx || pvIdx == polymerMoleWeightIdx;
415 return pvIdx == polymerConcentrationIdx;
425 if (pvIdx == polymerConcentrationIdx) {
426 return "polymer_waterconcentration";
429 return "polymer_molecularweight";
438 return static_cast<Scalar
>(1.0);
443 if constexpr (enablePolymer) {
444 if constexpr (enablePolymerMolarWeight)
445 return eqIdx == contiPolymerEqIdx || eqIdx == contiPolymerMolarWeightEqIdx;
447 return eqIdx == contiPolymerEqIdx;
453 static std::string
eqName(
unsigned eqIdx)
457 if (eqIdx == contiPolymerEqIdx)
458 return "conti^polymer";
460 return "conti^polymer_molecularweight";
463 static Scalar
eqWeight([[maybe_unused]]
unsigned eqIdx)
468 return static_cast<Scalar
>(1.0);
472 template <
class LhsEval>
473 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
474 const IntensiveQuantities& intQuants)
476 if constexpr (enablePolymer) {
477 const auto& fs = intQuants.fluidState();
479 LhsEval surfaceVolumeWater =
480 Toolbox::template decay<LhsEval>(fs.saturation(waterPhaseIdx))
481 * Toolbox::template decay<LhsEval>(fs.invB(waterPhaseIdx))
482 * Toolbox::template decay<LhsEval>(intQuants.porosity());
485 surfaceVolumeWater = max(surfaceVolumeWater, 1e-10);
488 const LhsEval massPolymer = surfaceVolumeWater
489 * Toolbox::template decay<LhsEval>(intQuants.polymerConcentration())
490 * (1.0 - Toolbox::template decay<LhsEval>(intQuants.polymerDeadPoreVolume()));
493 const LhsEval adsorptionPolymer =
494 Toolbox::template decay<LhsEval>(1.0 - intQuants.porosity())
495 * Toolbox::template decay<LhsEval>(intQuants.polymerRockDensity())
496 * Toolbox::template decay<LhsEval>(intQuants.polymerAdsorption());
498 LhsEval accumulationPolymer = massPolymer + adsorptionPolymer;
500 storage[contiPolymerEqIdx] += accumulationPolymer;
503 if constexpr (enablePolymerMolarWeight) {
504 accumulationPolymer = max(accumulationPolymer, 1e-10);
506 storage[contiPolymerMolarWeightEqIdx] += accumulationPolymer
507 * Toolbox::template decay<LhsEval> (intQuants.polymerMoleWeight());
513 [[maybe_unused]]
const ElementContext& elemCtx,
514 [[maybe_unused]]
unsigned scvfIdx,
515 [[maybe_unused]]
unsigned timeIdx)
518 if constexpr (enablePolymer) {
519 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
521 const unsigned upIdx = extQuants.upstreamIndex(FluidSystem::waterPhaseIdx);
522 const unsigned inIdx = extQuants.interiorIndex();
523 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
524 const unsigned contiWaterEqIdx = Indices::conti0EqIdx + Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
526 if (upIdx == inIdx) {
527 flux[contiPolymerEqIdx] =
528 extQuants.volumeFlux(waterPhaseIdx)
529 *up.fluidState().invB(waterPhaseIdx)
530 *up.polymerViscosityCorrection()
531 /extQuants.polymerShearFactor()
532 *up.polymerConcentration();
535 flux[contiWaterEqIdx] /=
536 extQuants.waterShearFactor();
539 flux[contiPolymerEqIdx] =
540 extQuants.volumeFlux(waterPhaseIdx)
541 *decay<Scalar>(up.fluidState().invB(waterPhaseIdx))
542 *decay<Scalar>(up.polymerViscosityCorrection())
543 /decay<Scalar>(extQuants.polymerShearFactor())
544 *decay<Scalar>(up.polymerConcentration());
547 flux[contiWaterEqIdx] /=
548 decay<Scalar>(extQuants.waterShearFactor());
552 if constexpr (enablePolymerMolarWeight) {
554 flux[contiPolymerMolarWeightEqIdx] =
555 flux[contiPolymerEqIdx]*up.polymerMoleWeight();
557 flux[contiPolymerMolarWeightEqIdx] =
558 flux[contiPolymerEqIdx]*decay<Scalar>(up.polymerMoleWeight());
572 return static_cast<Scalar
>(0.0);
575 template <
class DofEntity>
576 static void serializeEntity(
const Model& model, std::ostream& outstream,
const DofEntity& dof)
578 if constexpr (enablePolymer) {
579 unsigned dofIdx = model.dofMapper().index(dof);
580 const PrimaryVariables& priVars = model.solution(0)[dofIdx];
581 outstream << priVars[polymerConcentrationIdx];
582 outstream << priVars[polymerMoleWeightIdx];
586 template <
class DofEntity>
589 if constexpr (enablePolymer) {
590 unsigned dofIdx = model.dofMapper().index(dof);
591 PrimaryVariables& priVars0 = model.solution(0)[dofIdx];
592 PrimaryVariables& priVars1 = model.solution(1)[dofIdx];
594 instream >> priVars0[polymerConcentrationIdx];
595 instream >> priVars0[polymerMoleWeightIdx];
598 priVars1[polymerConcentrationIdx] = priVars0[polymerConcentrationIdx];
599 priVars1[polymerMoleWeightIdx] = priVars0[polymerMoleWeightIdx];
607 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
608 return params_.plyrockDeadPoreVolume_[satnumRegionIdx];
615 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
616 return params_.plyrockResidualResistanceFactor_[satnumRegionIdx];
623 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
624 return params_.plyrockRockDensityFactor_[satnumRegionIdx];
631 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
632 return params_.plyrockAdsorbtionIndex_[satnumRegionIdx];
639 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
640 return params_.plyrockMaxAdsorbtion_[satnumRegionIdx];
647 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
648 return params_.plyadsAdsorbedPolymer_[satnumRegionIdx];
655 unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
656 return params_.plyviscViscosityMultiplierTable_[pvtnumRegionIdx];
661 return params_.plyviscViscosityMultiplierTable_[pvtnumRegionIdx];
668 unsigned polymerMixRegionIdx = elemCtx.problem().plmixnumRegionIndex(elemCtx, scvIdx, timeIdx);
669 return params_.plymaxMaxConcentration_[polymerMixRegionIdx];
676 unsigned polymerMixRegionIdx = elemCtx.problem().plmixnumRegionIndex(elemCtx, scvIdx, timeIdx);
677 return params_.plymixparToddLongstaff_[polymerMixRegionIdx];
682 const unsigned scvIdx,
683 const unsigned timeIdx)
685 const unsigned polymerMixRegionIdx = elemCtx.problem().plmixnumRegionIndex(elemCtx, scvIdx, timeIdx);
686 return params_.plyvmhCoefficients_[polymerMixRegionIdx];
691 return params_.hasPlyshlog_;
696 return params_.hasShrate_;
699 static const Scalar
shrate(
unsigned pvtnumRegionIdx)
701 return params_.shrate_[pvtnumRegionIdx];
710 template <
class Evaluation>
712 unsigned pvtnumRegionIdx,
713 const Evaluation& v0)
715 using ToolboxLocal = MathToolbox<Evaluation>;
717 const auto& viscosityMultiplierTable = params_.plyviscViscosityMultiplierTable_[pvtnumRegionIdx];
718 Scalar viscosityMultiplier = viscosityMultiplierTable.eval(scalarValue(polymerConcentration),
true);
720 const Scalar eps = 1e-14;
722 if (std::abs((viscosityMultiplier - 1.0)) < eps)
723 return ToolboxLocal::createConstant(v0, 1.0);
725 const std::vector<Scalar>& shearEffectRefLogVelocity = params_.plyshlogShearEffectRefLogVelocity_[pvtnumRegionIdx];
726 auto v0AbsLog = log(abs(v0));
728 if (v0AbsLog < shearEffectRefLogVelocity[0])
729 return ToolboxLocal::createConstant(v0, 1.0);
735 const std::vector<Scalar>& shearEffectRefMultiplier = params_.plyshlogShearEffectRefMultiplier_[pvtnumRegionIdx];
736 size_t numTableEntries = shearEffectRefLogVelocity.size();
737 assert(shearEffectRefMultiplier.size() == numTableEntries);
739 std::vector<Scalar> shearEffectMultiplier(numTableEntries, 1.0);
740 for (
size_t i = 0; i < numTableEntries; ++i) {
741 shearEffectMultiplier[i] = (1.0 + (viscosityMultiplier - 1.0)*shearEffectRefMultiplier[i]) / viscosityMultiplier;
742 shearEffectMultiplier[i] = log(shearEffectMultiplier[i]);
746 TabulatedFunction logShearEffectMultiplier = TabulatedFunction(numTableEntries, shearEffectRefLogVelocity, shearEffectMultiplier,
false);
753 auto F = [&logShearEffectMultiplier, &v0AbsLog](
const Evaluation& u) {
754 return u + logShearEffectMultiplier.eval(u,
true) - v0AbsLog;
757 auto dF = [&logShearEffectMultiplier](
const Evaluation& u) {
758 return 1 + logShearEffectMultiplier.evalDerivative(u,
true);
764 bool converged =
false;
766 for (
int i = 0; i < 20; ++i) {
770 if (std::abs(scalarValue(f)) < 1e-12) {
776 throw std::runtime_error(
"Not able to compute shear velocity. \n");
780 return exp(logShearEffectMultiplier.eval(u,
true));
794template <
class TypeTag,
bool enablePolymerV>
795BlackOilPolymerParams<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
796BlackOilPolymerModule<TypeTag, enablePolymerV>::params_;
805template <class TypeTag, bool enablePolymerV = getPropValue<TypeTag, Properties::EnablePolymer>()>
820 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
821 static constexpr int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
822 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
823 static constexpr bool enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>();
824 static constexpr int polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
838 const auto linearizationType = elemCtx.linearizationType();
839 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
841 if constexpr (enablePolymerMolarWeight) {
842 polymerMoleWeight_ = priVars.makeEvaluation(polymerMoleWeightIdx, timeIdx, linearizationType);
850 const Scalar& maxPolymerAdsorption = elemCtx.problem().maxPolymerAdsorption(elemCtx, dofIdx, timeIdx);
856 const Evaluation resistanceFactor = 1.0 + (residualResistanceFactor - 1.0) *
polymerAdsorption_ / maxAdsorbtion;
859 if constexpr (!enablePolymerMolarWeight) {
861 const auto& fs =
asImp_().fluidState_;
862 const Evaluation& muWater = fs.viscosity(waterPhaseIdx);
868 const Evaluation viscosityPolymer = viscosityMultiplier.eval(cmax,
true) * muWater;
869 const Evaluation viscosityPolymerEffective = pow(viscosityMixture, plymixparToddLongstaff) * pow(viscosityPolymer, 1.0 - plymixparToddLongstaff);
870 const Evaluation viscosityWaterEffective = pow(viscosityMixture, plymixparToddLongstaff) * pow(muWater, 1.0 - plymixparToddLongstaff);
880 const Scalar k_mh = plyvmhCoefficients.k_mh;
881 const Scalar a_mh = plyvmhCoefficients.a_mh;
882 const Scalar gamma = plyvmhCoefficients.gamma;
883 const Scalar kappa = plyvmhCoefficients.kappa;
906 if constexpr (enablePolymerMolarWeight)
909 throw std::logic_error(
"polymerMoleWeight() is called but polymer milecular weight is disabled");
932 {
return *
static_cast<Implementation*
>(
this); }
946template <
class TypeTag>
960 {
throw std::logic_error(
"polymerMoleWeight() called but polymer molecular weight is disabled"); }
963 {
throw std::runtime_error(
"polymerConcentration() called but polymers are disabled"); }
966 {
throw std::runtime_error(
"polymerDeadPoreVolume() called but polymers are disabled"); }
969 {
throw std::runtime_error(
"polymerAdsorption() called but polymers are disabled"); }
972 {
throw std::runtime_error(
"polymerRockDensity() called but polymers are disabled"); }
975 {
throw std::runtime_error(
"polymerViscosityCorrection() called but polymers are disabled"); }
978 {
throw std::runtime_error(
"waterViscosityCorrection() called but polymers are disabled"); }
989template <class TypeTag, bool enablePolymerV = getPropValue<TypeTag, Properties::EnablePolymer>()>
1002 static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
1003 static constexpr int dimWorld = GridView::dimensionworld;
1004 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
1006 using Toolbox = MathToolbox<Evaluation>;
1008 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
1009 using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
1018 template <
class Dummy =
bool>
1024 throw std::runtime_error(
"The extension of the blackoil model for polymers is not yet "
1025 "implemented for problems specified using permeabilities.");
1034 template <
class Dummy =
bool>
1041 waterShearFactor_ = 1.0;
1042 polymerShearFactor_ = 1.0;
1047 const ExtensiveQuantities& extQuants = asImp_();
1048 unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
1049 unsigned interiorDofIdx = extQuants.interiorIndex();
1050 unsigned exteriorDofIdx = extQuants.exteriorIndex();
1051 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
1052 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
1053 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
1056 Evaluation poroAvg = intQuantsIn.porosity()*0.5 + intQuantsEx.porosity()*0.5;
1057 unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvfIdx, timeIdx);
1058 const Evaluation& Sw = up.fluidState().saturation(waterPhaseIdx);
1059 unsigned cellIdx = elemCtx.globalSpaceIndex(scvfIdx, timeIdx);
1060 const auto& materialLawManager = elemCtx.problem().materialLawManager();
1061 const auto& scaledDrainageInfo =
1062 materialLawManager->oilWaterScaledEpsInfoDrainage(cellIdx);
1063 const Scalar& Swcr = scaledDrainageInfo.Swcr;
1066 Evaluation denom = max(poroAvg * (Sw - Swcr), 1e-12);
1067 Evaluation waterVolumeVelocity = extQuants.volumeFlux(waterPhaseIdx) / denom;
1071 const Evaluation& relWater = up.relativePermeability(waterPhaseIdx);
1072 Scalar trans = elemCtx.problem().transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
1074 Scalar faceArea = elemCtx.stencil(timeIdx).interiorFace(scvfIdx).area();
1075 auto dist = elemCtx.pos(interiorDofIdx, timeIdx) - elemCtx.pos(exteriorDofIdx, timeIdx);
1077 Scalar absPerm = trans / faceArea * dist.two_norm();
1078 waterVolumeVelocity *=
1080 assert(isfinite(waterVolumeVelocity));
1088 waterVolumeVelocity);
1089 polymerShearFactor_ =
1092 waterVolumeVelocity*up.polymerViscosityCorrection());
1097 {
return polymerShearFactor_; }
1100 {
return waterShearFactor_; }
1104 Implementation& asImp_()
1105 {
return *
static_cast<Implementation*
>(
this); }
1107 Evaluation polymerShearFactor_;
1108 Evaluation waterShearFactor_;
1112template <
class TypeTag>
1130 {
throw std::runtime_error(
"polymerShearFactor() called but polymers are disabled"); }
1133 {
throw std::runtime_error(
"waterShearFactor() called but polymers are disabled"); }
Contains the parameters required to extend the black-oil model by polymer.
Declares the properties required by the black oil model.
const Evaluation & polymerShearFactor() const
Definition: blackoilpolymermodules.hh:1129
void updateShearMultipliersPerm(const ElementContext &, unsigned, unsigned)
Definition: blackoilpolymermodules.hh:1124
void updateShearMultipliers(const ElementContext &, unsigned, unsigned)
Definition: blackoilpolymermodules.hh:1119
const Evaluation & waterShearFactor() const
Definition: blackoilpolymermodules.hh:1132
Provides the polymer specific extensive quantities to the generic black-oil module's extensive quanti...
Definition: blackoilpolymermodules.hh:991
void updateShearMultipliers(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Method which calculates the shear factor based on flow velocity.
Definition: blackoilpolymermodules.hh:1036
const Evaluation & waterShearFactor() const
Definition: blackoilpolymermodules.hh:1099
void updateShearMultipliersPerm(const ElementContext &, unsigned, unsigned)
Method which calculates the shear factor based on flow velocity.
Definition: blackoilpolymermodules.hh:1020
const Evaluation & polymerShearFactor() const
Definition: blackoilpolymermodules.hh:1096
const Evaluation & polymerConcentration() const
Definition: blackoilpolymermodules.hh:962
const Evaluation & waterViscosityCorrection() const
Definition: blackoilpolymermodules.hh:977
const Evaluation & polymerRockDensity() const
Definition: blackoilpolymermodules.hh:971
const Evaluation & polymerDeadPoreVolume() const
Definition: blackoilpolymermodules.hh:965
const Evaluation & polymerMoleWeight() const
Definition: blackoilpolymermodules.hh:959
void polymerPropertiesUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilpolymermodules.hh:954
const Evaluation & polymerAdsorption() const
Definition: blackoilpolymermodules.hh:968
const Evaluation & polymerViscosityCorrection() const
Definition: blackoilpolymermodules.hh:974
Provides the volumetric quantities required for the equations needed by the polymers extension of the...
Definition: blackoilpolymermodules.hh:807
const Evaluation & polymerAdsorption() const
Definition: blackoilpolymermodules.hh:915
Scalar polymerRockDensity_
Definition: blackoilpolymermodules.hh:938
const Evaluation & polymerMoleWeight() const
Definition: blackoilpolymermodules.hh:904
const Evaluation & polymerViscosityCorrection() const
Definition: blackoilpolymermodules.hh:922
const Scalar & polymerRockDensity() const
Definition: blackoilpolymermodules.hh:918
void polymerPropertiesUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the intensive properties needed to handle polymers from the primary variables.
Definition: blackoilpolymermodules.hh:834
Evaluation waterViscosityCorrection_
Definition: blackoilpolymermodules.hh:941
Scalar polymerDeadPoreVolume_
Definition: blackoilpolymermodules.hh:937
Evaluation polymerAdsorption_
Definition: blackoilpolymermodules.hh:939
Evaluation polymerMoleWeight_
Definition: blackoilpolymermodules.hh:936
const Evaluation & waterViscosityCorrection() const
Definition: blackoilpolymermodules.hh:926
Evaluation polymerConcentration_
Definition: blackoilpolymermodules.hh:934
const Evaluation & polymerConcentration() const
Definition: blackoilpolymermodules.hh:901
Implementation & asImp_()
Definition: blackoilpolymermodules.hh:931
const Scalar & polymerDeadPoreVolume() const
Definition: blackoilpolymermodules.hh:912
Evaluation polymerViscosityCorrection_
Definition: blackoilpolymermodules.hh:940
Contains the high level supplements required to extend the black oil model by polymer.
Definition: blackoilpolymermodules.hh:61
static const Scalar plyrockAdsorbtionIndex(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:627
static const TabulatedFunction & plyviscViscosityMultiplierTable(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:651
static std::string primaryVarName(unsigned pvIdx)
Definition: blackoilpolymermodules.hh:421
static const TabulatedFunction & plyadsAdsorbedPolymer(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:643
const Scalar molarMass() const
Definition: blackoilpolymermodules.hh:784
static TabulatedTwoDFunction & getPlymwinjTable(const int tableNumber)
get the PLYMWINJ table
Definition: blackoilpolymermodules.hh:350
static void serializeEntity(const Model &model, std::ostream &outstream, const DofEntity &dof)
Definition: blackoilpolymermodules.hh:576
static Scalar eqWeight(unsigned eqIdx)
Definition: blackoilpolymermodules.hh:463
static void registerParameters()
Register all run-time parameters for the black-oil polymer module.
Definition: blackoilpolymermodules.hh:393
static bool eqApplies(unsigned eqIdx)
Definition: blackoilpolymermodules.hh:441
static Scalar primaryVarWeight(unsigned pvIdx)
Definition: blackoilpolymermodules.hh:433
static BlackOilPolymerParams< Scalar >::SkprpolyTable & getSkprpolyTable(const int tableNumber)
get the SKPRPOLY table
Definition: blackoilpolymermodules.hh:379
static bool hasShrate()
Definition: blackoilpolymermodules.hh:694
static const TabulatedFunction & plyviscViscosityMultiplierTable(unsigned pvtnumRegionIdx)
Definition: blackoilpolymermodules.hh:659
static const Scalar plyrockResidualResistanceFactor(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:611
static void registerOutputModules(Model &model, Simulator &simulator)
Register all polymer specific VTK and ECL output modules.
Definition: blackoilpolymermodules.hh:402
static Scalar computeUpdateError(const PrimaryVariables &, const EqVector &)
Return how much a Newton-Raphson update is considered an error.
Definition: blackoilpolymermodules.hh:566
static const BlackOilPolymerParams< Scalar >::PlyvmhCoefficients & plyvmhCoefficients(const ElementContext &elemCtx, const unsigned scvIdx, const unsigned timeIdx)
Definition: blackoilpolymermodules.hh:681
static const Scalar shrate(unsigned pvtnumRegionIdx)
Definition: blackoilpolymermodules.hh:699
static bool primaryVarApplies(unsigned pvIdx)
Definition: blackoilpolymermodules.hh:409
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilpolymermodules.hh:473
static const Scalar plyrockRockDensityFactor(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:619
static const Scalar plymaxMaxConcentration(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:664
static bool hasPlyshlog()
Definition: blackoilpolymermodules.hh:689
static std::string eqName(unsigned eqIdx)
Definition: blackoilpolymermodules.hh:453
static const Scalar plyrockDeadPoreVolume(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:603
static const Scalar plymixparToddLongstaff(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:672
static TabulatedTwoDFunction & getSkprwatTable(const int tableNumber)
get the SKPRWAT table
Definition: blackoilpolymermodules.hh:364
static void deserializeEntity(Model &model, std::istream &instream, const DofEntity &dof)
Definition: blackoilpolymermodules.hh:587
static const Scalar plyrockMaxAdsorbtion(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:635
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:512
static Evaluation computeShearFactor(const Evaluation &polymerConcentration, unsigned pvtnumRegionIdx, const Evaluation &v0)
Computes the shear factor.
Definition: blackoilpolymermodules.hh:711
VTK output module for the black oil model's polymer related quantities.
Definition: vtkblackoilpolymermodule.hh:71
Definition: blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:235
Definition: blackoilpolymerparams.hh:99
Definition: blackoilpolymerparams.hh:106
Struct holding the parameters for the BlackOilPolymerModule class.
Definition: blackoilpolymerparams.hh:41
AdsorptionBehaviour
Definition: blackoilpolymerparams.hh:45
IntervalTabulated2DFunction< Scalar > TabulatedTwoDFunction
Definition: blackoilpolymerparams.hh:43
Tabulated1DFunction< Scalar > TabulatedFunction
Definition: blackoilpolymerparams.hh:42