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