| 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 ExtendedRichardsModel | ||
| 10 | * \brief Volume averaged quantities required by the extended Richards model. | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_RICHARDSEXTENDED_VOLUME_VARIABLES_HH | ||
| 14 | #define DUMUX_RICHARDSEXTENDED_VOLUME_VARIABLES_HH | ||
| 15 | |||
| 16 | #include <cassert> | ||
| 17 | |||
| 18 | #include <dune/common/exceptions.hh> | ||
| 19 | |||
| 20 | #include <dumux/porousmediumflow/volumevariables.hh> | ||
| 21 | #include <dumux/porousmediumflow/nonisothermal/volumevariables.hh> | ||
| 22 | #include <dumux/material/idealgas.hh> | ||
| 23 | #include <dumux/material/constants.hh> | ||
| 24 | #include <dumux/material/solidstates/updatesolidvolumefractions.hh> | ||
| 25 | |||
| 26 | #include "primaryvariableswitch.hh" | ||
| 27 | |||
| 28 | namespace Dumux { | ||
| 29 | |||
| 30 | /*! | ||
| 31 | * \ingroup ExtendedRichardsModel | ||
| 32 | * \brief Volume averaged quantities required by the extended Richards model. | ||
| 33 | * | ||
| 34 | * This contains the quantities which are are constant within a finite | ||
| 35 | * volume in the Richards model. | ||
| 36 | */ | ||
| 37 | template <class Traits> | ||
| 38 | 257826 | class ExtendedRichardsVolumeVariables | |
| 39 | : public PorousMediumFlowVolumeVariables<Traits> | ||
| 40 | , public EnergyVolumeVariables<Traits, ExtendedRichardsVolumeVariables<Traits> > | ||
| 41 | { | ||
| 42 | using ParentType = PorousMediumFlowVolumeVariables<Traits>; | ||
| 43 | using EnergyVolVars = EnergyVolumeVariables<Traits, ExtendedRichardsVolumeVariables<Traits> >; | ||
| 44 | using Scalar = typename Traits::PrimaryVariables::value_type; | ||
| 45 | using PermeabilityType = typename Traits::PermeabilityType; | ||
| 46 | using ModelTraits = typename Traits::ModelTraits; | ||
| 47 | |||
| 48 | static constexpr int numFluidComps = ParentType::numFluidComponents(); | ||
| 49 | static constexpr int numPhases = ParentType::numFluidPhases(); | ||
| 50 | |||
| 51 | using EffDiffModel = typename Traits::EffectiveDiffusivityModel; | ||
| 52 | |||
| 53 | // checks if the fluid system uses the Richards model index convention | ||
| 54 | static constexpr auto fsCheck = ModelTraits::checkFluidSystem(typename Traits::FluidSystem{}); | ||
| 55 | |||
| 56 | public: | ||
| 57 | //! Export type of the fluid system | ||
| 58 | using FluidSystem = typename Traits::FluidSystem; | ||
| 59 | //! Export type of the fluid state | ||
| 60 | using FluidState = typename Traits::FluidState; | ||
| 61 | //! Export type of the fluid state | ||
| 62 | //! Export type of solid state | ||
| 63 | using SolidState = typename Traits::SolidState; | ||
| 64 | //! Export type of solid system | ||
| 65 | using SolidSystem = typename Traits::SolidSystem; | ||
| 66 | using Indices = typename Traits::ModelTraits::Indices; | ||
| 67 | using PrimaryVariableSwitch = ExtendedRichardsPrimaryVariableSwitch; | ||
| 68 | |||
| 69 | //! Export phase indices | ||
| 70 | static constexpr auto liquidPhaseIdx = Traits::FluidSystem::phase0Idx; | ||
| 71 | static constexpr auto gasPhaseIdx = Traits::FluidSystem::phase1Idx; | ||
| 72 | |||
| 73 | /*! | ||
| 74 | * \brief Updates all quantities for a given control volume. | ||
| 75 | * | ||
| 76 | * \param elemSol A vector containing all primary variables connected to the element | ||
| 77 | * \param problem The object specifying the problem which ought to | ||
| 78 | * be simulated | ||
| 79 | * \param element An element which contains part of the control volume | ||
| 80 | * \param scv The sub-control volume | ||
| 81 | */ | ||
| 82 | template<class ElemSol, class Problem, class Element, class Scv> | ||
| 83 | 538662 | void update(const ElemSol &elemSol, | |
| 84 | const Problem &problem, | ||
| 85 | const Element &element, | ||
| 86 | const Scv& scv) | ||
| 87 | { | ||
| 88 | 538662 | ParentType::update(elemSol, problem, element, scv); | |
| 89 | |||
| 90 | 1077324 | const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol); | |
| 91 | |||
| 92 | 538662 | const auto& priVars = elemSol[scv.localDofIndex()]; | |
| 93 | 538662 | const auto phasePresence = priVars.state(); | |
| 94 | |||
| 95 | // precompute the minimum capillary pressure (entry pressure) | ||
| 96 | // needed to make sure we don't compute unphysical capillary pressures and thus saturations | ||
| 97 |
2/2✓ Branch 0 taken 3571 times.
✓ Branch 1 taken 535091 times.
|
538662 | minPc_ = fluidMatrixInteraction.endPointPc(); |
| 98 | |||
| 99 | //update porosity before calculating the effective properties depending on it | ||
| 100 |
2/2✓ Branch 0 taken 3571 times.
✓ Branch 1 taken 535091 times.
|
538662 | updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps); |
| 101 | |||
| 102 | typename FluidSystem::ParameterCache paramCache; | ||
| 103 | 538662 | auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx) | |
| 104 | { | ||
| 105 | ✗ | return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx); | |
| 106 | }; | ||
| 107 | |||
| 108 |
2/2✓ Branch 0 taken 3571 times.
✓ Branch 1 taken 535091 times.
|
538662 | if (phasePresence == Indices::gasPhaseOnly) |
| 109 | { | ||
| 110 | 3571 | moleFraction_[liquidPhaseIdx] = 1.0; | |
| 111 | 3571 | massFraction_[liquidPhaseIdx] = 1.0; | |
| 112 | 3571 | moleFraction_[gasPhaseIdx] = priVars[Indices::switchIdx]; | |
| 113 | |||
| 114 | 7142 | const auto averageMolarMassGasPhase = (moleFraction_[gasPhaseIdx]*FluidSystem::molarMass(liquidPhaseIdx)) + | |
| 115 | 3571 | ((1-moleFraction_[gasPhaseIdx])*FluidSystem::molarMass(gasPhaseIdx)); | |
| 116 | |||
| 117 | //X_w = x_w* M_w/ M_avg | ||
| 118 | 3571 | massFraction_[gasPhaseIdx] = priVars[Indices::switchIdx]*FluidSystem::molarMass(liquidPhaseIdx)/averageMolarMassGasPhase; | |
| 119 | |||
| 120 | 3571 | fluidState_.setSaturation(liquidPhaseIdx, 0.0); | |
| 121 | 3571 | fluidState_.setSaturation(gasPhaseIdx, 1.0); | |
| 122 | |||
| 123 | 3571 | EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState_, solidState_); | |
| 124 | |||
| 125 | // get pc for sw = 0.0 | ||
| 126 | 3571 | const Scalar pc = fluidMatrixInteraction.pc(0.0); | |
| 127 | |||
| 128 | // set the wetting pressure | ||
| 129 | 3571 | fluidState_.setPressure(liquidPhaseIdx, problem.nonwettingReferencePressure() - pc); | |
| 130 | 3571 | fluidState_.setPressure(gasPhaseIdx, problem.nonwettingReferencePressure()); | |
| 131 | |||
| 132 | 14284 | molarDensity_[liquidPhaseIdx] = FluidSystem::H2O::liquidDensity(temperature(), pressure(liquidPhaseIdx))/FluidSystem::H2O::molarMass(); | |
| 133 | 10713 | molarDensity_[gasPhaseIdx] = IdealGas<Scalar>::molarDensity(temperature(), problem.nonwettingReferencePressure()); | |
| 134 | |||
| 135 | // density and viscosity | ||
| 136 | 3571 | paramCache.updateAll(fluidState_); | |
| 137 | 3571 | fluidState_.setDensity(liquidPhaseIdx, FluidSystem::density(fluidState_, paramCache, liquidPhaseIdx)); | |
| 138 | 3571 | fluidState_.setDensity(gasPhaseIdx, FluidSystem::density(fluidState_, paramCache, gasPhaseIdx)); | |
| 139 | 3571 | fluidState_.setViscosity(liquidPhaseIdx, FluidSystem::viscosity(fluidState_, paramCache, liquidPhaseIdx)); | |
| 140 | |||
| 141 | // compute and set the enthalpy | ||
| 142 | 3571 | fluidState_.setEnthalpy(liquidPhaseIdx, EnergyVolVars::enthalpy(fluidState_, paramCache, liquidPhaseIdx)); | |
| 143 | 3571 | fluidState_.setEnthalpy(gasPhaseIdx, EnergyVolVars::enthalpy(fluidState_, paramCache, gasPhaseIdx)); | |
| 144 | |||
| 145 | //binary diffusion coefficients | ||
| 146 | 3571 | paramCache.updateAll(fluidState_); | |
| 147 | 3571 | effectiveDiffCoeff_ = getEffectiveDiffusionCoefficient(gasPhaseIdx, | |
| 148 | FluidSystem::comp1Idx, | ||
| 149 | FluidSystem::comp0Idx); | ||
| 150 | } | ||
| 151 |
1/2✓ Branch 0 taken 535091 times.
✗ Branch 1 not taken.
|
535091 | else if (phasePresence == Indices::bothPhases) |
| 152 | { | ||
| 153 | 535091 | completeFluidState_(elemSol, problem, element, scv, fluidState_, solidState_); | |
| 154 | |||
| 155 | // we want to account for diffusion in the air phase | ||
| 156 | // use Raoult to compute the water mole fraction in air | ||
| 157 | 2140364 | molarDensity_[liquidPhaseIdx] = FluidSystem::H2O::liquidDensity(temperature(), pressure(liquidPhaseIdx))/FluidSystem::H2O::molarMass(); | |
| 158 | 1605273 | molarDensity_[gasPhaseIdx] = IdealGas<Scalar>::molarDensity(temperature(), problem.nonwettingReferencePressure()); | |
| 159 | 535091 | moleFraction_[liquidPhaseIdx] = 1.0; | |
| 160 | |||
| 161 | 1605273 | moleFraction_[gasPhaseIdx] = FluidSystem::H2O::vaporPressure(temperature()) / problem.nonwettingReferencePressure(); | |
| 162 | |||
| 163 | 1070182 | const auto averageMolarMassGasPhase = (moleFraction_[gasPhaseIdx]*FluidSystem::molarMass(liquidPhaseIdx)) + | |
| 164 | 535091 | ((1-moleFraction_[gasPhaseIdx])*FluidSystem::molarMass(gasPhaseIdx)); | |
| 165 | |||
| 166 | //X_w = x_w* M_w/ M_avg | ||
| 167 | 535091 | massFraction_[gasPhaseIdx] = moleFraction_[gasPhaseIdx]*FluidSystem::molarMass(liquidPhaseIdx)/averageMolarMassGasPhase; | |
| 168 | |||
| 169 | // binary diffusion coefficients | ||
| 170 | 535091 | paramCache.updateAll(fluidState_); | |
| 171 | 535091 | effectiveDiffCoeff_ = getEffectiveDiffusionCoefficient(gasPhaseIdx, | |
| 172 | FluidSystem::comp1Idx, | ||
| 173 | FluidSystem::comp0Idx); | ||
| 174 | } | ||
| 175 | ✗ | else if (phasePresence == Indices::liquidPhaseOnly) | |
| 176 | { | ||
| 177 | ✗ | completeFluidState_(elemSol, problem, element, scv, fluidState_, solidState_); | |
| 178 | |||
| 179 | ✗ | molarDensity_[liquidPhaseIdx] = FluidSystem::H2O::liquidDensity(temperature(), pressure(liquidPhaseIdx))/FluidSystem::H2O::molarMass(); | |
| 180 | ✗ | molarDensity_[gasPhaseIdx] = IdealGas<Scalar>::molarDensity(temperature(), problem.nonwettingReferencePressure()); | |
| 181 | ✗ | moleFraction_[liquidPhaseIdx] = 1.0; | |
| 182 | ✗ | moleFraction_[gasPhaseIdx] = 0.0; | |
| 183 | ✗ | massFraction_[liquidPhaseIdx] = 1.0; | |
| 184 | ✗ | massFraction_[gasPhaseIdx] = 0.0; | |
| 185 | |||
| 186 | // binary diffusion coefficients (none required for liquid phase only) | ||
| 187 | ✗ | effectiveDiffCoeff_ = 0.0; | |
| 188 | } | ||
| 189 | |||
| 190 | ////////// | ||
| 191 | // specify the other parameters | ||
| 192 | ////////// | ||
| 193 | 1077324 | relativePermeabilityWetting_ = fluidMatrixInteraction.krw(fluidState_.saturation(liquidPhaseIdx)); | |
| 194 | 538662 | EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_); | |
| 195 | 1077324 | permeability_ = problem.spatialParams().permeability(element, scv, elemSol); | |
| 196 | 538662 | EnergyVolVars::updateEffectiveThermalConductivity(); | |
| 197 | 538662 | } | |
| 198 | |||
| 199 | /*! | ||
| 200 | * \brief Returns the fluid configuration at the given primary | ||
| 201 | * variables. | ||
| 202 | */ | ||
| 203 | const FluidState &fluidState() const | ||
| 204 |
2/4✓ Branch 0 taken 478494 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 447336 times.
✗ Branch 3 not taken.
|
3627192 | { return fluidState_; } |
| 205 | |||
| 206 | /*! | ||
| 207 | * \brief Returns the phase state for the control volume. | ||
| 208 | */ | ||
| 209 | const SolidState &solidState() const | ||
| 210 | 538662 | { return solidState_; } | |
| 211 | |||
| 212 | /*! | ||
| 213 | * \brief Returns the temperature. | ||
| 214 | */ | ||
| 215 | Scalar temperature() const | ||
| 216 |
1/2✓ Branch 0 taken 364230 times.
✗ Branch 1 not taken.
|
4060031 | { return fluidState_.temperature(); } |
| 217 | |||
| 218 | /*! | ||
| 219 | * \brief Returns the average porosity [] within the control volume. | ||
| 220 | * | ||
| 221 | * The porosity is defined as the ratio of the pore space to the | ||
| 222 | * total volume, i.e. \f[ \Phi := \frac{V_{pore}}{V_{pore} + V_{rock}} \f] | ||
| 223 | */ | ||
| 224 | Scalar porosity() const | ||
| 225 |
2/4✗ Branch 4 not taken.
✓ Branch 5 taken 538662 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 538662 times.
|
9509448 | { return solidState_.porosity(); } |
| 226 | |||
| 227 | /*! | ||
| 228 | * \brief Returns the permeability within the control volume in \f$[m^2]\f$. | ||
| 229 | */ | ||
| 230 | const PermeabilityType& permeability() const | ||
| 231 | ✗ | { return permeability_; } | |
| 232 | |||
| 233 | /*! | ||
| 234 | * \brief Returns the average absolute saturation [] of a given | ||
| 235 | * fluid phase within the finite volume. | ||
| 236 | * | ||
| 237 | * The saturation of a fluid phase is defined as the fraction of | ||
| 238 | * the pore volume filled by it, i.e. | ||
| 239 | * \f[ S_\alpha := \frac{V_\alpha}{V_{pore}} = \phi \frac{V_\alpha}{V} \f] | ||
| 240 | * | ||
| 241 | * \param phaseIdx The index of the fluid phase | ||
| 242 | */ | ||
| 243 | Scalar saturation(const int phaseIdx = liquidPhaseIdx) const | ||
| 244 |
6/8✗ Branch 2 not taken.
✓ Branch 3 taken 538662 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 538662 times.
✓ Branch 6 taken 11 times.
✓ Branch 7 taken 35461 times.
✓ Branch 8 taken 11 times.
✓ Branch 9 taken 35461 times.
|
5902992 | { return fluidState_.saturation(phaseIdx); } |
| 245 | |||
| 246 | /*! | ||
| 247 | * \brief Returns the average mass density \f$\mathrm{[kg/m^3]}\f$ of a given | ||
| 248 | * fluid phase within the control volume. | ||
| 249 | * | ||
| 250 | * \param phaseIdx The index of the fluid phase | ||
| 251 | */ | ||
| 252 | Scalar density(const int phaseIdx = liquidPhaseIdx) const | ||
| 253 |
8/20✓ Branch 0 taken 478494 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 478494 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 478494 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 478494 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 447336 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 447336 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 447336 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 447336 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
17055900 | { return fluidState_.density(phaseIdx); } |
| 254 | |||
| 255 | /*! | ||
| 256 | * \brief Returns the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within | ||
| 257 | * the control volume. | ||
| 258 | * | ||
| 259 | * For the nonwetting phase (i.e. the gas phase), we assume | ||
| 260 | * infinite mobility, which implies that the nonwetting phase | ||
| 261 | * pressure is equal to the finite volume's reference pressure | ||
| 262 | * defined by the problem. | ||
| 263 | * | ||
| 264 | * \param phaseIdx The index of the fluid phase | ||
| 265 | */ | ||
| 266 | Scalar pressure(const int phaseIdx = liquidPhaseIdx) const | ||
| 267 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 348 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 348 times.
|
7024169 | { return fluidState_.pressure(phaseIdx); } |
| 268 | |||
| 269 | /*! | ||
| 270 | * \brief Returns the effective mobility \f$\mathrm{[1/(Pa*s)]}\f$ of a given phase within | ||
| 271 | * the control volume. | ||
| 272 | * | ||
| 273 | * The mobility of a fluid phase is defined as the relative | ||
| 274 | * permeability of the phase (given by the chosen material law) | ||
| 275 | * divided by the dynamic viscosity of the fluid, i.e. | ||
| 276 | * \f[ \lambda_\alpha := \frac{k_{r\alpha}}{\mu_\alpha} \f] | ||
| 277 | * | ||
| 278 | * \param phaseIdx The index of the fluid phase | ||
| 279 | */ | ||
| 280 | Scalar mobility(const int phaseIdx = liquidPhaseIdx) const | ||
| 281 |
4/16✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 478494 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 478494 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 447336 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 447336 times.
✗ Branch 15 not taken.
|
5554980 | { return relativePermeability(phaseIdx)/fluidState_.viscosity(phaseIdx); } |
| 282 | |||
| 283 | /*! | ||
| 284 | * \brief Returns the dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of a given phase within | ||
| 285 | * the control volume. | ||
| 286 | * | ||
| 287 | * \param phaseIdx The index of the fluid phase | ||
| 288 | * \note The nonwetting phase is infinitely mobile | ||
| 289 | */ | ||
| 290 | Scalar viscosity(const int phaseIdx = liquidPhaseIdx) const | ||
| 291 | { return phaseIdx == liquidPhaseIdx ? fluidState_.viscosity(liquidPhaseIdx) : 0.0; } | ||
| 292 | |||
| 293 | /*! | ||
| 294 | * \brief Returns relative permeability [-] of a given phase within | ||
| 295 | * the control volume. | ||
| 296 | * | ||
| 297 | * \param phaseIdx The index of the fluid phase | ||
| 298 | */ | ||
| 299 | ✗ | Scalar relativePermeability(const int phaseIdx = liquidPhaseIdx) const | |
| 300 |
4/22✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 478494 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 478494 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 447336 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 447336 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
2777490 | { return phaseIdx == liquidPhaseIdx ? relativePermeabilityWetting_ : 1.0; } |
| 301 | |||
| 302 | /*! | ||
| 303 | * \brief Returns the effective capillary pressure \f$\mathrm{[Pa]}\f$ within the | ||
| 304 | * control volume. | ||
| 305 | * | ||
| 306 | * The capillary pressure is defined as the difference in | ||
| 307 | * pressures of the nonwetting and the wetting phase, i.e. | ||
| 308 | * \f[ p_c = p_n - p_w \f] | ||
| 309 | * | ||
| 310 | * \note Capillary pressures are always larger than the entry pressure | ||
| 311 | * This regularization doesn't affect the residual in which pc is not needed. | ||
| 312 | */ | ||
| 313 | Scalar capillaryPressure() const | ||
| 314 | { | ||
| 315 | using std::max; | ||
| 316 |
3/12✓ Branch 0 taken 33000 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 33000 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 33000 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
99000 | return max(minPc_, pressure(gasPhaseIdx) - pressure(liquidPhaseIdx)); |
| 317 | } | ||
| 318 | |||
| 319 | /*! | ||
| 320 | * \brief Returns the pressureHead \f$\mathrm{[cm]}\f$ of a given phase within | ||
| 321 | * the control volume. | ||
| 322 | * | ||
| 323 | * For the nonwetting phase (i.e. the gas phase), we assume | ||
| 324 | * infinite mobility, which implies that the nonwetting phase | ||
| 325 | * pressure is equal to the finite volume's reference pressure | ||
| 326 | * defined by the problem. | ||
| 327 | * | ||
| 328 | * \param phaseIdx The index of the fluid phase | ||
| 329 | * \note this function is here as a convenience to the user to not have to | ||
| 330 | * manually do a conversion. It is not correct if the density is not constant | ||
| 331 | * or the gravity different | ||
| 332 | */ | ||
| 333 | Scalar pressureHead(const int phaseIdx = liquidPhaseIdx) const | ||
| 334 | ✗ | { return 100.0 *(pressure(phaseIdx) - pressure(gasPhaseIdx))/density(phaseIdx)/9.81; } | |
| 335 | |||
| 336 | /*! | ||
| 337 | * \brief Returns the water content of a fluid phase within the finite volume. | ||
| 338 | * | ||
| 339 | * The water content is defined as the fraction of | ||
| 340 | * the saturation divided by the porosity. | ||
| 341 | |||
| 342 | * \param phaseIdx The index of the fluid phase | ||
| 343 | * \note this function is here as a convenience to the user to not have to | ||
| 344 | * manually do a conversion. | ||
| 345 | */ | ||
| 346 | Scalar waterContent(const int phaseIdx = liquidPhaseIdx) const | ||
| 347 | 99000 | { return saturation(phaseIdx) * solidState_.porosity(); } | |
| 348 | |||
| 349 | /*! | ||
| 350 | * \brief Returns the mole fraction of a given component in a | ||
| 351 | * given phase within the control volume in \f$[-]\f$. | ||
| 352 | * | ||
| 353 | * \param phaseIdx The phase index | ||
| 354 | * \param compIdx The component index | ||
| 355 | */ | ||
| 356 | 1258800 | Scalar moleFraction(const int phaseIdx, const int compIdx) const | |
| 357 | { | ||
| 358 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1258800 times.
|
1258800 | if (compIdx != FluidSystem::comp0Idx) |
| 359 | ✗ | DUNE_THROW(Dune::InvalidStateException, "There is only one component for Richards!"); | |
| 360 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 11 times.
|
1259159 | return moleFraction_[phaseIdx]; |
| 361 | } | ||
| 362 | |||
| 363 | /*! | ||
| 364 | * \brief Returns the mole fraction of a given component in a | ||
| 365 | * given phase within the control volume in \f$[-]\f$. | ||
| 366 | * | ||
| 367 | * \param phaseIdx The phase index | ||
| 368 | * \param compIdx The component index | ||
| 369 | */ | ||
| 370 | 2985930 | Scalar massFraction(const int phaseIdx, const int compIdx) const | |
| 371 | { | ||
| 372 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2985930 times.
|
2985930 | if (compIdx != FluidSystem::comp0Idx) |
| 373 | ✗ | DUNE_THROW(Dune::InvalidStateException, "There is only one component for Richards!"); | |
| 374 | 2985930 | return massFraction_[phaseIdx]; | |
| 375 | } | ||
| 376 | |||
| 377 | /*! | ||
| 378 | * \brief Returns the mass density of a given phase within the | ||
| 379 | * control volume in \f$[mol/m^3]\f$. | ||
| 380 | * | ||
| 381 | * \param phaseIdx The phase index | ||
| 382 | */ | ||
| 383 | Scalar molarDensity(const int phaseIdx) const | ||
| 384 | { | ||
| 385 | 1225800 | return molarDensity_[phaseIdx]; | |
| 386 | } | ||
| 387 | |||
| 388 | /*! | ||
| 389 | * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$. | ||
| 390 | */ | ||
| 391 | 538662 | Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const | |
| 392 | { | ||
| 393 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 538662 times.
|
538662 | assert(phaseIdx == gasPhaseIdx); |
| 394 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 538662 times.
|
538662 | assert(compIIdx != compJIdx); |
| 395 | typename FluidSystem::ParameterCache paramCache; | ||
| 396 | 538662 | paramCache.updatePhase(fluidState_, phaseIdx); | |
| 397 | 538662 | return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx); | |
| 398 | } | ||
| 399 | |||
| 400 | /*! | ||
| 401 | * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$. | ||
| 402 | */ | ||
| 403 | ✗ | Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const | |
| 404 | { | ||
| 405 | ✗ | assert(phaseIdx == gasPhaseIdx); | |
| 406 | ✗ | assert(compIIdx != compJIdx); | |
| 407 | 11070 | return effectiveDiffCoeff_; | |
| 408 | } | ||
| 409 | private: | ||
| 410 | /*! | ||
| 411 | * \brief Fills the fluid state according to the primary variables. | ||
| 412 | * | ||
| 413 | * Taking the information from the primary variables, | ||
| 414 | * the fluid state is filled with every information that is | ||
| 415 | * necessary to evaluate the model's local residual. | ||
| 416 | * | ||
| 417 | * \param elemSol A vector containing all primary variables connected to the element | ||
| 418 | * \param problem The problem at hand. | ||
| 419 | * \param element The current element. | ||
| 420 | * \param scv The subcontrol volume. | ||
| 421 | * \param fluidState The fluid state to fill. | ||
| 422 | * \param solidState The solid state to fill. | ||
| 423 | */ | ||
| 424 | template<class ElemSol, class Problem, class Element, class Scv> | ||
| 425 | 535091 | void completeFluidState_(const ElemSol& elemSol, | |
| 426 | const Problem& problem, | ||
| 427 | const Element& element, | ||
| 428 | const Scv& scv, | ||
| 429 | FluidState& fluidState, | ||
| 430 | SolidState& solidState) | ||
| 431 | { | ||
| 432 | 535091 | EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState, solidState); | |
| 433 | |||
| 434 | 1070182 | const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol); | |
| 435 | |||
| 436 | 535091 | const auto& priVars = elemSol[scv.localDofIndex()]; | |
| 437 | |||
| 438 | // set the wetting pressure | ||
| 439 | using std::max; | ||
| 440 | 535091 | Scalar minPc = fluidMatrixInteraction.pc(1.0); | |
| 441 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 535091 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 535091 times.
|
1070182 | fluidState.setPressure(liquidPhaseIdx, priVars[Indices::pressureIdx]); |
| 442 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 535091 times.
✓ Branch 2 taken 535091 times.
✗ Branch 3 not taken.
|
535091 | fluidState.setPressure(gasPhaseIdx, max(problem.nonwettingReferencePressure(), fluidState.pressure(liquidPhaseIdx) + minPc)); |
| 443 | |||
| 444 | // compute the capillary pressure to compute the saturation | ||
| 445 | // make sure that we the capillary pressure is not smaller than the minimum pc | ||
| 446 | // this would possibly return unphysical values from regularized material laws | ||
| 447 | using std::max; | ||
| 448 |
1/2✓ Branch 0 taken 535091 times.
✗ Branch 1 not taken.
|
535091 | const Scalar pc = max(fluidMatrixInteraction.endPointPc(), |
| 449 |
1/2✓ Branch 0 taken 535091 times.
✗ Branch 1 not taken.
|
535091 | problem.nonwettingReferencePressure() - fluidState.pressure(liquidPhaseIdx)); |
| 450 | 535091 | const Scalar sw = fluidMatrixInteraction.sw(pc); | |
| 451 | 535091 | fluidState.setSaturation(liquidPhaseIdx, sw); | |
| 452 | 535091 | fluidState.setSaturation(gasPhaseIdx, 1.0-sw); | |
| 453 | |||
| 454 | // density and viscosity | ||
| 455 | typename FluidSystem::ParameterCache paramCache; | ||
| 456 | 535091 | paramCache.updateAll(fluidState); | |
| 457 | 535091 | fluidState.setDensity(liquidPhaseIdx, | |
| 458 | FluidSystem::density(fluidState, paramCache, liquidPhaseIdx)); | ||
| 459 | 535091 | fluidState.setDensity(gasPhaseIdx, | |
| 460 | FluidSystem::density(fluidState, paramCache, gasPhaseIdx)); | ||
| 461 | |||
| 462 | 535091 | fluidState.setViscosity(liquidPhaseIdx, | |
| 463 | FluidSystem::viscosity(fluidState, paramCache, liquidPhaseIdx)); | ||
| 464 | |||
| 465 | // compute and set the enthalpy | ||
| 466 | 535091 | fluidState.setEnthalpy(liquidPhaseIdx, EnergyVolVars::enthalpy(fluidState, paramCache, liquidPhaseIdx)); | |
| 467 | 535091 | fluidState.setEnthalpy(gasPhaseIdx, EnergyVolVars::enthalpy(fluidState, paramCache, gasPhaseIdx)); | |
| 468 | 535091 | } | |
| 469 | |||
| 470 | protected: | ||
| 471 | FluidState fluidState_; //!< the fluid state | ||
| 472 | SolidState solidState_; | ||
| 473 | Scalar relativePermeabilityWetting_; //!< the relative permeability of the wetting phase | ||
| 474 | PermeabilityType permeability_; //!< the intrinsic permeability | ||
| 475 | Scalar minPc_; //!< the minimum capillary pressure (entry pressure) | ||
| 476 | Scalar moleFraction_[numPhases]; //!< The water mole fractions in water and air | ||
| 477 | Scalar massFraction_[numPhases]; //!< The water mass fractions in water and air | ||
| 478 | Scalar molarDensity_[numPhases]; //!< The molar density of water and air | ||
| 479 | |||
| 480 | // Effective diffusion coefficients for the phases | ||
| 481 | Scalar effectiveDiffCoeff_; | ||
| 482 | |||
| 483 | }; | ||
| 484 | } // end namespace Dumux | ||
| 485 | |||
| 486 | #endif | ||
| 487 |