63 using Toolbox = MathToolbox<Evaluation>;
65 static constexpr unsigned microbialConcentrationIdx = Indices::microbialConcentrationIdx;
66 static constexpr unsigned oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
67 static constexpr unsigned ureaConcentrationIdx = Indices::ureaConcentrationIdx;
68 static constexpr unsigned biofilmConcentrationIdx = Indices::biofilmConcentrationIdx;
69 static constexpr unsigned calciteConcentrationIdx = Indices::calciteConcentrationIdx;
70 static constexpr unsigned contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
71 static constexpr unsigned contiOxygenEqIdx = Indices::contiOxygenEqIdx;
72 static constexpr unsigned contiUreaEqIdx = Indices::contiUreaEqIdx;
73 static constexpr unsigned contiBiofilmEqIdx = Indices::contiBiofilmEqIdx;
74 static constexpr unsigned contiCalciteEqIdx = Indices::contiCalciteEqIdx;
75 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
77 static constexpr unsigned enableMICP = enableMICPV;
97 const PrimaryVariables& priVars = model.solution(1)[dofIdx];
98 if (phi - priVars[biofilmConcentrationIdx] - priVars[calciteConcentrationIdx] < toleranceBeforeClogging())
99 throw std::logic_error(
"Clogging has been (almost) reached in at least one cell\n");
118 Simulator& simulator)
127 static bool eqApplies(
unsigned eqIdx)
133 return eqIdx == contiMicrobialEqIdx || eqIdx == contiOxygenEqIdx || eqIdx == contiUreaEqIdx || eqIdx == contiBiofilmEqIdx || eqIdx == contiCalciteEqIdx;
136 static Scalar eqWeight([[maybe_unused]]
unsigned eqIdx)
138 assert(eqApplies(eqIdx));
141 return static_cast<Scalar
>(1.0);
145 template <
class LhsEval>
146 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
147 const IntensiveQuantities& intQuants)
152 LhsEval surfaceVolumeWater = Toolbox::template decay<LhsEval>(intQuants.porosity());
154 surfaceVolumeWater = max(surfaceVolumeWater, 1e-10);
156 const LhsEval massMicrobes = surfaceVolumeWater * Toolbox::template decay<LhsEval>(intQuants.microbialConcentration());
157 LhsEval accumulationMicrobes = massMicrobes;
158 storage[contiMicrobialEqIdx] += accumulationMicrobes;
160 const LhsEval massOxygen = surfaceVolumeWater * Toolbox::template decay<LhsEval>(intQuants.oxygenConcentration());
161 LhsEval accumulationOxygen = massOxygen;
162 storage[contiOxygenEqIdx] += accumulationOxygen;
164 const LhsEval massUrea = surfaceVolumeWater * Toolbox::template decay<LhsEval>(intQuants.ureaConcentration());
165 LhsEval accumulationUrea = massUrea;
166 storage[contiUreaEqIdx] += accumulationUrea;
168 const LhsEval massBiofilm = Toolbox::template decay<LhsEval>(intQuants.biofilmConcentration());
169 LhsEval accumulationBiofilm = massBiofilm;
170 storage[contiBiofilmEqIdx] += accumulationBiofilm;
172 const LhsEval massCalcite = Toolbox::template decay<LhsEval>(intQuants.calciteConcentration());
173 LhsEval accumulationCalcite = massCalcite;
174 storage[contiCalciteEqIdx] += accumulationCalcite;
177 static void computeFlux(RateVector& flux,
178 const ElementContext& elemCtx,
186 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
188 const unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
189 const unsigned inIdx = extQuants.interiorIndex();
190 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
193 if (upIdx == inIdx) {
194 flux[contiMicrobialEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * up.microbialConcentration();
195 flux[contiOxygenEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * up.oxygenConcentration();
196 flux[contiUreaEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * up.ureaConcentration();
199 flux[contiMicrobialEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * decay<Scalar>(up.microbialConcentration());
200 flux[contiOxygenEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * decay<Scalar>(up.oxygenConcentration());
201 flux[contiUreaEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * decay<Scalar>(up.ureaConcentration());
206 static void addSource(RateVector& source,
207 const ElementContext& elemCtx,
216 const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
217 const auto& K = elemCtx.problem().intrinsicPermeability(elemCtx, dofIdx, 0);
218 size_t numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
220 for (
unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
221 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
222 unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
223 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
224 const Evaluation& mobWater = up.mobility(waterPhaseIdx);
227 Evaluation waterVolumeVelocity = extQuants.volumeFlux(waterPhaseIdx) / (K[0][0] * mobWater);
228 dpW = std::max(dpW, abs(waterVolumeVelocity));
232 Scalar k_a = microbialAttachmentRate();
233 Scalar k_d = microbialDeathRate();
234 Scalar rho_b = densityBiofilm();
235 Scalar rho_c = densityCalcite();
236 Scalar k_str = detachmentRate();
237 Scalar k_o = halfVelocityOxygen();
238 Scalar k_u = halfVelocityUrea() / 10.0;
239 Scalar mu = maximumGrowthRate();
240 Scalar mu_u = maximumUreaUtilization() / 10.0;
241 Scalar Y_sb = yieldGrowthCoefficient();
242 Scalar F = oxygenConsumptionFactor();
243 Scalar Y_uc = 1.67 * 10;
246 source[Indices::contiMicrobialEqIdx] += intQuants.microbialConcentration() * intQuants.porosity() *
247 (Y_sb * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()) - k_d - k_a)
248 + rho_b * intQuants.biofilmConcentration() * k_str * pow(intQuants.porosity() * dpW, 0.58);
250 source[Indices::contiOxygenEqIdx] -= (intQuants.microbialConcentration() * intQuants.porosity() + rho_b * intQuants.biofilmConcentration()) *
251 F * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration());
253 source[Indices::contiUreaEqIdx] -= rho_b * intQuants.biofilmConcentration() * mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration());
255 source[Indices::contiBiofilmEqIdx] += intQuants.biofilmConcentration() * (Y_sb * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()) - k_d
256 - k_str * pow(intQuants.porosity() * dpW, 0.58) - Y_uc * (rho_b / rho_c) * intQuants.biofilmConcentration() * mu_u *
257 (intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration())) / (intQuants.porosity() + intQuants.biofilmConcentration()))
258 + k_a * intQuants.microbialConcentration() * intQuants.porosity() / rho_b;
260 source[Indices::contiCalciteEqIdx] += (rho_b / rho_c) * intQuants.biofilmConcentration() * Y_uc * mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration());
263 static const Scalar densityBiofilm()
265 return params_.densityBiofilm_;
268 static const Scalar densityCalcite()
270 return params_.densityCalcite_;
273 static const Scalar detachmentRate()
275 return params_.detachmentRate_;
278 static const Scalar criticalPorosity()
280 return params_.criticalPorosity_;
283 static const Scalar fittingFactor()
285 return params_.fittingFactor_;
288 static const Scalar halfVelocityOxygen()
290 return params_.halfVelocityOxygen_;
293 static const Scalar halfVelocityUrea()
295 return params_.halfVelocityUrea_;
298 static const Scalar maximumGrowthRate()
300 return params_.maximumGrowthRate_;
303 static const Scalar maximumOxygenConcentration()
305 return params_.maximumOxygenConcentration_;
308 static const Scalar maximumUreaConcentration()
310 return params_.maximumUreaConcentration_ / 10.0;
313 static const Scalar maximumUreaUtilization()
315 return params_.maximumUreaUtilization_;
318 static const Scalar microbialAttachmentRate()
320 return params_.microbialAttachmentRate_;
323 static const Scalar microbialDeathRate()
325 return params_.microbialDeathRate_;
328 static const Scalar minimumPermeability()
330 return params_.minimumPermeability_;
333 static const Scalar oxygenConsumptionFactor()
335 return params_.oxygenConsumptionFactor_;
338 static const Scalar toleranceBeforeClogging()
340 return params_.toleranceBeforeClogging_;
343 static const Scalar yieldGrowthCoefficient()
345 return params_.yieldGrowthCoefficient_;
348 static const std::vector<Scalar> phi()
354 static BlackOilMICPParams<Scalar> params_;
402 const auto linearizationType = elemCtx.linearizationType();
403 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
404 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
405 const auto& K = elemCtx.problem().intrinsicPermeability(elemCtx, dofIdx, timeIdx);
406 Scalar referencePorosity_ = elemCtx.problem().porosity(elemCtx, dofIdx, timeIdx);
407 Scalar eta = MICPModule::fittingFactor();
408 Scalar k_min = MICPModule::minimumPermeability();
409 Scalar phi_crit = MICPModule::criticalPorosity();
411 microbialConcentration_ = priVars.makeEvaluation(microbialConcentrationIdx, timeIdx, linearizationType);
412 oxygenConcentration_ = priVars.makeEvaluation(oxygenConcentrationIdx, timeIdx, linearizationType);
413 ureaConcentration_ = priVars.makeEvaluation(ureaConcentrationIdx, timeIdx, linearizationType);
414 biofilmConcentration_ = priVars.makeEvaluation(biofilmConcentrationIdx, timeIdx, linearizationType);
415 calciteConcentration_ = priVars.makeEvaluation(calciteConcentrationIdx, timeIdx, linearizationType);
418 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]);