68 using Toolbox = MathToolbox<Evaluation>;
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;
80 static constexpr unsigned enablePolymer = enablePolymerV;
84 static constexpr unsigned numPhases = FluidSystem::numPhases;
97 const auto iterTable = params_.plymwinjTables_.find(tableNumber);
98 if (iterTable != params_.plymwinjTables_.end()) {
99 return iterTable->second;
102 throw std::runtime_error(
" the PLYMWINJ table " + std::to_string(tableNumber) +
" does not exist\n");
111 const auto iterTable = params_.skprwatTables_.find(tableNumber);
112 if (iterTable != params_.skprwatTables_.end()) {
113 return iterTable->second;
116 throw std::runtime_error(
" the SKPRWAT table " + std::to_string(tableNumber) +
" does not exist\n");
126 const auto iterTable = params_.skprpolyTables_.find(tableNumber);
127 if (iterTable != params_.skprpolyTables_.end()) {
128 return iterTable->second;
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)
168 assert(primaryVarApplies(pvIdx));
170 if (pvIdx == polymerConcentrationIdx) {
171 return "polymer_waterconcentration";
174 return "polymer_molecularweight";
178 static Scalar primaryVarWeight([[maybe_unused]]
unsigned pvIdx)
180 assert(primaryVarApplies(pvIdx));
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)
200 assert(eqApplies(eqIdx));
202 if (eqIdx == contiPolymerEqIdx)
203 return "conti^polymer";
205 return "conti^polymer_molecularweight";
208 static Scalar eqWeight([[maybe_unused]]
unsigned eqIdx)
210 assert(eqApplies(eqIdx));
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();
224 LhsEval surfaceVolumeWater =
225 Toolbox::template decay<LhsEval>(fs.saturation(waterPhaseIdx))
226 * Toolbox::template decay<LhsEval>(fs.invB(waterPhaseIdx))
227 * Toolbox::template decay<LhsEval>(intQuants.porosity());
230 surfaceVolumeWater = max(surfaceVolumeWater, 1e-10);
233 const LhsEval massPolymer = surfaceVolumeWater
234 * Toolbox::template decay<LhsEval>(intQuants.polymerConcentration())
235 * (1.0 - Toolbox::template decay<LhsEval>(intQuants.polymerDeadPoreVolume()));
238 const LhsEval adsorptionPolymer =
239 Toolbox::template decay<LhsEval>(1.0 - intQuants.porosity())
240 * Toolbox::template decay<LhsEval>(intQuants.polymerRockDensity())
241 * Toolbox::template decay<LhsEval>(intQuants.polymerAdsorption());
243 LhsEval accumulationPolymer = massPolymer + adsorptionPolymer;
245 storage[contiPolymerEqIdx] += accumulationPolymer;
248 if constexpr (enablePolymerMolarWeight) {
249 accumulationPolymer = max(accumulationPolymer, 1e-10);
251 storage[contiPolymerMolarWeightEqIdx] += accumulationPolymer
252 * Toolbox::template decay<LhsEval> (intQuants.polymerMoleWeight());
257 static void computeFlux([[maybe_unused]] RateVector& flux,
258 [[maybe_unused]]
const ElementContext& elemCtx,
259 [[maybe_unused]]
unsigned scvfIdx,
260 [[maybe_unused]]
unsigned timeIdx)
263 if constexpr (enablePolymer) {
264 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
266 const unsigned upIdx = extQuants.upstreamIndex(FluidSystem::waterPhaseIdx);
267 const unsigned inIdx = extQuants.interiorIndex();
268 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
269 const unsigned contiWaterEqIdx = Indices::conti0EqIdx + Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
271 if (upIdx == inIdx) {
272 flux[contiPolymerEqIdx] =
273 extQuants.volumeFlux(waterPhaseIdx)
274 *up.fluidState().invB(waterPhaseIdx)
275 *up.polymerViscosityCorrection()
276 /extQuants.polymerShearFactor()
277 *up.polymerConcentration();
280 flux[contiWaterEqIdx] /=
281 extQuants.waterShearFactor();
284 flux[contiPolymerEqIdx] =
285 extQuants.volumeFlux(waterPhaseIdx)
286 *decay<Scalar>(up.fluidState().invB(waterPhaseIdx))
287 *decay<Scalar>(up.polymerViscosityCorrection())
288 /decay<Scalar>(extQuants.polymerShearFactor())
289 *decay<Scalar>(up.polymerConcentration());
292 flux[contiWaterEqIdx] /=
293 decay<Scalar>(extQuants.waterShearFactor());
297 if constexpr (enablePolymerMolarWeight) {
299 flux[contiPolymerMolarWeightEqIdx] =
300 flux[contiPolymerEqIdx]*up.polymerMoleWeight();
302 flux[contiPolymerMolarWeightEqIdx] =
303 flux[contiPolymerEqIdx]*decay<Scalar>(up.polymerMoleWeight());
317 return static_cast<Scalar
>(0.0);
320 template <
class DofEntity>
321 static void serializeEntity(
const Model& model, std::ostream& outstream,
const DofEntity& dof)
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>
332 static void deserializeEntity(Model& model, std::istream& instream,
const DofEntity& dof)
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];
339 instream >> priVars0[polymerConcentrationIdx];
340 instream >> priVars0[polymerMoleWeightIdx];
343 priVars1[polymerConcentrationIdx] = priVars0[polymerConcentrationIdx];
344 priVars1[polymerMoleWeightIdx] = priVars0[polymerMoleWeightIdx];
348 static const Scalar plyrockDeadPoreVolume(
const ElementContext& elemCtx,
352 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
353 return params_.plyrockDeadPoreVolume_[satnumRegionIdx];
356 static const Scalar plyrockResidualResistanceFactor(
const ElementContext& elemCtx,
360 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
361 return params_.plyrockResidualResistanceFactor_[satnumRegionIdx];
364 static const Scalar plyrockRockDensityFactor(
const ElementContext& elemCtx,
368 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
369 return params_.plyrockRockDensityFactor_[satnumRegionIdx];
372 static const Scalar plyrockAdsorbtionIndex(
const ElementContext& elemCtx,
376 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
377 return params_.plyrockAdsorbtionIndex_[satnumRegionIdx];
380 static const Scalar plyrockMaxAdsorbtion(
const ElementContext& elemCtx,
384 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
385 return params_.plyrockMaxAdsorbtion_[satnumRegionIdx];
388 static const TabulatedFunction& plyadsAdsorbedPolymer(
const ElementContext& elemCtx,
392 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
393 return params_.plyadsAdsorbedPolymer_[satnumRegionIdx];
396 static const TabulatedFunction& plyviscViscosityMultiplierTable(
const ElementContext& elemCtx,
400 unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
401 return params_.plyviscViscosityMultiplierTable_[pvtnumRegionIdx];
404 static const TabulatedFunction& plyviscViscosityMultiplierTable(
unsigned pvtnumRegionIdx)
406 return params_.plyviscViscosityMultiplierTable_[pvtnumRegionIdx];
409 static const Scalar plymaxMaxConcentration(
const ElementContext& elemCtx,
413 unsigned polymerMixRegionIdx = elemCtx.problem().plmixnumRegionIndex(elemCtx, scvIdx, timeIdx);
414 return params_.plymaxMaxConcentration_[polymerMixRegionIdx];
417 static const Scalar plymixparToddLongstaff(
const ElementContext& elemCtx,
421 unsigned polymerMixRegionIdx = elemCtx.problem().plmixnumRegionIndex(elemCtx, scvIdx, timeIdx);
422 return params_.plymixparToddLongstaff_[polymerMixRegionIdx];
425 static const typename BlackOilPolymerParams<Scalar>::PlyvmhCoefficients&
426 plyvmhCoefficients(
const ElementContext& elemCtx,
427 const unsigned scvIdx,
428 const unsigned timeIdx)
430 const unsigned polymerMixRegionIdx = elemCtx.problem().plmixnumRegionIndex(elemCtx, scvIdx, timeIdx);
431 return params_.plyvmhCoefficients_[polymerMixRegionIdx];
434 static bool hasPlyshlog()
436 return params_.hasPlyshlog_;
439 static bool hasShrate()
441 return params_.hasShrate_;
444 static const Scalar shrate(
unsigned pvtnumRegionIdx)
446 return params_.shrate_[pvtnumRegionIdx];
455 template <
class Evaluation>
457 unsigned pvtnumRegionIdx,
458 const Evaluation& v0)
460 using ToolboxLocal = MathToolbox<Evaluation>;
462 const auto& viscosityMultiplierTable = params_.plyviscViscosityMultiplierTable_[pvtnumRegionIdx];
463 Scalar viscosityMultiplier = viscosityMultiplierTable.eval(scalarValue(polymerConcentration),
true);
465 const Scalar eps = 1e-14;
467 if (std::abs((viscosityMultiplier - 1.0)) < eps)
468 return ToolboxLocal::createConstant(v0, 1.0);
470 const std::vector<Scalar>& shearEffectRefLogVelocity = params_.plyshlogShearEffectRefLogVelocity_[pvtnumRegionIdx];
471 auto v0AbsLog = log(abs(v0));
473 if (v0AbsLog < shearEffectRefLogVelocity[0])
474 return ToolboxLocal::createConstant(v0, 1.0);
480 const std::vector<Scalar>& shearEffectRefMultiplier = params_.plyshlogShearEffectRefMultiplier_[pvtnumRegionIdx];
481 size_t numTableEntries = shearEffectRefLogVelocity.size();
482 assert(shearEffectRefMultiplier.size() == numTableEntries);
484 std::vector<Scalar> shearEffectMultiplier(numTableEntries, 1.0);
485 for (
size_t i = 0; i < numTableEntries; ++i) {
486 shearEffectMultiplier[i] = (1.0 + (viscosityMultiplier - 1.0)*shearEffectRefMultiplier[i]) / viscosityMultiplier;
487 shearEffectMultiplier[i] = log(shearEffectMultiplier[i]);
491 TabulatedFunction logShearEffectMultiplier = TabulatedFunction(numTableEntries, shearEffectRefLogVelocity, shearEffectMultiplier,
false);
498 auto F = [&logShearEffectMultiplier, &v0AbsLog](
const Evaluation& u) {
499 return u + logShearEffectMultiplier.eval(u,
true) - v0AbsLog;
502 auto dF = [&logShearEffectMultiplier](
const Evaluation& u) {
503 return 1 + logShearEffectMultiplier.evalDerivative(u,
true);
509 bool converged =
false;
511 for (
int i = 0; i < 20; ++i) {
515 if (std::abs(scalarValue(f)) < 1e-12) {
521 throw std::runtime_error(
"Not able to compute shear velocity. \n");
525 return exp(logShearEffectMultiplier.eval(u,
true));
529 const Scalar molarMass()
const
535 static BlackOilPolymerParams<Scalar> params_;
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);
591 const Scalar& maxAdsorbtion = PolymerModule::plyrockMaxAdsorbtion(elemCtx, dofIdx, timeIdx);
592 const auto& plyadsAdsorbedPolymer = PolymerModule::plyadsAdsorbedPolymer(elemCtx, dofIdx, timeIdx);
593 polymerAdsorption_ = plyadsAdsorbedPolymer.eval(polymerConcentration_,
true);
595 const Scalar& maxPolymerAdsorption = elemCtx.problem().maxPolymerAdsorption(elemCtx, dofIdx, timeIdx);
596 polymerAdsorption_ = std::max(Evaluation(maxPolymerAdsorption) , polymerAdsorption_);
600 const Scalar& residualResistanceFactor = PolymerModule::plyrockResidualResistanceFactor(elemCtx, dofIdx, timeIdx);
601 const Evaluation resistanceFactor = 1.0 + (residualResistanceFactor - 1.0) * polymerAdsorption_ / maxAdsorbtion;
604 if constexpr (!enablePolymerMolarWeight) {
605 const Scalar cmax = PolymerModule::plymaxMaxConcentration(elemCtx, dofIdx, timeIdx);
606 const auto& fs = asImp_().fluidState_;
607 const Evaluation& muWater = fs.viscosity(waterPhaseIdx);
608 const auto& viscosityMultiplier = PolymerModule::plyviscViscosityMultiplierTable(elemCtx, dofIdx, timeIdx);
609 const Evaluation viscosityMixture = viscosityMultiplier.eval(polymerConcentration_,
true) * muWater;
612 const Scalar plymixparToddLongstaff = PolymerModule::plymixparToddLongstaff(elemCtx, dofIdx, timeIdx);
613 const Evaluation viscosityPolymer = viscosityMultiplier.eval(cmax,
true) * muWater;
614 const Evaluation viscosityPolymerEffective = pow(viscosityMixture, plymixparToddLongstaff) * pow(viscosityPolymer, 1.0 - plymixparToddLongstaff);
615 const Evaluation viscosityWaterEffective = pow(viscosityMixture, plymixparToddLongstaff) * pow(muWater, 1.0 - plymixparToddLongstaff);
617 const Evaluation cbar = polymerConcentration_ / cmax;
619 waterViscosityCorrection_ = muWater * ((1.0 - cbar) / viscosityWaterEffective + cbar / viscosityPolymerEffective);
621 polymerViscosityCorrection_ = (muWater / waterViscosityCorrection_) / viscosityPolymerEffective;
624 const auto& plyvmhCoefficients = PolymerModule::plyvmhCoefficients(elemCtx, dofIdx, timeIdx);
625 const Scalar k_mh = plyvmhCoefficients.k_mh;
626 const Scalar a_mh = plyvmhCoefficients.a_mh;
627 const Scalar gamma = plyvmhCoefficients.gamma;
628 const Scalar kappa = plyvmhCoefficients.kappa;
632 const Evaluation intrinsicViscosity = k_mh * pow(polymerMoleWeight_ * 1000., a_mh);
633 const Evaluation x = polymerConcentration_ * intrinsicViscosity;
634 waterViscosityCorrection_ = 1.0 / (1.0 + gamma * (x + kappa * x * x));
635 polymerViscosityCorrection_ = 1.0;
639 asImp_().mobility_[waterPhaseIdx] *= waterViscosityCorrection_ / resistanceFactor;
642 polymerDeadPoreVolume_ = PolymerModule::plyrockDeadPoreVolume(elemCtx, dofIdx, timeIdx);
643 polymerRockDensity_ = PolymerModule::plyrockRockDensityFactor(elemCtx, dofIdx, timeIdx);
786 waterShearFactor_ = 1.0;
787 polymerShearFactor_ = 1.0;
789 if (!PolymerModule::hasPlyshlog())
792 const ExtensiveQuantities& extQuants = asImp_();
793 unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
794 unsigned interiorDofIdx = extQuants.interiorIndex();
795 unsigned exteriorDofIdx = extQuants.exteriorIndex();
796 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
797 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
798 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
801 Evaluation poroAvg = intQuantsIn.porosity()*0.5 + intQuantsEx.porosity()*0.5;
802 unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvfIdx, timeIdx);
803 const Evaluation& Sw = up.fluidState().saturation(waterPhaseIdx);
804 unsigned cellIdx = elemCtx.globalSpaceIndex(scvfIdx, timeIdx);
805 const auto& materialLawManager = elemCtx.problem().materialLawManager();
806 const auto& scaledDrainageInfo =
807 materialLawManager->oilWaterScaledEpsInfoDrainage(cellIdx);
808 const Scalar& Swcr = scaledDrainageInfo.Swcr;
811 Evaluation denom = max(poroAvg * (Sw - Swcr), 1e-12);
812 Evaluation waterVolumeVelocity = extQuants.volumeFlux(waterPhaseIdx) / denom;
815 if (PolymerModule::hasShrate()) {
816 const Evaluation& relWater = up.relativePermeability(waterPhaseIdx);
817 Scalar trans = elemCtx.problem().transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
819 Scalar faceArea = elemCtx.stencil(timeIdx).interiorFace(scvfIdx).area();
820 auto dist = elemCtx.pos(interiorDofIdx, timeIdx) - elemCtx.pos(exteriorDofIdx, timeIdx);
822 Scalar absPerm = trans / faceArea * dist.two_norm();
823 waterVolumeVelocity *=
824 PolymerModule::shrate(pvtnumRegionIdx)*sqrt(poroAvg*Sw / (relWater*absPerm));
825 assert(isfinite(waterVolumeVelocity));
833 waterVolumeVelocity);
834 polymerShearFactor_ =
837 waterVolumeVelocity*up.polymerViscosityCorrection());