GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/material/fluidstates/nonequilibrium.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 56 77 72.7%
Functions: 9 25 36.0%
Branches: 79 114 69.3%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup FluidStates
10 * \brief Represents all relevant thermodynamic quantities of a
11 * multi-phase, multi-component fluid system without using
12 * any assumptions.
13 */
14 #ifndef DUMUX_NONEQUILIBRIUM_FLUID_STATE_HH
15 #define DUMUX_NONEQUILIBRIUM_FLUID_STATE_HH
16
17 #include <cmath>
18 #include <algorithm>
19 #include <dune/common/exceptions.hh>
20
21 namespace Dumux {
22
23 /*!
24 * \ingroup FluidStates
25 * \brief Represents all relevant thermodynamic quantities of a
26 * multi-phase, multi-component fluid system without using
27 * any assumptions.
28 */
29 template <class ScalarType, class FluidSystem>
30 535095 class NonEquilibriumFluidState
31 {
32 public:
33 static constexpr int numPhases = FluidSystem::numPhases;
34 static constexpr int numComponents = FluidSystem::numComponents;
35
36 //! export the scalar type
37 using Scalar = ScalarType;
38
39 /*****************************************************
40 * Generic access to fluid properties (No assumptions
41 * on thermodynamic equilibrium required)
42 *****************************************************/
43 /*!
44 * \brief Returns the index of the wetting phase in the
45 * fluid-solid configuration (for porous medium systems).
46 */
47 ✗ int wettingPhase() const { return wPhaseIdx_; }
48 /*!
49 * \brief Returns the saturation \f$S_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[-]}\f$.
50 *
51 * The saturation is defined as the pore space occupied by the fluid divided by the total pore space:
52 * \f[S_\alpha := \frac{\phi \mathcal{V}_\alpha}{\phi \mathcal{V}}\f]
53 *
54 * \param phaseIdx the index of the phase
55 */
56 Scalar saturation(int phaseIdx) const
57
5/14
✗ Branch 0 not taken.
✓ Branch 1 taken 1627501 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 4863180 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 397768 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 795536 times.
 48740666 { return saturation_[phaseIdx]; }
58
59 /*!
60 * \brief Returns the molar fraction \f$x^\kappa_\alpha\f$ of the component \f$\kappa\f$ in fluid phase \f$\alpha\f$ in \f$\mathrm{[-]}\f$.
61 *
62 * The molar fraction \f$x^\kappa_\alpha\f$ is defined as the ratio of the number of molecules
63 * of component \f$\kappa\f$ and the total number of molecules of the phase \f$\alpha\f$.
64 *
65 * \param phaseIdx the index of the phase
66 * \param compIdx the index of the component
67 */
68 Scalar moleFraction(int phaseIdx, int compIdx) const
69 91704541 { return moleFraction_[phaseIdx][compIdx]; }
70
71
72 /*!
73 * \brief Returns the mass fraction \f$X^\kappa_\alpha\f$ of component \f$\kappa\f$ in fluid phase \f$\alpha\f$ in \f$\mathrm{[-]}\f$.
74 *
75 * The mass fraction \f$X^\kappa_\alpha\f$ is defined as the weight of all molecules of a
76 * component divided by the total mass of the fluid phase. It is related with the component's mole fraction by means of the relation
77 *
78 * \f[X^\kappa_\alpha = x^\kappa_\alpha \frac{M^\kappa}{\overline M_\alpha}\;,\f]
79 *
80 * where \f$M^\kappa\f$ is the molar mass of component \f$\kappa\f$ and
81 * \f$\overline M_\alpha\f$ is the mean molar mass of a molecule of phase
82 * \f$\alpha\f$.
83 *
84 * \param phaseIdx the index of the phase
85 * \param compIdx the index of the component
86 */
87 70696174 Scalar massFraction(int phaseIdx, int compIdx) const
88 {
89 using std::abs;
90 using std::max;
91
2/2
✓ Branch 0 taken 21405440 times.
✓ Branch 1 taken 42480598 times.
 70696174 return abs(sumMoleFractions_[phaseIdx])
92 70696174 * moleFraction_[phaseIdx][compIdx]
93
2/2
✓ Branch 0 taken 21405440 times.
✓ Branch 1 taken 42480598 times.
 70696174 * FluidSystem::molarMass(compIdx)
94
4/4
✓ Branch 0 taken 70696172 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 70696172 times.
✓ Branch 3 taken 2 times.
 141392348 / max(1e-40, abs(averageMolarMass_[phaseIdx]));
95 }
96
97 /*!
98 * \brief The average molar mass \f$\overline M_\alpha\f$ of phase \f$\alpha\f$ in \f$\mathrm{[kg/mol]}\f$
99 *
100 * The average molar mass is the mean mass of a mole of the
101 * fluid at current composition. It is defined as the sum of the
102 * component's molar masses weighted by the current mole fraction:
103 * \f[\mathrm{ \overline M_\alpha = \sum_\kappa M^\kappa x_\alpha^\kappa}\f]
104 */
105 Scalar averageMolarMass(int phaseIdx) const
106
3/5
✓ Branch 0 taken 89466 times.
✓ Branch 1 taken 121213 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 3041968 { return averageMolarMass_[phaseIdx]; }
107
108 /*!
109 * \brief The molar concentration \f$c^\kappa_\alpha\f$ of component \f$\kappa\f$ in fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
110 *
111 * This quantity is usually called "molar concentration" or just
112 * "concentration", but there are many other (though less common)
113 * measures for concentration.
114 *
115 * http://en.wikipedia.org/wiki/Concentration
116 */
117 Scalar molarity(int phaseIdx, int compIdx) const
118 4 { return molarDensity(phaseIdx)*moleFraction(phaseIdx, compIdx); }
119
120 /*!
121 * \brief The fugacity coefficient \f$\Phi^\kappa_\alpha\f$ of component \f$\kappa\f$ in fluid phase \f$\alpha\f$ in \f$\mathrm{[-]}\f$
122 */
123 Scalar fugacityCoefficient(int phaseIdx, int compIdx) const
124 5665768 { return fugacityCoefficient_[phaseIdx][compIdx]; }
125
126 /*!
127 * \brief The fugacity \f$f^\kappa_\alpha\f$ of component \f$\kappa\f$
128 * in fluid phase \f$\alpha\f$ in \f$\mathrm{[Pa]}\f$
129 *
130 * The fugacity is defined as:
131 * \f$f_\alpha^\kappa := \Phi^\kappa_\alpha x^\kappa_\alpha p_\alpha \;,\f$
132 * where \f$\Phi^\kappa_\alpha\f$ is the fugacity coefficient \cite reid1987 .
133 * The physical meaning of fugacity becomes clear from the equation:
134 * \f[f_\alpha^\kappa = p_\alpha \exp\left\{\frac{\zeta^\kappa_\alpha}{R T_\alpha} \right\} \;,\f]
135 * where \f$\zeta^\kappa_\alpha\f$ represents the \f$\kappa\f$'s chemical
136 * potential in phase \f$\alpha\f$, \f$R\f$ stands for the ideal gas constant,
137 * and \f$T_\alpha\f$ for the absolute temperature of phase \f$\alpha\f$. Assuming thermal equilibrium,
138 * there is a one-to-one mapping between a component's chemical potential
139 * \f$\zeta^\kappa_\alpha\f$ and its fugacity \f$f^\kappa_\alpha\f$. In this
140 * case chemical equilibrium can thus be expressed by:
141 * \f[f^\kappa := f^\kappa_\alpha = f^\kappa_\beta\quad\forall \alpha, \beta\f]
142 */
143 Scalar fugacity(int phaseIdx, int compIdx) const
144 {
145 40320 return pressure_[phaseIdx]
146
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 40322 *fugacityCoefficient_[phaseIdx][compIdx]
147 40320 *moleFraction_[phaseIdx][compIdx];
148 }
149
150 Scalar fugacity(int compIdx) const
151 80640 { return fugacity(0, compIdx); }
152
153 /*!
154 * \brief The molar volume \f$v_{mol,\alpha}\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[m^3/mol]}\f$
155 *
156 * This quantity is the inverse of the molar density.
157 */
158 Scalar molarVolume(int phaseIdx) const
159
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 { return 1.0/molarDensity(phaseIdx); }
160
161 /*!
162 * \brief The mass density \f$\rho_\alpha\f$ of the fluid phase
163 * \f$\alpha\f$ in \f$\mathrm{[kg/m^3]}\f$
164 */
165 ✗ Scalar density(int phaseIdx) const
166
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 92748658 { return density_[phaseIdx]; }
167
168 /*!
169 * \brief The molar density \f$\rho_{mol,\alpha}\f$
170 * of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
171 *
172 * The molar density is defined by the mass density \f$\rho_\alpha\f$ and the mean molar mass \f$\overline M_\alpha\f$:
173 *
174 * \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{\overline M_\alpha} \;.\f]
175 */
176 ✗ Scalar molarDensity(int phaseIdx) const
177
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 57161022 { return molarDensity_[phaseIdx]; }
178
179 /*!
180 * \brief The absolute temperature\f$T_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[K]}\f$
181 */
182 ✗ Scalar temperature(const int phaseIdx) const
183
22/23
✓ Branch 0 taken 5161380 times.
✓ Branch 1 taken 4863180 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 275268 times.
✓ Branch 4 taken 272468 times.
✓ Branch 5 taken 1414846 times.
✓ Branch 6 taken 1414846 times.
✓ Branch 7 taken 1476636 times.
✓ Branch 8 taken 1476636 times.
✓ Branch 9 taken 4832576 times.
✓ Branch 10 taken 4832576 times.
✓ Branch 11 taken 5377512 times.
✓ Branch 12 taken 5377512 times.
✓ Branch 13 taken 501200 times.
✓ Branch 14 taken 2326274 times.
✓ Branch 15 taken 2308220 times.
✓ Branch 16 taken 2097542 times.
✓ Branch 17 taken 2161052 times.
✓ Branch 18 taken 3775448 times.
✓ Branch 19 taken 4783084 times.
✓ Branch 20 taken 1147636 times.
✓ Branch 21 taken 292600 times.
✓ Branch 22 taken 295400 times.
 85008130 { return temperature_[phaseIdx]; }
184
185 /*!
186 * \brief Get the equilibrium temperature \f$\mathrm{[K]}\f$ of the fluid phases.
187 */
188 Scalar temperature() const
189 { return temperatureEquil_ ; }
190
191 /*!
192 * \brief The pressure \f$p_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[Pa]}\f$
193 */
194 ✗ Scalar pressure(int phaseIdx) const
195
18/20
✓ Branch 0 taken 4848480 times.
✓ Branch 1 taken 4865980 times.
✓ Branch 2 taken 397768 times.
✓ Branch 3 taken 1414846 times.
✓ Branch 4 taken 1414846 times.
✓ Branch 5 taken 1204168 times.
✓ Branch 6 taken 1601936 times.
✓ Branch 7 taken 4832576 times.
✓ Branch 8 taken 4832576 times.
✓ Branch 9 taken 4832576 times.
✓ Branch 10 taken 5628112 times.
✓ Branch 11 taken 1591072 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 2222842 times.
✓ Branch 15 taken 2433520 times.
✓ Branch 16 taken 3250610 times.
✓ Branch 17 taken 3478788 times.
✓ Branch 18 taken 3650148 times.
✓ Branch 19 taken 3650148 times.
 91412958 { return pressure_[phaseIdx]; }
196
197 /*!
198 * \brief The partial pressure of a component in a phase \f$\mathrm{[Pa]}\f$
199 * \todo is this needed?
200 */
201 Scalar partialPressure(int phaseIdx, int compIdx) const
202 {
203 assert(FluidSystem::isGas(phaseIdx));
204 return moleFraction(phaseIdx, compIdx) * pressure(phaseIdx);
205 }
206
207 /*!
208 * \brief The specific enthalpy \f$h_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[J/kg]}\f$
209 */
210 ✗ Scalar enthalpy(int phaseIdx) const
211
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 18192332 { return enthalpy_[phaseIdx]; }
212
213 /*!
214 * \brief The specific internal energy \f$u_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[J/kg]}\f$
215 *
216 * The specific internal energy is defined by the relation:
217 *
218 * \f[u_\alpha = h_\alpha - \frac{p_\alpha}{\rho_\alpha}\f]
219 */
220 Scalar internalEnergy(int phaseIdx) const
221
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 17227362 { return enthalpy_[phaseIdx] - pressure(phaseIdx)/density(phaseIdx); }
222
223 /*!
224 * \brief The dynamic viscosity \f$\mu_\alpha\f$ of fluid phase \f$\alpha\f$ in \f$\mathrm{[Pa s]}\f$
225 */
226 ✗ Scalar viscosity(int phaseIdx) const
227
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 42869738 { return viscosity_[phaseIdx]; }
228
229 /*****************************************************
230 * Setter methods. Note that these are not part of the
231 * generic FluidState interface but specific for each
232 * implementation...
233 *****************************************************/
234 /*!
235 * \brief Set the temperature \f$\mathrm{[K]}\f$ of a fluid phase
236 */
237 ✗ void setTemperature(int phaseIdx, Scalar value)
238 4224707 { temperature_[phaseIdx] = value; }
239
240 /*!
241 * \brief Set the temperature \f$\mathrm{[K]}\f$ of all fluid phases.
242 */
243 void setTemperature(Scalar value)
244 { temperatureEquil_ = value; }
245
246 /*!
247 * \brief Set the fluid pressure of a phase \f$\mathrm{[Pa]}\f$
248 */
249 ✗ void setPressure(int phaseIdx, Scalar value)
250 4224707 { pressure_[phaseIdx] = value; }
251
252 /*!
253 * \brief Set the saturation of a phase \f$\mathrm{[-]}\f$
254 */
255 ✗ void setSaturation(int phaseIdx, Scalar value)
256 3803351 { saturation_[phaseIdx] = value; }
257
258 /*!
259 * \brief Set the mole fraction of a component in a phase \f$\mathrm{[-]}\f$
260 */
261 13690634 void setMoleFraction(int phaseIdx, int compIdx, Scalar value)
262 {
263 using std::isfinite;
264
2/4
✓ Branch 0 taken 13690634 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 13690634 times.
✗ Branch 3 not taken.
 27381268 if (isfinite(averageMolarMass_[phaseIdx]))
265 {
266 13690634 Scalar delta = value - moleFraction_[phaseIdx][compIdx];
267
268 13690634 moleFraction_[phaseIdx][compIdx] = value;
269
270 13690634 sumMoleFractions_[phaseIdx] += delta;
271
2/2
✓ Branch 0 taken 1391823 times.
✓ Branch 1 taken 9116667 times.
 16474280 averageMolarMass_[phaseIdx] += delta*FluidSystem::molarMass(compIdx);
272 }
273 else
274 {
275 ✗ moleFraction_[phaseIdx][compIdx] = value;
276
277 // re-calculate the mean molar mass
278 ✗ sumMoleFractions_[phaseIdx] = 0.0;
279 ✗ averageMolarMass_[phaseIdx] = 0.0;
280 ✗ for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
281 ✗ sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compJIdx];
282 ✗ averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compJIdx]*FluidSystem::molarMass(compJIdx);
283 }
284 }
285 13690634 }
286
287 /*!
288 * \brief Set the mass fraction of a component in a phase \f$\mathrm{[-]}\f$
289 * and update the average molar mass \f$\mathrm{[kg/mol]}\f$ according
290 * to the current composition of the phase
291 */
292 void setMassFraction(int phaseIdx, int compIdx, Scalar value)
293 {
294 if (numComponents != 2)
295 DUNE_THROW(Dune::NotImplemented, "This currently only works for 2 components.");
296 else
297 {
298 // calculate average molar mass of the gas phase
299 Scalar M1 = FluidSystem::molarMass(compIdx);
300 Scalar M2 = FluidSystem::molarMass(1-compIdx);
301 Scalar X2 = 1.0-value;
302 Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
303
304 moleFraction_[phaseIdx][compIdx] = value * avgMolarMass / M1;
305 moleFraction_[phaseIdx][1-compIdx] = 1.0-moleFraction_[phaseIdx][compIdx];
306
307 // re-calculate the mean molar mass
308 sumMoleFractions_[phaseIdx] = 0.0;
309 averageMolarMass_[phaseIdx] = 0.0;
310 for (int compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
311 sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compJIdx];
312 averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compJIdx]*FluidSystem::molarMass(compJIdx);
313 }
314 }
315 }
316
317 /*!
318 * \brief Set the fugacity of a component in a phase \f$\mathrm{[-]}\f$
319 */
320 void setFugacityCoefficient(int phaseIdx, int compIdx, Scalar value)
321 10482440 { fugacityCoefficient_[phaseIdx][compIdx] = value; }
322
323 /*!
324 * \brief Set the density of a phase \f$\mathrm{[kg / m^3]}\f$
325 */
326 ✗ void setDensity(int phaseIdx, Scalar value)
327 6843721 { density_[phaseIdx] = value; }
328
329 /*!
330 * \brief Set the molar density of a phase \f$\mathrm{[kg / m^3]}\f$
331 */
332 ✗ void setMolarDensity(int phaseIdx, Scalar value)
333 6843721 { molarDensity_[phaseIdx] = value; }
334
335 /*!
336 * \brief Set the specific enthalpy of a phase \f$\mathrm{[J/m^3]}\f$
337 */
338 ✗ void setEnthalpy(int phaseIdx, Scalar value)
339 4221515 { enthalpy_[phaseIdx] = value; }
340
341 /*!
342 * \brief Set the dynamic viscosity of a phase \f$\mathrm{[Pa s]}\f$
343 */
344 ✗ void setViscosity(int phaseIdx, Scalar value)
345 3800159 { viscosity_[phaseIdx] = value; }
346
347 /*!
348 * \brief Set the index of the wetting phase
349 */
350 ✗ void setWettingPhase(int phaseIdx)
351
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 397768 times.
 1204168 { wPhaseIdx_ = phaseIdx; }
