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