GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/material/fluidsystems/spe5.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 87 88 98.9%
Functions: 14 14 100.0%
Branches: 195 300 65.0%
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 * \brief @copybrief Dumux::FluidSystems::Spe5
11 */
12 #ifndef DUMUX_SPE5_FLUID_SYSTEM_HH
13 #define DUMUX_SPE5_FLUID_SYSTEM_HH
14
15 #include "spe5parametercache.hh"
16
17 #include <dumux/common/spline.hh>
18 #include <dumux/material/eos/pengrobinsonmixture.hh>
19 #include <dumux/io/name.hh>
20
21 namespace Dumux::FluidSystems {
22
23 /*!
24 * \ingroup FluidSystems
25 * \brief The fluid system for the SPE-5 benchmark problem.
26 *
27 * A three-phase fluid system featuring gas, oil and water as phases and the seven
28 * components distilled water, Methane \f$(\mathrm{C_1})\f$, Propane \f$(\mathrm{C_3})\f$,
29 * Pentane \f$(\mathrm{C_5})\f$, Heptane \f$(\mathrm{C_7})\f$, Decane
30 * \f$(\mathrm{C_{10}})\f$, Pentadecane \f$(\mathrm{C_{15}})\f$ and Icosane
31 * \f$(\mathrm{C_{20}})\f$. For the water phase the IAPWS-97 formulation is used as
32 * equation of state, while for the gas and oil phases a Peng-Robinson
33 * equation of state with slightly modified parameters is used. This fluid system is highly
34 * non-linear, and the gas and oil phases also cannot be considered ideal.
35 *
36 * See:
37 *
38 * J.E. Killough, et al.: Fifth Comparative Solution Project:
39 * Evaluation of Miscible Flood Simulators, Ninth SPE Symposium on
40 * Reservoir Simulation, 1987 \cite SPE5
41 */
42 template <class Scalar>
43 class Spe5
44 {
45 using ThisType = FluidSystems::Spe5<Scalar>;
46
47 using PengRobinsonMixture = Dumux::PengRobinsonMixture<Scalar, ThisType>;
48 using PengRobinson = Dumux::PengRobinson<Scalar>;
49
50 public:
51 using ParameterCache = Spe5ParameterCache<Scalar, ThisType>;
52
53 /****************************************
54 * Fluid phase parameters
55 ****************************************/
56
57 //! Number of phases in the fluid system
58 static const int numPhases = 3;
59
60 //! Index of the gas phase
61 static const int gPhaseIdx = 0;
62 //! Index of the water phase
63 static const int wPhaseIdx = 1;
64 //! Index of the oil phase
65 static const int oPhaseIdx = 2;
66
67 //! The component for pure water to be used
68 using H2O = Dumux::Components::H2O<Scalar>;
69
70 /*!
71 * \brief Return the human readable name of a fluid phase
72 * \param phaseIdx The index of the fluid phase to consider
73 */
74 3 static std::string phaseName(int phaseIdx)
75 {
76
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 3 assert(0 <= phaseIdx && phaseIdx < numPhases);
77
3/3
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 3 switch (phaseIdx)
78 {
79 1 case gPhaseIdx: return IOName::gaseousPhase();
80 1 case wPhaseIdx: return IOName::aqueousPhase();
81 1 case oPhaseIdx: return IOName::naplPhase();
82 }
83 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
84 }
85
86 /*!
87 * \brief Returns whether the fluids are miscible
88 */
89 static constexpr bool isMiscible()
90 { return true; }
91
92 /*!
93 * \brief Return whether a phase is gaseous
94 * \param phaseIdx The index of the fluid phase to consider
95 */
96 static constexpr bool isGas(int phaseIdx)
97 {
98 assert(0 <= phaseIdx && phaseIdx < numPhases);
99 return phaseIdx == gPhaseIdx;
100 }
101
102 /*!
103 * \brief Return whether a phase is compressible
104 * \param phaseIdx The index of the fluid phase to consider
105 *
106 * In the SPE-5 problems all fluids are compressible...
107 */
108 static constexpr bool isCompressible(int phaseIdx)
109 {
110 assert(0 <= phaseIdx && phaseIdx < numPhases);
111 return true;
112 }
113
114 /*!
115 * \brief Returns true if and only if a fluid phase is assumed to
116 * be an ideal mixture.
117 *
118 * We define an ideal mixture as a fluid phase where the fugacity
119 * coefficients of all components times the pressure of the phase
120 * are independent on the fluid composition. This assumption is true
121 * if Henry's law and Raoult's law apply. If you are unsure what
122 * this function should return, it is safe to return false. The
123 * only damage done will be (slightly) increased computation times
124 * in some cases.
125 * \param phaseIdx The index of the fluid phase to consider
126 */
127 static bool isIdealMixture(int phaseIdx)
128 {
129 // always use the reference oil for the fugacity coefficients,
130 // so they cannot be dependent on composition and they the
131 // phases thus always an ideal mixture
132 return phaseIdx == wPhaseIdx;
133 }
134
135 /*!
136 * \brief Returns true if and only if a fluid phase is assumed to
137 * be an ideal gas.
138 *
139 * \param phaseIdx The index of the fluid phase to consider
140 */
141 static bool isIdealGas(int phaseIdx)
142 {
143 assert(0 <= phaseIdx && phaseIdx < numPhases);
144 return false; // gas is not ideal here!
145 }
146
147 /****************************************
148 * Component related parameters
149 ****************************************/
150
151 //! Number of components in the fluid system
152 static const int numComponents = 7;
153
154 static const int H2OIdx = 0;
155 static const int C1Idx = 1;
156 static const int C3Idx = 2;
157 static const int C6Idx = 3;
158 static const int C10Idx = 4;
159 static const int C15Idx = 5;
160 static const int C20Idx = 6;
161
162 /*!
163 * \brief Return the human readable name of a component
164 * \param compIdx The index of the component to consider
165 */
166 7 static std::string componentName(int compIdx)
167 {
168
9/18
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 6 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 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 25 not taken.
✗ Branch 26 not taken.
 8 static std::string name[] = {
169 H2O::name(),
170
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 std::string("C1"),
171 std::string("C3"),
172 std::string("C6"),
173 std::string("C10"),
174 std::string("C15"),
175 std::string("C20")
176 };
177
178
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
 7 assert(0 <= compIdx && compIdx < numComponents);
179 7 return name[compIdx];
180 }
181
182 /*!
183 * \brief Return the molar mass of a component in \f$\mathrm{[kg/mol]}\f$.
184 * \param compIdx The index of the component to consider
185 */
186 506065 static Scalar molarMass(int compIdx)
187 {
188 static const Scalar M[] = {
189 H2O::molarMass(),
190 16.04e-3, // C1
191 44.10e-3, // C3
192 86.18e-3, // C6
193 142.29e-3, // C10
194 206.00e-3, // C15
195 282.00e-3 // C20
196 };
197
198 assert(0 <= compIdx && compIdx < numComponents);
199
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 506072 return M[compIdx];
200 }
201
202 /*!
203 * \brief Critical temperature of a component \f$\mathrm{[K]}\f$.
204 * \param compIdx The index of the component to consider
205 */
206 103215 static Scalar criticalTemperature(int compIdx)
207 {
208 static const Scalar Tcrit[] = {
209 H2O::criticalTemperature(), // H2O
210 343.0*5/9, // C1
211 665.7*5/9, // C3
212 913.4*5/9, // C6
213 1111.8*5/9, // C10
214 1270.0*5/9, // C15
215 1380.0*5/9 // C20
216 };
217
218 assert(0 <= compIdx && compIdx < numComponents);
219
4/4
✓ Branch 0 taken 36845 times.
✓ Branch 1 taken 14738 times.
✓ Branch 2 taken 36880 times.
✓ Branch 3 taken 14752 times.
 103215 return Tcrit[compIdx];
220 }
221
222 /*!
223 * \brief Critical pressure of a component \f$\mathrm{[Pa]}\f$.
224 * \param compIdx The index of the component to consider
225 */
226 103215 static Scalar criticalPressure(int compIdx)
227 {
228 static const Scalar pcrit[] = {
229 H2O::criticalPressure(), // H2O
230 667.8 * 6894.7573, // C1
231 616.3 * 6894.7573, // C3
232 436.9 * 6894.7573, // C6
233 304.0 * 6894.7573, // C10
234 200.0 * 6894.7573, // C15
235 162.0 * 6894.7573 // C20
236 };
237
238 assert(0 <= compIdx && compIdx < numComponents);
239
4/4
✓ Branch 0 taken 36845 times.
✓ Branch 1 taken 14738 times.
✓ Branch 2 taken 36880 times.
✓ Branch 3 taken 14752 times.
 103215 return pcrit[compIdx];
240 }
241
242 /*!
243 * \brief Molar volume of a component at the critical point \f$\mathrm{[m^3/mol]}\f$.
244 * \param compIdx The index of the component to consider
245 */
246 static Scalar criticalMolarVolume(int compIdx)
247 {
248 static const Scalar R = 8.314472;
249 static const Scalar vcrit[] = {
250 H2O::criticalMolarVolume(), // H2O
251 0.290*R*criticalTemperature(1)/criticalPressure(1), // C1
252 0.277*R*criticalTemperature(2)/criticalPressure(2), // C3
253 0.264*R*criticalTemperature(3)/criticalPressure(3), // C6
254 0.257*R*criticalTemperature(4)/criticalPressure(4), // C10
255 0.245*R*criticalTemperature(5)/criticalPressure(5), // C15
256 0.235*R*criticalTemperature(6)/criticalPressure(6) // C20
257 };
258
259 assert(0 <= compIdx && compIdx < numComponents);
260 return vcrit[compIdx];
261 }
262
263 /*!
264 * \brief The acentric factor of a component \f$\mathrm{[-]}\f$.
265 * \param compIdx The index of the component to consider
266 */
267 103215 static Scalar acentricFactor(int compIdx)
268 {
269 static const Scalar accFac[] = {
270 0.344, // H2O (from Reid, et al.)
271 0.0130, // C1
272 0.1524, // C3
273 0.3007, // C6
274 0.4885, // C10
275 0.6500, // C15
276 0.8500 // C20
277 };
278
279 assert(0 <= compIdx && compIdx < numComponents);
280
4/4
✓ Branch 0 taken 36845 times.
✓ Branch 1 taken 14738 times.
✓ Branch 2 taken 36880 times.
✓ Branch 3 taken 14752 times.
 103215 return accFac[compIdx];
281 }
282
283 /*!
284 * \brief Returns the interaction coefficient for two components.
285 * \param comp1Idx The index of component 1 to consider
286 * \param comp2Idx The index of component 2 to consider
287 *
288 * The values are given by the SPE5 paper.
289 */
290
6/6
✓ Branch 0 taken 154749 times.
✓ Branch 1 taken 206332 times.
✓ Branch 2 taken 795858 times.
✓ Branch 3 taken 1061144 times.
✓ Branch 4 taken 154686 times.
✓ Branch 5 taken 206248 times.
 2579017 static Scalar interactionCoefficient(int comp1Idx, int comp2Idx)
291 {
292 using std::min;
293 using std::max;
294
6/6
✓ Branch 0 taken 154749 times.
✓ Branch 1 taken 206332 times.
✓ Branch 2 taken 795858 times.
✓ Branch 3 taken 1061144 times.
✓ Branch 4 taken 154686 times.
✓ Branch 5 taken 206248 times.
 2579017 int i = min(comp1Idx, comp2Idx);
295 2579017 int j = max(comp1Idx, comp2Idx);
296
12/12
✓ Branch 0 taken 81059 times.
✓ Branch 1 taken 280022 times.
✓ Branch 2 taken 51583 times.
✓ Branch 3 taken 29476 times.
✓ Branch 4 taken 416878 times.
✓ Branch 5 taken 1440124 times.
✓ Branch 6 taken 265286 times.
✓ Branch 7 taken 151592 times.
✓ Branch 8 taken 81026 times.
✓ Branch 9 taken 279908 times.
✓ Branch 10 taken 51562 times.
✓ Branch 11 taken 29464 times.
 2579017 if (i == C1Idx && (j == C15Idx || j == C20Idx))
297 return 0.05;
298
12/12
✓ Branch 0 taken 265284 times.
✓ Branch 1 taken 66321 times.
✓ Branch 2 taken 36845 times.
✓ Branch 3 taken 29476 times.
✓ Branch 4 taken 1364328 times.
✓ Branch 5 taken 341082 times.
✓ Branch 6 taken 189490 times.
✓ Branch 7 taken 151592 times.
✓ Branch 8 taken 265176 times.
✓ Branch 9 taken 66294 times.
✓ Branch 10 taken 36830 times.
✓ Branch 11 taken 29464 times.
 2368485 else if (i == C3Idx && (j == C15Idx || j == C20Idx))
299 return 0.005;
300 return 0;
301 }
302
303 /****************************************
304 * Methods which compute stuff
305 ****************************************/
306
307 /*!
308 * \brief Initialize the fluid system's static parameters
309 */
310 2 static void init()
311 {
312 PengRobinsonParamsMixture<Scalar, ThisType, gPhaseIdx, /*useSpe5=*/true> prParams;
313
314 // find envelopes of the 'a' and 'b' parameters for the range
315 // 273.15K <= T <= 373.15K and 10e3 Pa <= p <= 100e6 Pa. For
316 // this we take advantage of the fact that 'a' and 'b' for
317 // mixtures is just a convex combination of the attractive and
318 // repulsive parameters of the pure components
319 2 Scalar minT = 273.15, maxT = 373.15;
320 2 Scalar minP = 1e4, maxP = 100e6;
321
322 2 Scalar minA = 1e100, maxA = -1e100;
323 2 Scalar minB = 1e100, maxB = -1e100;
324
325 2 prParams.updatePure(minT, minP);
326 using std::min;
327 using std::max;
328
2/2
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 2 times.
 16 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
329
4/4
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 12 times.
 24 minA = min(prParams.pureParams(compIdx).a(), minA);
330
4/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 2 times.
 16 maxA = max(prParams.pureParams(compIdx).a(), maxA);
331
3/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 14 times.
 26 minB = min(prParams.pureParams(compIdx).b(), minB);
332
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14 times.
 14 maxB = max(prParams.pureParams(compIdx).b(), maxB);
333 }
334
335 2 prParams.updatePure(maxT, minP);
336
2/2
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 2 times.
 16 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
337
3/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 14 times.
✗ Branch 3 not taken.
 26 minA = min(prParams.pureParams(compIdx).a(), minA);
338
3/4
✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 2 times.
 28 maxA = max(prParams.pureParams(compIdx).a(), maxA);
339
4/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 2 times.
 26 minB = min(prParams.pureParams(compIdx).b(), minB);
340
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2 times.
 26 maxB = max(prParams.pureParams(compIdx).b(), maxB);
341 }
342
343 2 prParams.updatePure(minT, maxP);
344
2/2
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 2 times.
 16 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
345
3/4
✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 2 times.
 28 minA = min(prParams.pureParams(compIdx).a(), minA);
346
4/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 2 times.
 26 maxA = max(prParams.pureParams(compIdx).a(), maxA);
347
4/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 2 times.
 26 minB = min(prParams.pureParams(compIdx).b(), minB);
348
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2 times.
 26 maxB = max(prParams.pureParams(compIdx).b(), maxB);
349 }
350
351 2 prParams.updatePure(maxT, maxP);
352
2/2
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 2 times.
 16 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
