OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-25 18:50:19 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #3076 from totto82/drsdt_dyn

Browse Source 
Implement convective dissolution rate
This commit is contained in:
Tor Harald Sandve 2021-03-18 08:45:46 +01:00 committed by GitHub
commit 3df75f5ab6
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 141 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
ebos/eclbasevanguard.hh
View File

@ -735,7 +735,7 @@ public:
}
/*!
* \brief Returns the depth of an degree of freedom [m]
* \brief Returns the depth of a degree of freedom [m]
*
* For ECL problems this is defined as the average of the depth of an element and is
* thus slightly different from the depth of an element's centroid.
@ -745,6 +745,19 @@ public:
return cellCenterDepth_[globalSpaceIdx];
}
/*!
* \brief Returns the thickness of a degree of freedom [m]
*
* For ECL problems this is defined as the average of the depths of the top surface
* corners minus the average of the depths of the bottom surface corners
* The cell thickness is computed only when needed.
*/
Scalar cellThickness(unsigned globalSpaceIdx) const
{
assert(!cellThickness_.empty());
return cellThickness_[globalSpaceIdx];
}
/*!
* \brief Get the number of cells in the global leaf grid view.
* \warn This is a collective operation that needs to be called
@ -809,6 +822,38 @@ protected:
}
}
}
void updateCellThickness_()
{
bool drsdtcon = false;
auto schIt = this->schedule().begin();
const auto& schEndIt = this->schedule().end();
for(; schIt != schEndIt; ++schIt) {
const auto& oilVaporizationControl = schIt->oilvap();
if(oilVaporizationControl.getType() == Opm::OilVaporizationProperties::OilVaporization::DRSDTCON) {
drsdtcon = true;
break;
}
}
if (!drsdtcon)
return;
ElementMapper elemMapper(this->gridView(), Dune::mcmgElementLayout());
int numElements = this->gridView().size(/*codim=*/0);
cellThickness_.resize(numElements);
auto elemIt = this->gridView().template begin</*codim=*/0>();
const auto& elemEndIt = this->gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& element = *elemIt;
const unsigned int elemIdx = elemMapper.index(element);
cellThickness_[elemIdx] = asImp_().computeCellThickness(element);
}
}
Scalar computeCellThickness(const Element& element) const {
OPM_THROW(std::runtime_error, "cellThickness not implemented for this grid!");
}
private:
void updateOutputDir_()
@ -839,6 +884,8 @@ private:
ioConfig.setEclCompatibleRST(!EWOMS_GET_PARAM(TypeTag, bool, EnableOpmRstFile));
}
// computed from averaging cell corner depths
Scalar cellCenterDepth(const Element& element) const
{
typedef typename Element::Geometry Geometry;
@ -902,11 +949,14 @@ protected:
*/
std::vector<int> cartesianToCompressed_;
/*! \brief Cell center depths computed
* from averaging cell corner depths
/*! \brief Cell center depths
*/
std::vector<Scalar> cellCenterDepth_;
/*! \brief Cell thichness
*/
std::vector<Scalar> cellThickness_;
/*! \brief information about wells in parallel
*
* For each well in the model there is an entry with its name

27
ebos/eclcpgridvanguard.hh
View File

@ -99,6 +99,7 @@ public:
private:
typedef Dune::CartesianIndexMapper<Grid> CartesianIndexMapper;
using Element = typename GridView::template Codim<0>::Entity;
public:
EclCpGridVanguard(Simulator& simulator)
@ -269,7 +270,7 @@ public:
this->updateGridView_();
this->updateCartesianToCompressedMapping_();
this->updateCellDepths_();
this->updateCellThickness_();
#if HAVE_MPI
if (mpiSize > 1) {
@ -406,6 +407,30 @@ protected:
#endif
}
Scalar computeCellThickness(const Element& element) const
{
typedef typename Element::Geometry Geometry;
static constexpr int zCoord = Element::dimension - 1;
Scalar zz1 = 0.0;
Scalar zz2 = 0.0;
const Geometry geometry = element.geometry();
const int corners = geometry.corners();
// This code only works with CP-grid where the
// number of corners are 8 and
// also assumes that the first
// 4 corners are the top surface and
// the 4 next are the bottomn.
assert(corners == 8);
for (int i=0; i < 4; ++i){
zz1 += geometry.corner(i)[zCoord];
zz2 += geometry.corner(i+4)[zCoord];
}
zz1 /=4;
zz2 /=4;
return zz2-zz1;
}
std::unique_ptr<Grid> grid_;
std::unique_ptr<EquilGrid> equilGrid_;
std::unique_ptr<CartesianIndexMapper> cartesianIndexMapper_;

64
ebos/eclproblem.hh
View File

@ -896,13 +896,21 @@ public:
// deal with DRSDT
unsigned ntpvt = eclState.runspec().tabdims().getNumPVTTables();
size_t numDof = this->model().numGridDof();
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
//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);
dRsDtOnlyFreeGas_.resize(ntpvt, false);
lastRs_.resize(numDof, 0.0);
maxDRv_.resize(ntpvt, 1e30);
lastRv_.resize(numDof, 0.0);
maxOilSaturation_.resize(numDof, 0.0);
if (drsdtConvective_()) {
convectiveDrs_.resize(numDof, 1.0);
}
}
readRockParameters_();
@ -1974,10 +1982,15 @@ public:
if (!drsdtActive_() || maxDRs_[pvtRegionIdx] < 0.0)
return std::numeric_limits<Scalar>::max()/2.0;
Scalar scaling = 1.0;
if(drsdtConvective_()) {
scaling = convectiveDrs_[globalDofIdx];
}
// this is a bit hacky because it assumes that a time discretization with only
// two time indices is used.
if (timeIdx == 0)
return lastRs_[globalDofIdx] + maxDRs_[pvtRegionIdx];
return lastRs_[globalDofIdx] + maxDRs_[pvtRegionIdx] * scaling;
else
return lastRs_[globalDofIdx];
}
@ -2317,6 +2330,15 @@ private:
}
bool drsdtConvective_() const
{
const auto& simulator = this->simulator();
int episodeIdx = std::max(simulator.episodeIndex(), 0);
const auto& oilVaporizationControl = simulator.vanguard().schedule()[episodeIdx].oilvap();
return (oilVaporizationControl.drsdtConvective());
}
// update the parameters needed for DRSDT and DRVDT
void updateCompositionChangeLimits_()
{
@ -2326,6 +2348,43 @@ private:
int episodeIdx = std::max(simulator.episodeIndex(), 0);
const auto& oilVaporizationControl = simulator.vanguard().schedule()[episodeIdx].oilvap();
if (drsdtConvective_()) {
// 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
Scalar g = this->gravity_[dim - 1];
ElementContext elemCtx(simulator);
const auto& vanguard = simulator.vanguard();
auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element& elem = *elemIt;
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
unsigned compressedDofIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const DimMatrix& perm = intrinsicPermeability(compressedDofIdx);
const Scalar permz = perm[dim - 1][dim - 1]; // The Z permeability
Scalar distZ = vanguard.cellThickness(compressedDofIdx);
const auto& iq = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = iq.fluidState();
Scalar t = Opm::getValue(fs.temperature(FluidSystem::oilPhaseIdx));
Scalar p = Opm::getValue(fs.pressure(FluidSystem::oilPhaseIdx));
Scalar so = Opm::getValue(fs.saturation(FluidSystem::oilPhaseIdx));
Scalar rssat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(),t,p);
Scalar saturatedDensity = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(),t,p);
Scalar rsZero = 0.0;
Scalar pureDensity = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(),t,p,rsZero);
Scalar deltaDensity = saturatedDensity-pureDensity;
Scalar rs = Opm::getValue(fs.Rs());
Scalar visc = FluidSystem::oilPvt().viscosity(fs.pvtRegionIndex(),t,p,rs);
Scalar poro = Opm::getValue(iq.porosity());
// 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
convectiveDrs_[compressedDofIdx] = permz * rssat * Opm::max(0.0, deltaDensity) * g / ( so * visc * distZ * poro);
}
}
if (oilVaporizationControl.drsdtActive()) {
ElementContext elemCtx(simulator);
const auto& vanguard = simulator.vanguard();
@ -3384,6 +3443,7 @@ private:
std::vector<bool> dRsDtOnlyFreeGas_; // apply the DRSDT rate limit only to cells that exhibit free gas
std::vector<Scalar> lastRs_;
std::vector<Scalar> maxDRs_;
std::vector<Scalar> convectiveDrs_;
std::vector<Scalar> lastRv_;
std::vector<Scalar> maxDRv_;
constexpr static Scalar freeGasMinSaturation_ = 1e-7;