28#ifndef EWOMS_BLACK_OIL_SOLVENT_MODULE_HH
29#define EWOMS_BLACK_OIL_SOLVENT_MODULE_HH
31#include <opm/common/Exceptions.hpp>
33#include <opm/material/fluidsystems/blackoilpvt/SolventPvt.hpp>
43#include <opm/material/common/Valgrind.hpp>
45#include <dune/common/fvector.hh>
56template <class TypeTag, bool enableSolventV = getPropValue<TypeTag, Properties::EnableSolvent>()>
71 using Toolbox = MathToolbox<Evaluation>;
80 static constexpr unsigned solventSaturationIdx = Indices::solventSaturationIdx;
81 static constexpr unsigned contiSolventEqIdx = Indices::contiSolventEqIdx;
82 static constexpr unsigned enableSolvent = enableSolventV;
83 static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
84 static constexpr unsigned numPhases = FluidSystem::numPhases;
85 static constexpr bool blackoilConserveSurfaceVolume = getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();
86 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
99 { params_.solventPvt_ = value; }
110 if constexpr (enableSolvent)
118 Simulator& simulator)
120 if constexpr (enableSolvent)
126 if constexpr (enableSolvent)
127 return pvIdx == solventSaturationIdx;
136 return "saturation_solvent";
144 return static_cast<Scalar
>(1.0);
149 if constexpr (enableSolvent)
150 return eqIdx == contiSolventEqIdx;
155 static std::string
eqName([[maybe_unused]]
unsigned eqIdx)
159 return "conti^solvent";
162 static Scalar
eqWeight([[maybe_unused]]
unsigned eqIdx)
167 return static_cast<Scalar
>(1.0);
170 template <
class LhsEval>
171 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
172 const IntensiveQuantities& intQuants)
174 if constexpr (enableSolvent) {
175 if constexpr (blackoilConserveSurfaceVolume) {
176 storage[contiSolventEqIdx] +=
177 Toolbox::template decay<LhsEval>(intQuants.porosity())
178 * Toolbox::template decay<LhsEval>(intQuants.solventSaturation())
179 * Toolbox::template decay<LhsEval>(intQuants.solventInverseFormationVolumeFactor());
181 storage[contiSolventEqIdx] += Toolbox::template decay<LhsEval>(intQuants.porosity())
182 * Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(waterPhaseIdx))
183 * Toolbox::template decay<LhsEval>(intQuants.fluidState().invB(waterPhaseIdx))
184 * Toolbox::template decay<LhsEval>(intQuants.rsSolw());
188 storage[contiSolventEqIdx] +=
189 Toolbox::template decay<LhsEval>(intQuants.porosity())
190 * Toolbox::template decay<LhsEval>(intQuants.solventSaturation())
191 * Toolbox::template decay<LhsEval>(intQuants.solventDensity());
193 storage[contiSolventEqIdx] += Toolbox::template decay<LhsEval>(intQuants.porosity())
194 * Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(waterPhaseIdx))
195 * Toolbox::template decay<LhsEval>(intQuants.fluidState().density(waterPhaseIdx))
196 * Toolbox::template decay<LhsEval>(intQuants.rsSolw());
204 [[maybe_unused]]
const ElementContext& elemCtx,
205 [[maybe_unused]]
unsigned scvfIdx,
206 [[maybe_unused]]
unsigned timeIdx)
209 if constexpr (enableSolvent) {
210 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
212 unsigned upIdx = extQuants.solventUpstreamIndex();
213 unsigned inIdx = extQuants.interiorIndex();
214 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
216 if constexpr (blackoilConserveSurfaceVolume) {
218 flux[contiSolventEqIdx] =
219 extQuants.solventVolumeFlux()
220 *up.solventInverseFormationVolumeFactor();
222 flux[contiSolventEqIdx] =
223 extQuants.solventVolumeFlux()
224 *decay<Scalar>(up.solventInverseFormationVolumeFactor());
229 flux[contiSolventEqIdx] +=
230 extQuants.volumeFlux(waterPhaseIdx)
231 * up.fluidState().invB(waterPhaseIdx)
234 flux[contiSolventEqIdx] +=
235 extQuants.volumeFlux(waterPhaseIdx)
236 *decay<Scalar>(up.fluidState().invB(waterPhaseIdx))
237 *decay<Scalar>(up.rsSolw());
242 flux[contiSolventEqIdx] =
243 extQuants.solventVolumeFlux()
244 *up.solventDensity();
246 flux[contiSolventEqIdx] =
247 extQuants.solventVolumeFlux()
248 *decay<Scalar>(up.solventDensity());
253 flux[contiSolventEqIdx] +=
254 extQuants.volumeFlux(waterPhaseIdx)
255 * up.fluidState().density(waterPhaseIdx)
258 flux[contiSolventEqIdx] +=
259 extQuants.volumeFlux(waterPhaseIdx)
260 *decay<Scalar>(up.fluidState().density(waterPhaseIdx))
261 *decay<Scalar>(up.rsSolw());
271 Scalar solventSaturation,
274 if constexpr (!enableSolvent) {
275 priVars.setPrimaryVarsMeaningSolvent(PrimaryVariables::SolventMeaning::Disabled);
280 priVars.setPrimaryVarsMeaningSolvent(PrimaryVariables::SolventMeaning::Ss);
281 priVars[solventSaturationIdx] = solventSaturation;
283 priVars.setPrimaryVarsMeaningSolvent(PrimaryVariables::SolventMeaning::Rsolw);
284 priVars[solventSaturationIdx] = solventRsw;
292 const PrimaryVariables& oldPv,
293 const EqVector& delta)
295 if constexpr (enableSolvent)
297 newPv[solventSaturationIdx] = oldPv[solventSaturationIdx] - delta[solventSaturationIdx];
309 return static_cast<Scalar
>(0.0);
318 return std::abs(Toolbox::scalarValue(resid[contiSolventEqIdx]));
321 template <
class DofEntity>
322 static void serializeEntity(
const Model& model, std::ostream& outstream,
const DofEntity& dof)
324 if constexpr (enableSolvent) {
325 unsigned dofIdx = model.dofMapper().index(dof);
327 const PrimaryVariables& priVars = model.solution(0)[dofIdx];
328 outstream << priVars[solventSaturationIdx];
332 template <
class DofEntity>
335 if constexpr (enableSolvent) {
336 unsigned dofIdx = model.dofMapper().index(dof);
338 PrimaryVariables& priVars0 = model.solution(0)[dofIdx];
339 PrimaryVariables& priVars1 = model.solution(1)[dofIdx];
341 instream >> priVars0[solventSaturationIdx];
344 priVars1 = priVars0[solventSaturationIdx];
350 return params_.solventPvt_;
356 return params_.co2GasPvt_;
361 return params_.h2GasPvt_;
366 return params_.brineCo2Pvt_;
371 return params_.brineH2Pvt_;
374 static const TabulatedFunction&
ssfnKrg(
const ElementContext& elemCtx,
378 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
379 return params_.ssfnKrg_[satnumRegionIdx];
382 static const TabulatedFunction&
ssfnKrs(
const ElementContext& elemCtx,
386 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
387 return params_.ssfnKrs_[satnumRegionIdx];
390 static const TabulatedFunction&
sof2Krn(
const ElementContext& elemCtx,
394 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
395 return params_.sof2Krn_[satnumRegionIdx];
398 static const TabulatedFunction&
misc(
const ElementContext& elemCtx,
402 unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
403 return params_.misc_[miscnumRegionIdx];
406 static const TabulatedFunction&
pmisc(
const ElementContext& elemCtx,
410 unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
411 return params_.pmisc_[miscnumRegionIdx];
414 static const TabulatedFunction&
msfnKrsg(
const ElementContext& elemCtx,
418 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
419 return params_.msfnKrsg_[satnumRegionIdx];
422 static const TabulatedFunction&
msfnKro(
const ElementContext& elemCtx,
426 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
427 return params_.msfnKro_[satnumRegionIdx];
430 static const TabulatedFunction&
sorwmis(
const ElementContext& elemCtx,
434 unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
435 return params_.sorwmis_[miscnumRegionIdx];
438 static const TabulatedFunction&
sgcwmis(
const ElementContext& elemCtx,
442 unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
443 return params_.sgcwmis_[miscnumRegionIdx];
446 static const TabulatedFunction&
tlPMixTable(
const ElementContext& elemCtx,
450 unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
451 return params_.tlPMixTable_[miscnumRegionIdx];
458 unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
459 return params_.tlMixParamViscosity_[miscnumRegionIdx];
466 unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
467 return params_.tlMixParamDensity_[miscnumRegionIdx];
472 return params_.isMiscible_;
475 template <
class Value>
476 static const Value
solubilityLimit(
unsigned pvtIdx,
const Value& temperature,
const Value& pressure,
const Value& saltConcentration)
483 return brineCo2Pvt().saturatedGasDissolutionFactor(pvtIdx, temperature, pressure, saltConcentration);
485 return brineH2Pvt().saturatedGasDissolutionFactor(pvtIdx, temperature, pressure, saltConcentration);
490 return params_.rsSolw_active_;
495 return params_.co2sol_;
500 return params_.h2sol_;
507template <
class TypeTag,
bool enableSolventV>
508BlackOilSolventParams<typename BlackOilSolventModule<TypeTag, enableSolventV>::Scalar>
509BlackOilSolventModule<TypeTag, enableSolventV>::params_;
518template <class TypeTag, bool enableSolventV = getPropValue<TypeTag, Properties::EnableSolvent>()>
533 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
534 static constexpr int solventSaturationIdx = Indices::solventSaturationIdx;
535 static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
536 static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
537 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
538 static constexpr double cutOff = 1e-12;
552 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
554 auto& fs =
asImp_().fluidState_;
556 if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
557 solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx, elemCtx.linearizationType());
584 auto& fs =
asImp_().fluidState_;
590 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
591 if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
593 }
else if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Rsolw) {
594 rsSolw_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx, elemCtx.linearizationType());
605 const Evaluation& p = FluidSystem::phaseIsActive(oilPhaseIdx)? fs.pressure(oilPhaseIdx) : fs.pressure(gasPhaseIdx);
607 const Evaluation& pgImisc = fs.pressure(gasPhaseIdx);
610 const auto& problem = elemCtx.problem();
611 Evaluation pgMisc = 0.0;
612 std::array<Evaluation, numPhases> pC;
613 const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, timeIdx);
614 MaterialLaw::capillaryPressures(pC, materialParams, fs);
617 const auto linearizationType = elemCtx.linearizationType();
618 if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pg)
619 pgMisc = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, linearizationType);
621 const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, linearizationType);
622 pgMisc = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
625 fs.setPressure(gasPhaseIdx, pmisc * pgMisc + (1.0 - pmisc) * pgImisc);
631 if (gasSolventSat.value() < cutOff)
641 const Evaluation& p = FluidSystem::phaseIsActive(oilPhaseIdx)? fs.pressure(oilPhaseIdx) : fs.pressure(gasPhaseIdx);
642 const Evaluation miscibility = misc.eval(Fsolgas,
true) * pmisc.eval(p,
true);
645 unsigned cellIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
646 const auto& materialLawManager = elemCtx.problem().materialLawManager();
647 const auto& scaledDrainageInfo =
648 materialLawManager->oilWaterScaledEpsInfoDrainage(cellIdx);
650 const Scalar& sogcr = scaledDrainageInfo.Sogcr;
651 Evaluation sor = sogcr;
652 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
653 const Evaluation& sw = fs.saturation(waterPhaseIdx);
655 sor = miscibility * sorwmis.eval(sw,
true) + (1.0 - miscibility) * sogcr;
657 const Scalar& sgcr = scaledDrainageInfo.Sgcr;
658 Evaluation sgc = sgcr;
659 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
660 const Evaluation& sw = fs.saturation(waterPhaseIdx);
662 sgc = miscibility * sgcwmis.eval(sw,
true) + (1.0 - miscibility) * sgcr;
665 Evaluation oilGasSolventSat = gasSolventSat;
666 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
667 oilGasSolventSat += fs.saturation(oilPhaseIdx);
669 const Evaluation zero = 0.0;
670 const Evaluation oilGasSolventEffSat = std::max(oilGasSolventSat - sor - sgc, zero);
672 Evaluation F_totalGas = 0.0;
673 if (oilGasSolventEffSat.value() > cutOff) {
674 const Evaluation gasSolventEffSat = std::max(gasSolventSat - sgc, zero);
675 F_totalGas = gasSolventEffSat / oilGasSolventEffSat;
681 const Evaluation mkrgt = msfnKrsg.eval(F_totalGas,
true) * sof2Krn.eval(oilGasSolventSat,
true);
682 const Evaluation mkro = msfnKro.eval(F_totalGas,
true) * sof2Krn.eval(oilGasSolventSat,
true);
684 Evaluation& kro =
asImp_().mobility_[oilPhaseIdx];
685 Evaluation& krg =
asImp_().mobility_[gasPhaseIdx];
688 krg *= (1.0 - miscibility);
689 krg += miscibility * mkrgt;
690 kro *= (1.0 - miscibility);
691 kro += miscibility * mkro;
698 Evaluation& krg =
asImp_().mobility_[gasPhaseIdx];
700 krg *= ssfnKrg.eval(Fhydgas,
true);
714 const auto& iq =
asImp_();
715 unsigned pvtRegionIdx = iq.pvtRegionIndex();
716 const Evaluation& T = iq.fluidState().temperature(gasPhaseIdx);
717 const Evaluation& p = iq.fluidState().pressure(gasPhaseIdx);
719 const Evaluation rv = 0.0;
720 const Evaluation rvw = 0.0;
729 auto& fs =
asImp_().fluidState_;
730 const auto& bw = brineCo2Pvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p,
rsSolw());
732 const auto denw = bw*brineCo2Pvt.waterReferenceDensity(pvtRegionIdx)
733 +
rsSolw()*bw*brineCo2Pvt.gasReferenceDensity(pvtRegionIdx);
734 fs.setDensity(waterPhaseIdx, denw);
735 fs.setInvB(waterPhaseIdx, bw);
736 Evaluation& mobw =
asImp_().mobility_[waterPhaseIdx];
737 const auto& muWat = fs.viscosity(waterPhaseIdx);
738 const auto& muWatEff = brineCo2Pvt.viscosity(pvtRegionIdx, T, p,
rsSolw());
739 mobw *= muWat / muWatEff;
747 auto& fs =
asImp_().fluidState_;
748 const auto& bw = brineH2Pvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p,
rsSolw());
750 const auto denw = bw*brineH2Pvt.waterReferenceDensity(pvtRegionIdx)
751 +
rsSolw()*bw*brineH2Pvt.gasReferenceDensity(pvtRegionIdx);
752 fs.setDensity(waterPhaseIdx, denw);
753 fs.setInvB(waterPhaseIdx, bw);
754 Evaluation& mobw =
asImp_().mobility_[waterPhaseIdx];
755 const auto& muWat = fs.viscosity(waterPhaseIdx);
756 const auto& muWatEff = brineH2Pvt.viscosity(pvtRegionIdx, T, p,
rsSolw());
757 mobw *= muWat / muWatEff;
767 effectiveProperties(elemCtx, scvIdx, timeIdx);
799 void effectiveProperties(
const ElementContext& elemCtx,
813 auto& fs =
asImp_().fluidState_;
816 Evaluation oilEffSat = 0.0;
817 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
818 oilEffSat = fs.saturation(oilPhaseIdx);
820 Evaluation gasEffSat = fs.saturation(gasPhaseIdx);
822 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
825 const Evaluation zero = 0.0;
826 const Evaluation& sw = fs.saturation(waterPhaseIdx);
827 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
828 oilEffSat = std::max(oilEffSat - sorwmis.eval(sw,
true), zero);
830 gasEffSat = std::max(gasEffSat - sgcwmis.eval(sw,
true), zero);
831 solventEffSat = std::max(solventEffSat - sgcwmis.eval(sw,
true), zero);
833 const Evaluation oilGasSolventEffSat = oilEffSat + gasEffSat + solventEffSat;
834 const Evaluation oilSolventEffSat = oilEffSat + solventEffSat;
835 const Evaluation solventGasEffSat = solventEffSat + gasEffSat;
842 const Evaluation& p = FluidSystem::phaseIsActive(oilPhaseIdx)? fs.pressure(oilPhaseIdx) : fs.pressure(gasPhaseIdx);
845 const Evaluation pmisc = pmiscTable.eval(p,
true);
849 const Evaluation& muGas = fs.viscosity(gasPhaseIdx);
852 assert(muGas.value() > 0);
853 assert(muSolvent.value() > 0);
854 const Evaluation muGasPow = pow(muGas, 0.25);
855 const Evaluation muSolventPow = pow(muSolvent, 0.25);
857 Evaluation muMixSolventGas = muGas;
858 if (solventGasEffSat > cutOff)
859 muMixSolventGas *= muSolvent / pow(((gasEffSat / solventGasEffSat) * muSolventPow) + ((solventEffSat / solventGasEffSat) * muGasPow) , 4.0);
861 Evaluation muOil = 1.0;
862 Evaluation muOilPow = 1.0;
863 Evaluation muMixOilSolvent = 1.0;
864 Evaluation muOilEff = 1.0;
865 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
866 muOil = fs.viscosity(oilPhaseIdx);
867 assert(muOil.value() > 0);
868 muOilPow = pow(muOil, 0.25);
869 muMixOilSolvent = muOil;
870 if (oilSolventEffSat > cutOff)
871 muMixOilSolvent *= muSolvent / pow(((oilEffSat / oilSolventEffSat) * muSolventPow) + ((solventEffSat / oilSolventEffSat) * muOilPow) , 4.0);
873 muOilEff = pow(muOil,1.0 - tlMixParamMu) * pow(muMixOilSolvent, tlMixParamMu);
875 Evaluation muMixSolventGasOil = muOil;
876 if (oilGasSolventEffSat > cutOff)
877 muMixSolventGasOil *= muSolvent * muGas / pow(((oilEffSat / oilGasSolventEffSat) * muSolventPow * muGasPow)
878 + ((solventEffSat / oilGasSolventEffSat) * muOilPow * muGasPow) + ((gasEffSat / oilGasSolventEffSat) * muSolventPow * muOilPow), 4.0);
880 Evaluation muGasEff = pow(muGas,1.0 - tlMixParamMu) * pow(muMixSolventGas, tlMixParamMu);
881 Evaluation muSolventEff = pow(muSolvent,1.0 - tlMixParamMu) * pow(muMixSolventGasOil, tlMixParamMu);
884 const Evaluation& rhoGas = fs.density(gasPhaseIdx);
885 Evaluation rhoOil = 0.0;
886 if (FluidSystem::phaseIsActive(oilPhaseIdx))
887 rhoOil = fs.density(oilPhaseIdx);
898 const Evaluation muOilEffPow = pow(pow(muOil, 1.0 - tlMixParamRho) * pow(muMixOilSolvent, tlMixParamRho), 0.25);
899 const Evaluation muGasEffPow = pow(pow(muGas, 1.0 - tlMixParamRho) * pow(muMixSolventGas, tlMixParamRho), 0.25);
900 const Evaluation muSolventEffPow = pow(pow(muSolvent, 1.0 - tlMixParamRho) * pow(muMixSolventGasOil, tlMixParamRho), 0.25);
902 const Evaluation oilGasEffSaturation = oilEffSat + gasEffSat;
903 Evaluation sof = 0.0;
904 Evaluation sgf = 0.0;
905 if (oilGasEffSaturation.value() > cutOff) {
906 sof = oilEffSat / oilGasEffSaturation;
907 sgf = gasEffSat / oilGasEffSaturation;
910 const Evaluation muSolventOilGasPow = muSolventPow * ((sgf * muOilPow) + (sof * muGasPow));
912 Evaluation rhoMixSolventGasOil = 0.0;
913 if (oilGasSolventEffSat.value() > cutOff)
914 rhoMixSolventGasOil = (rhoOil * oilEffSat / oilGasSolventEffSat) + (rhoGas * gasEffSat / oilGasSolventEffSat) + (rhoSolvent * solventEffSat / oilGasSolventEffSat);
916 Evaluation rhoGasEff = 0.0;
917 if (std::abs(muSolventPow.value() - muGasPow.value()) < cutOff)
918 rhoGasEff = ((1.0 - tlMixParamRho) * rhoGas) + (tlMixParamRho * rhoMixSolventGasOil);
920 const Evaluation solventGasEffFraction = (muGasPow * (muSolventPow - muGasEffPow)) / (muGasEffPow * (muSolventPow - muGasPow));
921 rhoGasEff = (rhoGas * solventGasEffFraction) + (rhoSolvent * (1.0 - solventGasEffFraction));
924 Evaluation rhoSolventEff = 0.0;
925 if (std::abs((muSolventOilGasPow.value() - (muOilPow.value() * muGasPow.value()))) < cutOff)
926 rhoSolventEff = ((1.0 - tlMixParamRho) * rhoSolvent) + (tlMixParamRho * rhoMixSolventGasOil);
928 const Evaluation sfraction_se = (muSolventOilGasPow - (muOilPow * muGasPow * muSolventPow / muSolventEffPow)) / (muSolventOilGasPow - (muOilPow * muGasPow));
929 rhoSolventEff = (rhoSolvent * sfraction_se) + (rhoGas * sgf * (1.0 - sfraction_se)) + (rhoOil * sof * (1.0 - sfraction_se));
933 unsigned pvtRegionIdx =
asImp_().pvtRegionIndex();
934 Evaluation bGasEff = rhoGasEff;
935 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
936 bGasEff /= (FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) + FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) * fs.Rv());
938 bGasEff /= (FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx));
943 const Evaluation bg = fs.invB(gasPhaseIdx);
945 const Evaluation bg_eff = pmisc * bGasEff + (1.0 - pmisc) * bg;
946 const Evaluation bs_eff = pmisc * bSolventEff + (1.0 - pmisc) * bs;
949 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
950 fs.setDensity(gasPhaseIdx,
952 *(FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx)
953 + FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx)*fs.Rv()));
955 fs.setDensity(gasPhaseIdx,
957 *FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx));
962 Evaluation& mobg =
asImp_().mobility_[gasPhaseIdx];
963 muGasEff = bg_eff / (pmisc * bGasEff / muGasEff + (1.0 - pmisc) * bg / muGas);
964 mobg *= muGas / muGasEff;
967 solventViscosity_ = bs_eff / (pmisc * bSolventEff / muSolventEff + (1.0 - pmisc) * bs / muSolvent);
969 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
970 Evaluation rhoOilEff = 0.0;
971 if (std::abs(muOilPow.value() - muSolventPow.value()) < cutOff) {
972 rhoOilEff = ((1.0 - tlMixParamRho) * rhoOil) + (tlMixParamRho * rhoMixSolventGasOil);
975 const Evaluation solventOilEffFraction = (muOilPow * (muOilEffPow - muSolventPow)) / (muOilEffPow * (muOilPow - muSolventPow));
976 rhoOilEff = (rhoOil * solventOilEffFraction) + (rhoSolvent * (1.0 - solventOilEffFraction));
978 const Evaluation bOilEff = rhoOilEff / (FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) + FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) * fs.Rs());
979 const Evaluation bo = fs.invB(oilPhaseIdx);
980 const Evaluation bo_eff = pmisc * bOilEff + (1.0 - pmisc) * bo;
981 fs.setDensity(oilPhaseIdx,
983 *(FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx)
984 + FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx)*fs.Rs()));
987 Evaluation& mobo =
asImp_().mobility_[oilPhaseIdx];
988 muOilEff = bo_eff / (pmisc * bOilEff / muOilEff + (1.0 - pmisc) * bo / muOil);
989 mobo *= muOil / muOilEff;
995 {
return *
static_cast<Implementation*
>(
this); }
1008template <
class TypeTag>
1033 {
throw std::runtime_error(
"solventSaturation() called but solvents are disabled"); }
1036 {
throw std::runtime_error(
"rsSolw() called but solvents are disabled"); }
1039 {
throw std::runtime_error(
"solventDensity() called but solvents are disabled"); }
1042 {
throw std::runtime_error(
"solventViscosity() called but solvents are disabled"); }
1045 {
throw std::runtime_error(
"solventMobility() called but solvents are disabled"); }
1048 {
throw std::runtime_error(
"solventInverseFormationVolumeFactor() called but solvents are disabled"); }
1051 {
throw std::runtime_error(
"solventRefDensity() called but solvents are disabled"); }
1061template <class TypeTag, bool enableSolventV = getPropValue<TypeTag, Properties::EnableSolvent>()>
1074 using Toolbox = MathToolbox<Evaluation>;
1076 static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
1077 static constexpr int dimWorld = GridView::dimensionworld;
1079 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
1080 using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
1087 template <
class Dummy =
bool>
1093 const auto& gradCalc = elemCtx.gradientCalculator();
1096 const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
1097 const auto& faceNormal = scvf.normal();
1099 unsigned i = scvf.interiorIndex();
1100 unsigned j = scvf.exteriorIndex();
1103 DimEvalVector solventPGrad;
1105 gradCalc.calculateGradient(solventPGrad,
1109 Valgrind::CheckDefined(solventPGrad);
1112 if (Parameters::Get<Parameters::EnableGravity>()) {
1115 const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
1116 const auto& gEx = elemCtx.problem().gravity(elemCtx, j, timeIdx);
1118 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
1119 const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
1121 const auto& posIn = elemCtx.pos(i, timeIdx);
1122 const auto& posEx = elemCtx.pos(j, timeIdx);
1123 const auto& posFace = scvf.integrationPos();
1126 DimVector distVecIn(posIn);
1127 DimVector distVecEx(posEx);
1128 DimVector distVecTotal(posEx);
1130 distVecIn -= posFace;
1131 distVecEx -= posFace;
1132 distVecTotal -= posIn;
1133 Scalar absDistTotalSquared = distVecTotal.two_norm2();
1136 auto rhoIn = intQuantsIn.solventDensity();
1137 auto pStatIn = - rhoIn*(gIn*distVecIn);
1141 Scalar rhoEx = Toolbox::value(intQuantsEx.solventDensity());
1142 Scalar pStatEx = - rhoEx*(gEx*distVecEx);
1148 DimEvalVector f(distVecTotal);
1149 f *= (pStatEx - pStatIn)/absDistTotalSquared;
1152 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
1153 solventPGrad[dimIdx] += f[dimIdx];
1155 if (!isfinite(solventPGrad[dimIdx]))
1156 throw NumericalProblem(
"Non-finite potential gradient for solvent 'phase'");
1161 Evaluation solventPGradNormal = 0.0;
1162 for (
unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx)
1163 solventPGradNormal += solventPGrad[dimIdx]*faceNormal[dimIdx];
1165 if (solventPGradNormal > 0) {
1166 solventUpstreamDofIdx_ = j;
1167 solventDownstreamDofIdx_ = i;
1170 solventUpstreamDofIdx_ = i;
1171 solventDownstreamDofIdx_ = j;
1174 const auto& up = elemCtx.intensiveQuantities(solventUpstreamDofIdx_, timeIdx);
1180 if (solventUpstreamDofIdx_ == i)
1181 solventVolumeFlux_ = solventPGradNormal*up.solventMobility();
1183 solventVolumeFlux_ = solventPGradNormal*scalarValue(up.solventMobility());
1190 template <
class Dummy =
bool>
1196 const ExtensiveQuantities& extQuants = asImp_();
1198 unsigned interiorDofIdx = extQuants.interiorIndex();
1199 unsigned exteriorDofIdx = extQuants.exteriorIndex();
1200 assert(interiorDofIdx != exteriorDofIdx);
1202 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
1203 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
1205 unsigned I = elemCtx.globalSpaceIndex(interiorDofIdx, timeIdx);
1206 unsigned J = elemCtx.globalSpaceIndex(exteriorDofIdx, timeIdx);
1208 Scalar thpres = elemCtx.problem().thresholdPressure(I, J);
1209 Scalar trans = elemCtx.problem().transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
1210 Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
1212 Scalar zIn = elemCtx.problem().dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
1213 Scalar zEx = elemCtx.problem().dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
1214 Scalar distZ = zIn - zEx;
1216 const Evaluation& rhoIn = intQuantsIn.solventDensity();
1217 Scalar rhoEx = Toolbox::value(intQuantsEx.solventDensity());
1218 const Evaluation& rhoAvg = rhoIn*0.5 + rhoEx*0.5;
1220 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(gasPhaseIdx);
1221 Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(gasPhaseIdx));
1222 pressureExterior += distZ*g*rhoAvg;
1224 Evaluation pressureDiffSolvent = pressureExterior - pressureInterior;
1225 if (std::abs(scalarValue(pressureDiffSolvent)) > thpres) {
1226 if (pressureDiffSolvent < 0.0)
1227 pressureDiffSolvent += thpres;
1229 pressureDiffSolvent -= thpres;
1232 pressureDiffSolvent = 0.0;
1234 if (pressureDiffSolvent > 0.0) {
1235 solventUpstreamDofIdx_ = exteriorDofIdx;
1236 solventDownstreamDofIdx_ = interiorDofIdx;
1238 else if (pressureDiffSolvent < 0.0) {
1239 solventUpstreamDofIdx_ = interiorDofIdx;
1240 solventDownstreamDofIdx_ = exteriorDofIdx;
1246 solventUpstreamDofIdx_ = std::min(interiorDofIdx, exteriorDofIdx);
1247 solventDownstreamDofIdx_ = std::max(interiorDofIdx, exteriorDofIdx);
1248 solventVolumeFlux_ = 0.0;
1252 Scalar faceArea = elemCtx.stencil(timeIdx).interiorFace(scvfIdx).area();
1253 const IntensiveQuantities& up = elemCtx.intensiveQuantities(solventUpstreamDofIdx_, timeIdx);
1254 if (solventUpstreamDofIdx_ == interiorDofIdx)
1255 solventVolumeFlux_ =
1256 up.solventMobility()
1258 *pressureDiffSolvent;
1260 solventVolumeFlux_ =
1261 scalarValue(up.solventMobility())
1263 *pressureDiffSolvent;
1267 {
return solventUpstreamDofIdx_; }
1270 {
return solventDownstreamDofIdx_; }
1273 {
return solventVolumeFlux_; }
1279 Implementation& asImp_()
1280 {
return *
static_cast<Implementation*
>(
this); }
1282 Evaluation solventVolumeFlux_;
1283 unsigned solventUpstreamDofIdx_;
1284 unsigned solventDownstreamDofIdx_;
1287template <
class TypeTag>
1305 {
throw std::runtime_error(
"solventUpstreamIndex() called but solvents are disabled"); }
1308 {
throw std::runtime_error(
"solventDownstreamIndex() called but solvents are disabled"); }
1311 {
throw std::runtime_error(
"solventVolumeFlux() called but solvents are disabled"); }
1314 {
throw std::runtime_error(
"setSolventVolumeFlux() called but solvents are disabled"); }
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:1304
unsigned solventDownstreamIndex() const
Definition: blackoilsolventmodules.hh:1307
void setSolventVolumeFlux(const Evaluation &)
Definition: blackoilsolventmodules.hh:1313
const Evaluation & solventVolumeFlux() const
Definition: blackoilsolventmodules.hh:1310
void updateVolumeFluxPerm(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1294
void updateVolumeFluxTrans(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1299
Provides the solvent specific extensive quantities to the generic black-oil module's extensive quanti...
Definition: blackoilsolventmodules.hh:1063
unsigned solventUpstreamIndex() const
Definition: blackoilsolventmodules.hh:1266
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:1089
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:1192
unsigned solventDownstreamIndex() const
Definition: blackoilsolventmodules.hh:1269
const Evaluation & solventVolumeFlux() const
Definition: blackoilsolventmodules.hh:1272
void setSolventVolumeFlux(const Evaluation &solventVolumeFlux)
Definition: blackoilsolventmodules.hh:1275
const Scalar & solventRefDensity() const
Definition: blackoilsolventmodules.hh:1050
const Evaluation & solventDensity() const
Definition: blackoilsolventmodules.hh:1038
const Evaluation & solventInverseFormationVolumeFactor() const
Definition: blackoilsolventmodules.hh:1047
const Evaluation & rsSolw() const
Definition: blackoilsolventmodules.hh:1035
void solventPreSatFuncUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1017
const Evaluation & solventViscosity() const
Definition: blackoilsolventmodules.hh:1041
const Evaluation & solventMobility() const
Definition: blackoilsolventmodules.hh:1044
const Evaluation & solventSaturation() const
Definition: blackoilsolventmodules.hh:1032
void solventPvtUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1027
void solventPostSatFuncUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilsolventmodules.hh:1022
Provides the volumetric quantities required for the equations needed by the solvents extension of the...
Definition: blackoilsolventmodules.hh:520
void solventPreSatFuncUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Called before the saturation functions are doing their magic.
Definition: blackoilsolventmodules.hh:548
const Evaluation & solventMobility() const
Definition: blackoilsolventmodules.hh:786
Evaluation solventViscosity_
Definition: blackoilsolventmodules.hh:1001
Evaluation solventMobility_
Definition: blackoilsolventmodules.hh:1002
Scalar solventRefDensity_
Definition: blackoilsolventmodules.hh:1005
Evaluation solventDensity_
Definition: blackoilsolventmodules.hh:1000
Evaluation rsSolw_
Definition: blackoilsolventmodules.hh:999
Implementation & asImp_()
Definition: blackoilsolventmodules.hh:994
Evaluation solventSaturation_
Definition: blackoilsolventmodules.hh:998
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:710
const Evaluation & solventDensity() const
Definition: blackoilsolventmodules.hh:780
Evaluation solventInvFormationVolumeFactor_
Definition: blackoilsolventmodules.hh:1003
const Evaluation & rsSolw() const
Definition: blackoilsolventmodules.hh:777
const Evaluation & solventViscosity() const
Definition: blackoilsolventmodules.hh:783
void solventPostSatFuncUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Called after the saturation functions have been doing their magic.
Definition: blackoilsolventmodules.hh:578
const Evaluation & solventInverseFormationVolumeFactor() const
Definition: blackoilsolventmodules.hh:789
const Evaluation & solventSaturation() const
Definition: blackoilsolventmodules.hh:774
const Scalar & solventRefDensity() const
Definition: blackoilsolventmodules.hh:793
Evaluation hydrocarbonSaturation_
Definition: blackoilsolventmodules.hh:997
Contains the high level supplements required to extend the black oil model by solvents.
Definition: blackoilsolventmodules.hh:58
static void setSolventPvt(const SolventPvt &value)
Specify the solvent PVT of a all PVT regions.
Definition: blackoilsolventmodules.hh:98
static const TabulatedFunction & sof2Krn(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:390
static const TabulatedFunction & ssfnKrs(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:382
static const TabulatedFunction & pmisc(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:406
static const TabulatedFunction & msfnKrsg(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:414
static const TabulatedFunction & sgcwmis(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:438
static Scalar computeResidualError(const EqVector &resid)
Return how much a residual is considered an error.
Definition: blackoilsolventmodules.hh:315
static const BrineH2Pvt & brineH2Pvt()
Definition: blackoilsolventmodules.hh:369
static Scalar primaryVarWeight(unsigned pvIdx)
Definition: blackoilsolventmodules.hh:139
static void setParams(BlackOilSolventParams< Scalar > &¶ms)
Set parameters.
Definition: blackoilsolventmodules.hh:90
static std::string eqName(unsigned eqIdx)
Definition: blackoilsolventmodules.hh:155
static Scalar eqWeight(unsigned eqIdx)
Definition: blackoilsolventmodules.hh:162
static void serializeEntity(const Model &model, std::ostream &outstream, const DofEntity &dof)
Definition: blackoilsolventmodules.hh:322
static bool isSolubleInWater()
Definition: blackoilsolventmodules.hh:488
static bool isCO2Sol()
Definition: blackoilsolventmodules.hh:493
static const TabulatedFunction & msfnKro(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:422
static void registerOutputModules(Model &model, Simulator &simulator)
Register all solvent specific VTK and ECL output modules.
Definition: blackoilsolventmodules.hh:117
static const TabulatedFunction & ssfnKrg(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:374
static void setIsMiscible(const bool isMiscible)
Definition: blackoilsolventmodules.hh:102
static void registerParameters()
Register all run-time parameters for the black-oil solvent module.
Definition: blackoilsolventmodules.hh:108
static Scalar computeUpdateError(const PrimaryVariables &, const EqVector &)
Return how much a Newton-Raphson update is considered an error.
Definition: blackoilsolventmodules.hh:303
static const SolventPvt & solventPvt()
Definition: blackoilsolventmodules.hh:348
static bool isH2Sol()
Definition: blackoilsolventmodules.hh:498
static bool isMiscible()
Definition: blackoilsolventmodules.hh:470
static const H2GasPvt & h2GasPvt()
Definition: blackoilsolventmodules.hh:359
static const Scalar & tlMixParamDensity(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:462
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:291
static const TabulatedFunction & tlPMixTable(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:446
static const Scalar & tlMixParamViscosity(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:454
static void assignPrimaryVars(PrimaryVariables &priVars, Scalar solventSaturation, Scalar solventRsw)
Assign the solvent specific primary variables to a PrimaryVariables object.
Definition: blackoilsolventmodules.hh:270
static const TabulatedFunction & sorwmis(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:430
static std::string primaryVarName(unsigned pvIdx)
Definition: blackoilsolventmodules.hh:132
static bool eqApplies(unsigned eqIdx)
Definition: blackoilsolventmodules.hh:147
static const Value solubilityLimit(unsigned pvtIdx, const Value &temperature, const Value &pressure, const Value &saltConcentration)
Definition: blackoilsolventmodules.hh:476
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilsolventmodules.hh:171
static void deserializeEntity(Model &model, std::istream &instream, const DofEntity &dof)
Definition: blackoilsolventmodules.hh:333
static const Co2GasPvt & co2GasPvt()
Definition: blackoilsolventmodules.hh:354
static const BrineCo2Pvt & brineCo2Pvt()
Definition: blackoilsolventmodules.hh:364
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:203
static const TabulatedFunction & misc(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:398
static bool primaryVarApplies(unsigned pvIdx)
Definition: blackoilsolventmodules.hh:124
Callback class for a phase pressure.
Definition: quantitycallbacks.hh:84
void setPhaseIndex(unsigned phaseIdx)
Set the index of the fluid phase for which the pressure should be returned.
Definition: quantitycallbacks.hh:108
VTK output module for the black oil model's solvent related quantities.
Definition: vtkblackoilsolventmodule.hpp:54
Defines the common parameters for the porous medium multi-phase models.
Definition: blackoilboundaryratevector.hh:37
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:235
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