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