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