| 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 TwoPTwoCModel | ||
| 10 | * \brief Contains the quantities which are constant within a | ||
| 11 | * finite volume in the two-phase two-component model. | ||
| 12 | */ | ||
| 13 | |||
| 14 | #ifndef DUMUX_2P2C_VOLUME_VARIABLES_HH | ||
| 15 | #define DUMUX_2P2C_VOLUME_VARIABLES_HH | ||
| 16 | |||
| 17 | #include <dumux/material/fluidstates/compositional.hh> | ||
| 18 | #include <dumux/material/constraintsolvers/computefromreferencephase.hh> | ||
| 19 | #include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh> | ||
| 20 | |||
| 21 | #include <dumux/discretization/method.hh> | ||
| 22 | #include <dumux/porousmediumflow/volumevariables.hh> | ||
| 23 | #include <dumux/porousmediumflow/nonisothermal/volumevariables.hh> | ||
| 24 | #include <dumux/porousmediumflow/2p/formulation.hh> | ||
| 25 | #include <dumux/material/solidstates/updatesolidvolumefractions.hh> | ||
| 26 | #include <dumux/porousmediumflow/2pnc/primaryvariableswitch.hh> | ||
| 27 | |||
| 28 | namespace Dumux { | ||
| 29 | |||
| 30 | // forward declaration | ||
| 31 | template <class Traits, bool enableChemicalNonEquilibrium, bool useConstraintSolver> | ||
| 32 | class TwoPTwoCVolumeVariablesImplementation; | ||
| 33 | |||
| 34 | |||
| 35 | /*! | ||
| 36 | * \ingroup TwoPTwoCModel | ||
| 37 | * \brief Contains the quantities which are constant within a | ||
| 38 | * finite volume in the two-phase two-component model. | ||
| 39 | */ | ||
| 40 | template <class Traits, bool useConstraintSolver = true> | ||
| 41 | using TwoPTwoCVolumeVariables = TwoPTwoCVolumeVariablesImplementation<Traits, Traits::ModelTraits::enableChemicalNonEquilibrium(), useConstraintSolver>; | ||
| 42 | |||
| 43 | /*! | ||
| 44 | * \ingroup TwoPTwoCModel | ||
| 45 | * \brief Contains the quantities which are constant within a | ||
| 46 | * finite volume in the two-phase two-component model. | ||
| 47 | * This is the base class for a 2p2c model with and without chemical nonequilibrium | ||
| 48 | */ | ||
| 49 | template <class Traits, class Impl> | ||
| 50 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 8132 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13056 times.
✓ Branch 4 taken 2432 times.
✗ Branch 5 not taken.
|
1443508 | class TwoPTwoCVolumeVariablesBase |
| 51 | : public PorousMediumFlowVolumeVariables<Traits> | ||
| 52 | , public EnergyVolumeVariables<Traits, TwoPTwoCVolumeVariables<Traits> > | ||
| 53 | { | ||
| 54 | using ParentType = PorousMediumFlowVolumeVariables< Traits>; | ||
| 55 | using EnergyVolVars = EnergyVolumeVariables<Traits, TwoPTwoCVolumeVariables<Traits> >; | ||
| 56 | |||
| 57 | using Scalar = typename Traits::PrimaryVariables::value_type; | ||
| 58 | using ModelTraits = typename Traits::ModelTraits; | ||
| 59 | |||
| 60 | static constexpr int numFluidComps = ParentType::numFluidComponents(); | ||
| 61 | // component indices | ||
| 62 | enum | ||
| 63 | { | ||
| 64 | comp0Idx = Traits::FluidSystem::comp0Idx, | ||
| 65 | comp1Idx = Traits::FluidSystem::comp1Idx, | ||
| 66 | phase0Idx = Traits::FluidSystem::phase0Idx, | ||
| 67 | phase1Idx = Traits::FluidSystem::phase1Idx | ||
| 68 | }; | ||
| 69 | |||
| 70 | // phase presence indices | ||
| 71 | enum | ||
| 72 | { | ||
| 73 | firstPhaseOnly = ModelTraits::Indices::firstPhaseOnly, | ||
| 74 | secondPhaseOnly = ModelTraits::Indices::secondPhaseOnly, | ||
| 75 | bothPhases = ModelTraits::Indices::bothPhases | ||
| 76 | }; | ||
| 77 | |||
| 78 | // primary variable indices | ||
| 79 | enum | ||
| 80 | { | ||
| 81 | switchIdx = ModelTraits::Indices::switchIdx, | ||
| 82 | pressureIdx = ModelTraits::Indices::pressureIdx | ||
| 83 | }; | ||
| 84 | |||
| 85 | // formulations | ||
| 86 | static constexpr auto formulation = ModelTraits::priVarFormulation(); | ||
| 87 | |||
| 88 | using PermeabilityType = typename Traits::PermeabilityType; | ||
| 89 | using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, typename Traits::FluidSystem>; | ||
| 90 | using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition< Scalar, typename Traits::FluidSystem >; | ||
| 91 | using EffDiffModel = typename Traits::EffectiveDiffusivityModel; | ||
| 92 | using DiffusionCoefficients = typename Traits::DiffusionType::DiffusionCoefficientsContainer; | ||
| 93 | public: | ||
| 94 | //! The type of the object returned by the fluidState() method | ||
| 95 | using FluidState = typename Traits::FluidState; | ||
| 96 | //! The fluid system used here | ||
| 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 = TwoPNCPrimaryVariableSwitch; | ||
| 106 | |||
| 107 | //! Return whether moles or masses are balanced | ||
| 108 | static constexpr bool useMoles() { return ModelTraits::useMoles(); } | ||
| 109 | //! Return the two-phase formulation used here | ||
| 110 | static constexpr TwoPFormulation priVarFormulation() { return formulation; } | ||
| 111 | |||
| 112 | // check for permissive combinations | ||
| 113 | static_assert(ModelTraits::numFluidPhases() == 2, "NumPhases set in the model is not two!"); | ||
| 114 | static_assert(ModelTraits::numFluidComponents() == 2, "NumComponents set in the model is not two!"); | ||
| 115 | static_assert((formulation == TwoPFormulation::p0s1 || formulation == TwoPFormulation::p1s0), "Chosen TwoPFormulation not supported!"); | ||
| 116 | |||
| 117 | /*! | ||
| 118 | * \brief Updates all quantities for a given control volume. | ||
| 119 | * | ||
| 120 | * \param elemSol A vector containing all primary variables connected to the element | ||
| 121 | * \param problem The object specifying the problem which ought to | ||
| 122 | * be simulated | ||
| 123 | * \param element An element which contains part of the control volume | ||
| 124 | * \param scv The sub control volume | ||
| 125 | */ | ||
| 126 | template<class ElemSol, class Problem, class Element, class Scv> | ||
| 127 | 13851694 | void update(const ElemSol& elemSol, const Problem& problem, const Element& element, const Scv& scv) | |
| 128 | { | ||
| 129 | 13851694 | ParentType::update(elemSol, problem, element, scv); | |
| 130 | 13851694 | asImp_().completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_); | |
| 131 | |||
| 132 |
2/2✓ Branch 0 taken 581621 times.
✓ Branch 1 taken 670699 times.
|
13851694 | const auto& spatialParams = problem.spatialParams(); |
| 133 |
2/2✓ Branch 0 taken 581621 times.
✓ Branch 1 taken 670699 times.
|
13851694 | const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol); |
| 134 | |||
| 135 | 13851694 | const int wPhaseIdx = fluidState_.wettingPhase(); | |
| 136 | 13851694 | const int nPhaseIdx = 1 - wPhaseIdx; | |
| 137 | |||
| 138 | // relative permeabilities -> require wetting phase saturation as parameter! | ||
| 139 | 13851694 | relativePermeability_[wPhaseIdx] = fluidMatrixInteraction.krw(saturation(wPhaseIdx)); | |
| 140 |
2/2✓ Branch 0 taken 5848167 times.
✓ Branch 1 taken 7016643 times.
|
13851694 | relativePermeability_[nPhaseIdx] = fluidMatrixInteraction.krn(saturation(wPhaseIdx)); |
| 141 | |||
| 142 | // porosity & permeabilty | ||
| 143 |
2/2✓ Branch 0 taken 5848167 times.
✓ Branch 1 taken 7016643 times.
|
13851694 | updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps); |
| 144 | 13851694 | EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_); | |
| 145 |
2/2✓ Branch 0 taken 5848167 times.
✓ Branch 1 taken 7016643 times.
|
13851694 | permeability_ = spatialParams.permeability(element, scv, elemSol); |
| 146 | |||
| 147 | 13851694 | auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx) | |
| 148 | { | ||
| 149 | 27703388 | return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx); | |
| 150 | }; | ||
| 151 | |||
| 152 | 13851694 | effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient); | |
| 153 | |||
| 154 | 13851694 | EnergyVolVars::updateEffectiveThermalConductivity(); | |
| 155 | 13851694 | } | |
| 156 | |||
| 157 | /*! | ||
| 158 | * \brief Sets complete fluid state. | ||
| 159 | * | ||
| 160 | * \param elemSol A vector containing all primary variables connected to the element | ||
| 161 | * \param problem The object specifying the problem which ought to | ||
| 162 | * be simulated | ||
| 163 | * \param element An element which contains part of the control volume | ||
| 164 | * \param scv The sub-control volume | ||
| 165 | * \param fluidState A container with the current (physical) state of the fluid | ||
| 166 | * \param solidState A container with the current (physical) state of the solid | ||
| 167 | * | ||
| 168 | * Set temperature, saturations, capillary pressures, viscosities, densities and enthalpies. | ||
| 169 | */ | ||
| 170 | template<class ElemSol, class Problem, class Element, class Scv> | ||
| 171 | 13851694 | void completeFluidState(const ElemSol& elemSol, | |
| 172 | const Problem& problem, | ||
| 173 | const Element& element, | ||
| 174 | const Scv& scv, | ||
| 175 | FluidState& fluidState, | ||
| 176 | SolidState& solidState) | ||
| 177 | { | ||
| 178 | 13851694 | EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState, solidState); | |
| 179 | |||
| 180 | 13851694 | const auto& priVars = elemSol[scv.localDofIndex()]; | |
| 181 | 13851694 | const auto phasePresence = priVars.state(); | |
| 182 | |||
| 183 |
2/2✓ Branch 0 taken 12077574 times.
✓ Branch 1 taken 1774120 times.
|
13851694 | const auto& spatialParams = problem.spatialParams(); |
| 184 |
2/2✓ Branch 0 taken 12077574 times.
✓ Branch 1 taken 1774120 times.
|
13851694 | const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol); |
| 185 |
2/2✓ Branch 0 taken 12183592 times.
✓ Branch 1 taken 1668102 times.
|
13851694 | const auto wPhaseIdx = spatialParams.template wettingPhase<FluidSystem>(element, scv, elemSol); |
| 186 |
2/2✓ Branch 0 taken 12183592 times.
✓ Branch 1 taken 1668102 times.
|
13851694 | fluidState.setWettingPhase(wPhaseIdx); |
| 187 | |||
| 188 | // set the saturations | ||
| 189 |
2/2✓ Branch 0 taken 12183592 times.
✓ Branch 1 taken 1668102 times.
|
13851694 | if (phasePresence == firstPhaseOnly) |
| 190 | { | ||
| 191 | 12183592 | fluidState.setSaturation(phase0Idx, 1.0); | |
| 192 | 12183592 | fluidState.setSaturation(phase1Idx, 0.0); | |
| 193 | } | ||
| 194 |
2/2✓ Branch 0 taken 17974 times.
✓ Branch 1 taken 1650128 times.
|
1668102 | else if (phasePresence == secondPhaseOnly) |
| 195 | { | ||
| 196 | 17974 | fluidState.setSaturation(phase0Idx, 0.0); | |
| 197 | 17974 | fluidState.setSaturation(phase1Idx, 1.0); | |
| 198 | } | ||
| 199 |
1/2✓ Branch 0 taken 1650128 times.
✗ Branch 1 not taken.
|
1650128 | else if (phasePresence == bothPhases) |
| 200 | { | ||
| 201 | if (formulation == TwoPFormulation::p0s1) | ||
| 202 | { | ||
| 203 | 1747548 | fluidState.setSaturation(phase1Idx, priVars[switchIdx]); | |
| 204 | 1747548 | fluidState.setSaturation(phase0Idx, 1 - priVars[switchIdx]); | |
| 205 | } | ||
| 206 | else | ||
| 207 | { | ||
| 208 | 1552708 | fluidState.setSaturation(phase0Idx, priVars[switchIdx]); | |
| 209 | 1552708 | fluidState.setSaturation(phase1Idx, 1 - priVars[switchIdx]); | |
| 210 | } | ||
| 211 | } | ||
| 212 | else | ||
| 213 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase presence."); | |
| 214 | |||
| 215 | // set pressures of the fluid phases | ||
| 216 | 27703388 | pc_ = fluidMatrixInteraction.pc(fluidState.saturation(wPhaseIdx)); | |
| 217 | if (formulation == TwoPFormulation::p0s1) | ||
| 218 | { | ||
| 219 | 26114732 | fluidState.setPressure(phase0Idx, priVars[pressureIdx]); | |
| 220 | 26114732 | fluidState.setPressure(phase1Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] + pc_ | |
| 221 | : priVars[pressureIdx] - pc_); | ||
| 222 | } | ||
| 223 | else | ||
| 224 | { | ||
| 225 | 1588656 | fluidState.setPressure(phase1Idx, priVars[pressureIdx]); | |
| 226 | 1588656 | fluidState.setPressure(phase0Idx, (wPhaseIdx == phase0Idx) ? priVars[pressureIdx] - pc_ | |
| 227 | : priVars[pressureIdx] + pc_); | ||
| 228 | } | ||
| 229 | 13851694 | } | |
| 230 | |||
| 231 | /*! | ||
| 232 | * \brief Returns the phase state within the control volume. | ||
| 233 | */ | ||
| 234 | const FluidState &fluidState() const | ||
| 235 |
3/4✓ Branch 0 taken 12248 times.
✓ Branch 1 taken 72 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 48356 times.
|
318018162 | { return fluidState_; } |
| 236 | |||
| 237 | /*! | ||
| 238 | * \brief Returns the phase state for the control-volume. | ||
| 239 | */ | ||
| 240 | const SolidState &solidState() const | ||
| 241 | 12491118 | { return solidState_; } | |
| 242 | |||
| 243 | /*! | ||
| 244 | * \brief Returns the average molar mass \f$\mathrm{[kg/mol]}\f$ of the fluid phase. | ||
| 245 | * | ||
| 246 | * \param phaseIdx The phase index | ||
| 247 | */ | ||
| 248 | Scalar averageMolarMass(int phaseIdx) const | ||
| 249 | { return fluidState_.averageMolarMass(phaseIdx); } | ||
| 250 | |||
| 251 | /*! | ||
| 252 | * \brief Returns the saturation of a given phase within | ||
| 253 | * the control volume in \f$[-]\f$. | ||
| 254 | * | ||
| 255 | * \param phaseIdx The phase index | ||
| 256 | */ | ||
| 257 | Scalar saturation(const int phaseIdx) const | ||
| 258 |
16/26✓ Branch 0 taken 24 times.
✓ Branch 1 taken 2130962 times.
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 2130962 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 36482 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 36482 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 74770 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 74770 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 36482 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 36482 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 32 times.
✓ Branch 21 taken 664928 times.
✓ Branch 22 taken 48 times.
✓ Branch 23 taken 25647148 times.
✓ Branch 24 taken 16 times.
✓ Branch 25 taken 24982220 times.
|
919116616 | { return fluidState_.saturation(phaseIdx); } |
| 259 | |||
| 260 | /*! | ||
| 261 | * \brief Returns the mass fraction of a given component in a | ||
| 262 | * given phase within the control volume in \f$[-]\f$. | ||
| 263 | * | ||
| 264 | * \param phaseIdx The phase index | ||
| 265 | * \param compIdx The component index | ||
| 266 | */ | ||
| 267 | Scalar massFraction(const int phaseIdx, const int compIdx) const | ||
| 268 |
2/3✓ Branch 10 taken 4992 times.
✓ Branch 11 taken 5120 times.
✗ Branch 12 not taken.
|
298014012 | { return fluidState_.massFraction(phaseIdx, compIdx); } |
| 269 | |||
| 270 | /*! | ||
| 271 | * \brief Returns the mole fraction of a given component in a | ||
| 272 | * given phase within the control volume in \f$[-]\f$. | ||
| 273 | * | ||
| 274 | * \param phaseIdx The phase index | ||
| 275 | * \param compIdx The component index | ||
| 276 | */ | ||
| 277 | Scalar moleFraction(const int phaseIdx, const int compIdx) const | ||
| 278 |
6/8✓ Branch 0 taken 97072 times.
✓ Branch 1 taken 62076 times.
✓ Branch 2 taken 97072 times.
✓ Branch 3 taken 62076 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 170512 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 170512 times.
|
1295155592 | { return fluidState_.moleFraction(phaseIdx, compIdx); } |
| 279 | |||
| 280 | /*! | ||
| 281 | * \brief Returns the mass density of a given phase within the | ||
| 282 | * control volume in \f$[kg/m^3]\f$. | ||
| 283 | * | ||
| 284 | * \param phaseIdx The phase index | ||
| 285 | */ | ||
| 286 | Scalar density(const int phaseIdx) const | ||
| 287 |
9/12✓ Branch 0 taken 1885144 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1897392 times.
✓ Branch 3 taken 72 times.
✓ Branch 4 taken 64536 times.
✓ Branch 5 taken 86072 times.
✓ Branch 6 taken 52288 times.
✓ Branch 7 taken 86000 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 48356 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 48356 times.
|
1864689184 | { return fluidState_.density(phaseIdx); } |
| 288 | |||
| 289 | /*! | ||
| 290 | * \brief Returns the dynamic viscosity of the fluid within the | ||
| 291 | * control volume in \f$\mathrm{[Pa s]}\f$. | ||
| 292 | * | ||
| 293 | * \param phaseIdx The phase index | ||
| 294 | */ | ||
| 295 | Scalar viscosity(const int phaseIdx) const | ||
| 296 | { return fluidState_.viscosity(phaseIdx); } | ||
| 297 | |||
| 298 | /*! | ||
| 299 | * \brief Returns the mass density of a given phase within the | ||
| 300 | * control volume in \f$[mol/m^3]\f$. | ||
| 301 | * | ||
| 302 | * \param phaseIdx The phase index | ||
| 303 | */ | ||
| 304 | Scalar molarDensity(const int phaseIdx) const | ||
| 305 |
6/8✓ Branch 0 taken 43584 times.
✓ Branch 1 taken 256 times.
✓ Branch 2 taken 43584 times.
✓ Branch 3 taken 256 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 170512 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 170512 times.
|
1291209848 | { return fluidState_.molarDensity(phaseIdx); } |
| 306 | |||
| 307 | /*! | ||
| 308 | * \brief Returns the effective pressure of a given phase within | ||
| 309 | * the control volume in \f$[kg/(m*s^2)=N/m^2=Pa]\f$. | ||
| 310 | * | ||
| 311 | * \param phaseIdx The phase index | ||
| 312 | */ | ||
| 313 | Scalar pressure(const int phaseIdx) const | ||
| 314 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 85256 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 85256 times.
|
635481376 | { return fluidState_.pressure(phaseIdx); } |
| 315 | |||
| 316 | /*! | ||
| 317 | * \brief Returns temperature within the control volume in \f$[K]\f$. | ||
| 318 | * | ||
| 319 | * Note that we assume thermodynamic equilibrium, i.e. the | ||
| 320 | * temperature of the rock matrix and of all fluid phases are | ||
| 321 | * identical. | ||
| 322 | */ | ||
| 323 | Scalar temperature() const | ||
| 324 |
2/4✓ Branch 0 taken 9898304 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 9898304 times.
✗ Branch 3 not taken.
|
299810712 | { return fluidState_.temperature(/*phaseIdx=*/0); } |
| 325 | |||
| 326 | /*! | ||
| 327 | * \brief Returns the relative permeability of a given phase within | ||
| 328 | * the control volume in \f$[-]\f$. | ||
| 329 | * | ||
| 330 | * \param phaseIdx The phase index | ||
| 331 | */ | ||
| 332 | Scalar relativePermeability(const int phaseIdx) const | ||
| 333 | { return relativePermeability_[phaseIdx]; } | ||
| 334 | |||
| 335 | /*! | ||
| 336 | * \brief Returns the effective mobility of a given phase within | ||
| 337 | * the control volume in \f$[s*m/kg]\f$. | ||
| 338 | * | ||
| 339 | * \param phaseIdx The phase index | ||
| 340 | */ | ||
| 341 | Scalar mobility(const int phaseIdx) const | ||
| 342 | 1631188848 | { return relativePermeability_[phaseIdx]/fluidState_.viscosity(phaseIdx); } | |
| 343 | |||
| 344 | /*! | ||
| 345 | * \brief Returns the effective capillary pressure within the control volume | ||
| 346 | * in \f$[kg/(m*s^2)=N/m^2=Pa]\f$. | ||
| 347 | */ | ||
| 348 | ✗ | Scalar capillaryPressure() const | |
| 349 | ✗ | { return pc_; } | |
| 350 | |||
| 351 | /*! | ||
| 352 | * \brief Returns the average porosity within the control volume in \f$[-]\f$. | ||
| 353 | */ | ||
| 354 | Scalar porosity() const | ||
| 355 |
6/6✓ Branch 0 taken 32 times.
✓ Branch 1 taken 2721120 times.
✓ Branch 2 taken 48 times.
✓ Branch 3 taken 27703340 times.
✓ Branch 4 taken 16 times.
✓ Branch 5 taken 24982220 times.
|
891190676 | { return solidState_.porosity(); } |
| 356 | |||
| 357 | /*! | ||
| 358 | * \brief Returns the average permeability within the control volume in \f$[m^2]\f$. | ||
| 359 | */ | ||
| 360 | const PermeabilityType& permeability() const | ||
| 361 | ✗ | { return permeability_; } | |
| 362 | |||
| 363 | /*! | ||
| 364 | * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$. | ||
| 365 | */ | ||
| 366 | Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const | ||
| 367 | { | ||
| 368 | typename FluidSystem::ParameterCache paramCache; | ||
| 369 | 29002408 | paramCache.updatePhase(fluidState_, phaseIdx); | |
| 370 |
4/5✓ Branch 1 taken 48 times.
✓ Branch 2 taken 26197924 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 716880 times.
✓ Branch 5 taken 788536 times.
|
29002408 | return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx); |
| 371 | } | ||
| 372 | |||
| 373 | /*! | ||
| 374 | * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$. | ||
| 375 | */ | ||
| 376 | Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const | ||
| 377 | 122154812 | { return effectiveDiffCoeff_(phaseIdx, compIIdx, compJIdx); } | |
| 378 | |||
| 379 | /*! | ||
| 380 | * \brief Returns the wetting phase index | ||
| 381 | */ | ||
| 382 | int wettingPhase() const | ||
| 383 | { return fluidState_.wettingPhase(); } | ||
| 384 | |||
| 385 | private: | ||
| 386 | FluidState fluidState_; | ||
| 387 | SolidState solidState_; | ||
| 388 | |||
| 389 | Scalar pc_; // The capillary pressure | ||
| 390 | PermeabilityType permeability_; // Effective permeability within the control volume | ||
| 391 | |||
| 392 | // Relative permeability within the control volume | ||
| 393 | std::array<Scalar, ModelTraits::numFluidPhases()> relativePermeability_; | ||
| 394 | |||
| 395 | // Effective diffusion coefficients for the phases | ||
| 396 | DiffusionCoefficients effectiveDiffCoeff_; | ||
| 397 | |||
| 398 | protected: | ||
| 399 | const Impl &asImp_() const { return *static_cast<const Impl*>(this); } | ||
| 400 | Impl &asImp_() { return *static_cast<Impl*>(this); } | ||
| 401 | }; | ||
| 402 | /*! | ||
| 403 | * \ingroup TwoPTwoCModel | ||
| 404 | * \brief Contains the quantities which are constant within a | ||
| 405 | * finite volume in the two-phase two-component model. | ||
| 406 | * Specialization for chemical equilibrium | ||
| 407 | */ | ||
| 408 | template <class Traits, bool useConstraintSolver> | ||
| 409 |
7/12✗ Branch 0 not taken.
✓ Branch 1 taken 8132 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2040 times.
✓ Branch 4 taken 6092 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 13056 times.
✓ Branch 7 taken 2432 times.
✓ Branch 8 taken 13056 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2432 times.
✗ Branch 11 not taken.
|
1153320 | class TwoPTwoCVolumeVariablesImplementation<Traits, false, useConstraintSolver> |
| 410 | : public TwoPTwoCVolumeVariablesBase<Traits, TwoPTwoCVolumeVariablesImplementation<Traits, false, useConstraintSolver>> | ||
| 411 | { | ||
| 412 | using ParentType = TwoPTwoCVolumeVariablesBase< Traits, TwoPTwoCVolumeVariablesImplementation<Traits, false, useConstraintSolver>>; | ||
| 413 | using EnergyVolVars = EnergyVolumeVariables<Traits, TwoPTwoCVolumeVariables<Traits> >; | ||
| 414 | |||
| 415 | using Scalar = typename Traits::PrimaryVariables::value_type; | ||
| 416 | using ModelTraits = typename Traits::ModelTraits; | ||
| 417 | |||
| 418 | static constexpr int numFluidComps = ParentType::numFluidComponents(); | ||
| 419 | // component indices | ||
| 420 | enum | ||
| 421 | { | ||
| 422 | comp0Idx = Traits::FluidSystem::comp0Idx, | ||
| 423 | comp1Idx = Traits::FluidSystem::comp1Idx, | ||
| 424 | phase0Idx = Traits::FluidSystem::phase0Idx, | ||
| 425 | phase1Idx = Traits::FluidSystem::phase1Idx | ||
| 426 | }; | ||
| 427 | |||
| 428 | // phase presence indices | ||
| 429 | enum | ||
| 430 | { | ||
| 431 | firstPhaseOnly = ModelTraits::Indices::firstPhaseOnly, | ||
| 432 | secondPhaseOnly = ModelTraits::Indices::secondPhaseOnly, | ||
| 433 | bothPhases = ModelTraits::Indices::bothPhases | ||
| 434 | }; | ||
| 435 | |||
| 436 | // primary variable indices | ||
| 437 | enum | ||
| 438 | { | ||
| 439 | switchIdx = ModelTraits::Indices::switchIdx, | ||
| 440 | pressureIdx = ModelTraits::Indices::pressureIdx | ||
| 441 | }; | ||
| 442 | |||
| 443 | // formulations | ||
| 444 | static constexpr auto formulation = ModelTraits::priVarFormulation(); | ||
| 445 | |||
| 446 | using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, typename Traits::FluidSystem>; | ||
| 447 | using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition< Scalar, typename Traits::FluidSystem >; | ||
| 448 | public: | ||
| 449 | //! The type of the object returned by the fluidState() method | ||
| 450 | using FluidState = typename Traits::FluidState; | ||
| 451 | //! The fluid system used here | ||
| 452 | using FluidSystem = typename Traits::FluidSystem; | ||
| 453 | //! Export type of solid state | ||
| 454 | using SolidState = typename Traits::SolidState; | ||
| 455 | //! Export type of solid system | ||
| 456 | using SolidSystem = typename Traits::SolidSystem; | ||
| 457 | //! Export the primary variable switch | ||
| 458 | using PrimaryVariableSwitch = TwoPNCPrimaryVariableSwitch; | ||
| 459 | |||
| 460 | //! Return whether moles or masses are balanced | ||
| 461 | static constexpr bool useMoles() { return ModelTraits::useMoles(); } | ||
| 462 | //! Return the two-phase formulation used here | ||
| 463 | static constexpr TwoPFormulation priVarFormulation() { return formulation; } | ||
| 464 | |||
| 465 | // check for permissive combinations | ||
| 466 | static_assert(useMoles() || (!useMoles() && useConstraintSolver), "if !UseMoles, UseConstraintSolver has to be set to true"); | ||
| 467 | |||
| 468 | // The computations in the explicit composition update most probably assume a liquid-gas interface with | ||
| 469 | // liquid as first phase. TODO: is this really needed? The constraint solver does the job anyway, doesn't it? | ||
| 470 | static_assert(useConstraintSolver || (!FluidSystem::isGas(phase0Idx) && FluidSystem::isGas(phase1Idx)), | ||
| 471 | "Explicit composition calculation has to be re-checked for NON-liquid-gas equilibria"); | ||
| 472 | |||
| 473 | /*! | ||
| 474 | * \brief Sets complete fluid state. | ||
| 475 | * | ||
| 476 | * \param elemSol A vector containing all primary variables connected to the element | ||
| 477 | * \param problem The object specifying the problem which ought to | ||
| 478 | * be simulated | ||
| 479 | * \param element An element which contains part of the control volume | ||
| 480 | * \param scv The sub-control volume | ||
| 481 | * \param fluidState A container with the current (physical) state of the fluid | ||
| 482 | * \param solidState A container with the current (physical) state of the solid | ||
| 483 | * | ||
| 484 | * Set temperature, saturations, capillary pressures, viscosities, densities and enthalpies. | ||
| 485 | */ | ||
| 486 | template<class ElemSol, class Problem, class Element, class Scv> | ||
| 487 | 13457426 | void completeFluidState(const ElemSol& elemSol, | |
| 488 | const Problem& problem, | ||
| 489 | const Element& element, | ||
| 490 | const Scv& scv, | ||
| 491 | FluidState& fluidState, | ||
| 492 | SolidState& solidState) | ||
| 493 | { | ||
| 494 | 13457426 | ParentType::completeFluidState(elemSol, problem, element, scv, fluidState, solidState); | |
| 495 | |||
| 496 | 13457426 | const auto& priVars = elemSol[scv.localDofIndex()]; | |
| 497 | 13457426 | const auto phasePresence = priVars.state(); | |
| 498 | |||
| 499 | // calculate the phase compositions | ||
| 500 | typename FluidSystem::ParameterCache paramCache; | ||
| 501 | |||
| 502 | // If constraint solver is not used, get the phase pressures and set the fugacity coefficients here | ||
| 503 | if(!useConstraintSolver) | ||
| 504 | { | ||
| 505 | for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++ phaseIdx) | ||
| 506 | { | ||
| 507 | assert(FluidSystem::isIdealMixture(phaseIdx)); | ||
| 508 | for (int compIdx = 0; compIdx < ModelTraits::numFluidComponents(); ++ compIdx) { | ||
| 509 | Scalar phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx); | ||
| 510 | fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi); | ||
| 511 | } | ||
| 512 | } | ||
| 513 | } | ||
| 514 | |||
| 515 | // now comes the tricky part: calculate phase compositions | ||
| 516 |
2/2✓ Branch 0 taken 1255860 times.
✓ Branch 1 taken 12201566 times.
|
13457426 | const Scalar p0 = fluidState.pressure(phase0Idx); |
| 517 |
2/2✓ Branch 0 taken 1255860 times.
✓ Branch 1 taken 12201566 times.
|
13457426 | const Scalar p1 = fluidState.pressure(phase1Idx); |
| 518 |
2/2✓ Branch 0 taken 1255860 times.
✓ Branch 1 taken 12201566 times.
|
13457426 | if (phasePresence == bothPhases) |
| 519 | { | ||
| 520 | // both phases are present, phase compositions are a result | ||
| 521 | // of the equilibrium between the phases. This is the job | ||
| 522 | // of the "MiscibleMultiPhaseComposition" constraint solver | ||
| 523 | if(useConstraintSolver) | ||
| 524 | 1255860 | MiscibleMultiPhaseComposition::solve(fluidState, | |
| 525 | paramCache); | ||
| 526 | // ... or calculated explicitly this way ... | ||
| 527 | else | ||
| 528 | { | ||
| 529 | // get the partial pressure of the main component of the first phase within the | ||
| 530 | // second phase == vapor pressure due to equilibrium. Note that in this case the | ||
| 531 | // fugacityCoefficient * p is the vapor pressure (see implementation in respective fluidsystem) | ||
| 532 | const Scalar partPressLiquid = FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx)*p0; | ||
| 533 | |||
| 534 | // get the partial pressure of the main component of the gas phase | ||
| 535 | const Scalar partPressGas = p1 - partPressLiquid; | ||
| 536 | |||
| 537 | // calculate the mole fractions of the components within the nonwetting phase | ||
| 538 | const Scalar xnn = partPressGas / p1; | ||
| 539 | const Scalar xnw = partPressLiquid / p1; | ||
| 540 | |||
| 541 | // calculate the mole fractions of the components within the wetting phase | ||
| 542 | // note that in this case the fugacityCoefficient * p is the Henry Coefficient | ||
| 543 | // (see implementation in respective fluidsystem) | ||
| 544 | const Scalar xwn = partPressGas / (FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx)*p0); | ||
| 545 | const Scalar xww = 1.0 - xwn; | ||
| 546 | |||
| 547 | // set all mole fractions | ||
| 548 | fluidState.setMoleFraction(phase0Idx, comp0Idx, xww); | ||
| 549 | fluidState.setMoleFraction(phase0Idx, comp1Idx, xwn); | ||
| 550 | fluidState.setMoleFraction(phase1Idx, comp0Idx, xnw); | ||
| 551 | fluidState.setMoleFraction(phase1Idx, comp1Idx, xnn); | ||
| 552 | } | ||
| 553 | } | ||
| 554 |
2/2✓ Branch 0 taken 17974 times.
✓ Branch 1 taken 12183592 times.
|
12201566 | else if (phasePresence == secondPhaseOnly) |
| 555 | { | ||
| 556 | // only the second phase is present, composition is stored explicitly. | ||
| 557 | if( useMoles() ) | ||
| 558 | { | ||
| 559 | 35948 | fluidState.setMoleFraction(phase1Idx, comp1Idx, 1 - priVars[switchIdx]); | |
| 560 | 35948 | fluidState.setMoleFraction(phase1Idx, comp0Idx, priVars[switchIdx]); | |
| 561 | } | ||
| 562 | // setMassFraction() has only to be called 1-numComponents times | ||
| 563 | else | ||
| 564 | fluidState.setMassFraction(phase1Idx, comp0Idx, priVars[switchIdx]); | ||
| 565 | |||
| 566 | // calculate the composition of the remaining phases (as | ||
| 567 | // well as the densities of all phases). This is the job | ||
| 568 | // of the "ComputeFromReferencePhase" constraint solver | ||
| 569 | if (useConstraintSolver) | ||
| 570 | 17974 | ComputeFromReferencePhase::solve(fluidState, | |
| 571 | paramCache, | ||
| 572 | phase1Idx); | ||
| 573 | // ... or calculated explicitly this way ... | ||
| 574 | else | ||
| 575 | { | ||
| 576 | // note that the water phase is actually not existing! | ||
| 577 | // thus, this is used as phase switch criterion | ||
| 578 | const Scalar xnw = priVars[switchIdx]; | ||
| 579 | const Scalar xnn = 1.0 - xnw; | ||
| 580 | |||
| 581 | // first, xww: | ||
| 582 | // xnw * pn = "actual" (hypothetical) vapor pressure | ||
| 583 | // fugacityCoefficient * pw = vapor pressure given by thermodynamic conditions | ||
| 584 | // Here, xww is not actually the mole fraction of water in the wetting phase | ||
| 585 | // xww is only the ratio of "actual" vapor pressure / "thermodynamic" vapor pressure | ||
| 586 | // If xww > 1 : gas is over-saturated with water vapor, | ||
| 587 | // condensation takes place (see switch criterion in model) | ||
| 588 | const Scalar xww = xnw*p1/( FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx)*p0 ); | ||
| 589 | |||
| 590 | // second, xwn: | ||
| 591 | // partialPressure / xwn = Henry | ||
| 592 | // partialPressure = xnn * pn | ||
| 593 | // xwn = xnn * pn / Henry | ||
| 594 | // Henry = fugacityCoefficient * pw | ||
| 595 | const Scalar xwn = xnn*p1/( FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx)*p0 ); | ||
| 596 | |||
| 597 | fluidState.setMoleFraction(phase0Idx, comp0Idx, xww); | ||
| 598 | fluidState.setMoleFraction(phase0Idx, comp1Idx, xwn); | ||
| 599 | } | ||
| 600 | } | ||
| 601 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 12183592 times.
|
12183592 | else if (phasePresence == firstPhaseOnly) |
| 602 | { | ||
| 603 | // only the wetting phase is present, i.e. wetting phase | ||
| 604 | // composition is stored explicitly. | ||
| 605 | if( useMoles() ) // mole-fraction formulation | ||
| 606 | { | ||
| 607 | 24367184 | fluidState.setMoleFraction(phase0Idx, comp0Idx, 1-priVars[switchIdx]); | |
| 608 | 24367184 | fluidState.setMoleFraction(phase0Idx, comp1Idx, priVars[switchIdx]); | |
| 609 | } | ||
| 610 | // setMassFraction() has only to be called 1-numComponents times | ||
| 611 | else // mass-fraction formulation | ||
| 612 | fluidState.setMassFraction(phase0Idx, comp1Idx, priVars[switchIdx]); | ||
| 613 | |||
| 614 | // calculate the composition of the remaining phases (as | ||
| 615 | // well as the densities of all phases). This is the job | ||
| 616 | // of the "ComputeFromReferencePhase" constraint solver | ||
| 617 | if (useConstraintSolver) | ||
| 618 | 12183592 | ComputeFromReferencePhase::solve(fluidState, | |
| 619 | paramCache, | ||
| 620 | phase0Idx); | ||
| 621 | // ... or calculated explicitly this way ... | ||
| 622 | else | ||
| 623 | { | ||
| 624 | // note that the gas phase is actually not existing! | ||
| 625 | // thus, this is used as phase switch criterion | ||
| 626 | const Scalar xwn = priVars[switchIdx]; | ||
| 627 | |||
| 628 | // first, xnw: | ||
| 629 | // psteam = xnw * pn = partial pressure of water in gas phase | ||
| 630 | // psteam = fugacityCoefficient * pw | ||
| 631 | const Scalar xnw = ( FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp0Idx)*p0 )/p1; | ||
| 632 | |||
| 633 | // second, xnn: | ||
| 634 | // xwn = partialPressure / Henry | ||
| 635 | // partialPressure = pn * xnn | ||
| 636 | // xwn = pn * xnn / Henry | ||
| 637 | // xnn = xwn * Henry / pn | ||
| 638 | // Henry = fugacityCoefficient * pw | ||
| 639 | const Scalar xnn = xwn*( FluidSystem::fugacityCoefficient(fluidState, phase0Idx, comp1Idx)*p0 )/p1; | ||
| 640 | |||
| 641 | fluidState.setMoleFraction(phase1Idx, comp1Idx, xnn); | ||
| 642 | fluidState.setMoleFraction(phase1Idx, comp0Idx, xnw); | ||
| 643 | } | ||
| 644 | } | ||
| 645 | |||
| 646 |
2/2✓ Branch 0 taken 26914852 times.
✓ Branch 1 taken 13457426 times.
|
40372278 | for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx) |
| 647 | { | ||
| 648 | // set the viscosity and desity here if constraintsolver is not used | ||
| 649 | if(!useConstraintSolver) | ||
| 650 | { | ||
| 651 | paramCache.updateComposition(fluidState, phaseIdx); | ||
| 652 | const Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx); | ||
| 653 | fluidState.setDensity(phaseIdx, rho); | ||
| 654 | Scalar rhoMolar = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx); | ||
| 655 | fluidState.setMolarDensity(phaseIdx, rhoMolar); | ||
| 656 | } | ||
| 657 | |||
| 658 | // compute and set the enthalpy | ||
| 659 | 26914852 | const Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx); | |
| 660 | 26914852 | fluidState.setViscosity(phaseIdx,mu); | |
| 661 | 26914852 | Scalar h = EnergyVolVars::enthalpy(fluidState, paramCache, phaseIdx); | |
| 662 | 53829704 | fluidState.setEnthalpy(phaseIdx, h); | |
| 663 | } | ||
| 664 | 13457426 | } | |
| 665 | |||
| 666 | }; | ||
| 667 | /*! | ||
| 668 | * \ingroup TwoPTwoCModel | ||
| 669 | * \brief Contains the quantities which are constant within a | ||
| 670 | * finite volume in the two-phase two-component model. | ||
| 671 | * Specialization for chemical non-equilibrium. | ||
| 672 | * The equilibrium mole fraction is calculated using Henry's and Raoult's law | ||
| 673 | */ | ||
| 674 | template <class Traits, bool useConstraintSolver> | ||
| 675 | 173272 | class TwoPTwoCVolumeVariablesImplementation<Traits, true, useConstraintSolver> | |
| 676 | : public TwoPTwoCVolumeVariablesBase<Traits, TwoPTwoCVolumeVariablesImplementation<Traits, true, useConstraintSolver>> | ||
| 677 | { | ||
| 678 | using ParentType = TwoPTwoCVolumeVariablesBase< Traits, TwoPTwoCVolumeVariablesImplementation<Traits, true, useConstraintSolver>>; | ||
| 679 | using EnergyVolVars = EnergyVolumeVariables<Traits, TwoPTwoCVolumeVariables<Traits> >; | ||
| 680 | |||
| 681 | using Scalar = typename Traits::PrimaryVariables::value_type; | ||
| 682 | using ModelTraits = typename Traits::ModelTraits; | ||
| 683 | |||
| 684 | static constexpr int numFluidComps = ParentType::numFluidComponents(); | ||
| 685 | // component indices | ||
| 686 | enum | ||
| 687 | { | ||
| 688 | comp0Idx = Traits::FluidSystem::comp0Idx, | ||
| 689 | comp1Idx = Traits::FluidSystem::comp1Idx, | ||
| 690 | phase0Idx = Traits::FluidSystem::phase0Idx, | ||
| 691 | phase1Idx = Traits::FluidSystem::phase1Idx | ||
| 692 | }; | ||
| 693 | |||
| 694 | // phase presence indices | ||
| 695 | enum | ||
| 696 | { | ||
| 697 | firstPhaseOnly = ModelTraits::Indices::firstPhaseOnly, | ||
| 698 | secondPhaseOnly = ModelTraits::Indices::secondPhaseOnly, | ||
| 699 | bothPhases = ModelTraits::Indices::bothPhases | ||
| 700 | }; | ||
| 701 | |||
| 702 | // primary variable indices | ||
| 703 | enum | ||
| 704 | { | ||
| 705 | switchIdx = ModelTraits::Indices::switchIdx, | ||
| 706 | pressureIdx = ModelTraits::Indices::pressureIdx | ||
| 707 | }; | ||
| 708 | |||
| 709 | // formulations | ||
| 710 | static constexpr auto formulation = ModelTraits::priVarFormulation(); | ||
| 711 | |||
| 712 | using PermeabilityType = typename Traits::PermeabilityType; | ||
| 713 | using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, typename Traits::FluidSystem>; | ||
| 714 | using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition< Scalar, typename Traits::FluidSystem >; | ||
| 715 | |||
| 716 | public: | ||
| 717 | //! The type of the object returned by the fluidState() method | ||
| 718 | using FluidState = typename Traits::FluidState; | ||
| 719 | //! The fluid system used here | ||
| 720 | using FluidSystem = typename Traits::FluidSystem; | ||
| 721 | //! Export type of solid state | ||
| 722 | using SolidState = typename Traits::SolidState; | ||
| 723 | //! Export type of solid system | ||
| 724 | using SolidSystem = typename Traits::SolidSystem; | ||
| 725 | //! Export the primary variable switch | ||
| 726 | using PrimaryVariableSwitch = TwoPNCPrimaryVariableSwitch; | ||
| 727 | |||
| 728 | //! Return whether moles or masses are balanced | ||
| 729 | static constexpr bool useMoles() { return ModelTraits::useMoles(); } | ||
| 730 | //! Return the two-phase formulation used here | ||
| 731 | static constexpr TwoPFormulation priVarFormulation() { return formulation; } | ||
| 732 | |||
| 733 | // check for permissive combinations | ||
| 734 | static_assert(useMoles() || (!useMoles() && useConstraintSolver), "if !UseMoles, UseConstraintSolver has to be set to true"); | ||
| 735 | static_assert((formulation == TwoPFormulation::p0s1 || formulation == TwoPFormulation::p1s0), "Chosen TwoPFormulation not supported!"); | ||
| 736 | |||
| 737 | // The computations in the explicit composition update most probably assume a liquid-gas interface with | ||
| 738 | // liquid as first phase. TODO: is this really needed? The constraint solver does the job anyway, doesn't it? | ||
| 739 | static_assert(useConstraintSolver || (!FluidSystem::isGas(phase0Idx) && FluidSystem::isGas(phase1Idx)), | ||
| 740 | "Explicit composition calculation has to be re-checked for NON-liquid-gas equilibria"); | ||
| 741 | |||
| 742 | /*! | ||
| 743 | * \brief Sets complete fluid state. | ||
| 744 | * | ||
| 745 | * \param elemSol A vector containing all primary variables connected to the element | ||
| 746 | * \param problem The object specifying the problem which ought to | ||
| 747 | * be simulated | ||
| 748 | * \param element An element which contains part of the control volume | ||
| 749 | * \param scv The sub-control volume | ||
| 750 | * \param fluidState A container with the current (physical) state of the fluid | ||
| 751 | * \param solidState A container with the current (physical) state of the solid | ||
| 752 | * | ||
| 753 | * Set temperature, saturations, capillary pressures, viscosities, densities and enthalpies. | ||
| 754 | */ | ||
| 755 | template<class ElemSol, class Problem, class Element, class Scv> | ||
| 756 | 394268 | void completeFluidState(const ElemSol& elemSol, | |
| 757 | const Problem& problem, | ||
| 758 | const Element& element, | ||
| 759 | const Scv& scv, | ||
| 760 | FluidState& fluidState, | ||
| 761 | SolidState& solidState) | ||
| 762 | { | ||
| 763 | 394268 | ParentType::completeFluidState(elemSol, problem, element, scv, fluidState, solidState); | |
| 764 | |||
| 765 | 394268 | const auto& priVars = elemSol[scv.localDofIndex()]; | |
| 766 | |||
| 767 | ///////////// | ||
| 768 | // set the fluid compositions | ||
| 769 | ///////////// | ||
| 770 | typename FluidSystem::ParameterCache paramCache; | ||
| 771 | 394268 | paramCache.updateAll(fluidState); | |
| 772 | |||
| 773 | 394268 | updateMoleFraction(fluidState, | |
| 774 | paramCache, | ||
| 775 | priVars); | ||
| 776 | |||
| 777 | |||
| 778 |
2/2✓ Branch 0 taken 788536 times.
✓ Branch 1 taken 394268 times.
|
1182804 | for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++phaseIdx) |
| 779 | { | ||
| 780 | // compute and set the enthalpy | ||
| 781 | 788536 | const Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx); | |
| 782 | 788536 | fluidState.setViscosity(phaseIdx,mu); | |
| 783 | 788536 | Scalar h = EnergyVolVars::enthalpy(fluidState, paramCache, phaseIdx); | |
| 784 | 1577072 | fluidState.setEnthalpy(phaseIdx, h); | |
| 785 | } | ||
| 786 | 394268 | } | |
| 787 | |||
| 788 | /*! | ||
| 789 | * \brief Updates composition of all phases from the primary variables. | ||
| 790 | * | ||
| 791 | * \param actualFluidState Container for all the secondary variables concerning the fluids | ||
| 792 | * \param paramCache Container for cache parameters | ||
| 793 | * \param priVars The primary Variables | ||
| 794 | */ | ||
| 795 | 394268 | void updateMoleFraction(FluidState &actualFluidState, | |
| 796 | typename Traits::FluidSystem::ParameterCache & paramCache, | ||
| 797 | const typename Traits::PrimaryVariables& priVars) | ||
| 798 | { | ||
| 799 | 394268 | const auto phasePresence = priVars.state(); | |
| 800 | |||
| 801 | 394268 | Scalar xwnNonEquil = 0.0; | |
| 802 | 394268 | Scalar xwwNonEquil = 0.0; | |
| 803 | 394268 | Scalar xnwNonEquil = 0.0; | |
| 804 | 394268 | Scalar xnnNonEquil = 0.0; | |
| 805 | |||
| 806 |
1/2✓ Branch 0 taken 394268 times.
✗ Branch 1 not taken.
|
394268 | if (phasePresence == bothPhases) |
| 807 | { | ||
| 808 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 394268 times.
|
394268 | xwnNonEquil = priVars[ModelTraits::numFluidPhases()]; |
| 809 | 394268 | xwwNonEquil = 1-xwnNonEquil; | |
| 810 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 394268 times.
|
394268 | xnwNonEquil = priVars[ModelTraits::numFluidPhases()+comp1Idx]; |
| 811 | |||
| 812 | //if the capillary pressure is high, we need to use Kelvin's equation to reduce the vapor pressure | ||
| 813 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 394268 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 394268 times.
|
788536 | if (actualFluidState.saturation(phase0Idx) < 0.01) |
| 814 | { | ||
| 815 | ✗ | const Scalar p1 = actualFluidState.pressure(phase1Idx); | |
| 816 | ✗ | const Scalar partPressLiquid = FluidSystem::vaporPressure(actualFluidState, comp0Idx); | |
| 817 | ✗ | xnwNonEquil =std::min(partPressLiquid/p1, xnwNonEquil); | |
| 818 | } | ||
| 819 | 394268 | xnnNonEquil = 1- xnwNonEquil; | |
| 820 | |||
| 821 | 394268 | actualFluidState.setMoleFraction(phase0Idx, comp0Idx, xwwNonEquil); | |
| 822 | 394268 | actualFluidState.setMoleFraction(phase0Idx, comp1Idx, xwnNonEquil); | |
| 823 | 394268 | actualFluidState.setMoleFraction(phase1Idx, comp0Idx, xnwNonEquil); | |
| 824 | 394268 | actualFluidState.setMoleFraction(phase1Idx, comp1Idx, xnnNonEquil); | |
| 825 | |||
| 826 | //check if we need this | ||
| 827 |
2/2✓ Branch 0 taken 788536 times.
✓ Branch 1 taken 394268 times.
|
1182804 | for(int phaseIdx=0; phaseIdx<ModelTraits::numFluidPhases(); ++phaseIdx) |
| 828 |
2/2✓ Branch 0 taken 1577072 times.
✓ Branch 1 taken 788536 times.
|
2365608 | for (int compIdx = 0; compIdx < numFluidComps; ++compIdx) |
| 829 | { | ||
| 830 | 1577072 | const Scalar phi = FluidSystem::fugacityCoefficient(actualFluidState, | |
| 831 | paramCache, | ||
| 832 | phaseIdx, | ||
| 833 | compIdx); | ||
| 834 | 3154144 | actualFluidState.setFugacityCoefficient(phaseIdx, | |
| 835 | compIdx, | ||
| 836 | phi); | ||
| 837 | } | ||
| 838 | |||
| 839 | 394268 | FluidState equilFluidState; // the fluidState *on the interface* i.e. chemical equilibrium | |
| 840 | 394268 | equilFluidState.assign(actualFluidState) ; | |
| 841 | // If constraint solver is not used, get the phase pressures and set the fugacity coefficients here | ||
| 842 | if(!useConstraintSolver) | ||
| 843 | { | ||
| 844 | for (int phaseIdx = 0; phaseIdx < ModelTraits::numFluidPhases(); ++ phaseIdx) | ||
| 845 | { | ||
| 846 | assert(FluidSystem::isIdealMixture(phaseIdx)); | ||
| 847 | for (int compIdx = 0; compIdx < ModelTraits::numFluidComponents(); ++ compIdx) { | ||
| 848 | Scalar phi = FluidSystem::fugacityCoefficient(equilFluidState, paramCache, phaseIdx, compIdx); | ||
| 849 | equilFluidState.setFugacityCoefficient(phaseIdx, compIdx, phi); | ||
| 850 | } | ||
| 851 | } | ||
| 852 | } | ||
| 853 | |||
| 854 | // now comes the tricky part: calculate phase compositions | ||
| 855 | 394268 | const Scalar p0 = equilFluidState.pressure(phase0Idx); | |
| 856 | 394268 | const Scalar p1 = equilFluidState.pressure(phase1Idx); | |
| 857 | // both phases are present, phase compositions are a result | ||
| 858 | // of the equilibrium between the phases. This is the job | ||
| 859 | // of the "MiscibleMultiPhaseComposition" constraint solver | ||
| 860 | if(useConstraintSolver) | ||
| 861 | 394268 | MiscibleMultiPhaseComposition::solve(equilFluidState, | |
| 862 | paramCache); | ||
| 863 | // ... or calculated explicitly this way ... | ||
| 864 | else | ||
| 865 | { | ||
| 866 | // get the partial pressure of the main component of the first phase within the | ||
| 867 | // second phase == vapor pressure due to equilibrium. Note that in this case the | ||
| 868 | // fugacityCoefficient * p is the vapor pressure (see implementation in respective fluidsystem) | ||
| 869 | const Scalar partPressLiquid = FluidSystem::fugacityCoefficient(equilFluidState, phase0Idx, comp0Idx)*p0; | ||
| 870 | |||
| 871 | // get the partial pressure of the main component of the gas phase | ||
| 872 | const Scalar partPressGas = p1 - partPressLiquid; | ||
| 873 | |||
| 874 | // calculate the mole fractions of the components within the nonwetting phase | ||
| 875 | const Scalar xnn = partPressGas / p1; | ||
| 876 | const Scalar xnw = partPressLiquid / p1; | ||
| 877 | |||
| 878 | // calculate the mole fractions of the components within the wetting phase | ||
| 879 | // note that in this case the fugacityCoefficient * p is the Henry Coefficient | ||
| 880 | // (see implementation in respective fluidsystem) | ||
| 881 | const Scalar xwn = partPressGas / (FluidSystem::fugacityCoefficient(equilFluidState, phase0Idx, comp1Idx)*p0); | ||
| 882 | const Scalar xww = 1.0 - xwn; | ||
| 883 | |||
| 884 | // set all mole fractions | ||
| 885 | equilFluidState.setMoleFraction(phase0Idx, comp0Idx, xww); | ||
| 886 | equilFluidState.setMoleFraction(phase0Idx, comp1Idx, xwn); | ||
| 887 | equilFluidState.setMoleFraction(phase1Idx, comp0Idx, xnw); | ||
| 888 | equilFluidState.setMoleFraction(phase1Idx, comp1Idx, xnn); | ||
| 889 | } | ||
| 890 | |||
| 891 | // Setting the equilibrium composition (in a kinetic model not necessarily the same as the actual mole fraction) | ||
| 892 |
2/2✓ Branch 0 taken 788536 times.
✓ Branch 1 taken 394268 times.
|
1182804 | for(int phaseIdx=0; phaseIdx<ModelTraits::numFluidPhases(); ++phaseIdx){ |
| 893 |
2/2✓ Branch 0 taken 1577072 times.
✓ Branch 1 taken 788536 times.
|
2365608 | for (int compIdx=0; compIdx< numFluidComps; ++ compIdx){ |
| 894 | 6308288 | xEquil_[phaseIdx][compIdx] = equilFluidState.moleFraction(phaseIdx, compIdx); | |
| 895 | } | ||
| 896 | } | ||
| 897 | } | ||
| 898 | else | ||
| 899 | { | ||
| 900 | ✗ | DUNE_THROW(Dune::InvalidStateException, "nonequilibrium is only possible for 2 phases present "); | |
| 901 | } | ||
| 902 | // compute densities of all phases | ||
| 903 |
2/2✓ Branch 0 taken 788536 times.
✓ Branch 1 taken 394268 times.
|
1182804 | for(int phaseIdx=0; phaseIdx<ModelTraits::numFluidPhases(); ++phaseIdx) |
| 904 | { | ||
| 905 | 788536 | const Scalar rho = FluidSystem::density(actualFluidState, paramCache, phaseIdx); | |
| 906 | 788536 | actualFluidState.setDensity(phaseIdx, rho); | |
| 907 | 788536 | const Scalar rhoMolar = FluidSystem::molarDensity(actualFluidState, paramCache, phaseIdx); | |
| 908 | 1577072 | actualFluidState.setMolarDensity(phaseIdx, rhoMolar); | |
| 909 | } | ||
| 910 | 394268 | } | |
| 911 | |||
| 912 | /*! | ||
| 913 | * \brief The mole fraction we would have in the case of chemical equilibrium / | ||
| 914 | * on the interface. | ||
| 915 | * | ||
| 916 | * \param phaseIdx The index of the fluid phase | ||
| 917 | * \param compIdx The local index of the component | ||
| 918 | */ | ||
| 919 | const Scalar xEquil(const unsigned int phaseIdx, const unsigned int compIdx) const | ||
| 920 | { | ||
| 921 | 3820800 | return xEquil_[phaseIdx][compIdx] ; | |
| 922 | } | ||
| 923 | |||
| 924 | private: | ||
| 925 | std::array<std::array<Scalar, numFluidComps>, ModelTraits::numFluidPhases()> xEquil_; | ||
| 926 | }; | ||
| 927 | |||
| 928 | } // end namespace Dumux | ||
| 929 | |||
| 930 | #endif | ||
| 931 |