From 9ebb3db5cc25a9dee056c506a58ac902ec666896 Mon Sep 17 00:00:00 2001 From: daavid00 Date: Wed, 6 Oct 2021 19:28:28 +0200 Subject: [PATCH] first version of micp implementation in flow --- .../blackoil/blackoilboundaryratevector.hh | 7 + .../blackoil/blackoilextensivequantities.hh | 2 + opm/models/blackoil/blackoilindices.hh | 62 +- .../blackoil/blackoilintensivequantities.hh | 15 +- opm/models/blackoil/blackoillocalresidual.hh | 11 + opm/models/blackoil/blackoilmicpmodules.hh | 675 ++++++++++++++++++ opm/models/blackoil/blackoilmodel.hh | 12 +- opm/models/blackoil/blackoilnewtonmethod.hh | 22 +- .../blackoil/blackoilonephaseindices.hh | 64 +- .../blackoil/blackoilprimaryvariables.hh | 69 +- opm/models/blackoil/blackoilproperties.hh | 3 + opm/models/blackoil/blackoilratevector.hh | 6 + .../blackoil/blackoiltwophaseindices.hh | 62 +- opm/models/io/vtkblackoilmicpmodule.hh | 264 +++++++ opm/models/io/vtkblackoilmodule.hh | 6 +- 15 files changed, 1234 insertions(+), 46 deletions(-) create mode 100644 opm/models/blackoil/blackoilmicpmodules.hh create mode 100644 opm/models/io/vtkblackoilmicpmodule.hh diff --git a/opm/models/blackoil/blackoilboundaryratevector.hh b/opm/models/blackoil/blackoilboundaryratevector.hh index ac4f43160..30526dabf 100644 --- a/opm/models/blackoil/blackoilboundaryratevector.hh +++ b/opm/models/blackoil/blackoilboundaryratevector.hh @@ -62,6 +62,7 @@ class BlackOilBoundaryRateVector : public GetPropType() }; + enum { enableMICP = getPropValue() }; static constexpr bool blackoilConserveSurfaceVolume = getPropValue(); @@ -177,6 +178,12 @@ public: (*this)[Indices::contiPolymerEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.polymerConcentration(); } + if (enableMICP) { + (*this)[Indices::contiMicrobialEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.microbialConcentration(); + (*this)[Indices::contiOxygenEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.oxygenConcentration(); + (*this)[Indices::contiUreaEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.ureaConcentration(); + } + // make sure that the right mass conservation quantities are used LocalResidual::adaptMassConservationQuantities_(*this, insideIntQuants.pvtRegionIndex()); diff --git a/opm/models/blackoil/blackoilextensivequantities.hh b/opm/models/blackoil/blackoilextensivequantities.hh index 9e741f9df..72f87ca6c 100644 --- a/opm/models/blackoil/blackoilextensivequantities.hh +++ b/opm/models/blackoil/blackoilextensivequantities.hh @@ -33,6 +33,7 @@ #include "blackoilpolymermodules.hh" #include "blackoilenergymodules.hh" #include "blackoildiffusionmodule.hh" +#include "blackoilmicpmodules.hh" #include namespace Opm { @@ -55,6 +56,7 @@ class BlackOilExtensiveQuantities , public BlackOilPolymerExtensiveQuantities , public BlackOilEnergyExtensiveQuantities , public BlackOilDiffusionExtensiveQuantities()> + , public BlackOilMICPExtensiveQuantities { using MultiPhaseParent = MultiPhaseBaseExtensiveQuantities; diff --git a/opm/models/blackoil/blackoilindices.hh b/opm/models/blackoil/blackoilindices.hh index addc41f08..6f015adc8 100644 --- a/opm/models/blackoil/blackoilindices.hh +++ b/opm/models/blackoil/blackoilindices.hh @@ -35,7 +35,7 @@ namespace Opm { * * \brief The primary variable and equation indices for the black-oil model. */ -template +template struct BlackOilIndices { //! Number of phases active at all times @@ -58,6 +58,9 @@ struct BlackOilIndices //! Shall energy be conserved? static const bool enableEnergy = numEnergyV > 0; + //! Is MICP involved? + static const bool enableMICP = numMICPsV > 0; + //! Number of solvent components to be considered static const int numSolvents = enableSolvent ? numSolventsV : 0; @@ -76,8 +79,11 @@ struct BlackOilIndices //! Number of salt equations to be considered static const int numBrine = enableBrine? 1 : 0; + //! Number of MICP components to be considered + static const int numMICPs = enableMICP ? numMICPsV : 0; + //! The number of equations - static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine; + static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine + numMICPs; //! \brief returns the index of "active" component static constexpr unsigned canonicalToActiveComponentIndex(unsigned compIdx) @@ -122,17 +128,37 @@ struct BlackOilIndices static const int polymerMoleWeightIdx = numPolymers > 1 ? polymerConcentrationIdx + 1 : -1000; + //! Index of the primary variable for the first MICP component + static const int microbialConcentrationIdx = + enableMICP ? PVOffset + numPhases + numSolvents : -1000; + + //! Index of the primary variable for the second MICP component + static const int oxygenConcentrationIdx = + numMICPs > 1 ? microbialConcentrationIdx + 1 : -1000; + + //! Index of the primary variable for the third MICP component + static const int ureaConcentrationIdx = + numMICPs > 2 ? oxygenConcentrationIdx + 1 : -1000; + + //! Index of the primary variable for the fourth MICP component + static const int biofilmConcentrationIdx = + numMICPs > 3 ? ureaConcentrationIdx + 1 : -1000; + + //! Index of the primary variable for the fifth MICP component + static const int calciteConcentrationIdx = + numMICPs > 4 ? biofilmConcentrationIdx + 1 : -1000; + //! Index of the primary variable for the foam static const int foamConcentrationIdx = - enableFoam ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers : -1000; + enableFoam ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs : -1000; //! Index of the primary variable for the brine static const int saltConcentrationIdx = - enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numExtbos + numPolymers + numFoam : -1000; + enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numExtbos + numPolymers + numMICPs + numFoam : -1000; //! Index of the primary variable for temperature static const int temperatureIdx = - enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : - 1000; + enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine : - 1000; //////// @@ -159,18 +185,38 @@ struct BlackOilIndices static const int contiPolymerMWEqIdx = numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000; + //! Index of the continuity equation for the first MICP component + static const int contiMicrobialEqIdx = + enableMICP ? PVOffset + numPhases + numSolvents + numExtbos : -1000; + + //! Index of the continuity equation for the second MICP component + static const int contiOxygenEqIdx = + numMICPs > 1 ? contiMicrobialEqIdx + 1 : -1000; + + //! Index of the continuity equation for the third MICP component + static const int contiUreaEqIdx = + numMICPs > 2 ? contiOxygenEqIdx + 1 : -1000; + + //! Index of the continuity equation for the fourth MICP component + static const int contiBiofilmEqIdx = + numMICPs > 3 ? contiUreaEqIdx + 1 : -1000; + + //! Index of the continuity equation for the fifth MICP component + static const int contiCalciteEqIdx = + numMICPs > 4 ? contiBiofilmEqIdx + 1 : -1000; + //! Index of the continuity equation for the foam component static const int contiFoamEqIdx = - enableFoam ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers : -1000; + enableFoam ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs : -1000; //! Index of the continuity equation for the salt water component static const int contiBrineEqIdx = - enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam : -1000; + enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam : -1000; //! Index of the continuity equation for energy static const int contiEnergyEqIdx = - enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine: -1000; + enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine: -1000; }; } // namespace Opm diff --git a/opm/models/blackoil/blackoilintensivequantities.hh b/opm/models/blackoil/blackoilintensivequantities.hh index c1a50c0ff..3e0495695 100644 --- a/opm/models/blackoil/blackoilintensivequantities.hh +++ b/opm/models/blackoil/blackoilintensivequantities.hh @@ -36,6 +36,7 @@ #include "blackoilbrinemodules.hh" #include "blackoilenergymodules.hh" #include "blackoildiffusionmodule.hh" +#include "blackoilmicpmodules.hh" #include #include #include @@ -62,6 +63,7 @@ class BlackOilIntensiveQuantities , public BlackOilFoamIntensiveQuantities , public BlackOilBrineIntensiveQuantities , public BlackOilEnergyIntensiveQuantities + , public BlackOilMICPIntensiveQuantities { using ParentType = GetPropType; using Implementation = GetPropType; @@ -85,6 +87,7 @@ class BlackOilIntensiveQuantities enum { enableTemperature = getPropValue() }; enum { enableEnergy = getPropValue() }; enum { enableDiffusion = getPropValue() }; + enum { enableMICP = getPropValue() }; enum { numPhases = getPropValue() }; enum { numComponents = getPropValue() }; enum { waterCompIdx = FluidSystem::waterCompIdx }; @@ -129,6 +132,7 @@ public: const auto& problem = elemCtx.problem(); const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx); + const auto& linearizationType = elemCtx.linearizationType(); asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx); @@ -170,7 +174,7 @@ public: } } if (gasEnabled && waterEnabled && !oilEnabled) { - Sg = 1.0 - Sw; + Sg = 1.0 - Sw; } Valgrind::CheckDefined(Sg); @@ -384,11 +388,19 @@ public: // deal with water induced rock compaction porosity_ *= problem.template rockCompPoroMultiplier(*this, globalSpaceIdx); + // the MICP processes change the porosity + if (enableMICP){ + Evaluation biofilm_ = priVars.makeEvaluation(Indices::biofilmConcentrationIdx, timeIdx, linearizationType); + Evaluation calcite_ = priVars.makeEvaluation(Indices::calciteConcentrationIdx, timeIdx, linearizationType); + porosity_ += - biofilm_ - calcite_; + } + asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx); asImp_().zPvtUpdate_(); asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx); asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache); asImp_().foamPropertiesUpdate_(elemCtx, dofIdx, timeIdx); + asImp_().MICPPropertiesUpdate_(elemCtx, dofIdx, timeIdx); // update the quantities which are required by the chosen // velocity model @@ -474,6 +486,7 @@ private: friend BlackOilEnergyIntensiveQuantities; friend BlackOilFoamIntensiveQuantities; friend BlackOilBrineIntensiveQuantities; + friend BlackOilMICPIntensiveQuantities; Implementation& asImp_() diff --git a/opm/models/blackoil/blackoillocalresidual.hh b/opm/models/blackoil/blackoillocalresidual.hh index c818c3422..4288c2183 100644 --- a/opm/models/blackoil/blackoillocalresidual.hh +++ b/opm/models/blackoil/blackoillocalresidual.hh @@ -36,6 +36,7 @@ #include "blackoilfoammodules.hh" #include "blackoilbrinemodules.hh" #include "blackoildiffusionmodule.hh" +#include "blackoilmicpmodules.hh" #include namespace Opm { @@ -88,6 +89,7 @@ class BlackOilLocalResidual : public GetPropType; using BrineModule = BlackOilBrineModule; using DiffusionModule = BlackOilDiffusionModule; + using MICPModule = BlackOilMICPModule; public: /*! @@ -161,6 +163,9 @@ public: // deal with salt (if present) BrineModule::addStorage(storage, intQuants); + + // deal with micp (if present) + MICPModule::addStorage(storage, intQuants); } /*! @@ -208,6 +213,9 @@ public: // deal with salt (if present) BrineModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx); + // deal with micp (if present) + MICPModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx); + DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx); } @@ -222,6 +230,9 @@ public: // retrieve the source term intrinsic to the problem elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx); + // deal with MICP (if present) + MICPModule::addSource(source, elemCtx, dofIdx, timeIdx); + // scale the source term of the energy equation if (enableEnergy) source[Indices::contiEnergyEqIdx] *= getPropValue(); diff --git a/opm/models/blackoil/blackoilmicpmodules.hh b/opm/models/blackoil/blackoilmicpmodules.hh new file mode 100644 index 000000000..e5cc9314e --- /dev/null +++ b/opm/models/blackoil/blackoilmicpmodules.hh @@ -0,0 +1,675 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief Contains the classes required to extend the black-oil model by MICP. + */ +#ifndef EWOMS_BLACK_OIL_MICP_MODULE_HH +#define EWOMS_BLACK_OIL_MICP_MODULE_HH + +#include "blackoilproperties.hh" +#include +#include + +#if HAVE_ECL_INPUT +#include +#include +#endif + +#include +#include + +#include + +#include + +namespace Opm { +/*! + * \ingroup BlackOil + * \brief Contains the high level supplements required to extend the black oil + * model by MICP. + */ +template ()> +class BlackOilMICPModule +{ + using Scalar = GetPropType; + using Evaluation = GetPropType; + using PrimaryVariables = GetPropType; + using IntensiveQuantities = GetPropType; + using ExtensiveQuantities = GetPropType; + using ElementContext = GetPropType; + using FluidSystem = GetPropType; + using Model = GetPropType; + using Simulator = GetPropType; + using EqVector = GetPropType; + using RateVector = GetPropType; + using Indices = GetPropType; + + using Toolbox = MathToolbox; + + static constexpr unsigned microbialConcentrationIdx = Indices::microbialConcentrationIdx; + static constexpr unsigned oxygenConcentrationIdx = Indices::oxygenConcentrationIdx; + static constexpr unsigned ureaConcentrationIdx = Indices::ureaConcentrationIdx; + static constexpr unsigned biofilmConcentrationIdx = Indices::biofilmConcentrationIdx; + static constexpr unsigned calciteConcentrationIdx = Indices::calciteConcentrationIdx; + static constexpr unsigned contiMicrobialEqIdx = Indices::contiMicrobialEqIdx; + static constexpr unsigned contiOxygenEqIdx = Indices::contiOxygenEqIdx; + static constexpr unsigned contiUreaEqIdx = Indices::contiUreaEqIdx; + static constexpr unsigned contiBiofilmEqIdx = Indices::contiBiofilmEqIdx; + static constexpr unsigned contiCalciteEqIdx = Indices::contiCalciteEqIdx; + static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx; + + static constexpr unsigned enableMICP = enableMICPV; + + static constexpr unsigned numEq = getPropValue(); + +public: + +#if HAVE_ECL_INPUT + // + //* \brief Initialize all internal data structures needed by the MICP module + // + static void initFromState(const EclipseState& eclState) + { + // some sanity checks: if MICP is enabled, the MICP keyword must be + // present, if MICP is disabled the keyword must not be present. + if (enableMICP && !eclState.runspec().micp()) { + throw std::runtime_error("Non-trivial MICP treatment requested at compile time, but " + "the deck does not contain the MICP keyword"); + } + else if (!enableMICP && eclState.runspec().micp()) { + throw std::runtime_error("MICP treatment disabled at compile time, but the deck " + "contains the MICP keyword"); + } + + if (!eclState.runspec().micp()) + return; // MICP treatment is supposed to be disabled*/ + + // initialize the objects which deal with the MICPpara keyword + const auto& MICPpara = eclState.getMICPpara(); + setMICPpara(MICPpara.getDensityBiofilm(), + MICPpara.getDensityCalcite(), + MICPpara.getDetachmentRate(), + MICPpara.getCriticalPorosity(), + MICPpara.getFittingFactor(), + MICPpara.getHalfVelocityOxygen(), + MICPpara.getHalfVelocityUrea(), + MICPpara.getMaximumGrowthRate(), + MICPpara.getMaximumUreaUtilization(), + MICPpara.getMicrobialAttachmentRate(), + MICPpara.getMicrobialDeathRate(), + MICPpara.getMinimumPermeability(), + MICPpara.getOxygenConsumptionFactor(), + MICPpara.getYieldGrowthCoefficient(), + MICPpara.getMaximumOxygenConcentration(), + MICPpara.getMaximumUreaConcentration(), + MICPpara.getToleranceBeforeClogging()); + // obtain the porosity for the clamp in the blackoilnewtonmethod + phi_ = eclState.fieldProps().get_double("PORO"); + } +#endif + + /*! + * \brief The simulator stops if "clogging" has been (almost) reached in any of the cells. + * + * I.e., porosity - biofilm - calcite < tol_clgg, where tol_clgg is a given tolerance. In the + * implemented model a permebaility-porosity relatonship is used where a minimum + * permeability value is reached if porosity - biofilm - calcite < phi_crit. + */ + static void checkCloggingMICP(const Model& model, const Scalar phi, unsigned dofIdx) + { + const PrimaryVariables& priVars = model.solution(/*timeIdx=*/1)[dofIdx]; + if (phi - priVars[biofilmConcentrationIdx] - priVars[calciteConcentrationIdx] < MICPparaToleranceBeforeClogging()) + throw std::logic_error("Clogging has been (almost) reached in at least one cell\n"); + } + + /*! + * \brief Specify the MICP properties a single region. + * + * The index of specified here must be in range [0, numSatRegions) + */ + static void setMICPpara(const Scalar& MICPparaDensityBiofilm, + const Scalar& MICPparaDensityCalcite, + const Scalar& MICPparaDetachmentRate, + const Scalar& MICPparaCriticalPorosity, + const Scalar& MICPparaFittingFactor, + const Scalar& MICPparaHalfVelocityOxygen, + const Scalar& MICPparaHalfVelocityUrea, + const Scalar& MICPparaMaximumGrowthRate, + const Scalar& MICPparaMaximumUreaUtilization, + const Scalar& MICPparaMicrobialAttachmentRate, + const Scalar& MICPparaMicrobialDeathRate, + const Scalar& MICPparaMinimumPermeability, + const Scalar& MICPparaOxygenConsumptionFactor, + const Scalar& MICPparaYieldGrowthCoefficient, + const Scalar& MICPparaMaximumOxygenConcentration, + const Scalar& MICPparaMaximumUreaConcentration, + const Scalar& MICPparaToleranceBeforeClogging) + { + MICPparaDensityBiofilm_ = MICPparaDensityBiofilm; + MICPparaDensityCalcite_ = MICPparaDensityCalcite; + MICPparaDetachmentRate_ = MICPparaDetachmentRate; + MICPparaCriticalPorosity_ = MICPparaCriticalPorosity; + MICPparaFittingFactor_ = MICPparaFittingFactor; + MICPparaHalfVelocityOxygen_ = MICPparaHalfVelocityOxygen; + MICPparaHalfVelocityUrea_ = MICPparaHalfVelocityUrea; + MICPparaMaximumGrowthRate_ = MICPparaMaximumGrowthRate; + MICPparaMaximumUreaUtilization_ = MICPparaMaximumUreaUtilization; + MICPparaMicrobialAttachmentRate_ = MICPparaMicrobialAttachmentRate; + MICPparaMicrobialDeathRate_ = MICPparaMicrobialDeathRate; + MICPparaMinimumPermeability_ = MICPparaMinimumPermeability; + MICPparaOxygenConsumptionFactor_ = MICPparaOxygenConsumptionFactor; + MICPparaYieldGrowthCoefficient_ = MICPparaYieldGrowthCoefficient; + MICPparaMaximumOxygenConcentration_ = MICPparaMaximumOxygenConcentration; + MICPparaMaximumUreaConcentration_ = MICPparaMaximumUreaConcentration; + MICPparaToleranceBeforeClogging_ = MICPparaToleranceBeforeClogging; + } + + /*! + * \brief Register all run-time parameters for the black-oil MICP module. + */ + static void registerParameters() + { + if (!enableMICP) + // MICP has been disabled at compile time + return; + + VtkBlackOilMICPModule::registerParameters(); + } + + /*! + * \brief Register all MICP specific VTK and ECL output modules. + */ + static void registerOutputModules(Model& model, + Simulator& simulator) + { + if (!enableMICP) + // MICP has been disabled at compile time + return; + + model.addOutputModule(new VtkBlackOilMICPModule(simulator)); + } + + static bool eqApplies(unsigned eqIdx) + { + if (!enableMICP) + return false; + + // All MICP components are true here + return eqIdx == contiMicrobialEqIdx || eqIdx == contiOxygenEqIdx || eqIdx == contiUreaEqIdx || eqIdx == contiBiofilmEqIdx || eqIdx == contiCalciteEqIdx; + } + + static Scalar eqWeight(unsigned eqIdx OPM_OPTIM_UNUSED) + { + assert(eqApplies(eqIdx)); + + // TODO: it may be beneficial to chose this differently. + return static_cast(1.0); + } + + // must be called after water storage is computed + template + static void addStorage(Dune::FieldVector& storage, + const IntensiveQuantities& intQuants) + { + if (!enableMICP) + return; + + LhsEval surfaceVolumeWater = Toolbox::template decay(intQuants.porosity()); + // avoid singular matrix if no water is present. + surfaceVolumeWater = max(surfaceVolumeWater, 1e-10); + // Suspended microbes in water phase + const LhsEval massMicrobes = surfaceVolumeWater * Toolbox::template decay(intQuants.microbialConcentration()); + LhsEval accumulationMicrobes = massMicrobes; + storage[contiMicrobialEqIdx] += accumulationMicrobes; + // Oxygen in water phase + const LhsEval massOxygen = surfaceVolumeWater * Toolbox::template decay(intQuants.oxygenConcentration()); + LhsEval accumulationOxygen = massOxygen; + storage[contiOxygenEqIdx] += accumulationOxygen; + // Urea in water phase + const LhsEval massUrea = surfaceVolumeWater * Toolbox::template decay(intQuants.ureaConcentration()); + LhsEval accumulationUrea = massUrea; + storage[contiUreaEqIdx] += accumulationUrea; + // Biofilm + const LhsEval massBiofilm = Toolbox::template decay(intQuants.biofilmConcentration()); + LhsEval accumulationBiofilm = massBiofilm; + storage[contiBiofilmEqIdx] += accumulationBiofilm; + // Calcite + const LhsEval massCalcite = Toolbox::template decay(intQuants.calciteConcentration()); + LhsEval accumulationCalcite = massCalcite; + storage[contiCalciteEqIdx] += accumulationCalcite; + } + + static void computeFlux(RateVector& flux, + const ElementContext& elemCtx, + unsigned scvfIdx, + unsigned timeIdx) + + { + if (!enableMICP) + return; + + const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx); + + const unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx); + const unsigned inIdx = extQuants.interiorIndex(); + const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx); + + + if (upIdx == inIdx) { + flux[contiMicrobialEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * up.microbialConcentration(); + flux[contiOxygenEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * up.oxygenConcentration(); + flux[contiUreaEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * up.ureaConcentration(); + } + else { + flux[contiMicrobialEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * decay(up.microbialConcentration()); + flux[contiOxygenEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * decay(up.oxygenConcentration()); + flux[contiUreaEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * decay(up.ureaConcentration()); + } + } + + // See https://doi.org/10.1016/j.ijggc.2021.103256 for the micp processes in the model. + static void addSource(RateVector& source, + const ElementContext& elemCtx, + unsigned dofIdx, + unsigned timeIdx) + + { + if (!enableMICP) + return; + + // compute dpW (max norm of the pressure gradient in the cell center) + const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx); + const auto& K = elemCtx.problem().intrinsicPermeability(elemCtx, dofIdx, 0); + size_t numInteriorFaces = elemCtx.numInteriorFaces(timeIdx); + Evaluation dpW = 0; + for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) { + const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx); + unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx); + const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx); + const Evaluation& mobWater = up.mobility(waterPhaseIdx); + + // compute water velocity from flux + Evaluation waterVolumeVelocity = extQuants.volumeFlux(waterPhaseIdx) / (K[0][0] * mobWater); + dpW = std::max(dpW, abs(waterVolumeVelocity)); + } + + // get the model parameters + Scalar k_a = MICPparaMicrobialAttachmentRate(); + Scalar k_d = MICPparaMicrobialDeathRate(); + Scalar rho_b = MICPparaDensityBiofilm(); + Scalar rho_c = MICPparaDensityCalcite(); + Scalar k_str = MICPparaDetachmentRate(); + Scalar k_o = MICPparaHalfVelocityOxygen(); + Scalar k_u = MICPparaHalfVelocityUrea() / 10.0;//Dividing by scaling factor 10 (see WellInterface_impl.hpp) + Scalar mu = MICPparaMaximumGrowthRate(); + Scalar mu_u = MICPparaMaximumUreaUtilization() / 10.0;//Dividing by scaling factor 10 (see WellInterface_impl.hpp) + Scalar Y_sb = MICPparaYieldGrowthCoefficient(); + Scalar F = MICPparaOxygenConsumptionFactor(); + Scalar Y_uc = 1.67 * 10; //Multiplying by scaling factor 10 (see WellInterface_impl.hpp) + + // compute the processes + source[Indices::contiMicrobialEqIdx] += intQuants.microbialConcentration() * intQuants.porosity() * + (Y_sb * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()) - k_d - k_a) + + rho_b * intQuants.biofilmConcentration() * k_str * pow(intQuants.porosity() * dpW, 0.58); + + source[Indices::contiOxygenEqIdx] -= (intQuants.microbialConcentration() * intQuants.porosity() + rho_b * intQuants.biofilmConcentration()) * + F * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()); + + source[Indices::contiUreaEqIdx] -= rho_b * intQuants.biofilmConcentration() * mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration()); + + source[Indices::contiBiofilmEqIdx] += intQuants.biofilmConcentration() * (Y_sb * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()) - k_d + - k_str * pow(intQuants.porosity() * dpW, 0.58) - Y_uc * (rho_b / rho_c) * intQuants.biofilmConcentration() * mu_u * + (intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration())) / (intQuants.porosity() + intQuants.biofilmConcentration())) + + k_a * intQuants.microbialConcentration() * intQuants.porosity() / rho_b; + + source[Indices::contiCalciteEqIdx] += (rho_b / rho_c) * intQuants.biofilmConcentration() * Y_uc * mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration()); + } + + static const Scalar MICPparaDensityBiofilm() + { + return MICPparaDensityBiofilm_; + } + + static const Scalar MICPparaDensityCalcite() + { + return MICPparaDensityCalcite_; + } + + static const Scalar MICPparaDetachmentRate() + { + return MICPparaDetachmentRate_; + } + + static const Scalar MICPparaCriticalPorosity() + { + return MICPparaCriticalPorosity_; + } + + static const Scalar MICPparaFittingFactor() + { + return MICPparaFittingFactor_; + } + + static const Scalar MICPparaHalfVelocityOxygen() + { + return MICPparaHalfVelocityOxygen_; + } + + static const Scalar MICPparaHalfVelocityUrea() + { + return MICPparaHalfVelocityUrea_; + } + + static const Scalar MICPparaMaximumGrowthRate() + { + return MICPparaMaximumGrowthRate_; + } + + static const Scalar MICPparaMaximumOxygenConcentration() + { + return MICPparaMaximumOxygenConcentration_; + } + + static const Scalar MICPparaMaximumUreaConcentration() + { + return MICPparaMaximumUreaConcentration_ / 10.0;//Dividing by scaling factor 10 (see WellInterface_impl.hpp); + } + + static const Scalar MICPparaMaximumUreaUtilization() + { + return MICPparaMaximumUreaUtilization_; + } + + static const Scalar MICPparaMicrobialAttachmentRate() + { + return MICPparaMicrobialAttachmentRate_; + } + + static const Scalar MICPparaMicrobialDeathRate() + { + return MICPparaMicrobialDeathRate_; + } + + static const Scalar MICPparaMinimumPermeability() + { + return MICPparaMinimumPermeability_; + } + + static const Scalar MICPparaOxygenConsumptionFactor() + { + return MICPparaOxygenConsumptionFactor_; + } + + static const Scalar MICPparaToleranceBeforeClogging() + { + return MICPparaToleranceBeforeClogging_; + } + + static const Scalar MICPparaYieldGrowthCoefficient() + { + return MICPparaYieldGrowthCoefficient_; + } + + static const std::vector phi() + { + return phi_; + } + +private: + static Scalar MICPparaDensityBiofilm_; + static Scalar MICPparaDensityCalcite_; + static Scalar MICPparaDetachmentRate_; + static Scalar MICPparaCriticalPorosity_; + static Scalar MICPparaFittingFactor_; + static Scalar MICPparaHalfVelocityOxygen_; + static Scalar MICPparaHalfVelocityUrea_; + static Scalar MICPparaMaximumGrowthRate_; + static Scalar MICPparaMaximumUreaUtilization_; + static Scalar MICPparaMicrobialAttachmentRate_; + static Scalar MICPparaMicrobialDeathRate_; + static Scalar MICPparaMinimumPermeability_; + static Scalar MICPparaOxygenConsumptionFactor_; + static Scalar MICPparaYieldGrowthCoefficient_; + static Scalar MICPparaMaximumOxygenConcentration_; + static Scalar MICPparaMaximumUreaConcentration_; + static Scalar MICPparaToleranceBeforeClogging_; + static std::vector phi_; +}; + + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaDensityBiofilm_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaDensityCalcite_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaDetachmentRate_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaCriticalPorosity_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaFittingFactor_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaHalfVelocityOxygen_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaHalfVelocityUrea_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaMaximumGrowthRate_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaMaximumUreaUtilization_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaMicrobialAttachmentRate_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaMicrobialDeathRate_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaMinimumPermeability_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaOxygenConsumptionFactor_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaYieldGrowthCoefficient_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaMaximumOxygenConcentration_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaMaximumUreaConcentration_; + +template +typename BlackOilMICPModule::Scalar +BlackOilMICPModule::MICPparaToleranceBeforeClogging_; + +template +std::vector::Scalar> +BlackOilMICPModule::phi_; + + +/*! + * \ingroup BlackOil + * \class Opm::BlackOilMICPIntensiveQuantities + * + * \brief Provides the volumetric quantities required for the equations needed by the + * MICP extension of the black-oil model. + */ +template ()> +class BlackOilMICPIntensiveQuantities +{ + using Implementation = GetPropType; + + using Scalar = GetPropType; + using Evaluation = GetPropType; + using PrimaryVariables = GetPropType; + using FluidSystem = GetPropType; + using Indices = GetPropType; + using ElementContext = GetPropType; + + using MICPModule = BlackOilMICPModule; + + static constexpr int microbialConcentrationIdx = Indices::microbialConcentrationIdx; + static constexpr int oxygenConcentrationIdx = Indices::oxygenConcentrationIdx; + static constexpr int ureaConcentrationIdx = Indices::ureaConcentrationIdx; + static constexpr int biofilmConcentrationIdx = Indices::biofilmConcentrationIdx; + static constexpr int calciteConcentrationIdx = Indices::calciteConcentrationIdx; + static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx; + + +public: + + /*! + * \brief Update the intensive properties needed to handle MICP from the + * primary variables + * + */ + void MICPPropertiesUpdate_(const ElementContext& elemCtx, + unsigned dofIdx, + unsigned timeIdx) + { + const auto linearizationType = elemCtx.linearizationType(); + const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx); + const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx); + const auto& K = elemCtx.problem().intrinsicPermeability(elemCtx, dofIdx, timeIdx); + Scalar referencePorosity_ = elemCtx.problem().porosity(elemCtx, dofIdx, timeIdx); + Scalar eta = MICPModule::MICPparaFittingFactor(); + Scalar k_min = MICPModule::MICPparaMinimumPermeability(); + Scalar phi_crit = MICPModule::MICPparaCriticalPorosity(); + + microbialConcentration_ = priVars.makeEvaluation(microbialConcentrationIdx, timeIdx, linearizationType); + oxygenConcentration_ = priVars.makeEvaluation(oxygenConcentrationIdx, timeIdx, linearizationType); + ureaConcentration_ = priVars.makeEvaluation(ureaConcentrationIdx, timeIdx, linearizationType); + biofilmConcentration_ = priVars.makeEvaluation(biofilmConcentrationIdx, timeIdx, linearizationType); + calciteConcentration_ = priVars.makeEvaluation(calciteConcentrationIdx, timeIdx, linearizationType); + + // Permeability reduction due to MICP, by adjusting the water mobility + asImp_().mobility_[waterPhaseIdx] *= max((pow((intQuants.porosity() - phi_crit) / (referencePorosity_ - phi_crit), eta) + k_min / K[0][0])/(1. + k_min / K[0][0]), k_min / K[0][0]); + + } + + const Evaluation& microbialConcentration() const + { return microbialConcentration_; } + + const Evaluation& oxygenConcentration() const + { return oxygenConcentration_; } + + const Evaluation& ureaConcentration() const + { return ureaConcentration_; } + + const Evaluation& biofilmConcentration() const + { return biofilmConcentration_; } + + const Evaluation& calciteConcentration() const + { return calciteConcentration_; } + + +protected: + Implementation& asImp_() + { return *static_cast(this); } + + Evaluation microbialConcentration_; + Evaluation oxygenConcentration_; + Evaluation ureaConcentration_; + Evaluation biofilmConcentration_; + Evaluation calciteConcentration_; + +}; + +template +class BlackOilMICPIntensiveQuantities +{ + using Evaluation = GetPropType; + using ElementContext = GetPropType; + using Scalar = GetPropType; + +public: + void MICPPropertiesUpdate_(const ElementContext&, + unsigned, + unsigned) + { } + + const Evaluation& microbialConcentration() const + { throw std::logic_error("microbialConcentration() called but MICP is disabled"); } + + const Evaluation& oxygenConcentration() const + { throw std::logic_error("oxygenConcentration() called but MICP is disabled"); } + + const Evaluation& ureaConcentration() const + { throw std::logic_error("ureaConcentration() called but MICP is disabled"); } + + const Evaluation& biofilmConcentration() const + { throw std::logic_error("biofilmConcentration() called but MICP is disabled"); } + + const Evaluation& calciteConcentration() const + { throw std::logic_error("calciteConcentration() called but MICP is disabled"); } +}; + +/*! + * \ingroup BlackOil + * \class Opm::BlackOilMICPExtensiveQuantities + * + * \brief Provides the MICP specific extensive quantities to the generic black-oil + * module's extensive quantities. + */ +template ()> +class BlackOilMICPExtensiveQuantities +{ + using Implementation = GetPropType; + +private: + Implementation& asImp_() + { return *static_cast(this); } + +}; + +template +class BlackOilMICPExtensiveQuantities{}; + +} // namespace Opm + +#endif diff --git a/opm/models/blackoil/blackoilmodel.hh b/opm/models/blackoil/blackoilmodel.hh index cfcf37d44..555373555 100644 --- a/opm/models/blackoil/blackoilmodel.hh +++ b/opm/models/blackoil/blackoilmodel.hh @@ -47,6 +47,7 @@ #include "blackoilbrinemodules.hh" #include "blackoilextbomodules.hh" #include "blackoildarcyfluxmodule.hh" +#include "blackoilmicpmodules.hh" #include #include @@ -78,7 +79,8 @@ struct BlackOilModel { using InheritsFrom = std::tuple; }; + MultiPhaseBaseModel, + VtkBlackOilMICP>; }; } // namespace TTag //! Set the local residual function @@ -131,7 +133,8 @@ struct Indices getPropValue(), getPropValue(), getPropValue(), - /*PVOffset=*/0>; }; + /*PVOffset=*/0, + getPropValue()>; }; //! Set the fluid system to the black-oil fluid system by default template @@ -158,6 +161,8 @@ template struct EnableFoam { static constexpr bool value = false; }; template struct EnableBrine { static constexpr bool value = false; }; +template +struct EnableMICP { static constexpr bool value = false; }; //! By default, the blackoil model is isothermal and does not conserve energy template @@ -287,6 +292,7 @@ class BlackOilModel using PolymerModule = BlackOilPolymerModule; using EnergyModule = BlackOilEnergyModule; using DiffusionModule = BlackOilDiffusionModule; + using MICPModule = BlackOilMICPModule; public: BlackOilModel(Simulator& simulator) @@ -305,6 +311,7 @@ public: PolymerModule::registerParameters(); EnergyModule::registerParameters(); DiffusionModule::registerParameters(); + MICPModule::registerParameters(); // register runtime parameters of the VTK output modules VtkBlackOilModule::registerParameters(); @@ -595,6 +602,7 @@ protected: SolventModule::registerOutputModules(*this, this->simulator_); PolymerModule::registerOutputModules(*this, this->simulator_); EnergyModule::registerOutputModules(*this, this->simulator_); + MICPModule::registerOutputModules(*this, this->simulator_); this->addOutputModule(new VtkBlackOilModule(this->simulator_)); this->addOutputModule(new VtkCompositionModule(this->simulator_)); diff --git a/opm/models/blackoil/blackoilnewtonmethod.hh b/opm/models/blackoil/blackoilnewtonmethod.hh index f7a554eb1..8041796be 100644 --- a/opm/models/blackoil/blackoilnewtonmethod.hh +++ b/opm/models/blackoil/blackoilnewtonmethod.hh @@ -32,6 +32,7 @@ #include #include +#include "blackoilmicpmodules.hh" #include @@ -111,8 +112,10 @@ class BlackOilNewtonMethod : public GetPropType; using EqVector = GetPropType; using Indices = GetPropType; + using FluidSystem = GetPropType; using Scalar = GetPropType; using Linearizer = GetPropType; + using MICPModule = BlackOilMICPModule; static const unsigned numEq = getPropValue(); @@ -248,6 +251,7 @@ protected: static constexpr bool enableFoam = Indices::foamConcentrationIdx >= 0; static constexpr bool enableBrine = Indices::saltConcentrationIdx >= 0; static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0; + static constexpr bool enableMICP = Indices::microbialConcentrationIdx >= 0; currentValue.checkDefined(); Valgrind::CheckDefined(update); @@ -259,7 +263,7 @@ protected: Scalar deltaSg = 0.0; Scalar deltaSs = 0.0; - if (Indices::waterEnabled) { + if (Indices::waterEnabled && FluidSystem::numActivePhases() > 1) { deltaSw = update[Indices::waterSaturationIdx]; deltaSo = -deltaSw; } @@ -371,6 +375,22 @@ protected: // keep the temperature within given values if (enableEnergy && pvIdx == Indices::temperatureIdx) nextValue[pvIdx] = std::clamp(nextValue[pvIdx], tempMin_, tempMax_); + + // Limit the variables to [0, cmax] values to improve the convergence. + // For the microorganisms we set this value equal to the biomass density value. + // For the oxygen and urea we set this value to the maximum injected + // concentration (the urea concentration has been scaled by 10). For + // the biofilm and calcite, we set this value equal to the porosity minus the clogging tolerance. + if (enableMICP && pvIdx == Indices::microbialConcentrationIdx) + nextValue[pvIdx] = std::clamp(nextValue[pvIdx], 0.0, MICPModule::MICPparaDensityBiofilm()); + if (enableMICP && pvIdx == Indices::oxygenConcentrationIdx) + nextValue[pvIdx] = std::clamp(nextValue[pvIdx], 0.0, MICPModule::MICPparaMaximumOxygenConcentration()); + if (enableMICP && pvIdx == Indices::ureaConcentrationIdx) + nextValue[pvIdx] = std::clamp(nextValue[pvIdx], 0.0, MICPModule::MICPparaMaximumUreaConcentration()); + if (enableMICP && pvIdx == Indices::biofilmConcentrationIdx) + nextValue[pvIdx] = std::clamp(nextValue[pvIdx], 0.0, MICPModule::phi()[globalDofIdx] - MICPModule::MICPparaToleranceBeforeClogging()); + if (enableMICP && pvIdx == Indices::calciteConcentrationIdx) + nextValue[pvIdx] = std::clamp(nextValue[pvIdx], 0.0, MICPModule::phi()[globalDofIdx] - MICPModule::MICPparaToleranceBeforeClogging()); } // switch the new primary variables to something which is physically meaningful. diff --git a/opm/models/blackoil/blackoilonephaseindices.hh b/opm/models/blackoil/blackoilonephaseindices.hh index 16c7b1549..2d6ab251b 100644 --- a/opm/models/blackoil/blackoilonephaseindices.hh +++ b/opm/models/blackoil/blackoilonephaseindices.hh @@ -23,7 +23,7 @@ /*! * \file * - * \copydoc Ewoms::BlackOilTwoPhaseIndices + * \copydoc Ewoms::BlackOilOnePhaseIndices */ #ifndef EWOMS_BLACK_OIL_ONE_PHASE_INDICES_HH #define EWOMS_BLACK_OIL_ONE_PHASE_INDICES_HH @@ -37,7 +37,7 @@ namespace Opm { * * \brief The primary variable and equation indices for the black-oil model. */ -template +template struct BlackOilOnePhaseIndices { //! Is phase enabled or not @@ -57,6 +57,9 @@ struct BlackOilOnePhaseIndices //! Shall energy be conserved? static const bool enableEnergy = numEnergyV > 0; + //! Is MICP involved? + static const bool enableMICP = numMICPsV > 0; + //! Number of solvent components to be considered static const int numSolvents = enableSolvent ? numSolventsV : 0; @@ -78,8 +81,11 @@ struct BlackOilOnePhaseIndices //! The number of fluid phases static const int numPhases = 1; + //! Number of MICP components to be considered + static const int numMICPs = enableMICP ? numMICPsV : 0; + //! The number of equations - static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine; + static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine + numMICPs; ////////////////////////////// // Primary variable indices @@ -115,17 +121,37 @@ struct BlackOilOnePhaseIndices static const int polymerMoleWeightIdx = numPolymers > 1 ? polymerConcentrationIdx + 1 : -1000; + //! Index of the primary variable for the first MICP component + static const int microbialConcentrationIdx = + enableMICP ? PVOffset + numPhases + numSolvents : -1000; + + //! Index of the primary variable for the second MICP component + static const int oxygenConcentrationIdx = + numMICPs > 1 ? microbialConcentrationIdx + 1 : -1000; + + //! Index of the primary variable for the third MICP component + static const int ureaConcentrationIdx = + numMICPs > 2 ? oxygenConcentrationIdx + 1 : -1000; + + //! Index of the primary variable for the fourth MICP component + static const int biofilmConcentrationIdx = + numMICPs > 3 ? ureaConcentrationIdx + 1 : -1000; + + //! Index of the primary variable for the fifth MICP component + static const int calciteConcentrationIdx = + numMICPs > 4 ? biofilmConcentrationIdx + 1 : -1000; + //! Index of the primary variable for the foam static const int foamConcentrationIdx = - enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000; + enableFoam ? PVOffset + numPhases + numSolvents + numPolymers + numMICPs : -1000; //! Index of the primary variable for the salt static const int saltConcentrationIdx = - enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam : -1000; + enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam : -1000; //! Index of the primary variable for temperature static const int temperatureIdx = - enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine: - 1000; + enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine: - 1000; ////////////////////// // Equation indices @@ -170,17 +196,37 @@ struct BlackOilOnePhaseIndices static const int contiPolymerMWEqIdx = numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000; + //! Index of the continuity equation for the first MICP component + static const int contiMicrobialEqIdx = + enableMICP ? PVOffset + numPhases + numSolvents : -1000; + + //! Index of the continuity equation for the second MICP component + static const int contiOxygenEqIdx = + numMICPs > 1 ? contiMicrobialEqIdx + 1 : -1000; + + //! Index of the continuity equation for the third MICP component + static const int contiUreaEqIdx = + numMICPs > 2 ? contiOxygenEqIdx + 1 : -1000; + + //! Index of the continuity equation for the fourth MICP component + static const int contiBiofilmEqIdx = + numMICPs > 3 ? contiUreaEqIdx + 1 : -1000; + + //! Index of the continuity equation for the fifth MICP component + static const int contiCalciteEqIdx = + numMICPs > 4 ? contiBiofilmEqIdx + 1 : -1000; + //! Index of the continuity equation for the foam component static const int contiFoamEqIdx = - enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000; + enableFoam ? PVOffset + numPhases + numSolvents + numPolymers + numMICPs : -1000; //! Index of the continuity equation for the salt component static const int contiBrineEqIdx = - enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam : -1000; + enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam : -1000; //! Index of the continuity equation for energy static const int contiEnergyEqIdx = - enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine: -1000; + enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine: -1000; }; } // namespace Opm diff --git a/opm/models/blackoil/blackoilprimaryvariables.hh b/opm/models/blackoil/blackoilprimaryvariables.hh index 0a4cee6a7..12971b8f2 100644 --- a/opm/models/blackoil/blackoilprimaryvariables.hh +++ b/opm/models/blackoil/blackoilprimaryvariables.hh @@ -2,20 +2,16 @@ // vi: set et ts=4 sw=4 sts=4: /* This file is part of the Open Porous Media project (OPM). - OPM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 2 of the License, or (at your option) any later version. - OPM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. - You should have received a copy of the GNU General Public License along with OPM. If not, see . - Consult the COPYING file in the top-level source directory of this module for the precise wording of the license and the list of copyright holders. @@ -35,6 +31,7 @@ #include "blackoilenergymodules.hh" #include "blackoilfoammodules.hh" #include "blackoilbrinemodules.hh" +#include "blackoilmicpmodules.hh" #include @@ -105,6 +102,7 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables enum { enableFoam = getPropValue() }; enum { enableBrine = getPropValue() }; enum { enableEnergy = getPropValue() }; + enum { enableMICP = getPropValue() }; enum { gasCompIdx = FluidSystem::gasCompIdx }; enum { waterCompIdx = FluidSystem::waterCompIdx }; enum { oilCompIdx = FluidSystem::oilCompIdx }; @@ -117,6 +115,7 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables using EnergyModule = BlackOilEnergyModule; using FoamModule = BlackOilFoamModule; using BrineModule = BlackOilBrineModule; + using MICPModule = BlackOilMICPModule; static_assert(numPhases == 3, "The black-oil model assumes three phases!"); static_assert(numComponents == 3, "The black-oil model assumes three components!"); @@ -312,13 +311,13 @@ public: } else { throw std::logic_error("For single-phase runs, only pure water is presently allowed."); } - + } else if (primaryVarsMeaning() == Sw_po_Sg) { if (waterEnabled) (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx)); - if (gasEnabled && waterEnabled && !oilEnabled) { - //-> water-gas system + if (gasEnabled && waterEnabled && !oilEnabled) { + //-> water-gas system (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(gasPhaseIdx)); } else if (oilEnabled) { @@ -603,13 +602,13 @@ public: Scalar Ssol = 0.0; if (enableSolvent) Ssol =(*this) [Indices::solventSaturationIdx]; - + Scalar So = 1.0 - Sw - Sg - Ssol; Sw = std::min(std::max(Sw,0.0),1.0); So = std::min(std::max(So,0.0),1.0); Sg = std::min(std::max(Sg,0.0),1.0); Ssol = std::min(std::max(Ssol,0.0),1.0); - + Scalar St = Sw + So + Sg + Ssol; Sw = Sw/St; Sg = Sg/St; @@ -622,8 +621,8 @@ public: if (enableSolvent) (*this)[Indices::solventSaturationIdx] = Ssol; return not(St==1); - - }else if (primaryVarsMeaning() == Sw_po_Rs) { + + }else if (primaryVarsMeaning() == Sw_po_Rs) { Scalar Ssol = 0.0; if (enableSolvent) Ssol = (*this)[Indices::solventSaturationIdx]; @@ -639,7 +638,7 @@ public: if (waterEnabled) (*this)[Indices::waterSaturationIdx]= Sw; if (enableSolvent) - (*this)[Indices::solventSaturationIdx] = Ssol; + (*this)[Indices::solventSaturationIdx] = Ssol; return not(St==1); }else{ assert(primaryVarsMeaning() == Sw_pg_Rv); @@ -660,11 +659,11 @@ public: (*this)[Indices::waterSaturationIdx]= Sw; if (enableSolvent) (*this)[Indices::solventSaturationIdx] = Ssol; - return not(St==1); + return not(St==1); } - assert(false); + assert(false); } - + BlackOilPrimaryVariables& operator=(const BlackOilPrimaryVariables& other) = default; BlackOilPrimaryVariables& operator=(Scalar value) { @@ -749,6 +748,46 @@ private: return (*this)[Indices::temperatureIdx]; } + Scalar microbialConcentration_() const + { + if (!enableMICP) + return 0.0; + + return (*this)[Indices::microbialConcentrationIdx]; + } + + Scalar oxygenConcentration_() const + { + if (!enableMICP) + return 0.0; + + return (*this)[Indices::oxygenConcentrationIdx]; + } + + Scalar ureaConcentration_() const + { + if (!enableMICP) + return 0.0; + + return (*this)[Indices::ureaConcentrationIdx]; + } + + Scalar biofilmConcentration_() const + { + if (!enableMICP) + return 0.0; + + return (*this)[Indices::biofilmConcentrationIdx]; + } + + Scalar calciteConcentration_() const + { + if (!enableMICP) + return 0.0; + + return (*this)[Indices::calciteConcentrationIdx]; + } + template void computeCapillaryPressures_(Container& result, Scalar So, diff --git a/opm/models/blackoil/blackoilproperties.hh b/opm/models/blackoil/blackoilproperties.hh index 3cde076f9..a4a5ecbe0 100644 --- a/opm/models/blackoil/blackoilproperties.hh +++ b/opm/models/blackoil/blackoilproperties.hh @@ -58,6 +58,9 @@ struct EnableFoam { using type = UndefinedProperty; }; //! Enable the ECL-blackoil extension for salt template struct EnableBrine { using type = UndefinedProperty; }; +//! Enable the ECL-blackoil extension for MICP. +template +struct EnableMICP { using type = UndefinedProperty; }; //! Allow the spatial and temporal domains to exhibit non-constant temperature diff --git a/opm/models/blackoil/blackoilratevector.hh b/opm/models/blackoil/blackoilratevector.hh index a037052d0..d302da3e7 100644 --- a/opm/models/blackoil/blackoilratevector.hh +++ b/opm/models/blackoil/blackoilratevector.hh @@ -60,6 +60,7 @@ class BlackOilRateVector using PolymerModule = BlackOilPolymerModule; using FoamModule = BlackOilFoamModule; using BrineModule = BlackOilBrineModule; + using MICPModule = BlackOilMICPModule; enum { numEq = getPropValue() }; enum { numComponents = getPropValue() }; @@ -71,6 +72,7 @@ class BlackOilRateVector enum { enablePolymerMolarWeight = getPropValue() }; enum { enableFoam = getPropValue() }; enum { enableBrine = getPropValue() }; + enum { enableMICP = getPropValue() }; using Toolbox = MathToolbox; using ParentType = Dune::FieldVector; @@ -144,6 +146,10 @@ public: throw std::logic_error("setMolarRate() not implemented for salt water"); } + if ( enableMICP ) { + throw std::logic_error("setMolarRate() not implemented for MICP"); + } + // convert to "surface volume" if requested if (getPropValue()) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { diff --git a/opm/models/blackoil/blackoiltwophaseindices.hh b/opm/models/blackoil/blackoiltwophaseindices.hh index a238e353f..4405161b6 100644 --- a/opm/models/blackoil/blackoiltwophaseindices.hh +++ b/opm/models/blackoil/blackoiltwophaseindices.hh @@ -37,7 +37,7 @@ namespace Opm { * * \brief The primary variable and equation indices for the black-oil model. */ -template +template struct BlackOilTwoPhaseIndices { //! Is phase enabled or not @@ -57,6 +57,9 @@ struct BlackOilTwoPhaseIndices //! Shall energy be conserved? static const bool enableEnergy = numEnergyV > 0; + //! Is MICP involved? + static const bool enableMICP = numMICPsV > 0; + //! Number of solvent components to be considered static const int numSolvents = enableSolvent ? numSolventsV : 0; @@ -78,8 +81,11 @@ struct BlackOilTwoPhaseIndices //! The number of fluid phases static const int numPhases = 2; + //! Number of MICP components to be considered + static const int numMICPs = enableMICP ? numMICPsV : 0; + //! The number of equations - static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine; + static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine + numMICPs; ////////////////////////////// // Primary variable indices @@ -115,17 +121,37 @@ struct BlackOilTwoPhaseIndices static const int polymerMoleWeightIdx = numPolymers > 1 ? polymerConcentrationIdx + 1 : -1000; + //! Index of the primary variable for the first MICP component + static const int microbialConcentrationIdx = + enableMICP ? PVOffset + numPhases + numSolvents : -1000; + + //! Index of the primary variable for the second MICP component + static const int oxygenConcentrationIdx = + numMICPs > 1 ? microbialConcentrationIdx + 1 : -1000; + + //! Index of the primary variable for the third MICP component + static const int ureaConcentrationIdx = + numMICPs > 2 ? oxygenConcentrationIdx + 1 : -1000; + + //! Index of the primary variable for the fourth MICP component + static const int biofilmConcentrationIdx = + numMICPs > 3 ? ureaConcentrationIdx + 1 : -1000; + + //! Index of the primary variable for the fifth MICP component + static const int calciteConcentrationIdx = + numMICPs > 4 ? biofilmConcentrationIdx + 1 : -1000; + //! Index of the primary variable for the foam static const int foamConcentrationIdx = - enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000; + enableFoam ? PVOffset + numPhases + numSolvents + numPolymers + numMICPs : -1000; //! Index of the primary variable for the salt static const int saltConcentrationIdx = - enableBrine ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : -1000; + enableBrine ? PVOffset + numPhases + numSolvents + numPolymers + numMICPs + numFoam : -1000; //! Index of the primary variable for temperature static const int temperatureIdx = - enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : - 1000; + enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine : - 1000; ////////////////////// // Equation indices @@ -188,17 +214,37 @@ struct BlackOilTwoPhaseIndices static const int contiPolymerMWEqIdx = numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000; + //! Index of the continuity equation for the first MICP component + static const int contiMicrobialEqIdx = + enableMICP ? PVOffset + numPhases + numSolvents : -1000; + + //! Index of the continuity equation for the second MICP component + static const int contiOxygenEqIdx = + numMICPs > 1 ? contiMicrobialEqIdx + 1 : -1000; + + //! Index of the continuity equation for the third MICP component + static const int contiUreaEqIdx = + numMICPs > 2 ? contiOxygenEqIdx + 1 : -1000; + + //! Index of the continuity equation for the fourth MICP component + static const int contiBiofilmEqIdx = + numMICPs > 3 ? contiUreaEqIdx + 1 : -1000; + + //! Index of the continuity equation for the fifth MICP component + static const int contiCalciteEqIdx = + numMICPs > 4 ? contiBiofilmEqIdx + 1 : -1000; + //! Index of the continuity equation for the foam component static const int contiFoamEqIdx = - enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000; + enableFoam ? PVOffset + numPhases + numSolvents + numPolymers + numMICPs : -1000; //! Index of the continuity equation for the salt component static const int contiBrineEqIdx = - enableBrine ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : -1000; + enableBrine ? PVOffset + numPhases + numSolvents + numPolymers + numMICPs + numFoam : -1000; //! Index of the continuity equation for energy static const int contiEnergyEqIdx = - enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : -1000; + enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine: -1000; }; } // namespace Opm diff --git a/opm/models/io/vtkblackoilmicpmodule.hh b/opm/models/io/vtkblackoilmicpmodule.hh new file mode 100644 index 000000000..b68f6e2da --- /dev/null +++ b/opm/models/io/vtkblackoilmicpmodule.hh @@ -0,0 +1,264 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::VtkBlackOilMICPModule + */ +#ifndef EWOMS_VTK_BLACK_OIL_MICP_MODULE_HH +#define EWOMS_VTK_BLACK_OIL_MICP_MODULE_HH + +#include + +#include "vtkmultiwriter.hh" +#include "baseoutputmodule.hh" + +#include +#include +#include + +#include + +#include + +namespace Opm::Properties { + +namespace TTag { + +// create new type tag for the VTK multi-phase output +struct VtkBlackOilMICP {}; + +} // namespace TTag + +// create the property tags needed for the MICP output module +template +struct VtkWriteMicrobialConcentration { using type = UndefinedProperty; }; +template +struct VtkWriteOxygenConcentration { using type = UndefinedProperty; }; +template +struct VtkWriteUreaConcentration { using type = UndefinedProperty; }; +template +struct VtkWriteBiofilmConcentration { using type = UndefinedProperty; }; +template +struct VtkWriteCalciteConcentration { using type = UndefinedProperty; }; + +// set default values for what quantities to output +template +struct VtkWriteMicrobialConcentration { static constexpr bool value = true; }; +template +struct VtkWriteOxygenConcentration { static constexpr bool value = true; }; +template +struct VtkWriteUreaConcentration { static constexpr bool value = true; }; +template +struct VtkWriteBiofilmConcentration { static constexpr bool value = true; }; +template +struct VtkWriteCalciteConcentration { static constexpr bool value = true; }; + +} // namespace Opm::Properties + +namespace Opm { +/*! + * \ingroup Vtk + * + * \brief VTK output module for the MICP model's related quantities. + */ +template +class VtkBlackOilMICPModule : public BaseOutputModule +{ + using ParentType = BaseOutputModule; + + using Simulator = GetPropType; + using GridView = GetPropType; + using Scalar = GetPropType; + using Evaluation = GetPropType; + using ElementContext = GetPropType; + + static const int vtkFormat = getPropValue(); + using VtkMultiWriter = ::Opm::VtkMultiWriter; + + enum { enableMICP = getPropValue() }; + + using ScalarBuffer = typename ParentType::ScalarBuffer; + +public: + VtkBlackOilMICPModule(const Simulator& simulator) + : ParentType(simulator) + { } + + /*! + * \brief Register all run-time parameters for the multi-phase VTK output + * module. + */ + static void registerParameters() + { + if (!enableMICP) + return; + + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteMicrobialConcentration, + "Include the concentration of the microbial component in the water phase " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteOxygenConcentration, + "Include the concentration of the oxygen component in the water phase " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteUreaConcentration, + "Include the concentration of the urea component in the water phase " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteBiofilmConcentration, + "Include the biofilm volume fraction " + "in the VTK output files"); + EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteCalciteConcentration, + "Include the calcite volume fraction " + "in the VTK output files"); + } + + /*! + * \brief Allocate memory for the scalar fields we would like to + * write to the VTK file. + */ + void allocBuffers() + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + if (!enableMICP) + return; + + if (microbialConcentrationOutput_()) + this->resizeScalarBuffer_(microbialConcentration_); + if (oxygenConcentrationOutput_()) + this->resizeScalarBuffer_(oxygenConcentration_); + if (ureaConcentrationOutput_()) + this->resizeScalarBuffer_(ureaConcentration_); + if (biofilmConcentrationOutput_()) + this->resizeScalarBuffer_(biofilmConcentration_); + if (calciteConcentrationOutput_()) + this->resizeScalarBuffer_(calciteConcentration_); + } + + /*! + * \brief Modify the internal buffers according to the intensive quantities relevant for + * an element + */ + void processElement(const ElementContext& elemCtx) + { + if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput)) + return; + + if (!enableMICP) + return; + + for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) { + const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0); + unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); + + if (microbialConcentrationOutput_()) + microbialConcentration_[globalDofIdx] = + scalarValue(intQuants.microbialConcentration()); + + if (oxygenConcentrationOutput_()) + oxygenConcentration_[globalDofIdx] = + scalarValue(intQuants.oxygenConcentration()); + + if (ureaConcentrationOutput_()) + ureaConcentration_[globalDofIdx] = + 10 * scalarValue(intQuants.ureaConcentration());//Multypliging by scaling factor 10 (see WellInterface_impl.hpp) + + if (biofilmConcentrationOutput_()) + biofilmConcentration_[globalDofIdx] = + scalarValue(intQuants.biofilmConcentration()); + + if (calciteConcentrationOutput_()) + calciteConcentration_[globalDofIdx] = + scalarValue(intQuants.calciteConcentration()); + + } + } + + /*! + * \brief Add all buffers to the VTK output writer. + */ + void commitBuffers(BaseOutputWriter& baseWriter) + { + VtkMultiWriter *vtkWriter = dynamic_cast(&baseWriter); + if (!vtkWriter) + return; + + if (!enableMICP) + return; + + if (microbialConcentrationOutput_()) + this->commitScalarBuffer_(baseWriter, "microbial concentration", microbialConcentration_); + + if (oxygenConcentrationOutput_()) + this->commitScalarBuffer_(baseWriter, "oxygen concentration", oxygenConcentration_); + + if (ureaConcentrationOutput_()) + this->commitScalarBuffer_(baseWriter, "urea concentration", ureaConcentration_); + + if (biofilmConcentrationOutput_()) + this->commitScalarBuffer_(baseWriter, "biofilm fraction", biofilmConcentration_); + + if (calciteConcentrationOutput_()) + this->commitScalarBuffer_(baseWriter, "calcite fraction", calciteConcentration_); + + } + +private: + static bool microbialConcentrationOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteMicrobialConcentration); + return val; + } + + static bool oxygenConcentrationOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteOxygenConcentration); + return val; + } + + static bool ureaConcentrationOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteUreaConcentration); + return val; + } + + static bool biofilmConcentrationOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteBiofilmConcentration); + return val; + } + + static bool calciteConcentrationOutput_() + { + static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteCalciteConcentration); + return val; + } + + ScalarBuffer microbialConcentration_; + ScalarBuffer oxygenConcentration_; + ScalarBuffer ureaConcentration_; + ScalarBuffer biofilmConcentration_; + ScalarBuffer calciteConcentration_; +}; +} // namespace Opm + +#endif diff --git a/opm/models/io/vtkblackoilmodule.hh b/opm/models/io/vtkblackoilmodule.hh index 3c4bb232a..325a97c78 100644 --- a/opm/models/io/vtkblackoilmodule.hh +++ b/opm/models/io/vtkblackoilmodule.hh @@ -225,8 +225,10 @@ public: const auto& primaryVars = elemCtx.primaryVars(dofIdx, /*timeIdx=*/0); unsigned pvtRegionIdx = elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex(); - Scalar SoMax = std::max(getValue(fs.saturation(oilPhaseIdx)), - elemCtx.problem().maxOilSaturation(globalDofIdx)); + Scalar SoMax = 0.0; + if (FluidSystem::phaseIsActive(oilPhaseIdx)) + SoMax = std::max(getValue(fs.saturation(oilPhaseIdx)), + elemCtx.problem().maxOilSaturation(globalDofIdx)); if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx)) { Scalar x_oG = getValue(fs.moleFraction(oilPhaseIdx, gasCompIdx));