132 const Problem &problem,
133 const Element &element,
137 const auto&
priVars = elemSol[scv.localDofIndex()];
138 const auto phasePresence =
priVars.state();
143 if (phasePresence == threePhases)
147 sg_ = 1. - sw_ - sn_;
149 else if (phasePresence == wPhaseOnly)
155 else if (phasePresence == gnPhaseOnly)
161 else if (phasePresence == wnPhaseOnly)
167 else if (phasePresence == gPhaseOnly)
173 else if (phasePresence == wgPhaseOnly)
179 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
185 const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
188 const Scalar pcgw = fluidMatrixInteraction.pcgw(sw_, sn_);
189 const Scalar pcnw = fluidMatrixInteraction.pcnw(sw_, sn_);
190 const Scalar pcgn = fluidMatrixInteraction.pcgn(sw_, sn_);
192 const Scalar pcAlpha = fluidMatrixInteraction.pcAlpha(sw_, sn_);
193 const Scalar pcNW1 = 0.0;
196 if (phasePresence == threePhases || phasePresence == gnPhaseOnly || phasePresence == gPhaseOnly || phasePresence == wgPhaseOnly)
199 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
200 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
202 else if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly)
205 pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
206 pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
208 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
215 if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly || phasePresence == gPhaseOnly)
219 else if (phasePresence == threePhases)
222 Scalar temp = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, wCompIdx);
223 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
228 Scalar defect = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
229 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
232 while(abs(defect) > 0.01)
234 Scalar deltaT = 1.e-8 * temp;
235 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
240 Scalar fUp = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
241 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
243 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
248 Scalar fDown = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
249 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
251 temp = temp - defect * 2. * deltaT / (fUp - fDown);
253 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
258 defect = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
259 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
263 else if (phasePresence == wgPhaseOnly)
266 temp_ = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, wCompIdx);
268 else if (phasePresence == gnPhaseOnly)
271 temp_ = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, nCompIdx);
273 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
275 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
282 if (phasePresence == threePhases) {
286 Scalar partPressH2O = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx);
287 Scalar partPressNAPL = pg_ - partPressH2O;
289 Scalar xgn = partPressNAPL/pg_;
290 Scalar xgw = partPressH2O/pg_;
293 Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(
fluidState_, wPhaseIdx,nCompIdx);
297 Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(
fluidState_, nPhaseIdx,wCompIdx);
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);
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);
321 else if (phasePresence == wPhaseOnly) {
326 Scalar xwn =
priVars[switch2Idx];
327 Scalar xww = 1 - xwn;
330 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
331 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
335 Scalar xgn = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx) / pg_;
336 Scalar xgw = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx) / pg_;
342 Scalar xnn = xwn * FluidSystem::henryCoefficient(
fluidState_, wPhaseIdx,nCompIdx) / (xgn * pg_);
343 Scalar xnw = xgw*pg_ / FluidSystem::henryCoefficient(
fluidState_, nPhaseIdx,wCompIdx);
345 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
346 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
347 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
348 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
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);
364 else if (phasePresence == gnPhaseOnly) {
368 Scalar xnw =
priVars[switch2Idx];
370 Scalar xgn = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx) / pg_;
375 Scalar xww = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx) / pg_;
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);
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);
399 else if (phasePresence == wnPhaseOnly) {
402 Scalar partPressH2O = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx);
403 Scalar partPressNAPL = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
405 Scalar xgn = partPressNAPL/pg_;
406 Scalar xgw = partPressH2O/pg_;
409 Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(
fluidState_, wPhaseIdx,nCompIdx);
412 Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(
fluidState_, nPhaseIdx,wCompIdx);
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);
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);
436 else if (phasePresence == gPhaseOnly) {
440 const Scalar xgn =
priVars[switch2Idx];
441 Scalar xgw = 1 - xgn;
444 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
445 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
449 Scalar xww = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx) / pg_;
450 Scalar xnn = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx) / pg_;
452 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
453 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
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);
469 else if (phasePresence == wgPhaseOnly) {
471 const Scalar xgn =
priVars[switch2Idx];
472 Scalar xgw = 1 - xgn;
475 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
476 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
479 Scalar xwn = xgn*pg_/FluidSystem::henryCoefficient(
fluidState_, wPhaseIdx,nCompIdx);
483 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
484 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
488 Scalar xnn = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx) / pg_;
490 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
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);
513 if (phasePresence == threePhases)
517 sg_ = 1. - sw_ - sn_;
519 else if (phasePresence == wnPhaseOnly)
525 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
531 const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
534 const Scalar pcgw = fluidMatrixInteraction.pcgw(sw_, sn_);
535 const Scalar pcnw = fluidMatrixInteraction.pcnw(sw_, sn_);
536 const Scalar pcgn = fluidMatrixInteraction.pcgn(sw_, sn_);
538 const Scalar pcAlpha = fluidMatrixInteraction.pcAlpha(sw_, sn_);
539 const Scalar pcNW1 = 0.0;
542 if (phasePresence == threePhases)
545 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
546 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
548 else if (phasePresence == wnPhaseOnly)
551 pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
552 pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
554 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
561 if (phasePresence == wnPhaseOnly)
565 else if (phasePresence == threePhases)
569 Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, wCompIdx);
570 temp_ = tempOnlyWater;
574 Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, nCompIdx);
575 temp_ = tempOnlyNAPL;
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)
584 fluidState_.setTemperature(phaseIdx, tempOnlyWater);
587 Scalar defect = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
588 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
590 Scalar temp = tempOnlyWater;
593 while(abs(defect) > 0.01)
595 Scalar deltaT = 1.e-6;
596 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
601 Scalar fUp = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
602 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
604 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
609 Scalar fDown = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
610 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
612 temp = temp - defect * 2. * deltaT / (fUp - fDown);
614 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
619 defect = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
620 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
622 if (counter>10)
break;
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);
631 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
633 for (
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
639 if (phasePresence == threePhases) {
643 Scalar partPressH2O = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx);
644 Scalar partPressNAPL = pg_ - partPressH2O;
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;
652 Scalar xgn = partPressNAPL/pg_;
653 Scalar xgw = partPressH2O/pg_;
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);
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);
684 else if (phasePresence == wnPhaseOnly) {
686 Scalar partPressH2O = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx);
687 Scalar partPressNAPL = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
689 Scalar xgn = partPressNAPL/pg_;
690 Scalar xgw = partPressH2O/pg_;
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);
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);
720 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
723 const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
724 for (
int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx)
732 const Scalar kr = fluidMatrixInteraction.kr(phaseIdx,
735 mobility_[phaseIdx] = kr / mu;
740 bulkDensTimesAdsorpCoeff_ = fluidMatrixInteraction.adsorptionModel().bulkDensTimesAdsorpCoeff();
744 EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv,
solidState_);
746 auto getEffectiveDiffusionCoefficient = [&](
int phaseIdx,
int compIIdx,
int compJIdx)
748 return EffDiffModel::effectiveDiffusionCoefficient(*
this, phaseIdx, compIIdx, compJIdx);
751 effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
754 permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
758 for (
int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx)
760 Scalar h = FluidSystem::enthalpy(
fluidState_, phaseIdx);
764 EnergyVolVars::updateEffectiveThermalConductivity();