My Project
Loading...
Searching...
No Matches
blackoilmicpmodules.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef OPM_BLACK_OIL_MICP_MODULE_HH
29#define OPM_BLACK_OIL_MICP_MODULE_HH
30
31#include <dune/common/fvector.hh>
32
35
37
38#include <cstddef>
39#include <stdexcept>
40
41namespace Opm {
49{
63
64 using Toolbox = MathToolbox<Evaluation>;
65
66 using TabulatedFunction = typename BlackOilMICPParams<Scalar>::TabulatedFunction;
67
68 enum { waterCompIdx = FluidSystem::waterCompIdx };
69
70 static constexpr unsigned microbialConcentrationIdx = Indices::microbialConcentrationIdx;
71 static constexpr unsigned oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
72 static constexpr unsigned ureaConcentrationIdx = Indices::ureaConcentrationIdx;
73 static constexpr unsigned biofilmConcentrationIdx = Indices::biofilmConcentrationIdx;
74 static constexpr unsigned calciteConcentrationIdx = Indices::calciteConcentrationIdx;
75 static constexpr unsigned contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
76 static constexpr unsigned contiOxygenEqIdx = Indices::contiOxygenEqIdx;
77 static constexpr unsigned contiUreaEqIdx = Indices::contiUreaEqIdx;
78 static constexpr unsigned contiBiofilmEqIdx = Indices::contiBiofilmEqIdx;
79 static constexpr unsigned contiCalciteEqIdx = Indices::contiCalciteEqIdx;
80 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
81
82 static constexpr unsigned enableMICP = enableMICPV;
83
84 static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
85
86public:
89 {
90 params_ = params;
91 }
92
96 static void registerParameters()
97 {
98 if constexpr (enableMICP)
100 }
101
105 static void registerOutputModules(Model& model,
106 Simulator& simulator)
107 {
108 if constexpr (enableMICP)
109 model.addOutputModule(new VtkBlackOilMICPModule<TypeTag>(simulator));
110 }
111
112 static bool eqApplies(unsigned eqIdx)
113 {
114 if constexpr (enableMICP)
115 return eqIdx == contiMicrobialEqIdx || eqIdx == contiOxygenEqIdx || eqIdx == contiUreaEqIdx
116 || eqIdx == contiBiofilmEqIdx || eqIdx == contiCalciteEqIdx;
117 else
118 return false;
119 }
120
121 static Scalar eqWeight([[maybe_unused]] unsigned eqIdx)
122 {
123 assert(eqApplies(eqIdx));
124
125 // TODO: it may be beneficial to chose this differently.
126 return static_cast<Scalar>(1.0);
127 }
128
129 // must be called after water storage is computed
130 template <class LhsEval>
131 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
132 const IntensiveQuantities& intQuants)
133 {
134 if constexpr (enableMICP) {
135 const auto& fs = intQuants.fluidState();
136 LhsEval surfaceVolumeWater = Toolbox::template decay<LhsEval>(fs.invB(waterPhaseIdx))
137 * Toolbox::template decay<LhsEval>(intQuants.porosity());
138 // avoid singular matrix if no water is present
140 // suspended microbes in water phase
141 const LhsEval massMicrobes = surfaceVolumeWater * Toolbox::template decay<LhsEval>(intQuants.microbialConcentration());
143 storage[contiMicrobialEqIdx] += accumulationMicrobes;
144 // oxygen in water phase
145 const LhsEval massOxygen = surfaceVolumeWater * Toolbox::template decay<LhsEval>(intQuants.oxygenConcentration());
147 storage[contiOxygenEqIdx] += accumulationOxygen;
148 // urea in water phase (applying the scaling factor for the urea equation)
149 const LhsEval massUrea = surfaceVolumeWater * Toolbox::template decay<LhsEval>(intQuants.ureaConcentration());
151 storage[contiUreaEqIdx] += accumulationUrea;
153 // biofilm
154 const LhsEval massBiofilm = Toolbox::template decay<LhsEval>(intQuants.biofilmConcentration());
156 storage[contiBiofilmEqIdx] += accumulationBiofilm;
157 // calcite
158 const LhsEval massCalcite = Toolbox::template decay<LhsEval>(intQuants.calciteConcentration());
160 storage[contiCalciteEqIdx] += accumulationCalcite;
161 }
162 }
163
164 template <class UpEval, class Eval, class IntensiveQuantities>
165 static void addMICPFluxes_(RateVector& flux,
166 const Eval& volumeFlux,
167 const IntensiveQuantities& upFs)
168 {
169 if constexpr (enableMICP) {
170 flux[contiMicrobialEqIdx] +=
171 decay<UpEval>(upFs.microbialConcentration())
172 * volumeFlux;
173 flux[contiOxygenEqIdx] +=
174 decay<UpEval>(upFs.oxygenConcentration())
175 * volumeFlux;
176 flux[contiUreaEqIdx] +=
177 decay<UpEval>(upFs.ureaConcentration())
178 * volumeFlux;
179 }
180 }
181
182 // since the urea concentration can be much larger than 1, then we apply a scaling factor
183 static void applyScaling(RateVector& flux)
184 {
185 if constexpr (enableMICP) {
187 }
188 }
189
190 static void computeFlux([[maybe_unused]] RateVector& flux,
191 [[maybe_unused]] const ElementContext& elemCtx,
192 [[maybe_unused]] unsigned scvfIdx,
193 [[maybe_unused]] unsigned timeIdx)
194
195 {
196 if constexpr (enableMICP) {
197 flux[contiMicrobialEqIdx] = 0.0;
198 flux[contiOxygenEqIdx] = 0.0;
199 flux[contiUreaEqIdx] = 0.0;
200 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
201 unsigned focusIdx = elemCtx.focusDofIndex();
202 unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
203 if (upIdx == focusIdx)
205 else
207 }
208 }
209
210 template <class UpstreamEval>
211 static void addMICPFluxes_(RateVector& flux,
212 const ElementContext& elemCtx,
213 unsigned scvfIdx,
214 unsigned timeIdx)
215 {
216 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
217 unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
218 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
219 const auto& fs = up.fluidState();
220 const auto& volFlux = extQuants.volumeFlux(waterPhaseIdx);
222 }
223
224 // see https://doi.org/10.1016/j.ijggc.2021.103256 for the MICP processes in the model
225 static void addSource(RateVector& source,
226 const Problem& problem,
227 const IntensiveQuantities& intQuants,
228 unsigned globalSpaceIdex)
229 {
230 if constexpr (enableMICP) {
231 const auto& velocityInf = problem.model().linearizer().getVelocityInfo();
233 Scalar normVelocityCell = 0.0;
234 for (const auto& velocityInfo : velocityInfos) {
235 normVelocityCell = max(normVelocityCell, std::abs(velocityInfo.velocity[waterPhaseIdx]));
236 }
237
238 // get the model parameters
239 const auto b = intQuants.fluidState().invB(waterPhaseIdx);
240 unsigned satnumIdx = problem.satnumRegionIndex(globalSpaceIdex);
241 Scalar k_a = microbialAttachmentRate(satnumIdx);
242 Scalar k_d = microbialDeathRate(satnumIdx);
243 Scalar rho_b = densityBiofilm(satnumIdx);
244 Scalar rho_c = densityCalcite(satnumIdx);
245 Scalar k_str = detachmentRate(satnumIdx);
246 Scalar eta = detachmentExponent(satnumIdx);
247 Scalar k_o = halfVelocityOxygen(satnumIdx);
248 Scalar k_u = halfVelocityUrea(satnumIdx);
249 Scalar mu = maximumGrowthRate(satnumIdx);
250 Scalar mu_u = maximumUreaUtilization(satnumIdx);
251 Scalar Y_sb = yieldGrowthCoefficient(satnumIdx);
252 Scalar F = oxygenConsumptionFactor(satnumIdx);
253 Scalar Y_uc = yieldUreaToCalciteCoefficient(satnumIdx);
254
255 // compute Monod terms (the negative region is replaced by a straight line)
256 // Schäfer et al (1998) https://doi.org/10.1016/S0169-7722(97)00060-0
257 Evaluation k_g = mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration());
258 Evaluation k_c = mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration());
259 if (intQuants.oxygenConcentration() < 0)
260 k_g = mu * intQuants.oxygenConcentration() / k_o;
261 if (intQuants.ureaConcentration() < 0)
262 k_c = mu_u * intQuants.ureaConcentration() / k_u;
263
264 // compute the processes
265 source[Indices::contiMicrobialEqIdx] += intQuants.microbialConcentration() * intQuants.porosity() *
266 b * (Y_sb * k_g - k_d - k_a) +
267 rho_b * intQuants.biofilmConcentration() * k_str * pow(normVelocityCell, eta);
268
269 source[Indices::contiOxygenEqIdx] -= (intQuants.microbialConcentration() * intQuants.porosity() *
270 b + rho_b * intQuants.biofilmConcentration()) * F * k_g;
271
272 source[Indices::contiUreaEqIdx] -= rho_b * intQuants.biofilmConcentration() * k_c;
273
274 source[Indices::contiBiofilmEqIdx] += intQuants.biofilmConcentration() * (Y_sb * k_g - k_d -
275 k_str * pow(normVelocityCell, eta) - Y_uc * (rho_b / rho_c) *
276 intQuants.biofilmConcentration() * k_c / (intQuants.porosity() +
277 intQuants.biofilmConcentration())) + k_a * intQuants.microbialConcentration() *
278 intQuants.porosity() * b / rho_b;
279
280 source[Indices::contiCalciteEqIdx] += (rho_b / rho_c) * intQuants.biofilmConcentration() * Y_uc * k_c;
281
282 // since the urea concentration can be much larger than 1, then we apply a scaling factor
283 source[Indices::contiUreaEqIdx] *= getPropValue<TypeTag, Properties::BlackOilUreaScalingFactor>();
284 }
285 }
286
287 static void addSource([[maybe_unused]] RateVector& source,
288 [[maybe_unused]] const ElementContext& elemCtx,
289 [[maybe_unused]] unsigned dofIdx,
290 [[maybe_unused]] unsigned timeIdx)
291 {
292 if constexpr (enableMICP) {
293 const auto& problem = elemCtx.problem();
294 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
295 addSource(source, problem, intQuants, dofIdx);
296 }
297 }
298
299 static const Scalar densityBiofilm(unsigned satnumRegionIdx)
300 {
301 return params_.densityBiofilm_[satnumRegionIdx];
302 }
303
304 static const Scalar densityCalcite(unsigned satnumRegionIdx)
305 {
306 return params_.densityCalcite_[satnumRegionIdx];
307 }
308
309 static const Scalar detachmentRate(unsigned satnumRegionIdx)
310 {
311 return params_.detachmentRate_[satnumRegionIdx];
312 }
313
314 static const Scalar detachmentExponent(unsigned satnumRegionIdx)
315 {
316 return params_.detachmentExponent_[satnumRegionIdx];
317 }
318
319 static const Scalar halfVelocityOxygen(unsigned satnumRegionIdx)
320 {
321 return params_.halfVelocityOxygen_[satnumRegionIdx];
322 }
323
324 static const Scalar halfVelocityUrea(unsigned satnumRegionIdx)
325 {
326 return params_.halfVelocityUrea_[satnumRegionIdx];
327 }
328
329 static const Scalar maximumGrowthRate(unsigned satnumRegionIdx)
330 {
331 return params_.maximumGrowthRate_[satnumRegionIdx];
332 }
333
334 static const Scalar maximumUreaUtilization(unsigned satnumRegionIdx)
335 {
336 return params_.maximumUreaUtilization_[satnumRegionIdx];
337 }
338
339 static const Scalar microbialAttachmentRate(unsigned satnumRegionIdx)
340 {
341 return params_.microbialAttachmentRate_[satnumRegionIdx];
342 }
343
344 static const Scalar microbialDeathRate(unsigned satnumRegionIdx)
345 {
346 return params_.microbialDeathRate_[satnumRegionIdx];
347 }
348
349 static const Scalar oxygenConsumptionFactor(unsigned satnumRegionIdx)
350 {
351 return params_.oxygenConsumptionFactor_[satnumRegionIdx];
352 }
353
354 static const Scalar yieldGrowthCoefficient(unsigned satnumRegionIdx)
355 {
356 return params_.yieldGrowthCoefficient_[satnumRegionIdx];
357 }
358
359 static const Scalar yieldUreaToCalciteCoefficient(unsigned satnumRegionIdx)
360 {
361 return params_.yieldUreaToCalciteCoefficient_[satnumRegionIdx];
362 }
363
364 static const Scalar microbialDiffusion(unsigned pvtRegionIdx)
365 {
366 return params_.microbialDiffusion_[pvtRegionIdx];
367 }
368
369 static const Scalar oxygenDiffusion(unsigned pvtRegionIdx)
370 {
371 return params_.oxygenDiffusion_[pvtRegionIdx];
372 }
373
374 static const Scalar ureaDiffusion(unsigned pvtRegionIdx)
375 {
376 return params_.ureaDiffusion_[pvtRegionIdx];
377 }
378
379 static const TabulatedFunction& permfactTable(const ElementContext& elemCtx,
380 unsigned scvIdx,
381 unsigned timeIdx)
382 {
383 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
384 return params_.permfactTable_[satnumRegionIdx];
385 }
386
387 static const TabulatedFunction& permfactTable(unsigned satnumRegionIdx)
388 {
389 return params_.permfactTable_[satnumRegionIdx];
390 }
391
392private:
393 static BlackOilMICPParams<Scalar> params_;
394};
395
396
397template <class TypeTag, bool enableMICPV>
399BlackOilMICPModule<TypeTag, enableMICPV>::params_;
400
410{
412
419
421
422 static constexpr int microbialConcentrationIdx = Indices::microbialConcentrationIdx;
423 static constexpr int oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
424 static constexpr int ureaConcentrationIdx = Indices::ureaConcentrationIdx;
425 static constexpr int biofilmConcentrationIdx = Indices::biofilmConcentrationIdx;
426 static constexpr int calciteConcentrationIdx = Indices::calciteConcentrationIdx;
427 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
428
429
430public:
431
437 void MICPPropertiesUpdate_(const ElementContext& elemCtx,
438 unsigned dofIdx,
439 unsigned timeIdx)
440 {
441 const auto linearizationType = elemCtx.linearizationType();
442 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
443 const Scalar referencePorosity_ = elemCtx.problem().referencePorosity(dofIdx, timeIdx);
444
445 microbialConcentration_ = priVars.makeEvaluation(microbialConcentrationIdx, timeIdx, linearizationType);
446 oxygenConcentration_ = priVars.makeEvaluation(oxygenConcentrationIdx, timeIdx, linearizationType);
447 ureaConcentration_ = priVars.makeEvaluation(ureaConcentrationIdx, timeIdx, linearizationType);
448 biofilmConcentration_ = priVars.makeEvaluation(biofilmConcentrationIdx, timeIdx, linearizationType);
449 calciteConcentration_ = priVars.makeEvaluation(calciteConcentrationIdx, timeIdx, linearizationType);
450
451 const Evaluation porosityFactor = min(1.0 - (biofilmConcentration_+calciteConcentration_)/(referencePorosity_+1e-8), 1.0); //phi/phi_0
452 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, dofIdx, timeIdx);
453 const auto& permfactTable = MICPModule::permfactTable(satnumRegionIdx);
454 permFactor_ = permfactTable.eval(porosityFactor, /*extrapolation=*/true);
455
456 biofilmMass_ = referencePorosity_*biofilmConcentration_*MICPModule::densityBiofilm(satnumRegionIdx);
457 calciteMass_ = referencePorosity_*calciteConcentration_*MICPModule::densityCalcite(satnumRegionIdx);
458 }
459
460 const Evaluation& microbialConcentration() const
461 { return microbialConcentration_; }
462
463 const Evaluation& oxygenConcentration() const
464 { return oxygenConcentration_; }
465
466 const Evaluation& ureaConcentration() const
467 { return ureaConcentration_; }
468
469 const Evaluation& biofilmConcentration() const
470 { return biofilmConcentration_; }
471
472 const Evaluation& calciteConcentration() const
473 { return calciteConcentration_; }
474
475 const Evaluation biofilmMass() const
476 { return biofilmMass_; }
477
478 const Evaluation calciteMass() const
479 { return calciteMass_; }
480
481 const Evaluation& permFactor() const
482 { return permFactor_; }
483
484protected:
485 Evaluation microbialConcentration_;
486 Evaluation oxygenConcentration_;
487 Evaluation ureaConcentration_;
488 Evaluation biofilmConcentration_;
489 Evaluation calciteConcentration_;
490 Evaluation biofilmMass_;
491 Evaluation calciteMass_;
492 Evaluation permFactor_;
493
494};
495
496template <class TypeTag>
498{
502
503public:
504 void MICPPropertiesUpdate_(const ElementContext&,
505 unsigned,
506 unsigned)
507 { }
508
509 const Evaluation& microbialConcentration() const
510 { throw std::logic_error("microbialConcentration() called but MICP is disabled"); }
511
512 const Evaluation& oxygenConcentration() const
513 { throw std::logic_error("oxygenConcentration() called but MICP is disabled"); }
514
515 const Evaluation& ureaConcentration() const
516 { throw std::logic_error("ureaConcentration() called but MICP is disabled"); }
517
518 const Evaluation& biofilmConcentration() const
519 { throw std::logic_error("biofilmConcentration() called but MICP is disabled"); }
520
521 const Evaluation& calciteConcentration() const
522 { throw std::logic_error("calciteConcentration() called but MICP is disabled"); }
523
524 const Evaluation& biofilmMass() const
525 { throw std::logic_error("biofilmMass() called but MICP is disabled"); }
526
527 const Evaluation& calciteMass() const
528 { throw std::logic_error("calciteMass() called but MICP is disabled"); }
529
530 const Evaluation& permFactor() const
531 { throw std::logic_error("permFactor() called but MICP is disabled"); }
532};
533
546
547template <class TypeTag>
549
550} // namespace Opm
551
552#endif
Contains the parameters required to extend the black-oil model by MICP.
Declares the properties required by the black oil model.
Provides the MICP specific extensive quantities to the generic black-oil module's extensive quantitie...
Definition blackoilmicpmodules.hh:543
Provides the volumetric quantities required for the equations needed by the MICP extension of the bla...
Definition blackoilmicpmodules.hh:410
void MICPPropertiesUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the intensive properties needed to handle MICP from the primary variables.
Definition blackoilmicpmodules.hh:437
Contains the high level supplements required to extend the black oil model by MICP.
Definition blackoilmicpmodules.hh:49
static void registerParameters()
Register all run-time parameters for the black-oil MICP module.
Definition blackoilmicpmodules.hh:96
static void registerOutputModules(Model &model, Simulator &simulator)
Register all MICP specific VTK and ECL output modules.
Definition blackoilmicpmodules.hh:105
static void setParams(BlackOilMICPParams< Scalar > &&params)
Set parameters.
Definition blackoilmicpmodules.hh:88
VTK output module for the MICP model's related quantities.
Definition vtkblackoilmicpmodule.hpp:53
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
Struct holding the parameters for the BlackOilMICPModule class.
Definition blackoilmicpparams.hpp:44
VTK output module for the MICP model's related quantities.