blackoilsolventmodules.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_SOLVENT_MODULE_HH
29#define EWOMS_BLACK_OIL_SOLVENT_MODULE_HH
30
31#include <dune/common/fvector.hh>
32
33#include <opm/common/Exceptions.hpp>
34
35#include <opm/material/common/MathToolbox.hpp>
36#include <opm/material/common/Valgrind.hpp>
37#include <opm/material/fluidsystems/blackoilpvt/SolventPvt.hpp>
38
41
44
46
48
49#include <algorithm>
50#include <array>
51#include <cassert>
52#include <cmath>
53#include <istream>
54#include <memory>
55#include <ostream>
56#include <stdexcept>
57#include <string>
58
59namespace Opm {
60
66template <class TypeTag, bool enableSolventV = getPropValue<TypeTag, Properties::EnableSolvent>()>
68{
81 using Toolbox = MathToolbox<Evaluation>;
82 using SolventPvt = typename BlackOilSolventParams<Scalar>::SolventPvt;
83 using Co2GasPvt = typename BlackOilSolventParams<Scalar>::Co2GasPvt;
84 using H2GasPvt = typename BlackOilSolventParams<Scalar>::H2GasPvt;
85 using BrineCo2Pvt = typename BlackOilSolventParams<Scalar>::BrineCo2Pvt;
86 using BrineH2Pvt = typename BlackOilSolventParams<Scalar>::BrineH2Pvt;
87
88 using TabulatedFunction = typename BlackOilSolventParams<Scalar>::TabulatedFunction;
89
90 static constexpr unsigned solventSaturationIdx = Indices::solventSaturationIdx;
91 static constexpr unsigned contiSolventEqIdx = Indices::contiSolventEqIdx;
92 static constexpr unsigned enableSolvent = enableSolventV;
93 static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
94 static constexpr unsigned numPhases = FluidSystem::numPhases;
95 static constexpr bool blackoilConserveSurfaceVolume =
96 getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();
97 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
98
99public:
102 { params_ = params; }
103
107 static void setSolventPvt(const SolventPvt& value)
108 { params_.solventPvt_ = value; }
109
110 static void setIsMiscible(const bool isMiscible)
111 { params_.isMiscible_ = isMiscible; }
112
116 static void registerParameters()
117 {
118 if constexpr (enableSolvent) {
120 }
121 }
122
126 static void registerOutputModules(Model& model,
127 Simulator& simulator)
128 {
129 if constexpr (enableSolvent) {
130 model.addOutputModule(std::make_unique<VtkBlackOilSolventModule<TypeTag>>(simulator));
131 }
132 }
133
134 static bool primaryVarApplies(unsigned pvIdx)
135 {
136 if constexpr (enableSolvent) {
137 return pvIdx == solventSaturationIdx;
138 }
139 else {
140 return false;
141 }
142 }
143
144 static std::string primaryVarName([[maybe_unused]] unsigned pvIdx)
145 {
146 assert(primaryVarApplies(pvIdx));
147
148 return "saturation_solvent";
149 }
150
151 static Scalar primaryVarWeight([[maybe_unused]] unsigned pvIdx)
152 {
153 assert(primaryVarApplies(pvIdx));
154
155 // TODO: it may be beneficial to chose this differently.
156 return static_cast<Scalar>(1.0);
157 }
158
159 static bool eqApplies(unsigned eqIdx)
160 {
161 if constexpr (enableSolvent) {
162 return eqIdx == contiSolventEqIdx;
163 }
164 else {
165 return false;
166 }
167 }
168
169 static std::string eqName([[maybe_unused]] unsigned eqIdx)
170 {
171 assert(eqApplies(eqIdx));
172
173 return "conti^solvent";
174 }
175
176 static Scalar eqWeight([[maybe_unused]] unsigned eqIdx)
177 {
178 assert(eqApplies(eqIdx));
179
180 // TODO: it may be beneficial to chose this differently.
181 return static_cast<Scalar>(1.0);
182 }
183
184 template <class LhsEval>
185 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
186 const IntensiveQuantities& intQuants)
187 {
188 if constexpr (enableSolvent) {
189 if constexpr (blackoilConserveSurfaceVolume) {
190 storage[contiSolventEqIdx] +=
191 Toolbox::template decay<LhsEval>(intQuants.porosity()) *
192 Toolbox::template decay<LhsEval>(intQuants.solventSaturation()) *
193 Toolbox::template decay<LhsEval>(intQuants.solventInverseFormationVolumeFactor());
194 if (isSolubleInWater()) {
195 storage[contiSolventEqIdx] +=
196 Toolbox::template decay<LhsEval>(intQuants.porosity()) *
197 Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(waterPhaseIdx)) *
198 Toolbox::template decay<LhsEval>(intQuants.fluidState().invB(waterPhaseIdx)) *
199 Toolbox::template decay<LhsEval>(intQuants.rsSolw());
200 }
201 }
202 else {
203 storage[contiSolventEqIdx] +=
204 Toolbox::template decay<LhsEval>(intQuants.porosity()) *
205 Toolbox::template decay<LhsEval>(intQuants.solventSaturation()) *
206 Toolbox::template decay<LhsEval>(intQuants.solventDensity());
207 if (isSolubleInWater()) {
208 storage[contiSolventEqIdx] +=
209 Toolbox::template decay<LhsEval>(intQuants.porosity()) *
210 Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(waterPhaseIdx)) *
211 Toolbox::template decay<LhsEval>(intQuants.fluidState().density(waterPhaseIdx)) *
212 Toolbox::template decay<LhsEval>(intQuants.rsSolw());
213 }
214 }
215 }
216 }
217
218 static void computeFlux([[maybe_unused]] RateVector& flux,
219 [[maybe_unused]] const ElementContext& elemCtx,
220 [[maybe_unused]] unsigned scvfIdx,
221 [[maybe_unused]] unsigned timeIdx)
222
223 {
224 if constexpr (enableSolvent) {
225 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
226
227 const unsigned upIdx = extQuants.solventUpstreamIndex();
228 const unsigned inIdx = extQuants.interiorIndex();
229 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
230
231 if constexpr (blackoilConserveSurfaceVolume) {
232 if (upIdx == inIdx) {
233 flux[contiSolventEqIdx] = extQuants.solventVolumeFlux() *
234 up.solventInverseFormationVolumeFactor();
235 }
236 else {
237 flux[contiSolventEqIdx] = extQuants.solventVolumeFlux() *
238 decay<Scalar>(up.solventInverseFormationVolumeFactor());
239 }
240
241
242 if (isSolubleInWater()) {
243 if (upIdx == inIdx) {
244 flux[contiSolventEqIdx] +=
245 extQuants.volumeFlux(waterPhaseIdx) *
246 up.fluidState().invB(waterPhaseIdx) *
247 up.rsSolw();
248 }
249 else {
250 flux[contiSolventEqIdx] +=
251 extQuants.volumeFlux(waterPhaseIdx) *
252 decay<Scalar>(up.fluidState().invB(waterPhaseIdx)) *
253 decay<Scalar>(up.rsSolw());
254 }
255 }
256 }
257 else {
258 if (upIdx == inIdx) {
259 flux[contiSolventEqIdx] = extQuants.solventVolumeFlux() * up.solventDensity();
260 }
261 else {
262 flux[contiSolventEqIdx] = extQuants.solventVolumeFlux() * decay<Scalar>(up.solventDensity());
263 }
264
265 if (isSolubleInWater()) {
266 if (upIdx == inIdx) {
267 flux[contiSolventEqIdx] +=
268 extQuants.volumeFlux(waterPhaseIdx) *
269 up.fluidState().density(waterPhaseIdx) *
270 up.rsSolw();
271 }
272 else {
273 flux[contiSolventEqIdx] +=
274 extQuants.volumeFlux(waterPhaseIdx) *
275 decay<Scalar>(up.fluidState().density(waterPhaseIdx)) *
276 decay<Scalar>(up.rsSolw());
277 }
278 }
279 }
280 }
281 }
282
286 static void assignPrimaryVars(PrimaryVariables& priVars,
287 Scalar solventSaturation,
288 Scalar solventRsw)
289 {
290 if constexpr (!enableSolvent) {
291 priVars.setPrimaryVarsMeaningSolvent(PrimaryVariables::SolventMeaning::Disabled);
292 return;
293 }
294 // Determine the meaning of the solvent primary variables
295 if (solventSaturation > 0 || !isSolubleInWater()) {
296 priVars.setPrimaryVarsMeaningSolvent(PrimaryVariables::SolventMeaning::Ss);
297 priVars[solventSaturationIdx] = solventSaturation;
298 } else {
299 priVars.setPrimaryVarsMeaningSolvent(PrimaryVariables::SolventMeaning::Rsolw);
300 priVars[solventSaturationIdx] = solventRsw;
301 }
302 }
303
307 static void updatePrimaryVars(PrimaryVariables& newPv,
308 const PrimaryVariables& oldPv,
309 const EqVector& delta)
310 {
311 if constexpr (enableSolvent) {
312 // do a plain unchopped Newton update
313 newPv[solventSaturationIdx] = oldPv[solventSaturationIdx] - delta[solventSaturationIdx];
314 }
315 }
316
320 static Scalar computeUpdateError(const PrimaryVariables&,
321 const EqVector&)
322 {
323 // do not consider consider the cange of solvent primary variables for
324 // convergence
325 // TODO: maybe this should be changed
326 return static_cast<Scalar>(0.0);
327 }
328
332 static Scalar computeResidualError(const EqVector& resid)
333 {
334 // do not weight the residual of solvents when it comes to convergence
335 return std::abs(Toolbox::scalarValue(resid[contiSolventEqIdx]));
336 }
337
338 template <class DofEntity>
339 static void serializeEntity(const Model& model, std::ostream& outstream, const DofEntity& dof)
340 {
341 if constexpr (enableSolvent) {
342 const unsigned dofIdx = model.dofMapper().index(dof);
343
344 const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
345 outstream << priVars[solventSaturationIdx];
346 }
347 }
348
349 template <class DofEntity>
350 static void deserializeEntity(Model& model, std::istream& instream, const DofEntity& dof)
351 {
352 if constexpr (enableSolvent) {
353 const unsigned dofIdx = model.dofMapper().index(dof);
354
355 PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
356 PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
357
358 instream >> priVars0[solventSaturationIdx];
359
360 // set the primary variables for the beginning of the current time step.
361 priVars1 = priVars0[solventSaturationIdx];
362 }
363 }
364
365 static const SolventPvt& solventPvt()
366 { return params_.solventPvt_; }
367
368 static const Co2GasPvt& co2GasPvt()
369 { return params_.co2GasPvt_; }
370
371 static const H2GasPvt& h2GasPvt()
372 { return params_.h2GasPvt_; }
373
374 static const BrineCo2Pvt& brineCo2Pvt()
375 { return params_.brineCo2Pvt_; }
376
377 static const BrineH2Pvt& brineH2Pvt()
378 { return params_.brineH2Pvt_; }
379
380 template <class ElemContext>
381 static const TabulatedFunction& ssfnKrg(const ElemContext& elemCtx,
382 unsigned scvIdx,
383 unsigned timeIdx)
384 {
385 const unsigned satnumRegionIdx =
386 elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
387 return params_.ssfnKrg_[satnumRegionIdx];
388 }
389
390 template <class ElemContext>
391 static const TabulatedFunction& ssfnKrs(const ElemContext& elemCtx,
392 unsigned scvIdx,
393 unsigned timeIdx)
394 {
395 const unsigned satnumRegionIdx =
396 elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
397 return params_.ssfnKrs_[satnumRegionIdx];
398 }
399
400 template <class ElemContext>
401 static const TabulatedFunction& sof2Krn(const ElemContext& elemCtx,
402 unsigned scvIdx,
403 unsigned timeIdx)
404 {
405 const unsigned satnumRegionIdx =
406 elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
407 return params_.sof2Krn_[satnumRegionIdx];
408 }
409
410 template <class ElemContext>
411 static const TabulatedFunction& misc(const ElemContext& elemCtx,
412 unsigned scvIdx,
413 unsigned timeIdx)
414 {
415 const unsigned miscnumRegionIdx =
416 elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
417 return params_.misc_[miscnumRegionIdx];
418 }
419
420 template <class ElemContext>
421 static const TabulatedFunction& pmisc(const ElemContext& elemCtx,
422 unsigned scvIdx,
423 unsigned timeIdx)
424 {
425 const unsigned miscnumRegionIdx =
426 elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
427 return params_.pmisc_[miscnumRegionIdx];
428 }
429
430 template <class ElemContext>
431 static const TabulatedFunction& msfnKrsg(const ElemContext& elemCtx,
432 unsigned scvIdx,
433 unsigned timeIdx)
434 {
435 const unsigned satnumRegionIdx =
436 elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
437 return params_.msfnKrsg_[satnumRegionIdx];
438 }
439
440 template <class ElemContext>
441 static const TabulatedFunction& msfnKro(const ElemContext& elemCtx,
442 unsigned scvIdx,
443 unsigned timeIdx)
444 {
445 const unsigned satnumRegionIdx =
446 elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
447 return params_.msfnKro_[satnumRegionIdx];
448 }
449
450 template <class ElemContext>
451 static const TabulatedFunction& sorwmis(const ElemContext& elemCtx,
452 unsigned scvIdx,
453 unsigned timeIdx)
454 {
455 const unsigned miscnumRegionIdx =
456 elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
457 return params_.sorwmis_[miscnumRegionIdx];
458 }
459
460 template <class ElemContext>
461 static const TabulatedFunction& sgcwmis(const ElemContext& elemCtx,
462 unsigned scvIdx,
463 unsigned timeIdx)
464 {
465 const unsigned miscnumRegionIdx =
466 elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
467 return params_.sgcwmis_[miscnumRegionIdx];
468 }
469
470 template <class ElemContext>
471 static const TabulatedFunction& tlPMixTable(const ElemContext& elemCtx,
472 unsigned scvIdx,
473 unsigned timeIdx)
474 {
475 const unsigned miscnumRegionIdx =
476 elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
477 return params_.tlPMixTable_[miscnumRegionIdx];
478 }
479
480 template <class ElemContext>
481 static Scalar tlMixParamViscosity(const ElemContext& elemCtx,
482 unsigned scvIdx,
483 unsigned timeIdx)
484 {
485 const unsigned miscnumRegionIdx =
486 elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
487 return params_.tlMixParamViscosity_[miscnumRegionIdx];
488 }
489
490
491 template <class ElemContext>
492 static Scalar tlMixParamDensity(const ElemContext& elemCtx,
493 unsigned scvIdx,
494 unsigned timeIdx)
495 {
496 const unsigned miscnumRegionIdx =
497 elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
498 return params_.tlMixParamDensity_[miscnumRegionIdx];
499 }
500
501 static bool isMiscible()
502 { return params_.isMiscible_; }
503
504 template <class Value>
505 static Value solubilityLimit(unsigned pvtIdx, const Value& temperature,
506 const Value& pressure, const Value& saltConcentration)
507 {
508 if (!isSolubleInWater()) {
509 return 0.0;
510 }
511
512 assert(isCO2Sol() || isH2Sol());
513 if (isCO2Sol()) {
514 return brineCo2Pvt().saturatedGasDissolutionFactor(pvtIdx, temperature,
515 pressure, saltConcentration);
516 }
517 else {
518 return brineH2Pvt().saturatedGasDissolutionFactor(pvtIdx, temperature,
519 pressure, saltConcentration);
520 }
521 }
522
523 static bool isSolubleInWater()
524 { return params_.rsSolw_active_; }
525
526 static bool isCO2Sol()
527 { return params_.co2sol_; }
528
529 static bool isH2Sol()
530 { return params_.h2sol_; }
531
532private:
533 static BlackOilSolventParams<Scalar> params_; // the krg(Fs) column of the SSFN table
534};
535
536template <class TypeTag, bool enableSolventV>
537BlackOilSolventParams<typename BlackOilSolventModule<TypeTag, enableSolventV>::Scalar>
538BlackOilSolventModule<TypeTag, enableSolventV>::params_;
539
540template <class TypeTag, bool enableSolventV>
542
550template <class TypeTag>
551class BlackOilSolventIntensiveQuantities<TypeTag, /*enableSolventV=*/true>
552{
554
563
565
566 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
567 static constexpr int solventSaturationIdx = Indices::solventSaturationIdx;
568 static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
569 static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
570 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
571 static constexpr double cutOff = 1e-12;
572
573public:
580 void solventPreSatFuncUpdate_(const ElementContext& elemCtx,
581 unsigned dofIdx,
582 unsigned timeIdx)
583 {
584 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
585 this->solventPreSatFuncUpdate_(priVars, timeIdx, elemCtx.linearizationType());
586 }
587
588 void solventPreSatFuncUpdate_(const PrimaryVariables& priVars,
589 const unsigned timeIdx,
590 const LinearizationType linearizationType)
591 {
592 auto& fs = asImp_().fluidState_;
593 solventSaturation_ = 0.0;
594 if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
595 solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx, linearizationType);
596 }
597
598 hydrocarbonSaturation_ = fs.saturation(gasPhaseIdx);
599
600 // apply a cut-off. Don't waste calculations if no solvent
601 if (solventSaturation().value() < cutOff) {
602 return;
603 }
604
605 // make the saturation of the gas phase which is used by the saturation functions
606 // the sum of the solvent "saturation" and the saturation the hydrocarbon gas.
607 fs.setSaturation(gasPhaseIdx, hydrocarbonSaturation_ + solventSaturation_);
608 }
609
617 void solventPostSatFuncUpdate_(const ElementContext& elemCtx,
618 unsigned dofIdx,
619 unsigned timeIdx)
620 {
621 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
622 this->solventPostSatFuncUpdate_(elemCtx.problem(), priVars, elemCtx.globalSpaceIndex(dofIdx, timeIdx), timeIdx, elemCtx.linearizationType());
623 }
624
625private:
626 // This class is a private implementation detail.
627 template <class Problem>
628 struct ProblemAndCellIndexOnlyContext
629 {
630 const Problem& problem_;
631 unsigned int index_;
632 const Problem& problem() const { return problem_; }
633 unsigned int globalSpaceIndex([[maybe_unused]] const unsigned int spaceIdx,
634 [[maybe_unused]] const unsigned int timeIdx) const
635 {
636 return index_;
637 }
638 };
639
640public:
641 void solventPostSatFuncUpdate_(const Problem& problem,
642 const PrimaryVariables& priVars,
643 const unsigned globalSpaceIdx,
644 const unsigned timeIdx,
645 const LinearizationType linearizationType)
646 {
647 // revert the gas "saturation" of the fluid state back to the saturation of the
648 // hydrocarbon gas.
649 auto& fs = asImp_().fluidState_;
650 fs.setSaturation(gasPhaseIdx, hydrocarbonSaturation_);
651
652 // update rsSolw. This needs to be done after the pressure is defined in the fluid state.
653 rsSolw_ = 0.0;
654 if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
655 rsSolw_ = SolventModule::solubilityLimit(asImp_().pvtRegionIndex(),
656 fs.temperature(waterPhaseIdx),
657 fs.pressure(waterPhaseIdx),
658 fs.saltConcentration());
659 } else if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Rsolw) {
660 rsSolw_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx, linearizationType);
661 }
662
663 solventMobility_ = 0.0;
664
665 // apply a cut-off. Don't waste calculations if no solvent
666 if (solventSaturation().value() < cutOff) {
667 return;
668 }
669
670 ProblemAndCellIndexOnlyContext<Problem> elemCtx{problem, globalSpaceIdx};
671 unsigned dofIdx = 0; // Dummy
672
673 // Pressure effects on capillary pressure miscibility
674 if (SolventModule::isMiscible()) {
675 const Evaluation& p =
676 FluidSystem::phaseIsActive(oilPhaseIdx)
677 ? fs.pressure(oilPhaseIdx)
678 : fs.pressure(gasPhaseIdx);
679 const Evaluation pmisc =
680 SolventModule::pmisc(elemCtx, dofIdx, timeIdx).eval(p, /*extrapolate=*/true);
681 const Evaluation& pgImisc = fs.pressure(gasPhaseIdx);
682
683 // compute capillary pressure for miscible fluid
684 Evaluation pgMisc = 0.0;
685 std::array<Evaluation, numPhases> pC;
686 const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, timeIdx);
687 MaterialLaw::capillaryPressures(pC, materialParams, fs);
688
689 // oil is the reference phase for pressure
690 if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pg) {
691 pgMisc = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx,
692 linearizationType);
693 }
694 else {
695 const Evaluation& po =
696 priVars.makeEvaluation(Indices::pressureSwitchIdx,
697 timeIdx, linearizationType);
698 pgMisc = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
699 }
700
701 fs.setPressure(gasPhaseIdx, pmisc * pgMisc + (1.0 - pmisc) * pgImisc);
702 }
703
704 const Evaluation gasSolventSat = hydrocarbonSaturation_ + solventSaturation_;
705
706 if (gasSolventSat.value() < cutOff) { // avoid division by zero
707 return;
708 }
709
710 const Evaluation Fhydgas = hydrocarbonSaturation_ / gasSolventSat;
711 const Evaluation Fsolgas = solventSaturation_ / gasSolventSat;
712
713 // account for miscibility of oil and solvent
714 if (SolventModule::isMiscible() && FluidSystem::phaseIsActive(oilPhaseIdx)) {
715 const auto& misc = SolventModule::misc(elemCtx, dofIdx, timeIdx);
716 const auto& pmisc = SolventModule::pmisc(elemCtx, dofIdx, timeIdx);
717 const Evaluation& p =
718 FluidSystem::phaseIsActive(oilPhaseIdx)
719 ? fs.pressure(oilPhaseIdx)
720 : fs.pressure(gasPhaseIdx);
721 const Evaluation miscibility =
722 misc.eval(Fsolgas, /*extrapolate=*/true) *
723 pmisc.eval(p, /*extrapolate=*/true);
724
725 // TODO adjust endpoints of sn and ssg
726 const unsigned cellIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
727 const auto& materialLawManager = elemCtx.problem().materialLawManager();
728 const auto& scaledDrainageInfo =
729 materialLawManager->oilWaterScaledEpsInfoDrainage(cellIdx);
730
731 const Scalar sogcr = scaledDrainageInfo.Sogcr;
732 Evaluation sor = sogcr;
733 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
734 const Evaluation& sw = fs.saturation(waterPhaseIdx);
735 const auto& sorwmis = SolventModule::sorwmis(elemCtx, dofIdx, timeIdx);
736 sor = miscibility *
737 sorwmis.eval(sw, /*extrapolate=*/true) + (1.0 - miscibility) * sogcr;
738 }
739 const Scalar sgcr = scaledDrainageInfo.Sgcr;
740 Evaluation sgc = sgcr;
741 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
742 const Evaluation& sw = fs.saturation(waterPhaseIdx);
743 const auto& sgcwmis = SolventModule::sgcwmis(elemCtx, dofIdx, timeIdx);
744 sgc = miscibility * sgcwmis.eval(sw, /*extrapolate=*/true) + (1.0 - miscibility) * sgcr;
745 }
746
747 Evaluation oilGasSolventSat = gasSolventSat;
748 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
749 oilGasSolventSat += fs.saturation(oilPhaseIdx);
750 }
751 const Evaluation zero = 0.0;
752 const Evaluation oilGasSolventEffSat = std::max(oilGasSolventSat - sor - sgc, zero);
753
754 Evaluation F_totalGas = 0.0;
755 if (oilGasSolventEffSat.value() > cutOff) {
756 const Evaluation gasSolventEffSat = std::max(gasSolventSat - sgc, zero);
757 F_totalGas = gasSolventEffSat / oilGasSolventEffSat;
758 }
759 const auto& msfnKro = SolventModule::msfnKro(elemCtx, dofIdx, timeIdx);
760 const auto& msfnKrsg = SolventModule::msfnKrsg(elemCtx, dofIdx, timeIdx);
761 const auto& sof2Krn = SolventModule::sof2Krn(elemCtx, dofIdx, timeIdx);
762
763 const Evaluation mkrgt = msfnKrsg.eval(F_totalGas, /*extrapolate=*/true) *
764 sof2Krn.eval(oilGasSolventSat, /*extrapolate=*/true);
765 const Evaluation mkro = msfnKro.eval(F_totalGas, /*extrapolate=*/true) *
766 sof2Krn.eval(oilGasSolventSat, /*extrapolate=*/true);
767
768 Evaluation& kro = asImp_().mobility_[oilPhaseIdx];
769 Evaluation& krg = asImp_().mobility_[gasPhaseIdx];
770
771 // combine immiscible and miscible part of the relperm
772 krg *= 1.0 - miscibility;
773 krg += miscibility * mkrgt;
774 kro *= 1.0 - miscibility;
775 kro += miscibility * mkro;
776 }
777
778 // compute the mobility of the solvent "phase" and modify the gas phase
779 const auto& ssfnKrg = SolventModule::ssfnKrg(elemCtx, dofIdx, timeIdx);
780 const auto& ssfnKrs = SolventModule::ssfnKrs(elemCtx, dofIdx, timeIdx);
781
782 Evaluation& krg = asImp_().mobility_[gasPhaseIdx];
783 solventMobility_ = krg * ssfnKrs.eval(Fsolgas, /*extrapolate=*/true);
784 krg *= ssfnKrg.eval(Fhydgas, /*extrapolate=*/true);
785 }
786
793 void solventPvtUpdate_(const ElementContext& elemCtx,
794 unsigned scvIdx,
795 unsigned timeIdx)
796 {
797 const auto& iq = asImp_();
798 const unsigned pvtRegionIdx = iq.pvtRegionIndex();
799 const Evaluation& T = iq.fluidState().temperature(gasPhaseIdx);
800 const Evaluation& p = iq.fluidState().pressure(gasPhaseIdx);
801
802 const Evaluation rv = 0.0;
803 const Evaluation rvw = 0.0;
804 if (SolventModule::isCO2Sol() || SolventModule::isH2Sol() ){
805 if (SolventModule::isCO2Sol()) {
806 const auto& co2gasPvt = SolventModule::co2GasPvt();
807 solventInvFormationVolumeFactor_ =
808 co2gasPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rv, rvw);
809 solventRefDensity_ = co2gasPvt.gasReferenceDensity(pvtRegionIdx);
810 solventViscosity_ = co2gasPvt.viscosity(pvtRegionIdx, T, p, rv, rvw);
811
812 const auto& brineCo2Pvt = SolventModule::brineCo2Pvt();
813 auto& fs = asImp_().fluidState_;
814 const auto& bw = brineCo2Pvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rsSolw());
815
816 const auto denw = bw * brineCo2Pvt.waterReferenceDensity(pvtRegionIdx) +
817 rsSolw() * bw * brineCo2Pvt.gasReferenceDensity(pvtRegionIdx);
818 fs.setDensity(waterPhaseIdx, denw);
819 fs.setInvB(waterPhaseIdx, bw);
820 Evaluation& mobw = asImp_().mobility_[waterPhaseIdx];
821 const auto& muWat = fs.viscosity(waterPhaseIdx);
822 const auto& muWatEff = brineCo2Pvt.viscosity(pvtRegionIdx, T, p, rsSolw());
823 mobw *= muWat / muWatEff;
824 } else {
825 const auto& h2gasPvt = SolventModule::h2GasPvt();
826 solventInvFormationVolumeFactor_ =
827 h2gasPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rv, rvw);
828 solventRefDensity_ = h2gasPvt.gasReferenceDensity(pvtRegionIdx);
829 solventViscosity_ = h2gasPvt.viscosity(pvtRegionIdx, T, p, rv, rvw);
830
831 const auto& brineH2Pvt = SolventModule::brineH2Pvt();
832 auto& fs = asImp_().fluidState_;
833 const auto& bw = brineH2Pvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rsSolw());
834
835 const auto denw = bw * brineH2Pvt.waterReferenceDensity(pvtRegionIdx) +
836 rsSolw() * bw * brineH2Pvt.gasReferenceDensity(pvtRegionIdx);
837 fs.setDensity(waterPhaseIdx, denw);
838 fs.setInvB(waterPhaseIdx, bw);
839 Evaluation& mobw = asImp_().mobility_[waterPhaseIdx];
840 const auto& muWat = fs.viscosity(waterPhaseIdx);
841 const auto& muWatEff = brineH2Pvt.viscosity(pvtRegionIdx, T, p, rsSolw());
842 mobw *= muWat / muWatEff;
843 }
844 } else {
845 const auto& solventPvt = SolventModule::solventPvt();
846 solventInvFormationVolumeFactor_ = solventPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
847 solventRefDensity_ = solventPvt.referenceDensity(pvtRegionIdx);
848 solventViscosity_ = solventPvt.viscosity(pvtRegionIdx, T, p);
849 }
850
851 solventDensity_ = solventInvFormationVolumeFactor_*solventRefDensity_;
852 effectiveProperties(elemCtx, scvIdx, timeIdx);
853
854 solventMobility_ /= solventViscosity_;
855 }
856
857 const Evaluation& solventSaturation() const
858 { return solventSaturation_; }
859
860 const Evaluation& rsSolw() const
861 { return rsSolw_; }
862
863 const Evaluation& solventDensity() const
864 { return solventDensity_; }
865
866 const Evaluation& solventViscosity() const
867 { return solventViscosity_; }
868
869 const Evaluation& solventMobility() const
870 { return solventMobility_; }
871
872 const Evaluation& solventInverseFormationVolumeFactor() const
873 { return solventInvFormationVolumeFactor_; }
874
875 // This could be stored pr pvtRegion instead
876 Scalar solventRefDensity() const
877 { return solventRefDensity_; }
878
879private:
880 // Computes the effective properties based on
881 // Todd-Longstaff mixing model.
882 void effectiveProperties(const ElementContext& elemCtx,
883 unsigned scvIdx,
884 unsigned timeIdx)
885 {
886 if (!SolventModule::isMiscible()) {
887 return;
888 }
889
890 // Don't waste calculations if no solvent
891 // Apply a cut-off for small and negative solvent saturations
892 if (solventSaturation() < cutOff) {
893 return;
894 }
895
896 // We plan to update the fluidstate with the effective
897 // properties
898 auto& fs = asImp_().fluidState_;
899
900 // Compute effective saturations
901 Evaluation oilEffSat = 0.0;
902 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
903 oilEffSat = fs.saturation(oilPhaseIdx);
904 }
905 Evaluation gasEffSat = fs.saturation(gasPhaseIdx);
906 Evaluation solventEffSat = solventSaturation();
907 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
908 const auto& sorwmis = SolventModule::sorwmis(elemCtx, scvIdx, timeIdx);
909 const auto& sgcwmis = SolventModule::sgcwmis(elemCtx, scvIdx, timeIdx);
910 const Evaluation zero = 0.0;
911 const Evaluation& sw = fs.saturation(waterPhaseIdx);
912 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
913 oilEffSat = std::max(oilEffSat - sorwmis.eval(sw, /*extrapolate=*/true), zero);
914 }
915 gasEffSat = std::max(gasEffSat - sgcwmis.eval(sw, /*extrapolate=*/true), zero);
916 solventEffSat = std::max(solventEffSat - sgcwmis.eval(sw, /*extrapolate=*/true), zero);
917 }
918 const Evaluation oilGasSolventEffSat = oilEffSat + gasEffSat + solventEffSat;
919 const Evaluation oilSolventEffSat = oilEffSat + solventEffSat;
920 const Evaluation solventGasEffSat = solventEffSat + gasEffSat;
921
922 // Compute effective viscosities
923
924 // Mixing parameter for viscosity
925 // The pressureMixingParameter represent the miscibility of the solvent while the mixingParameterViscosity the effect of the porous media.
926 // The pressureMixingParameter is not implemented in ecl100.
927 const Evaluation& p =
928 FluidSystem::phaseIsActive(oilPhaseIdx)
929 ? fs.pressure(oilPhaseIdx)
930 : fs.pressure(gasPhaseIdx);
931 // account for pressure effects
932 const auto& pmiscTable = SolventModule::pmisc(elemCtx, scvIdx, timeIdx);
933 const Evaluation pmisc = pmiscTable.eval(p, /*extrapolate=*/true);
934 const auto& tlPMixTable = SolventModule::tlPMixTable(elemCtx, scvIdx, timeIdx);
935 const Evaluation tlMixParamMu = SolventModule::tlMixParamViscosity(elemCtx, scvIdx, timeIdx) *
936 tlPMixTable.eval(p, /*extrapolate=*/true);
937
938 const Evaluation& muGas = fs.viscosity(gasPhaseIdx);
939 const Evaluation& muSolvent = solventViscosity_;
940
941 assert(muGas.value() > 0);
942 assert(muSolvent.value() > 0);
943 const Evaluation muGasPow = pow(muGas, 0.25);
944 const Evaluation muSolventPow = pow(muSolvent, 0.25);
945
946 Evaluation muMixSolventGas = muGas;
947 if (solventGasEffSat > cutOff) {
948 muMixSolventGas *= muSolvent / pow(((gasEffSat / solventGasEffSat) * muSolventPow) +
949 ((solventEffSat / solventGasEffSat) * muGasPow), 4.0);
950 }
951
952 Evaluation muOil = 1.0;
953 Evaluation muOilPow = 1.0;
954 Evaluation muMixOilSolvent = 1.0;
955 Evaluation muOilEff = 1.0;
956 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
957 muOil = fs.viscosity(oilPhaseIdx);
958 assert(muOil.value() > 0);
959 muOilPow = pow(muOil, 0.25);
960 muMixOilSolvent = muOil;
961 if (oilSolventEffSat > cutOff) {
962 muMixOilSolvent *= muSolvent / pow(((oilEffSat / oilSolventEffSat) * muSolventPow) +
963 ((solventEffSat / oilSolventEffSat) * muOilPow), 4.0);
964 }
965
966 muOilEff = pow(muOil,1.0 - tlMixParamMu) * pow(muMixOilSolvent, tlMixParamMu);
967 }
968 Evaluation muMixSolventGasOil = muOil;
969 if (oilGasSolventEffSat > cutOff) {
970 muMixSolventGasOil *= muSolvent * muGas / pow(((oilEffSat / oilGasSolventEffSat) *
971 muSolventPow * muGasPow) +
972 ((solventEffSat / oilGasSolventEffSat) *
973 muOilPow * muGasPow) +
974 ((gasEffSat / oilGasSolventEffSat) *
975 muSolventPow * muOilPow), 4.0);
976 }
977
978 Evaluation muGasEff = pow(muGas,1.0 - tlMixParamMu) * pow(muMixSolventGas, tlMixParamMu);
979 const Evaluation muSolventEff = pow(muSolvent,1.0 - tlMixParamMu) * pow(muMixSolventGasOil, tlMixParamMu);
980
981 // Compute effective densities
982 const Evaluation& rhoGas = fs.density(gasPhaseIdx);
983 Evaluation rhoOil = 0.0;
984 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
985 rhoOil = fs.density(oilPhaseIdx);
986 }
987
988 const Evaluation& rhoSolvent = solventDensity_;
989
990 // Mixing parameter for density
991 // The pressureMixingParameter represent the miscibility of the
992 // solvent while the mixingParameterDenisty the effect of the porous media.
993 // The pressureMixingParameter is not implemented in ecl100.
994 const Evaluation tlMixParamRho = SolventModule::tlMixParamDensity(elemCtx, scvIdx, timeIdx) *
995 tlPMixTable.eval(p, /*extrapolate=*/true);
996
997 // compute effective viscosities for density calculations. These have to
998 // be recomputed as a different mixing parameter may be used.
999 const Evaluation muOilEffPow = pow(pow(muOil, 1.0 - tlMixParamRho) *
1000 pow(muMixOilSolvent, tlMixParamRho), 0.25);
1001 const Evaluation muGasEffPow = pow(pow(muGas, 1.0 - tlMixParamRho) *
1002 pow(muMixSolventGas, tlMixParamRho), 0.25);
1003 const Evaluation muSolventEffPow = pow(pow(muSolvent, 1.0 - tlMixParamRho) *
1004 pow(muMixSolventGasOil, tlMixParamRho), 0.25);
1005
1006 const Evaluation oilGasEffSaturation = oilEffSat + gasEffSat;
1007 Evaluation sof = 0.0;
1008 Evaluation sgf = 0.0;
1009 if (oilGasEffSaturation.value() > cutOff) {
1010 sof = oilEffSat / oilGasEffSaturation;
1011 sgf = gasEffSat / oilGasEffSaturation;
1012 }
1013
1014 const Evaluation muSolventOilGasPow = muSolventPow * ((sgf * muOilPow) + (sof * muGasPow));
1015
1016 Evaluation rhoMixSolventGasOil = 0.0;
1017 if (oilGasSolventEffSat.value() > cutOff) {
1018 rhoMixSolventGasOil = (rhoOil * oilEffSat / oilGasSolventEffSat) +
1019 (rhoGas * gasEffSat / oilGasSolventEffSat) +
1020 (rhoSolvent * solventEffSat / oilGasSolventEffSat);
1021 }
1022
1023 Evaluation rhoGasEff = 0.0;
1024 if (std::abs(muSolventPow.value() - muGasPow.value()) < cutOff) {
1025 rhoGasEff = ((1.0 - tlMixParamRho) * rhoGas) +
1026 (tlMixParamRho * rhoMixSolventGasOil);
1027 }
1028 else {
1029 const Evaluation solventGasEffFraction = (muGasPow * (muSolventPow - muGasEffPow)) /
1030 (muGasEffPow * (muSolventPow - muGasPow));
1031 rhoGasEff = (rhoGas * solventGasEffFraction) + (rhoSolvent * (1.0 - solventGasEffFraction));
1032 }
1033
1034 Evaluation rhoSolventEff = 0.0;
1035 if (std::abs((muSolventOilGasPow.value() - (muOilPow.value() * muGasPow.value()))) < cutOff) {
1036 rhoSolventEff = ((1.0 - tlMixParamRho) * rhoSolvent) + (tlMixParamRho * rhoMixSolventGasOil);
1037 }
1038 else {
1039 const Evaluation sfraction_se = (muSolventOilGasPow -
1040 muOilPow * muGasPow * muSolventPow / muSolventEffPow) /
1041 (muSolventOilGasPow - (muOilPow * muGasPow));
1042 rhoSolventEff = (rhoSolvent * sfraction_se) +
1043 (rhoGas * sgf * (1.0 - sfraction_se)) +
1044 (rhoOil * sof * (1.0 - sfraction_se));
1045 }
1046
1047 // compute invB from densities.
1048 const unsigned pvtRegionIdx = asImp_().pvtRegionIndex();
1049 Evaluation bGasEff = rhoGasEff;
1050 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1051 bGasEff /= (FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) +
1052 FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) * fs.Rv());
1053 } else {
1054 bGasEff /= (FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx));
1055 }
1056 const Evaluation bSolventEff = rhoSolventEff / solventRefDensity();
1057
1058 // copy the unmodified invB factors
1059 const Evaluation bg = fs.invB(gasPhaseIdx);
1060 const Evaluation bs = solventInverseFormationVolumeFactor();
1061 const Evaluation bg_eff = pmisc * bGasEff + (1.0 - pmisc) * bg;
1062 const Evaluation bs_eff = pmisc * bSolventEff + (1.0 - pmisc) * bs;
1063
1064 // set the densities
1065 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1066 fs.setDensity(gasPhaseIdx,
1067 bg_eff *
1068 (FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) +
1069 FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx)*fs.Rv()));
1070 } else {
1071 fs.setDensity(gasPhaseIdx,
1072 bg_eff * FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx));
1073 }
1074 solventDensity_ = bs_eff * solventRefDensity();
1075
1076 // set the mobility
1077 Evaluation& mobg = asImp_().mobility_[gasPhaseIdx];
1078 muGasEff = bg_eff / (pmisc * bGasEff / muGasEff + (1.0 - pmisc) * bg / muGas);
1079 mobg *= muGas / muGasEff;
1080
1081 // Update viscosity of solvent
1082 solventViscosity_ = bs_eff / (pmisc * bSolventEff / muSolventEff +
1083 (1.0 - pmisc) * bs / muSolvent);
1084
1085 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1086 Evaluation rhoOilEff = 0.0;
1087 if (std::abs(muOilPow.value() - muSolventPow.value()) < cutOff) {
1088 rhoOilEff = ((1.0 - tlMixParamRho) * rhoOil) + (tlMixParamRho * rhoMixSolventGasOil);
1089 }
1090 else {
1091 const Evaluation solventOilEffFraction = (muOilPow * (muOilEffPow - muSolventPow)) /
1092 (muOilEffPow * (muOilPow - muSolventPow));
1093 rhoOilEff = (rhoOil * solventOilEffFraction) + (rhoSolvent * (1.0 - solventOilEffFraction));
1094 }
1095 const Evaluation bOilEff =
1096 rhoOilEff / (FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) +
1097 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) * fs.Rs());
1098 const Evaluation bo = fs.invB(oilPhaseIdx);
1099 const Evaluation bo_eff = pmisc * bOilEff + (1.0 - pmisc) * bo;
1100 fs.setDensity(oilPhaseIdx,
1101 bo_eff * (FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) +
1102 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) * fs.Rs()));
1103
1104 // keep the mu*b interpolation
1105 Evaluation& mobo = asImp_().mobility_[oilPhaseIdx];
1106 muOilEff = bo_eff / (pmisc * bOilEff / muOilEff + (1.0 - pmisc) * bo / muOil);
1107 mobo *= muOil / muOilEff;
1108 }
1109 }
1110
1111protected:
1112 Implementation& asImp_()
1113 { return *static_cast<Implementation*>(this); }
1114
1117 Evaluation rsSolw_;
1122
1124};
1125
1126template <class TypeTag>
1128{
1132
1133public:
1134 void solventPreSatFuncUpdate_(const ElementContext&,
1135 unsigned,
1136 unsigned)
1137 {}
1138
1139 void solventPostSatFuncUpdate_(const ElementContext&,
1140 unsigned,
1141 unsigned)
1142 {}
1143
1144 void solventPvtUpdate_(const ElementContext&,
1145 unsigned,
1146 unsigned)
1147 {}
1148
1149 const Evaluation& solventSaturation() const
1150 { throw std::runtime_error("solventSaturation() called but solvents are disabled"); }
1151
1152 const Evaluation& rsSolw() const
1153 { throw std::runtime_error("rsSolw() called but solvents are disabled"); }
1154
1155 const Evaluation& solventDensity() const
1156 { throw std::runtime_error("solventDensity() called but solvents are disabled"); }
1157
1158 const Evaluation& solventViscosity() const
1159 { throw std::runtime_error("solventViscosity() called but solvents are disabled"); }
1160
1161 const Evaluation& solventMobility() const
1162 { throw std::runtime_error("solventMobility() called but solvents are disabled"); }
1163
1164 const Evaluation& solventInverseFormationVolumeFactor() const
1165 { throw std::runtime_error("solventInverseFormationVolumeFactor() called but solvents are disabled"); }
1166
1167 Scalar solventRefDensity() const
1168 { throw std::runtime_error("solventRefDensity() called but solvents are disabled"); }
1169};
1170
1178template <class TypeTag, bool enableSolventV = getPropValue<TypeTag, Properties::EnableSolvent>()>
1180{
1182
1190
1191 using Toolbox = MathToolbox<Evaluation>;
1192
1193 static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
1194 static constexpr int dimWorld = GridView::dimensionworld;
1195
1196 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
1197 using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
1198
1199public:
1204 template <class Dummy = bool> // we need to make this method a template to avoid
1205 // compiler errors if it is not instantiated!
1206 void updateVolumeFluxPerm(const ElementContext& elemCtx,
1207 unsigned scvfIdx,
1208 unsigned timeIdx)
1209 {
1210 const auto& gradCalc = elemCtx.gradientCalculator();
1211 PressureCallback<TypeTag> pressureCallback(elemCtx);
1212
1213 const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
1214 const auto& faceNormal = scvf.normal();
1215
1216 const unsigned i = scvf.interiorIndex();
1217 const unsigned j = scvf.exteriorIndex();
1218
1219 // calculate the "raw" pressure gradient
1220 DimEvalVector solventPGrad;
1221 pressureCallback.setPhaseIndex(gasPhaseIdx);
1222 gradCalc.calculateGradient(solventPGrad,
1223 elemCtx,
1224 scvfIdx,
1225 pressureCallback);
1226 Valgrind::CheckDefined(solventPGrad);
1227
1228 // correct the pressure gradients by the gravitational acceleration
1229 if (Parameters::Get<Parameters::EnableGravity>()) {
1230 // estimate the gravitational acceleration at a given SCV face
1231 // using the arithmetic mean
1232 const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
1233 const auto& gEx = elemCtx.problem().gravity(elemCtx, j, timeIdx);
1234
1235 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
1236 const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
1237
1238 const auto& posIn = elemCtx.pos(i, timeIdx);
1239 const auto& posEx = elemCtx.pos(j, timeIdx);
1240 const auto& posFace = scvf.integrationPos();
1241
1242 // the distance between the centers of the control volumes
1243 DimVector distVecIn(posIn);
1244 DimVector distVecEx(posEx);
1245 DimVector distVecTotal(posEx);
1246
1247 distVecIn -= posFace;
1248 distVecEx -= posFace;
1249 distVecTotal -= posIn;
1250 const Scalar absDistTotalSquared = distVecTotal.two_norm2();
1251
1252 // calculate the hydrostatic pressure at the integration point of the face
1253 auto rhoIn = intQuantsIn.solventDensity();
1254 auto pStatIn = - rhoIn * (gIn * distVecIn);
1255
1256 // the quantities on the exterior side of the face do not influence the
1257 // result for the TPFA scheme, so they can be treated as scalar values.
1258 Scalar rhoEx = Toolbox::value(intQuantsEx.solventDensity());
1259 Scalar pStatEx = - rhoEx * (gEx * distVecEx);
1260
1261 // compute the hydrostatic gradient between the two control volumes (this
1262 // gradient exhibitis the same direction as the vector between the two
1263 // control volume centers and the length (pStaticExterior -
1264 // pStaticInterior)/distanceInteriorToExterior
1265 DimEvalVector f(distVecTotal);
1266 f *= (pStatEx - pStatIn) / absDistTotalSquared;
1267
1268 // calculate the final potential gradient
1269 for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
1270 solventPGrad[dimIdx] += f[dimIdx];
1271
1272 if (!isfinite(solventPGrad[dimIdx])) {
1273 throw NumericalProblem("Non-finite potential gradient for solvent 'phase'");
1274 }
1275 }
1276 }
1277
1278 // determine the upstream and downstream DOFs
1279 Evaluation solventPGradNormal = 0.0;
1280 for (unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx) {
1281 solventPGradNormal += solventPGrad[dimIdx]*faceNormal[dimIdx];
1282 }
1283
1284 if (solventPGradNormal > 0) {
1285 solventUpstreamDofIdx_ = j;
1286 solventDownstreamDofIdx_ = i;
1287 }
1288 else {
1289 solventUpstreamDofIdx_ = i;
1290 solventDownstreamDofIdx_ = j;
1291 }
1292
1293 const auto& up = elemCtx.intensiveQuantities(solventUpstreamDofIdx_, timeIdx);
1294
1295 // this is also slightly hacky because it assumes that the derivative of the
1296 // flux between two DOFs only depends on the primary variables in the
1297 // upstream direction. For non-TPFA flux approximation schemes, this is not
1298 // true...
1299 if (solventUpstreamDofIdx_ == i) {
1300 solventVolumeFlux_ = solventPGradNormal * up.solventMobility();
1301 }
1302 else {
1303 solventVolumeFlux_ = solventPGradNormal * scalarValue(up.solventMobility());
1304 }
1305 }
1306
1311 template <class Dummy = bool> // we need to make this method a template to avoid
1312 // compiler errors if it is not instantiated!
1313 void updateVolumeFluxTrans(const ElementContext& elemCtx,
1314 unsigned scvfIdx,
1315 unsigned timeIdx)
1316 {
1317 const ExtensiveQuantities& extQuants = asImp_();
1318
1319 const unsigned interiorDofIdx = extQuants.interiorIndex();
1320 const unsigned exteriorDofIdx = extQuants.exteriorIndex();
1321 assert(interiorDofIdx != exteriorDofIdx);
1322
1323 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
1324 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
1325
1326 const unsigned I = elemCtx.globalSpaceIndex(interiorDofIdx, timeIdx);
1327 const unsigned J = elemCtx.globalSpaceIndex(exteriorDofIdx, timeIdx);
1328
1329 const Scalar thpres = elemCtx.problem().thresholdPressure(I, J);
1330 const Scalar trans = elemCtx.problem().transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
1331 const Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
1332
1333 const Scalar zIn = elemCtx.problem().dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
1334 const Scalar zEx = elemCtx.problem().dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
1335 const Scalar distZ = zIn - zEx;
1336
1337 const Evaluation& rhoIn = intQuantsIn.solventDensity();
1338 const Scalar rhoEx = Toolbox::value(intQuantsEx.solventDensity());
1339 const Evaluation& rhoAvg = rhoIn * 0.5 + rhoEx * 0.5;
1340
1341 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(gasPhaseIdx);
1342 Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(gasPhaseIdx));
1343 pressureExterior += distZ * g * rhoAvg;
1344
1345 Evaluation pressureDiffSolvent = pressureExterior - pressureInterior;
1346 if (std::abs(scalarValue(pressureDiffSolvent)) > thpres) {
1347 if (pressureDiffSolvent < 0.0) {
1348 pressureDiffSolvent += thpres;
1349 }
1350 else {
1351 pressureDiffSolvent -= thpres;
1352 }
1353 }
1354 else {
1355 pressureDiffSolvent = 0.0;
1356 }
1357
1358 if (pressureDiffSolvent > 0.0) {
1359 solventUpstreamDofIdx_ = exteriorDofIdx;
1360 solventDownstreamDofIdx_ = interiorDofIdx;
1361 }
1362 else if (pressureDiffSolvent < 0.0) {
1363 solventUpstreamDofIdx_ = interiorDofIdx;
1364 solventDownstreamDofIdx_ = exteriorDofIdx;
1365 }
1366 else {
1367 // pressure potential gradient is zero; force consistent upstream and
1368 // downstream indices over the intersection regardless of the side which it
1369 // is looked at.
1370 solventUpstreamDofIdx_ = std::min(interiorDofIdx, exteriorDofIdx);
1371 solventDownstreamDofIdx_ = std::max(interiorDofIdx, exteriorDofIdx);
1372 solventVolumeFlux_ = 0.0;
1373 return;
1374 }
1375
1376 const Scalar faceArea = elemCtx.stencil(timeIdx).interiorFace(scvfIdx).area();
1377 const IntensiveQuantities& up = elemCtx.intensiveQuantities(solventUpstreamDofIdx_, timeIdx);
1378 if (solventUpstreamDofIdx_ == interiorDofIdx) {
1379 solventVolumeFlux_ = up.solventMobility() * (-trans / faceArea) * pressureDiffSolvent;
1380 }
1381 else {
1382 solventVolumeFlux_ =
1383 scalarValue(up.solventMobility()) * (-trans / faceArea) * pressureDiffSolvent;
1384 }
1385 }
1386
1387 unsigned solventUpstreamIndex() const
1388 { return solventUpstreamDofIdx_; }
1389
1390 unsigned solventDownstreamIndex() const
1391 { return solventDownstreamDofIdx_; }
1392
1393 const Evaluation& solventVolumeFlux() const
1394 { return solventVolumeFlux_; }
1395
1397 { solventVolumeFlux_ = solventVolumeFlux; }
1398
1399private:
1400 Implementation& asImp_()
1401 { return *static_cast<Implementation*>(this); }
1402
1403 Evaluation solventVolumeFlux_;
1404 unsigned solventUpstreamDofIdx_;
1405 unsigned solventDownstreamDofIdx_;
1406};
1407
1408template <class TypeTag>
1410{
1413
1414public:
1415 void updateVolumeFluxPerm(const ElementContext&,
1416 unsigned,
1417 unsigned)
1418 {}
1419
1420 void updateVolumeFluxTrans(const ElementContext&,
1421 unsigned,
1422 unsigned)
1423 {}
1424
1425 unsigned solventUpstreamIndex() const
1426 { throw std::runtime_error("solventUpstreamIndex() called but solvents are disabled"); }
1427
1428 unsigned solventDownstreamIndex() const
1429 { throw std::runtime_error("solventDownstreamIndex() called but solvents are disabled"); }
1430
1431 const Evaluation& solventVolumeFlux() const
1432 { throw std::runtime_error("solventVolumeFlux() called but solvents are disabled"); }
1433
1434 void setSolventVolumeFlux(const Evaluation& /* solventVolumeFlux */)
1435 { throw std::runtime_error("setSolventVolumeFlux() called but solvents are disabled"); }
1436};
1437
1438} // namespace Opm
1439
1440#endif
Declares the properties required by the black oil model.
Contains the parameters required to extend the black-oil model by solvents.
unsigned solventUpstreamIndex() const
Definition: blackoilsolventmodules.hh:1425
unsigned solventDownstreamIndex() const
Definition: blackoilsolventmodules.hh:1428
void setSolventVolumeFlux(const Evaluation &)
Definition: blackoilsolventmodules.hh:1434
const Evaluation & solventVolumeFlux() const
Definition: blackoilsolventmodules.hh:1431
void updateVolumeFluxPerm(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1415
void updateVolumeFluxTrans(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1420
Provides the solvent specific extensive quantities to the generic black-oil module's extensive quanti...
Definition: blackoilsolventmodules.hh:1180
unsigned solventUpstreamIndex() const
Definition: blackoilsolventmodules.hh:1387
void updateVolumeFluxPerm(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Method which calculates the volume flux of the polymer "phase" using the pressure potential gradient ...
Definition: blackoilsolventmodules.hh:1206
void updateVolumeFluxTrans(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Method which calculates the volume flux of the polymer "phase" using the gas pressure potential diffe...
Definition: blackoilsolventmodules.hh:1313
unsigned solventDownstreamIndex() const
Definition: blackoilsolventmodules.hh:1390
const Evaluation & solventVolumeFlux() const
Definition: blackoilsolventmodules.hh:1393
void setSolventVolumeFlux(const Evaluation &solventVolumeFlux)
Definition: blackoilsolventmodules.hh:1396
const Evaluation & solventDensity() const
Definition: blackoilsolventmodules.hh:1155
const Evaluation & solventInverseFormationVolumeFactor() const
Definition: blackoilsolventmodules.hh:1164
const Evaluation & rsSolw() const
Definition: blackoilsolventmodules.hh:1152
void solventPreSatFuncUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1134
const Evaluation & solventViscosity() const
Definition: blackoilsolventmodules.hh:1158
Scalar solventRefDensity() const
Definition: blackoilsolventmodules.hh:1167
const Evaluation & solventMobility() const
Definition: blackoilsolventmodules.hh:1161
const Evaluation & solventSaturation() const
Definition: blackoilsolventmodules.hh:1149
void solventPvtUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1144
void solventPostSatFuncUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1139
void solventPreSatFuncUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Called before the saturation functions are doing their magic.
Definition: blackoilsolventmodules.hh:580
Implementation & asImp_()
Definition: blackoilsolventmodules.hh:1112
Evaluation rsSolw_
Definition: blackoilsolventmodules.hh:1117
Evaluation solventViscosity_
Definition: blackoilsolventmodules.hh:1119
const Evaluation & solventViscosity() const
Definition: blackoilsolventmodules.hh:866
void solventPvtUpdate_(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Update the intensive PVT properties needed to handle solvents from the primary variables.
Definition: blackoilsolventmodules.hh:793
const Evaluation & rsSolw() const
Definition: blackoilsolventmodules.hh:860
const Evaluation & solventInverseFormationVolumeFactor() const
Definition: blackoilsolventmodules.hh:872
const Evaluation & solventDensity() const
Definition: blackoilsolventmodules.hh:863
void solventPostSatFuncUpdate_(const Problem &problem, const PrimaryVariables &priVars, const unsigned globalSpaceIdx, const unsigned timeIdx, const LinearizationType linearizationType)
Definition: blackoilsolventmodules.hh:641
const Evaluation & solventMobility() const
Definition: blackoilsolventmodules.hh:869
Evaluation hydrocarbonSaturation_
Definition: blackoilsolventmodules.hh:1115
void solventPreSatFuncUpdate_(const PrimaryVariables &priVars, const unsigned timeIdx, const LinearizationType linearizationType)
Definition: blackoilsolventmodules.hh:588
void solventPostSatFuncUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Called after the saturation functions have been doing their magic.
Definition: blackoilsolventmodules.hh:617
Evaluation solventSaturation_
Definition: blackoilsolventmodules.hh:1116
const Evaluation & solventSaturation() const
Definition: blackoilsolventmodules.hh:857
Evaluation solventDensity_
Definition: blackoilsolventmodules.hh:1118
Evaluation solventInvFormationVolumeFactor_
Definition: blackoilsolventmodules.hh:1121
Scalar solventRefDensity_
Definition: blackoilsolventmodules.hh:1123
Scalar solventRefDensity() const
Definition: blackoilsolventmodules.hh:876
Evaluation solventMobility_
Definition: blackoilsolventmodules.hh:1120
Provides the volumetric quantities required for the equations needed by the solvents extension of the...
Definition: blackoilsolventmodules.hh:541
Contains the high level supplements required to extend the black oil model by solvents.
Definition: blackoilsolventmodules.hh:68
static void setSolventPvt(const SolventPvt &value)
Specify the solvent PVT of a all PVT regions.
Definition: blackoilsolventmodules.hh:107
static const TabulatedFunction & misc(const ElemContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:411
static const TabulatedFunction & sof2Krn(const ElemContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:401
static Scalar computeResidualError(const EqVector &resid)
Return how much a residual is considered an error.
Definition: blackoilsolventmodules.hh:332
static Scalar tlMixParamDensity(const ElemContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:492
static const TabulatedFunction & ssfnKrg(const ElemContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:381
static const BrineH2Pvt & brineH2Pvt()
Definition: blackoilsolventmodules.hh:377
static Scalar primaryVarWeight(unsigned pvIdx)
Definition: blackoilsolventmodules.hh:151
static void setParams(BlackOilSolventParams< Scalar > &&params)
Set parameters.
Definition: blackoilsolventmodules.hh:101
static std::string eqName(unsigned eqIdx)
Definition: blackoilsolventmodules.hh:169
static Scalar eqWeight(unsigned eqIdx)
Definition: blackoilsolventmodules.hh:176
static void serializeEntity(const Model &model, std::ostream &outstream, const DofEntity &dof)
Definition: blackoilsolventmodules.hh:339
static bool isSolubleInWater()
Definition: blackoilsolventmodules.hh:523
static bool isCO2Sol()
Definition: blackoilsolventmodules.hh:526
static const TabulatedFunction & sgcwmis(const ElemContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:461
static Scalar tlMixParamViscosity(const ElemContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:481
static void registerOutputModules(Model &model, Simulator &simulator)
Register all solvent specific VTK and ECL output modules.
Definition: blackoilsolventmodules.hh:126
static const TabulatedFunction & tlPMixTable(const ElemContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:471
static const TabulatedFunction & sorwmis(const ElemContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:451
static Value solubilityLimit(unsigned pvtIdx, const Value &temperature, const Value &pressure, const Value &saltConcentration)
Definition: blackoilsolventmodules.hh:505
static const TabulatedFunction & pmisc(const ElemContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:421
static void setIsMiscible(const bool isMiscible)
Definition: blackoilsolventmodules.hh:110
static void registerParameters()
Register all run-time parameters for the black-oil solvent module.
Definition: blackoilsolventmodules.hh:116
static Scalar computeUpdateError(const PrimaryVariables &, const EqVector &)
Return how much a Newton-Raphson update is considered an error.
Definition: blackoilsolventmodules.hh:320
static const SolventPvt & solventPvt()
Definition: blackoilsolventmodules.hh:365
static bool isH2Sol()
Definition: blackoilsolventmodules.hh:529
static bool isMiscible()
Definition: blackoilsolventmodules.hh:501
static const H2GasPvt & h2GasPvt()
Definition: blackoilsolventmodules.hh:371
static void updatePrimaryVars(PrimaryVariables &newPv, const PrimaryVariables &oldPv, const EqVector &delta)
Do a Newton-Raphson update the primary variables of the solvents.
Definition: blackoilsolventmodules.hh:307
static void assignPrimaryVars(PrimaryVariables &priVars, Scalar solventSaturation, Scalar solventRsw)
Assign the solvent specific primary variables to a PrimaryVariables object.
Definition: blackoilsolventmodules.hh:286
static const TabulatedFunction & msfnKrsg(const ElemContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:431
static const TabulatedFunction & ssfnKrs(const ElemContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:391
static std::string primaryVarName(unsigned pvIdx)
Definition: blackoilsolventmodules.hh:144
static bool eqApplies(unsigned eqIdx)
Definition: blackoilsolventmodules.hh:159
static const TabulatedFunction & msfnKro(const ElemContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:441
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilsolventmodules.hh:185
static void deserializeEntity(Model &model, std::istream &instream, const DofEntity &dof)
Definition: blackoilsolventmodules.hh:350
static const Co2GasPvt & co2GasPvt()
Definition: blackoilsolventmodules.hh:368
static const BrineCo2Pvt & brineCo2Pvt()
Definition: blackoilsolventmodules.hh:374
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:218
static bool primaryVarApplies(unsigned pvIdx)
Definition: blackoilsolventmodules.hh:134
Callback class for a phase pressure.
Definition: quantitycallbacks.hh:85
void setPhaseIndex(unsigned phaseIdx)
Set the index of the fluid phase for which the pressure should be returned.
Definition: quantitycallbacks.hh:109
VTK output module for the black oil model's solvent related quantities.
Definition: vtkblackoilsolventmodule.hpp:54
static void registerParameters()
Register all run-time parameters for the multi-phase VTK output module.
Definition: vtkblackoilsolventmodule.hpp:84
Defines the common parameters for the porous medium multi-phase models.
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
This method contains all callback classes for quantities that are required by some extensive quantiti...
Struct holding the parameters for the BlackOilSolventModule class.
Definition: blackoilsolventparams.hpp:49
::Opm::SolventPvt< Scalar > SolventPvt
Definition: blackoilsolventparams.hpp:54
::Opm::H2GasPvt< Scalar > H2GasPvt
Definition: blackoilsolventparams.hpp:53
::Opm::BrineCo2Pvt< Scalar > BrineCo2Pvt
Definition: blackoilsolventparams.hpp:50
::Opm::Co2GasPvt< Scalar > Co2GasPvt
Definition: blackoilsolventparams.hpp:52
Tabulated1DFunction< Scalar > TabulatedFunction
Definition: blackoilsolventparams.hpp:55
::Opm::BrineH2Pvt< Scalar > BrineH2Pvt
Definition: blackoilsolventparams.hpp:51
Definition: linearizationtype.hh:34