119 const Problem &problem,
120 const Element &element,
124 const auto&
priVars = elemSol[scv.localDofIndex()];
125 const auto phasePresence =
priVars.state();
128 using MaterialLaw =
typename Problem::SpatialParams::MaterialLaw;
129 const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
134 if (phasePresence == threePhases)
138 sg_ = 1. - sw_ - sn_;
140 else if (phasePresence == wPhaseOnly)
146 else if (phasePresence == gnPhaseOnly)
152 else if (phasePresence == wnPhaseOnly)
158 else if (phasePresence == gPhaseOnly)
164 else if (phasePresence == wgPhaseOnly)
170 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
178 if (phasePresence == threePhases || phasePresence == gnPhaseOnly || phasePresence == gPhaseOnly || phasePresence == wgPhaseOnly)
183 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
184 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
185 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
187 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
190 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
191 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
193 else if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly)
198 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
199 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
200 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
202 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
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.");
216 if (phasePresence == wPhaseOnly || phasePresence == wnPhaseOnly || phasePresence == gPhaseOnly)
220 else if (phasePresence == threePhases)
223 Scalar temp = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, wCompIdx);
224 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
229 Scalar defect = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
230 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
233 while(abs(defect) > 0.01)
235 Scalar deltaT = 1.e-8 * temp;
236 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
241 Scalar fUp = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
242 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
244 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
249 Scalar fDown = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
250 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
252 temp = temp - defect * 2. * deltaT / (fUp - fDown);
254 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
259 defect = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
260 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
264 else if (phasePresence == wgPhaseOnly)
267 temp_ = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, wCompIdx);
269 else if (phasePresence == gnPhaseOnly)
272 temp_ = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, nCompIdx);
274 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
277 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
284 if (phasePresence == threePhases) {
288 Scalar partPressH2O = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx);
289 Scalar partPressNAPL = pg_ - partPressH2O;
291 Scalar xgn = partPressNAPL/pg_;
292 Scalar xgw = partPressH2O/pg_;
295 Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(
fluidState_, wPhaseIdx,nCompIdx);
299 Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(
fluidState_, nPhaseIdx,wCompIdx);
302 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
303 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
304 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
305 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
306 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
307 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
309 Scalar rhoW = FluidSystem::density(
fluidState_, wPhaseIdx);
310 Scalar rhoG = FluidSystem::density(
fluidState_, gPhaseIdx);
311 Scalar rhoN = FluidSystem::density(
fluidState_, nPhaseIdx);
312 Scalar rhoWMolar = FluidSystem::molarDensity(
fluidState_, wPhaseIdx);
313 Scalar rhoGMolar = FluidSystem::molarDensity(
fluidState_, gPhaseIdx);
314 Scalar rhoNMolar = FluidSystem::molarDensity(
fluidState_, nPhaseIdx);
323 else if (phasePresence == wPhaseOnly) {
328 Scalar xwn =
priVars[switch2Idx];
329 Scalar xww = 1 - xwn;
332 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
333 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
337 Scalar xgn = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx) / pg_;
338 Scalar xgw = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx) / pg_;
344 Scalar xnn = xwn * FluidSystem::henryCoefficient(
fluidState_, wPhaseIdx,nCompIdx) / (xgn * pg_);
345 Scalar xnw = xgw*pg_ / FluidSystem::henryCoefficient(
fluidState_, nPhaseIdx,wCompIdx);
347 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
348 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
349 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
350 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
352 Scalar rhoW = FluidSystem::density(
fluidState_, wPhaseIdx);
353 Scalar rhoG = FluidSystem::density(
fluidState_, gPhaseIdx);
354 Scalar rhoN = FluidSystem::density(
fluidState_, nPhaseIdx);
355 Scalar rhoWMolar = FluidSystem::molarDensity(
fluidState_, wPhaseIdx);
356 Scalar rhoGMolar = FluidSystem::molarDensity(
fluidState_, gPhaseIdx);
357 Scalar rhoNMolar = FluidSystem::molarDensity(
fluidState_, nPhaseIdx);
366 else if (phasePresence == gnPhaseOnly) {
370 Scalar xnw =
priVars[switch2Idx];
372 Scalar xgn = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx) / pg_;
377 Scalar xww = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx) / pg_;
380 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
381 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
382 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
383 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
384 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
386 Scalar rhoW = FluidSystem::density(
fluidState_, wPhaseIdx);
387 Scalar rhoG = FluidSystem::density(
fluidState_, gPhaseIdx);
388 Scalar rhoN = FluidSystem::density(
fluidState_, nPhaseIdx);
389 Scalar rhoWMolar = FluidSystem::molarDensity(
fluidState_, wPhaseIdx);
390 Scalar rhoGMolar = FluidSystem::molarDensity(
fluidState_, gPhaseIdx);
391 Scalar rhoNMolar = FluidSystem::molarDensity(
fluidState_, nPhaseIdx);
401 else if (phasePresence == wnPhaseOnly) {
404 Scalar partPressH2O = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx);
405 Scalar partPressNAPL = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
407 Scalar xgn = partPressNAPL/pg_;
408 Scalar xgw = partPressH2O/pg_;
411 Scalar xwn = partPressNAPL / FluidSystem::henryCoefficient(
fluidState_, wPhaseIdx,nCompIdx);
414 Scalar xnw = partPressH2O / FluidSystem::henryCoefficient(
fluidState_, nPhaseIdx,wCompIdx);
417 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
418 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
419 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
420 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
421 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
422 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
424 Scalar rhoW = FluidSystem::density(
fluidState_, wPhaseIdx);
425 Scalar rhoG = FluidSystem::density(
fluidState_, gPhaseIdx);
426 Scalar rhoN = FluidSystem::density(
fluidState_, nPhaseIdx);
427 Scalar rhoWMolar = FluidSystem::molarDensity(
fluidState_, wPhaseIdx);
428 Scalar rhoGMolar = FluidSystem::molarDensity(
fluidState_, gPhaseIdx);
429 Scalar rhoNMolar = FluidSystem::molarDensity(
fluidState_, nPhaseIdx);
438 else if (phasePresence == gPhaseOnly) {
442 const Scalar xgn =
priVars[switch2Idx];
443 Scalar xgw = 1 - xgn;
446 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
447 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
451 Scalar xww = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx) / pg_;
452 Scalar xnn = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx) / pg_;
454 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
455 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
457 Scalar rhoW = FluidSystem::density(
fluidState_, wPhaseIdx);
458 Scalar rhoG = FluidSystem::density(
fluidState_, gPhaseIdx);
459 Scalar rhoN = FluidSystem::density(
fluidState_, nPhaseIdx);
460 Scalar rhoWMolar = FluidSystem::molarDensity(
fluidState_, wPhaseIdx);
461 Scalar rhoGMolar = FluidSystem::molarDensity(
fluidState_, gPhaseIdx);
462 Scalar rhoNMolar = FluidSystem::molarDensity(
fluidState_, nPhaseIdx);
471 else if (phasePresence == wgPhaseOnly) {
473 const Scalar xgn =
priVars[switch2Idx];
474 Scalar xgw = 1 - xgn;
477 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
478 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
481 Scalar xwn = xgn*pg_/FluidSystem::henryCoefficient(
fluidState_, wPhaseIdx,nCompIdx);
485 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
486 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
490 Scalar xnn = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx) / pg_;
492 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
494 Scalar rhoW = FluidSystem::density(
fluidState_, wPhaseIdx);
495 Scalar rhoG = FluidSystem::density(
fluidState_, gPhaseIdx);
496 Scalar rhoN = FluidSystem::density(
fluidState_, nPhaseIdx);
497 Scalar rhoWMolar = FluidSystem::molarDensity(
fluidState_, wPhaseIdx);
498 Scalar rhoGMolar = FluidSystem::molarDensity(
fluidState_, gPhaseIdx);
499 Scalar rhoNMolar = FluidSystem::molarDensity(
fluidState_, nPhaseIdx);
514 if (phasePresence == threePhases)
518 sg_ = 1. - sw_ - sn_;
520 else if (phasePresence == wnPhaseOnly)
526 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
534 if (phasePresence == threePhases)
539 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
540 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
541 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
543 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
546 pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
547 pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
549 else if (phasePresence == wnPhaseOnly)
554 Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
555 Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
556 Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
558 Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
561 pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1;
562 pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1);
564 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
572 if (phasePresence == wnPhaseOnly)
576 else if (phasePresence == threePhases)
580 Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, wCompIdx);
581 temp_ = tempOnlyWater;
585 Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, nCompIdx);
586 temp_ = tempOnlyNAPL;
591 Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, nCompIdx);
592 Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(
fluidState_, gPhaseIdx, wCompIdx);
593 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
595 fluidState_.setTemperature(phaseIdx, tempOnlyWater);
598 Scalar defect = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
599 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
601 Scalar temp = tempOnlyWater;
604 while(abs(defect) > 0.01)
606 Scalar deltaT = 1.e-6;
607 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
612 Scalar fUp = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
613 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
615 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
620 Scalar fDown = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
621 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
623 temp = temp - defect * 2. * deltaT / (fUp - fDown);
625 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
630 defect = pg_ - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx)
631 - FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
633 if (counter>10)
break;
635 if ((sw_>1.e-10)&&(sw_<0.01))
636 temp = temp + (sw_ - 1.e-10) * (temp - tempOnlyNAPL) / (0.01 - 1.e-10);
637 if ((sn_>1.e-10)&&(sn_<0.01))
638 temp = temp + (sn_ - 1.e-10) * (temp - tempOnlyWater) / (0.01 - 1.e-10);
642 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
645 for(
int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
652 if (phasePresence == threePhases) {
656 Scalar partPressH2O = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx);
657 Scalar partPressNAPL = pg_ - partPressH2O;
660 if (sw_<0.02) partPressH2O *= sw_/0.02;
661 if (partPressH2O < 0.) partPressH2O = 0;
662 if (sn_<0.02) partPressNAPL *= sn_ / 0.02;
663 if (partPressNAPL < 0.) partPressNAPL = 0;
665 Scalar xgn = partPressNAPL/pg_;
666 Scalar xgw = partPressH2O/pg_;
676 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
677 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
678 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
679 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
680 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
681 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
683 Scalar rhoW = FluidSystem::density(
fluidState_, wPhaseIdx);
684 Scalar rhoG = FluidSystem::density(
fluidState_, gPhaseIdx);
685 Scalar rhoN = FluidSystem::density(
fluidState_, nPhaseIdx);
686 Scalar rhoWMolar = FluidSystem::molarDensity(
fluidState_, wPhaseIdx);
687 Scalar rhoGMolar = FluidSystem::molarDensity(
fluidState_, gPhaseIdx);
688 Scalar rhoNMolar = FluidSystem::molarDensity(
fluidState_, nPhaseIdx);
697 else if (phasePresence == wnPhaseOnly) {
699 Scalar partPressH2O = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, wCompIdx);
700 Scalar partPressNAPL = FluidSystem::partialPressureGas(
fluidState_, gPhaseIdx, nCompIdx);
702 Scalar xgn = partPressNAPL/pg_;
703 Scalar xgw = partPressH2O/pg_;
712 fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww);
713 fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn);
714 fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw);
715 fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn);
716 fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw);
717 fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn);
719 Scalar rhoW = FluidSystem::density(
fluidState_, wPhaseIdx);
720 Scalar rhoG = FluidSystem::density(
fluidState_, gPhaseIdx);
721 Scalar rhoN = FluidSystem::density(
fluidState_, nPhaseIdx);
722 Scalar rhoWMolar = FluidSystem::molarDensity(
fluidState_, wPhaseIdx);
723 Scalar rhoGMolar = FluidSystem::molarDensity(
fluidState_, gPhaseIdx);
724 Scalar rhoNMolar = FluidSystem::molarDensity(
fluidState_, nPhaseIdx);
733 else DUNE_THROW(Dune::InvalidStateException,
"phasePresence: " << phasePresence <<
" is invalid.");
736 for (
int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx) {
744 kr = MaterialLaw::kr(materialParams, phaseIdx,
748 mobility_[phaseIdx] = kr / mu;
753 bulkDensTimesAdsorpCoeff_ = MaterialLaw::bulkDensTimesAdsorpCoeff(materialParams);
757 EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv,
solidState_);
759 auto getEffectiveDiffusionCoefficient = [&](
int phaseIdx,
int compIIdx,
int compJIdx)
761 return EffDiffModel::effectiveDiffusionCoefficient(*
this, phaseIdx, compIIdx, compJIdx);
764 effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient);
767 permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
772 for (
int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx)
774 Scalar h = FluidSystem::enthalpy(
fluidState_, phaseIdx);
778 EnergyVolVars::updateEffectiveThermalConductivity();