| 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 MPNCModel | ||
| 10 | * \brief Contains the secondary variables (Quantities which are | ||
| 11 | * constant within a finite volume) of the MpNc model. | ||
| 12 | */ | ||
| 13 | |||
| 14 | #ifndef DUMUX_MPNC_VOLUME_VARIABLES_HH | ||
| 15 | #define DUMUX_MPNC_VOLUME_VARIABLES_HH | ||
| 16 | |||
| 17 | #include <dumux/porousmediumflow/volumevariables.hh> | ||
| 18 | #include <dumux/porousmediumflow/nonisothermal/volumevariables.hh> | ||
| 19 | |||
| 20 | #include <dumux/material/constraintsolvers/compositionfromfugacities.hh> | ||
| 21 | #include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh> | ||
| 22 | #include <dumux/material/solidstates/updatesolidvolumefractions.hh> | ||
| 23 | |||
| 24 | #include "pressureformulation.hh" | ||
| 25 | |||
| 26 | namespace Dumux { | ||
| 27 | |||
| 28 | // forward declaration | ||
| 29 | template <class Traits, bool enableChemicalNonEquilibrium> | ||
| 30 | class MPNCVolumeVariablesImplementation; | ||
| 31 | /*! | ||
| 32 | * \ingroup MPNCModel | ||
| 33 | * \brief Contains the quantities which are constant within a finite volume in the MpNc model. | ||
| 34 | * | ||
| 35 | * \tparam Traits Class encapsulating types to be used by the vol vars | ||
| 36 | */ | ||
| 37 | template <class Traits> | ||
| 38 | using MPNCVolumeVariables = MPNCVolumeVariablesImplementation<Traits, Traits::ModelTraits::enableChemicalNonEquilibrium()>; | ||
| 39 | |||
| 40 | template <class Traits> | ||
| 41 | 1438920 | class MPNCVolumeVariablesImplementation<Traits, false> | |
| 42 | : public PorousMediumFlowVolumeVariables<Traits> | ||
| 43 | , public EnergyVolumeVariables<Traits, MPNCVolumeVariables<Traits> > | ||
| 44 | { | ||
| 45 | using ParentType = PorousMediumFlowVolumeVariables<Traits>; | ||
| 46 | using EnergyVolVars = EnergyVolumeVariables<Traits, MPNCVolumeVariables<Traits> >; | ||
| 47 | |||
| 48 | using Scalar = typename Traits::PrimaryVariables::value_type; | ||
| 49 | using PermeabilityType = typename Traits::PermeabilityType; | ||
| 50 | |||
| 51 | using ModelTraits = typename Traits::ModelTraits; | ||
| 52 | static constexpr auto pressureFormulation = ModelTraits::pressureFormulation(); | ||
| 53 | |||
| 54 | static constexpr bool enableThermalNonEquilibrium = ModelTraits::enableThermalNonEquilibrium(); | ||
| 55 | static constexpr bool enableChemicalNonEquilibrium = ModelTraits::enableChemicalNonEquilibrium(); | ||
| 56 | static constexpr bool enableDiffusion = ModelTraits::enableMolecularDiffusion(); | ||
| 57 | |||
| 58 | using ComponentVector = Dune::FieldVector<Scalar, ModelTraits::numFluidComponents()>; | ||
| 59 | using CompositionFromFugacities = Dumux::CompositionFromFugacities<Scalar, typename Traits::FluidSystem>; | ||
| 60 | using EffDiffModel = typename Traits::EffectiveDiffusivityModel; | ||
| 61 | using DiffusionCoefficients = typename Traits::DiffusionType::DiffusionCoefficientsContainer; | ||
| 62 | |||
| 63 | public: | ||
| 64 | //! Export the type encapsulating primary variable indices | ||
| 65 | using Indices = typename Traits::ModelTraits::Indices; | ||
| 66 | //! Export the underlying fluid system | ||
| 67 | using FluidSystem = typename Traits::FluidSystem; | ||
| 68 | //! Export the fluid state type | ||
| 69 | using FluidState = typename Traits::FluidState; | ||
| 70 | //! Export type of solid state | ||
| 71 | using SolidState = typename Traits::SolidState; | ||
| 72 | //! Export type of solid system | ||
| 73 | using SolidSystem = typename Traits::SolidSystem; | ||
| 74 | |||
| 75 | //! Return number of phases considered by the model | ||
| 76 | static constexpr int numFluidPhases() { return ModelTraits::numFluidPhases(); } | ||
| 77 | //! Return number of components considered by the model | ||
| 78 | static constexpr int numFluidComps = ParentType::numFluidComponents(); | ||
| 79 | |||
| 80 | /*! | ||
| 81 | * \brief Updates all quantities for a given control volume. | ||
| 82 | * | ||
| 83 | * \param elemSol A vector containing all primary variables connected to the element | ||
| 84 | * \param problem The object specifying the problem which ought to | ||
| 85 | * be simulated | ||
| 86 | * \param element An element which contains part of the control volume | ||
| 87 | * \param scv The sub-control volume | ||
| 88 | */ | ||
| 89 | template<class ElemSol, class Problem, class Element, class Scv> | ||
| 90 | 1193800 | void update(const ElemSol& elemSol, | |
| 91 | const Problem& problem, | ||
| 92 | const Element& element, | ||
| 93 | const Scv& scv) | ||
| 94 | { | ||
| 95 | 1193800 | ParentType::update(elemSol, problem, element, scv); | |
| 96 | |||
| 97 | 1193800 | completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_); | |
| 98 | |||
| 99 | // calculate the remaining quantities | ||
| 100 | 1193800 | const auto& spatialParams = problem.spatialParams(); | |
| 101 | 1193800 | const int wPhaseIdx = spatialParams.template wettingPhase<FluidSystem>(element, scv, elemSol); | |
| 102 | 1193800 | const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol); | |
| 103 | 1193800 | const auto relPerm = fluidMatrixInteraction.relativePermeabilities(fluidState_, wPhaseIdx); | |
| 104 | 4775200 | std::copy(relPerm.begin(), relPerm.end(), relativePermeability_.begin()); | |
| 105 | |||
| 106 | typename FluidSystem::ParameterCache paramCache; | ||
| 107 | 1193800 | paramCache.updateAll(fluidState_); | |
| 108 | |||
| 109 | // porosity | ||
| 110 | 1193800 | updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps); | |
| 111 | |||
| 112 | if constexpr (enableDiffusion) | ||
| 113 | { | ||
| 114 | 1193800 | auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx) | |
| 115 | 2387600 | { return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx); }; | |
| 116 | |||
| 117 | 1193800 | effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient); | |
| 118 | } | ||
| 119 | |||
| 120 | 1193800 | EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_); | |
| 121 |
4/4✓ Branch 0 taken 809916 times.
✓ Branch 1 taken 383884 times.
✓ Branch 2 taken 809916 times.
✓ Branch 3 taken 383884 times.
|
2387600 | permeability_ = problem.spatialParams().permeability(element, scv, elemSol); |
| 122 | 1193800 | EnergyVolVars::updateEffectiveThermalConductivity(); | |
| 123 | 1193800 | } | |
| 124 | |||
| 125 | /*! | ||
| 126 | * \brief Sets complete fluid state. | ||
| 127 | * | ||
| 128 | * \param elemSol A vector containing all primary variables connected to the element | ||
| 129 | * \param problem The object specifying the problem which ought to | ||
| 130 | * be simulated | ||
| 131 | * \param element An element which contains part of the control volume | ||
| 132 | * \param scv The sub-control volume | ||
| 133 | * \param fluidState A container with the current (physical) state of the fluid | ||
| 134 | * \param solidState A container with the current (physical) state of the solid | ||
| 135 | */ | ||
| 136 | |||
| 137 | template<class ElemSol, class Problem, class Element, class Scv> | ||
| 138 | 1193800 | void completeFluidState(const ElemSol& elemSol, | |
| 139 | const Problem& problem, | ||
| 140 | const Element& element, | ||
| 141 | const Scv& scv, | ||
| 142 | FluidState& fluidState, | ||
| 143 | SolidState& solidState) | ||
| 144 | { | ||
| 145 | |||
| 146 | ///////////// | ||
| 147 | // set the fluid phase temperatures | ||
| 148 | ///////////// | ||
| 149 | 1193800 | EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState, solidState); | |
| 150 | |||
| 151 | ///////////// | ||
| 152 | // set the phase saturations | ||
| 153 | ///////////// | ||
| 154 | 1193800 | auto&& priVars = elemSol[scv.localDofIndex()]; | |
| 155 | 1193800 | Scalar sumSat = 0; | |
| 156 |
2/2✓ Branch 0 taken 1193800 times.
✓ Branch 1 taken 1193800 times.
|
2387600 | for (int phaseIdx = 0; phaseIdx < numFluidPhases() - 1; ++phaseIdx) { |
| 157 | 1193800 | sumSat += priVars[Indices::s0Idx + phaseIdx]; | |
| 158 | 3581400 | fluidState.setSaturation(phaseIdx, priVars[Indices::s0Idx + phaseIdx]); | |
| 159 | } | ||
| 160 | 1193800 | fluidState.setSaturation(numFluidPhases() - 1, 1.0 - sumSat); | |
| 161 | |||
| 162 | ///////////// | ||
| 163 | // set the phase pressures | ||
| 164 | ///////////// | ||
| 165 | // capillary pressure parameters | ||
| 166 | 1193800 | const auto& spatialParams = problem.spatialParams(); | |
| 167 | 1193800 | const int wPhaseIdx = spatialParams.template wettingPhase<FluidSystem>(element, scv, elemSol); | |
| 168 | 1193800 | fluidState.setWettingPhase(wPhaseIdx); | |
| 169 | // capillary pressures | ||
| 170 | |||
| 171 | 1193800 | const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol); | |
| 172 | 1193800 | const auto capPress = fluidMatrixInteraction.capillaryPressures(fluidState, wPhaseIdx); | |
| 173 | |||
| 174 | // add to the pressure of the first fluid phase | ||
| 175 | // depending on which pressure is stored in the primary variables | ||
| 176 | if(pressureFormulation == MpNcPressureFormulation::mostWettingFirst){ | ||
| 177 | // This means that the pressures are sorted from the most wetting to the least wetting-1 in the primary variables vector. | ||
| 178 | // For two phases this means that there is one pressure as primary variable: pw | ||
| 179 | 1193800 | const Scalar pw = priVars[Indices::p0Idx]; | |
| 180 |
2/2✓ Branch 0 taken 2387600 times.
✓ Branch 1 taken 1193800 times.
|
3581400 | for (int phaseIdx = 0; phaseIdx < numFluidPhases(); ++phaseIdx) |
| 181 | 9550400 | fluidState.setPressure(phaseIdx, pw - capPress[0] + capPress[phaseIdx]); | |
| 182 | } | ||
| 183 | else if(pressureFormulation == MpNcPressureFormulation::leastWettingFirst){ | ||
| 184 | // This means that the pressures are sorted from the least wetting to the most wetting-1 in the primary variables vector. | ||
| 185 | // For two phases this means that there is one pressure as primary variable: pn | ||
| 186 | const Scalar pn = priVars[Indices::p0Idx]; | ||
| 187 | for (int phaseIdx = numFluidPhases()-1; phaseIdx >= 0; --phaseIdx) | ||
| 188 | fluidState.setPressure(phaseIdx, pn - capPress[numFluidPhases()-1] + capPress[phaseIdx]); | ||
| 189 | } | ||
| 190 | else | ||
| 191 | DUNE_THROW(Dune::InvalidStateException, "MPNCVolumeVariables do not support the chosen pressure formulation"); | ||
| 192 | |||
| 193 | ///////////// | ||
| 194 | // set the fluid compositions | ||
| 195 | ///////////// | ||
| 196 | typename FluidSystem::ParameterCache paramCache; | ||
| 197 | 1193800 | paramCache.updateAll(fluidState); | |
| 198 | |||
| 199 | 1193800 | ComponentVector fug; | |
| 200 | // retrieve component fugacities | ||
| 201 |
2/2✓ Branch 0 taken 2387600 times.
✓ Branch 1 taken 1193800 times.
|
3581400 | for (int compIdx = 0; compIdx < numFluidComps; ++compIdx) |
| 202 | 7162800 | fug[compIdx] = priVars[Indices::fug0Idx + compIdx]; | |
| 203 | |||
| 204 | // calculate phase compositions | ||
| 205 |
2/2✓ Branch 0 taken 2387600 times.
✓ Branch 1 taken 1193800 times.
|
3581400 | for (int phaseIdx = 0; phaseIdx < numFluidPhases(); ++phaseIdx) { |
| 206 | // initial guess | ||
| 207 |
2/2✓ Branch 0 taken 4775200 times.
✓ Branch 1 taken 2387600 times.
|
7162800 | for (int compIdx = 0; compIdx < numFluidComps; ++compIdx) { |
| 208 | // sanity check to make sure we have non-zero fugacity coefficients | ||
| 209 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 4775200 times.
|
4775200 | if (FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx) == 0.0) |
| 210 | ✗ | DUNE_THROW(NumericalProblem, "MPNCVolumeVariables do not support fluidsystems with fugacity coefficients of 0"); | |
| 211 | |||
| 212 | // set initial guess of the component's mole fraction | ||
| 213 | 9550400 | fluidState.setMoleFraction(phaseIdx, compIdx, 1.0/numFluidComps); | |
| 214 | } | ||
| 215 | |||
| 216 | // calculate the phase composition from the component | ||
| 217 | // fugacities | ||
| 218 | 2387600 | CompositionFromFugacities::guessInitial(fluidState, paramCache, phaseIdx, fug); | |
| 219 | 2387600 | CompositionFromFugacities::solve(fluidState, paramCache, phaseIdx, fug); | |
| 220 | } | ||
| 221 | |||
| 222 | // dynamic viscosities | ||
| 223 |
2/2✓ Branch 0 taken 2387600 times.
✓ Branch 1 taken 1193800 times.
|
3581400 | for (int phaseIdx = 0; phaseIdx < numFluidPhases(); ++phaseIdx) { |
| 224 | // viscosities | ||
| 225 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
2387600 | Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx); |
| 226 | 2387600 | fluidState.setViscosity(phaseIdx, mu); | |
| 227 | |||
| 228 | // compute and set the enthalpy | ||
| 229 | 2387600 | Scalar h = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx); | |
| 230 | 4775200 | fluidState.setEnthalpy(phaseIdx, h); | |
| 231 | } | ||
| 232 | 1193800 | } | |
| 233 | |||
| 234 | /*! | ||
| 235 | * \brief Returns the fluid configuration at the given primary | ||
| 236 | * variables. | ||
| 237 | */ | ||
| 238 | const FluidState &fluidState() const | ||
| 239 |
0/8✗ 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.
|
4989696 | { return fluidState_; } |
| 240 | |||
| 241 | /*! | ||
| 242 | * \brief Returns the phase state for the control-volume. | ||
| 243 | */ | ||
| 244 | const SolidState &solidState() const | ||
| 245 | ✗ | { return solidState_; } | |
| 246 | |||
| 247 | /*! | ||
| 248 | * \brief Returns the saturation of a given phase within | ||
| 249 | * the control volume in \f$[-]\f$. | ||
| 250 | * | ||
| 251 | * \param phaseIdx The phase index | ||
| 252 | */ | ||
| 253 | Scalar saturation(int phaseIdx) const | ||
| 254 |
4/4✓ Branch 0 taken 443 times.
✓ Branch 1 taken 2387157 times.
✓ Branch 2 taken 443 times.
✓ Branch 3 taken 2387157 times.
|
4775200 | { return fluidState_.saturation(phaseIdx); } |
| 255 | |||
| 256 | /*! | ||
| 257 | * \brief Returns the mass fraction of a given component in a | ||
| 258 | * given phase within the control volume in \f$[-]\f$. | ||
| 259 | * | ||
| 260 | * \param phaseIdx The phase index | ||
| 261 | * \param compIdx The component index | ||
| 262 | */ | ||
| 263 | Scalar massFraction(const int phaseIdx, const int compIdx) const | ||
| 264 | 21265584 | { return fluidState_.massFraction(phaseIdx, compIdx); } | |
| 265 | |||
| 266 | /*! | ||
| 267 | * \brief Returns the mole fraction of a given component in a | ||
| 268 | * given phase within the control volume in \f$[-]\f$. | ||
| 269 | * | ||
| 270 | * \param phaseIdx The phase index | ||
| 271 | * \param compIdx The component index | ||
| 272 | */ | ||
| 273 | Scalar moleFraction(const int phaseIdx, const int compIdx) const | ||
| 274 | 57959808 | { return fluidState_.moleFraction(phaseIdx, compIdx); } | |
| 275 | |||
| 276 | /*! | ||
| 277 | * \brief Returns the concentration \f$\mathrm{[mol/m^3]}\f$ of a component in the phase. | ||
| 278 | * | ||
| 279 | * \param phaseIdx The phase index | ||
| 280 | * \param compIdx The index of the component | ||
| 281 | */ | ||
| 282 | Scalar molarity(const int phaseIdx, int compIdx) const | ||
| 283 | { return fluidState_.molarity(phaseIdx, compIdx); } | ||
| 284 | |||
| 285 | /*! | ||
| 286 | * \brief Returns the molar density \f$\mathrm{[mol/m^3]}\f$ the of the fluid phase. | ||
| 287 | * | ||
| 288 | * \param phaseIdx The phase index | ||
| 289 | */ | ||
| 290 | Scalar molarDensity(const int phaseIdx) const | ||
| 291 | 57959808 | { return fluidState_.molarDensity(phaseIdx);} | |
| 292 | |||
| 293 | /*! | ||
| 294 | * \brief Returns the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within | ||
| 295 | * the control volume. | ||
| 296 | * | ||
| 297 | * \param phaseIdx The phase index | ||
| 298 | */ | ||
| 299 | Scalar pressure(const int phaseIdx) const | ||
| 300 | 47571648 | { return fluidState_.pressure(phaseIdx); } | |
| 301 | |||
| 302 | /*! | ||
| 303 | * \brief Returns the density \f$\mathrm{[kg/m^3]}\f$ the of the fluid phase. | ||
| 304 | * | ||
| 305 | * \param phaseIdx The phase index | ||
| 306 | */ | ||
| 307 | Scalar density(const int phaseIdx) const | ||
| 308 | 95132576 | { return fluidState_.density(phaseIdx); } | |
| 309 | |||
| 310 | /*! | ||
| 311 | * \brief Returns the temperature inside the sub-control volume. | ||
| 312 | * | ||
| 313 | * Note that we assume thermodynamic equilibrium, i.e. the | ||
| 314 | * temperature of the rock matrix and of all fluid phases are | ||
| 315 | * identical. | ||
| 316 | */ | ||
| 317 | Scalar temperature() const | ||
| 318 | { return fluidState_.temperature(0/* phaseIdx*/); } | ||
| 319 | |||
| 320 | Scalar temperature(const int phaseIdx) const | ||
| 321 | { return fluidState_.temperature(phaseIdx); } | ||
| 322 | |||
| 323 | /*! | ||
| 324 | * \brief Return enthalpy \f$\mathrm{[kg/m^3]}\f$ the of the fluid phase. | ||
| 325 | */ | ||
| 326 | Scalar enthalpy(const int phaseIdx) const | ||
| 327 | ✗ | { return fluidState_.enthalpy(phaseIdx); } | |
| 328 | |||
| 329 | /*! | ||
| 330 | * \brief Returns the internal energy \f$\mathrm{[kg/m^3]}\f$ the of the fluid phase. | ||
| 331 | */ | ||
| 332 | Scalar internalEnergy(const int phaseIdx) const | ||
| 333 | ✗ | { return fluidState_.internalEnergy(phaseIdx); } | |
| 334 | |||
| 335 | /*! | ||
| 336 | * \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of a fluid phase in | ||
| 337 | * the sub-control volume. | ||
| 338 | */ | ||
| 339 | ✗ | Scalar fluidThermalConductivity(const int phaseIdx) const | |
| 340 | ✗ | { return FluidSystem::thermalConductivity(fluidState_, phaseIdx); } | |
| 341 | |||
| 342 | /*! | ||
| 343 | * \brief Returns the fugacity \f$\mathrm{[kg/m^3]}\f$ the of the component. | ||
| 344 | */ | ||
| 345 | Scalar fugacity(const int compIdx) const | ||
| 346 | 76800 | { return fluidState_.fugacity(compIdx); } | |
| 347 | |||
| 348 | /*! | ||
| 349 | * \brief Returns the average molar mass \f$\mathrm{[kg/m^3]}\f$ the of the phase. | ||
| 350 | */ | ||
| 351 | Scalar averageMolarMass(const int phaseIdx) const | ||
| 352 | { return fluidState_.averageMolarMass(phaseIdx); } | ||
| 353 | |||
| 354 | /*! | ||
| 355 | * \brief Returns the effective mobility of a given phase within | ||
| 356 | * the control volume. | ||
| 357 | * | ||
| 358 | * \param phaseIdx The local index of the phases | ||
| 359 | */ | ||
| 360 | Scalar mobility(const unsigned int phaseIdx) const | ||
| 361 | { | ||
| 362 | 28979904 | return relativePermeability(phaseIdx)/fluidState_.viscosity(phaseIdx); | |
| 363 | } | ||
| 364 | |||
| 365 | /*! | ||
| 366 | * \brief Returns the viscosity of a given phase within | ||
| 367 | * the control volume. | ||
| 368 | */ | ||
| 369 | Scalar viscosity(const unsigned int phaseIdx) const | ||
| 370 | { return fluidState_.viscosity(phaseIdx); } | ||
| 371 | |||
| 372 | /*! | ||
| 373 | * \brief Returns the relative permeability of a given phase within | ||
| 374 | * the control volume. | ||
| 375 | * | ||
| 376 | * \param phaseIdx The local index of the phases | ||
| 377 | */ | ||
| 378 | Scalar relativePermeability(const unsigned int phaseIdx) const | ||
| 379 | 57959808 | { return relativePermeability_[phaseIdx]; } | |
| 380 | |||
| 381 | /*! | ||
| 382 | * \brief Returns the average porosity within the control volume. | ||
| 383 | */ | ||
| 384 | Scalar porosity() const | ||
| 385 |
4/4✓ Branch 0 taken 443 times.
✓ Branch 1 taken 2387157 times.
✓ Branch 2 taken 443 times.
✓ Branch 3 taken 2387157 times.
|
4775200 | { return solidState_.porosity(); } |
| 386 | |||
| 387 | /*! | ||
| 388 | * \brief Returns the permeability within the control volume in \f$[m^2]\f$. | ||
| 389 | */ | ||
| 390 | const PermeabilityType& permeability() const | ||
| 391 | ✗ | { return permeability_; } | |
| 392 | |||
| 393 | /*! | ||
| 394 | * \brief Returns true if the fluid state is in the active set | ||
| 395 | * for a phase, | ||
| 396 | * | ||
| 397 | * \param phaseIdx The local index of the phases | ||
| 398 | */ | ||
| 399 | bool isPhaseActive(const unsigned int phaseIdx) const | ||
| 400 | { | ||
| 401 | return | ||
| 402 | phasePresentIneq(fluidState(), phaseIdx) - | ||
| 403 | phaseNotPresentIneq(fluidState(), phaseIdx) | ||
| 404 | >= 0; | ||
| 405 | } | ||
| 406 | |||
| 407 | /*! | ||
| 408 | * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$. | ||
| 409 | */ | ||
| 410 | Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const | ||
| 411 | { | ||
| 412 | typename FluidSystem::ParameterCache paramCache; | ||
| 413 | 2387600 | paramCache.updatePhase(fluidState_, phaseIdx); | |
| 414 |
2/2✓ Branch 1 taken 443 times.
✓ Branch 2 taken 2387157 times.
|
2387600 | return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx); |
| 415 | } | ||
| 416 | |||
| 417 | /*! | ||
| 418 | * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$. | ||
| 419 | */ | ||
| 420 | Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const | ||
| 421 | 9902496 | { return effectiveDiffCoeff_(phaseIdx, compIIdx, compJIdx); } | |
| 422 | |||
| 423 | /*! | ||
| 424 | * \brief Returns the value of the NCP-function for a phase. | ||
| 425 | * | ||
| 426 | * \param phaseIdx The local index of the phases | ||
| 427 | */ | ||
| 428 | Scalar phaseNcp(const unsigned int phaseIdx) const | ||
| 429 | { | ||
| 430 | 4989696 | Scalar aEval = phaseNotPresentIneq(fluidState(), phaseIdx); | |
| 431 | 14969088 | Scalar bEval = phasePresentIneq(fluidState(), phaseIdx); | |
| 432 | 4989696 | if (aEval > bEval) | |
| 433 | return phasePresentIneq(fluidState(), phaseIdx); | ||
| 434 | return phaseNotPresentIneq(fluidState(), phaseIdx); | ||
| 435 | }; | ||
| 436 | |||
| 437 | /*! | ||
| 438 | * \brief Returns the value of the inequality where a phase is | ||
| 439 | * present. | ||
| 440 | * | ||
| 441 | * \param phaseIdx The local index of the phases | ||
| 442 | * \param fluidState Container for all the secondary variables concerning the fluids | ||
| 443 | */ | ||
| 444 | ✗ | Scalar phasePresentIneq(const FluidState &fluidState, | |
| 445 | const unsigned int phaseIdx) const | ||
| 446 |
4/4✓ Branch 0 taken 3804253 times.
✓ Branch 1 taken 1185443 times.
✓ Branch 2 taken 3804253 times.
✓ Branch 3 taken 1185443 times.
|
9979392 | { return fluidState.saturation(phaseIdx); } |
| 447 | |||
| 448 | /*! | ||
| 449 | * \brief Returns the value of the inequality where a phase is not | ||
| 450 | * present. | ||
| 451 | * | ||
| 452 | * \param phaseIdx The local index of the phases | ||
| 453 | * \param fluidState Container for all the secondary variables concerning the fluids | ||
| 454 | */ | ||
| 455 | ✗ | Scalar phaseNotPresentIneq(const FluidState &fluidState, | |
| 456 | const unsigned int phaseIdx) const | ||
| 457 | { | ||
| 458 | // difference of sum of mole fractions in the phase from 100% | ||
| 459 | ✗ | Scalar a = 1; | |
| 460 |
4/6✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 9979392 times.
✓ Branch 3 taken 4989696 times.
✓ Branch 4 taken 7608506 times.
✓ Branch 5 taken 3804253 times.
|
26381847 | for (int compIdx = 0; compIdx < numFluidComps; ++compIdx) |
| 461 | 35175796 | a -= fluidState.moleFraction(phaseIdx, compIdx); | |
| 462 | ✗ | return a; | |
| 463 | } | ||
| 464 | |||
| 465 | protected: | ||
| 466 | Scalar porosity_; //!< Effective porosity within the control volume | ||
| 467 | std::array<Scalar, ModelTraits::numFluidPhases()> relativePermeability_; //!< Effective relative permeability within the control volume | ||
| 468 | PermeabilityType permeability_; | ||
| 469 | |||
| 470 | //! Mass fractions of each component within each phase | ||
| 471 | FluidState fluidState_; | ||
| 472 | SolidState solidState_; | ||
| 473 | |||
| 474 | // Effective diffusion coefficients for the phases | ||
| 475 | DiffusionCoefficients effectiveDiffCoeff_; | ||
| 476 | }; | ||
| 477 | |||
| 478 | template <class Traits> | ||
| 479 | class MPNCVolumeVariablesImplementation<Traits, true> | ||
| 480 | : public PorousMediumFlowVolumeVariables<Traits> | ||
| 481 | , public EnergyVolumeVariables<Traits, MPNCVolumeVariables<Traits> > | ||
| 482 | { | ||
| 483 | using ParentType = PorousMediumFlowVolumeVariables< Traits>; | ||
| 484 | using EnergyVolVars = EnergyVolumeVariables<Traits, MPNCVolumeVariables<Traits> >; | ||
| 485 | using Scalar = typename Traits::PrimaryVariables::value_type; | ||
| 486 | using PermeabilityType = typename Traits::PermeabilityType; | ||
| 487 | |||
| 488 | using ModelTraits = typename Traits::ModelTraits; | ||
| 489 | static constexpr auto pressureFormulation = ModelTraits::pressureFormulation(); | ||
| 490 | |||
| 491 | static constexpr bool enableThermalNonEquilibrium = ModelTraits::enableThermalNonEquilibrium(); | ||
| 492 | static constexpr bool enableChemicalNonEquilibrium = ModelTraits::enableChemicalNonEquilibrium(); | ||
| 493 | static constexpr bool enableDiffusion = ModelTraits::enableMolecularDiffusion(); | ||
| 494 | |||
| 495 | using Indices = typename ModelTraits::Indices; | ||
| 496 | using ComponentVector = Dune::FieldVector<Scalar, ModelTraits::numFluidComponents()>; | ||
| 497 | using CompositionFromFugacities = Dumux::CompositionFromFugacities<Scalar, typename Traits::FluidSystem>; | ||
| 498 | using ParameterCache = typename Traits::FluidSystem::ParameterCache; | ||
| 499 | using EffDiffModel = typename Traits::EffectiveDiffusivityModel; | ||
| 500 | using DiffusionCoefficients = typename Traits::DiffusionType::DiffusionCoefficientsContainer; | ||
| 501 | |||
| 502 | public: | ||
| 503 | //! Export the underlying fluid system | ||
| 504 | using FluidSystem = typename Traits::FluidSystem; | ||
| 505 | //! Export the fluid state type | ||
| 506 | using FluidState = typename Traits::FluidState; | ||
| 507 | //! Export type of solid state | ||
| 508 | using SolidState = typename Traits::SolidState; | ||
| 509 | //! Export type of solid system | ||
| 510 | using SolidSystem = typename Traits::SolidSystem; | ||
| 511 | |||
| 512 | //! Return number of phases considered by the model | ||
| 513 | static constexpr int numFluidPhases() { return ModelTraits::numFluidPhases(); } | ||
| 514 | //! Return number of components considered by the model | ||
| 515 | static constexpr int numFluidComps = ParentType::numFluidComponents(); | ||
| 516 | using ConstraintSolver = MiscibleMultiPhaseComposition<Scalar, FluidSystem>; | ||
| 517 | |||
| 518 | /*! | ||
| 519 | * \brief Updates all quantities for a given control volume. | ||
| 520 | * | ||
| 521 | * \param elemSol A vector containing all primary variables connected to the element | ||
| 522 | * \param problem The object specifying the problem which ought to | ||
| 523 | * be simulated | ||
| 524 | * \param element An element which contains part of the control volume | ||
| 525 | * \param scv The sub-control volume | ||
| 526 | */ | ||
| 527 | template<class ElemSol, class Problem, class Element, class Scv> | ||
| 528 | 806400 | void update(const ElemSol& elemSol, | |
| 529 | const Problem& problem, | ||
| 530 | const Element& element, | ||
| 531 | const Scv& scv) | ||
| 532 | { | ||
| 533 | 806400 | ParentType::update(elemSol, problem, element, scv); | |
| 534 | |||
| 535 | 806400 | completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_); | |
| 536 | |||
| 537 | // calculate the remaining quantities | ||
| 538 | 806400 | const auto& spatialParams = problem.spatialParams(); | |
| 539 | 806400 | const int wPhaseIdx = spatialParams.template wettingPhase<FluidSystem>(element, scv, elemSol); | |
| 540 | 806400 | const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol); | |
| 541 | 806400 | const auto relPerm = fluidMatrixInteraction.relativePermeabilities(fluidState_, wPhaseIdx); | |
| 542 | 3225600 | std::copy(relPerm.begin(), relPerm.end(), relativePermeability_.begin()); | |
| 543 | |||
| 544 | typename FluidSystem::ParameterCache paramCache; | ||
| 545 | 806400 | paramCache.updateAll(fluidState_); | |
| 546 | |||
| 547 | 806400 | updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps); | |
| 548 | |||
| 549 | if constexpr (enableDiffusion) | ||
| 550 | { | ||
| 551 | 806400 | auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx) | |
| 552 | 1612800 | { return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx); }; | |
| 553 | |||
| 554 | 806400 | effectiveDiffCoeff_.update(getEffectiveDiffusionCoefficient); | |
| 555 | } | ||
| 556 | |||
| 557 | 806400 | EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_); | |
| 558 | 806400 | permeability_ = spatialParams.permeability(element, scv, elemSol); | |
| 559 | 806400 | EnergyVolVars::updateEffectiveThermalConductivity(); | |
| 560 | 806400 | } | |
| 561 | |||
| 562 | /*! | ||
| 563 | * \brief Sets complete fluid state. | ||
| 564 | * | ||
| 565 | * \param elemSol A vector containing all primary variables connected to the element | ||
| 566 | * \param problem The object specifying the problem which ought to | ||
| 567 | * be simulated | ||
| 568 | * \param element An element which contains part of the control volume | ||
| 569 | * \param scv The sub-control volume | ||
| 570 | * \param fluidState A container with the current (physical) state of the fluid | ||
| 571 | * \param solidState A container with the current (physical) state of the solid | ||
| 572 | */ | ||
| 573 | template<class ElemSol, class Problem, class Element, class Scv> | ||
| 574 | 806400 | void completeFluidState(const ElemSol& elemSol, | |
| 575 | const Problem& problem, | ||
| 576 | const Element& element, | ||
| 577 | const Scv& scv, | ||
| 578 | FluidState& fluidState, | ||
| 579 | SolidState& solidState) | ||
| 580 | { | ||
| 581 | ///////////// | ||
| 582 | // set the fluid phase temperatures | ||
| 583 | ///////////// | ||
| 584 | 806400 | EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState, solidState); | |
| 585 | ///////////// | ||
| 586 | // set the phase saturations | ||
| 587 | ///////////// | ||
| 588 | 806400 | auto&& priVars = elemSol[scv.localDofIndex()]; | |
| 589 | 806400 | Scalar sumSat = 0; | |
| 590 |
2/2✓ Branch 0 taken 806400 times.
✓ Branch 1 taken 806400 times.
|
1612800 | for (int phaseIdx = 0; phaseIdx < numFluidPhases() - 1; ++phaseIdx) { |
| 591 | 806400 | sumSat += priVars[Indices::s0Idx + phaseIdx]; | |
| 592 | 2419200 | fluidState.setSaturation(phaseIdx, priVars[Indices::s0Idx + phaseIdx]); | |
| 593 | } | ||
| 594 | 806400 | fluidState.setSaturation(numFluidPhases() - 1, 1.0 - sumSat); | |
| 595 | |||
| 596 | ///////////// | ||
| 597 | // set the phase pressures | ||
| 598 | ///////////// | ||
| 599 | // capillary pressure parameters | ||
| 600 | 806400 | const auto& spatialParams = problem.spatialParams(); | |
| 601 | 806400 | const int wPhaseIdx = spatialParams.template wettingPhase<FluidSystem>(element, scv, elemSol); | |
| 602 | 806400 | fluidState.setWettingPhase(wPhaseIdx); | |
| 603 | 806400 | const auto fluidMatrixInteraction = spatialParams.fluidMatrixInteraction(element, scv, elemSol); | |
| 604 | 806400 | const auto capPress = fluidMatrixInteraction.capillaryPressures(fluidState, wPhaseIdx); | |
| 605 | |||
| 606 | // add to the pressure of the first fluid phase | ||
| 607 | // depending on which pressure is stored in the primary variables | ||
| 608 | if(pressureFormulation == MpNcPressureFormulation::mostWettingFirst){ | ||
| 609 | // This means that the pressures are sorted from the most wetting to the least wetting-1 in the primary variables vector. | ||
| 610 | // For two phases this means that there is one pressure as primary variable: pw | ||
| 611 | const Scalar pw = priVars[Indices::p0Idx]; | ||
| 612 | for (int phaseIdx = 0; phaseIdx < numFluidPhases(); ++phaseIdx) | ||
| 613 | fluidState.setPressure(phaseIdx, pw - capPress[0] + capPress[phaseIdx]); | ||
| 614 | } | ||
| 615 | else if(pressureFormulation == MpNcPressureFormulation::leastWettingFirst){ | ||
| 616 | // This means that the pressures are sorted from the least wetting to the most wetting-1 in the primary variables vector. | ||
| 617 | // For two phases this means that there is one pressure as primary variable: pn | ||
| 618 | 806400 | const Scalar pn = priVars[Indices::p0Idx]; | |
| 619 |
2/2✓ Branch 0 taken 1612800 times.
✓ Branch 1 taken 806400 times.
|
2419200 | for (int phaseIdx = numFluidPhases()-1; phaseIdx >= 0; --phaseIdx) |
| 620 | 6451200 | fluidState.setPressure(phaseIdx, pn - capPress[numFluidPhases()-1] + capPress[phaseIdx]); | |
| 621 | } | ||
| 622 | else | ||
| 623 | DUNE_THROW(Dune::InvalidStateException, "MPNCVolumeVariables do not support the chosen pressure formulation"); | ||
| 624 | |||
| 625 | |||
| 626 | ///////////// | ||
| 627 | // set the fluid compositions | ||
| 628 | ///////////// | ||
| 629 | typename FluidSystem::ParameterCache paramCache; | ||
| 630 | 806400 | paramCache.updateAll(fluidState); | |
| 631 | |||
| 632 | 806400 | updateMoleFraction(fluidState, paramCache, priVars); | |
| 633 | |||
| 634 | // dynamic viscosities | ||
| 635 |
2/2✓ Branch 0 taken 1612800 times.
✓ Branch 1 taken 806400 times.
|
2419200 | for (int phaseIdx = 0; phaseIdx < numFluidPhases(); ++phaseIdx) { |
| 636 | // viscosities | ||
| 637 | 1612800 | Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx); | |
| 638 | 1612800 | fluidState.setViscosity(phaseIdx, mu); | |
| 639 | |||
| 640 | // compute and set the enthalpy | ||
| 641 | 1612800 | Scalar h = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx); | |
| 642 | 3225600 | fluidState.setEnthalpy(phaseIdx, h); | |
| 643 | } | ||
| 644 | 806400 | } | |
| 645 | |||
| 646 | /*! | ||
| 647 | * \brief Updates composition of all phases in the mutable | ||
| 648 | * parameters from the primary variables. | ||
| 649 | * | ||
| 650 | * \param actualFluidState Container for all the secondary variables concerning the fluids | ||
| 651 | * \param paramCache Container for cache parameters | ||
| 652 | * \param priVars The primary Variables | ||
| 653 | */ | ||
| 654 | 806400 | void updateMoleFraction(FluidState & actualFluidState, | |
| 655 | ParameterCache & paramCache, | ||
| 656 | const typename Traits::PrimaryVariables& priVars) | ||
| 657 | { | ||
| 658 | // setting the mole fractions of the fluid state | ||
| 659 |
2/2✓ Branch 0 taken 1612800 times.
✓ Branch 1 taken 806400 times.
|
2419200 | for(int phaseIdx=0; phaseIdx<numFluidPhases(); ++phaseIdx) |
| 660 | { | ||
| 661 | // set the component mole fractions | ||
| 662 |
2/2✓ Branch 0 taken 3225600 times.
✓ Branch 1 taken 1612800 times.
|
4838400 | for (int compIdx = 0; compIdx < numFluidComps; ++compIdx) { |
| 663 | 3225600 | actualFluidState.setMoleFraction(phaseIdx, | |
| 664 | compIdx, | ||
| 665 | 3225600 | priVars[Indices::moleFrac00Idx + | |
| 666 | 3225600 | phaseIdx*numFluidComps + | |
| 667 | compIdx]); | ||
| 668 | } | ||
| 669 | } | ||
| 670 | |||
| 671 | // TODO: Can this be removed, are fugacity coefficients used in the constraint solver? | ||
| 672 |
2/2✓ Branch 0 taken 1612800 times.
✓ Branch 1 taken 806400 times.
|
2419200 | for(int phaseIdx=0; phaseIdx<numFluidPhases(); ++phaseIdx) |
| 673 |
2/2✓ Branch 0 taken 3225600 times.
✓ Branch 1 taken 1612800 times.
|
4838400 | for (int compIdx = 0; compIdx < numFluidComps; ++compIdx) { |
| 674 | 3225600 | const Scalar phi = FluidSystem::fugacityCoefficient(actualFluidState, | |
| 675 | paramCache, | ||
| 676 | phaseIdx, | ||
| 677 | compIdx); | ||
| 678 | 6451200 | actualFluidState.setFugacityCoefficient(phaseIdx, | |
| 679 | compIdx, | ||
| 680 | phi); | ||
| 681 | } | ||
| 682 | |||
| 683 | 806400 | FluidState equilFluidState; // the fluidState *on the interface* i.e. chemical equilibrium | |
| 684 | 806400 | equilFluidState.assign(actualFluidState) ; | |
| 685 | 806400 | ConstraintSolver::solve(equilFluidState, | |
| 686 | paramCache) ; | ||
| 687 | |||
| 688 | // Setting the equilibrium composition (in a kinetic model not necessarily the same as the actual mole fraction) | ||
| 689 |
2/2✓ Branch 0 taken 1612800 times.
✓ Branch 1 taken 806400 times.
|
2419200 | for(int phaseIdx=0; phaseIdx<numFluidPhases(); ++phaseIdx){ |
| 690 |
2/2✓ Branch 0 taken 3225600 times.
✓ Branch 1 taken 1612800 times.
|
4838400 | for (int compIdx=0; compIdx< numFluidComps; ++ compIdx){ |
| 691 | 12902400 | xEquil_[phaseIdx][compIdx] = equilFluidState.moleFraction(phaseIdx, compIdx); | |
| 692 | } | ||
| 693 | } | ||
| 694 | |||
| 695 | // compute densities of all phases | ||
| 696 |
2/2✓ Branch 0 taken 1612800 times.
✓ Branch 1 taken 806400 times.
|
2419200 | for(int phaseIdx=0; phaseIdx<numFluidPhases(); ++phaseIdx){ |
| 697 | 1612800 | const Scalar rho = FluidSystem::density(actualFluidState, paramCache, phaseIdx); | |
| 698 | 1612800 | actualFluidState.setDensity(phaseIdx, rho); | |
| 699 | 1612800 | const Scalar rhoMolar = FluidSystem::molarDensity(actualFluidState, paramCache, phaseIdx); | |
| 700 | 3225600 | actualFluidState.setMolarDensity(phaseIdx, rhoMolar); | |
| 701 | } | ||
| 702 | |||
| 703 | 806400 | } | |
| 704 | |||
| 705 | /*! | ||
| 706 | * \brief The mole fraction we would have in the case of chemical equilibrium / | ||
| 707 | * on the interface. | ||
| 708 | * | ||
| 709 | * \param phaseIdx The index of the fluid phase | ||
| 710 | * \param compIdx The local index of the component | ||
| 711 | */ | ||
| 712 | const Scalar xEquil(const unsigned int phaseIdx, const unsigned int compIdx) const | ||
| 713 | { | ||
| 714 | 14545440 | return xEquil_[phaseIdx][compIdx] ; | |
| 715 | } | ||
| 716 | |||
| 717 | /*! | ||
| 718 | * \brief Returns the fluid configuration at the given primary | ||
| 719 | * variables | ||
| 720 | */ | ||
| 721 | const FluidState &fluidState() const | ||
| 722 | 26661600 | { return fluidState_; } | |
| 723 | |||
| 724 | /*! | ||
| 725 | * \brief Returns the phase state for the control-volume. | ||
| 726 | */ | ||
| 727 | const SolidState &solidState() const | ||
| 728 |
1/2✓ Branch 1 taken 2424240 times.
✗ Branch 2 not taken.
|
14545440 | { return solidState_; } |
| 729 | |||
| 730 | /*! | ||
| 731 | * \brief Returns the saturation of a given phase within | ||
| 732 | * the control volume in \f$[-]\f$. | ||
| 733 | * | ||
| 734 | * \param phaseIdx The phase index | ||
| 735 | */ | ||
| 736 | Scalar saturation(int phaseIdx) const | ||
| 737 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 1612800 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1612800 times.
|
64632960 | { return fluidState_.saturation(phaseIdx); } |
| 738 | |||
| 739 | /*! | ||
| 740 | * \brief Returns the mass fraction of a given component in a | ||
| 741 | * given phase within the control volume in \f$[-]\f$. | ||
| 742 | * | ||
| 743 | * \param phaseIdx The phase index | ||
| 744 | * \param compIdx The component index | ||
| 745 | */ | ||
| 746 | Scalar massFraction(const int phaseIdx, const int compIdx) const | ||
| 747 | 38787840 | { return fluidState_.massFraction(phaseIdx, compIdx); } | |
| 748 | |||
| 749 | /*! | ||
| 750 | * \brief Returns the mole fraction of a given component in a | ||
| 751 | * given phase within the control volume in \f$[-]\f$. | ||
| 752 | * | ||
| 753 | * \param phaseIdx The phase index | ||
| 754 | * \param compIdx The component index | ||
| 755 | */ | ||
| 756 | Scalar moleFraction(const int phaseIdx, const int compIdx) const | ||
| 757 | 87272640 | { return fluidState_.moleFraction(phaseIdx, compIdx); } | |
| 758 | |||
| 759 | /*! | ||
| 760 | * \brief Returns the concentration \f$\mathrm{[mol/m^3]}\f$ of a component in the phase. | ||
| 761 | * | ||
| 762 | * \param phaseIdx The phase index | ||
| 763 | * \param compIdx The index of the component | ||
| 764 | */ | ||
| 765 | Scalar molarity(const int phaseIdx, int compIdx) const | ||
| 766 | { return fluidState_.molarity(phaseIdx, compIdx); } | ||
| 767 | |||
| 768 | /*! | ||
| 769 | * \brief Returns the molar density \f$\mathrm{[mol/m^3]}\f$ the of the fluid phase. | ||
| 770 | * | ||
| 771 | * \param phaseIdx The phase index | ||
| 772 | */ | ||
| 773 | Scalar molarDensity(const int phaseIdx) const | ||
| 774 | 87272640 | { return fluidState_.molarDensity(phaseIdx);} | |
| 775 | |||
| 776 | /*! | ||
| 777 | * \brief Returns the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within | ||
| 778 | * the control volume. | ||
| 779 | * | ||
| 780 | * \param phaseIdx The phase index | ||
| 781 | */ | ||
| 782 | Scalar pressure(const int phaseIdx) const | ||
| 783 | 40158720 | { return fluidState_.pressure(phaseIdx); } | |
| 784 | |||
| 785 | /*! | ||
| 786 | * \brief Returns the density \f$\mathrm{[kg/m^3]}\f$ the of the fluid phase. | ||
| 787 | * | ||
| 788 | * \param phaseIdx The phase index | ||
| 789 | */ | ||
| 790 | Scalar density(const int phaseIdx) const | ||
| 791 | 156522240 | { return fluidState_.density(phaseIdx); } | |
| 792 | |||
| 793 | /*! | ||
| 794 | * \brief Returns the temperature inside the sub-control volume. | ||
| 795 | * | ||
| 796 | * Note that we assume thermodynamic equilibrium, i.e. the | ||
| 797 | * temperature of the rock matrix and of all fluid phases are | ||
| 798 | * identical. | ||
| 799 | */ | ||
| 800 | Scalar temperature() const | ||
| 801 | { return fluidState_.temperature(0/* phaseIdx*/); } | ||
| 802 | |||
| 803 | /*! | ||
| 804 | * \brief Returns the enthalpy \f$\mathrm{[kg/m^3]}\f$ the of the fluid phase. | ||
| 805 | */ | ||
| 806 | Scalar enthalpy(const int phaseIdx) const | ||
| 807 | 29090880 | { return fluidState_.enthalpy(phaseIdx); } | |
| 808 | |||
| 809 | /*! | ||
| 810 | * \brief Returns the internal energy \f$\mathrm{[kg/m^3]}\f$ the of the fluid phase. | ||
| 811 | */ | ||
| 812 | Scalar internalEnergy(const int phaseIdx) const | ||
| 813 | 19393920 | { return fluidState_.internalEnergy(phaseIdx); } | |
| 814 | |||
| 815 | /*! | ||
| 816 | * \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of a fluid phase in | ||
| 817 | * the sub-control volume. | ||
| 818 | */ | ||
| 819 | Scalar fluidThermalConductivity(const int phaseIdx) const | ||
| 820 | 2424240 | { return FluidSystem::thermalConductivity(fluidState_, phaseIdx); } | |
| 821 | |||
| 822 | /*! | ||
| 823 | * \brief Returns the fugacity \f$\mathrm{[kg/m^3]}\f$ the of the component. | ||
| 824 | */ | ||
| 825 | Scalar fugacity(const int compIdx) const | ||
| 826 | 80640 | { return fluidState_.fugacity(compIdx); } | |
| 827 | |||
| 828 | /*! | ||
| 829 | * \brief Returns the average molar mass \f$\mathrm{[kg/m^3]}\f$ the of the phase. | ||
| 830 | */ | ||
| 831 | Scalar averageMolarMass(const int phaseIdx) const | ||
| 832 | { return fluidState_.averageMolarMass(phaseIdx); } | ||
| 833 | |||
| 834 | /*! | ||
| 835 | * \brief Returns the effective mobility of a given phase within | ||
| 836 | * the control volume. | ||
| 837 | * | ||
| 838 | * \param phaseIdx The local index of the phases | ||
| 839 | */ | ||
| 840 | Scalar mobility(const unsigned int phaseIdx) const | ||
| 841 | { | ||
| 842 | 29433600 | return relativePermeability(phaseIdx)/fluidState_.viscosity(phaseIdx); | |
| 843 | } | ||
| 844 | |||
| 845 | /*! | ||
| 846 | * \brief Returns the viscosity of a given phase within | ||
| 847 | * the control volume. | ||
| 848 | */ | ||
| 849 | Scalar viscosity(const unsigned int phaseIdx) const | ||
| 850 | { return fluidState_.viscosity(phaseIdx); } | ||
| 851 | |||
| 852 | /*! | ||
| 853 | * \brief Returns the relative permeability of a given phase within | ||
| 854 | * the control volume. | ||
| 855 | * | ||
| 856 | * \param phaseIdx The local index of the phases | ||
| 857 | */ | ||
| 858 | Scalar relativePermeability(const unsigned int phaseIdx) const | ||
| 859 | 58867200 | { return relativePermeability_[phaseIdx]; } | |
| 860 | |||
| 861 | /*! | ||
| 862 | * \brief Returns the average porosity within the control volume. | ||
| 863 | */ | ||
| 864 | Scalar porosity() const | ||
| 865 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 1612800 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1612800 times.
|
67858560 | { return solidState_.porosity(); } |
| 866 | |||
| 867 | /*! | ||
| 868 | * \brief Returns the permeability within the control volume in \f$[m^2]\f$. | ||
| 869 | */ | ||
| 870 | const PermeabilityType& permeability() const | ||
| 871 | { return permeability_; } | ||
| 872 | |||
| 873 | /*! | ||
| 874 | * \brief Returns true if the fluid state is in the active set | ||
| 875 | * for a phase, | ||
| 876 | * | ||
| 877 | * \param phaseIdx The local index of the phases | ||
| 878 | */ | ||
| 879 | bool isPhaseActive(const unsigned int phaseIdx) const | ||
| 880 | { | ||
| 881 | return | ||
| 882 | phasePresentIneq(fluidState(), phaseIdx) - | ||
| 883 | phaseNotPresentIneq(fluidState(), phaseIdx) | ||
| 884 | >= 0; | ||
| 885 | } | ||
| 886 | |||
| 887 | /*! | ||
| 888 | * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$. | ||
| 889 | */ | ||
| 890 | Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const | ||
| 891 | { | ||
| 892 | typename FluidSystem::ParameterCache paramCache; | ||
| 893 | 6461280 | paramCache.updatePhase(fluidState_, phaseIdx); | |
| 894 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 1612800 times.
|
6461280 | return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx); |
| 895 | } | ||
| 896 | |||
| 897 | /*! | ||
| 898 | * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$. | ||
| 899 | */ | ||
| 900 | Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const | ||
| 901 | 9696960 | { return effectiveDiffCoeff_(phaseIdx, compIIdx, compJIdx); } | |
| 902 | |||
| 903 | /*! | ||
| 904 | * \brief Returns the value of the NCP-function for a phase. | ||
| 905 | * | ||
| 906 | * \param phaseIdx The local index of the phases | ||
| 907 | */ | ||
| 908 | Scalar phaseNcp(const unsigned int phaseIdx) const | ||
| 909 | { | ||
| 910 | 4848480 | Scalar aEval = phaseNotPresentIneq(fluidState(), phaseIdx); | |
| 911 | 14545440 | Scalar bEval = phasePresentIneq(fluidState(), phaseIdx); | |
| 912 | 4848480 | if (aEval > bEval) | |
| 913 | return phasePresentIneq(fluidState(), phaseIdx); | ||
| 914 | return phaseNotPresentIneq(fluidState(), phaseIdx); | ||
| 915 | }; | ||
| 916 | |||
| 917 | /*! | ||
| 918 | * \brief Returns the value of the inequality where a phase is | ||
| 919 | * present. | ||
| 920 | * | ||
| 921 | * \param phaseIdx The local index of the phases | ||
| 922 | * \param fluidState Container for all the secondary variables concerning the fluids | ||
| 923 | */ | ||
| 924 | ✗ | Scalar phasePresentIneq(const FluidState &fluidState, | |
| 925 | const unsigned int phaseIdx) const | ||
| 926 |
2/4✓ Branch 0 taken 4848480 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4848480 times.
✗ Branch 3 not taken.
|
9696960 | { return fluidState.saturation(phaseIdx); } |
| 927 | |||
| 928 | /*! | ||
| 929 | * \brief Returns the value of the inequality where a phase is not | ||
| 930 | * present. | ||
| 931 | * | ||
| 932 | * \param phaseIdx The local index of the phases | ||
| 933 | * \param fluidState Container for all the secondary variables concerning the fluids | ||
| 934 | */ | ||
| 935 | ✗ | Scalar phaseNotPresentIneq(const FluidState &fluidState, | |
| 936 | const unsigned int phaseIdx) const | ||
| 937 | { | ||
| 938 | // difference of sum of mole fractions in the phase from 100% | ||
| 939 | ✗ | Scalar a = 1; | |
| 940 |
4/6✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 9696960 times.
✓ Branch 3 taken 4848480 times.
✓ Branch 4 taken 9696960 times.
✓ Branch 5 taken 4848480 times.
|
29090880 | for (int compIdx = 0; compIdx < numFluidComps; ++compIdx) |
| 941 | 38787840 | a -= fluidState.moleFraction(phaseIdx, compIdx); | |
| 942 | ✗ | return a; | |
| 943 | } | ||
| 944 | |||
| 945 | protected: | ||
| 946 | std::array<Scalar, ModelTraits::numFluidPhases()> relativePermeability_; //!< Effective relative permeability within the control volume | ||
| 947 | PermeabilityType permeability_; | ||
| 948 | std::array<std::array<Scalar, numFluidComps>, numFluidPhases()> xEquil_; | ||
| 949 | |||
| 950 | //! Mass fractions of each component within each phase | ||
| 951 | FluidState fluidState_; | ||
| 952 | SolidState solidState_; | ||
| 953 | |||
| 954 | // Effective diffusion coefficients for the phases | ||
| 955 | DiffusionCoefficients effectiveDiffCoeff_; | ||
| 956 | }; | ||
| 957 | |||
| 958 | } // end namespace | ||
| 959 | |||
| 960 | #endif | ||
| 961 |