70 using TabulatedFunction =
typename BlackOilPolymerParams<Scalar>::TabulatedFunction;
71 using TabulatedTwoDFunction =
typename BlackOilPolymerParams<Scalar>::TabulatedTwoDFunction;
73 static constexpr unsigned polymerConcentrationIdx = Indices::polymerConcentrationIdx;
74 static constexpr unsigned polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
75 static constexpr unsigned contiPolymerEqIdx = Indices::contiPolymerEqIdx;
76 static constexpr unsigned contiPolymerMolarWeightEqIdx = Indices::contiPolymerMWEqIdx;
77 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
84 static constexpr unsigned numPhases = FluidSystem::numPhases;
98 if (
iterTable != params_.plymwinjTables_.end()) {
102 throw std::runtime_error(
" the PLYMWINJ table " + std::to_string(
tableNumber) +
" does not exist\n");
112 if (
iterTable != params_.skprwatTables_.end()) {
116 throw std::runtime_error(
" the SKPRWAT table " + std::to_string(
tableNumber) +
" does not exist\n");
127 if (
iterTable != params_.skprpolyTables_.end()) {
131 throw std::runtime_error(
" the SKPRPOLY table " + std::to_string(
tableNumber) +
" does not exist\n");
140 if constexpr (enablePolymer)
148 Simulator& simulator)
150 if constexpr (enablePolymer)
154 static bool primaryVarApplies(
unsigned pvIdx)
156 if constexpr (enablePolymer) {
157 if constexpr (enablePolymerMolarWeight)
158 return pvIdx == polymerConcentrationIdx ||
pvIdx == polymerMoleWeightIdx;
160 return pvIdx == polymerConcentrationIdx;
166 static std::string primaryVarName(
unsigned pvIdx)
170 if (
pvIdx == polymerConcentrationIdx) {
171 return "polymer_waterconcentration";
174 return "polymer_molecularweight";
183 return static_cast<Scalar
>(1.0);
186 static bool eqApplies(
unsigned eqIdx)
188 if constexpr (enablePolymer) {
189 if constexpr (enablePolymerMolarWeight)
190 return eqIdx == contiPolymerEqIdx ||
eqIdx == contiPolymerMolarWeightEqIdx;
192 return eqIdx == contiPolymerEqIdx;
198 static std::string eqName(
unsigned eqIdx)
202 if (
eqIdx == contiPolymerEqIdx)
203 return "conti^polymer";
205 return "conti^polymer_molecularweight";
213 return static_cast<Scalar
>(1.0);
217 template <
class LhsEval>
218 static void addStorage(Dune::FieldVector<LhsEval, numEq>&
storage,
219 const IntensiveQuantities& intQuants)
221 if constexpr (enablePolymer) {
222 const auto& fs = intQuants.fluidState();
234 * Toolbox::template
decay<LhsEval>(intQuants.polymerConcentration())
235 * (1.0 - Toolbox::template
decay<LhsEval>(intQuants.polymerDeadPoreVolume()));
240 * Toolbox::template
decay<LhsEval>(intQuants.polymerRockDensity())
241 * Toolbox::template
decay<LhsEval>(intQuants.polymerAdsorption());
248 if constexpr (enablePolymerMolarWeight) {
252 * Toolbox::template
decay<LhsEval> (intQuants.polymerMoleWeight());
263 if constexpr (enablePolymer) {
266 const unsigned upIdx =
extQuants.upstreamIndex(FluidSystem::waterPhaseIdx);
269 const unsigned contiWaterEqIdx = Indices::conti0EqIdx + Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
272 flux[contiPolymerEqIdx] =
274 *
up.fluidState().invB(waterPhaseIdx)
275 *
up.polymerViscosityCorrection()
277 *
up.polymerConcentration();
284 flux[contiPolymerEqIdx] =
297 if constexpr (enablePolymerMolarWeight) {
299 flux[contiPolymerMolarWeightEqIdx] =
300 flux[contiPolymerEqIdx]*
up.polymerMoleWeight();
302 flux[contiPolymerMolarWeightEqIdx] =
317 return static_cast<Scalar
>(0.0);
320 template <
class DofEntity>
323 if constexpr (enablePolymer) {
324 unsigned dofIdx = model.dofMapper().index(
dof);
325 const PrimaryVariables& priVars = model.solution(0)[dofIdx];
326 outstream << priVars[polymerConcentrationIdx];
327 outstream << priVars[polymerMoleWeightIdx];
331 template <
class DofEntity>
334 if constexpr (enablePolymer) {
335 unsigned dofIdx = model.dofMapper().index(
dof);
336 PrimaryVariables&
priVars0 = model.solution(0)[dofIdx];
337 PrimaryVariables&
priVars1 = model.solution(1)[dofIdx];
348 static const Scalar plyrockDeadPoreVolume(
const ElementContext& elemCtx,
356 static const Scalar plyrockResidualResistanceFactor(
const ElementContext& elemCtx,
364 static const Scalar plyrockRockDensityFactor(
const ElementContext& elemCtx,
372 static const Scalar plyrockAdsorbtionIndex(
const ElementContext& elemCtx,
380 static const Scalar plyrockMaxAdsorbtion(
const ElementContext& elemCtx,
388 static const TabulatedFunction& plyadsAdsorbedPolymer(
const ElementContext& elemCtx,
396 static const TabulatedFunction& plyviscViscosityMultiplierTable(
const ElementContext& elemCtx,
404 static const TabulatedFunction& plyviscViscosityMultiplierTable(
unsigned pvtnumRegionIdx)
409 static const Scalar plymaxMaxConcentration(
const ElementContext& elemCtx,
417 static const Scalar plymixparToddLongstaff(
const ElementContext& elemCtx,
425 static const typename BlackOilPolymerParams<Scalar>::PlyvmhCoefficients&
426 plyvmhCoefficients(
const ElementContext& elemCtx,
434 static bool hasPlyshlog()
436 return params_.hasPlyshlog_;
439 static bool hasShrate()
441 return params_.hasShrate_;
455 template <
class Evaluation>
458 const Evaluation&
v0)
465 const Scalar eps = 1
e-14;
468 return ToolboxLocal::createConstant(
v0, 1.0);
474 return ToolboxLocal::createConstant(
v0, 1.0);
509 bool converged =
false;
511 for (
int i = 0; i < 20; ++i) {
521 throw std::runtime_error(
"Not able to compute shear velocity. \n");
529 const Scalar molarMass()
const
566 static constexpr int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
567 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
569 static constexpr int polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
583 const auto linearizationType = elemCtx.linearizationType();
584 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx,
timeIdx);
585 polymerConcentration_ = priVars.makeEvaluation(polymerConcentrationIdx,
timeIdx, linearizationType);
586 if constexpr (enablePolymerMolarWeight) {
587 polymerMoleWeight_ = priVars.makeEvaluation(polymerMoleWeightIdx,
timeIdx, linearizationType);
592 const auto& plyadsAdsorbedPolymer = PolymerModule::plyadsAdsorbedPolymer(elemCtx, dofIdx,
timeIdx);
593 polymerAdsorption_ = plyadsAdsorbedPolymer.eval(polymerConcentration_,
true);
594 if (
static_cast<int>(PolymerModule::plyrockAdsorbtionIndex(elemCtx, dofIdx,
timeIdx)) ==
597 const Scalar& maxPolymerAdsorption = elemCtx.problem().maxPolymerAdsorption(elemCtx, dofIdx,
timeIdx);
598 polymerAdsorption_ = std::max(Evaluation(maxPolymerAdsorption) , polymerAdsorption_);
606 if constexpr (!enablePolymerMolarWeight) {
607 const Scalar
cmax = PolymerModule::plymaxMaxConcentration(elemCtx, dofIdx,
timeIdx);
608 const auto& fs = asImp_().fluidState_;
609 const Evaluation&
muWater = fs.viscosity(waterPhaseIdx);
614 const Scalar plymixparToddLongstaff = PolymerModule::plymixparToddLongstaff(elemCtx, dofIdx,
timeIdx);
619 const Evaluation
cbar = polymerConcentration_ /
cmax;
626 const auto& plyvmhCoefficients = PolymerModule::plyvmhCoefficients(elemCtx, dofIdx,
timeIdx);
627 const Scalar k_mh = plyvmhCoefficients.k_mh;
628 const Scalar a_mh = plyvmhCoefficients.a_mh;
629 const Scalar gamma = plyvmhCoefficients.gamma;
630 const Scalar kappa = plyvmhCoefficients.kappa;
636 waterViscosityCorrection_ = 1.0 / (1.0 + gamma * (x + kappa * x * x));
637 polymerViscosityCorrection_ = 1.0;
641 asImp_().mobility_[waterPhaseIdx] *= waterViscosityCorrection_ /
resistanceFactor;
644 polymerDeadPoreVolume_ = PolymerModule::plyrockDeadPoreVolume(elemCtx, dofIdx,
timeIdx);
645 polymerRockDensity_ = PolymerModule::plyrockRockDensityFactor(elemCtx, dofIdx,
timeIdx);
648 const Evaluation& polymerConcentration()
const
649 {
return polymerConcentration_; }
651 const Evaluation& polymerMoleWeight()
const
653 if constexpr (enablePolymerMolarWeight)
654 return polymerMoleWeight_;
656 throw std::logic_error(
"polymerMoleWeight() is called but polymer milecular weight is disabled");
659 const Scalar& polymerDeadPoreVolume()
const
660 {
return polymerDeadPoreVolume_; }
662 const Evaluation& polymerAdsorption()
const
663 {
return polymerAdsorption_; }
665 const Scalar& polymerRockDensity()
const
666 {
return polymerRockDensity_; }
669 const Evaluation& polymerViscosityCorrection()
const
670 {
return polymerViscosityCorrection_; }
673 const Evaluation& waterViscosityCorrection()
const
674 {
return waterViscosityCorrection_; }
678 Implementation& asImp_()
679 {
return *
static_cast<Implementation*
>(
this); }
681 Evaluation polymerConcentration_;
683 Evaluation polymerMoleWeight_;
684 Scalar polymerDeadPoreVolume_;
685 Scalar polymerRockDensity_;
686 Evaluation polymerAdsorption_;
687 Evaluation polymerViscosityCorrection_;
688 Evaluation waterViscosityCorrection_;