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