3.1-git
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
52template <class Traits>
55, public EnergyVolumeVariables<Traits, ThreePWaterOilVolumeVariables<Traits> >
56{
59 using Scalar = typename Traits::PrimaryVariables::value_type;
60 using ModelTraits = typename Traits::ModelTraits;
61 using FS = typename Traits::FluidSystem;
62 static constexpr int numFluidComps = ParentType::numFluidComponents();
63
64 enum {
66
67 wCompIdx = FS::wCompIdx,
68 nCompIdx = FS::nCompIdx,
69
70 wPhaseIdx = FS::wPhaseIdx,
71 gPhaseIdx = FS::gPhaseIdx,
72 nPhaseIdx = FS::nPhaseIdx,
73
74 switch1Idx = ModelTraits::Indices::switch1Idx,
75 switch2Idx = ModelTraits::Indices::switch2Idx,
76 pressureIdx = ModelTraits::Indices::pressureIdx
77 };
78
79 // present phases
80 enum {
81 threePhases = ModelTraits::Indices::threePhases,
82 wPhaseOnly = ModelTraits::Indices::wPhaseOnly,
83 gnPhaseOnly = ModelTraits::Indices::gnPhaseOnly,
84 wnPhaseOnly = ModelTraits::Indices::wnPhaseOnly,
85 gPhaseOnly = ModelTraits::Indices::gPhaseOnly,
86 wgPhaseOnly = ModelTraits::Indices::wgPhaseOnly
87 };
88
89public:
91 using FluidState = typename Traits::FluidState;
93 using FluidSystem = typename Traits::FluidSystem;
95 using Indices = typename ModelTraits::Indices;
97 using SolidState = typename Traits::SolidState;
99 using SolidSystem = typename Traits::SolidSystem;
103 static constexpr bool onlyGasPhaseCanDisappear()
104 { return Traits::ModelTraits::onlyGasPhaseCanDisappear(); }
105
114 template<class ElemSol, class Problem, class Element, class Scv>
115 void update(const ElemSol &elemSol,
116 const Problem &problem,
117 const Element &element,
118 const Scv& scv)
119 {
120 ParentType::update(elemSol, problem, element, scv);
121 const auto& priVars = elemSol[scv.localDofIndex()];
122 const auto phasePresence = priVars.state();
123
124 // capillary pressure parameters
125 using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
126 const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
127
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.");
169
170 fluidState_.setSaturation(wPhaseIdx, sw_);
171 fluidState_.setSaturation(gPhaseIdx, sg_);
172 fluidState_.setSaturation(nPhaseIdx, sn_);
173
174 /* now the pressures */
175 if (phasePresence == threePhases || phasePresence == gnPhaseOnly || phasePresence == gPhaseOnly || phasePresence == wgPhaseOnly)
176 {
177 pg_ = priVars[pressureIdx];
178
179 // calculate capillary pressures
180 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
181 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
182 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
183
184 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
185 Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
186
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
194 // calculate capillary pressures
195 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
196 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
197 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
198
199 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
200 Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
201
202 pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
203 pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
204 }
205 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
207
208 fluidState_.setPressure(wPhaseIdx, pw_);
209 fluidState_.setPressure(gPhaseIdx, pg_);
210 fluidState_.setPressure(nPhaseIdx, pn_);
211
212 /* now the temperature */
213 if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly || phasePresence == gPhaseOnly)
214 {
215 temp_ = priVars[switch1Idx];
216 }
217 else if (phasePresence == threePhases)
218 {
219 // temp from inverse pwsat and pnsat which have to sum up to pg
220 Scalar temp = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx); // initial guess
221 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
222 {
223 fluidState_.setTemperature(phaseIdx, temp);
224 }
225 solidState_.setTemperature(temp);
226 Scalar defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
227 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
228
229 using std::abs;
230 while(abs(defect) > 0.01) // simply a small number chosen ...
231 {
232 Scalar deltaT = 1.e-8 * temp;
233 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
234 {
235 fluidState_.setTemperature(phaseIdx, temp+deltaT);
236 }
237 solidState_.setTemperature(temp+deltaT);
238 Scalar fUp = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
239 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
240
241 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
242 {
243 fluidState_.setTemperature(phaseIdx, temp-deltaT);
244 }
245 solidState_.setTemperature(temp-deltaT);
246 Scalar fDown = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
247 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
248
249 temp = temp - defect * 2. * deltaT / (fUp - fDown);
250
251 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
252 {
253 fluidState_.setTemperature(phaseIdx, temp);
254 }
255 solidState_.setTemperature(temp);
256 defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
257 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
258 }
259 temp_ = temp;
260 }
261 else if (phasePresence == wgPhaseOnly)
262 {
263 // temp from inverse pwsat
264 temp_ = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
265 }
266 else if (phasePresence == gnPhaseOnly)
267 {
268 // temp from inverse pnsat
269 temp_ = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
270 }
271 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
273
274 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
275 {
276 fluidState_.setTemperature(phaseIdx, temp_);
277 }
278 solidState_.setTemperature(temp_);
279
280 // now comes the tricky part: calculate phase composition
281 if (phasePresence == threePhases) {
282
283 // all phases are present, phase compositions are a
284 // result of the the gas <-> liquid equilibrium.
285 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
286 Scalar partPressNAPL = pg_ - partPressH2O;
287
288 Scalar xgn = partPressNAPL/pg_;
289 Scalar xgw = partPressH2O/pg_;
290
291 // Henry
292 Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
293 Scalar xww = 1.-xwn;
294
295 // Not yet filled with real numbers for the NAPL phase
296 Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
297 Scalar xnn = 1.-xnw;
298
299 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
300 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
301 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
302 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
303 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
304 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
305
306 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
307 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
308 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
309 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
310 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
311 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
312
313 fluidState_.setDensity(wPhaseIdx, rhoW);
314 fluidState_.setDensity(gPhaseIdx, rhoG);
315 fluidState_.setDensity(nPhaseIdx, rhoN);
316 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
317 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
318 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
319 }
320 else if (phasePresence == wPhaseOnly) {
321 // only the water phase is present, water phase composition is
322 // stored explicitly.
323
324 // extract mole fractions in the water phase
325 Scalar xwn = priVars[switch2Idx];
326 Scalar xww = 1 - xwn;
327
328 // write water mole fractions in the fluid state
329 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
330 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
331
332 // note that the gas phase is actually not existing!
333 // thus, this is used as phase switch criterion
334 Scalar xgn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
335 Scalar xgw = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
336
337
338 // note that the NAPL phase is actually not existing!
339 // thus, this is used as phase switch criterion
340 // maybe solubility would be better than this approach via Henry
341 Scalar xnn = xwn * FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx) / (xgn * pg_);
342 Scalar xnw = xgw*pg_ / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
343
344 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
345 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
346 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
347 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
348
349 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
350 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
351 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
352 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
353 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
354 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
355
356 fluidState_.setDensity(wPhaseIdx, rhoW);
357 fluidState_.setDensity(gPhaseIdx, rhoG);
358 fluidState_.setDensity(nPhaseIdx, rhoN);
359 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
360 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
361 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
362 }
363 else if (phasePresence == gnPhaseOnly) {
364
365 // only gas and NAPL phases are present
366
367 Scalar xnw = priVars[switch2Idx];
368 Scalar xnn = 1.-xnw;
369 Scalar xgn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
370 Scalar xgw = 1.-xgn;
371
372 // note that the water phase is actually not present
373 // the values are used as switching criteria
374 Scalar xww = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
375
376 // write mole fractions in the fluid state
377 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
378 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
379 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
380 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
381 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
382
383 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
384 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
385 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
386 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
387 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
388 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
389
390 fluidState_.setDensity(wPhaseIdx, rhoW);
391 fluidState_.setDensity(gPhaseIdx, rhoG);
392 fluidState_.setDensity(nPhaseIdx, rhoN);
393 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
394 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
395 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
396
397 }
398 else if (phasePresence == wnPhaseOnly) {
399 // water and NAPL are present, phase compositions are a
400 // mole fractions of non-existing gas phase are used as switching criteria
401 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
402 Scalar partPressNAPL = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
403
404 Scalar xgn = partPressNAPL/pg_;
405 Scalar xgw = partPressH2O/pg_;
406
407 // Henry
408 Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
409 Scalar xww = 1.-xwn;
410
411 Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(fluidState_, nPhaseIdx,wCompIdx);
412 Scalar xnn = 1.-xnw;
413
414 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
415 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
416 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
417 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
418 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
419 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
420
421 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
422 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
423 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
424 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
425 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
426 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
427
428 fluidState_.setDensity(wPhaseIdx, rhoW);
429 fluidState_.setDensity(gPhaseIdx, rhoG);
430 fluidState_.setDensity(nPhaseIdx, rhoN);
431 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
432 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
433 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
434 }
435 else if (phasePresence == gPhaseOnly) {
436 // only the gas phase is present, gas phase composition is
437 // stored explicitly here below.
438
439 const Scalar xgn = priVars[switch2Idx];
440 Scalar xgw = 1 - xgn;
441
442 // write mole fractions in the fluid state
443 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
444 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
445
446 // note that the water and NAPL phase is actually not present
447 // the values are used as switching criteria
448 Scalar xww = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) / pg_;
449 Scalar xnn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
450
451 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
452 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
453
454 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
455 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
456 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
457 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
458 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
459 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
460
461 fluidState_.setDensity(wPhaseIdx, rhoW);
462 fluidState_.setDensity(gPhaseIdx, rhoG);
463 fluidState_.setDensity(nPhaseIdx, rhoN);
464 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
465 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
466 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
467 }
468 else if (phasePresence == wgPhaseOnly) {
469 // only water and gas phases are present
470 const Scalar xgn = priVars[switch2Idx];
471 Scalar xgw = 1 - xgn;
472
473 // write mole fractions in the fluid state
474 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
475 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
476
477
478 Scalar xwn = xgn*pg_/FluidSystem::henryCoefficient(fluidState_, wPhaseIdx,nCompIdx);
479 Scalar xww = 1.-xwn;
480
481 // write mole fractions in the fluid state
482 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
483 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
484
485 // note that the NAPL phase is actually not existing!
486 // thus, this is used as phase switch criterion
487 Scalar xnn = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx) / pg_;
488
489 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
490
491 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
492 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
493 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
494 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
495 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
496 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
497
498 fluidState_.setDensity(wPhaseIdx, rhoW);
499 fluidState_.setDensity(gPhaseIdx, rhoG);
500 fluidState_.setDensity(nPhaseIdx, rhoN);
501 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
502 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
503 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
504 }
505 else
506 assert(false); // unhandled phase state
507 } // end of if(!UseSimpleModel), i.e. the more complex version with six phase states
508 else // use the simpler model with only two phase states
509 {
510 /* first the saturations */
511 if (phasePresence == threePhases)
512 {
513 sw_ = priVars[switch1Idx];
514 sn_ = priVars[switch2Idx];
515 sg_ = 1. - sw_ - sn_;
516 }
517 else if (phasePresence == wnPhaseOnly)
518 {
519 sn_ = priVars[switch2Idx];
520 sw_ = 1. - sn_;
521 sg_ = 0.;
522 }
523 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
525
526 fluidState_.setSaturation(wPhaseIdx, sw_);
527 fluidState_.setSaturation(gPhaseIdx, sg_);
528 fluidState_.setSaturation(nPhaseIdx, sn_);
529
530 /* now the pressures */
531 if (phasePresence == threePhases)
532 {
533 pg_ = priVars[pressureIdx];
534
535 // calculate capillary pressures
536 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
537 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
538 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
539
540 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
541 Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
542
543 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
544 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
545 }
546 else if (phasePresence == wnPhaseOnly)
547 {
548 pw_ = priVars[pressureIdx];
549
550 // calculate capillary pressures
551 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
552 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
553 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
554
555 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
556 Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
557
558 pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
559 pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
560 }
561 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
563
564 fluidState_.setPressure(wPhaseIdx, pw_);
565 fluidState_.setPressure(gPhaseIdx, pg_);
566 fluidState_.setPressure(nPhaseIdx, pn_);
567
568 /* now the temperature */
569 if (phasePresence == wnPhaseOnly)
570 {
571 temp_ = priVars[switch1Idx];
572 }
573 else if (phasePresence == threePhases)
574 {
575 if(sn_<=1.e-10) // this threshold values is chosen arbitrarily as a small number
576 {
577 Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
578 temp_ = tempOnlyWater;
579 }
580 if(sw_<=1.e-10) // this threshold values is chosen arbitrarily as a small number
581 {
582 Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
583 temp_ = tempOnlyNAPL;
584 }
585 else
586 {
587 // temp from inverse pwsat and pnsat which have to sum up to pg
588 Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx);
589 Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx);
590 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
591 {
592 fluidState_.setTemperature(phaseIdx, tempOnlyWater);
593 }
594 solidState_.setTemperature(tempOnlyWater);
595 Scalar defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
596 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
597
598 Scalar temp = tempOnlyWater; // initial guess
599 int counter = 0;
600 using std::abs;
601 while(abs(defect) > 0.01) // simply a small number chosen ...
602 {
603 Scalar deltaT = 1.e-6; // fixed number, but T should always be in the order of a few hundred Kelvin
604 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
605 {
606 fluidState_.setTemperature(phaseIdx, temp+deltaT);
607 }
608 solidState_.setTemperature(temp+deltaT);
609 Scalar fUp = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
610 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
611
612 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
613 {
614 fluidState_.setTemperature(phaseIdx, temp-deltaT);
615 }
616 solidState_.setTemperature(temp-deltaT);
617 Scalar fDown = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
618 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
619
620 temp = temp - defect * 2. * deltaT / (fUp - fDown);
621
622 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
623 {
624 fluidState_.setTemperature(phaseIdx, temp);
625 }
626 solidState_.setTemperature(temp);
627 defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx)
628 - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
629 counter +=1;
630 if (counter>10) break;
631 }
632 if ((sw_>1.e-10)&&(sw_<0.01))
633 temp = temp + (sw_ - 1.e-10) * (temp - tempOnlyNAPL) / (0.01 - 1.e-10);
634 if ((sn_>1.e-10)&&(sn_<0.01))
635 temp = temp + (sn_ - 1.e-10) * (temp - tempOnlyWater) / (0.01 - 1.e-10);
636 temp_ = temp;
637 }
638 }
639 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
641
642 for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
643 {
644 fluidState_.setTemperature(phaseIdx, temp_);
645 }
646 solidState_.setTemperature(temp_);
647
648 // now comes the tricky part: calculate phase composition
649 if (phasePresence == threePhases) {
650
651 // all phases are present, phase compositions are a
652 // result of the the gas <-> liquid equilibrium.
653 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
654 Scalar partPressNAPL = pg_ - partPressH2O;
655 // regularized evaporation for small liquid phase saturations
656 // avoids negative saturations of liquid phases
657 if (sw_<0.02) partPressH2O *= sw_/0.02;
658 if (partPressH2O < 0.) partPressH2O = 0;
659 if (sn_<0.02) partPressNAPL *= sn_ / 0.02;
660 if (partPressNAPL < 0.) partPressNAPL = 0;
661
662 Scalar xgn = partPressNAPL/pg_;
663 Scalar xgw = partPressH2O/pg_;
664
665 // Immiscible liquid phases, mole fractions are just dummy values
666 Scalar xwn = 0;
667 Scalar xww = 1.-xwn;
668
669 // Not yet filled with real numbers for the NAPL phase
670 Scalar xnw = 0;
671 Scalar xnn = 1.-xnw;
672
673 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
674 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
675 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
676 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
677 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
678 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
679
680 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
681 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
682 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
683 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
684 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
685 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
686
687 fluidState_.setDensity(wPhaseIdx, rhoW);
688 fluidState_.setDensity(gPhaseIdx, rhoG);
689 fluidState_.setDensity(nPhaseIdx, rhoN);
690 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
691 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
692 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
693 }
694 else if (phasePresence == wnPhaseOnly) {
695 // mole fractions of non-existing gas phase are used as switching criteria
696 Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx);
697 Scalar partPressNAPL = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx);
698
699 Scalar xgn = partPressNAPL/pg_;
700 Scalar xgw = partPressH2O/pg_;
701
702 // immiscible liquid phases, mole fractions are just dummy values
703 Scalar xwn = 0;
704 Scalar xww = 1.-xwn;
705
706 Scalar xnw = 0;
707 Scalar xnn = 1.-xnw;
708
709 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
710 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
711 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
712 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
713 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
714 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
715
716 Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
717 Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
718 Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
719 Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx);
720 Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx);
721 Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx);
722
723 fluidState_.setDensity(wPhaseIdx, rhoW);
724 fluidState_.setDensity(gPhaseIdx, rhoG);
725 fluidState_.setDensity(nPhaseIdx, rhoN);
726 fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar);
727 fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar);
728 fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar);
729 }
730 else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
731 }
732
733 for (int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx) {
734 // Mobilities
735 const Scalar mu =
737 phaseIdx);
738 fluidState_.setViscosity(phaseIdx,mu);
739
740 Scalar kr;
741 kr = MaterialLaw::kr(materialParams, phaseIdx,
742 fluidState_.saturation(wPhaseIdx),
743 fluidState_.saturation(nPhaseIdx),
744 fluidState_.saturation(gPhaseIdx));
745 mobility_[phaseIdx] = kr / mu;
746 Valgrind::CheckDefined(mobility_[phaseIdx]);
747 }
748
749 // material dependent parameters for NAPL adsorption
750 bulkDensTimesAdsorpCoeff_ =
751 MaterialLaw::bulkDensTimesAdsorpCoeff(materialParams);
752
753 /* ATTENTION: The conversion to effective diffusion parameters
754 * for the porous media happens at another place!
755 */
756
757 diffusionCoefficient_[gPhaseIdx] = FluidSystem::diffusionCoefficient(fluidState_, gPhaseIdx);
758 diffusionCoefficient_[wPhaseIdx] = FluidSystem::diffusionCoefficient(fluidState_, wPhaseIdx);
759 /* no diffusion in NAPL phase considered at the moment, dummy values */
760 diffusionCoefficient_[nPhaseIdx] = 1.e-10;
761
762
763 // porosity
764 updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
765 EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
766
767 // permeability
768 permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
769 Valgrind::CheckDefined(permeability_);
770
771// fluidState_.setTemperature(temp_);
772 // the enthalpies (internal energies are directly calculated in the fluidstate
773 for (int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx) {
774 Scalar h = FluidSystem::enthalpy(fluidState_, phaseIdx);
775 fluidState_.setEnthalpy(phaseIdx, h);
776
777 }
778 }
779
783 const FluidState &fluidState() const
784 { return fluidState_; }
785
789 const SolidState &solidState() const
790 { return solidState_; }
791
797 Scalar averageMolarMass(int phaseIdx) const
798 { return fluidState_.averageMolarMass(phaseIdx); }
799
806 Scalar saturation(const int phaseIdx) const
807 { return fluidState_.saturation(phaseIdx); }
808
816 Scalar massFraction(const int phaseIdx, const int compIdx) const
817 { return fluidState_.massFraction(phaseIdx, compIdx); }
818
826 Scalar moleFraction(const int phaseIdx, const int compIdx) const
827 { return fluidState_.moleFraction(phaseIdx, compIdx); }
828
835 Scalar density(const int phaseIdx) const
836 { return fluidState_.density(phaseIdx); }
837
844 Scalar molarDensity(const int phaseIdx) const
845 { return fluidState_.molarDensity(phaseIdx); }
846
853 Scalar pressure(const int phaseIdx) const
854 { return fluidState_.pressure(phaseIdx); }
855
863 Scalar temperature() const
864 { return fluidState_.temperature(/*phaseIdx=*/0); }
865
872 Scalar mobility(const int phaseIdx) const
873 { return mobility_[phaseIdx]; }
874
875 Scalar viscosity(const int phaseIdx) const
876 { return fluidState_.viscosity(phaseIdx); }
877
881 Scalar capillaryPressure() const
882 { return fluidState_.capillaryPressure(); }
883
887 Scalar porosity() const
888 { return solidState_.porosity(); }
889
893 Scalar permeability() const
894 { return permeability_; }
895
899 Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
900 {
901 return diffusionCoefficient_[phaseIdx];
902 }
903
908 { return bulkDensTimesAdsorpCoeff_; }
909
916 Scalar internalEnergy(int phaseIdx) const
917 { return fluidState_.internalEnergy(phaseIdx); };
918
925 Scalar enthalpy(int phaseIdx) const
926 { return fluidState_.enthalpy(phaseIdx); };
927
928protected:
931
932private:
933 Scalar sw_, sg_, sn_, pg_, pw_, pn_, temp_;
934
935 Scalar moleFrac_[numPs][numFluidComps];
936 Scalar massFrac_[numPs][numFluidComps];
937
938 Scalar permeability_;
939 Scalar mobility_[numPs];
940 Scalar bulkDensTimesAdsorpCoeff_;
941 /* We need a tensor here !! */
943 Dune::FieldVector<Scalar, numPs> diffusionCoefficient_;
944 std::array<std::array<Scalar, numFluidComps-1>, numPs> diffCoefficient_;
945
946};
947} // end namespace Dumux
948
949#endif
Some exceptions thrown in DuMux
Define some often used mathematical functions.
Some templates to wrap the valgrind macros.
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
bool CheckDefined(const T &value)
Make valgrind complain if the object occupied by an object is undefined.
Definition: valgrind.hh:72
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Property tag Indices
Definition: porousmediumflow/sequential/properties.hh:59
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
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:56
const FluidState & fluidState() const
Returns the phase state for the control-volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:783
Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
Returns the diffusion coefficient.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:899
SolidState solidState_
Definition: porousmediumflow/3pwateroil/volumevariables.hh:930
Scalar molarDensity(const int phaseIdx) const
Returns the molar density of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:844
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:816
typename Traits::FluidState FluidState
The type of the object returned by the fluidState() method.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:91
Scalar permeability() const
Returns the permeability within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:893
typename Traits::FluidSystem FluidSystem
The type of the fluid system.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:93
Scalar viscosity(const int phaseIdx) const
Definition: porousmediumflow/3pwateroil/volumevariables.hh:875
Scalar bulkDensTimesAdsorpCoeff() const
Returns the adsorption information.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:907
typename Traits::SolidState SolidState
Export type of solid state.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:97
Scalar porosity() const
Returns the average porosity within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:887
Scalar enthalpy(int phaseIdx) const
Returns the total enthalpy of a phase in the sub-control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:925
Scalar mobility(const int phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:872
Scalar capillaryPressure() const
Returns the effective capillary pressure within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:881
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:826
Scalar saturation(const int phaseIdx) const
Returns the effective saturation of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:806
Scalar density(const int phaseIdx) const
Returns the mass density of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:835
typename Traits::SolidSystem SolidSystem
Export type of solid system.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:99
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:115
static constexpr bool onlyGasPhaseCanDisappear()
State if only the gas phase is allowed to disappear.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:103
Scalar temperature() const
Returns temperature inside the sub-control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:863
Scalar internalEnergy(int phaseIdx) const
Returns the total internal energy of a phase in the sub-control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:916
typename ModelTraits::Indices Indices
Export the indices.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:95
const SolidState & solidState() const
Returns the phase state for the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:789
Scalar averageMolarMass(int phaseIdx) const
Returns the average molar mass of the fluid phase.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:797
FluidState fluidState_
Definition: porousmediumflow/3pwateroil/volumevariables.hh:926
Scalar pressure(const int phaseIdx) const
Returns the effective pressure of a given phase within the control volume.
Definition: porousmediumflow/3pwateroil/volumevariables.hh:853
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.