353
3/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 14 times.
✗ Branch 3 not taken.
 26 minA = min(prParams.pureParams(compIdx).a(), minA);
354
3/4
✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 2 times.
 28 maxA = max(prParams.pureParams(compIdx).a(), maxA);
355
4/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 2 times.
 26 minB = min(prParams.pureParams(compIdx).b(), minB);
356
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 2 times.
 26 maxB = max(prParams.pureParams(compIdx).b(), maxB);
357 }
358
359 2 PengRobinson::init(/*aMin=*/minA, /*aMax=*/maxA, /*na=*/100,
360 /*bMin=*/minB, /*bMax=*/maxB, /*nb=*/200);
361 2 }
362
363 /*!
364 * \brief Calculate the density \f$\mathrm{[kg/m^3]}\f$ of a fluid phase
365 *
366 *
367 * \param fluidState An arbitrary fluid state
368 * \param paramCache Container for cache parameters
369 * \param phaseIdx The index of the fluid phase to consider
370 */
371 template <class FluidState>
372 56830 static Scalar density(const FluidState &fluidState,
373 const ParameterCache &paramCache,
374 int phaseIdx)
375 {
376
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 56830 times.
 56830 assert(0 <= phaseIdx && phaseIdx < numPhases);
377
378
2/3
✓ Branch 1 taken 56827 times.
✓ Branch 2 taken 3 times.
✗ Branch 0 not taken.
 56830 return fluidState.averageMolarMass(phaseIdx)/paramCache.molarVolume(phaseIdx);
379 }
380
381 /*!
382 * \brief The molar density \f$\rho_{mol,\alpha}\f$
383 * of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
384 *
385 * The molar density is defined by the
386 * mass density \f$\rho_\alpha\f$ and the mean molar mass \f$\overline M_\alpha\f$:
387 *
388 * \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{\overline M_\alpha} \;.\f]
389 */
390 template <class FluidState>
391
5/10
✗ Branch 0 not taken.
✓ Branch 1 taken 4989 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9972 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 34902 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 6966 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
 56830 static Scalar molarDensity(const FluidState &fluidState, const ParameterCache &paramCache, int phaseIdx)
