From 5c2b577854310cb081189a82cfb1bc0619b7902d Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 10 Nov 2023 16:02:17 +0100 Subject: [PATCH] Add H2SOL and CO2SOL model. That allows for dissolution of solvent in water --- .../blackoil/blackoilintensivequantities.hh | 10 +- opm/models/blackoil/blackoilnewtonmethod.hh | 16 +- .../blackoil/blackoilprimaryvariables.hh | 64 ++++- opm/models/blackoil/blackoilsolventmodules.hh | 221 ++++++++++++++++-- opm/models/blackoil/blackoilsolventparams.hh | 21 ++ opm/models/io/vtkblackoilsolventmodule.hh | 23 ++ 6 files changed, 323 insertions(+), 32 deletions(-) diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index 7912b48cc..e6a625a19 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -231,10 +231,12 @@ public: // deal with solvent if constexpr (enableSolvent) { - if (FluidSystem::phaseIsActive(oilPhaseIdx)) { - So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx); - } else if (FluidSystem::phaseIsActive(gasPhaseIdx)) { - Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx); + if(priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) { + if (FluidSystem::phaseIsActive(oilPhaseIdx)) { + So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx); + } else if (FluidSystem::phaseIsActive(gasPhaseIdx)) { + Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx); + } } } diff --git a/opm/models/blackoil/blackoilnewtonmethod.hh b/opm/models/blackoil/blackoilnewtonmethod.hh index 05eed406a..bb5e84a75 100644 --- a/opm/models/blackoil/blackoilnewtonmethod.hh +++ b/opm/models/blackoil/blackoilnewtonmethod.hh @@ -310,7 +310,7 @@ protected: deltaSg = update[Indices::compositionSwitchIdx]; deltaSo -= deltaSg; } - if (enableSolvent) { + if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) { deltaSs = update[Indices::solventSaturationIdx]; deltaSo -= deltaSs; } @@ -365,7 +365,13 @@ protected: } else if (enableSolvent && pvIdx == Indices::solventSaturationIdx) { // solvent saturation updates are also subject to the Appleyard chop - delta *= satAlpha; + if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) { + delta *= satAlpha; + } else { + //Ensure Rssolw factor does not become negative + if (delta > currentValue[Indices::solventSaturationIdx]) + delta = currentValue[Indices::solventSaturationIdx]; + } } else if (enableExtbo && pvIdx == Indices::zFractionIdx) { // z fraction updates are also subject to the Appleyard chop @@ -396,8 +402,10 @@ protected: nextValue[pvIdx] = currentValue[pvIdx] - delta; // keep the solvent saturation between 0 and 1 - if (enableSolvent && pvIdx == Indices::solventSaturationIdx) - nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], 0.0), 1.0); + if (enableSolvent && pvIdx == Indices::solventSaturationIdx) { + if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss ) + nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], 0.0), 1.0); + } // keep the z fraction between 0 and 1 if (enableExtbo && pvIdx == Indices::zFractionIdx) diff --git a/opm/models/blackoil/blackoilprimaryvariables.hh b/opm/models/blackoil/blackoilprimaryvariables.hh index a48ea8522..85e1c988b 100644 --- a/opm/models/blackoil/blackoilprimaryvariables.hh +++ b/opm/models/blackoil/blackoilprimaryvariables.hh @@ -83,6 +83,7 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables enum { pressureSwitchIdx = Indices::pressureSwitchIdx }; enum { compositionSwitchIdx = Indices::compositionSwitchIdx }; enum { saltConcentrationIdx = Indices::saltConcentrationIdx }; + enum { solventSaturationIdx = Indices::solventSaturationIdx }; static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0; static constexpr bool waterEnabled = Indices::waterEnabled; @@ -150,6 +151,12 @@ public: Disabled, // The primary variable is not used }; + enum class SolventMeaning { + Ss, // solvent saturation + Rsolw, // dissolved solvent in water + Disabled, // The primary variable is not used + }; + BlackOilPrimaryVariables() : ParentType() { @@ -167,6 +174,7 @@ public: Valgrind::SetUndefined(primaryVarsMeaningGas_); Valgrind::SetUndefined(primaryVarsMeaningPressure_); Valgrind::SetUndefined(primaryVarsMeaningBrine_); + Valgrind::SetUndefined(primaryVarsMeaningSolvent_); pvtRegionIdx_ = 0; } @@ -184,6 +192,7 @@ public: result.primaryVarsMeaningGas_ = GasMeaning::Rv; result.primaryVarsMeaningPressure_ = PressureMeaning::Pg; result.primaryVarsMeaningWater_ = WaterMeaning::Rsw; + result.primaryVarsMeaningSolvent_ = SolventMeaning::Ss; for (size_t i = 0; i < result.size(); ++i) { result[i] = i+1; } @@ -260,6 +269,18 @@ public: void setPrimaryVarsMeaningBrine(BrineMeaning newMeaning) { primaryVarsMeaningBrine_ = newMeaning; } + + SolventMeaning primaryVarsMeaningSolvent() const + { return primaryVarsMeaningSolvent_; } + + /*! + * \brief Set the interpretation which should be applied to the switching primary + * variables. + */ + + void setPrimaryVarsMeaningSolvent(SolventMeaning newMeaning) + { primaryVarsMeaningSolvent_ = newMeaning; } + /*! * \copydoc ImmisciblePrimaryVariables::assignMassConservative */ @@ -480,7 +501,6 @@ public: default: throw std::logic_error("No valid primary variable selected for composision"); } - checkDefined(); } /*! @@ -542,6 +562,31 @@ public: } } + // if solvent saturation disappeares: Ss (Solvent saturation) -> Rsolw (solvent dissolved in water) + // if solvent saturation appears: Rsolw (solvent dissolved in water) -> Ss (Solvent saturation) + // Scalar rsolw = 0.0; // not needed at the moment since we dont allow for vapwat in combination with rsolw + if constexpr (enableSolvent) { + if (SolventModule::isSolubleInWater()) { + Scalar p = (*this)[pressureSwitchIdx]; // cap-pressure? + Scalar solLimit = SolventModule::solubilityLimit(pvtRegionIndex(), T , p, saltConcentration); + if (primaryVarsMeaningSolvent() == SolventMeaning::Ss) { + Scalar solSat = (*this)[solventSaturationIdx]; + if (solSat < -eps){ //solvent dissappears + setPrimaryVarsMeaningSolvent(SolventMeaning::Rsolw); + (*this)[solventSaturationIdx] = solLimit; //set rsolw to solubility limit + + } + } + else if (primaryVarsMeaningSolvent() == SolventMeaning::Rsolw) { + Scalar rsolw = (*this)[solventSaturationIdx]; + if (rsolw > solLimit + eps){ //solvent appears as phase + setPrimaryVarsMeaningSolvent(SolventMeaning::Ss); + (*this)[solventSaturationIdx] = 0.0; + } + } + } + } + // keep track if any meaning has changed bool changed = false; @@ -826,7 +871,7 @@ public: sg = (*this)[Indices::compositionSwitchIdx]; Scalar ssol = 0.0; - if constexpr (enableSolvent) + if (primaryVarsMeaningSolvent() == SolventMeaning::Ss) ssol =(*this) [Indices::solventSaturationIdx]; Scalar so = 1.0 - sw - sg - ssol; @@ -843,7 +888,7 @@ public: (*this)[Indices::waterSwitchIdx] = sw; if (primaryVarsMeaningGas() == GasMeaning::Sg) (*this)[Indices::compositionSwitchIdx] = sg; - if constexpr (enableSolvent) + if (primaryVarsMeaningSolvent() == SolventMeaning::Ss) (*this) [Indices::solventSaturationIdx] = ssol; return !(st==1); @@ -877,6 +922,7 @@ public: Valgrind::CheckDefined(primaryVarsMeaningGas_); Valgrind::CheckDefined(primaryVarsMeaningPressure_); Valgrind::CheckDefined(primaryVarsMeaningBrine_); + Valgrind::CheckDefined(primaryVarsMeaningSolvent_); Valgrind::CheckDefined(pvtRegionIdx_); #endif // NDEBUG @@ -891,6 +937,7 @@ public: serializer(primaryVarsMeaningPressure_); serializer(primaryVarsMeaningGas_); serializer(primaryVarsMeaningBrine_); + serializer(primaryVarsMeaningSolvent_); serializer(pvtRegionIdx_); } @@ -901,6 +948,7 @@ public: this->primaryVarsMeaningPressure_ == rhs.primaryVarsMeaningPressure_ && this->primaryVarsMeaningGas_ == rhs.primaryVarsMeaningGas_ && this->primaryVarsMeaningBrine_ == rhs.primaryVarsMeaningBrine_ && + this->primaryVarsMeaningSolvent_ == rhs.primaryVarsMeaningSolvent_ && this->pvtRegionIdx_ == rhs.pvtRegionIdx_; } @@ -913,10 +961,11 @@ private: Scalar solventSaturation_() const { - if constexpr (enableSolvent) - return (*this)[Indices::solventSaturationIdx]; - else - return 0.0; + if constexpr (enableSolvent) { + if ( primaryVarsMeaningSolvent() == SolventMeaning::Ss) + return (*this)[Indices::solventSaturationIdx]; + } + return 0.0; } Scalar zFraction_() const @@ -1033,6 +1082,7 @@ private: PressureMeaning primaryVarsMeaningPressure_; GasMeaning primaryVarsMeaningGas_; BrineMeaning primaryVarsMeaningBrine_; + SolventMeaning primaryVarsMeaningSolvent_; unsigned short pvtRegionIdx_; }; diff --git a/opm/models/blackoil/blackoilsolventmodules.hh b/opm/models/blackoil/blackoilsolventmodules.hh index 701c2e6d5..787d73190 100644 --- a/opm/models/blackoil/blackoilsolventmodules.hh +++ b/opm/models/blackoil/blackoilsolventmodules.hh @@ -78,9 +78,12 @@ class BlackOilSolventModule using EqVector = GetPropType; using RateVector = GetPropType; using Indices = GetPropType; - using Toolbox = MathToolbox; using SolventPvt = typename BlackOilSolventParams::SolventPvt; + using Co2GasPvt = typename BlackOilSolventParams::Co2GasPvt; + using H2GasPvt = typename BlackOilSolventParams::H2GasPvt; + using BrineCo2Pvt = typename BlackOilSolventParams::BrineCo2Pvt; + using BrineH2Pvt = typename BlackOilSolventParams::BrineH2Pvt; using TabulatedFunction = typename BlackOilSolventParams::TabulatedFunction; @@ -90,6 +93,7 @@ class BlackOilSolventModule static constexpr unsigned numEq = getPropValue(); static constexpr unsigned numPhases = FluidSystem::numPhases; static constexpr bool blackoilConserveSurfaceVolume = getPropValue(); + static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx; public: @@ -111,7 +115,26 @@ public: if (!eclState.runspec().phases().active(Phase::SOLVENT)) return; // solvent treatment is supposed to be disabled - params_.solventPvt_.initFromState(eclState, schedule); + params_.co2sol_ = eclState.runspec().co2Sol(); + params_.h2sol_ = eclState.runspec().h2Sol(); + + if (isCO2Sol() && isH2Sol()) { + throw std::runtime_error("CO2SOL and H2SOL can not be used together"); + } + + if (isCO2Sol() || isH2Sol() ) { + if (isCO2Sol()) { + params_.co2GasPvt_.initFromState(eclState, schedule); + params_.brineCo2Pvt_.initFromState(eclState, schedule); + } else { + params_.h2GasPvt_.initFromState(eclState, schedule); + params_.brineH2Pvt_.initFromState(eclState, schedule); + } + if (eclState.getSimulationConfig().hasDISGASW()) { + params_.rsSolw_active_ = true; + } + } else + params_.solventPvt_.initFromState(eclState, schedule); const auto& tableManager = eclState.getTableManager(); // initialize the objects which deal with the SSFN keyword @@ -413,12 +436,25 @@ public: Toolbox::template decay(intQuants.porosity()) * Toolbox::template decay(intQuants.solventSaturation()) * Toolbox::template decay(intQuants.solventInverseFormationVolumeFactor()); + if (isSolubleInWater()) { + storage[contiSolventEqIdx] += Toolbox::template decay(intQuants.porosity()) + * Toolbox::template decay(intQuants.fluidState().saturation(waterPhaseIdx)) + * Toolbox::template decay(intQuants.fluidState().invB(waterPhaseIdx)) + * Toolbox::template decay(intQuants.rsSolw()); + } } else { storage[contiSolventEqIdx] += Toolbox::template decay(intQuants.porosity()) * Toolbox::template decay(intQuants.solventSaturation()) * Toolbox::template decay(intQuants.solventDensity()); + if (isSolubleInWater()) { + storage[contiSolventEqIdx] += Toolbox::template decay(intQuants.porosity()) + * Toolbox::template decay(intQuants.fluidState().saturation(waterPhaseIdx)) + * Toolbox::template decay(intQuants.fluidState().density(waterPhaseIdx)) + * Toolbox::template decay(intQuants.rsSolw()); + } + } } } @@ -437,7 +473,7 @@ public: const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx); if constexpr (blackoilConserveSurfaceVolume) { - if (upIdx == inIdx) + if (upIdx == inIdx) flux[contiSolventEqIdx] = extQuants.solventVolumeFlux() *up.solventInverseFormationVolumeFactor(); @@ -445,6 +481,20 @@ public: flux[contiSolventEqIdx] = extQuants.solventVolumeFlux() *decay(up.solventInverseFormationVolumeFactor()); + + + if (isSolubleInWater()) { + if (upIdx == inIdx) + flux[contiSolventEqIdx] = + extQuants.volumeFlux(waterPhaseIdx) + * up.fluidState().invB(waterPhaseIdx) + * up.rsSolw(); + else + flux[contiSolventEqIdx] = + extQuants.volumeFlux(waterPhaseIdx) + *decay(up.fluidState().invB(waterPhaseIdx)) + *decay(up.rsSolw()); + } } else { if (upIdx == inIdx) @@ -455,6 +505,20 @@ public: flux[contiSolventEqIdx] = extQuants.solventVolumeFlux() *decay(up.solventDensity()); + + + if (isSolubleInWater()) { + if (upIdx == inIdx) + flux[contiSolventEqIdx] = + extQuants.volumeFlux(waterPhaseIdx) + * up.fluidState().density(waterPhaseIdx) + * up.rsSolw(); + else + flux[contiSolventEqIdx] = + extQuants.volumeFlux(waterPhaseIdx) + *decay(up.fluidState().density(waterPhaseIdx)) + *decay(up.rsSolw()); + } } } } @@ -463,10 +527,21 @@ public: * \brief Assign the solvent specific primary variables to a PrimaryVariables object */ static void assignPrimaryVars(PrimaryVariables& priVars, - Scalar solventSaturation) + Scalar solventSaturation, + Scalar solventRsw) { - if constexpr (enableSolvent) + if constexpr (!enableSolvent) { + priVars.setPrimaryVarsMeaningSolvent(PrimaryVariables::SolventMeaning::Disabled); + return; + } + // Determine the meaning of the solvent primary variables + if (solventSaturation > 0 || !isSolubleInWater()) { + priVars.setPrimaryVarsMeaningSolvent(PrimaryVariables::SolventMeaning::Ss); priVars[solventSaturationIdx] = solventSaturation; + } else { + priVars.setPrimaryVarsMeaningSolvent(PrimaryVariables::SolventMeaning::Rsolw); + priVars[solventSaturationIdx] = solventRsw; + } } /*! @@ -530,7 +605,30 @@ public: } static const SolventPvt& solventPvt() - { return params_.solventPvt_; } + { + return params_.solventPvt_; + } + + + static const Co2GasPvt& co2GasPvt() + { + return params_.co2GasPvt_; + } + + static const H2GasPvt& h2GasPvt() + { + return params_.h2GasPvt_; + } + + static const BrineCo2Pvt& brineCo2Pvt() + { + return params_.brineCo2Pvt_; + } + + static const BrineH2Pvt& brineH2Pvt() + { + return params_.brineH2Pvt_; + } static const TabulatedFunction& ssfnKrg(const ElementContext& elemCtx, unsigned scvIdx, @@ -633,6 +731,34 @@ public: return params_.isMiscible_; } + template + static const Value solubilityLimit(unsigned pvtIdx, const Value& temperature, const Value& pressure, const Value& saltConcentration) + { + if (!isSolubleInWater()) + return 0.0; + + assert(isCO2Sol() || isH2Sol()); + if (isCO2Sol()) + return brineCo2Pvt().saturatedGasDissolutionFactor(pvtIdx, temperature, pressure, saltConcentration); + else + return brineH2Pvt().saturatedGasDissolutionFactor(pvtIdx, temperature, pressure, saltConcentration); + } + + static bool isSolubleInWater() + { + return params_.rsSolw_active_; + } + + static bool isCO2Sol() + { + return params_.co2sol_; + } + + static bool isH2Sol() + { + return params_.h2sol_; + } + private: static BlackOilSolventParams params_; // the krg(Fs) column of the SSFN table }; @@ -683,8 +809,13 @@ public: unsigned timeIdx) { const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx); + auto& fs = asImp_().fluidState_; - solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx, elemCtx.linearizationType()); + solventSaturation_ = 0.0; + if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) { + solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx, elemCtx.linearizationType()); + } + hydrocarbonSaturation_ = fs.saturation(gasPhaseIdx); // apply a cut-off. Don't waste calculations if no solvent @@ -712,6 +843,16 @@ public: auto& fs = asImp_().fluidState_; fs.setSaturation(gasPhaseIdx, hydrocarbonSaturation_); + + // update rsSolw. This needs to be done after the pressure is defined in the fluid state. + rsSolw_ = 0.0; + const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx); + if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) { + rsSolw_ = SolventModule::solubilityLimit(asImp_().pvtRegionIndex(), fs.temperature(waterPhaseIdx), fs.pressure(waterPhaseIdx), fs.saltConcentration()); + } else if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Rsolw) { + rsSolw_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx, elemCtx.linearizationType()); + } + solventMobility_ = 0.0; // apply a cut-off. Don't waste calculations if no solvent @@ -726,7 +867,6 @@ public: // compute capillary pressure for miscible fluid const auto& problem = elemCtx.problem(); - const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx); Evaluation pgMisc = 0.0; std::array pC; const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, timeIdx); @@ -824,18 +964,58 @@ public: unsigned timeIdx) { const auto& iq = asImp_(); - const auto& fs = iq.fluidState(); - const auto& solventPvt = SolventModule::solventPvt(); - unsigned pvtRegionIdx = iq.pvtRegionIndex(); - solventRefDensity_ = solventPvt.referenceDensity(pvtRegionIdx); - const Evaluation& T = fs.temperature(gasPhaseIdx); - const Evaluation& p = fs.pressure(gasPhaseIdx); - solventInvFormationVolumeFactor_ = solventPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p); + const Evaluation& T = iq.fluidState().temperature(gasPhaseIdx); + const Evaluation& p = iq.fluidState().pressure(gasPhaseIdx); + + const Evaluation rv = 0.0; + const Evaluation rvw = 0.0; + if (SolventModule::isCO2Sol() || SolventModule::isH2Sol() ){ + if (SolventModule::isCO2Sol()) { + const auto& co2gasPvt = SolventModule::co2GasPvt(); + solventInvFormationVolumeFactor_ = co2gasPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rv, rvw); + solventRefDensity_ = co2gasPvt.gasReferenceDensity(pvtRegionIdx); + solventViscosity_ = co2gasPvt.viscosity(pvtRegionIdx, T, p, rv, rvw); + + const auto& brineCo2Pvt = SolventModule::brineCo2Pvt(); + auto& fs = asImp_().fluidState_; + const auto& bw = brineCo2Pvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rsSolw()); + + const auto denw = bw*brineCo2Pvt.waterReferenceDensity(pvtRegionIdx) + + rsSolw()*bw*brineCo2Pvt.gasReferenceDensity(pvtRegionIdx); + fs.setDensity(waterPhaseIdx, denw); + fs.setInvB(waterPhaseIdx, bw); + Evaluation& mobw = asImp_().mobility_[waterPhaseIdx]; + const auto& muWat = fs.viscosity(waterPhaseIdx); + const auto& muWatEff = brineCo2Pvt.viscosity(pvtRegionIdx, T, p, rsSolw()); + mobw *= muWat / muWatEff; + } else { + const auto& h2gasPvt = SolventModule::h2GasPvt(); + solventInvFormationVolumeFactor_ = h2gasPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rv, rvw); + solventRefDensity_ = h2gasPvt.gasReferenceDensity(pvtRegionIdx); + solventViscosity_ = h2gasPvt.viscosity(pvtRegionIdx, T, p, rv, rvw); + + const auto& brineH2Pvt = SolventModule::brineH2Pvt(); + auto& fs = asImp_().fluidState_; + const auto& bw = brineH2Pvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rsSolw()); + + const auto denw = bw*brineH2Pvt.waterReferenceDensity(pvtRegionIdx) + + rsSolw()*bw*brineH2Pvt.gasReferenceDensity(pvtRegionIdx); + fs.setDensity(waterPhaseIdx, denw); + fs.setInvB(waterPhaseIdx, bw); + Evaluation& mobw = asImp_().mobility_[waterPhaseIdx]; + const auto& muWat = fs.viscosity(waterPhaseIdx); + const auto& muWatEff = brineH2Pvt.viscosity(pvtRegionIdx, T, p, rsSolw()); + mobw *= muWat / muWatEff; + } + } else { + const auto& solventPvt = SolventModule::solventPvt(); + solventInvFormationVolumeFactor_ = solventPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p); + solventRefDensity_ = solventPvt.referenceDensity(pvtRegionIdx); + solventViscosity_ = solventPvt.viscosity(pvtRegionIdx, T, p); + } solventDensity_ = solventInvFormationVolumeFactor_*solventRefDensity_; - solventViscosity_ = solventPvt.viscosity(pvtRegionIdx, T, p); - effectiveProperties(elemCtx, scvIdx, timeIdx); solventMobility_ /= solventViscosity_; @@ -846,6 +1026,9 @@ public: const Evaluation& solventSaturation() const { return solventSaturation_; } + const Evaluation& rsSolw() const + { return rsSolw_; } + const Evaluation& solventDensity() const { return solventDensity_; } @@ -1036,6 +1219,7 @@ protected: Evaluation hydrocarbonSaturation_; Evaluation solventSaturation_; + Evaluation rsSolw_; Evaluation solventDensity_; Evaluation solventViscosity_; Evaluation solventMobility_; @@ -1071,6 +1255,9 @@ public: const Evaluation& solventSaturation() const { throw std::runtime_error("solventSaturation() called but solvents are disabled"); } + const Evaluation& rsSolw() const + { throw std::runtime_error("rsSolw() called but solvents are disabled"); } + const Evaluation& solventDensity() const { throw std::runtime_error("solventDensity() called but solvents are disabled"); } diff --git a/opm/models/blackoil/blackoilsolventparams.hh b/opm/models/blackoil/blackoilsolventparams.hh index 4dc8b6ab6..c96c91d38 100644 --- a/opm/models/blackoil/blackoilsolventparams.hh +++ b/opm/models/blackoil/blackoilsolventparams.hh @@ -29,6 +29,11 @@ #define EWOMS_BLACK_OIL_SOLVENT_PARAMS_HH #include +#include +#include +#include +#include + #include namespace Opm { @@ -40,6 +45,19 @@ struct BlackOilSolventParams { using SolventPvt = ::Opm::SolventPvt; SolventPvt solventPvt_; + + using Co2GasPvt = ::Opm::Co2GasPvt; + Co2GasPvt co2GasPvt_; + + using H2GasPvt = ::Opm::H2GasPvt; + H2GasPvt h2GasPvt_; + + using BrineCo2Pvt = ::Opm::BrineCo2Pvt; + BrineCo2Pvt brineCo2Pvt_; + + using BrineH2Pvt = ::Opm::BrineH2Pvt; + BrineH2Pvt brineH2Pvt_; + std::vector ssfnKrg_; // the krg(Fs) column of the SSFN table std::vector ssfnKrs_; // the krs(Fs) column of the SSFN table std::vector sof2Krn_; // the krn(Sn) column of the SOF2 table @@ -55,6 +73,9 @@ struct BlackOilSolventParams { std::vector tlPMixTable_; // the tlpmixpa(Po) column of the TLPMIXPA table bool isMiscible_; + bool rsSolw_active_ = false; + bool co2sol_; + bool h2sol_; /*! * \brief Specify the number of satuation regions. diff --git a/opm/models/io/vtkblackoilsolventmodule.hh b/opm/models/io/vtkblackoilsolventmodule.hh index 6ef0311b8..3b6daeacb 100644 --- a/opm/models/io/vtkblackoilsolventmodule.hh +++ b/opm/models/io/vtkblackoilsolventmodule.hh @@ -53,6 +53,8 @@ struct VtkBlackOilSolvent {}; template struct VtkWriteSolventSaturation { using type = UndefinedProperty; }; template +struct VtkWriteSolventRsw { using type = UndefinedProperty; }; +template struct VtkWriteSolventDensity { using type = UndefinedProperty; }; template struct VtkWriteSolventViscosity { using type = UndefinedProperty; }; @@ -63,6 +65,8 @@ struct VtkWriteSolventMobility { using type = UndefinedProperty; }; template struct VtkWriteSolventSaturation { static constexpr bool value = true; }; template +struct VtkWriteSolventRsw { static constexpr bool value = true; }; +template struct VtkWriteSolventDensity { static constexpr bool value = true; }; template struct VtkWriteSolventViscosity { static constexpr bool value = true; }; @@ -112,6 +116,9 @@ public: EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSolventSaturation, "Include the \"saturation\" of the solvent component " "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSolventRsw, + "Include the \"dissolved volume in water\" of the solvent component " + "in the VTK output files"); EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSolventDensity, "Include the \"density\" of the solvent component " "in the VTK output files"); @@ -137,6 +144,8 @@ public: if (solventSaturationOutput_()) this->resizeScalarBuffer_(solventSaturation_); + if (solventRswOutput_()) + this->resizeScalarBuffer_(solventRsw_); if (solventDensityOutput_()) this->resizeScalarBuffer_(solventDensity_); if (solventViscosityOutput_()) @@ -166,6 +175,10 @@ public: solventSaturation_[globalDofIdx] = Toolbox::scalarValue(intQuants.solventSaturation()); + if (solventRswOutput_()) + solventRsw_[globalDofIdx] = + Toolbox::scalarValue(intQuants.rsSolw()); + if (solventDensityOutput_()) solventDensity_[globalDofIdx] = Toolbox::scalarValue(intQuants.solventDensity()); @@ -195,6 +208,9 @@ public: if (solventSaturationOutput_()) this->commitScalarBuffer_(baseWriter, "saturation_solvent", solventSaturation_); + if (solventRswOutput_()) + this->commitScalarBuffer_(baseWriter, "dissolved_solvent", solventRsw_); + if (solventDensityOutput_()) this->commitScalarBuffer_(baseWriter, "density_solvent", solventDensity_); @@ -212,6 +228,12 @@ private: return val; } + static bool solventRswOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteSolventRsw); + return val; + } + static bool solventDensityOutput_() { static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteSolventDensity); @@ -231,6 +253,7 @@ private: } ScalarBuffer solventSaturation_; + ScalarBuffer solventRsw_; ScalarBuffer solventDensity_; ScalarBuffer solventViscosity_; ScalarBuffer solventMobility_;