3.2-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
37#include "primaryvariableswitch.hh"
38
39namespace Dumux {
40
46template <class Traits>
47class ThreePThreeCVolumeVariables
48: public PorousMediumFlowVolumeVariables<Traits>
49, public EnergyVolumeVariables<Traits, ThreePThreeCVolumeVariables<Traits> >
50{
51 using ParentType = PorousMediumFlowVolumeVariables<Traits>;
52 using EnergyVolVars = EnergyVolumeVariables<Traits, ThreePThreeCVolumeVariables<Traits> >;
53
54 using Scalar = typename Traits::PrimaryVariables::value_type;
55 using PermeabilityType = typename Traits::PermeabilityType;
56
57 using FS = typename Traits::FluidSystem;
58 using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FS>;
59 using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, FS>;
60
61 using ModelTraits = typename Traits::ModelTraits;
62 using Idx = typename ModelTraits::Indices;
63 static constexpr int numFluidComps = ParentType::numFluidComponents();
64 enum {
65 wCompIdx = FS::wCompIdx,
66 gCompIdx = FS::gCompIdx,
67 nCompIdx = FS::nCompIdx,
68
69 wPhaseIdx = FS::wPhaseIdx,
70 gPhaseIdx = FS::gPhaseIdx,
71 nPhaseIdx = FS::nPhaseIdx,
72
73 switch1Idx = Idx::switch1Idx,
74 switch2Idx = Idx::switch2Idx,
75 pressureIdx = Idx::pressureIdx
76 };
77
78 // present phases
79 enum {
80 threePhases = Idx::threePhases,
81 wPhaseOnly = Idx::wPhaseOnly,
82 gnPhaseOnly = Idx::gnPhaseOnly,
83 wnPhaseOnly = Idx::wnPhaseOnly,
84 gPhaseOnly = Idx::gPhaseOnly,
85 wgPhaseOnly = Idx::wgPhaseOnly
86 };
87
88 using EffDiffModel = typename Traits::EffectiveDiffusivityModel;
89 using DiffusionCoefficients = typename Traits::DiffusionType::DiffusionCoefficientsContainer;
90
91public:
93 using FluidState = typename Traits::FluidState;
95 using FluidSystem = typename Traits::FluidSystem;
97 using Indices = typename ModelTraits::Indices;
99 using SolidState = typename Traits::SolidState;
101 using SolidSystem = typename Traits::SolidSystem;
103 using PrimaryVariableSwitch = ThreePThreeCPrimaryVariableSwitch;
104
114 template<class ElemSol, class Problem, class Element, class Scv>
115 void update(const ElemSol &elemSol,
116 const Problem &problem,
117 const Element &element,
118 const Scv& scv)
119 {
120 ParentType::update(elemSol, problem, element, scv);
121 const auto& priVars = elemSol[scv.localDofIndex()];
122 const auto phasePresence = priVars.state();
123
124 constexpr bool useConstraintSolver = ModelTraits::useConstraintSolver();
125
126 // capillary pressure parameters
127 const auto &materialParams =
128 problem.spatialParams().materialLawParams(element, scv, elemSol);
129
130
131 EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState_, solidState_);
132
133 /* first the saturations */
134 if (phasePresence == threePhases)
135 {
136 sw_ = priVars[switch1Idx];
137 sn_ = priVars[switch2Idx];
138 sg_ = 1. - sw_ - sn_;
139 }
140 else if (phasePresence == wPhaseOnly)
141 {
142 sw_ = 1.;
143 sn_ = 0.;
144 sg_ = 0.;
145 }
146 else if (phasePresence == gnPhaseOnly)
147 {
148 sw_ = 0.;
149 sn_ = priVars[switch2Idx];
150 sg_ = 1. - sn_;
151 }
152 else if (phasePresence == wnPhaseOnly)
153 {
154 sn_ = priVars[switch2Idx];
155 sw_ = 1. - sn_;
156 sg_ = 0.;
157 }
158 else if (phasePresence == gPhaseOnly)
159 {
160 sw_ = 0.;
161 sn_ = 0.;
162 sg_ = 1.;
163 }
164 else if (phasePresence == wgPhaseOnly)
165 {
166 sw_ = priVars[switch1Idx];
167 sn_ = 0.;
168 sg_ = 1. - sw_;
169 }
170 else
171 DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
172 Valgrind::CheckDefined(sg_);
173
174 fluidState_.setSaturation(wPhaseIdx, sw_);
175 fluidState_.setSaturation(gPhaseIdx, sg_);
176 fluidState_.setSaturation(nPhaseIdx, sn_);
177
178 /* now the pressures */
179 pg_ = priVars[pressureIdx];
180
181 // calculate capillary pressures
182 using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
183 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
184 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
185 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
186
187 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
188 Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
189
190 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
191 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
192
193 fluidState_.setPressure(wPhaseIdx, pw_);
194 fluidState_.setPressure(gPhaseIdx, pg_);
195 fluidState_.setPressure(nPhaseIdx, pn_);
196
197 // calculate and set all fugacity coefficients. this is
198 // possible because we require all phases to be an ideal
199 // mixture, i.e. fugacity coefficients are not supposed to
200 // depend on composition!
201 typename FluidSystem::ParameterCache paramCache;
202 // assert(FluidSystem::isIdealGas(gPhaseIdx));
203 for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++ phaseIdx) {
204 assert(FluidSystem::isIdealMixture(phaseIdx));
205
206 for (int compIdx = 0; compIdx < ModelTraits::numFluidComponents(); ++ compIdx) {
207 Scalar phi = FluidSystem::fugacityCoefficient(fluidState_, paramCache, phaseIdx, compIdx);
208 fluidState_.setFugacityCoefficient(phaseIdx, compIdx, phi);
209 }
210 }
211
212 // now comes the tricky part: calculate phase composition
213 if (phasePresence == threePhases) {
214 // all phases are present, phase compositions are a
215 // result of the the gas <-> liquid equilibrium. This is
216 // the job of the "MiscibleMultiPhaseComposition"
217 // constraint solver ...
218 if (useConstraintSolver) {
219 MiscibleMultiPhaseComposition::solve(fluidState_,
220 paramCache);
221 }
222 // ... or calculated explicitly this way ...
223 // please note that we experienced some problems with un-regularized
224 // partial pressures due to their calculation from fugacity coefficients -
225 // that's why they are regularized below "within physically meaningful bounds"
226 else {
227 Scalar partPressH2O = FluidSystem::fugacityCoefficient(fluidState_,
228 wPhaseIdx,
229 wCompIdx) * pw_;
230 if (partPressH2O > pg_) partPressH2O = pg_;
231 Scalar partPressNAPL = FluidSystem::fugacityCoefficient(fluidState_,
232 nPhaseIdx,
233 nCompIdx) * pn_;
234 if (partPressNAPL > pg_) partPressNAPL = pg_;
235 Scalar partPressAir = pg_ - partPressH2O - partPressNAPL;
236
237 Scalar xgn = partPressNAPL/pg_;
238 Scalar xgw = partPressH2O/pg_;
239 Scalar xgg = partPressAir/pg_;
240
241 // actually, it's nothing else than Henry coefficient
242 Scalar xwn = partPressNAPL
243 / (FluidSystem::fugacityCoefficient(fluidState_,
244 wPhaseIdx,nCompIdx)
245 * pw_);
246 Scalar xwg = partPressAir
247 / (FluidSystem::fugacityCoefficient(fluidState_,
248 wPhaseIdx,gCompIdx)
249 * pw_);
250 Scalar xww = 1.-xwg-xwn;
251
252 Scalar xnn = 1.-2.e-10;
253 Scalar xna = 1.e-10;
254 Scalar xnw = 1.e-10;
255
256 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
257 fluidState_.setMoleFraction(wPhaseIdx, gCompIdx, xwg);
258 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
259 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
260 fluidState_.setMoleFraction(gPhaseIdx, gCompIdx, xgg);
261 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
262 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
263 fluidState_.setMoleFraction(nPhaseIdx, gCompIdx, xna);
264 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
265
266 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
267 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
268 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
269 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
270 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
271 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
272
273 fluidState_.setDensity(wPhaseIdx, rhoW);
274 fluidState_.setDensity(gPhaseIdx, rhoG);
275 fluidState_.setDensity(nPhaseIdx, rhoN);
276 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
277 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
278 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
279 }
280 }
281 else if (phasePresence == wPhaseOnly) {
282 // only the water phase is present, water phase composition is
283 // stored explicitly.
284
285 // extract mole fractions in the water phase
286 Scalar xwg = priVars[switch1Idx];
287 Scalar xwn = priVars[switch2Idx];
288 Scalar xww = 1 - xwg - xwn;
289
290 // write water mole fractions in the fluid state
291 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
292 fluidState_.setMoleFraction(wPhaseIdx, gCompIdx, xwg);
293 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
294
295 // calculate the composition of the remaining phases (as
296 // well as the densities of all phases). this is the job
297 // of the "ComputeFromReferencePhase" constraint solver ...
298 if (useConstraintSolver)
299 {
300 ComputeFromReferencePhase::solve(fluidState_,
301 paramCache,
302 wPhaseIdx);
303 }
304 // ... or calculated explicitly this way ...
305 else {
306 // note that the gas phase is actually not existing!
307 // thus, this is used as phase switch criterion
308 Scalar xgg = xwg * FluidSystem::fugacityCoefficient(fluidState_,
309 wPhaseIdx,gCompIdx)
310 * pw_ / pg_;
311 Scalar xgn = xwn * FluidSystem::fugacityCoefficient(fluidState_,
312 wPhaseIdx,nCompIdx)
313 * pw_ / pg_;
314 Scalar xgw = FluidSystem::fugacityCoefficient(fluidState_,
315 wPhaseIdx,wCompIdx)
316 * pw_ / pg_;
317
318
319 // note that the gas phase is actually not existing!
320 // thus, this is used as phase switch criterion
321 Scalar xnn = xwn * FluidSystem::fugacityCoefficient(fluidState_,
322 wPhaseIdx,nCompIdx)
323 * pw_;
324 Scalar xna = 1.e-10;
325 Scalar xnw = 1.e-10;
326
327 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
328 fluidState_.setMoleFraction(gPhaseIdx, gCompIdx, xgg);
329 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
330 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
331 fluidState_.setMoleFraction(nPhaseIdx, gCompIdx, xna);
332 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
333
334 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
335 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
336 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
337 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
338 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
339 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
340
341 fluidState_.setDensity(wPhaseIdx, rhoW);
342 fluidState_.setDensity(gPhaseIdx, rhoG);
343 fluidState_.setDensity(nPhaseIdx, rhoN);
344 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
345 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
346 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
347 }
348 }
349 else if (phasePresence == gnPhaseOnly) {
350 // only gas and NAPL phases are present
351 // we have all (partly hypothetical) phase pressures
352 // and temperature and the mole fraction of water in
353 // the gas phase
354
355 // we have all (partly hypothetical) phase pressures
356 // and temperature and the mole fraction of water in
357 // the gas phase
358 Scalar partPressNAPL = fluidState_.fugacityCoefficient(nPhaseIdx, nCompIdx)*pn_;
359 if (partPressNAPL > pg_) partPressNAPL = pg_;
360
361 Scalar xgw = priVars[switch1Idx];
362 Scalar xgn = partPressNAPL/pg_;
363 Scalar xgg = 1.-xgw-xgn;
364
365 // write mole fractions in the fluid state
366 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
367 fluidState_.setMoleFraction(gPhaseIdx, gCompIdx, xgg);
368 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
369
370 // calculate the composition of the remaining phases (as
371 // well as the densities of all phases). this is the job
372 // of the "ComputeFromReferencePhase" constraint solver
373 ComputeFromReferencePhase::solve(fluidState_,
374 paramCache,
375 gPhaseIdx);
376 }
377 else if (phasePresence == wnPhaseOnly) {
378 // only water and NAPL phases are present
379 Scalar partPressNAPL = fluidState_.fugacityCoefficient(nPhaseIdx,nCompIdx)*pn_;
380 if (partPressNAPL > pg_) partPressNAPL = pg_;
381 Scalar henryC = fluidState_.fugacityCoefficient(wPhaseIdx,nCompIdx)*pw_;
382
383 Scalar xwg = priVars[switch1Idx];
384 Scalar xwn = partPressNAPL/henryC;
385 Scalar xww = 1.-xwg-xwn;
386
387 // write mole fractions in the fluid state
388 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
389 fluidState_.setMoleFraction(wPhaseIdx, gCompIdx, xwg);
390 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
391
392 // calculate the composition of the remaining phases (as
393 // well as the densities of all phases). this is the job
394 // of the "ComputeFromReferencePhase" constraint solver
395 ComputeFromReferencePhase::solve(fluidState_,
396 paramCache,
397 wPhaseIdx);
398 }
399 else if (phasePresence == gPhaseOnly) {
400 // only the gas phase is present, gas phase composition is
401 // stored explicitly here below.
402
403 const Scalar xgw = priVars[switch1Idx];
404 const Scalar xgn = priVars[switch2Idx];
405 Scalar xgg = 1 - xgw - xgn;
406
407 // write mole fractions in the fluid state
408 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
409 fluidState_.setMoleFraction(gPhaseIdx, gCompIdx, xgg);
410 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
411
412 // calculate the composition of the remaining phases (as
413 // well as the densities of all phases). this is the job
414 // of the "ComputeFromReferencePhase" constraint solver ...
415 if (useConstraintSolver)
416 {
417 ComputeFromReferencePhase::solve(fluidState_,
418 paramCache,
419 gPhaseIdx);
420 }
421 // ... or calculated explicitly this way ...
422 else {
423
424 // note that the water phase is actually not existing!
425 // thus, this is used as phase switch criterion
426 Scalar xww = xgw * pg_
427 / (FluidSystem::fugacityCoefficient(fluidState_,
428 wPhaseIdx,wCompIdx)
429 * pw_);
430 Scalar xwn = 1.e-10;
431 Scalar xwg = 1.e-10;
432
433 // note that the NAPL phase is actually not existing!
434 // thus, this is used as phase switch criterion
435 Scalar xnn = xgn * pg_
436 / (FluidSystem::fugacityCoefficient(fluidState_,
437 nPhaseIdx,nCompIdx)
438 * pn_);
439 Scalar xna = 1.e-10;
440 Scalar xnw = 1.e-10;
441
442 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
443 fluidState_.setMoleFraction(wPhaseIdx, gCompIdx, xwg);
444 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
445 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
446 fluidState_.setMoleFraction(nPhaseIdx, gCompIdx, xna);
447 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
448
449 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
450 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
451 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
452 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
453 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
454 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
455
456 fluidState_.setDensity(wPhaseIdx, rhoW);
457 fluidState_.setDensity(gPhaseIdx, rhoG);
458 fluidState_.setDensity(nPhaseIdx, rhoN);
459 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
460 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
461 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
462 }
463 }
464 else if (phasePresence == wgPhaseOnly) {
465 // only water and gas phases are present
466 Scalar xgn = priVars[switch2Idx];
467 Scalar partPressH2O = fluidState_.fugacityCoefficient(wPhaseIdx, wCompIdx)*pw_;
468 if (partPressH2O > pg_) partPressH2O = pg_;
469
470 Scalar xgw = partPressH2O/pg_;
471 Scalar xgg = 1.-xgn-xgw;
472
473 // write mole fractions in the fluid state
474 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
475 fluidState_.setMoleFraction(gPhaseIdx, gCompIdx, xgg);
476 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
477
478 // calculate the composition of the remaining phases (as
479 // well as the densities of all phases). this is the job
480 // of the "ComputeFromReferencePhase" constraint solver ...
481 if (useConstraintSolver)
482 {
483 ComputeFromReferencePhase::solve(fluidState_,
484 paramCache,
485 gPhaseIdx);
486 }
487 // ... or calculated explicitly this way ...
488 else {
489 // actually, it's nothing else than Henry coefficient
490 Scalar xwn = xgn * pg_
491 / (FluidSystem::fugacityCoefficient(fluidState_,
492 wPhaseIdx,nCompIdx)
493 * pw_);
494 Scalar xwg = xgg * pg_
495 / (FluidSystem::fugacityCoefficient(fluidState_,
496 wPhaseIdx,gCompIdx)
497 * pw_);
498 Scalar xww = 1.-xwg-xwn;
499
500 // note that the NAPL phase is actually not existing!
501 // thus, this is used as phase switch criterion
502 Scalar xnn = xgn * pg_
503 / (FluidSystem::fugacityCoefficient(fluidState_,
504 nPhaseIdx,nCompIdx)
505 * pn_);
506 Scalar xna = 1.e-10;
507 Scalar xnw = 1.e-10;
508
509 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
510 fluidState_.setMoleFraction(wPhaseIdx, gCompIdx, xwg);
511 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
512 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
513 fluidState_.setMoleFraction(nPhaseIdx, gCompIdx, xna);
514 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
515
516 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
517 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
518 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
519 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
520 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
521 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
522
523 fluidState_.setDensity(wPhaseIdx, rhoW);
524 fluidState_.setDensity(gPhaseIdx, rhoG);
525 fluidState_.setDensity(nPhaseIdx, rhoN);
526 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
527 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
528 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
529 }
530 }
531 else
532 {
533 DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
534 }
535
536 for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx) {
537 // Mobilities
538 const Scalar mu =
539 FluidSystem::viscosity(fluidState_,
540 paramCache,
541 phaseIdx);
542 fluidState_.setViscosity(phaseIdx,mu);
543
544 Scalar kr;
545 kr = MaterialLaw::kr(materialParams, phaseIdx,
546 fluidState_.saturation(wPhaseIdx),
547 fluidState_.saturation(nPhaseIdx),
548 fluidState_.saturation(gPhaseIdx));
549 mobility_[phaseIdx] = kr / mu;
550 Valgrind::CheckDefined(mobility_[phaseIdx]);
551 }
552
553 // material dependent parameters for NAPL adsorption
554 bulkDensTimesAdsorpCoeff_ =
555 MaterialLaw::bulkDensTimesAdsorpCoeff(materialParams);
556
557 /* compute the diffusion coefficient
558 * \note This is the part of the diffusion coefficient determined by the fluid state, e.g.
559 * important if they are tabularized. In the diffusive flux computation (e.g. Fick's law)
560 * this gets converted into an effecient coefficient depending on saturation and porosity.
561 * We can then add a normalized tensorial component
562 * e.g. obtained from DTI from the spatial params (currently not implemented)
563 */
564
565 auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
566 {
567 return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
568 };
569
570 // porosity & permeabilty
571 updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
572
573 effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
574
575 EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
576 permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
577 EnergyVolVars::updateEffectiveThermalConductivity();
578
579 // compute and set the enthalpy
580 for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx)
581 {
582 Scalar h = EnergyVolVars::enthalpy(fluidState_, paramCache, phaseIdx);
583 fluidState_.setEnthalpy(phaseIdx, h);
584 }
585 }
586
590 const FluidState &fluidState() const
591 { return fluidState_; }
592
596 const SolidState &solidState() const
597 { return solidState_; }
598
604 Scalar averageMolarMass(int phaseIdx) const
605 { return fluidState_.averageMolarMass(phaseIdx); }
606
613 Scalar saturation(const int phaseIdx) const
614 { return fluidState_.saturation(phaseIdx); }
615
623 Scalar massFraction(const int phaseIdx, const int compIdx) const
624 { return fluidState_.massFraction(phaseIdx, compIdx); }
625
633 Scalar moleFraction(const int phaseIdx, const int compIdx) const
634 { return fluidState_.moleFraction(phaseIdx, compIdx); }
635
642 Scalar density(const int phaseIdx) const
643 { return fluidState_.density(phaseIdx); }
644
651 Scalar molarDensity(const int phaseIdx) const
652 { return fluidState_.molarDensity(phaseIdx); }
653
660 Scalar pressure(const int phaseIdx) const
661 { return fluidState_.pressure(phaseIdx); }
662
670 Scalar temperature() const
671 { return fluidState_.temperature(/*phaseIdx=*/0); }
672
679 Scalar mobility(const int phaseIdx) const
680 {
681 return mobility_[phaseIdx];
682 }
683
687 Scalar capillaryPressure() const
688 { return fluidState_.capillaryPressure(); }
689
693 Scalar porosity() const
694 { return solidState_.porosity(); }
695
699 Scalar bulkDensTimesAdsorpCoeff() const
700 { return bulkDensTimesAdsorpCoeff_; }
701
705 const PermeabilityType& permeability() const
706 { return permeability_; }
707
711 [[deprecated("Will be removed after release 3.2. Use diffusionCoefficient(phaseIdx, compIIdx, compJIdx)!")]]
712 Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
713 {
714 if (compIdx != phaseIdx)
715 return diffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
716 else
717 DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for phaseIdx = compIdx");
718 }
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 Scalar bulkDensTimesAdsorpCoeff_;
750
751 DiffusionCoefficients effectiveDiffCoeff_;
752};
753
754} // end namespace Dumux
755
756#endif
constants.hh
A central place for various physical constants occuring in some equations.
compositional.hh
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
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...
updatesolidvolumefractions.hh
Update the solid volume fractions (inert and reacitve) and set them in the solidstate.
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
Valgrind::CheckDefined
bool CheckDefined(const T &value)
Make valgrind complain if the object occupied by an object is undefined.
Definition: valgrind.hh:72
Dumux
Definition: adapt.hh:29
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::ComputeFromReferencePhase
Computes all quantities of a generic fluid state if a reference phase has been specified.
Definition: computefromreferencephase.hh:77
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:112
Dumux::MiscibleMultiPhaseComposition
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: misciblemultiphasecomposition.hh:60
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:82
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:50
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:623
Dumux::ThreePThreeCVolumeVariables::Indices
typename ModelTraits::Indices Indices
export the indices
Definition: porousmediumflow/3p3c/volumevariables.hh:97
Dumux::ThreePThreeCVolumeVariables::averageMolarMass
Scalar averageMolarMass(int phaseIdx) const
Returns the average molar mass of the fluid phase.
Definition: porousmediumflow/3p3c/volumevariables.hh:604
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:115
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:613
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:633
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:651
Dumux::ThreePThreeCVolumeVariables::SolidState
typename Traits::SolidState SolidState
export type of solid state
Definition: porousmediumflow/3p3c/volumevariables.hh:99
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:660
Dumux::ThreePThreeCVolumeVariables::solidState
const SolidState & solidState() const
Returns the phase state for the control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:596
Dumux::ThreePThreeCVolumeVariables::permeability
const PermeabilityType & permeability() const
Returns the average permeability within the control volume in .
Definition: porousmediumflow/3p3c/volumevariables.hh:705
Dumux::ThreePThreeCVolumeVariables::FluidState
typename Traits::FluidState FluidState
export fluid state type
Definition: porousmediumflow/3p3c/volumevariables.hh:93
Dumux::ThreePThreeCVolumeVariables::temperature
Scalar temperature() const
Returns temperature inside the sub-control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:670
Dumux::ThreePThreeCVolumeVariables::capillaryPressure
Scalar capillaryPressure() const
Returns the effective capillary pressure within the control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:687
Dumux::ThreePThreeCVolumeVariables::diffusionCoefficient
Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
Returns the diffusion coefficient.
Definition: porousmediumflow/3p3c/volumevariables.hh:712
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:679
Dumux::ThreePThreeCVolumeVariables::fluidState
const FluidState & fluidState() const
Returns the phase state for the control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:590
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:642
Dumux::ThreePThreeCVolumeVariables::FluidSystem
typename Traits::FluidSystem FluidSystem
export fluid system type
Definition: porousmediumflow/3p3c/volumevariables.hh:95
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:699
Dumux::ThreePThreeCVolumeVariables::porosity
Scalar porosity() const
Returns the average porosity within the control volume.
Definition: porousmediumflow/3p3c/volumevariables.hh:693
Dumux::ThreePThreeCVolumeVariables::SolidSystem
typename Traits::SolidSystem SolidSystem
export type of solid system
Definition: porousmediumflow/3p3c/volumevariables.hh:101
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