72 using SolventPvt =
typename BlackOilSolventParams<Scalar>::SolventPvt;
73 using Co2GasPvt =
typename BlackOilSolventParams<Scalar>::Co2GasPvt;
74 using H2GasPvt =
typename BlackOilSolventParams<Scalar>::H2GasPvt;
75 using BrineCo2Pvt =
typename BlackOilSolventParams<Scalar>::BrineCo2Pvt;
76 using BrineH2Pvt =
typename BlackOilSolventParams<Scalar>::BrineH2Pvt;
78 using TabulatedFunction =
typename BlackOilSolventParams<Scalar>::TabulatedFunction;
80 static constexpr unsigned solventSaturationIdx = Indices::solventSaturationIdx;
81 static constexpr unsigned contiSolventEqIdx = Indices::contiSolventEqIdx;
84 static constexpr unsigned numPhases = FluidSystem::numPhases;
86 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
99 { params_.solventPvt_ = value; }
102 static void setIsMiscible(
const bool isMiscible)
103 { params_.isMiscible_ = isMiscible; }
110 if constexpr (enableSolvent)
118 Simulator& simulator)
120 if constexpr (enableSolvent)
124 static bool primaryVarApplies(
unsigned pvIdx)
126 if constexpr (enableSolvent)
127 return pvIdx == solventSaturationIdx;
136 return "saturation_solvent";
144 return static_cast<Scalar
>(1.0);
147 static bool eqApplies(
unsigned eqIdx)
149 if constexpr (enableSolvent)
150 return eqIdx == contiSolventEqIdx;
159 return "conti^solvent";
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) {
179 * Toolbox::template
decay<LhsEval>(intQuants.solventInverseFormationVolumeFactor());
180 if (isSolubleInWater()) {
182 * Toolbox::template
decay<LhsEval>(intQuants.fluidState().saturation(waterPhaseIdx))
183 * Toolbox::template
decay<LhsEval>(intQuants.fluidState().invB(waterPhaseIdx))
192 if (isSolubleInWater()) {
194 * Toolbox::template
decay<LhsEval>(intQuants.fluidState().saturation(waterPhaseIdx))
195 * Toolbox::template
decay<LhsEval>(intQuants.fluidState().density(waterPhaseIdx))
209 if constexpr (enableSolvent) {
216 if constexpr (blackoilConserveSurfaceVolume) {
218 flux[contiSolventEqIdx] =
220 *
up.solventInverseFormationVolumeFactor();
222 flux[contiSolventEqIdx] =
227 if (isSolubleInWater()) {
229 flux[contiSolventEqIdx] +=
231 *
up.fluidState().invB(waterPhaseIdx)
234 flux[contiSolventEqIdx] +=
242 flux[contiSolventEqIdx] =
244 *
up.solventDensity();
246 flux[contiSolventEqIdx] =
251 if (isSolubleInWater()) {
253 flux[contiSolventEqIdx] +=
255 *
up.fluidState().density(waterPhaseIdx)
258 flux[contiSolventEqIdx] +=
271 Scalar solventSaturation,
274 if constexpr (!enableSolvent) {
275 priVars.setPrimaryVarsMeaningSolvent(PrimaryVariables::SolventMeaning::Disabled);
279 if (solventSaturation > 0 || !isSolubleInWater()) {
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>
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];
348 static const SolventPvt& solventPvt()
350 return params_.solventPvt_;
354 static const Co2GasPvt& co2GasPvt()
356 return params_.co2GasPvt_;
359 static const H2GasPvt& h2GasPvt()
361 return params_.h2GasPvt_;
364 static const BrineCo2Pvt& brineCo2Pvt()
366 return params_.brineCo2Pvt_;
369 static const BrineH2Pvt& brineH2Pvt()
371 return params_.brineH2Pvt_;
374 static const TabulatedFunction& ssfnKrg(
const ElementContext& elemCtx,
382 static const TabulatedFunction& ssfnKrs(
const ElementContext& elemCtx,
390 static const TabulatedFunction& sof2Krn(
const ElementContext& elemCtx,
398 static const TabulatedFunction& misc(
const ElementContext& elemCtx,
406 static const TabulatedFunction& pmisc(
const ElementContext& elemCtx,
414 static const TabulatedFunction& msfnKrsg(
const ElementContext& elemCtx,
422 static const TabulatedFunction& msfnKro(
const ElementContext& elemCtx,
430 static const TabulatedFunction& sorwmis(
const ElementContext& elemCtx,
438 static const TabulatedFunction& sgcwmis(
const ElementContext& elemCtx,
446 static const TabulatedFunction& tlPMixTable(
const ElementContext& elemCtx,
454 static const Scalar& tlMixParamViscosity(
const ElementContext& elemCtx,
462 static const Scalar& tlMixParamDensity(
const ElementContext& elemCtx,
470 static bool isMiscible()
472 return params_.isMiscible_;
475 template <
class Value>
476 static const Value solubilityLimit(
unsigned pvtIdx,
const Value& temperature,
const Value& pressure,
const Value& saltConcentration)
478 if (!isSolubleInWater())
481 assert(isCO2Sol() || isH2Sol());
483 return brineCo2Pvt().saturatedGasDissolutionFactor(pvtIdx, temperature, pressure, saltConcentration);
485 return brineH2Pvt().saturatedGasDissolutionFactor(pvtIdx, temperature, pressure, saltConcentration);
488 static bool isSolubleInWater()
490 return params_.rsSolw_active_;
493 static bool isCO2Sol()
495 return params_.co2sol_;
498 static bool isH2Sol()
500 return params_.h2sol_;
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 = 1
e-12;
552 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx,
timeIdx);
554 auto& fs = asImp_().fluidState_;
555 solventSaturation_ = 0.0;
556 if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
557 solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx,
timeIdx, elemCtx.linearizationType());
560 hydrocarbonSaturation_ = fs.saturation(gasPhaseIdx);
563 if (solventSaturation().value() < cutOff)
568 fs.setSaturation(gasPhaseIdx, hydrocarbonSaturation_ + solventSaturation_);
584 auto& fs = asImp_().fluidState_;
585 fs.setSaturation(gasPhaseIdx, hydrocarbonSaturation_);
590 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx,
timeIdx);
591 if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
592 rsSolw_ = SolventModule::solubilityLimit(asImp_().pvtRegionIndex(), fs.temperature(waterPhaseIdx), fs.pressure(waterPhaseIdx), fs.saltConcentration());
593 }
else if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Rsolw) {
594 rsSolw_ = priVars.makeEvaluation(solventSaturationIdx,
timeIdx, elemCtx.linearizationType());
597 solventMobility_ = 0.0;
600 if (solventSaturation().value() < cutOff)
604 if (SolventModule::isMiscible()) {
605 const Evaluation&
p = FluidSystem::phaseIsActive(oilPhaseIdx)? fs.pressure(oilPhaseIdx) : fs.pressure(gasPhaseIdx);
606 const Evaluation pmisc = SolventModule::pmisc(elemCtx, dofIdx,
timeIdx).eval(
p,
true);
607 const Evaluation&
pgImisc = fs.pressure(gasPhaseIdx);
610 const auto& problem = elemCtx.problem();
612 std::array<Evaluation, numPhases>
pC;
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);
629 Evaluation
gasSolventSat = hydrocarbonSaturation_ + solventSaturation_;
638 if (SolventModule::isMiscible() && FluidSystem::phaseIsActive(oilPhaseIdx)) {
639 const auto& misc = SolventModule::misc(elemCtx, dofIdx,
timeIdx);
640 const auto& pmisc = SolventModule::pmisc(elemCtx, dofIdx,
timeIdx);
641 const Evaluation&
p = FluidSystem::phaseIsActive(oilPhaseIdx)? fs.pressure(oilPhaseIdx) : fs.pressure(gasPhaseIdx);
646 const auto& materialLawManager = elemCtx.problem().materialLawManager();
648 materialLawManager->oilWaterScaledEpsInfoDrainage(
cellIdx);
652 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
653 const Evaluation&
sw = fs.saturation(waterPhaseIdx);
654 const auto& sorwmis = SolventModule::sorwmis(elemCtx, dofIdx,
timeIdx);
658 Evaluation
sgc = sgcr;
659 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
660 const Evaluation&
sw = fs.saturation(waterPhaseIdx);
661 const auto& sgcwmis = SolventModule::sgcwmis(elemCtx, dofIdx,
timeIdx);
666 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
669 const Evaluation
zero = 0.0;
677 const auto& msfnKro = SolventModule::msfnKro(elemCtx, dofIdx,
timeIdx);
678 const auto& msfnKrsg = SolventModule::msfnKrsg(elemCtx, dofIdx,
timeIdx);
679 const auto& sof2Krn = SolventModule::sof2Krn(elemCtx, dofIdx,
timeIdx);
684 Evaluation&
kro = asImp_().mobility_[oilPhaseIdx];
685 Evaluation&
krg = asImp_().mobility_[gasPhaseIdx];
695 const auto& ssfnKrg = SolventModule::ssfnKrg(elemCtx, dofIdx,
timeIdx);
696 const auto& ssfnKrs = SolventModule::ssfnKrs(elemCtx, dofIdx,
timeIdx);
698 Evaluation&
krg = asImp_().mobility_[gasPhaseIdx];
699 solventMobility_ =
krg * ssfnKrs.eval(
Fsolgas,
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;
721 if (SolventModule::isCO2Sol() || SolventModule::isH2Sol() ){
722 if (SolventModule::isCO2Sol()) {
723 const auto&
co2gasPvt = SolventModule::co2GasPvt();
724 solventInvFormationVolumeFactor_ =
co2gasPvt.inverseFormationVolumeFactor(pvtRegionIdx, T,
p, rv, rvw);
725 solventRefDensity_ =
co2gasPvt.gasReferenceDensity(pvtRegionIdx);
726 solventViscosity_ =
co2gasPvt.viscosity(pvtRegionIdx, T,
p, rv, rvw);
728 const auto& brineCo2Pvt = SolventModule::brineCo2Pvt();
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());
741 const auto&
h2gasPvt = SolventModule::h2GasPvt();
742 solventInvFormationVolumeFactor_ =
h2gasPvt.inverseFormationVolumeFactor(pvtRegionIdx, T,
p, rv, rvw);
743 solventRefDensity_ =
h2gasPvt.gasReferenceDensity(pvtRegionIdx);
744 solventViscosity_ =
h2gasPvt.viscosity(pvtRegionIdx, T,
p, rv, rvw);
746 const auto& brineH2Pvt = SolventModule::brineH2Pvt();
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());
760 const auto& solventPvt = SolventModule::solventPvt();
761 solventInvFormationVolumeFactor_ = solventPvt.inverseFormationVolumeFactor(pvtRegionIdx, T,
p);
762 solventRefDensity_ = solventPvt.referenceDensity(pvtRegionIdx);
763 solventViscosity_ = solventPvt.viscosity(pvtRegionIdx, T,
p);
766 solventDensity_ = solventInvFormationVolumeFactor_*solventRefDensity_;
769 solventMobility_ /= solventViscosity_;
774 const Evaluation& solventSaturation()
const
775 {
return solventSaturation_; }
777 const Evaluation& rsSolw()
const
780 const Evaluation& solventDensity()
const
781 {
return solventDensity_; }
783 const Evaluation& solventViscosity()
const
784 {
return solventViscosity_; }
786 const Evaluation& solventMobility()
const
787 {
return solventMobility_; }
789 const Evaluation& solventInverseFormationVolumeFactor()
const
790 {
return solventInvFormationVolumeFactor_; }
793 const Scalar& solventRefDensity()
const
794 {
return solventRefDensity_; }
799 void effectiveProperties(
const ElementContext& elemCtx,
803 if (!SolventModule::isMiscible())
808 if (solventSaturation() < cutOff)
813 auto& fs = asImp_().fluidState_;
817 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
820 Evaluation
gasEffSat = fs.saturation(gasPhaseIdx);
822 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
823 const auto& sorwmis = SolventModule::sorwmis(elemCtx,
scvIdx,
timeIdx);
824 const auto& sgcwmis = SolventModule::sgcwmis(elemCtx,
scvIdx,
timeIdx);
825 const Evaluation
zero = 0.0;
826 const Evaluation&
sw = fs.saturation(waterPhaseIdx);
827 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
842 const Evaluation&
p = FluidSystem::phaseIsActive(oilPhaseIdx)? fs.pressure(oilPhaseIdx) : fs.pressure(gasPhaseIdx);
846 const auto& tlPMixTable = SolventModule::tlPMixTable(elemCtx,
scvIdx,
timeIdx);
849 const Evaluation&
muGas = fs.viscosity(gasPhaseIdx);
850 const Evaluation&
muSolvent = solventViscosity_;
861 Evaluation
muOil = 1.0;
865 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
866 muOil = fs.viscosity(oilPhaseIdx);
884 const Evaluation&
rhoGas = fs.density(gasPhaseIdx);
886 if (FluidSystem::phaseIsActive(oilPhaseIdx))
887 rhoOil = fs.density(oilPhaseIdx);
889 const Evaluation&
rhoSolvent = solventDensity_;
903 Evaluation
sof = 0.0;
904 Evaluation
sgf = 0.0;
933 unsigned pvtRegionIdx = asImp_().pvtRegionIndex();
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);
944 const Evaluation
bs = solventInverseFormationVolumeFactor();
945 const Evaluation
bg_eff = pmisc *
bGasEff + (1.0 - pmisc) * bg;
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));
959 solventDensity_ =
bs_eff*solventRefDensity();
962 Evaluation&
mobg = asImp_().mobility_[gasPhaseIdx];
969 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
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];
994 Implementation& asImp_()
995 {
return *
static_cast<Implementation*
>(
this); }
997 Evaluation hydrocarbonSaturation_;
998 Evaluation solventSaturation_;
1000 Evaluation solventDensity_;
1001 Evaluation solventViscosity_;
1002 Evaluation solventMobility_;
1003 Evaluation solventInvFormationVolumeFactor_;
1005 Scalar solventRefDensity_;
Provides the solvent specific extensive quantities to the generic black-oil module's extensive quanti...
Definition blackoilsolventmodules.hh:1063
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