| 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-FileCopyrightText: 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 FluidSystems | ||
| 10 | * \copybrief Dumux::FluidSystems::H2OAirXylene | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_H2O_AIR_XYLENE_FLUID_SYSTEM_HH | ||
| 13 | #define DUMUX_H2O_AIR_XYLENE_FLUID_SYSTEM_HH | ||
| 14 | |||
| 15 | #include <dumux/material/idealgas.hh> | ||
| 16 | #include <dumux/material/components/air.hh> | ||
| 17 | #include <dumux/material/components/h2o.hh> | ||
| 18 | #include <dumux/material/components/xylene.hh> | ||
| 19 | #include <dumux/material/components/tabulatedcomponent.hh> | ||
| 20 | |||
| 21 | #include <dumux/material/binarycoefficients/h2o_air.hh> | ||
| 22 | #include <dumux/material/binarycoefficients/h2o_xylene.hh> | ||
| 23 | #include <dumux/material/binarycoefficients/air_xylene.hh> | ||
| 24 | |||
| 25 | #include <dumux/material/fluidsystems/base.hh> | ||
| 26 | |||
| 27 | #include <dumux/io/name.hh> | ||
| 28 | |||
| 29 | namespace Dumux::FluidSystems { | ||
| 30 | |||
| 31 | /*! | ||
| 32 | * \ingroup FluidSystems | ||
| 33 | * \brief A three-phase fluid system featuring gas, NAPL and water as phases and | ||
| 34 | * distilled water \f$(\mathrm{H_2O})\f$ and air (Pseudo component composed of | ||
| 35 | * \f$\mathrm{79\%\;N_2}\f$, \f$\mathrm{20\%\;O_2}\f$ and Mesitylene \f$(\mathrm{C_8H_{10}})\f$ as components. | ||
| 36 | * | ||
| 37 | * \note This fluid system assumes all phases to be ideal mixtures. | ||
| 38 | */ | ||
| 39 | template <class Scalar, | ||
| 40 | class H2OType = Components::TabulatedComponent<Components::H2O<Scalar> > > | ||
| 41 | class H2OAirXylene | ||
| 42 | : public Base<Scalar, H2OAirXylene<Scalar, H2OType> > | ||
| 43 | { | ||
| 44 | using ThisType = H2OAirXylene<Scalar, H2OType>; | ||
| 45 | |||
| 46 | public: | ||
| 47 | using H2O = H2OType; | ||
| 48 | using NAPL = Components::Xylene<Scalar>; | ||
| 49 | using Air = Dumux::Components::Air<Scalar>; | ||
| 50 | |||
| 51 | static const int numPhases = 3; | ||
| 52 | static const int numComponents = 3; | ||
| 53 | |||
| 54 | static const int wPhaseIdx = 0; // index of the water phase | ||
| 55 | static const int nPhaseIdx = 1; // index of the NAPL phase | ||
| 56 | static const int gPhaseIdx = 2; // index of the gas phase | ||
| 57 | |||
| 58 | static const int H2OIdx = 0; | ||
| 59 | static const int NAPLIdx = 1; | ||
| 60 | static const int AirIdx = 2; | ||
| 61 | |||
| 62 | // export component indices to indicate the main component | ||
| 63 | // of the corresponding phase at atmospheric pressure 1 bar | ||
| 64 | // and room temperature 20°C: | ||
| 65 | static const int wCompIdx = H2OIdx; | ||
| 66 | static const int nCompIdx = NAPLIdx; | ||
| 67 | static const int gCompIdx = AirIdx; | ||
| 68 | |||
| 69 | /*! | ||
| 70 | * \brief Initialize the fluid system's static parameters generically | ||
| 71 | * | ||
| 72 | * If a tabulated H2O component is used, we do our best to create | ||
| 73 | * tables that always work. | ||
| 74 | */ | ||
| 75 | 3 | static void init() | |
| 76 | { | ||
| 77 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | init(/*tempMin=*/273.15, |
| 78 | /*tempMax=*/623.15, | ||
| 79 | /*numTemp=*/100, | ||
| 80 | /*pMin=*/0.0, | ||
| 81 | /*pMax=*/20e6, | ||
| 82 | /*numP=*/200); | ||
| 83 | 2 | } | |
| 84 | |||
| 85 | /*! | ||
| 86 | * \brief Initialize the fluid system's static parameters using | ||
| 87 | * problem specific temperature and pressure ranges | ||
| 88 | * | ||
| 89 | * \param tempMin The minimum temperature used for tabulation of water \f$\mathrm{[K]}\f$ | ||
| 90 | * \param tempMax The maximum temperature used for tabulation of water \f$\mathrm{[K]}\f$ | ||
| 91 | * \param nTemp The number of ticks on the temperature axis of the table of water | ||
| 92 | * \param pressMin The minimum pressure used for tabulation of water \f$\mathrm{[Pa]}\f$ | ||
| 93 | * \param pressMax The maximum pressure used for tabulation of water \f$\mathrm{[Pa]}\f$ | ||
| 94 | * \param nPress The number of ticks on the pressure axis of the table of water | ||
| 95 | */ | ||
| 96 | 3 | static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp, | |
| 97 | Scalar pressMin, Scalar pressMax, unsigned nPress) | ||
| 98 | { | ||
| 99 | if (H2O::isTabulated) | ||
| 100 | { | ||
| 101 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
3 | H2O::init(tempMin, tempMax, nTemp, |
| 102 | pressMin, pressMax, nPress); | ||
| 103 | } | ||
| 104 | } | ||
| 105 | |||
| 106 | /*! | ||
| 107 | * \brief Returns whether the fluids are miscible | ||
| 108 | */ | ||
| 109 | static constexpr bool isMiscible() | ||
| 110 | { return true; } | ||
| 111 | |||
| 112 | /*! | ||
| 113 | * \brief Return whether a phase is gaseous | ||
| 114 | * | ||
| 115 | * \param phaseIdx The index of the fluid phase to consider | ||
| 116 | */ | ||
| 117 | 3 | static constexpr bool isGas(int phaseIdx) | |
| 118 | { | ||
| 119 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
| 120 | return phaseIdx == gPhaseIdx; | ||
| 121 | } | ||
| 122 | |||
| 123 | /*! | ||
| 124 | * \brief Returns true if and only if a fluid phase is assumed to | ||
| 125 | * be an ideal gas. | ||
| 126 | * | ||
| 127 | * \param phaseIdx The index of the fluid phase to consider | ||
| 128 | */ | ||
| 129 | static bool isIdealGas(int phaseIdx) | ||
| 130 | { return phaseIdx == gPhaseIdx && H2O::gasIsIdeal() && Air::gasIsIdeal() && NAPL::gasIsIdeal(); } | ||
| 131 | |||
| 132 | /*! | ||
| 133 | * \brief Returns true if and only if a fluid phase is assumed to | ||
| 134 | * be an ideal mixture. | ||
| 135 | * | ||
| 136 | * We define an ideal mixture as a fluid phase where the fugacity | ||
| 137 | * coefficients of all components times the pressure of the phase | ||
| 138 | * are independent on the fluid composition. This assumption is true | ||
| 139 | * if Henry's law and Raoult's law apply. If you are unsure what | ||
| 140 | * this function should return, it is safe to return false. The | ||
| 141 | * only damage done will be (slightly) increased computation times | ||
| 142 | * in some cases. | ||
| 143 | * | ||
| 144 | * \param phaseIdx The index of the fluid phase to consider | ||
| 145 | */ | ||
| 146 | ✗ | static bool isIdealMixture(int phaseIdx) | |
| 147 | { | ||
| 148 | ✗ | assert(0 <= phaseIdx && phaseIdx < numPhases); | |
| 149 | // we assume Henry's and Raoult's laws for the water phase and | ||
| 150 | // and no interaction between gas molecules of different | ||
| 151 | // components, so all phases are ideal mixtures! | ||
| 152 | return true; | ||
| 153 | } | ||
| 154 | |||
| 155 | /*! | ||
| 156 | * \brief Returns true if and only if a fluid phase is assumed to | ||
| 157 | * be compressible. | ||
| 158 | * | ||
| 159 | * Compressible means that the partial derivative of the density | ||
| 160 | * to the fluid pressure is always larger than zero. | ||
| 161 | * | ||
| 162 | * \param phaseIdx The index of the fluid phase to consider | ||
| 163 | */ | ||
| 164 | 3 | static constexpr bool isCompressible(int phaseIdx) | |
| 165 | { | ||
| 166 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
| 167 | // gases are always compressible | ||
| 168 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
|
3 | if (phaseIdx == gPhaseIdx) |
| 169 | return true; | ||
| 170 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
|
2 | else if (phaseIdx == wPhaseIdx) |
| 171 | // the water component decides for the water phase... | ||
| 172 | return H2O::liquidIsCompressible(); | ||
| 173 | |||
| 174 | // the NAPL component decides for the napl phase... | ||
| 175 | return NAPL::liquidIsCompressible(); | ||
| 176 | } | ||
| 177 | |||
| 178 | /*! | ||
| 179 | * \brief Return the human readable name of a phase (used in indices) | ||
| 180 | * \param phaseIdx The index of the fluid phase to consider | ||
| 181 | */ | ||
| 182 | 39 | static std::string phaseName(int phaseIdx) | |
| 183 | { | ||
| 184 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 39 times.
|
39 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
| 185 |
3/3✓ Branch 0 taken 13 times.
✓ Branch 1 taken 13 times.
✓ Branch 2 taken 13 times.
|
39 | switch (phaseIdx) |
| 186 | { | ||
| 187 | 13 | case wPhaseIdx: return IOName::aqueousPhase(); | |
| 188 | 13 | case nPhaseIdx: return IOName::naplPhase(); | |
| 189 | 13 | case gPhaseIdx: return IOName::gaseousPhase(); | |
| 190 | } | ||
| 191 | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | ||
| 192 | } | ||
| 193 | |||
| 194 | /*! | ||
| 195 | * \brief Return the human readable name of a component (used in indices) | ||
| 196 | * \param compIdx The index of the component to consider | ||
| 197 | */ | ||
| 198 | 21 | static std::string componentName(int compIdx) | |
| 199 | { | ||
| 200 |
3/4✓ Branch 0 taken 7 times.
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
|
21 | switch (compIdx) { |
| 201 | 7 | case H2OIdx: return H2O::name(); | |
| 202 | 7 | case AirIdx: return Air::name(); | |
| 203 | 7 | case NAPLIdx: return NAPL::name(); | |
| 204 | } | ||
| 205 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); | |
| 206 | } | ||
| 207 | |||
| 208 | /*! | ||
| 209 | * \brief Return the molar mass of a component in \f$\mathrm{[kg/mol]}\f$ | ||
| 210 | * \param compIdx The index of the component to consider | ||
| 211 | */ | ||
| 212 |
2/3✓ Branch 0 taken 33 times.
✓ Branch 1 taken 237080022 times.
✗ Branch 2 not taken.
|
246194994 | static Scalar molarMass(int compIdx) |
| 213 | { | ||
| 214 | switch (compIdx) { | ||
| 215 | case H2OIdx: return H2O::molarMass(); | ||
| 216 | case AirIdx: return Air::molarMass(); | ||
| 217 | case NAPLIdx: return NAPL::molarMass(); | ||
| 218 | } | ||
| 219 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); | |
| 220 | } | ||
| 221 | |||
| 222 | using Base<Scalar, ThisType>::density; | ||
| 223 | /*! | ||
| 224 | * \brief Given a phase's composition, temperature, pressure, and | ||
| 225 | * the partial pressures of all components, return its | ||
| 226 | * density \f$\mathrm{[kg/m^3]}\f$. | ||
| 227 | * | ||
| 228 | * We apply Eq. (7) | ||
| 229 | * in Class et al. (2002a) \cite A3:class:2002b <BR> | ||
| 230 | * for the water density. | ||
| 231 | * | ||
| 232 | * \param fluidState The fluid state | ||
| 233 | * \param phaseIdx The index of the phase | ||
| 234 | */ | ||
| 235 | template <class FluidState> | ||
| 236 | 9889974 | static Scalar density(const FluidState &fluidState, int phaseIdx) | |
| 237 | { | ||
| 238 |
2/2✓ Branch 0 taken 3296658 times.
✓ Branch 1 taken 3296659 times.
|
6593317 | if (phaseIdx == wPhaseIdx) { |
| 239 | // This assumes each gas molecule displaces exactly one | ||
| 240 | // molecule in the liquid. | ||
| 241 | 3296658 | return H2O::liquidMolarDensity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)) | |
| 242 | 3296658 | * (H2O::molarMass()*fluidState.moleFraction(wPhaseIdx, H2OIdx) | |
| 243 | 3296658 | + Air::molarMass()*fluidState.moleFraction(wPhaseIdx, AirIdx) | |
| 244 | 3296658 | + NAPL::molarMass()*fluidState.moleFraction(wPhaseIdx, NAPLIdx)); | |
| 245 | } | ||
| 246 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3296658 times.
|
3296659 | else if (phaseIdx == nPhaseIdx) { |
| 247 | // assume pure NAPL for the NAPL phase | ||
| 248 | 3296658 | Scalar pressure = NAPL::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100; | |
| 249 | 3296658 | return NAPL::liquidDensity(fluidState.temperature(phaseIdx), pressure); | |
| 250 | } | ||
| 251 | |||
| 252 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3296658 times.
|
3296658 | assert (phaseIdx == gPhaseIdx); |
| 253 | 3296658 | Scalar pH2O = | |
| 254 | 3296658 | fluidState.moleFraction(gPhaseIdx, H2OIdx) * | |
| 255 | 3296658 | fluidState.pressure(gPhaseIdx); | |
| 256 | 3296658 | Scalar pAir = | |
| 257 | 3296658 | fluidState.moleFraction(gPhaseIdx, AirIdx) * | |
| 258 | 3296658 | fluidState.pressure(gPhaseIdx); | |
| 259 | 3296658 | Scalar pNAPL = | |
| 260 | 3296658 | fluidState.moleFraction(gPhaseIdx, NAPLIdx) * | |
| 261 | 3296658 | fluidState.pressure(gPhaseIdx); | |
| 262 | 3296658 | return H2O::gasDensity(fluidState.temperature(phaseIdx), pH2O) | |
| 263 | 3296658 | + Air::gasDensity(fluidState.temperature(phaseIdx), pAir) | |
| 264 | 3296658 | + NAPL::gasDensity(fluidState.temperature(phaseIdx), pNAPL); | |
| 265 | } | ||
| 266 | |||
| 267 | using Base<Scalar, ThisType>::molarDensity; | ||
| 268 | //! \copydoc Base<Scalar,ThisType>::molarDensity(const FluidState&,int) | ||
| 269 | template <class FluidState> | ||
| 270 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3296657 times.
|
9889974 | static Scalar molarDensity(const FluidState &fluidState, int phaseIdx) |
| 271 | { | ||
| 272 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3296657 times.
|
9889974 | Scalar temperature = fluidState.temperature(phaseIdx); |
| 273 | 3296660 | Scalar pressure = fluidState.pressure(phaseIdx); | |
| 274 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3296659 times.
|
3296660 | if (phaseIdx == nPhaseIdx) |
| 275 | { | ||
| 276 | 3296658 | return NAPL::liquidMolarDensity(temperature, pressure); | |
| 277 | } | ||
| 278 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3296658 times.
|
3296659 | else if (phaseIdx == wPhaseIdx) |
| 279 | { // This assumes each gas molecule displaces exactly one | ||
| 280 | // molecule in the liquid. | ||
| 281 | 3296658 | return H2O::liquidMolarDensity(temperature, pressure); | |
| 282 | } | ||
| 283 | else | ||
| 284 | { | ||
| 285 | 3296658 | return H2O::gasMolarDensity(temperature, fluidState.partialPressure(gPhaseIdx, H2OIdx)) | |
| 286 | 3296658 | + NAPL::gasMolarDensity(temperature, fluidState.partialPressure(gPhaseIdx, NAPLIdx)) | |
| 287 | 3296658 | + Air::gasMolarDensity(temperature, fluidState.partialPressure(gPhaseIdx, AirIdx)); | |
| 288 | } | ||
| 289 | } | ||
| 290 | |||
| 291 | using Base<Scalar, ThisType>::viscosity; | ||
| 292 | /*! | ||
| 293 | * \brief Return the viscosity of a phase \f$\mathrm{[Pa s]}\f$. | ||
| 294 | * \param fluidState The fluid state | ||
| 295 | * \param phaseIdx The index of the phase to consider | ||
| 296 | * | ||
| 297 | * \todo Check the parameter phiCAW for the xylene case and give a physical meaningful name | ||
| 298 | */ | ||
| 299 | template <class FluidState> | ||
| 300 | 9889974 | static Scalar viscosity(const FluidState &fluidState, | |
| 301 | int phaseIdx) | ||
| 302 | { | ||
| 303 |
2/2✓ Branch 0 taken 3296658 times.
✓ Branch 1 taken 6593316 times.
|
9889974 | if (phaseIdx == wPhaseIdx) { |
| 304 | // assume pure water viscosity | ||
| 305 | 3296658 | return H2O::liquidViscosity(fluidState.temperature(phaseIdx), | |
| 306 | fluidState.pressure(phaseIdx)); | ||
| 307 | } | ||
| 308 |
2/2✓ Branch 0 taken 3296658 times.
✓ Branch 1 taken 3296658 times.
|
6593316 | else if (phaseIdx == nPhaseIdx) { |
| 309 | // assume pure NAPL viscosity | ||
| 310 | 3296658 | return NAPL::liquidViscosity(fluidState.temperature(phaseIdx), | |
| 311 | 3296658 | fluidState.pressure(phaseIdx)); | |
| 312 | } | ||
| 313 | |||
| 314 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3296658 times.
|
3296658 | assert (phaseIdx == gPhaseIdx); |
| 315 | |||
| 316 | /* Wilke method. See: | ||
| 317 | * | ||
| 318 | * See: R. Reid, et al.: The Properties of Gases and Liquids, | ||
| 319 | * 4th edition, McGraw-Hill, 1987, 407-410 | ||
| 320 | * 5th edition, McGraw-Hill, 20001, p. 9.21/22 | ||
| 321 | * | ||
| 322 | * in this case, we use a simplified version in order to avoid | ||
| 323 | * computationally costly evaluation of sqrt and pow functions and | ||
| 324 | * divisions | ||
| 325 | * -- compare e.g. with Promo Class p. 32/33 | ||
| 326 | */ | ||
| 327 | Scalar muResult; | ||
| 328 | 3296660 | const Scalar mu[numComponents] = { | |
| 329 | 3296658 | h2oGasViscosityInMixture(fluidState.temperature(phaseIdx),fluidState.pressure(phaseIdx)), | |
| 330 | 3296658 | Air::gasViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)), | |
| 331 | 3296658 | NAPL::gasViscosity(fluidState.temperature(phaseIdx), NAPL::vaporPressure(fluidState.temperature(phaseIdx))) | |
| 332 | }; | ||
| 333 | // molar masses | ||
| 334 | const Scalar M[numComponents] = { | ||
| 335 | H2O::molarMass(), | ||
| 336 | Air::molarMass(), | ||
| 337 | NAPL::molarMass() | ||
| 338 | }; | ||
| 339 | |||
| 340 | 3296658 | Scalar muAW = mu[AirIdx]*fluidState.moleFraction(gPhaseIdx, AirIdx) | |
| 341 | 3296658 | + mu[H2OIdx]*fluidState.moleFraction(gPhaseIdx, H2OIdx) | |
| 342 | 3296658 | / (fluidState.moleFraction(gPhaseIdx, AirIdx) | |
| 343 | 3296658 | + fluidState.moleFraction(gPhaseIdx, H2OIdx)); | |
| 344 | 3296658 | Scalar xAW = fluidState.moleFraction(gPhaseIdx, AirIdx) | |
| 345 | 3296657 | + fluidState.moleFraction(gPhaseIdx, H2OIdx); | |
| 346 | |||
| 347 | 3296658 | Scalar MAW = (fluidState.moleFraction(gPhaseIdx, AirIdx)*Air::molarMass() | |
| 348 | 3296658 | + fluidState.moleFraction(gPhaseIdx, H2OIdx)*H2O::molarMass()) | |
| 349 | / xAW; | ||
| 350 | |||
| 351 | 3296658 | Scalar phiCAW = 0.3; // simplification for this particular system | |
| 352 | /* actually like this | ||
| 353 | * using std::sqrt; | ||
| 354 | * using std::pow; | ||
| 355 | * Scalar phiCAW = pow(1.+sqrt(mu[NAPLIdx]/muAW)*pow(MAW/M[NAPLIdx],0.25),2) | ||
| 356 | * / sqrt(8.*(1.+M[NAPLIdx]/MAW)); | ||
| 357 | */ | ||
| 358 | 3296658 | Scalar phiAWC = phiCAW * muAW*M[NAPLIdx]/(mu[NAPLIdx]*MAW); | |
| 359 | |||
| 360 | 3296658 | muResult = (xAW*muAW)/(xAW+fluidState.moleFraction(gPhaseIdx, NAPLIdx)*phiAWC) | |
| 361 | 3296658 | + (fluidState.moleFraction(gPhaseIdx, NAPLIdx) * mu[NAPLIdx]) | |
| 362 | 3296658 | / (fluidState.moleFraction(gPhaseIdx, NAPLIdx) + xAW*phiCAW); | |
| 363 | 3296658 | return muResult; | |
| 364 | } | ||
| 365 | |||
| 366 | |||
| 367 | using Base<Scalar, ThisType>::diffusionCoefficient; | ||
| 368 | //! \copydoc Base<Scalar,ThisType>::diffusionCoefficient(const FluidState&,int,int) | ||
| 369 | template <class FluidState> | ||
| 370 |
2/2✓ Branch 0 taken 6593314 times.
✓ Branch 1 taken 13186628 times.
|
19779951 | static Scalar diffusionCoefficient(const FluidState &fluidState, |
| 371 | int phaseIdx, | ||
| 372 | int compIdx) | ||
| 373 | { | ||
| 374 |
2/2✓ Branch 0 taken 6593314 times.
✓ Branch 1 taken 13186628 times.
|
19779951 | Scalar temperature = fluidState.temperature(phaseIdx); |
| 375 | 19779951 | Scalar pressure = fluidState.pressure(phaseIdx); | |
| 376 |
2/2✓ Branch 0 taken 6593317 times.
✓ Branch 1 taken 13186634 times.
|
19779951 | if (phaseIdx==gPhaseIdx) { |
| 377 | 6593317 | Scalar diffAC = BinaryCoeff::Air_Xylene::gasDiffCoeff(temperature, pressure); | |
| 378 |
2/2✓ Branch 1 taken 3296657 times.
✓ Branch 2 taken 3296657 times.
|
6593317 | Scalar diffWC = BinaryCoeff::H2O_Xylene::gasDiffCoeff(temperature, pressure); |
| 379 |
2/2✓ Branch 0 taken 3296657 times.
✓ Branch 1 taken 3296657 times.
|
6593317 | Scalar diffAW = BinaryCoeff::H2O_Air::gasDiffCoeff(temperature, pressure); |
| 380 | |||
| 381 | 6593317 | const Scalar xga = fluidState.moleFraction(gPhaseIdx, AirIdx); | |
| 382 | 6593317 | const Scalar xgw = fluidState.moleFraction(gPhaseIdx, H2OIdx); | |
| 383 | 6593317 | const Scalar xgc = fluidState.moleFraction(gPhaseIdx, NAPLIdx); | |
| 384 | |||
| 385 |
2/2✓ Branch 0 taken 3296658 times.
✓ Branch 1 taken 3296659 times.
|
6593317 | if (compIdx==NAPLIdx) |
| 386 | 3296658 | return (1.- xgw)/(xga/diffAW + xgc/diffWC); | |
| 387 |
2/2✓ Branch 0 taken 3296658 times.
✓ Branch 1 taken 1 times.
|
3296659 | else if (compIdx==H2OIdx) |
| 388 | 3296658 | return (1.- xgc)/(xgw/diffWC + xga/diffAC); | |
| 389 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
1 | else if (compIdx==AirIdx) |
| 390 |
10/20✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 1 times.
✗ Branch 29 not taken.
|
5 | DUNE_THROW(Dune::InvalidStateException, |
| 391 | "Diffusivity of Air in the gas phase " | ||
| 392 | "is constraint by sum of diffusive fluxes = 0 !\n"); | ||
| 393 |
2/2✓ Branch 0 taken 6593317 times.
✓ Branch 1 taken 6593317 times.
|
13186634 | } else if (phaseIdx==wPhaseIdx){ |
| 394 | 6593317 | Scalar diffACl = BinaryCoeff::Air_Xylene::liquidDiffCoeff(temperature, pressure); | |
| 395 |
2/4✓ Branch 0 taken 3296657 times.
✓ Branch 1 taken 3296657 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
6593317 | Scalar diffWCl = BinaryCoeff::H2O_Xylene::liquidDiffCoeff(temperature, pressure); |
| 396 |
2/4✓ Branch 0 taken 3296657 times.
✓ Branch 1 taken 3296657 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
6593317 | Scalar diffAWl = BinaryCoeff::H2O_Air::liquidDiffCoeff(temperature, pressure); |
| 397 | |||
| 398 | 6593317 | Scalar xwa = fluidState.moleFraction(wPhaseIdx, AirIdx); | |
| 399 | 6593317 | Scalar xww = fluidState.moleFraction(wPhaseIdx, H2OIdx); | |
| 400 | 6593317 | Scalar xwc = fluidState.moleFraction(wPhaseIdx, NAPLIdx); | |
| 401 | |||
| 402 | Scalar diffCont; | ||
| 403 | |||
| 404 |
3/4✓ Branch 0 taken 3296658 times.
✓ Branch 1 taken 3296658 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
6593317 | switch (compIdx) { |
| 405 | 3296658 | case NAPLIdx: | |
| 406 | 3296658 | diffCont = (1.- xww)/(xwa/diffAWl + xwc/diffWCl); | |
| 407 | 3296658 | return diffCont; | |
| 408 | 3296658 | case AirIdx: | |
| 409 | 3296658 | diffCont = (1.- xwc)/(xww/diffWCl + xwa/diffACl); | |
| 410 | 3296658 | return diffCont; | |
| 411 | 1 | case H2OIdx: | |
| 412 |
10/20✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 1 times.
✗ Branch 29 not taken.
|
4 | DUNE_THROW(Dune::InvalidStateException, |
| 413 | "Diffusivity of water in the water phase " | ||
| 414 | "is constraint by sum of diffusive fluxes = 0 !\n"); | ||
| 415 | } | ||
| 416 | } else if (phaseIdx==nPhaseIdx) { | ||
| 417 | return 0.0; | ||
| 418 | } | ||
| 419 | return 0; | ||
| 420 | } | ||
| 421 | |||
| 422 | using Base<Scalar, ThisType>::binaryDiffusionCoefficient; | ||
| 423 | //! \copydoc Base<Scalar,ThisType>::binaryDiffusionCoefficient(const FluidState&,int,int,int) | ||
| 424 | template <class FluidState> | ||
| 425 | 27 | static Scalar binaryDiffusionCoefficient(const FluidState &fluidState, | |
| 426 | int phaseIdx, | ||
| 427 | int compIIdx, | ||
| 428 | int compJIdx) | ||
| 429 | { | ||
| 430 |
10/20✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 27 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 27 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 27 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 27 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 27 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 27 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 27 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 27 times.
✗ Branch 29 not taken.
|
108 | DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OAirXylene::binaryDiffusionCoefficient()"); |
| 431 | } | ||
| 432 | |||
| 433 | using Base<Scalar, ThisType>::fugacityCoefficient; | ||
| 434 | /*! | ||
| 435 | * \brief Returns the fugacity coefficient \f$\mathrm{[-]}\f$ of a component in a | ||
| 436 | * phase. | ||
| 437 | * \param fluidState The fluid state | ||
| 438 | * \param phaseIdx The index of the phase to consider | ||
| 439 | * \param compIdx The index of the component to consider | ||
| 440 | * | ||
| 441 | * In this case, things are actually pretty simple. We have an ideal | ||
| 442 | * solution. Thus, the fugacity coefficient is 1 in the gas phase | ||
| 443 | * (fugacity equals the partial pressure of the component in the gas phase) | ||
| 444 | * respectively in the liquid phases it is the Henry coefficients | ||
| 445 | * divided by pressure. | ||
| 446 | */ | ||
| 447 | template <class FluidState> | ||
| 448 | 6587092 | static Scalar fugacityCoefficient(const FluidState &fluidState, | |
| 449 | int phaseIdx, | ||
| 450 | int compIdx) | ||
| 451 | { | ||
| 452 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 36263236 times.
|
36263236 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
| 453 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 36263236 times.
|
36263236 | assert(0 <= compIdx && compIdx < numComponents); |
| 454 | |||
| 455 |
2/2✓ Branch 0 taken 13186628 times.
✓ Branch 1 taken 23076599 times.
|
42850319 | Scalar T = fluidState.temperature(phaseIdx); |
| 456 | 3296666 | Scalar p = fluidState.pressure(phaseIdx); | |
| 457 | |||
| 458 |
2/2✓ Branch 0 taken 13186631 times.
✓ Branch 1 taken 23076605 times.
|
36263236 | if (phaseIdx == wPhaseIdx) { |
| 459 |
2/2✓ Branch 0 taken 3296658 times.
✓ Branch 1 taken 9889973 times.
|
13186631 | if (compIdx == H2OIdx) |
| 460 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 3290426 times.
|
6587084 | return H2O::vaporPressure(T)/p; |
| 461 |
2/2✓ Branch 0 taken 6593315 times.
✓ Branch 1 taken 3296658 times.
|
9889973 | else if (compIdx == AirIdx) |
| 462 | 6593315 | return BinaryCoeff::H2O_Air::henry(T)/p; | |
| 463 | else if (compIdx == NAPLIdx) | ||
| 464 | 6593315 | return BinaryCoeff::H2O_Xylene::henry(T)/p; | |
| 465 | } | ||
| 466 | |||
| 467 | // for the NAPL phase, we assume currently that nothing is | ||
| 468 | // dissolved. this means that the affinity of the NAPL | ||
| 469 | // component to the NAPL phase is much higher than for the | ||
| 470 | // other components, i.e. the fugacity coefficient is much | ||
| 471 | // smaller. | ||
| 472 |
2/2✓ Branch 0 taken 13186631 times.
✓ Branch 1 taken 9889974 times.
|
23076605 | if (phaseIdx == nPhaseIdx) { |
| 473 | 13186631 | Scalar phiNapl = NAPL::vaporPressure(T)/p; | |
| 474 |
2/2✓ Branch 0 taken 6593316 times.
✓ Branch 1 taken 6593315 times.
|
13186631 | if (compIdx == NAPLIdx) |
| 475 | return phiNapl; | ||
| 476 |
2/2✓ Branch 0 taken 3296658 times.
✓ Branch 1 taken 3296658 times.
|
6593316 | else if (compIdx == AirIdx) |
| 477 | 3296658 | return 1e6*phiNapl; | |
| 478 | else if (compIdx == H2OIdx) | ||
| 479 | 3296658 | return 1e6*phiNapl; | |
| 480 | } | ||
| 481 | |||
| 482 | // for the gas phase, assume an ideal gas when it comes to | ||
| 483 | // fugacity (-> fugacity == partial pressure) | ||
| 484 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9889974 times.
|
9889974 | assert(phaseIdx == gPhaseIdx); |
| 485 | return 1.0; | ||
| 486 | } | ||
| 487 | |||
| 488 | template <class FluidState> | ||
| 489 | static Scalar kelvinVaporPressure(const FluidState &fluidState, | ||
| 490 | const int phaseIdx, | ||
| 491 | const int compIdx) | ||
| 492 | { | ||
| 493 | DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OAirXylene::kelvinVaporPressure()"); | ||
| 494 | } | ||
| 495 | |||
| 496 | using Base<Scalar, ThisType>::enthalpy; | ||
| 497 | /*! | ||
| 498 | * \brief Given all mole fractions in a phase, return the specific | ||
| 499 | * phase enthalpy \f$\mathrm{[J/kg]}\f$. | ||
| 500 | * \param fluidState The fluid state | ||
| 501 | * \param phaseIdx The index of the phase to consider | ||
| 502 | * | ||
| 503 | * \note This system neglects the contribution of gas-molecules in the liquid phase. | ||
| 504 | * This contribution is probably not big. Somebody would have to find out the enthalpy of solution for this system. ... | ||
| 505 | */ | ||
| 506 | template <class FluidState> | ||
| 507 | 9889974 | static Scalar enthalpy(const FluidState &fluidState, | |
| 508 | int phaseIdx) | ||
| 509 | { | ||
| 510 |
2/2✓ Branch 0 taken 3296658 times.
✓ Branch 1 taken 6593316 times.
|
9889974 | if (phaseIdx == wPhaseIdx) { |
| 511 | 3296658 | return H2O::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); | |
| 512 | } | ||
| 513 |
2/2✓ Branch 0 taken 3296658 times.
✓ Branch 1 taken 3296658 times.
|
6593316 | else if (phaseIdx == nPhaseIdx) { |
| 514 | 3296658 | return NAPL::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); | |
| 515 | } | ||
| 516 |
1/2✓ Branch 0 taken 3296658 times.
✗ Branch 1 not taken.
|
3296658 | else if (phaseIdx == gPhaseIdx) { // gas phase enthalpy depends strongly on composition |
| 517 | 3296658 | Scalar hgc = NAPL::gasEnthalpy(fluidState.temperature(phaseIdx), | |
| 518 | fluidState.pressure(phaseIdx)); | ||
| 519 | 3296658 | Scalar hgw = H2O::gasEnthalpy(fluidState.temperature(phaseIdx), | |
| 520 | fluidState.pressure(phaseIdx)); | ||
| 521 | // pressure is only a dummy here (not dependent on pressure, just temperature) | ||
| 522 | 3296658 | Scalar hga = Air::gasEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); | |
| 523 | |||
| 524 | 3296658 | Scalar result = 0; | |
| 525 | 3296658 | result += hgw * fluidState.massFraction(gPhaseIdx, H2OIdx); | |
| 526 | 3296658 | result += hga * fluidState.massFraction(gPhaseIdx, AirIdx); | |
| 527 | 3296658 | result += hgc * fluidState.massFraction(gPhaseIdx, NAPLIdx); | |
| 528 | |||
| 529 | 3296658 | return result; | |
| 530 | } | ||
| 531 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | |
| 532 | } | ||
| 533 | |||
| 534 | /*! | ||
| 535 | * \brief Returns the specific enthalpy \f$\mathrm{[J/kg]}\f$ of a component in a specific phase | ||
| 536 | * \param fluidState The fluid state | ||
| 537 | * \param phaseIdx The index of the phase | ||
| 538 | * \param componentIdx The index of the component | ||
| 539 | */ | ||
| 540 | template <class FluidState> | ||
| 541 | static Scalar componentEnthalpy(const FluidState& fluidState, int phaseIdx, int componentIdx) | ||
| 542 | { | ||
| 543 | const Scalar T = fluidState.temperature(phaseIdx); | ||
| 544 | const Scalar p = fluidState.pressure(phaseIdx); | ||
| 545 | |||
| 546 | if (phaseIdx == wPhaseIdx) | ||
| 547 | { | ||
| 548 | if (componentIdx == H2OIdx) | ||
| 549 | return H2O::liquidEnthalpy(T, p); | ||
| 550 | else if (componentIdx == NAPLIdx) | ||
| 551 | DUNE_THROW(Dune::NotImplemented, "The component enthalpy for NAPL in water is not implemented."); | ||
| 552 | else if (componentIdx == AirIdx) | ||
| 553 | DUNE_THROW(Dune::NotImplemented, "The component enthalpy for Air in water is not implemented."); | ||
| 554 | DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << componentIdx); | ||
| 555 | } | ||
| 556 | else if (phaseIdx == nPhaseIdx) | ||
| 557 | { | ||
| 558 | if (componentIdx == H2OIdx) | ||
| 559 | DUNE_THROW(Dune::NotImplemented, "The component enthalpy for water in NAPL is not implemented."); | ||
| 560 | else if (componentIdx == NAPLIdx) | ||
| 561 | return NAPL::liquidEnthalpy(T, p); | ||
| 562 | else if (componentIdx == AirIdx) | ||
| 563 | DUNE_THROW(Dune::NotImplemented, "The component enthalpy for air in NAPL is not implemented."); | ||
| 564 | DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << componentIdx); | ||
| 565 | } | ||
| 566 | else if (phaseIdx == gPhaseIdx) | ||
| 567 | { | ||
| 568 | if (componentIdx == H2OIdx) | ||
| 569 | return H2O::gasEnthalpy(T, p); | ||
| 570 | else if (componentIdx == NAPLIdx) | ||
| 571 | return NAPL::gasEnthalpy(T, p); | ||
| 572 | else if (componentIdx == AirIdx) | ||
| 573 | return Air::gasEnthalpy(T,p); | ||
| 574 | DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << componentIdx); | ||
| 575 | } | ||
| 576 | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | ||
| 577 | } | ||
| 578 | |||
| 579 | using Base<Scalar, ThisType>::heatCapacity; | ||
| 580 | //! \copydoc Base<Scalar,ThisType>::heatCapacity(const FluidState&,int) | ||
| 581 | template <class FluidState> | ||
| 582 | 3 | static Scalar heatCapacity(const FluidState &fluidState, | |
| 583 | int phaseIdx) | ||
| 584 | { | ||
| 585 |
10/20✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 3 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 3 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 3 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 3 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 3 times.
✗ Branch 29 not taken.
|
12 | DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OAirXylene::heatCapacity()"); |
| 586 | } | ||
| 587 | |||
| 588 | using Base<Scalar, ThisType>::thermalConductivity; | ||
| 589 | //! \copydoc Base<Scalar,ThisType>::thermalConductivity(const FluidState&,int) | ||
| 590 | template <class FluidState> | ||
| 591 |
2/2✓ Branch 0 taken 3296657 times.
✓ Branch 1 taken 6593314 times.
|
9889974 | static Scalar thermalConductivity(const FluidState &fluidState, |
| 592 | int phaseIdx) | ||
| 593 | { | ||
| 594 |
2/2✓ Branch 0 taken 3296657 times.
✓ Branch 1 taken 6593314 times.
|
9889974 | const Scalar temperature = fluidState.temperature(phaseIdx) ; |
| 595 | 9889974 | const Scalar pressure = fluidState.pressure(phaseIdx); | |
| 596 |
2/2✓ Branch 0 taken 3296658 times.
✓ Branch 1 taken 6593316 times.
|
9889974 | if (phaseIdx == wPhaseIdx) |
| 597 | { | ||
| 598 | 3296658 | return H2O::liquidThermalConductivity(temperature, pressure); | |
| 599 | } | ||
| 600 |
2/2✓ Branch 0 taken 3296658 times.
✓ Branch 1 taken 3296658 times.
|
6593316 | else if (phaseIdx == nPhaseIdx) |
| 601 | { | ||
| 602 | return NAPL::liquidThermalConductivity(temperature, pressure); | ||
| 603 | } | ||
| 604 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3296658 times.
|
3296658 | else if (phaseIdx == gPhaseIdx) |
| 605 | { | ||
| 606 | return Air::gasThermalConductivity(temperature, pressure); | ||
| 607 | } | ||
| 608 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | |
| 609 | } | ||
| 610 | |||
| 611 | private: | ||
| 612 | static Scalar waterPhaseDensity_(Scalar T, Scalar pw, Scalar xww, Scalar xwa, Scalar xwc) | ||
| 613 | { | ||
| 614 | Scalar rholH2O = H2O::liquidDensity(T, pw); | ||
| 615 | Scalar clH2O = rholH2O/H2O::molarMass(); | ||
| 616 | |||
| 617 | // this assumes each dissolved molecule displaces exactly one | ||
| 618 | // water molecule in the liquid | ||
| 619 | return | ||
| 620 | clH2O*(xww*H2O::molarMass() + xwa*Air::molarMass() + xwc*NAPL::molarMass()); | ||
| 621 | } | ||
| 622 | |||
| 623 | static Scalar gasPhaseDensity_(Scalar T, Scalar pg, Scalar xgw, Scalar xga, Scalar xgc) | ||
| 624 | { | ||
| 625 | return H2O::gasDensity(T, pg*xgw) + Air::gasDensity(T, pg*xga) + NAPL::gasDensity(T, pg*xgc); | ||
| 626 | } | ||
| 627 | |||
| 628 | static Scalar NAPLPhaseDensity_(Scalar T, Scalar pn) | ||
| 629 | { | ||
| 630 | return | ||
| 631 | NAPL::liquidDensity(T, pn); | ||
| 632 | } | ||
| 633 | |||
| 634 | }; | ||
| 635 | |||
| 636 | } // end namespace Dumux::FluidSystems | ||
| 637 | |||
| 638 | #endif | ||
| 639 |