My Project
Loading...
Searching...
No Matches
blackoilfoammodules.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 EWOMS_BLACK_OIL_FOAM_MODULE_HH
29#define EWOMS_BLACK_OIL_FOAM_MODULE_HH
30
31#include <dune/common/fvector.hh>
32
33#include <opm/common/OpmLog/OpmLog.hpp>
34
35#include <opm/input/eclipse/EclipseState/Phase.hpp>
36
39
42
43#if HAVE_ECL_INPUT
44#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
45#include <opm/input/eclipse/EclipseState/Tables/FoamadsTable.hpp>
46#include <opm/input/eclipse/EclipseState/Tables/FoammobTable.hpp>
47#endif
48
49#include <string>
50
51namespace Opm {
52
60{
73
74 using Toolbox = MathToolbox<Evaluation>;
75
76 using TabulatedFunction = typename BlackOilFoamParams<Scalar>::TabulatedFunction;
77
78 static constexpr unsigned foamConcentrationIdx = Indices::foamConcentrationIdx;
79 static constexpr unsigned contiFoamEqIdx = Indices::contiFoamEqIdx;
80 static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
81 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
82
83 static constexpr unsigned enableFoam = enableFoamV;
84
85 static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
86 static constexpr unsigned numPhases = FluidSystem::numPhases;
87
88 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
89
90public:
93 {
94 params_ = params;
95 }
96
100 static void registerParameters()
101 {
102 }
103
107 static void registerOutputModules(Model&,
108 Simulator&)
109 {
110 if constexpr (enableFoam) {
111 if (Parameters::Get<Parameters::EnableVtkOutput>()) {
112 OpmLog::warning("VTK output requested, currently unsupported by the foam module.");
113 }
114 }
115 //model.addOutputModule(new VtkBlackOilFoamModule<TypeTag>(simulator));
116 }
117
118 static bool primaryVarApplies(unsigned pvIdx)
119 {
120 if constexpr (enableFoam)
121 return pvIdx == foamConcentrationIdx;
122 else
123 return false;
124 }
125
126 static std::string primaryVarName([[maybe_unused]] unsigned pvIdx)
127 {
128 assert(primaryVarApplies(pvIdx));
129 return "foam_concentration";
130 }
131
132 static Scalar primaryVarWeight([[maybe_unused]] unsigned pvIdx)
133 {
134 assert(primaryVarApplies(pvIdx));
135
136 // TODO: it may be beneficial to chose this differently.
137 return static_cast<Scalar>(1.0);
138 }
139
140 static bool eqApplies(unsigned eqIdx)
141 {
142 if constexpr (enableFoam)
143 return eqIdx == contiFoamEqIdx;
144 else
145 return false;
146
147 }
148
149 static std::string eqName([[maybe_unused]] unsigned eqIdx)
150 {
151 assert(eqApplies(eqIdx));
152
153 return "conti^foam";
154 }
155
156 static Scalar eqWeight([[maybe_unused]] unsigned eqIdx)
157 {
158 assert(eqApplies(eqIdx));
159
160 // TODO: it may be beneficial to chose this differently.
161 return static_cast<Scalar>(1.0);
162 }
163
164 // must be called after water storage is computed
165 template <class LhsEval>
166 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
167 const IntensiveQuantities& intQuants)
168 {
169 if constexpr (enableFoam) {
170 const auto& fs = intQuants.fluidState();
171
172 LhsEval surfaceVolume = Toolbox::template decay<LhsEval>(intQuants.porosity());
173 if (params_.transport_phase_ == Phase::WATER) {
174 surfaceVolume *= (Toolbox::template decay<LhsEval>(fs.saturation(waterPhaseIdx))
175 * Toolbox::template decay<LhsEval>(fs.invB(waterPhaseIdx)));
176 } else if (params_.transport_phase_ == Phase::GAS) {
177 surfaceVolume *= (Toolbox::template decay<LhsEval>(fs.saturation(gasPhaseIdx))
178 * Toolbox::template decay<LhsEval>(fs.invB(gasPhaseIdx)));
179 } else if (params_.transport_phase_ == Phase::SOLVENT) {
180 if constexpr (enableSolvent) {
181 surfaceVolume *= (Toolbox::template decay<LhsEval>( intQuants.solventSaturation())
182 * Toolbox::template decay<LhsEval>(intQuants.solventInverseFormationVolumeFactor()));
183 }
184 } else {
185 throw std::runtime_error("Transport phase is GAS/WATER/SOLVENT");
186 }
187
188 // Avoid singular matrix if no gas is present.
190
191 // Foam/surfactant in free phase.
193 * Toolbox::template decay<LhsEval>(intQuants.foamConcentration());
194
195 // Adsorbed foam/surfactant.
196 const LhsEval adsorbedFoam =
197 Toolbox::template decay<LhsEval>(1.0 - intQuants.porosity())
198 * Toolbox::template decay<LhsEval>(intQuants.foamRockDensity())
199 * Toolbox::template decay<LhsEval>(intQuants.foamAdsorbed());
200
202 storage[contiFoamEqIdx] += accumulationFoam;
203 }
204 }
205
206 static void computeFlux([[maybe_unused]] RateVector& flux,
207 [[maybe_unused]] const ElementContext& elemCtx,
208 [[maybe_unused]] unsigned scvfIdx,
209 [[maybe_unused]] unsigned timeIdx)
210
211 {
212 if constexpr (enableFoam) {
213 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
214 const unsigned inIdx = extQuants.interiorIndex();
215
216 // The effect of the mobility reduction factor is
217 // incorporated in the mobility for the relevant phase,
218 // so fluxes do not need modification here.
219 switch (transportPhase()) {
220 case Phase::WATER: {
221 const unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
222 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
223 if (upIdx == inIdx) {
224 flux[contiFoamEqIdx] =
225 extQuants.volumeFlux(waterPhaseIdx)
226 *up.fluidState().invB(waterPhaseIdx)
227 *up.foamConcentration();
228 } else {
229 flux[contiFoamEqIdx] =
230 extQuants.volumeFlux(waterPhaseIdx)
231 *decay<Scalar>(up.fluidState().invB(waterPhaseIdx))
232 *decay<Scalar>(up.foamConcentration());
233 }
234 break;
235 }
236 case Phase::GAS: {
237 const unsigned upIdx = extQuants.upstreamIndex(gasPhaseIdx);
238 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
239 if (upIdx == inIdx) {
240 flux[contiFoamEqIdx] =
241 extQuants.volumeFlux(gasPhaseIdx)
242 *up.fluidState().invB(gasPhaseIdx)
243 *up.foamConcentration();
244 } else {
245 flux[contiFoamEqIdx] =
246 extQuants.volumeFlux(gasPhaseIdx)
247 *decay<Scalar>(up.fluidState().invB(gasPhaseIdx))
248 *decay<Scalar>(up.foamConcentration());
249 }
250 break;
251 }
252 case Phase::SOLVENT: {
253 if constexpr (enableSolvent) {
254 const unsigned upIdx = extQuants.solventUpstreamIndex();
255 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
256 if (upIdx == inIdx) {
257 flux[contiFoamEqIdx] =
258 extQuants.solventVolumeFlux()
259 *up.solventInverseFormationVolumeFactor()
260 *up.foamConcentration();
261 } else {
262 flux[contiFoamEqIdx] =
263 extQuants.solventVolumeFlux()
264 *decay<Scalar>(up.solventInverseFormationVolumeFactor())
265 *decay<Scalar>(up.foamConcentration());
266 }
267 } else {
268 throw std::runtime_error("Foam transport phase is SOLVENT but SOLVENT is not activated.");
269 }
270 break;
271 }
272 default: {
273 throw std::runtime_error("Foam transport phase must be GAS/WATER/SOLVENT.");
274 }
275 }
276 }
277 }
278
282 static Scalar computeUpdateError(const PrimaryVariables&,
283 const EqVector&)
284 {
285 // do not consider the change of foam primary variables for convergence
286 // TODO: maybe this should be changed
287 return static_cast<Scalar>(0.0);
288 }
289
290 template <class DofEntity>
291 static void serializeEntity([[maybe_unused]] const Model& model,
292 [[maybe_unused]] std::ostream& outstream,
293 [[maybe_unused]] const DofEntity& dof)
294 {
295 if constexpr (enableFoam) {
296 unsigned dofIdx = model.dofMapper().index(dof);
297 const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
298 outstream << priVars[foamConcentrationIdx];
299 }
300 }
301
302 template <class DofEntity>
303 static void deserializeEntity([[maybe_unused]] Model& model,
304 [[maybe_unused]] std::istream& instream,
305 [[maybe_unused]] const DofEntity& dof)
306 {
307 if constexpr (enableFoam) {
308 unsigned dofIdx = model.dofMapper().index(dof);
309 PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
310 PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
311
312 instream >> priVars0[foamConcentrationIdx];
313
314 // set the primary variables for the beginning of the current time step.
315 priVars1[foamConcentrationIdx] = priVars0[foamConcentrationIdx];
316 }
317 }
318
319 static const Scalar foamRockDensity(const ElementContext& elemCtx,
320 unsigned scvIdx,
321 unsigned timeIdx)
322 {
323 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
324 return params_.foamRockDensity_[satnumRegionIdx];
325 }
326
327 static bool foamAllowDesorption(const ElementContext& elemCtx,
328 unsigned scvIdx,
329 unsigned timeIdx)
330 {
331 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
332 return params_.foamAllowDesorption_[satnumRegionIdx];
333 }
334
335 static const TabulatedFunction& adsorbedFoamTable(const ElementContext& elemCtx,
336 unsigned scvIdx,
337 unsigned timeIdx)
338 {
339 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
340 return params_.adsorbedFoamTable_[satnumRegionIdx];
341 }
342
343 static const TabulatedFunction& gasMobilityMultiplierTable(const ElementContext& elemCtx,
344 unsigned scvIdx,
345 unsigned timeIdx)
346 {
347 unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
348 return params_.gasMobilityMultiplierTable_[pvtnumRegionIdx];
349 }
350
351 static const typename BlackOilFoamParams<Scalar>::FoamCoefficients&
352 foamCoefficients(const ElementContext& elemCtx,
353 const unsigned scvIdx,
354 const unsigned timeIdx)
355 {
356 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
357 return params_.foamCoefficients_[satnumRegionIdx];
358 }
359
360 static Phase transportPhase() {
361 return params_.transport_phase_;
362 }
363
364private:
365 static BlackOilFoamParams<Scalar> params_;
366};
367
368template <class TypeTag, bool enableFoam>
370BlackOilFoamModule<TypeTag, enableFoam>::params_;
371
381{
383
391
393
395 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
396
397 static constexpr int foamConcentrationIdx = Indices::foamConcentrationIdx;
398 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
399 static constexpr unsigned oilPhaseIdx = FluidSystem::oilPhaseIdx;
400 static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
401
402public:
403
409 void foamPropertiesUpdate_(const ElementContext& elemCtx,
410 unsigned dofIdx,
411 unsigned timeIdx)
412 {
413 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
414 foamConcentration_ = priVars.makeEvaluation(foamConcentrationIdx, timeIdx);
415 const auto& fs = asImp_().fluidState_;
416
417 // Compute gas mobility reduction factor
418 Evaluation mobilityReductionFactor = 1.0;
419 if (false) {
420 // The functional model is used.
421 // TODO: allow this model.
422 // In order to do this we must allow transport to be in the water phase, not just the gas phase.
423 const auto& foamCoefficients = FoamModule::foamCoefficients(elemCtx, dofIdx, timeIdx);
424
425 const Scalar fm_mob = foamCoefficients.fm_mob;
426
427 const Scalar fm_surf = foamCoefficients.fm_surf;
428 const Scalar ep_surf = foamCoefficients.ep_surf;
429
430 const Scalar fm_oil = foamCoefficients.fm_oil;
431 const Scalar fl_oil = foamCoefficients.fl_oil;
432 const Scalar ep_oil = foamCoefficients.ep_oil;
433
434 const Scalar fm_dry = foamCoefficients.fm_dry;
435 const Scalar ep_dry = foamCoefficients.ep_dry;
436
437 const Scalar fm_cap = foamCoefficients.fm_cap;
438 const Scalar ep_cap = foamCoefficients.ep_cap;
439
440 const Evaluation C_surf = foamConcentration_;
441 const Evaluation Ca = 1e10; // TODO: replace with proper capillary number.
442 const Evaluation S_o = fs.saturation(oilPhaseIdx);
443 const Evaluation S_w = fs.saturation(waterPhaseIdx);
444
445 Evaluation F1 = pow(C_surf/fm_surf, ep_surf);
446 Evaluation F2 = pow((fm_oil-S_o)/(fm_oil-fl_oil), ep_oil);
447 Evaluation F3 = pow(fm_cap/Ca, ep_cap);
448 Evaluation F7 = 0.5 + atan(ep_dry*(S_w-fm_dry))/M_PI;
449
450 mobilityReductionFactor = 1./(1. + fm_mob*F1*F2*F3*F7);
451 } else {
452 // The tabular model is used.
453 // Note that the current implementation only includes the effect of foam concentration (FOAMMOB),
454 // and not the optional pressure dependence (FOAMMOBP) or shear dependence (FOAMMOBS).
455 const auto& gasMobilityMultiplier = FoamModule::gasMobilityMultiplierTable(elemCtx, dofIdx, timeIdx);
456 mobilityReductionFactor = gasMobilityMultiplier.eval(foamConcentration_, /* extrapolate = */ true);
457 }
458
459 // adjust mobility
460 switch (FoamModule::transportPhase()) {
461 case Phase::WATER: {
462 asImp_().mobility_[waterPhaseIdx] *= mobilityReductionFactor;
463 break;
464 }
465 case Phase::GAS: {
466 asImp_().mobility_[gasPhaseIdx] *= mobilityReductionFactor;
467 break;
468 }
469 case Phase::SOLVENT: {
470 if constexpr (enableSolvent) {
471 asImp_().solventMobility_ *= mobilityReductionFactor;
472 } else {
473 throw std::runtime_error("Foam transport phase is SOLVENT but SOLVENT is not activated.");
474 }
475 break;
476 }
477 default: {
478 throw std::runtime_error("Foam transport phase must be GAS/WATER/SOLVENT.");
479 }
480 }
481
482 foamRockDensity_ = FoamModule::foamRockDensity(elemCtx, dofIdx, timeIdx);
483
484 const auto& adsorbedFoamTable = FoamModule::adsorbedFoamTable(elemCtx, dofIdx, timeIdx);
485 foamAdsorbed_ = adsorbedFoamTable.eval(foamConcentration_, /*extrapolate=*/true);
486 if (!FoamModule::foamAllowDesorption(elemCtx, dofIdx, timeIdx)) {
487 throw std::runtime_error("Foam module does not support the 'no desorption' option.");
488 }
489 }
490
491 const Evaluation& foamConcentration() const
492 { return foamConcentration_; }
493
494 Scalar foamRockDensity() const
495 { return foamRockDensity_; }
496
497 const Evaluation& foamAdsorbed() const
498 { return foamAdsorbed_; }
499
500protected:
501 Implementation& asImp_()
502 { return *static_cast<Implementation*>(this); }
503
504 Evaluation foamConcentration_;
505 Scalar foamRockDensity_;
506 Evaluation foamAdsorbed_;
507};
508
509template <class TypeTag>
511{
515
516public:
517 void foamPropertiesUpdate_(const ElementContext&,
518 unsigned,
519 unsigned)
520 { }
521
522
523 const Evaluation& foamConcentration() const
524 { throw std::runtime_error("foamConcentration() called but foam is disabled"); }
525
526 Scalar foamRockDensity() const
527 { throw std::runtime_error("foamRockDensity() called but foam is disabled"); }
528
529 Scalar foamAdsorbed() const
530 { throw std::runtime_error("foamAdsorbed() called but foam is disabled"); }
531};
532
533} // namespace Opm
534
535#endif
Contains the parameters to extend the black-oil model to include the effects of foam.
Declares the properties required by the black oil model.
Provides the volumetric quantities required for the equations needed by the polymers extension of the...
Definition blackoilfoammodules.hh:381
void foamPropertiesUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the intensive properties needed to handle polymers from the primary variables.
Definition blackoilfoammodules.hh:409
Contains the high level supplements required to extend the black oil model to include the effects of ...
Definition blackoilfoammodules.hh:60
static void registerOutputModules(Model &, Simulator &)
Register all foam specific VTK and ECL output modules.
Definition blackoilfoammodules.hh:107
static void registerParameters()
Register all run-time parameters for the black-oil foam module.
Definition blackoilfoammodules.hh:100
static void setParams(BlackOilFoamParams< Scalar > &&params)
Set parameters.
Definition blackoilfoammodules.hh:92
static Scalar computeUpdateError(const PrimaryVariables &, const EqVector &)
Return how much a Newton-Raphson update is considered an error.
Definition blackoilfoammodules.hh:282
Declare the properties used by the infrastructure code of the finite volume discretizations.
Declare the properties used by the infrastructure code of the finite volume discretizations.
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 BlackoilFoamModule class.
Definition blackoilfoamparams.hpp:46