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