blackoilpolymermodules.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef EWOMS_BLACK_OIL_POLYMER_MODULE_HH
29#define EWOMS_BLACK_OIL_POLYMER_MODULE_HH
30
31#include <dune/common/fvector.hh>
32
33#include <opm/common/OpmLog/OpmLog.hpp>
34
35#include <opm/material/common/MathToolbox.hpp>
36
39
41
43
44#include <algorithm>
45#include <cassert>
46#include <cmath>
47#include <cstddef>
48#include <istream>
49#include <memory>
50#include <ostream>
51#include <stdexcept>
52#include <string>
53#include <vector>
54
55namespace Opm {
56
62template <class TypeTag, bool enablePolymerV = getPropValue<TypeTag, Properties::EnablePolymer>()>
64{
77
78 using Toolbox = MathToolbox<Evaluation>;
79
80 using TabulatedFunction = typename BlackOilPolymerParams<Scalar>::TabulatedFunction;
81 using TabulatedTwoDFunction = typename BlackOilPolymerParams<Scalar>::TabulatedTwoDFunction;
82
83 static constexpr unsigned polymerConcentrationIdx = Indices::polymerConcentrationIdx;
84 static constexpr unsigned polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
85 static constexpr unsigned contiPolymerEqIdx = Indices::contiPolymerEqIdx;
86 static constexpr unsigned contiPolymerMolarWeightEqIdx = Indices::contiPolymerMWEqIdx;
87 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
88
89 static constexpr unsigned enablePolymer = enablePolymerV;
90 static constexpr bool enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>();
91
92 static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
93 static constexpr unsigned numPhases = FluidSystem::numPhases;
94
95public:
98 {
99 params_ = params;
100 }
101
105 static TabulatedTwoDFunction& getPlymwinjTable(const int tableNumber)
106 {
107 const auto iterTable = params_.plymwinjTables_.find(tableNumber);
108 if (iterTable != params_.plymwinjTables_.end()) {
109 return iterTable->second;
110 }
111 else {
112 throw std::runtime_error(" the PLYMWINJ table " + std::to_string(tableNumber) + " does not exist\n");
113 }
114 }
115
119 static TabulatedTwoDFunction& getSkprwatTable(const int tableNumber)
120 {
121 const auto iterTable = params_.skprwatTables_.find(tableNumber);
122 if (iterTable != params_.skprwatTables_.end()) {
123 return iterTable->second;
124 }
125 else {
126 throw std::runtime_error(" the SKPRWAT table " + std::to_string(tableNumber) + " does not exist\n");
127 }
128 }
129
134 getSkprpolyTable(const int tableNumber)
135 {
136 const auto iterTable = params_.skprpolyTables_.find(tableNumber);
137 if (iterTable != params_.skprpolyTables_.end()) {
138 return iterTable->second;
139 }
140 else {
141 throw std::runtime_error(" the SKPRPOLY table " + std::to_string(tableNumber) + " does not exist\n");
142 }
143 }
144
148 static void registerParameters()
149 {
150 if constexpr (enablePolymer) {
152 }
153 }
154
158 static void registerOutputModules(Model& model,
159 Simulator& simulator)
160 {
161 if constexpr (enablePolymer) {
162 model.addOutputModule(std::make_unique<VtkBlackOilPolymerModule<TypeTag>>(simulator));
163 }
164 }
165
166 static bool primaryVarApplies(unsigned pvIdx)
167 {
168 if constexpr (enablePolymer) {
169 if constexpr (enablePolymerMolarWeight) {
170 return pvIdx == polymerConcentrationIdx || pvIdx == polymerMoleWeightIdx;
171 }
172 else {
173 return pvIdx == polymerConcentrationIdx;
174 }
175 }
176 else {
177 return false;
178 }
179 }
180
181 static std::string primaryVarName(unsigned pvIdx)
182 {
183 assert(primaryVarApplies(pvIdx));
184
185 if (pvIdx == polymerConcentrationIdx) {
186 return "polymer_waterconcentration";
187 }
188 else {
189 return "polymer_molecularweight";
190 }
191 }
192
193 static Scalar primaryVarWeight([[maybe_unused]] unsigned pvIdx)
194 {
195 assert(primaryVarApplies(pvIdx));
196
197 // TODO: it may be beneficial to chose this differently.
198 return static_cast<Scalar>(1.0);
199 }
200
201 static bool eqApplies(unsigned eqIdx)
202 {
203 if constexpr (enablePolymer) {
204 if constexpr (enablePolymerMolarWeight) {
205 return eqIdx == contiPolymerEqIdx || eqIdx == contiPolymerMolarWeightEqIdx;
206 }
207 else {
208 return eqIdx == contiPolymerEqIdx;
209 }
210 }
211 else {
212 return false;
213 }
214 }
215
216 static std::string eqName(unsigned eqIdx)
217 {
218 assert(eqApplies(eqIdx));
219
220 if (eqIdx == contiPolymerEqIdx) {
221 return "conti^polymer";
222 }
223 else {
224 return "conti^polymer_molecularweight";
225 }
226 }
227
228 static Scalar eqWeight([[maybe_unused]] unsigned eqIdx)
229 {
230 assert(eqApplies(eqIdx));
231
232 // TODO: it may be beneficial to chose this differently.
233 return static_cast<Scalar>(1.0);
234 }
235
236 // must be called after water storage is computed
237 template <class LhsEval>
238 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
239 const IntensiveQuantities& intQuants)
240 {
241 if constexpr (enablePolymer) {
242 const auto& fs = intQuants.fluidState();
243
244 // avoid singular matrix if no water is present.
245 const LhsEval surfaceVolumeWater =
246 max(Toolbox::template decay<LhsEval>(fs.saturation(waterPhaseIdx)) *
247 Toolbox::template decay<LhsEval>(fs.invB(waterPhaseIdx)) *
248 Toolbox::template decay<LhsEval>(intQuants.porosity()),
249 1e-10);
250
251 // polymer in water phase
252 const LhsEval massPolymer =
253 surfaceVolumeWater *
254 Toolbox::template decay<LhsEval>(intQuants.polymerConcentration()) *
255 (1.0 - Toolbox::template decay<LhsEval>(intQuants.polymerDeadPoreVolume()));
256
257 // polymer in solid phase
258 const LhsEval adsorptionPolymer =
259 Toolbox::template decay<LhsEval>(1.0 - intQuants.porosity()) *
260 Toolbox::template decay<LhsEval>(intQuants.polymerRockDensity()) *
261 Toolbox::template decay<LhsEval>(intQuants.polymerAdsorption());
262
263 LhsEval accumulationPolymer = massPolymer + adsorptionPolymer;
264
265 storage[contiPolymerEqIdx] += accumulationPolymer;
266
267 // tracking the polymer molecular weight
268 if constexpr (enablePolymerMolarWeight) {
269 accumulationPolymer = max(accumulationPolymer, 1e-10);
270
271 storage[contiPolymerMolarWeightEqIdx] +=
272 accumulationPolymer * Toolbox::template decay<LhsEval>(intQuants.polymerMoleWeight());
273 }
274 }
275 }
276
277 static void computeFlux([[maybe_unused]] RateVector& flux,
278 [[maybe_unused]] const ElementContext& elemCtx,
279 [[maybe_unused]] unsigned scvfIdx,
280 [[maybe_unused]] unsigned timeIdx)
281 {
282 if constexpr (enablePolymer) {
283 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
284
285 const unsigned upIdx = extQuants.upstreamIndex(FluidSystem::waterPhaseIdx);
286 const unsigned inIdx = extQuants.interiorIndex();
287 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
288 const unsigned contiWaterEqIdx =
289 Indices::conti0EqIdx + Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
290
291 if (upIdx == inIdx) {
292 flux[contiPolymerEqIdx] =
293 extQuants.volumeFlux(waterPhaseIdx) *
294 up.fluidState().invB(waterPhaseIdx) *
295 up.polymerViscosityCorrection() /
296 extQuants.polymerShearFactor() *
297 up.polymerConcentration();
298
299 // modify water
300 flux[contiWaterEqIdx] /= extQuants.waterShearFactor();
301 }
302 else {
303 flux[contiPolymerEqIdx] =
304 extQuants.volumeFlux(waterPhaseIdx) *
305 decay<Scalar>(up.fluidState().invB(waterPhaseIdx)) *
306 decay<Scalar>(up.polymerViscosityCorrection()) /
307 decay<Scalar>(extQuants.polymerShearFactor()) *
308 decay<Scalar>(up.polymerConcentration());
309
310 // modify water
311 flux[contiWaterEqIdx] /= decay<Scalar>(extQuants.waterShearFactor());
312 }
313
314 // flux related to transport of polymer molecular weight
315 if constexpr (enablePolymerMolarWeight) {
316 if (upIdx == inIdx) {
317 flux[contiPolymerMolarWeightEqIdx] =
318 flux[contiPolymerEqIdx] * up.polymerMoleWeight();
319 }
320 else {
321 flux[contiPolymerMolarWeightEqIdx] =
322 flux[contiPolymerEqIdx] * decay<Scalar>(up.polymerMoleWeight());
323 }
324 }
325 }
326 }
327
331 static Scalar computeUpdateError(const PrimaryVariables&,
332 const EqVector&)
333 {
334 // do not consider consider the cange of polymer primary variables for
335 // convergence
336 // TODO: maybe this should be changed
337 return static_cast<Scalar>(0.0);
338 }
339
340 template <class DofEntity>
341 static void serializeEntity(const Model& model, std::ostream& outstream, const DofEntity& dof)
342 {
343 if constexpr (enablePolymer) {
344 const unsigned dofIdx = model.dofMapper().index(dof);
345 const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
346 outstream << priVars[polymerConcentrationIdx];
347 outstream << priVars[polymerMoleWeightIdx];
348 }
349 }
350
351 template <class DofEntity>
352 static void deserializeEntity(Model& model, std::istream& instream, const DofEntity& dof)
353 {
354 if constexpr (enablePolymer) {
355 const unsigned dofIdx = model.dofMapper().index(dof);
356 PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
357 PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
358
359 instream >> priVars0[polymerConcentrationIdx];
360 instream >> priVars0[polymerMoleWeightIdx];
361
362 // set the primary variables for the beginning of the current time step.
363 priVars1[polymerConcentrationIdx] = priVars0[polymerConcentrationIdx];
364 priVars1[polymerMoleWeightIdx] = priVars0[polymerMoleWeightIdx];
365 }
366 }
367
368 static Scalar plyrockDeadPoreVolume(const ElementContext& elemCtx,
369 unsigned scvIdx,
370 unsigned timeIdx)
371 {
372 const unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
373 return params_.plyrockDeadPoreVolume_[satnumRegionIdx];
374 }
375
376 static Scalar plyrockResidualResistanceFactor(const ElementContext& elemCtx,
377 unsigned scvIdx,
378 unsigned timeIdx)
379 {
380 const unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
381 return params_.plyrockResidualResistanceFactor_[satnumRegionIdx];
382 }
383
384 static Scalar plyrockRockDensityFactor(const ElementContext& elemCtx,
385 unsigned scvIdx,
386 unsigned timeIdx)
387 {
388 const unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
389 return params_.plyrockRockDensityFactor_[satnumRegionIdx];
390 }
391
392 static Scalar plyrockAdsorbtionIndex(const ElementContext& elemCtx,
393 unsigned scvIdx,
394 unsigned timeIdx)
395 {
396 const unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
397 return params_.plyrockAdsorbtionIndex_[satnumRegionIdx];
398 }
399
400 static Scalar plyrockMaxAdsorbtion(const ElementContext& elemCtx,
401 unsigned scvIdx,
402 unsigned timeIdx)
403 {
404 const unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
405 return params_.plyrockMaxAdsorbtion_[satnumRegionIdx];
406 }
407
408 static const TabulatedFunction&
409 plyadsAdsorbedPolymer(const ElementContext& elemCtx,
410 unsigned scvIdx,
411 unsigned timeIdx)
412 {
413 const unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
414 return params_.plyadsAdsorbedPolymer_[satnumRegionIdx];
415 }
416
417 static const TabulatedFunction&
418 plyviscViscosityMultiplierTable(const ElementContext& elemCtx,
419 unsigned scvIdx,
420 unsigned timeIdx)
421 {
422 unsigned pvtnumRegionIdx =
423 elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
424 return params_.plyviscViscosityMultiplierTable_[pvtnumRegionIdx];
425 }
426
427 static const TabulatedFunction&
428 plyviscViscosityMultiplierTable(unsigned pvtnumRegionIdx)
429 { return params_.plyviscViscosityMultiplierTable_[pvtnumRegionIdx]; }
430
431 static Scalar plymaxMaxConcentration(const ElementContext& elemCtx,
432 unsigned scvIdx,
433 unsigned timeIdx)
434 {
435 const unsigned polymerMixRegionIdx =
436 elemCtx.problem().plmixnumRegionIndex(elemCtx, scvIdx, timeIdx);
437 return params_.plymaxMaxConcentration_[polymerMixRegionIdx];
438 }
439
440 static Scalar plymixparToddLongstaff(const ElementContext& elemCtx,
441 unsigned scvIdx,
442 unsigned timeIdx)
443 {
444 const unsigned polymerMixRegionIdx =
445 elemCtx.problem().plmixnumRegionIndex(elemCtx, scvIdx, timeIdx);
446 return params_.plymixparToddLongstaff_[polymerMixRegionIdx];
447 }
448
450 plyvmhCoefficients(const ElementContext& elemCtx,
451 const unsigned scvIdx,
452 const unsigned timeIdx)
453 {
454 const unsigned polymerMixRegionIdx =
455 elemCtx.problem().plmixnumRegionIndex(elemCtx, scvIdx, timeIdx);
456 return params_.plyvmhCoefficients_[polymerMixRegionIdx];
457 }
458
459 static bool hasPlyshlog()
460 { return params_.hasPlyshlog_; }
461
462 static bool hasShrate()
463 { return params_.hasShrate_; }
464
465 static Scalar shrate(unsigned pvtnumRegionIdx)
466 { return params_.shrate_[pvtnumRegionIdx]; }
467
474 template <class Evaluation>
475 static Evaluation computeShearFactor(const Evaluation& polymerConcentration,
476 unsigned pvtnumRegionIdx,
477 const Evaluation& v0)
478 {
479 using ToolboxLocal = MathToolbox<Evaluation>;
480
481 const auto& viscosityMultiplierTable = params_.plyviscViscosityMultiplierTable_[pvtnumRegionIdx];
482 const Scalar viscosityMultiplier =
483 viscosityMultiplierTable.eval(scalarValue(polymerConcentration), /*extrapolate=*/true);
484
485 const Scalar eps = 1e-14;
486 // return 1.0 if the polymer has no effect on the water.
487 if (std::abs((viscosityMultiplier - 1.0)) < eps) {
488 return ToolboxLocal::createConstant(v0, 1.0);
489 }
490
491 const std::vector<Scalar>& shearEffectRefLogVelocity =
492 params_.plyshlogShearEffectRefLogVelocity_[pvtnumRegionIdx];
493 const auto v0AbsLog = log(abs(v0));
494 // return 1.0 if the velocity /sharte is smaller than the first velocity entry.
495 if (v0AbsLog < shearEffectRefLogVelocity[0]) {
496 return ToolboxLocal::createConstant(v0, 1.0);
497 }
498
499 // compute shear factor from input
500 // Z = (1 + (P - 1) * M(v)) / P
501 // where M(v) is computed from user input
502 // and P = viscosityMultiplier
503 const std::vector<Scalar>& shearEffectRefMultiplier =
504 params_.plyshlogShearEffectRefMultiplier_[pvtnumRegionIdx];
505 const std::size_t numTableEntries = shearEffectRefLogVelocity.size();
506 assert(shearEffectRefMultiplier.size() == numTableEntries);
507
508 std::vector<Scalar> shearEffectMultiplier(numTableEntries, 1.0);
509 for (std::size_t i = 0; i < numTableEntries; ++i) {
510 shearEffectMultiplier[i] = (1.0 + (viscosityMultiplier - 1.0) *
511 shearEffectRefMultiplier[i]) / viscosityMultiplier;
512 shearEffectMultiplier[i] = log(shearEffectMultiplier[i]);
513 }
514 // store the logarithmic velocity and logarithmic multipliers in a table for easy look up and
515 // linear interpolation in the logarithmic space.
516 const TabulatedFunction logShearEffectMultiplier =
517 TabulatedFunction(numTableEntries, shearEffectRefLogVelocity,
518 shearEffectMultiplier, /*bool sortInputs =*/ false);
519
520 // Find sheared velocity (v) that satisfies
521 // F = log(v) + log (Z) - log(v0) = 0;
522
523 // Set up the function
524 // u = log(v)
525 auto F = [&logShearEffectMultiplier, &v0AbsLog](const Evaluation& u) {
526 return u + logShearEffectMultiplier.eval(u, true) - v0AbsLog;
527 };
528 // and its derivative
529 auto dF = [&logShearEffectMultiplier](const Evaluation& u) {
530 return 1 + logShearEffectMultiplier.evalDerivative(u, true);
531 };
532
533 // Solve F = 0 using Newton
534 // Use log(v0) as initial value for u
535 auto u = v0AbsLog;
536 bool converged = false;
537 // TODO make this into parameters
538 for (int i = 0; i < 20; ++i) {
539 const auto f = F(u);
540 const auto df = dF(u);
541 u -= f / df;
542 if (std::abs(scalarValue(f)) < 1e-12) {
543 converged = true;
544 break;
545 }
546 }
547 if (!converged) {
548 throw std::runtime_error("Not able to compute shear velocity. \n");
549 }
550
551 // return the shear factor
552 return exp(logShearEffectMultiplier.eval(u, /*extrapolate=*/true));
553 }
554
555 Scalar molarMass() const
556 { return 0.25; } // kg/mol
557
558private:
559 static BlackOilPolymerParams<Scalar> params_;
560};
561
562template <class TypeTag, bool enablePolymerV>
563BlackOilPolymerParams<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
564BlackOilPolymerModule<TypeTag, enablePolymerV>::params_;
565
566template <class TypeTag, bool enablePolymerV>
568
576template <class TypeTag>
577class BlackOilPolymerIntensiveQuantities<TypeTag, /*enablePolymerV=*/true>
578{
580
588
590
591 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
592 static constexpr int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
593 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
594 static constexpr bool enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>();
595 static constexpr int polymerMoleWeightIdx = Indices::polymerMoleWeightIdx;
596
597public:
603 void polymerPropertiesUpdate_(const ElementContext& elemCtx,
604 unsigned dofIdx,
605 unsigned timeIdx)
606 {
607 const auto linearizationType = elemCtx.linearizationType();
608 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
609 polymerConcentration_ = priVars.makeEvaluation(polymerConcentrationIdx, timeIdx, linearizationType);
610 if constexpr (enablePolymerMolarWeight) {
611 polymerMoleWeight_ = priVars.makeEvaluation(polymerMoleWeightIdx, timeIdx, linearizationType);
612 }
613
614 // permeability reduction due to polymer
615 const Scalar& maxAdsorbtion = PolymerModule::plyrockMaxAdsorbtion(elemCtx, dofIdx, timeIdx);
616 const auto& plyadsAdsorbedPolymer = PolymerModule::plyadsAdsorbedPolymer(elemCtx, dofIdx, timeIdx);
617 polymerAdsorption_ = plyadsAdsorbedPolymer.eval(polymerConcentration_, /*extrapolate=*/true);
618 if (static_cast<int>(PolymerModule::plyrockAdsorbtionIndex(elemCtx, dofIdx, timeIdx)) ==
620 {
621 const auto maxPolymerAdsorption =
622 elemCtx.problem().maxPolymerAdsorption(elemCtx, dofIdx, timeIdx);
623 polymerAdsorption_ = std::max(Evaluation(maxPolymerAdsorption), polymerAdsorption_);
624 }
625
626 // compute resitanceFactor
627 const Scalar& residualResistanceFactor =
628 PolymerModule::plyrockResidualResistanceFactor(elemCtx, dofIdx, timeIdx);
629 const Evaluation resistanceFactor = 1.0 + (residualResistanceFactor - 1.0) *
630 polymerAdsorption_ / maxAdsorbtion;
631
632 // compute effective viscosities
633 if constexpr (!enablePolymerMolarWeight) {
634 const Scalar cmax = PolymerModule::plymaxMaxConcentration(elemCtx, dofIdx, timeIdx);
635 const auto& fs = asImp_().fluidState_;
636 const Evaluation& muWater = fs.viscosity(waterPhaseIdx);
637 const auto& viscosityMultiplier =
638 PolymerModule::plyviscViscosityMultiplierTable(elemCtx, dofIdx, timeIdx);
639 const Evaluation viscosityMixture =
640 viscosityMultiplier.eval(polymerConcentration_, /*extrapolate=*/true) * muWater;
641
642 // Do the Todd-Longstaff mixing
643 const Scalar plymixparToddLongstaff = PolymerModule::plymixparToddLongstaff(elemCtx, dofIdx, timeIdx);
644 const Evaluation viscosityPolymer = viscosityMultiplier.eval(cmax, /*extrapolate=*/true) * muWater;
645 const Evaluation viscosityPolymerEffective =
646 pow(viscosityMixture, plymixparToddLongstaff) * pow(viscosityPolymer, 1.0 - plymixparToddLongstaff);
647 const Evaluation viscosityWaterEffective =
648 pow(viscosityMixture, plymixparToddLongstaff) * pow(muWater, 1.0 - plymixparToddLongstaff);
649
650 const Evaluation cbar = polymerConcentration_ / cmax;
651 // waterViscosity / effectiveWaterViscosity
652 waterViscosityCorrection_ = muWater * ((1.0 - cbar) / viscosityWaterEffective + cbar / viscosityPolymerEffective);
653 // effectiveWaterViscosity / effectivePolymerViscosity
654 polymerViscosityCorrection_ = (muWater / waterViscosityCorrection_) / viscosityPolymerEffective;
655 }
656 else { // based on PLYVMH
657 const auto& plyvmhCoefficients = PolymerModule::plyvmhCoefficients(elemCtx, dofIdx, timeIdx);
658 const Scalar k_mh = plyvmhCoefficients.k_mh;
659 const Scalar a_mh = plyvmhCoefficients.a_mh;
660 const Scalar gamma = plyvmhCoefficients.gamma;
661 const Scalar kappa = plyvmhCoefficients.kappa;
662
663 // viscosity model based on Mark-Houwink equation and Huggins equation
664 // 1000 is a emperical constant, most likely related to unit conversion
665 const Evaluation intrinsicViscosity = k_mh * pow(polymerMoleWeight_ * 1000., a_mh);
666 const Evaluation x = polymerConcentration_ * intrinsicViscosity;
667 waterViscosityCorrection_ = 1.0 / (1.0 + gamma * (x + kappa * x * x));
668 polymerViscosityCorrection_ = 1.0;
669 }
670
671 // adjust water mobility
672 asImp_().mobility_[waterPhaseIdx] *= waterViscosityCorrection_ / resistanceFactor;
673
674 // update rock properties
675 polymerDeadPoreVolume_ = PolymerModule::plyrockDeadPoreVolume(elemCtx, dofIdx, timeIdx);
676 polymerRockDensity_ = PolymerModule::plyrockRockDensityFactor(elemCtx, dofIdx, timeIdx);
677 }
678
679 const Evaluation& polymerConcentration() const
680 { return polymerConcentration_; }
681
682 const Evaluation& polymerMoleWeight() const
683 {
684 if constexpr (enablePolymerMolarWeight) {
685 return polymerMoleWeight_;
686 }
687 else {
688 throw std::logic_error("polymerMoleWeight() is called but polymer milecular weight is disabled");
689 }
690 }
691
693 { return polymerDeadPoreVolume_; }
694
695 const Evaluation& polymerAdsorption() const
696 { return polymerAdsorption_; }
697
698 Scalar polymerRockDensity() const
699 { return polymerRockDensity_; }
700
701 // effectiveWaterViscosity / effectivePolymerViscosity
702 const Evaluation& polymerViscosityCorrection() const
703 { return polymerViscosityCorrection_; }
704
705 // waterViscosity / effectiveWaterViscosity
706 const Evaluation& waterViscosityCorrection() const
707 { return waterViscosityCorrection_; }
708
709protected:
710 Implementation& asImp_()
711 { return *static_cast<Implementation*>(this); }
712
714 // polymer molecular weight
721};
722
723template <class TypeTag>
725{
729
730public:
731 void polymerPropertiesUpdate_(const ElementContext&,
732 unsigned,
733 unsigned)
734 {}
735
736 const Evaluation& polymerMoleWeight() const
737 { throw std::logic_error("polymerMoleWeight() called but polymer molecular weight is disabled"); }
738
739 const Evaluation& polymerConcentration() const
740 { throw std::runtime_error("polymerConcentration() called but polymers are disabled"); }
741
742 const Evaluation& polymerDeadPoreVolume() const
743 { throw std::runtime_error("polymerDeadPoreVolume() called but polymers are disabled"); }
744
745 const Evaluation& polymerAdsorption() const
746 { throw std::runtime_error("polymerAdsorption() called but polymers are disabled"); }
747
748 const Evaluation& polymerRockDensity() const
749 { throw std::runtime_error("polymerRockDensity() called but polymers are disabled"); }
750
751 const Evaluation& polymerViscosityCorrection() const
752 { throw std::runtime_error("polymerViscosityCorrection() called but polymers are disabled"); }
753
754 const Evaluation& waterViscosityCorrection() const
755 { throw std::runtime_error("waterViscosityCorrection() called but polymers are disabled"); }
756};
757
758template <class TypeTag, bool enablePolymerV>
760
768template <class TypeTag>
769class BlackOilPolymerExtensiveQuantities<TypeTag, /*enablePolymerV=*/true>
770{
772
780
781 static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
782 static constexpr int dimWorld = GridView::dimensionworld;
783 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
784
785 using Toolbox = MathToolbox<Evaluation>;
787 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
788 using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
789
790public:
797 template <class Dummy = bool> // we need to make this method a template to avoid
798 // compiler errors if it is not instantiated!
799 void updateShearMultipliersPerm(const ElementContext&,
800 unsigned,
801 unsigned)
802 {
803 throw std::runtime_error("The extension of the blackoil model for polymers is not yet "
804 "implemented for problems specified using permeabilities.");
805 }
806
813 template <class Dummy = bool> // we need to make this method a template to avoid
814 // compiler errors if it is not instantiated!
815 void updateShearMultipliers(const ElementContext& elemCtx,
816 unsigned scvfIdx,
817 unsigned timeIdx)
818 {
819 waterShearFactor_ = 1.0;
820 polymerShearFactor_ = 1.0;
821
822 if (!PolymerModule::hasPlyshlog()) {
823 return;
824 }
825
826 const ExtensiveQuantities& extQuants = asImp_();
827 const unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
828 const unsigned interiorDofIdx = extQuants.interiorIndex();
829 const unsigned exteriorDofIdx = extQuants.exteriorIndex();
830 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
831 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
832 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
833
834 // compute water velocity from flux
835 const Evaluation poroAvg = intQuantsIn.porosity() * 0.5 + intQuantsEx.porosity() * 0.5;
836 const unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvfIdx, timeIdx);
837 const Evaluation& Sw = up.fluidState().saturation(waterPhaseIdx);
838 const unsigned cellIdx = elemCtx.globalSpaceIndex(scvfIdx, timeIdx);
839 const auto& materialLawManager = elemCtx.problem().materialLawManager();
840 const auto& scaledDrainageInfo =
841 materialLawManager->oilWaterScaledEpsInfoDrainage(cellIdx);
842 const Scalar& Swcr = scaledDrainageInfo.Swcr;
843
844 // guard against zero porosity and no mobile water
845 const Evaluation denom = max(poroAvg * (Sw - Swcr), 1e-12);
846 Evaluation waterVolumeVelocity = extQuants.volumeFlux(waterPhaseIdx) / denom;
847
848 // if shrate is specified. Compute shrate based on the water velocity
849 if (PolymerModule::hasShrate()) {
850 const Evaluation& relWater = up.relativePermeability(waterPhaseIdx);
851 const Scalar trans = elemCtx.problem().transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
852 if (trans > 0.0) {
853 const Scalar faceArea = elemCtx.stencil(timeIdx).interiorFace(scvfIdx).area();
854 const auto dist = elemCtx.pos(interiorDofIdx, timeIdx) - elemCtx.pos(exteriorDofIdx, timeIdx);
855 // compute permeability from transmissibility.
856 const Scalar absPerm = trans / faceArea * dist.two_norm();
857 waterVolumeVelocity *=
858 PolymerModule::shrate(pvtnumRegionIdx) * sqrt(poroAvg * Sw / (relWater * absPerm));
859 assert(isfinite(waterVolumeVelocity));
860 }
861 }
862
863 // compute share factors for water and polymer
864 waterShearFactor_ =
865 PolymerModule::computeShearFactor(up.polymerConcentration(),
866 pvtnumRegionIdx,
867 waterVolumeVelocity);
868 polymerShearFactor_ =
869 PolymerModule::computeShearFactor(up.polymerConcentration(),
870 pvtnumRegionIdx,
871 waterVolumeVelocity * up.polymerViscosityCorrection());
872 }
873
874 const Evaluation& polymerShearFactor() const
875 { return polymerShearFactor_; }
876
877 const Evaluation& waterShearFactor() const
878 { return waterShearFactor_; }
879
880private:
881 Implementation& asImp_()
882 { return *static_cast<Implementation*>(this); }
883
884 Evaluation polymerShearFactor_;
885 Evaluation waterShearFactor_;
886};
887
888template <class TypeTag>
890{
893
894public:
895 void updateShearMultipliers(const ElementContext&,
896 unsigned,
897 unsigned)
898 {}
899
900 void updateShearMultipliersPerm(const ElementContext&,
901 unsigned,
902 unsigned)
903 {}
904
905 const Evaluation& polymerShearFactor() const
906 { throw std::runtime_error("polymerShearFactor() called but polymers are disabled"); }
907
908 const Evaluation& waterShearFactor() const
909 { throw std::runtime_error("waterShearFactor() called but polymers are disabled"); }
910};
911
912} // namespace Opm
913
914#endif
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:905
void updateShearMultipliersPerm(const ElementContext &, unsigned, unsigned)
Definition: blackoilpolymermodules.hh:900
void updateShearMultipliers(const ElementContext &, unsigned, unsigned)
Definition: blackoilpolymermodules.hh:895
const Evaluation & waterShearFactor() const
Definition: blackoilpolymermodules.hh:908
const Evaluation & polymerShearFactor() const
Definition: blackoilpolymermodules.hh:874
void updateShearMultipliersPerm(const ElementContext &, unsigned, unsigned)
Method which calculates the shear factor based on flow velocity.
Definition: blackoilpolymermodules.hh:799
const Evaluation & waterShearFactor() const
Definition: blackoilpolymermodules.hh:877
void updateShearMultipliers(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Method which calculates the shear factor based on flow velocity.
Definition: blackoilpolymermodules.hh:815
Provides the polymer specific extensive quantities to the generic black-oil module's extensive quanti...
Definition: blackoilpolymermodules.hh:759
const Evaluation & polymerConcentration() const
Definition: blackoilpolymermodules.hh:739
const Evaluation & waterViscosityCorrection() const
Definition: blackoilpolymermodules.hh:754
const Evaluation & polymerRockDensity() const
Definition: blackoilpolymermodules.hh:748
const Evaluation & polymerDeadPoreVolume() const
Definition: blackoilpolymermodules.hh:742
const Evaluation & polymerMoleWeight() const
Definition: blackoilpolymermodules.hh:736
void polymerPropertiesUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilpolymermodules.hh:731
const Evaluation & polymerAdsorption() const
Definition: blackoilpolymermodules.hh:745
const Evaluation & polymerViscosityCorrection() const
Definition: blackoilpolymermodules.hh:751
Evaluation polymerViscosityCorrection_
Definition: blackoilpolymermodules.hh:719
Scalar polymerDeadPoreVolume_
Definition: blackoilpolymermodules.hh:716
void polymerPropertiesUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the intensive properties needed to handle polymers from the primary variables.
Definition: blackoilpolymermodules.hh:603
Implementation & asImp_()
Definition: blackoilpolymermodules.hh:710
Scalar polymerRockDensity_
Definition: blackoilpolymermodules.hh:717
const Evaluation & polymerMoleWeight() const
Definition: blackoilpolymermodules.hh:682
Evaluation polymerConcentration_
Definition: blackoilpolymermodules.hh:713
Scalar polymerDeadPoreVolume() const
Definition: blackoilpolymermodules.hh:692
Evaluation polymerMoleWeight_
Definition: blackoilpolymermodules.hh:715
const Evaluation & waterViscosityCorrection() const
Definition: blackoilpolymermodules.hh:706
const Evaluation & polymerAdsorption() const
Definition: blackoilpolymermodules.hh:695
const Evaluation & polymerConcentration() const
Definition: blackoilpolymermodules.hh:679
const Evaluation & polymerViscosityCorrection() const
Definition: blackoilpolymermodules.hh:702
Evaluation waterViscosityCorrection_
Definition: blackoilpolymermodules.hh:720
Scalar polymerRockDensity() const
Definition: blackoilpolymermodules.hh:698
Evaluation polymerAdsorption_
Definition: blackoilpolymermodules.hh:718
Provides the volumetric quantities required for the equations needed by the polymers extension of the...
Definition: blackoilpolymermodules.hh:567
Contains the high level supplements required to extend the black oil model by polymer.
Definition: blackoilpolymermodules.hh:64
static const TabulatedFunction & plyviscViscosityMultiplierTable(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:418
static std::string primaryVarName(unsigned pvIdx)
Definition: blackoilpolymermodules.hh:181
static const TabulatedFunction & plyadsAdsorbedPolymer(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:409
static TabulatedTwoDFunction & getPlymwinjTable(const int tableNumber)
get the PLYMWINJ table
Definition: blackoilpolymermodules.hh:105
static Scalar plyrockAdsorbtionIndex(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:392
static void serializeEntity(const Model &model, std::ostream &outstream, const DofEntity &dof)
Definition: blackoilpolymermodules.hh:341
static Scalar eqWeight(unsigned eqIdx)
Definition: blackoilpolymermodules.hh:228
static void registerParameters()
Register all run-time parameters for the black-oil polymer module.
Definition: blackoilpolymermodules.hh:148
static bool eqApplies(unsigned eqIdx)
Definition: blackoilpolymermodules.hh:201
static Scalar primaryVarWeight(unsigned pvIdx)
Definition: blackoilpolymermodules.hh:193
static BlackOilPolymerParams< Scalar >::SkprpolyTable & getSkprpolyTable(const int tableNumber)
get the SKPRPOLY table
Definition: blackoilpolymermodules.hh:134
static bool hasShrate()
Definition: blackoilpolymermodules.hh:462
static Scalar plyrockRockDensityFactor(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:384
static Scalar plyrockDeadPoreVolume(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:368
static const TabulatedFunction & plyviscViscosityMultiplierTable(unsigned pvtnumRegionIdx)
Definition: blackoilpolymermodules.hh:428
static void registerOutputModules(Model &model, Simulator &simulator)
Register all polymer specific VTK and ECL output modules.
Definition: blackoilpolymermodules.hh:158
static Scalar plymaxMaxConcentration(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:431
static Scalar computeUpdateError(const PrimaryVariables &, const EqVector &)
Return how much a Newton-Raphson update is considered an error.
Definition: blackoilpolymermodules.hh:331
Scalar molarMass() const
Definition: blackoilpolymermodules.hh:555
static const BlackOilPolymerParams< Scalar >::PlyvmhCoefficients & plyvmhCoefficients(const ElementContext &elemCtx, const unsigned scvIdx, const unsigned timeIdx)
Definition: blackoilpolymermodules.hh:450
static bool primaryVarApplies(unsigned pvIdx)
Definition: blackoilpolymermodules.hh:166
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilpolymermodules.hh:238
static bool hasPlyshlog()
Definition: blackoilpolymermodules.hh:459
static Scalar shrate(unsigned pvtnumRegionIdx)
Definition: blackoilpolymermodules.hh:465
static std::string eqName(unsigned eqIdx)
Definition: blackoilpolymermodules.hh:216
static void setParams(BlackOilPolymerParams< Scalar > &&params)
Set parameters.
Definition: blackoilpolymermodules.hh:97
static TabulatedTwoDFunction & getSkprwatTable(const int tableNumber)
get the SKPRWAT table
Definition: blackoilpolymermodules.hh:119
static void deserializeEntity(Model &model, std::istream &instream, const DofEntity &dof)
Definition: blackoilpolymermodules.hh:352
static Scalar plymixparToddLongstaff(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:440
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:277
static Scalar plyrockMaxAdsorbtion(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:400
static Scalar plyrockResidualResistanceFactor(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:376
static Evaluation computeShearFactor(const Evaluation &polymerConcentration, unsigned pvtnumRegionIdx, const Evaluation &v0)
Computes the shear factor.
Definition: blackoilpolymermodules.hh:475
VTK output module for the black oil model's polymer related quantities.
Definition: vtkblackoilpolymermodule.hpp:61
static void registerParameters()
Register all run-time parameters for the multi-phase VTK output module.
Definition: vtkblackoilpolymermodule.hpp:91
Definition: blackoilboundaryratevector.hh:39
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:233
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
The Opm property system, traits with inheritance.
Definition: blackoilpolymerparams.hpp:85
Definition: blackoilpolymerparams.hpp:92
Struct holding the parameters for the BlackOilPolymerModule class.
Definition: blackoilpolymerparams.hpp:45
IntervalTabulated2DFunction< Scalar > TabulatedTwoDFunction
Definition: blackoilpolymerparams.hpp:47
Tabulated1DFunction< Scalar > TabulatedFunction
Definition: blackoilpolymerparams.hpp:46