3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
porousmediumflow/3p3c/volumevariables.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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program 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 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
25#ifndef DUMUX_3P3C_VOLUME_VARIABLES_HH
26#define DUMUX_3P3C_VOLUME_VARIABLES_HH
27
28#include <dumux/material/constants.hh>
29#include <dumux/material/fluidstates/compositional.hh>
30#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
31#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
32
33#include <dumux/porousmediumflow/volumevariables.hh>
34#include <dumux/porousmediumflow/nonisothermal/volumevariables.hh>
35#include <dumux/material/solidstates/updatesolidvolumefractions.hh>
36#include <dumux/common/optionalscalar.hh>
37
38#include "primaryvariableswitch.hh"
39
40#include <dumux/common/deprecated.hh>
41
42namespace Dumux {
43
44namespace Detail {
45// helper struct and function detecting if the fluid matrix interaction features a adsorptionModel() function
46template <class FluidMatrixInteraction>
47using AdsorptionModelDetector = decltype(std::declval<FluidMatrixInteraction>().adsorptionModel());
48
49template<class FluidMatrixInteraction>
50static constexpr bool hasAdsorptionModel()
51{ return Dune::Std::is_detected<AdsorptionModelDetector, FluidMatrixInteraction>::value; }
52
53}
54
60template <class Traits>
61class ThreePThreeCVolumeVariables
62: public PorousMediumFlowVolumeVariables<Traits>
63, public EnergyVolumeVariables<Traits, ThreePThreeCVolumeVariables<Traits> >
64{
65 using ParentType = PorousMediumFlowVolumeVariables<Traits>;
66 using EnergyVolVars = EnergyVolumeVariables<Traits, ThreePThreeCVolumeVariables<Traits> >;
67
68 using Scalar = typename Traits::PrimaryVariables::value_type;
69 using PermeabilityType = typename Traits::PermeabilityType;
70
71 using FS = typename Traits::FluidSystem;
72 using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FS>;
73 using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, FS>;
74
75 using ModelTraits = typename Traits::ModelTraits;
76 using Idx = typename ModelTraits::Indices;
77 static constexpr int numFluidComps = ParentType::numFluidComponents();
78 enum {
79 wCompIdx = FS::wCompIdx,
80 gCompIdx = FS::gCompIdx,
81 nCompIdx = FS::nCompIdx,
82
83 wPhaseIdx = FS::wPhaseIdx,
84 gPhaseIdx = FS::gPhaseIdx,
85 nPhaseIdx = FS::nPhaseIdx,
86
87 switch1Idx = Idx::switch1Idx,
88 switch2Idx = Idx::switch2Idx,
89 pressureIdx = Idx::pressureIdx
90 };
91
92 // present phases
93 enum {
94 threePhases = Idx::threePhases,
95 wPhaseOnly = Idx::wPhaseOnly,
96 gnPhaseOnly = Idx::gnPhaseOnly,
97 wnPhaseOnly = Idx::wnPhaseOnly,
98 gPhaseOnly = Idx::gPhaseOnly,
99 wgPhaseOnly = Idx::wgPhaseOnly
100 };
101
102 using EffDiffModel = typename Traits::EffectiveDiffusivityModel;
103 using DiffusionCoefficients = typename Traits::DiffusionType::DiffusionCoefficientsContainer;
104
105public:
107 using FluidState = typename Traits::FluidState;
109 using FluidSystem = typename Traits::FluidSystem;
111 using Indices = typename ModelTraits::Indices;
113 using SolidState = typename Traits::SolidState;
115 using SolidSystem = typename Traits::SolidSystem;
117 using PrimaryVariableSwitch = ThreePThreeCPrimaryVariableSwitch;
118
128 template<class ElemSol, class Problem, class Element, class Scv>
129 void update(const ElemSol &elemSol,
130 const Problem &problem,
131 const Element &element,
132 const Scv& scv)
133 {
134 ParentType::update(elemSol, problem, element, scv);
135 const auto& priVars = elemSol[scv.localDofIndex()];
136 const auto phasePresence = priVars.state();
137
138 constexpr bool useConstraintSolver = ModelTraits::useConstraintSolver();
139
140 EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState_, solidState_);
141
142 /* first the saturations */
143 if (phasePresence == threePhases)
144 {
145 sw_ = priVars[switch1Idx];
146 sn_ = priVars[switch2Idx];
147 sg_ = 1. - sw_ - sn_;
148 }
149 else if (phasePresence == wPhaseOnly)
150 {
151 sw_ = 1.;
152 sn_ = 0.;
153 sg_ = 0.;
154 }
155 else if (phasePresence == gnPhaseOnly)
156 {
157 sw_ = 0.;
158 sn_ = priVars[switch2Idx];
159 sg_ = 1. - sn_;
160 }
161 else if (phasePresence == wnPhaseOnly)
162 {
163 sn_ = priVars[switch2Idx];
164 sw_ = 1. - sn_;
165 sg_ = 0.;
166 }
167 else if (phasePresence == gPhaseOnly)
168 {
169 sw_ = 0.;
170 sn_ = 0.;
171 sg_ = 1.;
172 }
173 else if (phasePresence == wgPhaseOnly)
174 {
175 sw_ = priVars[switch1Idx];
176 sn_ = 0.;
177 sg_ = 1. - sw_;
178 }
179 else
180 DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
181
182 fluidState_.setSaturation(wPhaseIdx, sw_);
183 fluidState_.setSaturation(gPhaseIdx, sg_);
184 fluidState_.setSaturation(nPhaseIdx, sn_);
185
186 /* now the pressures */
187 pg_ = priVars[pressureIdx];
188
189 // calculate capillary pressures
190
191 // old material law interface is deprecated: Replace this by
192 // const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
193 // after the release of 3.3, when the deprecated interface is no longer supported
194 const auto fluidMatrixInteraction = Deprecated::makePcKrSw<3>(Scalar{}, problem.spatialParams(), element, scv, elemSol);
195
196 Scalar pcgw = fluidMatrixInteraction.pcgw(sw_, sn_);
197 Scalar pcnw = fluidMatrixInteraction.pcnw(sw_, sn_);
198 Scalar pcgn = fluidMatrixInteraction.pcgn(sw_, sn_);
199
200 const Scalar pcAlpha = fluidMatrixInteraction.pcAlpha(sw_, sn_);
201 const Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
202
203 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
204 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
205
206 fluidState_.setPressure(wPhaseIdx, pw_);
207 fluidState_.setPressure(gPhaseIdx, pg_);
208 fluidState_.setPressure(nPhaseIdx, pn_);
209
210 // calculate and set all fugacity coefficients. this is
211 // possible because we require all phases to be an ideal
212 // mixture, i.e. fugacity coefficients are not supposed to
213 // depend on composition!
214 typename FluidSystem::ParameterCache paramCache;
215 // assert(FluidSystem::isIdealGas(gPhaseIdx));
216 for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++ phaseIdx) {
217 assert(FluidSystem::isIdealMixture(phaseIdx));
218
219 for (int compIdx = 0; compIdx < ModelTraits::numFluidComponents(); ++ compIdx) {
220 Scalar phi = FluidSystem::fugacityCoefficient(fluidState_, paramCache, phaseIdx, compIdx);
221 fluidState_.setFugacityCoefficient(phaseIdx, compIdx, phi);
222 }
223 }
224
225 // now comes the tricky part: calculate phase composition
226 if (phasePresence == threePhases) {
227 // all phases are present, phase compositions are a
228 // result of the the gas <-> liquid equilibrium. This is
229 // the job of the "MiscibleMultiPhaseComposition"
230 // constraint solver ...
231 if (useConstraintSolver) {
232 MiscibleMultiPhaseComposition::solve(fluidState_,
233 paramCache);
234 }
235 // ... or calculated explicitly this way ...
236 // please note that we experienced some problems with un-regularized
237 // partial pressures due to their calculation from fugacity coefficients -
238 // that's why they are regularized below "within physically meaningful bounds"
239 else {
240 Scalar partPressH2O = FluidSystem::fugacityCoefficient(fluidState_,
241 wPhaseIdx,
242 wCompIdx) * pw_;
243 if (partPressH2O > pg_) partPressH2O = pg_;
244 Scalar partPressNAPL = FluidSystem::fugacityCoefficient(fluidState_,
245 nPhaseIdx,
246 nCompIdx) * pn_;
247 if (partPressNAPL > pg_) partPressNAPL = pg_;
248 Scalar partPressAir = pg_ - partPressH2O - partPressNAPL;
249
250 Scalar xgn = partPressNAPL/pg_;
251 Scalar xgw = partPressH2O/pg_;
252 Scalar xgg = partPressAir/pg_;
253
254 // actually, it's nothing else than Henry coefficient
255 Scalar xwn = partPressNAPL
256 / (FluidSystem::fugacityCoefficient(fluidState_,
257 wPhaseIdx,nCompIdx)
258 * pw_);
259 Scalar xwg = partPressAir
260 / (FluidSystem::fugacityCoefficient(fluidState_,
261 wPhaseIdx,gCompIdx)
262 * pw_);
263 Scalar xww = 1.-xwg-xwn;
264
265 Scalar xnn = 1.-2.e-10;
266 Scalar xna = 1.e-10;
267 Scalar xnw = 1.e-10;
268
269 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
270 fluidState_.setMoleFraction(wPhaseIdx, gCompIdx, xwg);
271 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
272 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
273 fluidState_.setMoleFraction(gPhaseIdx, gCompIdx, xgg);
274 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
275 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
276 fluidState_.setMoleFraction(nPhaseIdx, gCompIdx, xna);
277 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
278
279 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
280 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
281 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
282 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
283 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
284 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
285
286 fluidState_.setDensity(wPhaseIdx, rhoW);
287 fluidState_.setDensity(gPhaseIdx, rhoG);
288 fluidState_.setDensity(nPhaseIdx, rhoN);
289 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
290 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
291 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
292 }
293 }
294 else if (phasePresence == wPhaseOnly) {
295 // only the water phase is present, water phase composition is
296 // stored explicitly.
297
298 // extract mole fractions in the water phase
299 Scalar xwg = priVars[switch1Idx];
300 Scalar xwn = priVars[switch2Idx];
301 Scalar xww = 1 - xwg - xwn;
302
303 // write water mole fractions in the fluid state
304 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
305 fluidState_.setMoleFraction(wPhaseIdx, gCompIdx, xwg);
306 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
307
308 // calculate the composition of the remaining phases (as
309 // well as the densities of all phases). this is the job
310 // of the "ComputeFromReferencePhase" constraint solver ...
311 if (useConstraintSolver)
312 {
313 ComputeFromReferencePhase::solve(fluidState_,
314 paramCache,
315 wPhaseIdx);
316 }
317 // ... or calculated explicitly this way ...
318 else {
319 // note that the gas phase is actually not existing!
320 // thus, this is used as phase switch criterion
321 Scalar xgg = xwg * FluidSystem::fugacityCoefficient(fluidState_,
322 wPhaseIdx,gCompIdx)
323 * pw_ / pg_;
324 Scalar xgn = xwn * FluidSystem::fugacityCoefficient(fluidState_,
325 wPhaseIdx,nCompIdx)
326 * pw_ / pg_;
327 Scalar xgw = FluidSystem::fugacityCoefficient(fluidState_,
328 wPhaseIdx,wCompIdx)
329 * pw_ / pg_;
330
331
332 // note that the gas phase is actually not existing!
333 // thus, this is used as phase switch criterion
334 Scalar xnn = xwn * FluidSystem::fugacityCoefficient(fluidState_,
335 wPhaseIdx,nCompIdx)
336 * pw_;
337 Scalar xna = 1.e-10;
338 Scalar xnw = 1.e-10;
339
340 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
341 fluidState_.setMoleFraction(gPhaseIdx, gCompIdx, xgg);
342 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
343 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
344 fluidState_.setMoleFraction(nPhaseIdx, gCompIdx, xna);
345 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
346
347 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
348 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
349 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
350 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
351 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
352 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
353
354 fluidState_.setDensity(wPhaseIdx, rhoW);
355 fluidState_.setDensity(gPhaseIdx, rhoG);
356 fluidState_.setDensity(nPhaseIdx, rhoN);
357 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
358 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
359 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
360 }
361 }
362 else if (phasePresence == gnPhaseOnly) {
363 // only gas and NAPL phases are present
364 // we have all (partly hypothetical) phase pressures
365 // and temperature and the mole fraction of water in
366 // the gas phase
367
368 // we have all (partly hypothetical) phase pressures
369 // and temperature and the mole fraction of water in
370 // the gas phase
371 Scalar partPressNAPL = fluidState_.fugacityCoefficient(nPhaseIdx, nCompIdx)*pn_;
372 if (partPressNAPL > pg_) partPressNAPL = pg_;
373
374 Scalar xgw = priVars[switch1Idx];
375 Scalar xgn = partPressNAPL/pg_;
376 Scalar xgg = 1.-xgw-xgn;
377
378 // write mole fractions in the fluid state
379 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
380 fluidState_.setMoleFraction(gPhaseIdx, gCompIdx, xgg);
381 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
382
383 // calculate the composition of the remaining phases (as
384 // well as the densities of all phases). this is the job
385 // of the "ComputeFromReferencePhase" constraint solver
386 ComputeFromReferencePhase::solve(fluidState_,
387 paramCache,
388 gPhaseIdx);
389 }
390 else if (phasePresence == wnPhaseOnly) {
391 // only water and NAPL phases are present
392 Scalar partPressNAPL = fluidState_.fugacityCoefficient(nPhaseIdx,nCompIdx)*pn_;
393 if (partPressNAPL > pg_) partPressNAPL = pg_;
394 Scalar henryC = fluidState_.fugacityCoefficient(wPhaseIdx,nCompIdx)*pw_;
395
396 Scalar xwg = priVars[switch1Idx];
397 Scalar xwn = partPressNAPL/henryC;
398 Scalar xww = 1.-xwg-xwn;
399
400 // write mole fractions in the fluid state
401 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
402 fluidState_.setMoleFraction(wPhaseIdx, gCompIdx, xwg);
403 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
404
405 // calculate the composition of the remaining phases (as
406 // well as the densities of all phases). this is the job
407 // of the "ComputeFromReferencePhase" constraint solver
408 ComputeFromReferencePhase::solve(fluidState_,
409 paramCache,
410 wPhaseIdx);
411 }
412 else if (phasePresence == gPhaseOnly) {
413 // only the gas phase is present, gas phase composition is
414 // stored explicitly here below.
415
416 const Scalar xgw = priVars[switch1Idx];
417 const Scalar xgn = priVars[switch2Idx];
418 Scalar xgg = 1 - xgw - xgn;
419
420 // write mole fractions in the fluid state
421 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
422 fluidState_.setMoleFraction(gPhaseIdx, gCompIdx, xgg);
423 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
424
425 // calculate the composition of the remaining phases (as
426 // well as the densities of all phases). this is the job
427 // of the "ComputeFromReferencePhase" constraint solver ...
428 if (useConstraintSolver)
429 {
430 ComputeFromReferencePhase::solve(fluidState_,
431 paramCache,
432 gPhaseIdx);
433 }
434 // ... or calculated explicitly this way ...
435 else {
436
437 // note that the water phase is actually not existing!
438 // thus, this is used as phase switch criterion
439 Scalar xww = xgw * pg_
440 / (FluidSystem::fugacityCoefficient(fluidState_,
441 wPhaseIdx,wCompIdx)
442 * pw_);
443 Scalar xwn = 1.e-10;
444 Scalar xwg = 1.e-10;
445
446 // note that the NAPL phase is actually not existing!
447 // thus, this is used as phase switch criterion
448 Scalar xnn = xgn * pg_
449 / (FluidSystem::fugacityCoefficient(fluidState_,
450 nPhaseIdx,nCompIdx)
451 * pn_);
452 Scalar xna = 1.e-10;
453 Scalar xnw = 1.e-10;
454
455 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
456 fluidState_.setMoleFraction(wPhaseIdx, gCompIdx, xwg);
457 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
458 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
459 fluidState_.setMoleFraction(nPhaseIdx, gCompIdx, xna);
460 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
461
462 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
463 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
464 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
465 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
466 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
467 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
468
469 fluidState_.setDensity(wPhaseIdx, rhoW);
470 fluidState_.setDensity(gPhaseIdx, rhoG);
471 fluidState_.setDensity(nPhaseIdx, rhoN);
472 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
473 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
474 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
475 }
476 }
477 else if (phasePresence == wgPhaseOnly) {
478 // only water and gas phases are present
479 Scalar xgn = priVars[switch2Idx];
480 Scalar partPressH2O = fluidState_.fugacityCoefficient(wPhaseIdx, wCompIdx)*pw_;
481 if (partPressH2O > pg_) partPressH2O = pg_;
482
483 Scalar xgw = partPressH2O/pg_;
484 Scalar xgg = 1.-xgn-xgw;
485
486 // write mole fractions in the fluid state
487 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
488 fluidState_.setMoleFraction(gPhaseIdx, gCompIdx, xgg);
489 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
490
491 // calculate the composition of the remaining phases (as
492 // well as the densities of all phases). this is the job
493 // of the "ComputeFromReferencePhase" constraint solver ...
494 if (useConstraintSolver)
495 {
496 ComputeFromReferencePhase::solve(fluidState_,
497 paramCache,
498 gPhaseIdx);
499 }
500 // ... or calculated explicitly this way ...
501 else {
502 // actually, it's nothing else than Henry coefficient
503 Scalar xwn = xgn * pg_
504 / (FluidSystem::fugacityCoefficient(fluidState_,
505 wPhaseIdx,nCompIdx)
506 * pw_);
507 Scalar xwg = xgg * pg_
508 / (FluidSystem::fugacityCoefficient(fluidState_,
509 wPhaseIdx,gCompIdx)
510 * pw_);
511 Scalar xww = 1.-xwg-xwn;
512
513 // note that the NAPL phase is actually not existing!
514 // thus, this is used as phase switch criterion
515 Scalar xnn = xgn * pg_
516 / (FluidSystem::fugacityCoefficient(fluidState_,
517 nPhaseIdx,nCompIdx)
518 * pn_);
519 Scalar xna = 1.e-10;
520 Scalar xnw = 1.e-10;
521
522 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
523 fluidState_.setMoleFraction(wPhaseIdx, gCompIdx, xwg);
524 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
525 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
526 fluidState_.setMoleFraction(nPhaseIdx, gCompIdx, xna);
527 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
528
529 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
530 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
531 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
532 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
533 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
534 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
535
536 fluidState_.setDensity(wPhaseIdx, rhoW);
537 fluidState_.setDensity(gPhaseIdx, rhoG);
538 fluidState_.setDensity(nPhaseIdx, rhoN);
539 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
540 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
541 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
542 }
543 }
544 else
545 DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
546
547 for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
548 {
549 // mobilities
550 const Scalar mu =
551 FluidSystem::viscosity(fluidState_,
552 paramCache,
553 phaseIdx);
554 fluidState_.setViscosity(phaseIdx,mu);
555
556 const Scalar kr = fluidMatrixInteraction.kr(phaseIdx,
557 fluidState_.saturation(wPhaseIdx),
558 fluidState_.saturation(nPhaseIdx));
559 mobility_[phaseIdx] = kr / mu;
560 }
561
562 // material dependent parameters for NAPL adsorption (only if law is provided)
563 if constexpr (Detail::hasAdsorptionModel<std::decay_t<decltype(fluidMatrixInteraction)>>())
564 bulkDensTimesAdsorpCoeff_ = fluidMatrixInteraction.adsorptionModel().bulkDensTimesAdsorpCoeff();
565
566 /* compute the diffusion coefficient
567 * \note This is the part of the diffusion coefficient determined by the fluid state, e.g.
568 * important if they are tabularized. In the diffusive flux computation (e.g. Fick's law)
569 * this gets converted into an effecient coefficient depending on saturation and porosity.
570 * We can then add a normalized tensorial component
571 * e.g. obtained from DTI from the spatial params (currently not implemented)
572 */
573
574 auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
575 {
576 return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
577 };
578
579 // porosity & permeabilty
580 updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
581
582 effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
583
584 EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
585 permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
586 EnergyVolVars::updateEffectiveThermalConductivity();
587
588 // compute and set the enthalpy
589 for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
590 {
591 Scalar h = EnergyVolVars::enthalpy(fluidState_, paramCache, phaseIdx);
592 fluidState_.setEnthalpy(phaseIdx, h);
593 }
594 }
595
599 const FluidState &fluidState() const
600 { return fluidState_; }
601
605 const SolidState &solidState() const
606 { return solidState_; }
607
613 Scalar averageMolarMass(int phaseIdx) const
614 { return fluidState_.averageMolarMass(phaseIdx); }
615
622 Scalar saturation(const int phaseIdx) const
623 { return fluidState_.saturation(phaseIdx); }
624
632 Scalar massFraction(const int phaseIdx, const int compIdx) const
633 { return fluidState_.massFraction(phaseIdx, compIdx); }
634
642 Scalar moleFraction(const int phaseIdx, const int compIdx) const
643 { return fluidState_.moleFraction(phaseIdx, compIdx); }
644
651 Scalar density(const int phaseIdx) const
652 { return fluidState_.density(phaseIdx); }
653
660 Scalar molarDensity(const int phaseIdx) const
661 { return fluidState_.molarDensity(phaseIdx); }
662
669 Scalar pressure(const int phaseIdx) const
670 { return fluidState_.pressure(phaseIdx); }
671
679 Scalar temperature() const
680 { return fluidState_.temperature(/*phaseIdx=*/0); }
681
688 Scalar mobility(const int phaseIdx) const
689 { return mobility_[phaseIdx]; }
690
694 Scalar capillaryPressure() const
695 { return fluidState_.capillaryPressure(); }
696
700 Scalar porosity() const
701 { return solidState_.porosity(); }
702
706 Scalar bulkDensTimesAdsorpCoeff() const
707 {
708 if (bulkDensTimesAdsorpCoeff_)
709 return bulkDensTimesAdsorpCoeff_.value();
710 else
711 DUNE_THROW(Dune::NotImplemented, "Your spatialParams do not provide an adsorption model");
712 }
713
717 const PermeabilityType& permeability() const
718 { return permeability_; }
719
720 /*
721 * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
722 */
723 Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
724 {
725 typename FluidSystem::ParameterCache paramCache;
726 paramCache.updatePhase(fluidState_, phaseIdx);
727 return FluidSystem::diffusionCoefficient(fluidState_, paramCache, phaseIdx, compJIdx);
728 }
729
733 Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
734 { return effectiveDiffCoeff_(phaseIdx, compIIdx, compJIdx); }
735
736protected:
737 FluidState fluidState_;
738 SolidState solidState_;
739
740
741private:
742 Scalar sw_, sg_, sn_, pg_, pw_, pn_;
743
744 Scalar moleFrac_[ModelTraits::numFluidPhases()][ModelTraits::numFluidComponents()];
745 Scalar massFrac_[ModelTraits::numFluidPhases()][ModelTraits::numFluidComponents()];
746
747 PermeabilityType permeability_;
748 Scalar mobility_[ModelTraits::numFluidPhases()];
749 OptionalScalar<Scalar> bulkDensTimesAdsorpCoeff_;
750
751 DiffusionCoefficients effectiveDiffCoeff_;
752};
753
754} // end namespace Dumux
755
756#endif
deprecated.hh
Helpers for deprecation.
optionalscalar.hh
A wrapper that can either contain a valid Scalar or NaN.
computefromreferencephase.hh
Computes all quantities of a generic fluid state if a reference phase has been specified.
misciblemultiphasecomposition.hh
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
compositional.hh
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
updatesolidvolumefractions.hh
Update the solid volume fractions (inert and reacitve) and set them in the solidstate.
constants.hh
A central place for various physical constants occuring in some equations.
Dumux::updateSolidVolumeFractions
void updateSolidVolumeFractions(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv, SolidState &solidState, const int solidVolFracOffset)
update the solid volume fractions (inert and reacitve) and set them in the solidstate
Definition: updatesolidvolumefractions.hh:36
Dumux
Definition: adapt.hh:29
Dumux::Detail::AdsorptionModelDetector
decltype(std::declval< FluidMatrixInteraction >().adsorptionModel()) AdsorptionModelDetector
Definition: porousmediumflow/3p3c/volumevariables.hh:47
Dumux::Detail::hasAdsorptionModel
static constexpr bool hasAdsorptionModel()
Definition: porousmediumflow/3p3c/volumevariables.hh:50
Dumux::IOName::phasePresence
std::string phasePresence() noexcept
I/O name of phase presence.
Definition: name.hh:147
Dumux::IOName::viscosity
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
Dumux::IOName::molarDensity
std::string molarDensity(int phaseIdx) noexcept
I/O name of molar density for multiphase systems.
Definition: name.hh:83
Dumux::IOName::density
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Dumux::OptionalScalar::value
T value() const
Definition: optionalscalar.hh:48
Dumux::ComputeFromReferencePhase
Computes all quantities of a generic fluid state if a reference phase has been specified.
Definition: computefromreferencephase.hh:76
Dumux::ComputeFromReferencePhase::solve
static void solve(FluidState &fluidState, ParameterCache &paramCache, int refPhaseIdx)
Computes all quantities of a generic fluid state if a reference phase has been specified.
Definition: computefromreferencephase.hh:111
Dumux::MiscibleMultiPhaseComposition
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: misciblemultiphasecomposition.hh:59
Dumux::MiscibleMultiPhaseComposition::solve
static void solve(FluidState &fluidState, ParameterCache &paramCache, int knownPhaseIdx=0)
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: misciblemultiphasecomposition.hh:81
Dumux::ThreePThreeCPrimaryVariableSwitch
The primary variable switch controlling the phase presence state variable.
Definition: 3p3c/primaryvariableswitch.hh:40
Dumux::ThreePThreeCVolumeVariables
Contains the quantities which are are constant within a finite volume in the three-phase three-compon...
Definition: porousmediumflow/3p3c/volumevariables.hh:64
Dumux::ThreePThreeCVolumeVariables::massFraction
Scalar massFraction(const int phaseIdx, const int compIdx) const
Returns the mass fraction of a given component in a given phase within the control volume in .
Definition: porousmediumflow/3p3c/volumevariables.hh:632
Dumux::ThreePThreeCVolumeVariables::Indices
typename ModelTraits::Indices Indices
export the indices
Definition: porousmediumflow/3p3c/volumevariables.hh:111
Dumux::ThreePThreeCVolumeVariables::averageMolarMass
Scalar averageMolarMass(int phaseIdx) const
Returns the average molar mass of the fluid phase.
Definition: porousmediumflow/3p3c/volumevariables.hh:613
Dumux::ThreePThreeCVolumeVariables::fluidState_
FluidState fluidState_
Definition: porousmediumflow/3p3c/volumevariables.hh:737
Dumux::ThreePThreeCVolumeVariables::update
void update(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv)
Update all quantities for a given control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:129
Dumux::ThreePThreeCVolumeVariables::solidState_
SolidState solidState_
Definition: porousmediumflow/3p3c/volumevariables.hh:738
Dumux::ThreePThreeCVolumeVariables::saturation
Scalar saturation(const int phaseIdx) const
Returns the effective saturation of a given phase within the control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:622
Dumux::ThreePThreeCVolumeVariables::moleFraction
Scalar moleFraction(const int phaseIdx, const int compIdx) const
Returns the mole fraction of a given component in a given phase within the control volume in .
Definition: porousmediumflow/3p3c/volumevariables.hh:642
Dumux::ThreePThreeCVolumeVariables::molarDensity
Scalar molarDensity(const int phaseIdx) const
Returns the molar density of a given phase within the control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:660
Dumux::ThreePThreeCVolumeVariables::SolidState
typename Traits::SolidState SolidState
export type of solid state
Definition: porousmediumflow/3p3c/volumevariables.hh:113
Dumux::ThreePThreeCVolumeVariables::pressure
Scalar pressure(const int phaseIdx) const
Returns the effective pressure of a given phase within the control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:669
Dumux::ThreePThreeCVolumeVariables::solidState
const SolidState & solidState() const
Returns the phase state for the control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:605
Dumux::ThreePThreeCVolumeVariables::permeability
const PermeabilityType & permeability() const
Returns the average permeability within the control volume in .
Definition: porousmediumflow/3p3c/volumevariables.hh:717
Dumux::ThreePThreeCVolumeVariables::FluidState
typename Traits::FluidState FluidState
export fluid state type
Definition: porousmediumflow/3p3c/volumevariables.hh:107
Dumux::ThreePThreeCVolumeVariables::temperature
Scalar temperature() const
Returns temperature inside the sub-control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:679
Dumux::ThreePThreeCVolumeVariables::capillaryPressure
Scalar capillaryPressure() const
Returns the effective capillary pressure within the control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:694
Dumux::ThreePThreeCVolumeVariables::mobility
Scalar mobility(const int phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:688
Dumux::ThreePThreeCVolumeVariables::fluidState
const FluidState & fluidState() const
Returns the phase state for the control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:599
Dumux::ThreePThreeCVolumeVariables::density
Scalar density(const int phaseIdx) const
Returns the mass density of a given phase within the control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:651
Dumux::ThreePThreeCVolumeVariables::FluidSystem
typename Traits::FluidSystem FluidSystem
export fluid system type
Definition: porousmediumflow/3p3c/volumevariables.hh:109
Dumux::ThreePThreeCVolumeVariables::effectiveDiffusionCoefficient
Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
Returns the effective diffusion coefficients for a phase in .
Definition: porousmediumflow/3p3c/volumevariables.hh:733
Dumux::ThreePThreeCVolumeVariables::diffusionCoefficient
Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
Definition: porousmediumflow/3p3c/volumevariables.hh:723
Dumux::ThreePThreeCVolumeVariables::bulkDensTimesAdsorpCoeff
Scalar bulkDensTimesAdsorpCoeff() const
Returns the adsorption information.
Definition: porousmediumflow/3p3c/volumevariables.hh:706
Dumux::ThreePThreeCVolumeVariables::porosity
Scalar porosity() const
Returns the average porosity within the control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:700
Dumux::ThreePThreeCVolumeVariables::SolidSystem
typename Traits::SolidSystem SolidSystem
export type of solid system
Definition: porousmediumflow/3p3c/volumevariables.hh:115
Dumux::EnergyVolumeVariablesImplementation
Definition: porousmediumflow/nonisothermal/volumevariables.hh:75
Dumux::PorousMediumFlowVolumeVariables
The isothermal base class.
Definition: porousmediumflow/volumevariables.hh:40
Dumux::PorousMediumFlowVolumeVariables::numFluidComponents
static constexpr int numFluidComponents()
Return number of components considered by the model.
Definition: porousmediumflow/volumevariables.hh:52
Dumux::PorousMediumFlowVolumeVariables::priVars
const PrimaryVariables & priVars() const
Returns the vector of primary variables.
Definition: porousmediumflow/volumevariables.hh:76
Dumux::PorousMediumFlowVolumeVariables::update
void update(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv)
Updates all quantities for a given control volume.
Definition: porousmediumflow/volumevariables.hh:64
volumevariables.hh
Base class for the model specific class which provides access to all volume averaged quantities.
volumevariables.hh
Base class for the model specific class which provides access to all volume averaged quantities.
primaryvariableswitch.hh
The primary variable switch for the extended Richards model.
Generated by doxygen 1.9.3