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