OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

extending drsdtcon with regimes, and option for GAS/WATER

Browse Source
This commit is contained in:
Trine Mykkeltvedt
2024-06-07 08:58:27 +02:00
committed by Tor Harald Sandve
parent ad595fed5e
commit ab2bd924db
4 changed files with 104 additions and 25 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
opm/simulators/flow/EquilInitializer.hpp
View File

@@ -174,6 +174,7 @@ public:
const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, regionIdx);
fluidState.setEnthalpy(phaseIdx, h);
}
}
// set the salt concentration

18
opm/simulators/flow/FlowProblem.hpp
View File

@@ -57,7 +57,9 @@
#include <opm/models/common/directionalmobility.hh>
#include <opm/models/utils/pffgridvector.hh>
#include <opm/models/blackoil/blackoilmodel.hh>
#include <opm/models/blackoil/blackoilconvectivemixingmodule.hh>
#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
#include <opm/output/eclipse/EclipseIO.hpp>
@@ -180,6 +182,8 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
using MICPModule = BlackOilMICPModule<TypeTag>;
using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag>;
using ModuleParams = typename BlackOilLocalResidualTPFA<TypeTag>::ModuleParams;
using InitialFluidState = typename EquilInitializer<TypeTag>::ScalarFluidState;
@@ -534,6 +538,8 @@ public:
const auto& schedule = simulator.vanguard().schedule();
const auto& events = schedule[episodeIdx].events();
if (episodeIdx >= 0 && events.hasEvent(ScheduleEvents::GEO_MODIFIER)) {
// bring the contents of the keywords to the current state of the SCHEDULE
// section.
@@ -589,6 +595,8 @@ public:
FluidSystem::setVapPars(0.0, 0.0);
}
}
ConvectiveMixingModule::beginEpisode(simulator.vanguard().eclState(), simulator.vanguard().schedule(), episodeIdx, moduleParams_.convectiveMixingModuleParam);
}
/*!
@@ -1855,6 +1863,10 @@ public:
eclWriter_->mutableOutputModule().setCnvData(data);
}
const ModuleParams& moduleParams() const {
return moduleParams_;
}
template<class Serializer>
void serializeOp(Serializer& serializer)
{
@@ -1953,6 +1965,7 @@ protected:
return true;
}
bool updateMaxOilSaturation_()
{
OPM_TIMEBLOCK(updateMaxOilSaturation);
@@ -2597,7 +2610,6 @@ protected:
return this->rockCompTransMultVal_[dofIdx];
}
private:
struct PffDofData_
{
@@ -2883,6 +2895,10 @@ private:
BCData<int> bcindex_;
bool nonTrivialBoundaryConditions_ = false;
bool explicitRockCompaction_ = false;
ModuleParams moduleParams_;
};
} // namespace Opm

66
opm/simulators/flow/MixingRateControls.cpp
View File

@@ -101,8 +101,7 @@ init(std::size_t numDof, int episodeIdx, const unsigned ntpvt)
//TODO We may want to only allocate these properties only if active.
//But since they may be activated at later time we need some more
//intrastructure to handle it
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
maxDRv_.resize(ntpvt, 1e30);
lastRv_.resize(numDof, 0.0);
maxDRs_.resize(ntpvt, 1e30);
@@ -113,7 +112,6 @@ init(std::size_t numDof, int episodeIdx, const unsigned ntpvt)
if (this->drsdtConvective(episodeIdx)) {
convectiveDrs_.resize(numDof, 1.0);
}
}
}
template<class FluidSystem>
@@ -123,7 +121,7 @@ drsdtActive(int episodeIdx) const
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
return (oilVaporizationControl.drsdtActive() && bothOilGasActive);
return (oilVaporizationControl.drsdtActive());
}
template<class FluidSystem>
@@ -133,7 +131,7 @@ drvdtActive(int episodeIdx) const
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
return (oilVaporizationControl.drvdtActive() && bothOilGasActive);
return (oilVaporizationControl.drvdtActive());
}
template<class FluidSystem>
@@ -143,7 +141,7 @@ drsdtConvective(int episodeIdx) const
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
return (oilVaporizationControl.drsdtConvective() && bothOilGasActive);
return (oilVaporizationControl.drsdtConvective());
}
template<class FluidSystem>
@@ -273,29 +271,61 @@ updateConvectiveDRsDt_(const unsigned compressedDofIdx,
const Scalar p,
const Scalar rs,
const Scalar so,
const Scalar sg,
const Scalar poro,
const Scalar permz,
const Scalar distZ,
const Scalar gravity,
const Scalar salt,
const Scalar Xhi,
const Scalar Psi,
const Scalar omegainn,
const int pvtRegionIndex)
{
const Scalar rssat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIndex, t, p);
const Scalar saturatedInvB
= FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIndex, t, p);
const Scalar rssat = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
FluidSystem::waterPvt().saturatedGasDissolutionFactor(pvtRegionIndex, t, p, salt) :
FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIndex, t, p);
const Scalar saturatedInvB = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
FluidSystem::waterPvt().saturatedInverseFormationVolumeFactor(pvtRegionIndex, t, p, salt) :
FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIndex, t, p);
const Scalar rsZero = 0.0;
const Scalar pureDensity
= FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIndex, t, p, rsZero)
* FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex);
const Scalar saturatedDensity = saturatedInvB
* (FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex)
+ rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIndex));
const Scalar deltaDensity = saturatedDensity - pureDensity;
const Scalar visc = FluidSystem::oilPvt().viscosity(pvtRegionIndex, t, p, rs);
const Scalar sg_max = 1.0;
const Scalar pureDensity = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
(FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIndex, t, p, rsZero, salt)
* FluidSystem::waterPvt().waterReferenceDensity(pvtRegionIndex)) :
(FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIndex, t, p, rsZero)
* FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex));
const Scalar saturatedDensity = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
(saturatedInvB * (FluidSystem::waterPvt().waterReferenceDensity(pvtRegionIndex)
+ rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIndex))) :
(saturatedInvB * (FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIndex)
+ rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIndex)));
Scalar deltaDensity = saturatedDensity - pureDensity;
const Scalar visc = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
FluidSystem::waterPvt().viscosity(pvtRegionIndex, t, p, rs, salt) :
FluidSystem::oilPvt().viscosity(pvtRegionIndex, t, p, rs);
// Note that for so = 0 this gives no limits (inf) for the dissolution rate
// Also we restrict the effect of convective mixing to positive density differences
// i.e. we only allow for fingers moving downward
Scalar co2Density = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIndex,
t,p,0.0 /*=Rv*/, 0.0 /*=Rvw*/) * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIndex);
Scalar factor = 1.0;
Scalar X = (rs - rssat * sg) / (rssat * ( 1.0 - sg));
Scalar omega = 0.0;
if ((rs >= (rssat * sg))){
if(X > Psi){
factor = 0.0;
omega = omegainn;
}
} else {
factor /= Xhi;
deltaDensity = (saturatedDensity - co2Density);
}
convectiveDrs_[compressedDofIdx]
= permz * rssat * max(0.0, deltaDensity) * gravity / (so * visc * distZ * poro);
= factor * permz * rssat * max(0.0, deltaDensity) * gravity / ( std::max(sg_max - sg, 0.0) * visc * distZ * poro) + (omega/Xhi);
}
template class MixingRateControls<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>>;

