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 | /*! | ||
9 | * \file | ||
10 | * \ingroup FluidSystems | ||
11 | * \brief @copybrief Dumux::FluidSystems::H2OAir | ||
12 | */ | ||
13 | #ifndef DUMUX_H2O_AIR_SYSTEM_HH | ||
14 | #define DUMUX_H2O_AIR_SYSTEM_HH | ||
15 | |||
16 | #include <cassert> | ||
17 | #include <iomanip> | ||
18 | |||
19 | #include <dumux/material/idealgas.hh> | ||
20 | #include <dumux/material/fluidsystems/base.hh> | ||
21 | |||
22 | #include <dumux/material/binarycoefficients/h2o_air.hh> | ||
23 | #include <dumux/material/components/air.hh> | ||
24 | #include <dumux/material/components/tabulatedcomponent.hh> | ||
25 | #include <dumux/material/components/h2o.hh> | ||
26 | |||
27 | #include <dumux/common/exceptions.hh> | ||
28 | |||
29 | #include <dumux/io/name.hh> | ||
30 | |||
31 | namespace Dumux { | ||
32 | namespace FluidSystems { | ||
33 | /*! | ||
34 | * \ingroup FluidSystems | ||
35 | * \brief Policy for the H2O-air fluid system | ||
36 | */ | ||
37 | template<bool fastButSimplifiedRelations = false> | ||
38 | struct H2OAirDefaultPolicy | ||
39 | { | ||
40 | static constexpr bool useH2ODensityAsLiquidMixtureDensity() { return fastButSimplifiedRelations; } | ||
41 | static constexpr bool useIdealGasDensity() { return fastButSimplifiedRelations; } | ||
42 | static constexpr bool useAirViscosityAsGasMixtureViscosity() { return fastButSimplifiedRelations; } | ||
43 | }; | ||
44 | |||
45 | /*! | ||
46 | * \ingroup FluidSystems | ||
47 | * | ||
48 | * \brief A compositional two-phase fluid system with water and air as | ||
49 | * components in both, the liquid and the gas phase. | ||
50 | * | ||
51 | * This fluidsystem features gas and liquid phases of distilled water | ||
52 | * \f$(\mathrm{H_2O})\f$) and air (Pseudo component composed of \f$\mathrm{79\%\;N_2}\f$, | ||
53 | * \f$\mathrm{20\%\;O_2}\f$ and \f$\mathrm{1\%\;Ar}\f$) as components. It is applied by | ||
54 | * default with the tabulated version of water of the IAPWS-formulation. | ||
55 | */ | ||
56 | template <class Scalar, | ||
57 | class H2Otype = Components::TabulatedComponent<Components::H2O<Scalar> >, | ||
58 | class Policy = H2OAirDefaultPolicy<>, | ||
59 | bool useKelvinVaporPressure = false> | ||
60 | class H2OAir | ||
61 | : public Base<Scalar, H2OAir<Scalar, H2Otype, Policy> > | ||
62 | { | ||
63 | using ThisType = H2OAir<Scalar,H2Otype, Policy>; | ||
64 | using IdealGas = Dumux::IdealGas<Scalar>; | ||
65 | |||
66 | public: | ||
67 | using H2O = H2Otype; | ||
68 | using Air = Dumux::Components::Air<Scalar>; | ||
69 | |||
70 | static constexpr int numPhases = 2; //!< Number of phases in the fluid system | ||
71 | static constexpr int numComponents = 2; //!< Number of components in the fluid system | ||
72 | |||
73 | static constexpr int liquidPhaseIdx = 0; //!< index of the liquid phase | ||
74 | static constexpr int gasPhaseIdx = 1; //!< index of the gas phase | ||
75 | static constexpr int phase0Idx = liquidPhaseIdx; //!< index of the first phase | ||
76 | static constexpr int phase1Idx = gasPhaseIdx; //!< index of the second phase | ||
77 | |||
78 | static constexpr int H2OIdx = 0; //!< index of the first component | ||
79 | static constexpr int AirIdx = 1; //!< index of the second component | ||
80 | static constexpr int comp0Idx = H2OIdx; //!< index of the first component | ||
81 | static constexpr int comp1Idx = AirIdx; //!< index of the second component | ||
82 | static constexpr int liquidCompIdx = H2OIdx; //!< index of the liquid component | ||
83 | static constexpr int gasCompIdx = AirIdx; //!< index of the gas component | ||
84 | |||
85 | /*! | ||
86 | * \brief Return the human readable name of a phase | ||
87 | * | ||
88 | * \param phaseIdx index of the phase | ||
89 | */ | ||
90 | 1835 | static std::string phaseName(int phaseIdx) | |
91 | { | ||
92 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1822 times.
|
1835 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
93 |
2/3✓ Branch 0 taken 817 times.
✓ Branch 1 taken 1005 times.
✗ Branch 2 not taken.
|
1835 | switch (phaseIdx) |
94 | { | ||
95 | 824 | case liquidPhaseIdx: return IOName::liquidPhase(); | |
96 | 1011 | case gasPhaseIdx: return IOName::gaseousPhase(); | |
97 | } | ||
98 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | |
99 | } | ||
100 | |||
101 | /*! | ||
102 | * \brief Returns whether the fluids are miscible | ||
103 | */ | ||
104 | static constexpr bool isMiscible() | ||
105 | { return true; } | ||
106 | |||
107 | /*! | ||
108 | * \brief Return whether a phase is gaseous | ||
109 | * | ||
110 | * \param phaseIdx The index of the fluid phase to consider | ||
111 | */ | ||
112 | static constexpr bool isGas(int phaseIdx) | ||
113 | { | ||
114 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 162488276 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 4 times.
|
162488288 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
115 | return phaseIdx == gasPhaseIdx; | ||
116 | } | ||
117 | |||
118 | /*! | ||
119 | * \brief Returns true if and only if a fluid phase is assumed to | ||
120 | * be an ideal mixture. | ||
121 | * | ||
122 | * We define an ideal mixture as a fluid phase where the fugacity | ||
123 | * coefficients of all components times the pressure of the phase | ||
124 | * are independent on the fluid composition. This assumption is true | ||
125 | * if Henry's law and Raoult's law apply. If you are unsure what | ||
126 | * this function should return, it is safe to return false. The | ||
127 | * only damage done will be (slightly) increased computation times | ||
128 | * in some cases. | ||
129 | * | ||
130 | * \param phaseIdx The index of the fluid phase to consider | ||
131 | */ | ||
132 | static constexpr bool isIdealMixture(int phaseIdx) | ||
133 | { | ||
134 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 25191 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 25191 times.
|
50382 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
135 | // we assume Henry's and Raoult's laws for the water phase and | ||
136 | // and no interaction between gas molecules of different | ||
137 | // components, so all phases are ideal mixtures! | ||
138 | return true; | ||
139 | } | ||
140 | |||
141 | /*! | ||
142 | * \brief Returns true if and only if a fluid phase is assumed to | ||
143 | * be compressible. | ||
144 | * | ||
145 | * Compressible means that the partial derivative of the density | ||
146 | * to the fluid pressure is always larger than zero. | ||
147 | * | ||
148 | * \param phaseIdx The index of the fluid phase to consider | ||
149 | */ | ||
150 | static constexpr bool isCompressible(int phaseIdx) | ||
151 | { | ||
152 | assert(0 <= phaseIdx && phaseIdx < numPhases); | ||
153 | // ideal gases are always compressible | ||
154 |
4/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 1 times.
|
4 | if (phaseIdx == gasPhaseIdx) |
155 | return true; | ||
156 | // the water component decides for the liquid phase... | ||
157 | return H2O::liquidIsCompressible(); | ||
158 | } | ||
159 | |||
160 | /*! | ||
161 | * \brief Returns true if and only if a fluid phase is assumed to | ||
162 | * have a constant viscosity. | ||
163 | * | ||
164 | * \param phaseIdx The index of the fluid phase to consider | ||
165 | */ | ||
166 | static constexpr bool viscosityIsConstant(int phaseIdx) | ||
167 | { | ||
168 | // water decides for the liquid phase | ||
169 | if (phaseIdx == liquidPhaseIdx) | ||
170 | return H2O::liquidViscosityIsConstant(); | ||
171 | // air decides if policy is enabled | ||
172 | else if (phaseIdx == gasPhaseIdx && Policy::useAirViscosityAsGasMixtureViscosity()) | ||
173 | return Air::gasViscosityIsConstant(); | ||
174 | // in general it depends on the mixture | ||
175 | else | ||
176 | return false; | ||
177 | } | ||
178 | |||
179 | /*! | ||
180 | * \brief Returns true if and only if a fluid phase is assumed to | ||
181 | * be an ideal gas. | ||
182 | * | ||
183 | * \param phaseIdx The index of the fluid phase to consider | ||
184 | */ | ||
185 | static constexpr bool isIdealGas(int phaseIdx) | ||
186 | { | ||
187 | assert(0 <= phaseIdx && phaseIdx < numPhases); | ||
188 | |||
189 | // let the fluids decide | ||
190 | if (phaseIdx == gasPhaseIdx) | ||
191 | return H2O::gasIsIdeal() && Air::gasIsIdeal(); | ||
192 | return false; // not a gas | ||
193 | } | ||
194 | |||
195 | /**************************************** | ||
196 | * Component related static parameters | ||
197 | ****************************************/ | ||
198 | /*! | ||
199 | * \brief Return the human readable name of a component | ||
200 | * | ||
201 | * \param compIdx index of the component | ||
202 | */ | ||
203 | 673154 | static std::string componentName(int compIdx) | |
204 | { | ||
205 |
2/3✓ Branch 0 taken 399031 times.
✓ Branch 1 taken 274109 times.
✗ Branch 2 not taken.
|
673154 | switch (compIdx) |
206 | { | ||
207 | 399038 | case H2OIdx: return H2O::name(); | |
208 | 274116 | case AirIdx: return Air::name(); | |
209 | } | ||
210 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); | |
211 | } | ||
212 | |||
213 | /*! | ||
214 | * \brief Return the molar mass of a component \f$\mathrm{[kg/mol]}\f$. | ||
215 | * | ||
216 | * \param compIdx index of the component | ||
217 | */ | ||
218 | 868048102 | static Scalar molarMass(int compIdx) | |
219 | { | ||
220 |
2/3✓ Branch 0 taken 420231135 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 447816889 times.
|
868048102 | switch (compIdx) |
221 | { | ||
222 | case H2OIdx: return H2O::molarMass(); | ||
223 | 420231174 | case AirIdx: return Air::molarMass(); | |
224 | } | ||
225 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx); | |
226 | } | ||
227 | |||
228 | /*! | ||
229 | * \brief Critical temperature of a component \f$\mathrm{[K]}\f$. | ||
230 | * | ||
231 | * \param compIdx The index of the component to consider | ||
232 | */ | ||
233 | static Scalar criticalTemperature(int compIdx) | ||
234 | { | ||
235 | static const Scalar TCrit[] = { | ||
236 | H2O::criticalTemperature(), | ||
237 | Air::criticalTemperature() | ||
238 | }; | ||
239 | |||
240 | assert(0 <= compIdx && compIdx < numComponents); | ||
241 | return TCrit[compIdx]; | ||
242 | } | ||
243 | |||
244 | /*! | ||
245 | * \brief Critical pressure of a component \f$\mathrm{[Pa]}\f$. | ||
246 | * | ||
247 | * \param compIdx The index of the component to consider | ||
248 | */ | ||
249 | static Scalar criticalPressure(int compIdx) | ||
250 | { | ||
251 | static const Scalar pCrit[] = { | ||
252 | H2O::criticalPressure(), | ||
253 | Air::criticalPressure() | ||
254 | }; | ||
255 | |||
256 | assert(0 <= compIdx && compIdx < numComponents); | ||
257 | return pCrit[compIdx]; | ||
258 | } | ||
259 | |||
260 | /*! | ||
261 | * \brief Vapor pressure of a component \f$\mathrm{[Pa]}\f$. | ||
262 | * | ||
263 | * \param fluidState The fluid state | ||
264 | * \param compIdx The index of the component to consider | ||
265 | */ | ||
266 | template <class FluidState> | ||
267 | 1852922 | static Scalar vaporPressure(const FluidState &fluidState, int compIdx) | |
268 | { | ||
269 |
1/2✓ Branch 0 taken 1064379 times.
✗ Branch 1 not taken.
|
1852922 | if (compIdx == H2OIdx) |
270 | { | ||
271 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
1852922 | const auto t = fluidState.temperature(H2OIdx); |
272 | // cppcheck-suppress internalAstError | ||
273 | if constexpr (!useKelvinVaporPressure) | ||
274 | 1852922 | return H2O::vaporPressure(t); | |
275 | else | ||
276 | { | ||
277 | ✗ | const auto pc = (fluidState.wettingPhase() == (int) H2OIdx) | |
278 | ✗ | ? fluidState.pressure(AirIdx)-fluidState.pressure(H2OIdx) | |
279 | ✗ | : fluidState.pressure(H2OIdx)-fluidState.pressure(AirIdx); | |
280 | ✗ | return H2O::vaporPressure(t)*exp( -pc * molarMass(H2OIdx) | |
281 | ✗ | / density(fluidState, H2OIdx) | |
282 | ✗ | / (Dumux::Constants<Scalar>::R*t) ); | |
283 | } | ||
284 | } | ||
285 | ✗ | else if (compIdx == AirIdx) | |
286 | // return Air::vaporPressure(fluidState.temperature(AirIdx)); | ||
287 | ✗ | DUNE_THROW(Dune::NotImplemented, "Air::vaporPressure(t)"); | |
288 | else | ||
289 | ✗ | DUNE_THROW(Dune::NotImplemented, "Invalid component index " << compIdx); | |
290 | } | ||
291 | |||
292 | /*! | ||
293 | * \brief Molar volume of a component at the critical point \f$\mathrm{[m^3/mol]}\f$. | ||
294 | * | ||
295 | * \param compIdx The index of the component to consider | ||
296 | */ | ||
297 | static Scalar criticalMolarVolume(int compIdx) | ||
298 | { | ||
299 | DUNE_THROW(Dune::NotImplemented, | ||
300 | "H2OAirFluidSystem::criticalMolarVolume()"); | ||
301 | } | ||
302 | |||
303 | /*! | ||
304 | * \brief The acentric factor of a component \f$\mathrm{[-]}\f$. | ||
305 | * | ||
306 | * \param compIdx The index of the component to consider | ||
307 | */ | ||
308 | static Scalar acentricFactor(int compIdx) | ||
309 | { | ||
310 | static const Scalar accFac[] = { | ||
311 | H2O::acentricFactor(), | ||
312 | Air::acentricFactor() | ||
313 | }; | ||
314 | |||
315 | assert(0 <= compIdx && compIdx < numComponents); | ||
316 | return accFac[compIdx]; | ||
317 | } | ||
318 | |||
319 | /**************************************** | ||
320 | * thermodynamic relations | ||
321 | ****************************************/ | ||
322 | |||
323 | /*! | ||
324 | * \brief Initialize the fluid system's static parameters generically | ||
325 | * | ||
326 | * If a tabulated H2O component is used, we do our best to create | ||
327 | * tables that always work. | ||
328 | */ | ||
329 | static void init() | ||
330 | { | ||
331 |
2/4✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
|
41 | init(/*tempMin=*/273.15, |
332 | /*tempMax=*/623.15, | ||
333 | /*numTemp=*/100, | ||
334 | /*pMin=*/-10., | ||
335 | /*pMax=*/20e6, | ||
336 | /*numP=*/200); | ||
337 | } | ||
338 | |||
339 | /*! | ||
340 | * \brief Initialize the fluid system's static parameters using | ||
341 | * problem specific temperature and pressure ranges | ||
342 | * | ||
343 | * \param tempMin The minimum temperature used for tabulation of water \f$\mathrm{[K]}\f$ | ||
344 | * \param tempMax The maximum temperature used for tabulation of water\f$\mathrm{[K]}\f$ | ||
345 | * \param nTemp The number of ticks on the temperature axis of the table of water | ||
346 | * \param pressMin The minimum pressure used for tabulation of water \f$\mathrm{[Pa]}\f$ | ||
347 | * \param pressMax The maximum pressure used for tabulation of water \f$\mathrm{[Pa]}\f$ | ||
348 | * \param nPress The number of ticks on the pressure axis of the table of water | ||
349 | */ | ||
350 | 50 | static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp, | |
351 | Scalar pressMin, Scalar pressMax, unsigned nPress) | ||
352 | { | ||
353 | 50 | std::cout << "The H2O-air fluid system was configured with the following policy:\n"; | |
354 | 100 | std::cout << " - use H2O density as liquid mixture density: " << std::boolalpha << Policy::useH2ODensityAsLiquidMixtureDensity() << "\n"; | |
355 | 100 | std::cout << " - use ideal gas density: " << std::boolalpha << Policy::useIdealGasDensity() << "\n"; | |
356 | 100 | std::cout << " - use air viscosity as gas mixture viscosity: " << std::boolalpha << Policy::useAirViscosityAsGasMixtureViscosity() << std::endl; | |
357 | |||
358 | if (H2O::isTabulated) | ||
359 | { | ||
360 | 42 | H2O::init(tempMin, tempMax, nTemp, | |
361 | pressMin, pressMax, nPress); | ||
362 | } | ||
363 | 50 | } | |
364 | |||
365 | using Base<Scalar, ThisType>::density; | ||
366 | /*! | ||
367 | * \brief Given a phase's composition, temperature, pressure, and | ||
368 | * the partial pressures of all components, return its | ||
369 | * density \f$\mathrm{[kg/m^3]}\f$. | ||
370 | * | ||
371 | * If Policy::useH2ODensityAsLiquidMixtureDensity() == false, we apply Eq. (7) | ||
372 | * in Class et al. (2002a) \cite A3:class:2002b <BR> | ||
373 | * for the liquid density. | ||
374 | * | ||
375 | * \param phaseIdx index of the phase | ||
376 | * \param fluidState the fluid state | ||
377 | * | ||
378 | */ | ||
379 | template <class FluidState> | ||
380 | 186888440 | static Scalar density(const FluidState &fluidState, | |
381 | const int phaseIdx) | ||
382 | { | ||
383 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 177262673 times.
|
186888440 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
384 | |||
385 |
2/2✓ Branch 0 taken 76898896 times.
✓ Branch 1 taken 100363764 times.
|
186888440 | const Scalar T = fluidState.temperature(phaseIdx); |
386 |
2/2✓ Branch 0 taken 76898896 times.
✓ Branch 1 taken 100363764 times.
|
186888440 | const Scalar p = fluidState.pressure(phaseIdx); |
387 | |||
388 |
2/2✓ Branch 0 taken 76898903 times.
✓ Branch 1 taken 100363770 times.
|
186888440 | if (phaseIdx == phase0Idx) |
389 | { | ||
390 | if (Policy::useH2ODensityAsLiquidMixtureDensity()) | ||
391 | // assume pure water | ||
392 | 1577076 | return H2O::liquidDensity(T, p); | |
393 | else | ||
394 | { | ||
395 | // See: Eq. (7) in Class et al. (2002a) | ||
396 | // This assumes each gas molecule displaces exactly one | ||
397 | // molecule in the liquid. | ||
398 | 19103036 | return H2O::liquidMolarDensity(T, p) | |
399 | 19103036 | * (H2O::molarMass()*fluidState.moleFraction(liquidPhaseIdx, H2OIdx) | |
400 | 38206064 | + Air::molarMass()*fluidState.moleFraction(liquidPhaseIdx, AirIdx)); | |
401 | } | ||
402 | } | ||
403 | else if (phaseIdx == gasPhaseIdx) | ||
404 | { | ||
405 | if (Policy::useIdealGasDensity()) | ||
406 | // for the gas phase assume an ideal gas | ||
407 | { | ||
408 | 57328848 | const Scalar averageMolarMass = fluidState.averageMolarMass(gasPhaseIdx); | |
409 | 114657696 | return IdealGas::density(averageMolarMass, T, p); | |
410 | } | ||
411 | |||
412 | 52147982 | return H2O::gasDensity(T, fluidState.partialPressure(gasPhaseIdx, H2OIdx)) | |
413 | 55637470 | + Air::gasDensity(T, fluidState.partialPressure(gasPhaseIdx, AirIdx)); | |
414 | } | ||
415 | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | ||
416 | } | ||
417 | |||
418 | using Base<Scalar, ThisType>::molarDensity; | ||
419 | //! \copydoc Base<Scalar,ThisType>::molarDensity(const FluidState&,int) | ||
420 | template <class FluidState> | ||
421 | 65085106 | static Scalar molarDensity(const FluidState &fluidState, int phaseIdx) | |
422 | { | ||
423 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 59908027 times.
|
65085106 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
424 | |||
425 |
2/2✓ Branch 0 taken 19657388 times.
✓ Branch 1 taken 40250626 times.
|
65085106 | const Scalar T = fluidState.temperature(phaseIdx); |
426 |
2/2✓ Branch 0 taken 19657388 times.
✓ Branch 1 taken 40250626 times.
|
65085106 | const Scalar p = fluidState.pressure(phaseIdx); |
427 | |||
428 |
2/2✓ Branch 0 taken 19657395 times.
✓ Branch 1 taken 40250632 times.
|
65085106 | if (phaseIdx == phase0Idx) |
429 | { | ||
430 | // assume pure water or that each gas molecule displaces exactly one | ||
431 | // molecule in the liquid. | ||
432 | 19891574 | return H2O::liquidMolarDensity(T, p); | |
433 | } | ||
434 | else if (phaseIdx == phase1Idx) | ||
435 | { | ||
436 | if (Policy::useIdealGasDensity()) | ||
437 | // for the gas phase assume an ideal gas | ||
438 | 1577084 | { return IdealGas::molarDensity(T, p); } | |
439 | |||
440 | 44914998 | return H2O::gasMolarDensity(T, fluidState.partialPressure(phase1Idx, H2OIdx)) | |
441 | 45425010 | + Air::gasMolarDensity(T, fluidState.partialPressure(phase1Idx, AirIdx)); | |
442 | } | ||
443 | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | ||
444 | } | ||
445 | |||
446 | using Base<Scalar, ThisType>::viscosity; | ||
447 | /*! | ||
448 | * \brief Calculate the dynamic viscosity of a fluid phase \f$\mathrm{[Pa*s]}\f$ | ||
449 | * | ||
450 | * Compositional effects in the gas phase are accounted by the Wilke method. | ||
451 | * See Reid et al. (1987) \cite reid1987 <BR> | ||
452 | * 4th edition, McGraw-Hill, 1987, 407-410 | ||
453 | * 5th edition, McGraw-Hill, 2001, p. 9.21/22 | ||
454 | * \note Compositional effects for a liquid mixture have to be implemented. | ||
455 | * | ||
456 | * \param fluidState An arbitrary fluid state | ||
457 | * \param phaseIdx The index of the fluid phase to consider | ||
458 | */ | ||
459 | template <class FluidState> | ||
460 | 116896258 | static Scalar viscosity(const FluidState &fluidState, | |
461 | int phaseIdx) | ||
462 | { | ||
463 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 114784965 times.
|
116896258 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
464 | |||
465 |
2/2✓ Branch 0 taken 20752858 times.
✓ Branch 1 taken 94032094 times.
|
116896258 | Scalar T = fluidState.temperature(phaseIdx); |
466 |
2/2✓ Branch 0 taken 20752858 times.
✓ Branch 1 taken 94032094 times.
|
116896258 | Scalar p = fluidState.pressure(phaseIdx); |
467 | |||
468 |
2/2✓ Branch 0 taken 20752865 times.
✓ Branch 1 taken 94032100 times.
|
116896258 | if (phaseIdx == liquidPhaseIdx) |
469 | { | ||
470 | // assume pure water for the liquid phase | ||
471 | 19497306 | return H2O::liquidViscosity(T, p); | |
472 | } | ||
473 | else if (phaseIdx == gasPhaseIdx) | ||
474 | { | ||
475 | if(Policy::useAirViscosityAsGasMixtureViscosity()){ | ||
476 | 394274 | return Air::gasViscosity(T, p); | |
477 | } | ||
478 | else //using a complicated version of this fluid system | ||
479 | { | ||
480 | // Wilke method (Reid et al.): | ||
481 | 39763166 | Scalar muResult = 0; | |
482 | 79526332 | const Scalar mu[numComponents] = { | |
483 | 39763166 | h2oGasViscosityInMixture(T, p), | |
484 | 39763166 | Air::gasViscosity(T, p) | |
485 | }; | ||
486 | |||
487 | // molar masses | ||
488 | 39763166 | const Scalar M[numComponents] = { | |
489 | H2O::molarMass(), | ||
490 | Air::molarMass() | ||
491 | }; | ||
492 | |||
493 |
2/2✓ Branch 0 taken 75772118 times.
✓ Branch 1 taken 37886059 times.
|
119289498 | for (int i = 0; i < numComponents; ++i) |
494 | { | ||
495 | Scalar divisor = 0; | ||
496 | using std::sqrt; | ||
497 |
2/2✓ Branch 0 taken 151544236 times.
✓ Branch 1 taken 75772118 times.
|
238578996 | for (int j = 0; j < numComponents; ++j) |
498 | { | ||
499 | // 1 + (mu[i]/mu[j]^1/2 * (M[i]/M[j])^1/4) | ||
500 | 159052664 | Scalar phiIJ = 1 + sqrt(mu[i]/mu[j] * sqrt(M[j]/M[i])); | |
501 | 159052664 | phiIJ = phiIJ * phiIJ / sqrt(8*(1 + M[i]/M[j])); | |
502 |
2/2✓ Branch 0 taken 2979476 times.
✓ Branch 1 taken 2979476 times.
|
315125828 | divisor += fluidState.moleFraction(phaseIdx, j)*phiIJ; |
503 | } | ||
504 |
2/2✓ Branch 0 taken 1489738 times.
✓ Branch 1 taken 1489738 times.
|
157562914 | muResult += fluidState.moleFraction(phaseIdx, i)*mu[i] / divisor; |
505 | } | ||
506 | return muResult; | ||
507 | } | ||
508 | } | ||
509 | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | ||
510 | } | ||
511 | |||
512 | using Base<Scalar, ThisType>::fugacityCoefficient; | ||
513 | /*! | ||
514 | * \brief Returns the fugacity coefficient \f$\mathrm{[-]}\f$ of a component in a | ||
515 | * phase. | ||
516 | * | ||
517 | * The fugacity coefficient \f$\phi^\kappa_\alpha\f$ of | ||
518 | * component \f$\kappa\f$ in phase \f$\alpha\f$ is connected to | ||
519 | * the fugacity \f$f^\kappa_\alpha\f$ and the component's mole | ||
520 | * fraction \f$x^\kappa_\alpha\f$ by means of the relation | ||
521 | * | ||
522 | * \f[ | ||
523 | f^\kappa_\alpha = \phi^\kappa_\alpha\;x^\kappa_\alpha\;p_\alpha | ||
524 | \f] | ||
525 | * where \f$p_\alpha\f$ is the pressure of the fluid phase. | ||
526 | * | ||
527 | * For liquids with very low miscibility this boils down to the | ||
528 | * Henry constant for the solutes and the saturated vapor pressure | ||
529 | * both divided by phase pressure. | ||
530 | */ | ||
531 | template <class FluidState> | ||
532 | 4257540 | static Scalar fugacityCoefficient(const FluidState &fluidState, | |
533 | int phaseIdx, | ||
534 | int compIdx) | ||
535 | { | ||
536 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4257514 times.
|
4257540 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
537 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4257514 times.
|
4257540 | assert(0 <= compIdx && compIdx < numComponents); |
538 | |||
539 |
2/2✓ Branch 0 taken 2128744 times.
✓ Branch 1 taken 2128744 times.
|
4257540 | Scalar T = fluidState.temperature(phaseIdx); |
540 |
2/2✓ Branch 0 taken 2128744 times.
✓ Branch 1 taken 2128744 times.
|
4257540 | Scalar p = fluidState.pressure(phaseIdx); |
541 | |||
542 |
2/2✓ Branch 0 taken 2128758 times.
✓ Branch 1 taken 2128756 times.
|
4257540 | if (phaseIdx == liquidPhaseIdx) { |
543 |
2/2✓ Branch 0 taken 1064379 times.
✓ Branch 1 taken 1064379 times.
|
2128772 | if (compIdx == H2OIdx) |
544 | 1064386 | return vaporPressure(fluidState, compIdx)/p; | |
545 | 2128772 | return BinaryCoeff::H2O_Air::henry(T)/p; | |
546 | } | ||
547 | |||
548 | // for the gas phase, assume an ideal gas when it comes to | ||
549 | // fugacity (-> fugacity == partial pressure) | ||
550 | return 1.0; | ||
551 | } | ||
552 | |||
553 | /*! | ||
554 | * \brief Returns the relative humidity of the gas phase. | ||
555 | * | ||
556 | * The relative humidity is the ratio of the partial pressure of water vapor | ||
557 | * to the equilibrium vapor pressure of water at a given temperature. | ||
558 | */ | ||
559 | template <class FluidState> | ||
560 | static Scalar relativeHumidity(const FluidState &fluidState) | ||
561 | { | ||
562 | return fluidState.partialPressure(gasPhaseIdx, comp0Idx) | ||
563 | / H2O::vaporPressure(fluidState.temperature(gasPhaseIdx)); | ||
564 | } | ||
565 | |||
566 | using Base<Scalar, ThisType>::diffusionCoefficient; | ||
567 | //! \copydoc Base<Scalar,ThisType>::diffusionCoefficient(const FluidState&,int,int) | ||
568 | template <class FluidState> | ||
569 | 52 | static Scalar diffusionCoefficient(const FluidState &fluidState, | |
570 | int phaseIdx, | ||
571 | int compIdx) | ||
572 | { | ||
573 |
7/16✓ Branch 2 taken 26 times.
✗ Branch 3 not taken.
✓ Branch 11 taken 26 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 26 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 26 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 26 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 26 times.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 26 times.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
|
572 | DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OAir::diffusionCoefficient()"); |
574 | } | ||
575 | |||
576 | using Base<Scalar, ThisType>::binaryDiffusionCoefficient; | ||
577 | //! \copydoc Base<Scalar,ThisType>::binaryDiffusionCoefficient(const FluidState&,int,int,int) | ||
578 | template <class FluidState> | ||
579 | 60298216 | static Scalar binaryDiffusionCoefficient(const FluidState &fluidState, | |
580 | int phaseIdx, | ||
581 | int compIIdx, | ||
582 | int compJIdx) | ||
583 | { | ||
584 | using std::swap; | ||
585 |
2/2✓ Branch 0 taken 37471597 times.
✓ Branch 1 taken 20715287 times.
|
60298216 | if (compIIdx > compJIdx) |
586 | 39348714 | swap(compIIdx, compJIdx); | |
587 | |||
588 |
2/2✓ Branch 0 taken 20320980 times.
✓ Branch 1 taken 37865852 times.
|
60298216 | Scalar T = fluidState.temperature(phaseIdx); |
589 |
2/2✓ Branch 0 taken 20320980 times.
✓ Branch 1 taken 37865852 times.
|
60298216 | Scalar p = fluidState.pressure(phaseIdx); |
590 | |||
591 | // we are in the liquid phase | ||
592 |
2/2✓ Branch 0 taken 20321008 times.
✓ Branch 1 taken 37865876 times.
|
60298216 | if (phaseIdx == liquidPhaseIdx) |
593 | { | ||
594 |
4/4✓ Branch 0 taken 20321001 times.
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 20320994 times.
✓ Branch 3 taken 7 times.
|
20555212 | if (compIIdx == H2OIdx && compJIdx == AirIdx) |
595 | 41110368 | return BinaryCoeff::H2O_Air::liquidDiffCoeff(T, p); | |
596 | else | ||
597 |
10/22✓ Branch 2 taken 14 times.
✗ Branch 3 not taken.
✓ Branch 11 taken 14 times.
✗ Branch 12 not taken.
✓ Branch 16 taken 14 times.
✗ Branch 17 not taken.
✓ Branch 20 taken 14 times.
✗ Branch 21 not taken.
✓ Branch 24 taken 14 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 14 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 14 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 14 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 14 times.
✗ Branch 36 not taken.
✗ Branch 38 not taken.
✓ Branch 39 taken 14 times.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
|
392 | DUNE_THROW(Dune::InvalidStateException, |
598 | "Binary diffusion coefficient of components " | ||
599 | << compIIdx << " and " << compJIdx | ||
600 | << " in phase " << phaseIdx << " is undefined!\n"); | ||
601 | } | ||
602 | |||
603 | // we are in the gas phase | ||
604 |
1/2✓ Branch 0 taken 37865876 times.
✗ Branch 1 not taken.
|
39743004 | else if (phaseIdx == gasPhaseIdx) |
605 | { | ||
606 |
4/4✓ Branch 0 taken 37865870 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 37865864 times.
✓ Branch 3 taken 6 times.
|
39743004 | if (compIIdx == H2OIdx && compJIdx == AirIdx) |
607 | 79485960 | return BinaryCoeff::H2O_Air::gasDiffCoeff(T, p); | |
608 | else | ||
609 |
10/22✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 11 taken 12 times.
✗ Branch 12 not taken.
✓ Branch 16 taken 12 times.
✗ Branch 17 not taken.
✓ Branch 20 taken 12 times.
✗ Branch 21 not taken.
✓ Branch 24 taken 12 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 12 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 12 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 12 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 12 times.
✗ Branch 36 not taken.
✗ Branch 38 not taken.
✓ Branch 39 taken 12 times.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
|
336 | DUNE_THROW(Dune::InvalidStateException, |
610 | "Binary diffusion coefficient of components " | ||
611 | << compIIdx << " and " << compJIdx | ||
612 | << " in phase " << phaseIdx << " is undefined!\n"); | ||
613 | } | ||
614 | |||
615 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | |
616 | } | ||
617 | |||
618 | using Base<Scalar, ThisType>::enthalpy; | ||
619 | /*! | ||
620 | * \brief Given a phase's composition, temperature and pressure, | ||
621 | * return its specific enthalpy \f$\mathrm{[J/kg]}\f$. | ||
622 | * \param fluidState An arbitrary fluid state | ||
623 | * \param phaseIdx The index of the fluid phase to consider | ||
624 | * | ||
625 | * See: | ||
626 | * Class 2001: | ||
627 | * Theorie und numerische Modellierung nichtisothermer Mehrphasenprozesse in NAPL-kontaminierten porösen Medien | ||
628 | * Chapter 2.1.13 Innere Energie, Wäremekapazität, Enthalpie \cite A3:class:2001 <BR> | ||
629 | * | ||
630 | * Formula (2.42): | ||
631 | * the specific enthalpy of a gasphase result from the sum of (enthalpies*mass fraction) of the components | ||
632 | * | ||
633 | * \todo This system neglects the contribution of gas-molecules in the liquid phase. | ||
634 | * This contribution is probably not big. Somebody would have to find out the enthalpy of solution for this system. ... | ||
635 | */ | ||
636 | template <class FluidState> | ||
637 | 32738809 | static Scalar enthalpy(const FluidState &fluidState, | |
638 | int phaseIdx) | ||
639 | { | ||
640 |
2/2✓ Branch 0 taken 6224201 times.
✓ Branch 1 taken 22774130 times.
|
32738809 | const Scalar T = fluidState.temperature(phaseIdx); |
641 |
2/2✓ Branch 0 taken 6224201 times.
✓ Branch 1 taken 22774130 times.
|
32738809 | const Scalar p = fluidState.pressure(phaseIdx); |
642 | |||
643 |
2/2✓ Branch 0 taken 6224208 times.
✓ Branch 1 taken 22774136 times.
|
32738809 | if (phaseIdx == liquidPhaseIdx) |
644 | 6744403 | return H2O::liquidEnthalpy(T, p); | |
645 | |||
646 |
1/2✓ Branch 0 taken 22774136 times.
✗ Branch 1 not taken.
|
25994406 | else if (phaseIdx == gasPhaseIdx) |
647 | 25994406 | return H2O::gasEnthalpy(T, p)*fluidState.massFraction(gasPhaseIdx, H2OIdx) | |
648 | 51988812 | + Air::gasEnthalpy(T, p)*fluidState.massFraction(gasPhaseIdx, AirIdx); | |
649 | |||
650 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | |
651 | } | ||
652 | |||
653 | /*! | ||
654 | * \brief Returns the specific enthalpy \f$\mathrm{[J/kg]}\f$ of a component in a specific phase | ||
655 | * \param fluidState An arbitrary fluid state | ||
656 | * \param phaseIdx The index of the fluid phase to consider | ||
657 | * \param componentIdx The index of the component to consider | ||
658 | * | ||
659 | */ | ||
660 | template <class FluidState> | ||
661 | 61368936 | static Scalar componentEnthalpy(const FluidState &fluidState, | |
662 | int phaseIdx, | ||
663 | int componentIdx) | ||
664 | { | ||
665 |
2/2✓ Branch 0 taken 7884800 times.
✓ Branch 1 taken 50344288 times.
|
61368936 | const Scalar T = fluidState.temperature(phaseIdx); |
666 |
2/2✓ Branch 0 taken 7884800 times.
✓ Branch 1 taken 50344288 times.
|
61368936 | const Scalar p = fluidState.pressure(phaseIdx); |
667 | |||
668 |
2/2✓ Branch 0 taken 7884800 times.
✓ Branch 1 taken 50344288 times.
|
61368936 | if (phaseIdx == liquidPhaseIdx) |
669 | { | ||
670 | // the liquid enthalpy is constant | ||
671 | 7884800 | return H2O::liquidEnthalpy(T, p); | |
672 | } | ||
673 |
1/2✓ Branch 0 taken 50344288 times.
✗ Branch 1 not taken.
|
53484136 | else if (phaseIdx == gasPhaseIdx) |
674 | { | ||
675 |
2/2✓ Branch 0 taken 25172144 times.
✓ Branch 1 taken 25172144 times.
|
53484136 | if (componentIdx == H2OIdx) |
676 | { | ||
677 | 26742068 | return H2O::gasEnthalpy(T, p); | |
678 | } | ||
679 |
1/2✓ Branch 0 taken 25172144 times.
✗ Branch 1 not taken.
|
26742068 | else if (componentIdx == AirIdx) |
680 | { | ||
681 | 53484136 | return Air::gasEnthalpy(T, p); | |
682 | } | ||
683 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << componentIdx); | |
684 | } | ||
685 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | |
686 | } | ||
687 | |||
688 | using Base<Scalar, ThisType>::thermalConductivity; | ||
689 | /*! | ||
690 | * \brief Thermal conductivity of a fluid phase \f$\mathrm{[W/(m K)]}\f$. | ||
691 | * \param fluidState An arbitrary fluid state | ||
692 | * \param phaseIdx The index of the fluid phase to consider | ||
693 | * | ||
694 | * Use the conductivity of air and water as a first approximation. | ||
695 | * Source: | ||
696 | * http://en.wikipedia.org/wiki/List_of_thermal_conductivities | ||
697 | */ | ||
698 | template <class FluidState> | ||
699 | 5878176 | static Scalar thermalConductivity(const FluidState &fluidState, | |
700 | int phaseIdx) | ||
701 | { | ||
702 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5878163 times.
|
5878176 | assert(0 <= phaseIdx && phaseIdx < numPhases); |
703 | |||
704 |
3/6✓ Branch 0 taken 5827341 times.
✓ Branch 1 taken 1063513 times.
✓ Branch 2 taken 1012704 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
57614508 | const Scalar temperature = fluidState.temperature(phaseIdx) ; |
705 |
4/8✓ Branch 0 taken 5827341 times.
✓ Branch 1 taken 1063513 times.
✓ Branch 2 taken 1012704 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1012704 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
108338136 | const Scalar pressure = fluidState.pressure(phaseIdx); |
706 |
2/4✓ Branch 0 taken 5827348 times.
✓ Branch 1 taken 1063519 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
6890880 | if (phaseIdx == liquidPhaseIdx) |
707 | { | ||
708 | 5320999 | return H2O::liquidThermalConductivity(temperature, pressure); | |
709 | } | ||
710 | else if (phaseIdx == gasPhaseIdx) | ||
711 | { | ||
712 | return Air::gasThermalConductivity(temperature, pressure); | ||
713 | } | ||
714 | else | ||
715 | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | ||
716 | } | ||
717 | |||
718 | using Base<Scalar, ThisType>::heatCapacity; | ||
719 | /*! | ||
720 | * \brief Specific isobaric heat capacity of a fluid phase. | ||
721 | * \f$\mathrm{[J/(kg*K)}\f$. | ||
722 | * | ||
723 | * \todo Check whether the gas phase enthalpy is a linear mixture of the component | ||
724 | * enthalpies and the mole fractions is a good assumption. | ||
725 | * | ||
726 | * \param fluidState An arbitrary fluid state | ||
727 | * \param phaseIdx for which phase to give back the heat capacity | ||
728 | */ | ||
729 | template <class FluidState> | ||
730 | 17368803 | static Scalar heatCapacity(const FluidState &fluidState, | |
731 | int phaseIdx) | ||
732 | { | ||
733 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 17368790 times.
|
17368803 | const Scalar temperature = fluidState.temperature(phaseIdx); |
734 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 17368790 times.
|
17368803 | const Scalar pressure = fluidState.pressure(phaseIdx); |
735 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 17368790 times.
|
17368803 | if (phaseIdx == liquidPhaseIdx) |
736 | { | ||
737 | // influence of air is neglected | ||
738 | 5 | return H2O::liquidHeatCapacity(temperature, pressure); | |
739 | } | ||
740 |
1/2✓ Branch 0 taken 17368790 times.
✗ Branch 1 not taken.
|
17368796 | else if (phaseIdx == gasPhaseIdx) |
741 | { | ||
742 | 17368796 | return Air::gasHeatCapacity(temperature, pressure) * fluidState.moleFraction(gasPhaseIdx, AirIdx) | |
743 | 34737586 | + H2O::gasHeatCapacity(temperature, pressure) * fluidState.moleFraction(gasPhaseIdx, H2OIdx); | |
744 | } | ||
745 | else | ||
746 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx); | |
747 | } | ||
748 | }; | ||
749 | |||
750 | } // end namespace FluidSystems | ||
751 | } // end namespace Dumux | ||
752 | |||
753 | #endif | ||
754 |