392 {
393
2/2
✓ Branch 0 taken 23268 times.
✓ Branch 1 taken 11634 times.
 56830 return 1.0/paramCache.molarVolume(phaseIdx);
394 }
395
396 /*!
397 * \brief Calculate the dynamic viscosity of a fluid phase \f$\mathrm{[Pa*s]}\f$
398 * \param fs An arbitrary fluid state
399 * \param paramCache Container for cache parameters
400 * \param phaseIdx The index of the fluid phase to consider
401 */
402 template <class FluidState>
403 3 static Scalar viscosity(const FluidState &fs,
404 const ParameterCache &paramCache,
405 int phaseIdx)
406 {
407
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 3 assert(0 <= phaseIdx && phaseIdx <= numPhases);
408
409
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 if (phaseIdx == gPhaseIdx) {
410 // given by SPE-5 in table on page 64. we use a constant
411 // viscosity, though...
412 return 0.0170e-2 * 0.1;
413 }
414
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 else if (phaseIdx == wPhaseIdx)
415 // given by SPE-5: 0.7 centi-Poise = 0.0007 Pa s
416 return 0.7e-2 * 0.1;
417 else {
418
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 assert(phaseIdx == oPhaseIdx);
419 // given by SPE-5 in table on page 64. we use a constant
420 // viscosity, though...
421 return 0.208e-2 * 0.1;
422 }
423 }
424
425 /*!
426 * \brief Calculate the fugacity coefficient \f$\mathrm{[-]}\f$ of an individual
427 * component in a fluid phase
428 *
429 * The fugacity coefficient \f$\mathrm{\phi^\kappa_\alpha}\f$ is connected
430 * to the fugacity \f$\mathrm{f^\kappa_\alpha}\f$ and the component's mole
431 * fraction in a phase \f$\mathrm{x^\kappa_\alpha}\f$ by means of the
432 * relation
433 *
434 * \f[ f^\kappa_\alpha = \phi^\kappa_\alpha \cdot x^\kappa_\alpha \cdot p_alpha \f]
435 * \param fs An arbitrary fluid state
436 * \param paramCache Container for cache parameters
437 * \param phaseIdx The index of the fluid phase to consider
438 * \param compIdx The index of the component to consider
439 */
440 template <class FluidState>
441 316386 static Scalar fugacityCoefficient(const FluidState &fs,
442 const ParameterCache &paramCache,
443 int phaseIdx,
444 int compIdx)
445 {
446
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 316386 times.
 316386 assert(0 <= phaseIdx && phaseIdx <= numPhases);
447
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 316386 times.
 316386 assert(0 <= compIdx && compIdx <= numComponents);
448
449
2/2
✓ Branch 0 taken 265216 times.
✓ Branch 1 taken 51170 times.
 316386 if (phaseIdx == oPhaseIdx || phaseIdx == gPhaseIdx)
450 265216 return PengRobinsonMixture::computeFugacityCoefficient(fs,
451 paramCache,
452 phaseIdx,
453 265216 compIdx);
454 else {
455
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 51170 times.
 51170 assert(phaseIdx == wPhaseIdx);
456 51170 return henryCoeffWater_(compIdx, fs.temperature(wPhaseIdx))
457 51170 / fs.pressure(wPhaseIdx);
458 }
459 }
460
461
462 /*!
463 * \brief Calculate the binary molecular diffusion coefficient for
464 * a component in a fluid phase \f$\mathrm{[mol^2 * s / (kg*m^3)]}\f$
465 *
466 * Molecular diffusion of a component \f$\mathrm{\kappa}\f$ is caused by a
467 * gradient of the chemical potential and follows the law
468 *
469 * \f[ J = - D grad \mu_\kappa \f]
470 *
471 * where \f$\mathrm{\mu_\kappa}\f$ is the component's chemical potential,
472 * \f$\mathrm{D}\f$ is the diffusion coefficient and \f$\mathrm{J}\f$ is the
473 * diffusive flux. \f$\mathrm{\mu_\kappa}\f$ is connected to the component's
474 * fugacity \f$\mathrm{f_\kappa}\f$ by the relation
475 *
476 * \f[ \mu_\kappa = R T_\alpha \mathrm{ln} \frac{f_\kappa}{p_\alpha} \f]
477 *
478 * where \f$\mathrm{p_\alpha}\f$ and \f$\mathrm{T_\alpha}\f$ are the fluid phase'
479 * pressure and temperature.
480 * \param fs An arbitrary fluid state
481 * \param paramCache Container for cache parameters
482 * \param phaseIdx The index of the fluid phase to consider
483 * \param compIdx The index of the component to consider
484 */
485 template <class FluidState>
486 21 static Scalar diffusionCoefficient(const FluidState &fs,
487 const ParameterCache &paramCache,
488 int phaseIdx,
489 int compIdx)
490
10/20
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 21 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 21 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 21 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 21 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 21 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 21 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 21 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 21 times.
✗ Branch 29 not taken.
 84 { DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients"); }
491
492 /*!
493 * \brief Given a phase's composition, temperature and pressure,
494 * return the binary diffusion coefficient \f$\mathrm{[m^2/s]}\f$ for components
495 * \f$i\f$ and \f$j\f$ in this phase.
496 *
497 * \param fluidState An arbitrary fluid state
498 * \param paramCache Container for cache parameters
499 * \param phaseIdx The index of the fluid phase to consider
500 * \param compIIdx The index of the first component to consider
501 * \param compJIdx The index of the second component to consider
502 */
503 template <class FluidState>
504 147 static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
505 const ParameterCache &paramCache,
506 int phaseIdx,
507 int compIIdx,
508 int compJIdx)
509
10/20
✓ Branch 1 taken 147 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 147 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 147 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 147 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 147 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 147 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 147 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 147 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 147 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 147 times.
✗ Branch 29 not taken.
 588 { DUNE_THROW(Dune::NotImplemented, "Binary diffusion coefficients"); }
