116 const Problem &problem,
117 const Element &element,
121 const auto&
priVars = elemSol[scv.localDofIndex()];
122 const auto phasePresence =
priVars.state();
125 using MaterialLaw =
typename Problem::SpatialParams::MaterialLaw;
126 const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
131 if (phasePresence == threePhases)
135 sg_ = 1. - sw_ - sn_;
137 else if (phasePresence == wPhaseOnly)
143 else if (phasePresence == gnPhaseOnly)
149 else if (phasePresence == wnPhaseOnly)
155 else if (phasePresence == gPhaseOnly)
161 else if (phasePresence == wgPhaseOnly)
167 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
175 if (phasePresence == threePhases || phasePresence == gnPhaseOnly || phasePresence == gPhaseOnly || phasePresence == wgPhaseOnly)
180 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
181 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
182 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
184 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
187 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
188 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
190 else if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly)
195 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
196 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
197 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
199 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
202 pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
203 pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
205 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
213 if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly || phasePresence == gPhaseOnly)
217 else if (phasePresence == threePhases)
220 Scalar temp = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, wCompIdx);
221 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
226 Scalar defect = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
227 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
230 while(abs(defect) > 0.01)
232 Scalar deltaT = 1.e-8 * temp;
233 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
238 Scalar fUp = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
239 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
241 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
246 Scalar fDown = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
247 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
249 temp = temp - defect * 2. * deltaT / (fUp - fDown);
251 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
256 defect = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
257 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
261 else if (phasePresence == wgPhaseOnly)
264 temp_ = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, wCompIdx);
266 else if (phasePresence == gnPhaseOnly)
269 temp_ = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, nCompIdx);
271 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
274 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
281 if (phasePresence == threePhases) {
285 Scalar partPressH2O = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx);
286 Scalar partPressNAPL = pg_ - partPressH2O;
288 Scalar xgn = partPressNAPL/pg_;
289 Scalar xgw = partPressH2O/pg_;
292 Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(
fluidState_, wPhaseIdx,nCompIdx);
296 Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(
fluidState_, nPhaseIdx,wCompIdx);
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);
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);
320 else if (phasePresence == wPhaseOnly) {
325 Scalar xwn =
priVars[switch2Idx];
326 Scalar xww = 1 - xwn;
329 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
330 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
334 Scalar xgn = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx) / pg_;
335 Scalar xgw = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx) / pg_;
341 Scalar xnn = xwn * FluidSystem::henryCoefficient(
fluidState_, wPhaseIdx,nCompIdx) / (xgn * pg_);
342 Scalar xnw = xgw*pg_ / FluidSystem::henryCoefficient(
fluidState_, nPhaseIdx,wCompIdx);
344 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
345 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
346 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
347 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
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);
363 else if (phasePresence == gnPhaseOnly) {
367 Scalar xnw =
priVars[switch2Idx];
369 Scalar xgn = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx) / pg_;
374 Scalar xww = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx) / pg_;
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);
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);
398 else if (phasePresence == wnPhaseOnly) {
401 Scalar partPressH2O = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx);
402 Scalar partPressNAPL = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
404 Scalar xgn = partPressNAPL/pg_;
405 Scalar xgw = partPressH2O/pg_;
408 Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(
fluidState_, wPhaseIdx,nCompIdx);
411 Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(
fluidState_, nPhaseIdx,wCompIdx);
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);
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);
435 else if (phasePresence == gPhaseOnly) {
439 const Scalar xgn =
priVars[switch2Idx];
440 Scalar xgw = 1 - xgn;
443 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
444 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
448 Scalar xww = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx) / pg_;
449 Scalar xnn = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx) / pg_;
451 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
452 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
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);
468 else if (phasePresence == wgPhaseOnly) {
470 const Scalar xgn =
priVars[switch2Idx];
471 Scalar xgw = 1 - xgn;
474 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
475 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
478 Scalar xwn = xgn*pg_/FluidSystem::henryCoefficient(
fluidState_, wPhaseIdx,nCompIdx);
482 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
483 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
487 Scalar xnn = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx) / pg_;
489 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
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);
511 if (phasePresence == threePhases)
515 sg_ = 1. - sw_ - sn_;
517 else if (phasePresence == wnPhaseOnly)
523 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
531 if (phasePresence == threePhases)
536 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
537 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
538 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
540 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
543 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
544 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
546 else if (phasePresence == wnPhaseOnly)
551 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
552 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
553 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
555 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
558 pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
559 pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
561 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
569 if (phasePresence == wnPhaseOnly)
573 else if (phasePresence == threePhases)
577 Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, wCompIdx);
578 temp_ = tempOnlyWater;
582 Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, nCompIdx);
583 temp_ = tempOnlyNAPL;
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)
592 fluidState_.setTemperature(phaseIdx, tempOnlyWater);
595 Scalar defect = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
596 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
598 Scalar temp = tempOnlyWater;
601 while(abs(defect) > 0.01)
603 Scalar deltaT = 1.e-6;
604 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
609 Scalar fUp = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
610 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
612 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
617 Scalar fDown = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
618 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
620 temp = temp - defect * 2. * deltaT / (fUp - fDown);
622 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
627 defect = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
628 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
630 if (counter>10)
break;
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);
639 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
642 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
649 if (phasePresence == threePhases) {
653 Scalar partPressH2O = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx);
654 Scalar partPressNAPL = pg_ - partPressH2O;
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;
662 Scalar xgn = partPressNAPL/pg_;
663 Scalar xgw = partPressH2O/pg_;
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);
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);
694 else if (phasePresence == wnPhaseOnly) {
696 Scalar partPressH2O = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx);
697 Scalar partPressNAPL = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
699 Scalar xgn = partPressNAPL/pg_;
700 Scalar xgw = partPressH2O/pg_;
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);
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);
730 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
733 for (
int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx) {
741 kr = MaterialLaw::kr(materialParams, phaseIdx,
745 mobility_[phaseIdx] = kr / mu;
750 bulkDensTimesAdsorpCoeff_ =
751 MaterialLaw::bulkDensTimesAdsorpCoeff(materialParams);
757 diffusionCoefficient_[gPhaseIdx] = FluidSystem::diffusionCoefficient(
fluidState_, gPhaseIdx);
758 diffusionCoefficient_[wPhaseIdx] = FluidSystem::diffusionCoefficient(
fluidState_, wPhaseIdx);
760 diffusionCoefficient_[nPhaseIdx] = 1.e-10;
765 EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv,
solidState_);
768 permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
773 for (
int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx) {
774 Scalar h = FluidSystem::enthalpy(
fluidState_, phaseIdx);