44
opm/simulators/flow/MixingRateControls.hpp
View File

@@ -117,20 +117,46 @@ public:
const int pvtRegionIdx,
const std::array<bool,3>& active)
{
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
if (active[0]) {
// This implements the convective DRSDT as described in
// Sandve et al. "Convective dissolution in field scale CO2 storage simulations using the OPM Flow
// simulator" Submitted to TCCS 11, 2021
// modification and introduction of regimes following Mykkeltvedt et al. Submitted to TIMP 2024
const auto& fs = iq.fluidState();
const auto& temperature = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
getValue(fs.temperature(FluidSystem::waterPhaseIdx)) :
getValue(fs.temperature(FluidSystem::oilPhaseIdx));
const auto& pressure = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
getValue(fs.pressure(FluidSystem::waterPhaseIdx)) :
getValue(fs.pressure(FluidSystem::oilPhaseIdx));
const auto& sLiquid = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
getValue(fs.saturation(FluidSystem::waterPhaseIdx)) :
getValue(fs.saturation(FluidSystem::oilPhaseIdx));
const auto& rs = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
getValue(fs.Rsw()) :
getValue(fs.Rs());
const auto& salt = getValue(fs.saltSaturation());
this->updateConvectiveDRsDt_(compressedDofIdx,
getValue(fs.temperature(FluidSystem::oilPhaseIdx)),
getValue(fs.pressure(FluidSystem::oilPhaseIdx)),
getValue(fs.Rs()),
getValue(fs.saturation(FluidSystem::oilPhaseIdx)),
temperature,
pressure,
rs,
sLiquid,
getValue(fs.saturation(FluidSystem::gasPhaseIdx)),
getValue(iq.porosity()),
permZ,
distZ,
gravity,
salt,
oilVaporizationControl.getMaxDRSDT(fs.pvtRegionIndex()),
oilVaporizationControl.getPsi(fs.pvtRegionIndex()),
oilVaporizationControl.getOmega(fs.pvtRegionIndex()),
fs.pvtRegionIndex());
}
@@ -138,13 +164,14 @@ public:
const auto& fs = iq.fluidState();
using FluidState = typename std::decay<decltype(fs)>::type;
const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
constexpr Scalar freeGasMinSaturation_ = 1e-7;
if (oilVaporizationControl.getOption(pvtRegionIdx) ||
fs.saturation(FluidSystem::gasPhaseIdx) > freeGasMinSaturation_) {
lastRs_[compressedDofIdx]
= BlackOil::template getRs_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex());
= (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
BlackOil::template getRsw_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex()) :
BlackOil::template getRs_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex());
}
else
lastRs_[compressedDofIdx] = std::numeric_limits<Scalar>::infinity();
@@ -164,10 +191,15 @@ private:
const Scalar p,
const Scalar rs,
const Scalar so,
const Scalar sg,
const Scalar poro,
const Scalar permz,
const Scalar distZ,
const Scalar gravity,
const Scalar salt,
const Scalar Xhi,
const Scalar Psi,
const Scalar omegainn,
const int pvtRegionIndex);
std::vector<Scalar> lastRv_;