3.4
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
porousmediumflow/3pwateroil/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 *****************************************************************************/
26#ifndef DUMUX_3P2CNI_VOLUME_VARIABLES_HH
27#define DUMUX_3P2CNI_VOLUME_VARIABLES_HH
28
29#include <vector>
30#include <iostream>
31
32#include <dumux/common/math.hh>
35
39
42
44
45namespace Dumux {
46
47namespace Detail {
48// helper struct and function detecting if the fluid matrix interaction features a adsorptionModel() function
49#ifndef DOXYGEN // hide from doxygen
50template <class FluidMatrixInteraction>
51using AdsorptionModelDetector = decltype(std::declval<FluidMatrixInteraction>().adsorptionModel());
52#endif // DOXYGEN
53
54template<class FluidMatrixInteraction>
55static constexpr bool hasAdsorptionModel()
56{ return Dune::Std::is_detected<AdsorptionModelDetector, FluidMatrixInteraction>::value; }
57
58}
59
65template <class Traits>
68, public EnergyVolumeVariables<Traits, ThreePWaterOilVolumeVariables<Traits> >
69{
72 using Scalar = typename Traits::PrimaryVariables::value_type;
73 using ModelTraits = typename Traits::ModelTraits;
74 using FS = typename Traits::FluidSystem;
75 static constexpr int numFluidComps = ParentType::numFluidComponents();
76
77 enum {
79
80 wCompIdx = FS::wCompIdx,
81 nCompIdx = FS::nCompIdx,
82
83 wPhaseIdx = FS::wPhaseIdx,
84 gPhaseIdx = FS::gPhaseIdx,
85 nPhaseIdx = FS::nPhaseIdx,
86
87 switch1Idx = ModelTraits::Indices::switch1Idx,
88 switch2Idx = ModelTraits::Indices::switch2Idx,
89 pressureIdx = ModelTraits::Indices::pressureIdx
90 };
91
92 // present phases
93 enum {
94 threePhases = ModelTraits::Indices::threePhases,
95 wPhaseOnly = ModelTraits::Indices::wPhaseOnly,
96 gnPhaseOnly = ModelTraits::Indices::gnPhaseOnly,
97 wnPhaseOnly = ModelTraits::Indices::wnPhaseOnly,
98 gPhaseOnly = ModelTraits::Indices::gPhaseOnly,
99 wgPhaseOnly = ModelTraits::Indices::wgPhaseOnly
100 };
101
102 using EffDiffModel = typename Traits::EffectiveDiffusivityModel;
103 using DiffusionCoefficients = typename Traits::DiffusionType::DiffusionCoefficientsContainer;
104
105public:
107 using FluidState = typename Traits::FluidState;
109 using FluidSystem = typename Traits::FluidSystem;
111 using Indices = typename ModelTraits::Indices;
113 using SolidState = typename Traits::SolidState;
115 using SolidSystem = typename Traits::SolidSystem;
119 static constexpr bool onlyGasPhaseCanDisappear()
120 { return Traits::ModelTraits::onlyGasPhaseCanDisappear(); }
121
130 template<class ElemSol, class Problem, class Element, class Scv>
131 void update(const ElemSol &elemSol,
132 const Problem &problem,
133 const Element &element,
134 const Scv& scv)
135 {
136 ParentType::update(elemSol, problem, element, scv);
137 const auto& priVars = elemSol[scv.localDofIndex()];
138 const auto phasePresence = priVars.state();
139
140 if constexpr (!onlyGasPhaseCanDisappear())
141 {
142 /* first the saturations */
143 if (phasePresence == threePhases)
144 {
145 sw_ = priVars[switch1Idx];
146 sn_ = priVars[switch2Idx];
147 sg_ = 1. - sw_ - sn_;
148 }
149 else if (phasePresence == wPhaseOnly)
150 {
151 sw_ = 1.;
152 sn_ = 0.;
153 sg_ = 0.;
154 }
155 else if (phasePresence == gnPhaseOnly)
156 {
157 sw_ = 0.;
158 sn_ = priVars[switch1Idx];
159 sg_ = 1. - sn_;
160 }
161 else if (phasePresence == wnPhaseOnly)
162 {
163 sn_ = priVars[switch2Idx];
164 sw_ = 1. - sn_;
165 sg_ = 0.;
166 }
167 else if (phasePresence == gPhaseOnly)
168 {
169 sw_ = 0.;
170 sn_ = 0.;
171 sg_ = 1.;
172 }
173 else if (phasePresence == wgPhaseOnly)
174 {
175 sw_ = priVars[switch1Idx];
176 sn_ = 0.;
177 sg_ = 1. - sw_;
178 }
179 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
180
181 fluidState_.setSaturation(wPhaseIdx, sw_);
182 fluidState_.setSaturation(gPhaseIdx, sg_);
183 fluidState_.setSaturation(nPhaseIdx, sn_);
184
185 const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
186
187 // calculate capillary pressures
188 const Scalar pcgw = fluidMatrixInteraction.pcgw(sw_, sn_);
189 const Scalar pcnw = fluidMatrixInteraction.pcnw(sw_, sn_);
190 const Scalar pcgn = fluidMatrixInteraction.pcgn(sw_, sn_);
191
192 const Scalar pcAlpha = fluidMatrixInteraction.pcAlpha(sw_, sn_);
193 const Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
194
195 /* now the pressures */
196 if (phasePresence == threePhases || phasePresence == gnPhaseOnly || phasePresence == gPhaseOnly || phasePresence == wgPhaseOnly)
197 {
198 pg_ = priVars[pressureIdx];
199 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
200 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
201 }
202 else if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly)
203 {
204 pw_ = priVars[pressureIdx];
205 pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
206 pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
207 }
208 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
209
210 fluidState_.setPressure(wPhaseIdx, pw_);
211 fluidState_.setPressure(gPhaseIdx, pg_);
212 fluidState_.setPressure(nPhaseIdx, pn_);
213
214 /* now the temperature */
215 if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly || phasePresence == gPhaseOnly)
216 {
217 temp_ = priVars[switch1Idx];
218 }
219 else if (phasePresence == threePhases)
220 {
221 // temp from inverse pwsat and pnsat which have to sum up to pg
222 Scalar temp = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx); // initial guess
223 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
224 {
225 fluidState_.setTemperature(phaseIdx, temp);
226 }
227 solidState_.setTemperature(temp);
228 Scalar defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
229 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
230
231 using std::abs;
232 while(abs(defect) > 0.01) // simply a small number chosen ...
233 {
234 Scalar deltaT = 1.e-8 * temp;
235 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
236 {
237 fluidState_.setTemperature(phaseIdx, temp+deltaT);
238 }
239 solidState_.setTemperature(temp+deltaT);
240 Scalar fUp = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
241 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
242
243 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
244 {
245 fluidState_.setTemperature(phaseIdx, temp-deltaT);
246 }
247 solidState_.setTemperature(temp-deltaT);
248 Scalar fDown = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
249 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
250
251 temp = temp - defect * 2. * deltaT / (fUp - fDown);
252
253 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
254 {
255 fluidState_.setTemperature(phaseIdx, temp);
256 }
257 solidState_.setTemperature(temp);
258 defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
259 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
260 }
261 temp_ = temp;
262 }
263 else if (phasePresence == wgPhaseOnly)
264 {
265 // temp from inverse pwsat
266 temp_ = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
267 }
268 else if (phasePresence == gnPhaseOnly)
269 {
270 // temp from inverse pnsat
271 temp_ = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
272 }
273 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
274
275 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
276 {
277 fluidState_.setTemperature(phaseIdx, temp_);
278 }
279 solidState_.setTemperature(temp_);
280
281 // now comes the tricky part: calculate phase composition
282 if (phasePresence == threePhases) {
283
284 // all phases are present, phase compositions are a
285 // result of the the gas <-> liquid equilibrium.
286 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
287 Scalar partPressNAPL = pg_ - partPressH2O;
288
289 Scalar xgn = partPressNAPL/pg_;
290 Scalar xgw = partPressH2O/pg_;
291
292 // Henry
293 Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
294 Scalar xww = 1.-xwn;
295
296 // Not yet filled with real numbers for the NAPL phase
297 Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
298 Scalar xnn = 1.-xnw;
299
300 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
301 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
302 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
303 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
304 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
305 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
306
307 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
308 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
309 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
310 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
311 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
312 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
313
314 fluidState_.setDensity(wPhaseIdx, rhoW);
315 fluidState_.setDensity(gPhaseIdx, rhoG);
316 fluidState_.setDensity(nPhaseIdx, rhoN);
317 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
318 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
319 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
320 }
321 else if (phasePresence == wPhaseOnly) {
322 // only the water phase is present, water phase composition is
323 // stored explicitly.
324
325 // extract mole fractions in the water phase
326 Scalar xwn = priVars[switch2Idx];
327 Scalar xww = 1 - xwn;
328
329 // write water mole fractions in the fluid state
330 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
331 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
332
333 // note that the gas phase is actually not existing!
334 // thus, this is used as phase switch criterion
335 Scalar xgn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
336 Scalar xgw = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
337
338
339 // note that the NAPL phase is actually not existing!
340 // thus, this is used as phase switch criterion
341 // maybe solubility would be better than this approach via Henry
342 Scalar xnn = xwn * FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx) / (xgn * pg_);
343 Scalar xnw = xgw*pg_ / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
344
345 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
346 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
347 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
348 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
349
350 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
351 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
352 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
353 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
354 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
355 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
356
357 fluidState_.setDensity(wPhaseIdx, rhoW);
358 fluidState_.setDensity(gPhaseIdx, rhoG);
359 fluidState_.setDensity(nPhaseIdx, rhoN);
360 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
361 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
362 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
363 }
364 else if (phasePresence == gnPhaseOnly) {
365
366 // only gas and NAPL phases are present
367
368 Scalar xnw = priVars[switch2Idx];
369 Scalar xnn = 1.-xnw;
370 Scalar xgn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
371 Scalar xgw = 1.-xgn;
372
373 // note that the water phase is actually not present
374 // the values are used as switching criteria
375 Scalar xww = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
376
377 // write mole fractions in the fluid state
378 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
379 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
380 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
381 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
382 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
383
384 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
385 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
386 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
387 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
388 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
389 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
390
391 fluidState_.setDensity(wPhaseIdx, rhoW);
392 fluidState_.setDensity(gPhaseIdx, rhoG);
393 fluidState_.setDensity(nPhaseIdx, rhoN);
394 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
395 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
396 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
397
398 }
399 else if (phasePresence == wnPhaseOnly) {
400 // water and NAPL are present, phase compositions are a
401 // mole fractions of non-existing gas phase are used as switching criteria
402 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
403 Scalar partPressNAPL = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
404
405 Scalar xgn = partPressNAPL/pg_;
406 Scalar xgw = partPressH2O/pg_;
407
408 // Henry
409 Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
410 Scalar xww = 1.-xwn;
411
412 Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
413 Scalar xnn = 1.-xnw;
414
415 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
416 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
417 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
418 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
419 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
420 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
421
422 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
423 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
424 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
425 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
426 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
427 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
428
429 fluidState_.setDensity(wPhaseIdx, rhoW);
430 fluidState_.setDensity(gPhaseIdx, rhoG);
431 fluidState_.setDensity(nPhaseIdx, rhoN);
432 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
433 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
434 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
435 }
436 else if (phasePresence == gPhaseOnly) {
437 // only the gas phase is present, gas phase composition is
438 // stored explicitly here below.
439
440 const Scalar xgn = priVars[switch2Idx];
441 Scalar xgw = 1 - xgn;
442
443 // write mole fractions in the fluid state
444 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
445 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
446
447 // note that the water and NAPL phase is actually not present
448 // the values are used as switching criteria
449 Scalar xww = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
450 Scalar xnn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
451
452 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
453 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
454
455 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
456 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
457 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
458 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
459 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
460 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
461
462 fluidState_.setDensity(wPhaseIdx, rhoW);
463 fluidState_.setDensity(gPhaseIdx, rhoG);
464 fluidState_.setDensity(nPhaseIdx, rhoN);
465 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
466 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
467 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
468 }
469 else if (phasePresence == wgPhaseOnly) {
470 // only water and gas phases are present
471 const Scalar xgn = priVars[switch2Idx];
472 Scalar xgw = 1 - xgn;
473
474 // write mole fractions in the fluid state
475 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
476 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
477
478
479 Scalar xwn = xgn*pg_/FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
480 Scalar xww = 1.-xwn;
481
482 // write mole fractions in the fluid state
483 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
484 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
485
486 // note that the NAPL phase is actually not existing!
487 // thus, this is used as phase switch criterion
488 Scalar xnn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
489
490 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
491
492 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
493 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
494 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
495 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
496 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
497 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
498
499 fluidState_.setDensity(wPhaseIdx, rhoW);
500 fluidState_.setDensity(gPhaseIdx, rhoG);
501 fluidState_.setDensity(nPhaseIdx, rhoN);
502 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
503 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
504 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
505 }
506 else
507 assert(false); // unhandled phase state
508 } // end of if(!UseSimpleModel), i.e. the more complex version with six phase states
509
510 else // use the simpler model with only two phase states
511 {
512 /* first the saturations */
513 if (phasePresence == threePhases)
514 {
515 sw_ = priVars[switch1Idx];
516 sn_ = priVars[switch2Idx];
517 sg_ = 1. - sw_ - sn_;
518 }
519 else if (phasePresence == wnPhaseOnly)
520 {
521 sn_ = priVars[switch2Idx];
522 sw_ = 1. - sn_;
523 sg_ = 0.;
524 }
525 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
526
527 fluidState_.setSaturation(wPhaseIdx, sw_);
528 fluidState_.setSaturation(gPhaseIdx, sg_);
529 fluidState_.setSaturation(nPhaseIdx, sn_);
530
531 const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
532
533 // calculate capillary pressures
534 const Scalar pcgw = fluidMatrixInteraction.pcgw(sw_, sn_);
535 const Scalar pcnw = fluidMatrixInteraction.pcnw(sw_, sn_);
536 const Scalar pcgn = fluidMatrixInteraction.pcgn(sw_, sn_);
537
538 const Scalar pcAlpha = fluidMatrixInteraction.pcAlpha(sw_, sn_);
539 const Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
540
541 /* now the pressures */
542 if (phasePresence == threePhases)
543 {
544 pg_ = priVars[pressureIdx];
545 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
546 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
547 }
548 else if (phasePresence == wnPhaseOnly)
549 {
550 pw_ = priVars[pressureIdx];
551 pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
552 pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
553 }
554 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
555
556 fluidState_.setPressure(wPhaseIdx, pw_);
557 fluidState_.setPressure(gPhaseIdx, pg_);
558 fluidState_.setPressure(nPhaseIdx, pn_);
559
560 /* now the temperature */
561 if (phasePresence == wnPhaseOnly)
562 {
563 temp_ = priVars[switch1Idx];
564 }
565 else if (phasePresence == threePhases)
566 {
567 if(sn_<=1.e-10) // this threshold values is chosen arbitrarily as a small number
568 {
569 Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
570 temp_ = tempOnlyWater;
571 }
572 if(sw_<=1.e-10) // this threshold values is chosen arbitrarily as a small number
573 {
574 Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
575 temp_ = tempOnlyNAPL;
576 }
577 else
578 {
579 // temp from inverse pwsat and pnsat which have to sum up to pg
580 Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
581 Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
582 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
583 {
584 fluidState_.setTemperature(phaseIdx, tempOnlyWater);
585 }
586 solidState_.setTemperature(tempOnlyWater);
587 Scalar defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
588 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
589
590 Scalar temp = tempOnlyWater; // initial guess
591 int counter = 0;
592 using std::abs;
593 while(abs(defect) > 0.01) // simply a small number chosen ...
594 {
595 Scalar deltaT = 1.e-6; // fixed number, but T should always be in the order of a few hundred Kelvin
596 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
597 {
598 fluidState_.setTemperature(phaseIdx, temp+deltaT);
599 }
600 solidState_.setTemperature(temp+deltaT);
601 Scalar fUp = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
602 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
603
604 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
605 {
606 fluidState_.setTemperature(phaseIdx, temp-deltaT);
607 }
608 solidState_.setTemperature(temp-deltaT);
609 Scalar fDown = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
610 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
611
612 temp = temp - defect * 2. * deltaT / (fUp - fDown);
613
614 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
615 {
616 fluidState_.setTemperature(phaseIdx, temp);
617 }
618 solidState_.setTemperature(temp);
619 defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
620 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
621 counter +=1;
622 if (counter>10) break;
623 }
624 if ((sw_>1.e-10)&&(sw_<0.01))
625 temp = temp + (sw_ - 1.e-10) * (temp - tempOnlyNAPL) / (0.01 - 1.e-10);
626 if ((sn_>1.e-10)&&(sn_<0.01))
627 temp = temp + (sn_ - 1.e-10) * (temp - tempOnlyWater) / (0.01 - 1.e-10);
628 temp_ = temp;
629 }
630 }
631 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
632
633 for (int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
634 fluidState_.setTemperature(phaseIdx, temp_);
635
636 solidState_.setTemperature(temp_);
637
638 // now comes the tricky part: calculate phase composition
639 if (phasePresence == threePhases) {
640
641 // all phases are present, phase compositions are a
642 // result of the the gas <-> liquid equilibrium.
643 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
644 Scalar partPressNAPL = pg_ - partPressH2O;
645 // regularized evaporation for small liquid phase saturations
646 // avoids negative saturations of liquid phases
647 if (sw_<0.02) partPressH2O *= sw_/0.02;
648 if (partPressH2O < 0.) partPressH2O = 0;
649 if (sn_<0.02) partPressNAPL *= sn_ / 0.02;
650 if (partPressNAPL < 0.) partPressNAPL = 0;
651
652 Scalar xgn = partPressNAPL/pg_;
653 Scalar xgw = partPressH2O/pg_;
654
655 // Immiscible liquid phases, mole fractions are just dummy values
656 Scalar xwn = 0;
657 Scalar xww = 1.-xwn;
658
659 // Not yet filled with real numbers for the NAPL phase
660 Scalar xnw = 0;
661 Scalar xnn = 1.-xnw;
662
663 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
664 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
665 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
666 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
667 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
668 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
669
670 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
671 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
672 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
673 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
674 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
675 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
676
677 fluidState_.setDensity(wPhaseIdx, rhoW);
678 fluidState_.setDensity(gPhaseIdx, rhoG);
679 fluidState_.setDensity(nPhaseIdx, rhoN);
680 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
681 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
682 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
683 }
684 else if (phasePresence == wnPhaseOnly) {
685 // mole fractions of non-existing gas phase are used as switching criteria
686 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
687 Scalar partPressNAPL = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
688
689 Scalar xgn = partPressNAPL/pg_;
690 Scalar xgw = partPressH2O/pg_;
691
692 // immiscible liquid phases, mole fractions are just dummy values
693 Scalar xwn = 0;
694 Scalar xww = 1.-xwn;
695
696 Scalar xnw = 0;
697 Scalar xnn = 1.-xnw;
698
699 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
700 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
701 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
702 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
703 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
704 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
705
706 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
707 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
708 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
709 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
710 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
711 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
712
713 fluidState_.setDensity(wPhaseIdx, rhoW);
714 fluidState_.setDensity(gPhaseIdx, rhoG);
715 fluidState_.setDensity(nPhaseIdx, rhoN);
716 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
717 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
718 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
719 }
720 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
721 }
722
723 const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
724 for (int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx)
725 {
726 // Mobilities
727 const Scalar mu =
729 phaseIdx);
730 fluidState_.setViscosity(phaseIdx,mu);
731
732 const Scalar kr = fluidMatrixInteraction.kr(phaseIdx,
733 fluidState_.saturation(wPhaseIdx),
734 fluidState_.saturation(nPhaseIdx));
735 mobility_[phaseIdx] = kr / mu;
736 }
737
738 // material dependent parameters for NAPL adsorption (only if law is provided)
739 if constexpr (Detail::hasAdsorptionModel<std::decay_t<decltype(fluidMatrixInteraction)>>())
740 bulkDensTimesAdsorpCoeff_ = fluidMatrixInteraction.adsorptionModel().bulkDensTimesAdsorpCoeff();
741
742 // porosity
743 updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
744 EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
745
746 auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
747 {
748 return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
749 };
750
751 effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
752
753 // permeability
754 permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
755
756 fluidState_.setTemperature(temp_);
757 // the enthalpies (internal energies are directly calculated in the fluidstate
758 for (int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx)
759 {
760 Scalar h = FluidSystem::enthalpy(fluidState_, phaseIdx);
761 fluidState_.setEnthalpy(phaseIdx, h);
762 }
763
764 EnergyVolVars::updateEffectiveThermalConductivity();
765 }
766
770 const FluidState &fluidState() const
771 { return fluidState_; }
772
776 const SolidState &solidState() const
777 { return solidState_; }
778
784 Scalar averageMolarMass(int phaseIdx) const
785 { return fluidState_.averageMolarMass(phaseIdx); }
786
793 Scalar saturation(const int phaseIdx) const
794 { return fluidState_.saturation(phaseIdx); }
795
803 Scalar massFraction(const int phaseIdx, const int compIdx) const
804 { return fluidState_.massFraction(phaseIdx, compIdx); }
805
813 Scalar moleFraction(const int phaseIdx, const int compIdx) const
814 { return fluidState_.moleFraction(phaseIdx, compIdx); }
815
822 Scalar density(const int phaseIdx) const
823 { return fluidState_.density(phaseIdx); }
824
831 Scalar molarDensity(const int phaseIdx) const
832 { return fluidState_.molarDensity(phaseIdx); }
833
840 Scalar pressure(const int phaseIdx) const
841 { return fluidState_.pressure(phaseIdx); }
842
850 Scalar temperature() const
851 { return fluidState_.temperature(/*phaseIdx=*/0); }
852
859 Scalar mobility(const int phaseIdx) const
860 { return mobility_[phaseIdx]; }
861
862 Scalar viscosity(const int phaseIdx) const
863 { return fluidState_.viscosity(phaseIdx); }
864
868 Scalar capillaryPressure() const
869 { return fluidState_.capillaryPressure(); }
870
874 Scalar porosity() const
875 { return solidState_.porosity(); }
876
880 Scalar permeability() const
881 { return permeability_; }
882
883 /*
884 * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
885 */
886 Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
887 {
888 if (phaseIdx != nPhaseIdx)
889 return FluidSystem::diffusionCoefficient(fluidState_, phaseIdx);
890 else
891 return 1.e-10;
892 }
893
897 Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
898 { return effectiveDiffCoeff_(phaseIdx, compIIdx, compJIdx); }
899
904 {
905 if (bulkDensTimesAdsorpCoeff_)
906 return bulkDensTimesAdsorpCoeff_.value();
907 else
908 DUNE_THROW(Dune::NotImplemented, "Your spatialParams do not provide an adsorption model");
909 }
910
917 Scalar internalEnergy(int phaseIdx) const
918 { return fluidState_.internalEnergy(phaseIdx); };
919
926 Scalar enthalpy(int phaseIdx) const
927 { return fluidState_.enthalpy(phaseIdx); };
928
929protected:
932
933private:
934 Scalar sw_, sg_, sn_, pg_, pw_, pn_, temp_;
935
936 Scalar moleFrac_[numPs][numFluidComps];
937 Scalar massFrac_[numPs][numFluidComps];
938
939 Scalar permeability_;
940 Scalar mobility_[numPs];
941 OptionalScalar<Scalar> bulkDensTimesAdsorpCoeff_;
942
944 DiffusionCoefficients effectiveDiffCoeff_;
945
946};
947} // end namespace Dumux
948
949#endif
Some exceptions thrown in DuMux
Define some often used mathematical functions.
A wrapper that can either contain a valid Scalar or NaN.
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Update the solid volume fractions (inert and reacitve) and set them in the solidstate.
A central place for various physical constants occuring in some equations.
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
Definition: adapt.hh:29
decltype(std::declval< FluidMatrixInteraction >().adsorptionModel()) AdsorptionModelDetector
Definition: porousmediumflow/3p3c/volumevariables.hh:45
static constexpr bool hasAdsorptionModel()
Definition: porousmediumflow/3p3c/volumevariables.hh:48
std::string phasePresence() noexcept
I/O name of phase presence.
Definition: name.hh:147
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
std::string molarDensity(int phaseIdx) noexcept
I/O name of molar density for multiphase systems.
Definition: name.hh:83
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
T value() const
Definition: optionalscalar.hh:48
The primary variable switch controlling the phase presence state variable.
Definition: 3pwateroil/primaryvariableswitch.hh:39
Contains the quantities which are are constant within a finite volume in the three-phase,...
Definition: porousmediumflow/3pwateroil/volumevariables.hh:69
const FluidState & fluidState() const
Returns the phase state for the control-volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:770
SolidState solidState_
Definition: porousmediumflow/3pwateroil/volumevariables.hh:931
Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
Definition: porousmediumflow/3pwateroil/volumevariables.hh:886
Scalar molarDensity(const int phaseIdx) const
Returns the molar density of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:831
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/3pwateroil/volumevariables.hh:803
typename Traits::FluidState FluidState
The type of the object returned by the fluidState() method.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:107
Scalar permeability() const
Returns the permeability within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:880
typename Traits::FluidSystem FluidSystem
The type of the fluid system.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:109
Scalar viscosity(const int phaseIdx) const
Definition: porousmediumflow/3pwateroil/volumevariables.hh:862
Scalar bulkDensTimesAdsorpCoeff() const
Returns the adsorption information.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:903
typename Traits::SolidState SolidState
Export type of solid state.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:113
Scalar porosity() const
Returns the average porosity within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:874
Scalar enthalpy(int phaseIdx) const
Returns the total enthalpy of a phase in the sub-control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:926
Scalar mobility(const int phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:859
Scalar capillaryPressure() const
Returns the effective capillary pressure within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:868
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/3pwateroil/volumevariables.hh:813
Scalar saturation(const int phaseIdx) const
Returns the effective saturation of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:793
Scalar density(const int phaseIdx) const
Returns the mass density of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:822
typename Traits::SolidSystem SolidSystem
Export type of solid system.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:115
void update(const ElemSol &elemSol, const Problem &problem, const Element &element, const Scv &scv)
Updates all quantities for a given control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:131
Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
Returns the effective diffusion coefficients for a phase in .
Definition: porousmediumflow/3pwateroil/volumevariables.hh:897
static constexpr bool onlyGasPhaseCanDisappear()
State if only the gas phase is allowed to disappear.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:119
Scalar temperature() const
Returns temperature inside the sub-control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:850
Scalar internalEnergy(int phaseIdx) const
Returns the total internal energy of a phase in the sub-control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:917
typename ModelTraits::Indices Indices
Export the indices.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:111
const SolidState & solidState() const
Returns the phase state for the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:776
Scalar averageMolarMass(int phaseIdx) const
Returns the average molar mass of the fluid phase.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:784
FluidState fluidState_
Definition: porousmediumflow/3pwateroil/volumevariables.hh:927
Scalar pressure(const int phaseIdx) const
Returns the effective pressure of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:840
Definition: porousmediumflow/nonisothermal/volumevariables.hh:75
The isothermal base class.
Definition: porousmediumflow/volumevariables.hh:40
static constexpr int numFluidComponents()
Return number of components considered by the model.
Definition: porousmediumflow/volumevariables.hh:52
const PrimaryVariables & priVars() const
Returns the vector of primary variables.
Definition: porousmediumflow/volumevariables.hh:76
static constexpr int numFluidPhases()
Return number of phases considered by the model.
Definition: porousmediumflow/volumevariables.hh:50
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
Base class for the model specific class which provides access to all volume averaged quantities.
Base class for the model specific class which provides access to all volume averaged quantities.
The primary variable switch for the extended Richards model.