From 3dd7fd0b3a7dd3fe7b44b55aae286965afc06d93 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Thu, 1 Jun 2017 07:51:44 +0200 Subject: [PATCH] Add polymer model to ebos Adds a conservation equation for polymer. Polymer concentration in the water phase is used as primary variable The polymer influences the viscosity of the water, and leaves gas and oil uneffected. A shear multiplier is computed if PLYSHLOG and/or SHRATE is specified based on either velocity or shrate. The shear multiplier effects the water and polymer viscosity. Tested and verified on the test cases in polymer_test_suite --- ebos/eclfluxmodule.hh | 3 + ebos/eclproblem.hh | 128 +++++++++++++++++++++++++++++++++++++++++- 2 files changed, 129 insertions(+), 2 deletions(-) diff --git a/ebos/eclfluxmodule.hh b/ebos/eclfluxmodule.hh index 22d101770..2695092c6 100644 --- a/ebos/eclfluxmodule.hh +++ b/ebos/eclfluxmodule.hh @@ -209,6 +209,9 @@ protected: void updateSolvent(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx) { asImp_().updateVolumeFluxTrans(elemCtx, scvfIdx, timeIdx); } + void updatePolymer(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx) + { asImp_().updateShearMultipliers(elemCtx, scvfIdx, timeIdx); } + /*! * \brief Update the required gradients for interior faces */ diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index dae5311fd..96b863932 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -249,6 +249,7 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) enum { numPhases = FluidSystem::numPhases }; enum { numComponents = FluidSystem::numComponents }; enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) }; + enum { enablePolymer = GET_PROP_VALUE(TypeTag, EnablePolymer) }; enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; @@ -267,8 +268,12 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem) typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams; typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + typedef BlackOilSolventModule SolventModule; + typedef BlackOilPolymerModule PolymerModule; + typedef Opm::CompositionalFluidState ScalarFluidState; typedef Opm::MathToolbox Toolbox; typedef Ewoms::EclSummaryWriter EclSummaryWriter; @@ -321,6 +326,8 @@ public: // Tell the solvent module to initialize its internal data structures const auto& gridManager = simulator.gridManager(); SolventModule::initFromDeck(gridManager.deck(), gridManager.eclState()); + PolymerModule::initFromDeck(gridManager.deck(), gridManager.eclState()); + } /*! @@ -360,6 +367,15 @@ public: simulator.setTimeStepSize(0.0); updatePffDofData_(); + + if (GET_PROP_VALUE(TypeTag, EnablePolymer)) { + const auto& gridManager = this->simulator().gridManager(); + const auto& gridView = gridManager.gridView(); + int numElements = gridView.size(/*codim=*/0); + maxPolymerAdsorption_.resize(numElements, 0.0); + } + + } void prefetch(const Element& elem) const @@ -453,7 +469,11 @@ public: if (updateHysteresis_()) this->model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0); + this->model().updateMaxOilSaturations(); + + if (GET_PROP_VALUE(TypeTag, EnablePolymer)) + updateMaxPolymerAdsorption_(); if (!GET_PROP_VALUE(TypeTag, DisableWells)) // set up the wells @@ -739,6 +759,30 @@ public: std::shared_ptr materialLawManager() { return materialLawManager_; } + + + /*! + * \brief Returns the initial solvent saturation for a given a cell index + */ + Scalar solventSaturation(unsigned elemIdx) const + { + if (solventSaturation_.empty()) + return 0; + + return solventSaturation_[elemIdx]; + } + + /*! + * \brief Returns the initial polymer concentration for a given a cell index + */ + Scalar polymerConcentration(unsigned elemIdx) const + { + if (polymerConcentration_.empty()) + return 0; + + return polymerConcentration_[elemIdx]; + } + /*! * \brief Returns the index of the relevant region for thermodynmic properties */ @@ -796,6 +840,24 @@ public: return miscnum_[elemIdx]; } + /*! + * \brief Returns the max polymer adsorption value + */ + template + Scalar maxPolymerAdsorption(const Context& context, unsigned spaceIdx, unsigned timeIdx) const + { return maxPolymerAdsorption(context.globalSpaceIndex(spaceIdx, timeIdx)); } + + /*! + * \brief Returns the max polymer adsorption value + */ + Scalar maxPolymerAdsorption(unsigned elemIdx) const + { + if (maxPolymerAdsorption_.empty()) + return 0; + + return maxPolymerAdsorption_[elemIdx]; + } + /*! @@ -848,8 +910,11 @@ public: else values.assignNaive(initialFluidStates_[globalDofIdx]); - //if (enableSolvent) - // values[Indices::solventSaturationIdx] = whatever; + if (enableSolvent) + values[Indices::solventSaturationIdx] = solventSaturation_[globalDofIdx]; + + if (enablePolymer) + values[Indices::polymerConcentrationIdx] = polymerConcentration_[globalDofIdx]; } /*! @@ -1164,8 +1229,11 @@ private: readExplicitInitialCondition_(); else readEquilInitialCondition_(); + + readBlackoilExtentionsInitialConditions_(); } + void readEquilInitialCondition_() { const auto& deck = this->simulator().gridManager().deck(); @@ -1381,6 +1449,36 @@ private: } } } + void readBlackoilExtentionsInitialConditions_() + { + const auto& gridManager = this->simulator().gridManager(); + const auto& eclState = gridManager.eclState(); + size_t numDof = this->model().numGridDof(); + + if (enableSolvent) { + const std::vector& solventSaturationData = eclState.get3DProperties().getDoubleGridProperty("SSOL").getData(); + solventSaturation_.resize(numDof,0.0); + for (size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) { + size_t cartesianDofIdx = gridManager.cartesianIndex(dofIdx); + assert(0 <= cartesianDofIdx); + assert(cartesianDofIdx <= solventSaturationData.size()); + solventSaturation_[dofIdx] = solventSaturationData[cartesianDofIdx]; + } + } + + if (enablePolymer) { + const std::vector& polyConcentrationData = eclState.get3DProperties().getDoubleGridProperty("SPOLY").getData(); + polymerConcentration_.resize(numDof,0.0); + for (size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) { + size_t cartesianDofIdx = gridManager.cartesianIndex(dofIdx); + assert(0 <= cartesianDofIdx); + assert(cartesianDofIdx <= polyConcentrationData.size()); + polymerConcentration_[dofIdx] = polyConcentrationData[cartesianDofIdx]; + } + } + } + + // update the hysteresis parameters of the material laws for the whole grid bool updateHysteresis_() @@ -1407,6 +1505,26 @@ private: return true; } + void updateMaxPolymerAdsorption_() + { + // we need to update the max polymer adsoption data for all elements + ElementContext elemCtx(this->simulator()); + const auto& gridManager = this->simulator().gridManager(); + auto elemIt = gridManager.gridView().template begin(); + const auto& elemEndIt = gridManager.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 auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + + maxPolymerAdsorption_[compressedDofIdx] = std::max(maxPolymerAdsorption_[compressedDofIdx] , intQuants.polymerAdsorption().value()); + } + } + void updatePvtnum_() { const auto& eclState = this->simulator().gridManager().eclState(); @@ -1513,9 +1631,15 @@ private: std::vector rockTableIdx_; std::vector rockParams_; + std::vector maxPolymerAdsorption_; + bool useMassConservativeInitialCondition_; std::vector initialFluidStates_; + std::vector polymerConcentration_; + std::vector solventSaturation_; + + EclWellManager wellManager_; EclDeckUnits deckUnits_;