| Line | Branch | Exec | Source |
|---|---|---|---|
| 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 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup ThreePWaterOilModel | ||
| 10 | * \brief Contains the quantities which are constant within a | ||
| 11 | * finite volume in the three-phase, two-component model. | ||
| 12 | */ | ||
| 13 | |||
| 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> | ||
| 21 | #include <dumux/porousmediumflow/volumevariables.hh> | ||
| 22 | #include <dumux/porousmediumflow/nonisothermal/volumevariables.hh> | ||
| 23 | |||
| 24 | #include <dumux/material/constants.hh> | ||
| 25 | #include <dumux/material/fluidstates/compositional.hh> | ||
| 26 | #include <dumux/material/solidstates/updatesolidvolumefractions.hh> | ||
| 27 | |||
| 28 | #include <dumux/common/optionalscalar.hh> | ||
| 29 | #include <dumux/common/exceptions.hh> | ||
| 30 | |||
| 31 | #include "primaryvariableswitch.hh" | ||
| 32 | |||
| 33 | namespace Dumux { | ||
| 34 | |||
| 35 | namespace Detail { | ||
| 36 | // helper struct and function detecting if the fluid matrix interaction features a adsorptionModel() function | ||
| 37 | #ifndef DOXYGEN // hide from doxygen | ||
| 38 | template <class FluidMatrixInteraction> | ||
| 39 | using AdsorptionModelDetector = decltype(std::declval<FluidMatrixInteraction>().adsorptionModel()); | ||
| 40 | #endif // DOXYGEN | ||
| 41 | |||
| 42 | template<class FluidMatrixInteraction> | ||
| 43 | static constexpr bool hasAdsorptionModel() | ||
| 44 | { return Dune::Std::is_detected<AdsorptionModelDetector, FluidMatrixInteraction>::value; } | ||
| 45 | |||
| 46 | } | ||
| 47 | |||
| 48 | /*! | ||
| 49 | * \ingroup ThreePWaterOilModel | ||
| 50 | * \brief Contains the quantities which are are constant within a | ||
| 51 | * finite volume in the three-phase, two-component model. | ||
| 52 | */ | ||
| 53 | template <class Traits> | ||
| 54 | 2480080 | class ThreePWaterOilVolumeVariables | |
| 55 | : public PorousMediumFlowVolumeVariables<Traits> | ||
| 56 | , public EnergyVolumeVariables<Traits, ThreePWaterOilVolumeVariables<Traits> > | ||
| 57 | { | ||
| 58 | using ParentType = PorousMediumFlowVolumeVariables<Traits>; | ||
| 59 | using EnergyVolVars = EnergyVolumeVariables<Traits, ThreePWaterOilVolumeVariables<Traits> >; | ||
| 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 { | ||
| 66 | numPs = ParentType::numFluidPhases(), | ||
| 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 | |||
| 93 | public: | ||
| 94 | //! The type of the object returned by the fluidState() method | ||
| 95 | using FluidState = typename Traits::FluidState; | ||
| 96 | //! The type of the fluid system | ||
| 97 | using FluidSystem = typename Traits::FluidSystem; | ||
| 98 | //! Export the indices | ||
| 99 | using Indices = typename ModelTraits::Indices; | ||
| 100 | //! Export type of solid state | ||
| 101 | using SolidState = typename Traits::SolidState; | ||
| 102 | //! Export type of solid system | ||
| 103 | using SolidSystem = typename Traits::SolidSystem; | ||
| 104 | //! Export the primary variable switch | ||
| 105 | using PrimaryVariableSwitch = ThreePWaterOilPrimaryVariableSwitch; | ||
| 106 | //! State if only the gas phase is allowed to disappear | ||
| 107 | static constexpr bool onlyGasPhaseCanDisappear() | ||
| 108 | { return Traits::ModelTraits::onlyGasPhaseCanDisappear(); } | ||
| 109 | |||
| 110 | /*! | ||
| 111 | * \brief Updates all quantities for a given control volume. | ||
| 112 | * | ||
| 113 | * \param elemSol A vector containing all primary variables connected to the element | ||
| 114 | * \param problem The object specifying the problem which ought to be simulated | ||
| 115 | * \param element An element which contains part of the control volume | ||
| 116 | * \param scv The sub control volume | ||
| 117 | */ | ||
| 118 | template<class ElemSol, class Problem, class Element, class Scv> | ||
| 119 | 3854948 | void update(const ElemSol &elemSol, | |
| 120 | const Problem &problem, | ||
| 121 | const Element &element, | ||
| 122 | const Scv& scv) | ||
| 123 | { | ||
| 124 | 3854948 | ParentType::update(elemSol, problem, element, scv); | |
| 125 | 3854948 | const auto& priVars = elemSol[scv.localDofIndex()]; | |
| 126 | 3854948 | 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 |
2/2✓ Branch 0 taken 370 times.
✓ Branch 1 taken 3854578 times.
|
3854948 | if (phasePresence == threePhases) |
| 502 | { | ||
| 503 | 370 | sw_ = priVars[switch1Idx]; | |
| 504 | 370 | sn_ = priVars[switch2Idx]; | |
| 505 | 370 | sg_ = 1. - sw_ - sn_; | |
| 506 | } | ||
| 507 |
1/2✓ Branch 0 taken 3854578 times.
✗ Branch 1 not taken.
|
3854578 | else if (phasePresence == wnPhaseOnly) |
| 508 | { | ||
| 509 | 3854578 | sn_ = priVars[switch2Idx]; | |
| 510 | 3854578 | sw_ = 1. - sn_; | |
| 511 | 3854578 | sg_ = 0.; | |
| 512 | } | ||
| 513 | ✗ | else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid."); | |
| 514 | |||
| 515 | 3854948 | fluidState_.setSaturation(wPhaseIdx, sw_); | |
| 516 | 3854948 | fluidState_.setSaturation(gPhaseIdx, sg_); | |
| 517 | 3854948 | fluidState_.setSaturation(nPhaseIdx, sn_); | |
| 518 | |||
| 519 | 7709896 | const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol); | |
| 520 | |||
| 521 | // calculate capillary pressures | ||
| 522 | 3854948 | const Scalar pcgw = fluidMatrixInteraction.pcgw(sw_, sn_); | |
| 523 | 3854948 | const Scalar pcnw = fluidMatrixInteraction.pcnw(sw_, sn_); | |
| 524 | 3854948 | const Scalar pcgn = fluidMatrixInteraction.pcgn(sw_, sn_); | |
| 525 | |||
| 526 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 3854936 times.
|
3854948 | const Scalar pcAlpha = fluidMatrixInteraction.pcAlpha(sw_, sn_); |
| 527 | 3854948 | const Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file | |
| 528 | |||
| 529 | /* now the pressures */ | ||
| 530 |
2/2✓ Branch 0 taken 370 times.
✓ Branch 1 taken 3854578 times.
|
3854948 | if (phasePresence == threePhases) |
| 531 | { | ||
| 532 | 370 | pg_ = priVars[pressureIdx]; | |
| 533 | 370 | pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1); | |
| 534 | 370 | pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1; | |
| 535 | } | ||
| 536 |
1/2✓ Branch 0 taken 3854578 times.
✗ Branch 1 not taken.
|
3854578 | else if (phasePresence == wnPhaseOnly) |
| 537 | { | ||
| 538 | 3854578 | pw_ = priVars[pressureIdx]; | |
| 539 | 3854578 | pn_ = pw_ + pcAlpha * pcnw + (1.-pcAlpha)*pcNW1; | |
| 540 | 3854578 | pg_ = pn_ + pcAlpha * pcgn + (1.-pcAlpha)*(pcgw - pcNW1); | |
| 541 | } | ||
| 542 | ✗ | else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid."); | |
| 543 | |||
| 544 |
2/2✓ Branch 0 taken 3854578 times.
✓ Branch 1 taken 370 times.
|
3854948 | fluidState_.setPressure(wPhaseIdx, pw_); |
| 545 |
2/2✓ Branch 0 taken 3854578 times.
✓ Branch 1 taken 370 times.
|
3854948 | fluidState_.setPressure(gPhaseIdx, pg_); |
| 546 |
2/2✓ Branch 0 taken 3854578 times.
✓ Branch 1 taken 370 times.
|
3854948 | fluidState_.setPressure(nPhaseIdx, pn_); |
| 547 | |||
| 548 | /* now the temperature */ | ||
| 549 |
2/2✓ Branch 0 taken 3854578 times.
✓ Branch 1 taken 370 times.
|
3854948 | if (phasePresence == wnPhaseOnly) |
| 550 | { | ||
| 551 | 7709156 | temp_ = priVars[switch1Idx]; | |
| 552 | } | ||
| 553 |
1/2✓ Branch 0 taken 370 times.
✗ Branch 1 not taken.
|
370 | else if (phasePresence == threePhases) |
| 554 | { | ||
| 555 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 358 times.
|
370 | if(sn_<=1.e-10) // this threshold values is chosen arbitrarily as a small number |
| 556 | { | ||
| 557 | 12 | Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx); | |
| 558 | 12 | temp_ = tempOnlyWater; | |
| 559 | } | ||
| 560 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 368 times.
|
370 | if(sw_<=1.e-10) // this threshold values is chosen arbitrarily as a small number |
| 561 | { | ||
| 562 | 2 | Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx); | |
| 563 | 2 | temp_ = tempOnlyNAPL; | |
| 564 | } | ||
| 565 | else | ||
| 566 | { | ||
| 567 | // temp from inverse pwsat and pnsat which have to sum up to pg | ||
| 568 | 368 | Scalar tempOnlyNAPL = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, nCompIdx); | |
| 569 | 368 | Scalar tempOnlyWater = FluidSystem::inverseVaporPressureCurve(fluidState_, gPhaseIdx, wCompIdx); | |
| 570 |
2/2✓ Branch 1 taken 1104 times.
✓ Branch 2 taken 368 times.
|
1472 | for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) |
| 571 | { | ||
| 572 | 2208 | fluidState_.setTemperature(phaseIdx, tempOnlyWater); | |
| 573 | } | ||
| 574 | 368 | solidState_.setTemperature(tempOnlyWater); | |
| 575 | 368 | Scalar defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) | |
| 576 | 368 | - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx); | |
| 577 | |||
| 578 | 368 | Scalar temp = tempOnlyWater; // initial guess | |
| 579 | 368 | int counter = 0; | |
| 580 | using std::abs; | ||
| 581 |
4/4✓ Branch 0 taken 368 times.
✓ Branch 1 taken 368 times.
✓ Branch 2 taken 368 times.
✓ Branch 3 taken 368 times.
|
1472 | 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 |
2/2✓ Branch 0 taken 1104 times.
✓ Branch 1 taken 368 times.
|
1472 | for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) |
| 585 | { | ||
| 586 | 2208 | fluidState_.setTemperature(phaseIdx, temp+deltaT); | |
| 587 | } | ||
| 588 | 368 | solidState_.setTemperature(temp+deltaT); | |
| 589 | 368 | Scalar fUp = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) | |
| 590 | 368 | - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx); | |
| 591 | |||
| 592 |
2/2✓ Branch 0 taken 1104 times.
✓ Branch 1 taken 368 times.
|
1472 | for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) |
| 593 | { | ||
| 594 | 2208 | fluidState_.setTemperature(phaseIdx, temp-deltaT); | |
| 595 | } | ||
| 596 | 368 | solidState_.setTemperature(temp-deltaT); | |
| 597 | 368 | Scalar fDown = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) | |
| 598 | 368 | - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx); | |
| 599 | |||
| 600 | 368 | temp = temp - defect * 2. * deltaT / (fUp - fDown); | |
| 601 | |||
| 602 |
2/2✓ Branch 0 taken 1104 times.
✓ Branch 1 taken 368 times.
|
1472 | for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) |
| 603 | { | ||
| 604 | 2208 | fluidState_.setTemperature(phaseIdx, temp); | |
| 605 | } | ||
| 606 | 368 | solidState_.setTemperature(temp); | |
| 607 | 368 | defect = pg_ - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx) | |
| 608 | 368 | - FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx); | |
| 609 | 368 | counter +=1; | |
| 610 |
1/2✓ Branch 0 taken 368 times.
✗ Branch 1 not taken.
|
368 | if (counter>10) break; |
| 611 | } | ||
| 612 |
2/4✓ Branch 0 taken 368 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 368 times.
|
368 | if ((sw_>1.e-10)&&(sw_<0.01)) |
| 613 | ✗ | temp = temp + (sw_ - 1.e-10) * (temp - tempOnlyNAPL) / (0.01 - 1.e-10); | |
| 614 |
3/4✓ Branch 0 taken 357 times.
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 357 times.
|
368 | if ((sn_>1.e-10)&&(sn_<0.01)) |
| 615 | ✗ | temp = temp + (sn_ - 1.e-10) * (temp - tempOnlyWater) / (0.01 - 1.e-10); | |
| 616 | 368 | temp_ = temp; | |
| 617 | } | ||
| 618 | } | ||
| 619 | ✗ | else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid."); | |
| 620 | |||
| 621 |
2/2✓ Branch 0 taken 11564844 times.
✓ Branch 1 taken 3854948 times.
|
15419792 | for (int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) |
| 622 | 23129688 | fluidState_.setTemperature(phaseIdx, temp_); | |
| 623 | |||
| 624 |
2/2✓ Branch 0 taken 370 times.
✓ Branch 1 taken 3854578 times.
|
3854948 | solidState_.setTemperature(temp_); |
| 625 | |||
| 626 | // now comes the tricky part: calculate phase composition | ||
| 627 |
2/2✓ Branch 0 taken 370 times.
✓ Branch 1 taken 3854578 times.
|
3854948 | if (phasePresence == threePhases) { |
| 628 | |||
| 629 | // all phases are present, phase compositions are a | ||
| 630 | // result of the the gas <-> liquid equilibrium. | ||
| 631 | 370 | Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx); | |
| 632 | 370 | Scalar partPressNAPL = pg_ - partPressH2O; | |
| 633 | // regularized evaporation for small liquid phase saturations | ||
| 634 | // avoids negative saturations of liquid phases | ||
| 635 |
2/2✓ Branch 0 taken 13 times.
✓ Branch 1 taken 357 times.
|
370 | if (sw_<0.02) partPressH2O *= sw_/0.02; |
| 636 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 368 times.
|
370 | if (partPressH2O < 0.) partPressH2O = 0; |
| 637 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 358 times.
|
370 | if (sn_<0.02) partPressNAPL *= sn_ / 0.02; |
| 638 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 358 times.
|
370 | if (partPressNAPL < 0.) partPressNAPL = 0; |
| 639 | |||
| 640 | 370 | Scalar xgn = partPressNAPL/pg_; | |
| 641 | 370 | Scalar xgw = partPressH2O/pg_; | |
| 642 | |||
| 643 | // Immiscible liquid phases, mole fractions are just dummy values | ||
| 644 | 370 | Scalar xwn = 0; | |
| 645 | 370 | Scalar xww = 1.-xwn; | |
| 646 | |||
| 647 | // Not yet filled with real numbers for the NAPL phase | ||
| 648 | 370 | Scalar xnw = 0; | |
| 649 | 370 | Scalar xnn = 1.-xnw; | |
| 650 | |||
| 651 | 370 | fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww); | |
| 652 | 370 | fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn); | |
| 653 | 370 | fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw); | |
| 654 | 370 | fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn); | |
| 655 | 370 | fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw); | |
| 656 | 370 | fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn); | |
| 657 | |||
| 658 | 370 | Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx); | |
| 659 | 368 | Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx); | |
| 660 | 368 | Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx); | |
| 661 | 368 | Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx); | |
| 662 | 368 | Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx); | |
| 663 | 368 | Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx); | |
| 664 | |||
| 665 | 368 | fluidState_.setDensity(wPhaseIdx, rhoW); | |
| 666 | 368 | fluidState_.setDensity(gPhaseIdx, rhoG); | |
| 667 | 368 | fluidState_.setDensity(nPhaseIdx, rhoN); | |
| 668 | 368 | fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar); | |
| 669 | 368 | fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar); | |
| 670 | 368 | fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar); | |
| 671 | } | ||
| 672 |
1/2✓ Branch 0 taken 3854578 times.
✗ Branch 1 not taken.
|
3854578 | else if (phasePresence == wnPhaseOnly) { |
| 673 | // mole fractions of non-existing gas phase are used as switching criteria | ||
| 674 | 3854578 | Scalar partPressH2O = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, wCompIdx); | |
| 675 | 3854578 | Scalar partPressNAPL = FluidSystem::partialPressureGas(fluidState_, gPhaseIdx, nCompIdx); | |
| 676 | |||
| 677 | 3854578 | Scalar xgn = partPressNAPL/pg_; | |
| 678 | 3854578 | Scalar xgw = partPressH2O/pg_; | |
| 679 | |||
| 680 | // immiscible liquid phases, mole fractions are just dummy values | ||
| 681 | 3854578 | Scalar xwn = 0; | |
| 682 | 3854578 | Scalar xww = 1.-xwn; | |
| 683 | |||
| 684 | 3854578 | Scalar xnw = 0; | |
| 685 | 3854578 | Scalar xnn = 1.-xnw; | |
| 686 | |||
| 687 | 3854578 | fluidState_.setMoleFraction(wPhaseIdx, wCompIdx, xww); | |
| 688 | 3854578 | fluidState_.setMoleFraction(wPhaseIdx, nCompIdx, xwn); | |
| 689 | 3854578 | fluidState_.setMoleFraction(gPhaseIdx, wCompIdx, xgw); | |
| 690 | 3854578 | fluidState_.setMoleFraction(gPhaseIdx, nCompIdx, xgn); | |
| 691 | 3854578 | fluidState_.setMoleFraction(nPhaseIdx, wCompIdx, xnw); | |
| 692 | 3854578 | fluidState_.setMoleFraction(nPhaseIdx, nCompIdx, xnn); | |
| 693 | |||
| 694 | 3854578 | Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx); | |
| 695 | 3854578 | Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx); | |
| 696 | 3854578 | Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx); | |
| 697 | 3854578 | Scalar rhoWMolar = FluidSystem::molarDensity(fluidState_, wPhaseIdx); | |
| 698 | 3854578 | Scalar rhoGMolar = FluidSystem::molarDensity(fluidState_, gPhaseIdx); | |
| 699 | 3854578 | Scalar rhoNMolar = FluidSystem::molarDensity(fluidState_, nPhaseIdx); | |
| 700 | |||
| 701 | 3854578 | fluidState_.setDensity(wPhaseIdx, rhoW); | |
| 702 | 3854578 | fluidState_.setDensity(gPhaseIdx, rhoG); | |
| 703 | 3854578 | fluidState_.setDensity(nPhaseIdx, rhoN); | |
| 704 | 3854578 | fluidState_.setMolarDensity(wPhaseIdx, rhoWMolar); | |
| 705 | 3854578 | fluidState_.setMolarDensity(gPhaseIdx, rhoGMolar); | |
| 706 | 3854578 | fluidState_.setMolarDensity(nPhaseIdx, rhoNMolar); | |
| 707 | } | ||
| 708 | ✗ | else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid."); | |
| 709 | } | ||
| 710 | |||
| 711 | 7709892 | const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol); | |
| 712 |
2/2✓ Branch 0 taken 11564838 times.
✓ Branch 1 taken 3854946 times.
|
15419784 | for (int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx) |
| 713 | { | ||
| 714 | // Mobilities | ||
| 715 | const Scalar mu = | ||
| 716 | 11564838 | FluidSystem::viscosity(fluidState_, | |
| 717 | phaseIdx); | ||
| 718 | 11564838 | fluidState_.setViscosity(phaseIdx,mu); | |
| 719 | |||
| 720 | 34694514 | const Scalar kr = fluidMatrixInteraction.kr(phaseIdx, | |
| 721 | fluidState_.saturation(wPhaseIdx), | ||
| 722 | fluidState_.saturation(nPhaseIdx)); | ||
| 723 | 11564838 | 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 |
2/2✓ Branch 0 taken 483270 times.
✓ Branch 1 taken 3371676 times.
|
3854946 | updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps); |
| 732 | 3854946 | EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_); | |
| 733 | |||
| 734 | 3854946 | auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx) | |
| 735 | { | ||
| 736 | 11564838 | return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx); | |
| 737 | }; | ||
| 738 | |||
| 739 | 3854946 | effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient); | |
| 740 | |||
| 741 | // permeability | ||
| 742 |
4/4✓ Branch 0 taken 533126 times.
✓ Branch 1 taken 3321820 times.
✓ Branch 2 taken 533126 times.
✓ Branch 3 taken 3321820 times.
|
7709892 | permeability_ = problem.spatialParams().permeability(element, scv, elemSol); |
| 743 | |||
| 744 | 3854946 | fluidState_.setTemperature(temp_); | |
| 745 | // the enthalpies (internal energies are directly calculated in the fluidstate | ||
| 746 |
2/2✓ Branch 0 taken 11564838 times.
✓ Branch 1 taken 3854946 times.
|
15419784 | for (int phaseIdx = 0; phaseIdx < numPs; ++phaseIdx) |
| 747 | { | ||
| 748 | 11564838 | Scalar h = FluidSystem::enthalpy(fluidState_, phaseIdx); | |
| 749 | 23129676 | fluidState_.setEnthalpy(phaseIdx, h); | |
| 750 | } | ||
| 751 | |||
| 752 | 3854946 | EnergyVolVars::updateEffectiveThermalConductivity(); | |
| 753 | 3854946 | } | |
| 754 | |||
| 755 | /*! | ||
| 756 | * \brief Returns the phase state for the control-volume. | ||
| 757 | */ | ||
| 758 | const FluidState &fluidState() const | ||
| 759 |
6/8✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1664 times.
✓ Branch 5 taken 6 times.
✓ Branch 6 taken 153692 times.
✓ Branch 7 taken 6 times.
✓ Branch 8 taken 153692 times.
|
4165674 | { return fluidState_; } |
| 760 | |||
| 761 | /*! | ||
| 762 | * \brief Returns the phase state for the control volume. | ||
| 763 | */ | ||
| 764 | const SolidState &solidState() const | ||
| 765 | 3854946 | { return solidState_; } | |
| 766 | |||
| 767 | /*! | ||
| 768 | * \brief Returns the average molar mass \f$\mathrm{[kg/mol]}\f$ of the fluid phase. | ||
| 769 | * | ||
| 770 | * \param phaseIdx The phase index | ||
| 771 | */ | ||
| 772 | Scalar averageMolarMass(int phaseIdx) const | ||
| 773 | { return fluidState_.averageMolarMass(phaseIdx); } | ||
| 774 | |||
| 775 | /*! | ||
| 776 | * \brief Returns the effective saturation of a given phase within | ||
| 777 | * the control volume. | ||
| 778 | * | ||
| 779 | * \param phaseIdx The phase index | ||
| 780 | */ | ||
| 781 | Scalar saturation(const int phaseIdx) const | ||
| 782 |
4/8✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✓ Branch 9 taken 26 times.
✓ Branch 10 taken 2 times.
✓ Branch 11 taken 26 times.
|
302959052 | { return fluidState_.saturation(phaseIdx); } |
| 783 | |||
| 784 | /*! | ||
| 785 | * \brief Returns the mass fraction of a given component in a | ||
| 786 | * given phase within the control volume in \f$[-]\f$. | ||
| 787 | * | ||
| 788 | * \param phaseIdx The phase index | ||
| 789 | * \param compIdx The component index | ||
| 790 | */ | ||
| 791 | Scalar massFraction(const int phaseIdx, const int compIdx) const | ||
| 792 | 95846400 | { return fluidState_.massFraction(phaseIdx, compIdx); } | |
| 793 | |||
| 794 | /*! | ||
| 795 | * \brief Returns the mole fraction of a given component in a | ||
| 796 | * given phase within the control volume in \f$[-]\f$. | ||
| 797 | * | ||
| 798 | * \param phaseIdx The phase index | ||
| 799 | * \param compIdx The component index | ||
| 800 | */ | ||
| 801 | Scalar moleFraction(const int phaseIdx, const int compIdx) const | ||
| 802 | 383385600 | { return fluidState_.moleFraction(phaseIdx, compIdx); } | |
| 803 | |||
| 804 | /*! | ||
| 805 | * \brief Returns the mass density of a given phase within the | ||
| 806 | * control volume. | ||
| 807 | * | ||
| 808 | * \param phaseIdx The phase index | ||
| 809 | */ | ||
| 810 | Scalar density(const int phaseIdx) const | ||
| 811 | 575078400 | { return fluidState_.density(phaseIdx); } | |
| 812 | |||
| 813 | /*! | ||
| 814 | * \brief Returns the molar density of a given phase within the | ||
| 815 | * control volume. | ||
| 816 | * | ||
| 817 | * \param phaseIdx The phase index | ||
| 818 | */ | ||
| 819 | Scalar molarDensity(const int phaseIdx) const | ||
| 820 | 383385600 | { return fluidState_.molarDensity(phaseIdx); } | |
| 821 | |||
| 822 | /*! | ||
| 823 | * \brief Returns the effective pressure of a given phase within | ||
| 824 | * the control volume. | ||
| 825 | * | ||
| 826 | * \param phaseIdx The phase index | ||
| 827 | */ | ||
| 828 | Scalar pressure(const int phaseIdx) const | ||
| 829 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1664 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1664 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1664 times.
|
191699468 | { return fluidState_.pressure(phaseIdx); } |
| 830 | |||
| 831 | /*! | ||
| 832 | * \brief Returns temperature inside the sub-control volume. | ||
| 833 | * | ||
| 834 | * Note that we assume thermodynamic equilibrium, i.e. the | ||
| 835 | * temperatures of the rock matrix and of all fluid phases are | ||
| 836 | * identical. | ||
| 837 | */ | ||
| 838 | Scalar temperature() const | ||
| 839 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1664 times.
|
63900928 | { return fluidState_.temperature(/*phaseIdx=*/0); } |
| 840 | |||
| 841 | /*! | ||
| 842 | * \brief Returns the effective mobility of a given phase within | ||
| 843 | * the control volume. | ||
| 844 | * | ||
| 845 | * \param phaseIdx The phase index | ||
| 846 | */ | ||
| 847 | Scalar mobility(const int phaseIdx) const | ||
| 848 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1664 times.
|
143772928 | { return mobility_[phaseIdx]; } |
| 849 | |||
| 850 | Scalar viscosity(const int phaseIdx) const | ||
| 851 | 115200 | { return fluidState_.viscosity(phaseIdx); } | |
| 852 | |||
| 853 | /*! | ||
| 854 | * \brief Returns the effective capillary pressure within the control volume. | ||
| 855 | */ | ||
| 856 | Scalar capillaryPressure() const | ||
| 857 | { return fluidState_.capillaryPressure(); } | ||
| 858 | |||
| 859 | /*! | ||
| 860 | * \brief Returns the average porosity within the control volume. | ||
| 861 | */ | ||
| 862 | Scalar porosity() const | ||
| 863 |
0/4✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
295249092 | { return solidState_.porosity(); } |
| 864 | |||
| 865 | /*! | ||
| 866 | * \brief Returns the permeability within the control volume. | ||
| 867 | */ | ||
| 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 | |||
| 882 | /*! | ||
| 883 | * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$. | ||
| 884 | */ | ||
| 885 | Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const | ||
| 886 | 23961600 | { return effectiveDiffCoeff_(phaseIdx, compIIdx, compJIdx); } | |
| 887 | |||
| 888 | /*! | ||
| 889 | * \brief Returns the adsorption information. | ||
| 890 | */ | ||
| 891 | Scalar bulkDensTimesAdsorpCoeff() const | ||
| 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 | |||
| 899 | /*! | ||
| 900 | * \brief Returns the total internal energy of a phase in the | ||
| 901 | * sub-control volume. | ||
| 902 | * | ||
| 903 | * \param phaseIdx The phase index | ||
| 904 | */ | ||
| 905 | Scalar internalEnergy(int phaseIdx) const | ||
| 906 | 95846400 | { return fluidState_.internalEnergy(phaseIdx); }; | |
| 907 | |||
| 908 | /*! | ||
| 909 | * \brief Returns the total enthalpy of a phase in the sub-control | ||
| 910 | * volume. | ||
| 911 | * | ||
| 912 | * \param phaseIdx The phase index | ||
| 913 | */ | ||
| 914 | Scalar enthalpy(int phaseIdx) const | ||
| 915 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1664 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1664 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1664 times.
|
95853056 | { return fluidState_.enthalpy(phaseIdx); }; |
| 916 | |||
| 917 | protected: | ||
| 918 | FluidState fluidState_; | ||
| 919 | SolidState solidState_; | ||
| 920 | |||
| 921 | private: | ||
| 922 | Scalar sw_, sg_, sn_, pg_, pw_, pn_, temp_; | ||
| 923 | |||
| 924 | Scalar moleFrac_[numPs][numFluidComps]; | ||
| 925 | Scalar massFrac_[numPs][numFluidComps]; | ||
| 926 | |||
| 927 | Scalar permeability_; //!< Effective porosity within the control volume | ||
| 928 | Scalar mobility_[numPs]; //!< Effective mobility within the control volume | ||
| 929 | OptionalScalar<Scalar> bulkDensTimesAdsorpCoeff_; //!< the basis for calculating adsorbed NAPL | ||
| 930 | |||
| 931 | //!< Binary diffusion coefficients of the 3 components in the phases | ||
| 932 | DiffusionCoefficients effectiveDiffCoeff_; | ||
| 933 | |||
| 934 | }; | ||
| 935 | } // end namespace Dumux | ||
| 936 | |||
| 937 | #endif | ||
| 938 |