352
353 /*!
354 * \brief Retrieve all parameters from an arbitrary fluid
355 * state.
356 * \param fs Fluidstate
357 */
358 template <class FluidState>
359 1204168 void assign(const FluidState& fs)
360 {
361
2/2
✓ Branch 0 taken 2408336 times.
✓ Branch 1 taken 1204168 times.
 3612504 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
362 {
363 2408336 averageMolarMass_[phaseIdx] = 0;
364 2408336 sumMoleFractions_[phaseIdx] = 0;
365
2/2
✓ Branch 0 taken 4816672 times.
✓ Branch 1 taken 2408336 times.
 7225008 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
366 {
367 9633344 moleFraction_[phaseIdx][compIdx] = fs.moleFraction(phaseIdx, compIdx);
368 9633344 fugacityCoefficient_[phaseIdx][compIdx] = fs.fugacityCoefficient(phaseIdx, compIdx);
369 4816672 averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compIdx]*FluidSystem::molarMass(compIdx);
370 4816672 sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compIdx];
371 }
372 4816672 pressure_[phaseIdx] = fs.pressure(phaseIdx);
373 4816672 saturation_[phaseIdx] = fs.saturation(phaseIdx);
374 4816672 density_[phaseIdx] = fs.density(phaseIdx);
375 4816672 molarDensity_[phaseIdx] = fs.molarDensity(phaseIdx);
376 4816672 enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
377 4816672 viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
378 2408336 temperature_[phaseIdx] = fs.temperature(phaseIdx);
379 }
380 1204168 }
381
382 protected:
383 Scalar moleFraction_[numPhases][numComponents] = {};
384 Scalar fugacityCoefficient_[numPhases][numComponents] = {};
385 Scalar averageMolarMass_[numPhases] = {};
386 Scalar sumMoleFractions_[numPhases] = {};
387 Scalar pressure_[numPhases] = {};
388 Scalar saturation_[numPhases] = {};
389 Scalar density_[numPhases] = {};
390 Scalar molarDensity_[numPhases] = {};
391 Scalar enthalpy_[numPhases] = {};
392 Scalar viscosity_[numPhases] = {};
393 Scalar temperature_[numPhases] = {};
394 Scalar temperatureEquil_ = 0.0;
395
396 // For porous medium flow models, here we ... the index of the wetting
397 // phase (needed for vapor pressure evaluation if kelvin equation is used)
398 int wPhaseIdx_{0};
399 };
400
401 } // end namespace Dumux
402
403 #endif
404