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
47#include <algorithm>
48#include <array>
49#include <cassert>
50#include <cmath>
51#include <istream>
52#include <memory>
53#include <ostream>
54#include <stdexcept>
55#include <string>
56
57namespace Opm {
58
64template <class TypeTag, bool enableSolventV = getPropValue<TypeTag, Properties::EnableSolvent>()>
66{
79 using Toolbox = MathToolbox<Evaluation>;
80 using SolventPvt = typename BlackOilSolventParams<Scalar>::SolventPvt;
81 using Co2GasPvt = typename BlackOilSolventParams<Scalar>::Co2GasPvt;
82 using H2GasPvt = typename BlackOilSolventParams<Scalar>::H2GasPvt;
83 using BrineCo2Pvt = typename BlackOilSolventParams<Scalar>::BrineCo2Pvt;
84 using BrineH2Pvt = typename BlackOilSolventParams<Scalar>::BrineH2Pvt;
85
86 using TabulatedFunction = typename BlackOilSolventParams<Scalar>::TabulatedFunction;
87
88 static constexpr unsigned solventSaturationIdx = Indices::solventSaturationIdx;
89 static constexpr unsigned contiSolventEqIdx = Indices::contiSolventEqIdx;
90 static constexpr unsigned enableSolvent = enableSolventV;
91 static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
92 static constexpr unsigned numPhases = FluidSystem::numPhases;
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 static const TabulatedFunction& ssfnKrg(const ElementContext& elemCtx,
379 unsigned scvIdx,
380 unsigned timeIdx)
381 {
382 const unsigned satnumRegionIdx =
383 elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
384 return params_.ssfnKrg_[satnumRegionIdx];
385 }
386
387 static const TabulatedFunction& ssfnKrs(const ElementContext& elemCtx,
388 unsigned scvIdx,
389 unsigned timeIdx)
390 {
391 const unsigned satnumRegionIdx =
392 elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
393 return params_.ssfnKrs_[satnumRegionIdx];
394 }
395
396 static const TabulatedFunction& sof2Krn(const ElementContext& elemCtx,
397 unsigned scvIdx,
398 unsigned timeIdx)
399 {
400 const unsigned satnumRegionIdx =
401 elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
402 return params_.sof2Krn_[satnumRegionIdx];
403 }
404
405 static const TabulatedFunction& misc(const ElementContext& elemCtx,
406 unsigned scvIdx,
407 unsigned timeIdx)
408 {
409 const unsigned miscnumRegionIdx =
410 elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
411 return params_.misc_[miscnumRegionIdx];
412 }
413
414 static const TabulatedFunction& pmisc(const ElementContext& elemCtx,
415 unsigned scvIdx,
416 unsigned timeIdx)
417 {
418 const unsigned miscnumRegionIdx =
419 elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
420 return params_.pmisc_[miscnumRegionIdx];
421 }
422
423 static const TabulatedFunction& msfnKrsg(const ElementContext& elemCtx,
424 unsigned scvIdx,
425 unsigned timeIdx)
426 {
427 const unsigned satnumRegionIdx =
428 elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
429 return params_.msfnKrsg_[satnumRegionIdx];
430 }
431
432 static const TabulatedFunction& msfnKro(const ElementContext& elemCtx,
433 unsigned scvIdx,
434 unsigned timeIdx)
435 {
436 const unsigned satnumRegionIdx =
437 elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
438 return params_.msfnKro_[satnumRegionIdx];
439 }
440
441 static const TabulatedFunction& sorwmis(const ElementContext& elemCtx,
442 unsigned scvIdx,
443 unsigned timeIdx)
444 {
445 const unsigned miscnumRegionIdx =
446 elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
447 return params_.sorwmis_[miscnumRegionIdx];
448 }
449
450 static const TabulatedFunction& sgcwmis(const ElementContext& elemCtx,
451 unsigned scvIdx,
452 unsigned timeIdx)
453 {
454 const unsigned miscnumRegionIdx =
455 elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
456 return params_.sgcwmis_[miscnumRegionIdx];
457 }
458
459 static const TabulatedFunction& tlPMixTable(const ElementContext& elemCtx,
460 unsigned scvIdx,
461 unsigned timeIdx)
462 {
463 const unsigned miscnumRegionIdx =
464 elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
465 return params_.tlPMixTable_[miscnumRegionIdx];
466 }
467
468 static Scalar tlMixParamViscosity(const ElementContext& elemCtx,
469 unsigned scvIdx,
470 unsigned timeIdx)
471 {
472 const unsigned miscnumRegionIdx =
473 elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
474 return params_.tlMixParamViscosity_[miscnumRegionIdx];
475 }
476
477 static Scalar tlMixParamDensity(const ElementContext& elemCtx,
478 unsigned scvIdx,
479 unsigned timeIdx)
480 {
481 const unsigned miscnumRegionIdx =
482 elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
483 return params_.tlMixParamDensity_[miscnumRegionIdx];
484 }
485
486 static bool isMiscible()
487 { return params_.isMiscible_; }
488
489 template <class Value>
490 static Value solubilityLimit(unsigned pvtIdx, const Value& temperature,
491 const Value& pressure, const Value& saltConcentration)
492 {
493 if (!isSolubleInWater()) {
494 return 0.0;
495 }
496
497 assert(isCO2Sol() || isH2Sol());
498 if (isCO2Sol()) {
499 return brineCo2Pvt().saturatedGasDissolutionFactor(pvtIdx, temperature,
500 pressure, saltConcentration);
501 }
502 else {
503 return brineH2Pvt().saturatedGasDissolutionFactor(pvtIdx, temperature,
504 pressure, saltConcentration);
505 }
506 }
507
508 static bool isSolubleInWater()
509 { return params_.rsSolw_active_; }
510
511 static bool isCO2Sol()
512 { return params_.co2sol_; }
513
514 static bool isH2Sol()
515 { return params_.h2sol_; }
516
517private:
518 static BlackOilSolventParams<Scalar> params_; // the krg(Fs) column of the SSFN table
519};
520
521template <class TypeTag, bool enableSolventV>
522BlackOilSolventParams<typename BlackOilSolventModule<TypeTag, enableSolventV>::Scalar>
523BlackOilSolventModule<TypeTag, enableSolventV>::params_;
524
525template <class TypeTag, bool enableSolventV>
527
535template <class TypeTag>
536class BlackOilSolventIntensiveQuantities<TypeTag, /*enableSolventV=*/true>
537{
539
547
549
550 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
551 static constexpr int solventSaturationIdx = Indices::solventSaturationIdx;
552 static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
553 static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
554 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
555 static constexpr double cutOff = 1e-12;
556
557public:
564 void solventPreSatFuncUpdate_(const ElementContext& elemCtx,
565 unsigned dofIdx,
566 unsigned timeIdx)
567 {
568 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
569
570 auto& fs = asImp_().fluidState_;
571 solventSaturation_ = 0.0;
572 if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
573 solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx,
574 elemCtx.linearizationType());
575 }
576
577 hydrocarbonSaturation_ = fs.saturation(gasPhaseIdx);
578
579 // apply a cut-off. Don't waste calculations if no solvent
580 if (solventSaturation().value() < cutOff) {
581 return;
582 }
583
584 // make the saturation of the gas phase which is used by the saturation functions
585 // the sum of the solvent "saturation" and the saturation the hydrocarbon gas.
586 fs.setSaturation(gasPhaseIdx, hydrocarbonSaturation_ + solventSaturation_);
587 }
588
596 void solventPostSatFuncUpdate_(const ElementContext& elemCtx,
597 unsigned dofIdx,
598 unsigned timeIdx)
599 {
600 // revert the gas "saturation" of the fluid state back to the saturation of the
601 // hydrocarbon gas.
602 auto& fs = asImp_().fluidState_;
603 fs.setSaturation(gasPhaseIdx, hydrocarbonSaturation_);
604
605 // update rsSolw. This needs to be done after the pressure is defined in the fluid state.
606 rsSolw_ = 0.0;
607 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
608 if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
609 rsSolw_ = SolventModule::solubilityLimit(asImp_().pvtRegionIndex(),
610 fs.temperature(waterPhaseIdx),
611 fs.pressure(waterPhaseIdx),
612 fs.saltConcentration());
613 } else if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Rsolw) {
614 rsSolw_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx,
615 elemCtx.linearizationType());
616 }
617
618 solventMobility_ = 0.0;
619
620 // apply a cut-off. Don't waste calculations if no solvent
621 if (solventSaturation().value() < cutOff) {
622 return;
623 }
624
625 // Pressure effects on capillary pressure miscibility
626 if (SolventModule::isMiscible()) {
627 const Evaluation& p =
628 FluidSystem::phaseIsActive(oilPhaseIdx)
629 ? fs.pressure(oilPhaseIdx)
630 : fs.pressure(gasPhaseIdx);
631 const Evaluation pmisc =
632 SolventModule::pmisc(elemCtx, dofIdx, timeIdx).eval(p, /*extrapolate=*/true);
633 const Evaluation& pgImisc = fs.pressure(gasPhaseIdx);
634
635 // compute capillary pressure for miscible fluid
636 const auto& problem = elemCtx.problem();
637 Evaluation pgMisc = 0.0;
638 std::array<Evaluation, numPhases> pC;
639 const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, timeIdx);
640 MaterialLaw::capillaryPressures(pC, materialParams, fs);
641
642 // oil is the reference phase for pressure
643 const auto linearizationType = elemCtx.linearizationType();
644 if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pg) {
645 pgMisc = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx,
646 linearizationType);
647 }
648 else {
649 const Evaluation& po =
650 priVars.makeEvaluation(Indices::pressureSwitchIdx,
651 timeIdx, linearizationType);
652 pgMisc = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
653 }
654
655 fs.setPressure(gasPhaseIdx, pmisc * pgMisc + (1.0 - pmisc) * pgImisc);
656 }
657
658 const Evaluation gasSolventSat = hydrocarbonSaturation_ + solventSaturation_;
659
660 if (gasSolventSat.value() < cutOff) { // avoid division by zero
661 return;
662 }
663
664 const Evaluation Fhydgas = hydrocarbonSaturation_ / gasSolventSat;
665 const Evaluation Fsolgas = solventSaturation_ / gasSolventSat;
666
667 // account for miscibility of oil and solvent
668 if (SolventModule::isMiscible() && FluidSystem::phaseIsActive(oilPhaseIdx)) {
669 const auto& misc = SolventModule::misc(elemCtx, dofIdx, timeIdx);
670 const auto& pmisc = SolventModule::pmisc(elemCtx, dofIdx, timeIdx);
671 const Evaluation& p =
672 FluidSystem::phaseIsActive(oilPhaseIdx)
673 ? fs.pressure(oilPhaseIdx)
674 : fs.pressure(gasPhaseIdx);
675 const Evaluation miscibility =
676 misc.eval(Fsolgas, /*extrapolate=*/true) *
677 pmisc.eval(p, /*extrapolate=*/true);
678
679 // TODO adjust endpoints of sn and ssg
680 const unsigned cellIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
681 const auto& materialLawManager = elemCtx.problem().materialLawManager();
682 const auto& scaledDrainageInfo =
683 materialLawManager->oilWaterScaledEpsInfoDrainage(cellIdx);
684
685 const Scalar sogcr = scaledDrainageInfo.Sogcr;
686 Evaluation sor = sogcr;
687 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
688 const Evaluation& sw = fs.saturation(waterPhaseIdx);
689 const auto& sorwmis = SolventModule::sorwmis(elemCtx, dofIdx, timeIdx);
690 sor = miscibility *
691 sorwmis.eval(sw, /*extrapolate=*/true) + (1.0 - miscibility) * sogcr;
692 }
693 const Scalar sgcr = scaledDrainageInfo.Sgcr;
694 Evaluation sgc = sgcr;
695 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
696 const Evaluation& sw = fs.saturation(waterPhaseIdx);
697 const auto& sgcwmis = SolventModule::sgcwmis(elemCtx, dofIdx, timeIdx);
698 sgc = miscibility * sgcwmis.eval(sw, /*extrapolate=*/true) + (1.0 - miscibility) * sgcr;
699 }
700
701 Evaluation oilGasSolventSat = gasSolventSat;
702 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
703 oilGasSolventSat += fs.saturation(oilPhaseIdx);
704 }
705 const Evaluation zero = 0.0;
706 const Evaluation oilGasSolventEffSat = std::max(oilGasSolventSat - sor - sgc, zero);
707
708 Evaluation F_totalGas = 0.0;
709 if (oilGasSolventEffSat.value() > cutOff) {
710 const Evaluation gasSolventEffSat = std::max(gasSolventSat - sgc, zero);
711 F_totalGas = gasSolventEffSat / oilGasSolventEffSat;
712 }
713 const auto& msfnKro = SolventModule::msfnKro(elemCtx, dofIdx, timeIdx);
714 const auto& msfnKrsg = SolventModule::msfnKrsg(elemCtx, dofIdx, timeIdx);
715 const auto& sof2Krn = SolventModule::sof2Krn(elemCtx, dofIdx, timeIdx);
716
717 const Evaluation mkrgt = msfnKrsg.eval(F_totalGas, /*extrapolate=*/true) *
718 sof2Krn.eval(oilGasSolventSat, /*extrapolate=*/true);
719 const Evaluation mkro = msfnKro.eval(F_totalGas, /*extrapolate=*/true) *
720 sof2Krn.eval(oilGasSolventSat, /*extrapolate=*/true);
721
722 Evaluation& kro = asImp_().mobility_[oilPhaseIdx];
723 Evaluation& krg = asImp_().mobility_[gasPhaseIdx];
724
725 // combine immiscible and miscible part of the relperm
726 krg *= 1.0 - miscibility;
727 krg += miscibility * mkrgt;
728 kro *= 1.0 - miscibility;
729 kro += miscibility * mkro;
730 }
731
732 // compute the mobility of the solvent "phase" and modify the gas phase
733 const auto& ssfnKrg = SolventModule::ssfnKrg(elemCtx, dofIdx, timeIdx);
734 const auto& ssfnKrs = SolventModule::ssfnKrs(elemCtx, dofIdx, timeIdx);
735
736 Evaluation& krg = asImp_().mobility_[gasPhaseIdx];
737 solventMobility_ = krg * ssfnKrs.eval(Fsolgas, /*extrapolate=*/true);
738 krg *= ssfnKrg.eval(Fhydgas, /*extrapolate=*/true);
739 }
740
747 void solventPvtUpdate_(const ElementContext& elemCtx,
748 unsigned scvIdx,
749 unsigned timeIdx)
750 {
751 const auto& iq = asImp_();
752 const unsigned pvtRegionIdx = iq.pvtRegionIndex();
753 const Evaluation& T = iq.fluidState().temperature(gasPhaseIdx);
754 const Evaluation& p = iq.fluidState().pressure(gasPhaseIdx);
755
756 const Evaluation rv = 0.0;
757 const Evaluation rvw = 0.0;
758 if (SolventModule::isCO2Sol() || SolventModule::isH2Sol() ){
759 if (SolventModule::isCO2Sol()) {
760 const auto& co2gasPvt = SolventModule::co2GasPvt();
761 solventInvFormationVolumeFactor_ =
762 co2gasPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rv, rvw);
763 solventRefDensity_ = co2gasPvt.gasReferenceDensity(pvtRegionIdx);
764 solventViscosity_ = co2gasPvt.viscosity(pvtRegionIdx, T, p, rv, rvw);
765
766 const auto& brineCo2Pvt = SolventModule::brineCo2Pvt();
767 auto& fs = asImp_().fluidState_;
768 const auto& bw = brineCo2Pvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rsSolw());
769
770 const auto denw = bw * brineCo2Pvt.waterReferenceDensity(pvtRegionIdx) +
771 rsSolw() * bw * brineCo2Pvt.gasReferenceDensity(pvtRegionIdx);
772 fs.setDensity(waterPhaseIdx, denw);
773 fs.setInvB(waterPhaseIdx, bw);
774 Evaluation& mobw = asImp_().mobility_[waterPhaseIdx];
775 const auto& muWat = fs.viscosity(waterPhaseIdx);
776 const auto& muWatEff = brineCo2Pvt.viscosity(pvtRegionIdx, T, p, rsSolw());
777 mobw *= muWat / muWatEff;
778 } else {
779 const auto& h2gasPvt = SolventModule::h2GasPvt();
780 solventInvFormationVolumeFactor_ =
781 h2gasPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rv, rvw);
782 solventRefDensity_ = h2gasPvt.gasReferenceDensity(pvtRegionIdx);
783 solventViscosity_ = h2gasPvt.viscosity(pvtRegionIdx, T, p, rv, rvw);
784
785 const auto& brineH2Pvt = SolventModule::brineH2Pvt();
786 auto& fs = asImp_().fluidState_;
787 const auto& bw = brineH2Pvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rsSolw());
788
789 const auto denw = bw * brineH2Pvt.waterReferenceDensity(pvtRegionIdx) +
790 rsSolw() * bw * brineH2Pvt.gasReferenceDensity(pvtRegionIdx);
791 fs.setDensity(waterPhaseIdx, denw);
792 fs.setInvB(waterPhaseIdx, bw);
793 Evaluation& mobw = asImp_().mobility_[waterPhaseIdx];
794 const auto& muWat = fs.viscosity(waterPhaseIdx);
795 const auto& muWatEff = brineH2Pvt.viscosity(pvtRegionIdx, T, p, rsSolw());
796 mobw *= muWat / muWatEff;
797 }
798 } else {
799 const auto& solventPvt = SolventModule::solventPvt();
800 solventInvFormationVolumeFactor_ = solventPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
801 solventRefDensity_ = solventPvt.referenceDensity(pvtRegionIdx);
802 solventViscosity_ = solventPvt.viscosity(pvtRegionIdx, T, p);
803 }
804
805 solventDensity_ = solventInvFormationVolumeFactor_*solventRefDensity_;
806 effectiveProperties(elemCtx, scvIdx, timeIdx);
807
808 solventMobility_ /= solventViscosity_;
809 }
810
811 const Evaluation& solventSaturation() const
812 { return solventSaturation_; }
813
814 const Evaluation& rsSolw() const
815 { return rsSolw_; }
816
817 const Evaluation& solventDensity() const
818 { return solventDensity_; }
819
820 const Evaluation& solventViscosity() const
821 { return solventViscosity_; }
822
823 const Evaluation& solventMobility() const
824 { return solventMobility_; }
825
826 const Evaluation& solventInverseFormationVolumeFactor() const
827 { return solventInvFormationVolumeFactor_; }
828
829 // This could be stored pr pvtRegion instead
830 Scalar solventRefDensity() const
831 { return solventRefDensity_; }
832
833private:
834 // Computes the effective properties based on
835 // Todd-Longstaff mixing model.
836 void effectiveProperties(const ElementContext& elemCtx,
837 unsigned scvIdx,
838 unsigned timeIdx)
839 {
840 if (!SolventModule::isMiscible()) {
841 return;
842 }
843
844 // Don't waste calculations if no solvent
845 // Apply a cut-off for small and negative solvent saturations
846 if (solventSaturation() < cutOff) {
847 return;
848 }
849
850 // We plan to update the fluidstate with the effective
851 // properties
852 auto& fs = asImp_().fluidState_;
853
854 // Compute effective saturations
855 Evaluation oilEffSat = 0.0;
856 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
857 oilEffSat = fs.saturation(oilPhaseIdx);
858 }
859 Evaluation gasEffSat = fs.saturation(gasPhaseIdx);
860 Evaluation solventEffSat = solventSaturation();
861 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
862 const auto& sorwmis = SolventModule::sorwmis(elemCtx, scvIdx, timeIdx);
863 const auto& sgcwmis = SolventModule::sgcwmis(elemCtx, scvIdx, timeIdx);
864 const Evaluation zero = 0.0;
865 const Evaluation& sw = fs.saturation(waterPhaseIdx);
866 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
867 oilEffSat = std::max(oilEffSat - sorwmis.eval(sw, /*extrapolate=*/true), zero);
868 }
869 gasEffSat = std::max(gasEffSat - sgcwmis.eval(sw, /*extrapolate=*/true), zero);
870 solventEffSat = std::max(solventEffSat - sgcwmis.eval(sw, /*extrapolate=*/true), zero);
871 }
872 const Evaluation oilGasSolventEffSat = oilEffSat + gasEffSat + solventEffSat;
873 const Evaluation oilSolventEffSat = oilEffSat + solventEffSat;
874 const Evaluation solventGasEffSat = solventEffSat + gasEffSat;
875
876 // Compute effective viscosities
877
878 // Mixing parameter for viscosity
879 // The pressureMixingParameter represent the miscibility of the solvent while the mixingParameterViscosity the effect of the porous media.
880 // The pressureMixingParameter is not implemented in ecl100.
881 const Evaluation& p =
882 FluidSystem::phaseIsActive(oilPhaseIdx)
883 ? fs.pressure(oilPhaseIdx)
884 : fs.pressure(gasPhaseIdx);
885 // account for pressure effects
886 const auto& pmiscTable = SolventModule::pmisc(elemCtx, scvIdx, timeIdx);
887 const Evaluation pmisc = pmiscTable.eval(p, /*extrapolate=*/true);
888 const auto& tlPMixTable = SolventModule::tlPMixTable(elemCtx, scvIdx, timeIdx);
889 const Evaluation tlMixParamMu = SolventModule::tlMixParamViscosity(elemCtx, scvIdx, timeIdx) *
890 tlPMixTable.eval(p, /*extrapolate=*/true);
891
892 const Evaluation& muGas = fs.viscosity(gasPhaseIdx);
893 const Evaluation& muSolvent = solventViscosity_;
894
895 assert(muGas.value() > 0);
896 assert(muSolvent.value() > 0);
897 const Evaluation muGasPow = pow(muGas, 0.25);
898 const Evaluation muSolventPow = pow(muSolvent, 0.25);
899
900 Evaluation muMixSolventGas = muGas;
901 if (solventGasEffSat > cutOff) {
902 muMixSolventGas *= muSolvent / pow(((gasEffSat / solventGasEffSat) * muSolventPow) +
903 ((solventEffSat / solventGasEffSat) * muGasPow), 4.0);
904 }
905
906 Evaluation muOil = 1.0;
907 Evaluation muOilPow = 1.0;
908 Evaluation muMixOilSolvent = 1.0;
909 Evaluation muOilEff = 1.0;
910 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
911 muOil = fs.viscosity(oilPhaseIdx);
912 assert(muOil.value() > 0);
913 muOilPow = pow(muOil, 0.25);
914 muMixOilSolvent = muOil;
915 if (oilSolventEffSat > cutOff) {
916 muMixOilSolvent *= muSolvent / pow(((oilEffSat / oilSolventEffSat) * muSolventPow) +
917 ((solventEffSat / oilSolventEffSat) * muOilPow), 4.0);
918 }
919
920 muOilEff = pow(muOil,1.0 - tlMixParamMu) * pow(muMixOilSolvent, tlMixParamMu);
921 }
922 Evaluation muMixSolventGasOil = muOil;
923 if (oilGasSolventEffSat > cutOff) {
924 muMixSolventGasOil *= muSolvent * muGas / pow(((oilEffSat / oilGasSolventEffSat) *
925 muSolventPow * muGasPow) +
926 ((solventEffSat / oilGasSolventEffSat) *
927 muOilPow * muGasPow) +
928 ((gasEffSat / oilGasSolventEffSat) *
929 muSolventPow * muOilPow), 4.0);
930 }
931
932 Evaluation muGasEff = pow(muGas,1.0 - tlMixParamMu) * pow(muMixSolventGas, tlMixParamMu);
933 const Evaluation muSolventEff = pow(muSolvent,1.0 - tlMixParamMu) * pow(muMixSolventGasOil, tlMixParamMu);
934
935 // Compute effective densities
936 const Evaluation& rhoGas = fs.density(gasPhaseIdx);
937 Evaluation rhoOil = 0.0;
938 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
939 rhoOil = fs.density(oilPhaseIdx);
940 }
941
942 const Evaluation& rhoSolvent = solventDensity_;
943
944 // Mixing parameter for density
945 // The pressureMixingParameter represent the miscibility of the
946 // solvent while the mixingParameterDenisty the effect of the porous media.
947 // The pressureMixingParameter is not implemented in ecl100.
948 const Evaluation tlMixParamRho = SolventModule::tlMixParamDensity(elemCtx, scvIdx, timeIdx) *
949 tlPMixTable.eval(p, /*extrapolate=*/true);
950
951 // compute effective viscosities for density calculations. These have to
952 // be recomputed as a different mixing parameter may be used.
953 const Evaluation muOilEffPow = pow(pow(muOil, 1.0 - tlMixParamRho) *
954 pow(muMixOilSolvent, tlMixParamRho), 0.25);
955 const Evaluation muGasEffPow = pow(pow(muGas, 1.0 - tlMixParamRho) *
956 pow(muMixSolventGas, tlMixParamRho), 0.25);
957 const Evaluation muSolventEffPow = pow(pow(muSolvent, 1.0 - tlMixParamRho) *
958 pow(muMixSolventGasOil, tlMixParamRho), 0.25);
959
960 const Evaluation oilGasEffSaturation = oilEffSat + gasEffSat;
961 Evaluation sof = 0.0;
962 Evaluation sgf = 0.0;
963 if (oilGasEffSaturation.value() > cutOff) {
964 sof = oilEffSat / oilGasEffSaturation;
965 sgf = gasEffSat / oilGasEffSaturation;
966 }
967
968 const Evaluation muSolventOilGasPow = muSolventPow * ((sgf * muOilPow) + (sof * muGasPow));
969
970 Evaluation rhoMixSolventGasOil = 0.0;
971 if (oilGasSolventEffSat.value() > cutOff) {
972 rhoMixSolventGasOil = (rhoOil * oilEffSat / oilGasSolventEffSat) +
973 (rhoGas * gasEffSat / oilGasSolventEffSat) +
974 (rhoSolvent * solventEffSat / oilGasSolventEffSat);
975 }
976
977 Evaluation rhoGasEff = 0.0;
978 if (std::abs(muSolventPow.value() - muGasPow.value()) < cutOff) {
979 rhoGasEff = ((1.0 - tlMixParamRho) * rhoGas) +
980 (tlMixParamRho * rhoMixSolventGasOil);
981 }
982 else {
983 const Evaluation solventGasEffFraction = (muGasPow * (muSolventPow - muGasEffPow)) /
984 (muGasEffPow * (muSolventPow - muGasPow));
985 rhoGasEff = (rhoGas * solventGasEffFraction) + (rhoSolvent * (1.0 - solventGasEffFraction));
986 }
987
988 Evaluation rhoSolventEff = 0.0;
989 if (std::abs((muSolventOilGasPow.value() - (muOilPow.value() * muGasPow.value()))) < cutOff) {
990 rhoSolventEff = ((1.0 - tlMixParamRho) * rhoSolvent) + (tlMixParamRho * rhoMixSolventGasOil);
991 }
992 else {
993 const Evaluation sfraction_se = (muSolventOilGasPow -
994 muOilPow * muGasPow * muSolventPow / muSolventEffPow) /
995 (muSolventOilGasPow - (muOilPow * muGasPow));
996 rhoSolventEff = (rhoSolvent * sfraction_se) +
997 (rhoGas * sgf * (1.0 - sfraction_se)) +
998 (rhoOil * sof * (1.0 - sfraction_se));
999 }
1000
1001 // compute invB from densities.
1002 const unsigned pvtRegionIdx = asImp_().pvtRegionIndex();
1003 Evaluation bGasEff = rhoGasEff;
1004 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1005 bGasEff /= (FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) +
1006 FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) * fs.Rv());
1007 } else {
1008 bGasEff /= (FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx));
1009 }
1010 const Evaluation bSolventEff = rhoSolventEff / solventRefDensity();
1011
1012 // copy the unmodified invB factors
1013 const Evaluation bg = fs.invB(gasPhaseIdx);
1014 const Evaluation bs = solventInverseFormationVolumeFactor();
1015 const Evaluation bg_eff = pmisc * bGasEff + (1.0 - pmisc) * bg;
1016 const Evaluation bs_eff = pmisc * bSolventEff + (1.0 - pmisc) * bs;
1017
1018 // set the densities
1019 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1020 fs.setDensity(gasPhaseIdx,
1021 bg_eff *
1022 (FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) +
1023 FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx)*fs.Rv()));
1024 } else {
1025 fs.setDensity(gasPhaseIdx,
1026 bg_eff * FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx));
1027 }
1028 solventDensity_ = bs_eff * solventRefDensity();
1029
1030 // set the mobility
1031 Evaluation& mobg = asImp_().mobility_[gasPhaseIdx];
1032 muGasEff = bg_eff / (pmisc * bGasEff / muGasEff + (1.0 - pmisc) * bg / muGas);
1033 mobg *= muGas / muGasEff;
1034
1035 // Update viscosity of solvent
1036 solventViscosity_ = bs_eff / (pmisc * bSolventEff / muSolventEff +
1037 (1.0 - pmisc) * bs / muSolvent);
1038
1039 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1040 Evaluation rhoOilEff = 0.0;
1041 if (std::abs(muOilPow.value() - muSolventPow.value()) < cutOff) {
1042 rhoOilEff = ((1.0 - tlMixParamRho) * rhoOil) + (tlMixParamRho * rhoMixSolventGasOil);
1043 }
1044 else {
1045 const Evaluation solventOilEffFraction = (muOilPow * (muOilEffPow - muSolventPow)) /
1046 (muOilEffPow * (muOilPow - muSolventPow));
1047 rhoOilEff = (rhoOil * solventOilEffFraction) + (rhoSolvent * (1.0 - solventOilEffFraction));
1048 }
1049 const Evaluation bOilEff =
1050 rhoOilEff / (FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) +
1051 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) * fs.Rs());
1052 const Evaluation bo = fs.invB(oilPhaseIdx);
1053 const Evaluation bo_eff = pmisc * bOilEff + (1.0 - pmisc) * bo;
1054 fs.setDensity(oilPhaseIdx,
1055 bo_eff * (FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) +
1056 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) * fs.Rs()));
1057
1058 // keep the mu*b interpolation
1059 Evaluation& mobo = asImp_().mobility_[oilPhaseIdx];
1060 muOilEff = bo_eff / (pmisc * bOilEff / muOilEff + (1.0 - pmisc) * bo / muOil);
1061 mobo *= muOil / muOilEff;
1062 }
1063 }
1064
1065protected:
1066 Implementation& asImp_()
1067 { return *static_cast<Implementation*>(this); }
1068
1071 Evaluation rsSolw_;
1076
1078};
1079
1080template <class TypeTag>
1082{
1086
1087public:
1088 void solventPreSatFuncUpdate_(const ElementContext&,
1089 unsigned,
1090 unsigned)
1091 {}
1092
1093 void solventPostSatFuncUpdate_(const ElementContext&,
1094 unsigned,
1095 unsigned)
1096 {}
1097
1098 void solventPvtUpdate_(const ElementContext&,
1099 unsigned,
1100 unsigned)
1101 {}
1102
1103 const Evaluation& solventSaturation() const
1104 { throw std::runtime_error("solventSaturation() called but solvents are disabled"); }
1105
1106 const Evaluation& rsSolw() const
1107 { throw std::runtime_error("rsSolw() called but solvents are disabled"); }
1108
1109 const Evaluation& solventDensity() const
1110 { throw std::runtime_error("solventDensity() called but solvents are disabled"); }
1111
1112 const Evaluation& solventViscosity() const
1113 { throw std::runtime_error("solventViscosity() called but solvents are disabled"); }
1114
1115 const Evaluation& solventMobility() const
1116 { throw std::runtime_error("solventMobility() called but solvents are disabled"); }
1117
1118 const Evaluation& solventInverseFormationVolumeFactor() const
1119 { throw std::runtime_error("solventInverseFormationVolumeFactor() called but solvents are disabled"); }
1120
1121 Scalar solventRefDensity() const
1122 { throw std::runtime_error("solventRefDensity() called but solvents are disabled"); }
1123};
1124
1132template <class TypeTag, bool enableSolventV = getPropValue<TypeTag, Properties::EnableSolvent>()>
1134{
1136
1144
1145 using Toolbox = MathToolbox<Evaluation>;
1146
1147 static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
1148 static constexpr int dimWorld = GridView::dimensionworld;
1149
1150 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
1151 using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
1152
1153public:
1158 template <class Dummy = bool> // we need to make this method a template to avoid
1159 // compiler errors if it is not instantiated!
1160 void updateVolumeFluxPerm(const ElementContext& elemCtx,
1161 unsigned scvfIdx,
1162 unsigned timeIdx)
1163 {
1164 const auto& gradCalc = elemCtx.gradientCalculator();
1165 PressureCallback<TypeTag> pressureCallback(elemCtx);
1166
1167 const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
1168 const auto& faceNormal = scvf.normal();
1169
1170 const unsigned i = scvf.interiorIndex();
1171 const unsigned j = scvf.exteriorIndex();
1172
1173 // calculate the "raw" pressure gradient
1174 DimEvalVector solventPGrad;
1175 pressureCallback.setPhaseIndex(gasPhaseIdx);
1176 gradCalc.calculateGradient(solventPGrad,
1177 elemCtx,
1178 scvfIdx,
1179 pressureCallback);
1180 Valgrind::CheckDefined(solventPGrad);
1181
1182 // correct the pressure gradients by the gravitational acceleration
1183 if (Parameters::Get<Parameters::EnableGravity>()) {
1184 // estimate the gravitational acceleration at a given SCV face
1185 // using the arithmetic mean
1186 const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
1187 const auto& gEx = elemCtx.problem().gravity(elemCtx, j, timeIdx);
1188
1189 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
1190 const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
1191
1192 const auto& posIn = elemCtx.pos(i, timeIdx);
1193 const auto& posEx = elemCtx.pos(j, timeIdx);
1194 const auto& posFace = scvf.integrationPos();
1195
1196 // the distance between the centers of the control volumes
1197 DimVector distVecIn(posIn);
1198 DimVector distVecEx(posEx);
1199 DimVector distVecTotal(posEx);
1200
1201 distVecIn -= posFace;
1202 distVecEx -= posFace;
1203 distVecTotal -= posIn;
1204 const Scalar absDistTotalSquared = distVecTotal.two_norm2();
1205
1206 // calculate the hydrostatic pressure at the integration point of the face
1207 auto rhoIn = intQuantsIn.solventDensity();
1208 auto pStatIn = - rhoIn * (gIn * distVecIn);
1209
1210 // the quantities on the exterior side of the face do not influence the
1211 // result for the TPFA scheme, so they can be treated as scalar values.
1212 Scalar rhoEx = Toolbox::value(intQuantsEx.solventDensity());
1213 Scalar pStatEx = - rhoEx * (gEx * distVecEx);
1214
1215 // compute the hydrostatic gradient between the two control volumes (this
1216 // gradient exhibitis the same direction as the vector between the two
1217 // control volume centers and the length (pStaticExterior -
1218 // pStaticInterior)/distanceInteriorToExterior
1219 DimEvalVector f(distVecTotal);
1220 f *= (pStatEx - pStatIn) / absDistTotalSquared;
1221
1222 // calculate the final potential gradient
1223 for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
1224 solventPGrad[dimIdx] += f[dimIdx];
1225
1226 if (!isfinite(solventPGrad[dimIdx])) {
1227 throw NumericalProblem("Non-finite potential gradient for solvent 'phase'");
1228 }
1229 }
1230 }
1231
1232 // determine the upstream and downstream DOFs
1233 Evaluation solventPGradNormal = 0.0;
1234 for (unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx) {
1235 solventPGradNormal += solventPGrad[dimIdx]*faceNormal[dimIdx];
1236 }
1237
1238 if (solventPGradNormal > 0) {
1239 solventUpstreamDofIdx_ = j;
1240 solventDownstreamDofIdx_ = i;
1241 }
1242 else {
1243 solventUpstreamDofIdx_ = i;
1244 solventDownstreamDofIdx_ = j;
1245 }
1246
1247 const auto& up = elemCtx.intensiveQuantities(solventUpstreamDofIdx_, timeIdx);
1248
1249 // this is also slightly hacky because it assumes that the derivative of the
1250 // flux between two DOFs only depends on the primary variables in the
1251 // upstream direction. For non-TPFA flux approximation schemes, this is not
1252 // true...
1253 if (solventUpstreamDofIdx_ == i) {
1254 solventVolumeFlux_ = solventPGradNormal * up.solventMobility();
1255 }
1256 else {
1257 solventVolumeFlux_ = solventPGradNormal * scalarValue(up.solventMobility());
1258 }
1259 }
1260
1265 template <class Dummy = bool> // we need to make this method a template to avoid
1266 // compiler errors if it is not instantiated!
1267 void updateVolumeFluxTrans(const ElementContext& elemCtx,
1268 unsigned scvfIdx,
1269 unsigned timeIdx)
1270 {
1271 const ExtensiveQuantities& extQuants = asImp_();
1272
1273 const unsigned interiorDofIdx = extQuants.interiorIndex();
1274 const unsigned exteriorDofIdx = extQuants.exteriorIndex();
1275 assert(interiorDofIdx != exteriorDofIdx);
1276
1277 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
1278 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
1279
1280 const unsigned I = elemCtx.globalSpaceIndex(interiorDofIdx, timeIdx);
1281 const unsigned J = elemCtx.globalSpaceIndex(exteriorDofIdx, timeIdx);
1282
1283 const Scalar thpres = elemCtx.problem().thresholdPressure(I, J);
1284 const Scalar trans = elemCtx.problem().transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
1285 const Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
1286
1287 const Scalar zIn = elemCtx.problem().dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
1288 const Scalar zEx = elemCtx.problem().dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
1289 const Scalar distZ = zIn - zEx;
1290
1291 const Evaluation& rhoIn = intQuantsIn.solventDensity();
1292 const Scalar rhoEx = Toolbox::value(intQuantsEx.solventDensity());
1293 const Evaluation& rhoAvg = rhoIn * 0.5 + rhoEx * 0.5;
1294
1295 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(gasPhaseIdx);
1296 Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(gasPhaseIdx));
1297 pressureExterior += distZ * g * rhoAvg;
1298
1299 Evaluation pressureDiffSolvent = pressureExterior - pressureInterior;
1300 if (std::abs(scalarValue(pressureDiffSolvent)) > thpres) {
1301 if (pressureDiffSolvent < 0.0) {
1302 pressureDiffSolvent += thpres;
1303 }
1304 else {
1305 pressureDiffSolvent -= thpres;
1306 }
1307 }
1308 else {
1309 pressureDiffSolvent = 0.0;
1310 }
1311
1312 if (pressureDiffSolvent > 0.0) {
1313 solventUpstreamDofIdx_ = exteriorDofIdx;
1314 solventDownstreamDofIdx_ = interiorDofIdx;
1315 }
1316 else if (pressureDiffSolvent < 0.0) {
1317 solventUpstreamDofIdx_ = interiorDofIdx;
1318 solventDownstreamDofIdx_ = exteriorDofIdx;
1319 }
1320 else {
1321 // pressure potential gradient is zero; force consistent upstream and
1322 // downstream indices over the intersection regardless of the side which it
1323 // is looked at.
1324 solventUpstreamDofIdx_ = std::min(interiorDofIdx, exteriorDofIdx);
1325 solventDownstreamDofIdx_ = std::max(interiorDofIdx, exteriorDofIdx);
1326 solventVolumeFlux_ = 0.0;
1327 return;
1328 }
1329
1330 const Scalar faceArea = elemCtx.stencil(timeIdx).interiorFace(scvfIdx).area();
1331 const IntensiveQuantities& up = elemCtx.intensiveQuantities(solventUpstreamDofIdx_, timeIdx);
1332 if (solventUpstreamDofIdx_ == interiorDofIdx) {
1333 solventVolumeFlux_ = up.solventMobility() * (-trans / faceArea) * pressureDiffSolvent;
1334 }
1335 else {
1336 solventVolumeFlux_ =
1337 scalarValue(up.solventMobility()) * (-trans / faceArea) * pressureDiffSolvent;
1338 }
1339 }
1340
1341 unsigned solventUpstreamIndex() const
1342 { return solventUpstreamDofIdx_; }
1343
1344 unsigned solventDownstreamIndex() const
1345 { return solventDownstreamDofIdx_; }
1346
1347 const Evaluation& solventVolumeFlux() const
1348 { return solventVolumeFlux_; }
1349
1351 { solventVolumeFlux_ = solventVolumeFlux; }
1352
1353private:
1354 Implementation& asImp_()
1355 { return *static_cast<Implementation*>(this); }
1356
1357 Evaluation solventVolumeFlux_;
1358 unsigned solventUpstreamDofIdx_;
1359 unsigned solventDownstreamDofIdx_;
1360};
1361
1362template <class TypeTag>
1364{
1367
1368public:
1369 void updateVolumeFluxPerm(const ElementContext&,
1370 unsigned,
1371 unsigned)
1372 {}
1373
1374 void updateVolumeFluxTrans(const ElementContext&,
1375 unsigned,
1376 unsigned)
1377 {}
1378
1379 unsigned solventUpstreamIndex() const
1380 { throw std::runtime_error("solventUpstreamIndex() called but solvents are disabled"); }
1381
1382 unsigned solventDownstreamIndex() const
1383 { throw std::runtime_error("solventDownstreamIndex() called but solvents are disabled"); }
1384
1385 const Evaluation& solventVolumeFlux() const
1386 { throw std::runtime_error("solventVolumeFlux() called but solvents are disabled"); }
1387
1388 void setSolventVolumeFlux(const Evaluation& /* solventVolumeFlux */)
1389 { throw std::runtime_error("setSolventVolumeFlux() called but solvents are disabled"); }
1390};
1391
1392} // namespace Opm
1393
1394#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:1379
unsigned solventDownstreamIndex() const
Definition: blackoilsolventmodules.hh:1382
void setSolventVolumeFlux(const Evaluation &)
Definition: blackoilsolventmodules.hh:1388
const Evaluation & solventVolumeFlux() const
Definition: blackoilsolventmodules.hh:1385
void updateVolumeFluxPerm(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1369
void updateVolumeFluxTrans(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1374
Provides the solvent specific extensive quantities to the generic black-oil module's extensive quanti...
Definition: blackoilsolventmodules.hh:1134
unsigned solventUpstreamIndex() const
Definition: blackoilsolventmodules.hh:1341
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:1160
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:1267
unsigned solventDownstreamIndex() const
Definition: blackoilsolventmodules.hh:1344
const Evaluation & solventVolumeFlux() const
Definition: blackoilsolventmodules.hh:1347
void setSolventVolumeFlux(const Evaluation &solventVolumeFlux)
Definition: blackoilsolventmodules.hh:1350
const Evaluation & solventDensity() const
Definition: blackoilsolventmodules.hh:1109
const Evaluation & solventInverseFormationVolumeFactor() const
Definition: blackoilsolventmodules.hh:1118
const Evaluation & rsSolw() const
Definition: blackoilsolventmodules.hh:1106
void solventPreSatFuncUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1088
const Evaluation & solventViscosity() const
Definition: blackoilsolventmodules.hh:1112
Scalar solventRefDensity() const
Definition: blackoilsolventmodules.hh:1121
const Evaluation & solventMobility() const
Definition: blackoilsolventmodules.hh:1115
const Evaluation & solventSaturation() const
Definition: blackoilsolventmodules.hh:1103
void solventPvtUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1098
void solventPostSatFuncUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1093
void solventPreSatFuncUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Called before the saturation functions are doing their magic.
Definition: blackoilsolventmodules.hh:564
Implementation & asImp_()
Definition: blackoilsolventmodules.hh:1066
Evaluation rsSolw_
Definition: blackoilsolventmodules.hh:1071
Evaluation solventViscosity_
Definition: blackoilsolventmodules.hh:1073
const Evaluation & solventViscosity() const
Definition: blackoilsolventmodules.hh:820
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:747
const Evaluation & rsSolw() const
Definition: blackoilsolventmodules.hh:814
const Evaluation & solventInverseFormationVolumeFactor() const
Definition: blackoilsolventmodules.hh:826
const Evaluation & solventDensity() const
Definition: blackoilsolventmodules.hh:817
const Evaluation & solventMobility() const
Definition: blackoilsolventmodules.hh:823
Evaluation hydrocarbonSaturation_
Definition: blackoilsolventmodules.hh:1069
void solventPostSatFuncUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Called after the saturation functions have been doing their magic.
Definition: blackoilsolventmodules.hh:596
Evaluation solventSaturation_
Definition: blackoilsolventmodules.hh:1070
const Evaluation & solventSaturation() const
Definition: blackoilsolventmodules.hh:811
Evaluation solventDensity_
Definition: blackoilsolventmodules.hh:1072
Evaluation solventInvFormationVolumeFactor_
Definition: blackoilsolventmodules.hh:1075
Scalar solventRefDensity_
Definition: blackoilsolventmodules.hh:1077
Scalar solventRefDensity() const
Definition: blackoilsolventmodules.hh:830
Evaluation solventMobility_
Definition: blackoilsolventmodules.hh:1074
Provides the volumetric quantities required for the equations needed by the solvents extension of the...
Definition: blackoilsolventmodules.hh:526
Contains the high level supplements required to extend the black oil model by solvents.
Definition: blackoilsolventmodules.hh:66
static void setSolventPvt(const SolventPvt &value)
Specify the solvent PVT of a all PVT regions.
Definition: blackoilsolventmodules.hh:105
static const TabulatedFunction & sof2Krn(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:396
static const TabulatedFunction & ssfnKrs(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:387
static const TabulatedFunction & pmisc(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:414
static const TabulatedFunction & msfnKrsg(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:423
static const TabulatedFunction & sgcwmis(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:450
static Scalar computeResidualError(const EqVector &resid)
Return how much a residual is considered an error.
Definition: blackoilsolventmodules.hh:330
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:508
static bool isCO2Sol()
Definition: blackoilsolventmodules.hh:511
static const TabulatedFunction & msfnKro(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:432
static void registerOutputModules(Model &model, Simulator &simulator)
Register all solvent specific VTK and ECL output modules.
Definition: blackoilsolventmodules.hh:124
static Scalar tlMixParamDensity(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:477
static Value solubilityLimit(unsigned pvtIdx, const Value &temperature, const Value &pressure, const Value &saltConcentration)
Definition: blackoilsolventmodules.hh:490
static const TabulatedFunction & ssfnKrg(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:378
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:514
static bool isMiscible()
Definition: blackoilsolventmodules.hh:486
static const H2GasPvt & h2GasPvt()
Definition: blackoilsolventmodules.hh:369
static Scalar tlMixParamViscosity(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:468
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 const TabulatedFunction & tlPMixTable(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:459
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 & sorwmis(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:441
static std::string primaryVarName(unsigned pvIdx)
Definition: blackoilsolventmodules.hh:142
static bool eqApplies(unsigned eqIdx)
Definition: blackoilsolventmodules.hh:157
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 const TabulatedFunction & misc(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:405
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: 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