510
511 /*!
512 * \brief Given a phase's composition, temperature and pressure,
513 * calculate its specific enthalpy \f$\mathrm{[J/kg]}\f$.
514 * \param fs An arbitrary fluid state
515 * \param paramCache Container for cache parameters
516 * \param phaseIdx The index of the fluid phase to consider
517 */
518 template <class FluidState>
519 3 static Scalar enthalpy(const FluidState &fs,
520 const ParameterCache &paramCache,
521 int phaseIdx)
522
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, "Enthalpies"); }
523
524 /*!
525 * \brief Given a phase's composition, temperature and pressure,
526 * calculate its thermal conductivity \f$\mathrm{[W/(m K)]}\f$.
527 * \param fluidState An arbitrary fluid state
528 * \param paramCache Container for cache parameters
529 * \param phaseIdx The index of the fluid phase to consider
530 */
531 template <class FluidState>
532 3 static Scalar thermalConductivity(const FluidState &fluidState,
533 const ParameterCache &paramCache,
534 int phaseIdx)
535
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, "Thermal conductivities"); }
536
537 /*!
538 * \brief Given a phase's composition, temperature and pressure,
539 * calculate its heat capacity \f$\mathrm{[J/(kg K)]}\f$.
540 * \param fluidState An arbitrary fluid state
541 * \param paramCache Container for cache parameters
542 * \param phaseIdx The index of the fluid phase to consider
543 */
544 template <class FluidState>
545 3 static Scalar heatCapacity(const FluidState &fluidState,
546 const ParameterCache &paramCache,
547 int phaseIdx)
548
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, "Heat capacities"); }
549
550
551 private:
552 51170 static Scalar henryCoeffWater_(int compIdx, Scalar temperature)
553 {
554 // use henry's law for the solutes and the vapor pressure for
555 // the solvent.
556
3/4
✓ Branch 0 taken 7310 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 36550 times.
✓ Branch 3 taken 7310 times.
 51170 switch (compIdx) {
557 7310 case H2OIdx: return H2O::vaporPressure(temperature);
558
559 // the values of the Henry constant for the solutes have
560 // are faked so far...
561 case C1Idx: return 5.57601e+09;
562 case C3Idx: return 1e10;
563 case C6Idx: return 1e10;
564 case C10Idx: return 1e10;
565 case C15Idx: return 1e10;
566 case C20Idx: return 1e10;
567 ✗ default: DUNE_THROW(Dune::InvalidStateException, "Unknown component index " << compIdx);
568 }
569 }
570 };
571
572 } // end namespace Dumux::FluidSystems
573
574 #endif
575