From 9aa4c415ad68b6c218f913bdcc75f7c262c0772d Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 24 Sep 2020 21:51:59 +0200 Subject: [PATCH] Implement convective dissolution rate The convective DRSDT is activated with DRSDTCON --- ebos/eclbasevanguard.hh | 56 ++++++++++++++++++++++++++++++++-- ebos/eclcpgridvanguard.hh | 27 ++++++++++++++++- ebos/eclproblem.hh | 64 +++++++++++++++++++++++++++++++++++++-- 3 files changed, 141 insertions(+), 6 deletions(-) diff --git a/ebos/eclbasevanguard.hh b/ebos/eclbasevanguard.hh index 44c8fb7db..6d104db92 100644 --- a/ebos/eclbasevanguard.hh +++ b/ebos/eclbasevanguard.hh @@ -731,7 +731,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. @@ -741,6 +741,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 @@ -793,6 +806,38 @@ protected: cellCenterDepth_[elemIdx] = cellCenterDepth(element); } } + 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(); + const auto& elemEndIt = this->gridView().template end(); + 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_() @@ -823,6 +868,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; @@ -886,11 +933,14 @@ protected: */ std::vector cartesianToCompressed_; - /*! \brief Cell center depths computed - * from averaging cell corner depths + /*! \brief Cell center depths */ std::vector cellCenterDepth_; + /*! \brief Cell thichness + */ + std::vector cellThickness_; + /*! \brief information about wells in parallel * * For each well in the model there is an entry with its name diff --git a/ebos/eclcpgridvanguard.hh b/ebos/eclcpgridvanguard.hh index 8a8e296b1..7c00f0d4f 100644 --- a/ebos/eclcpgridvanguard.hh +++ b/ebos/eclcpgridvanguard.hh @@ -99,6 +99,7 @@ public: private: typedef Dune::CartesianIndexMapper 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) { @@ -398,6 +399,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_; std::unique_ptr equilGrid_; std::unique_ptr cartesianIndexMapper_; diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index eba3a2d3e..169530117 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -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_(); @@ -1973,10 +1981,15 @@ public: if (!drsdtActive_() || maxDRs_[pvtRegionIdx] < 0.0) return std::numeric_limits::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]; } @@ -2316,6 +2329,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_() { @@ -2325,6 +2347,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(); + const auto& elemEndIt = vanguard.gridView().template end(); + 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(); @@ -3385,6 +3444,7 @@ private: std::vector dRsDtOnlyFreeGas_; // apply the DRSDT rate limit only to cells that exhibit free gas std::vector lastRs_; std::vector maxDRs_; + std::vector convectiveDrs_; std::vector lastRv_; std::vector maxDRv_; constexpr static Scalar freeGasMinSaturation_ = 1e-7;