GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/examples/biomineralization/material/fluidsystems/icpcomplexsalinitybrine.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 37 38 97.4%
Functions: 3 3 100.0%
Branches: 8 32 25.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-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7
8 #ifndef DUMUX_ICP_COMPLEX_SALINITY_BRINE_FLUID_SYSTEM_HH
9 #define DUMUX_ICP_COMPLEX_SALINITY_BRINE_FLUID_SYSTEM_HH
10
11 // ## The brine adapter (`icpcomplexsalinitybrine.hh`)
12 //
13 // This file contains the brineadapter class__ which defines
14 // how to adapt values for communication between
15 // the "full" fluidsystem accounting for 4 components influencing the brine properties
16 // and the "brine" fluidsystem assuming water and a single salt determining the brine properties.
17 //
18 // [[content]]
19 //
20 // ### Include files
21 // [[details]] includes
22 // we include all necessary components
23 #include <dune/common/math.hh>
24
25 #include <dumux/material/fluidsystems/base.hh>
26 #include <dumux/material/constants.hh>
27 #include <dumux/material/components/h2o.hh>
28 #include <dumux/material/components/sodiumion.hh>
29 #include <dumux/material/components/chlorideion.hh>
30 #include <dumux/material/components/calciumion.hh>
31 #include <dumux/material/components/tabulatedcomponent.hh>
32
33 #include <dumux/common/exceptions.hh>
34
35 #include <dumux/io/name.hh>
36 // [[/details]]
37 //
38 // ### The brine adapter class
39 // In ICPComplexSalinityBrine we make the transition from 3 components additional to water contributing to salinity to the single component that is assumed in the brine fluid system.
40 // We enter the namespace Dumux. All Dumux functions and classes are in a namespace Dumux, to make sure they don`t clash with symbols from other libraries you may want to use in conjunction with Dumux.
41 // [[codeblock]]
42 namespace Dumux::FluidSystems {
43
44 template< class Scalar, class H2OType = Components::TabulatedComponent<Dumux::Components::H2O<Scalar>> >
45 class ICPComplexSalinityBrine : public Base< Scalar, ICPComplexSalinityBrine<Scalar, H2OType>>
46 {
47 using ThisType = ICPComplexSalinityBrine<Scalar, H2OType>;
48 using Base = Dumux::FluidSystems::Base<Scalar, ThisType>;
49
50 public:
51 //! export the involved components
52 using H2O = H2OType;
53 using Na = Components::SodiumIon<Scalar>;
54 using Cl = Components::ChlorideIon<Scalar>;
55 using Ca = Components::CalciumIon<Scalar>;
56 // [[/codeblock]]
57 //
58 // #### Component and phase definitions
59 // With the following function we define what phases and components will be used by the fluid system and define the indices used to distinguish phases and components in the course of the simulation.
60 // [[codeblock]]
61 // We use convenient declarations that we derive from the property system.
62 static constexpr int numPhases = 1; //!< Number of phases in the fluid system
63 static constexpr int numComponents = 4; //!< Number of components in the fluid system (H2O, Na, Cl, Ca)
64
65 static constexpr int phase0Idx = 0; //!< Index of the first (and only) phase
66 static constexpr int liquidPhaseIdx = phase0Idx; //!< The one considered phase is liquid
67
68 static constexpr int H2OIdx = 0; //!< index of the water component
69 static constexpr int NaIdx = 1; //!< index of the Na component
70 static constexpr int ClIdx = 2; //!< index of the Cl component
71 static constexpr int CaIdx = 3; //!< index of the Ca component
72 static constexpr int comp0Idx = H2OIdx; //!< index of the first component
73 static constexpr int comp1Idx = NaIdx; //!< index of the second component
74 static constexpr int comp2Idx = ClIdx; //!< index of the third component
75 static constexpr int comp3Idx = CaIdx; //!< index of the fourth component
76
77 // the fluid phase index
78 static std::string phaseName(int phaseIdx = liquidPhaseIdx)
79 {
80 assert(phaseIdx == liquidPhaseIdx);
81 return IOName::liquidPhase();
82 }
83
84 // the miscibility
85 static constexpr bool isMiscible()
86 {
87 return false;
88 }
89
90 // we do not have a gas phase
91 static constexpr bool isGas(int phaseIdx = liquidPhaseIdx)
92 {
93 assert(phaseIdx == liquidPhaseIdx);
94 return false;
95 }
96
97 // We need a ideal mixture
98 static constexpr bool isIdealMixture(int phaseIdx = liquidPhaseIdx)
99 {
100 assert(phaseIdx == liquidPhaseIdx);
101 return true;
102 }
103
104 // phase compressibility
105 static constexpr bool isCompressible(int phaseIdx = liquidPhaseIdx)
106 {
107 assert(phaseIdx == liquidPhaseIdx);
108 return H2O::liquidIsCompressible();
109 }
110
111 // no gas, thus no ideal gas
112 static constexpr bool isIdealGas(int phaseIdx = liquidPhaseIdx)
113 {
114 assert(phaseIdx == liquidPhaseIdx);
115 return false; /*we're a liquid!*/
116 }
117 // [[/codeblock]]
118 // [[exclude]]
119 // The component names
120 static std::string componentName(int compIdx)
121 {
122 switch (compIdx)
123 {
124 case H2OIdx: return H2O::name();
125 case CaIdx: return Ca::name();
126 case NaIdx: return Na::name();
127 case ClIdx: return Cl::name();
128 };
129 DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
130 }
131
132 // The component molar masses
133
1/2
✓ Branch 0 taken 20766176 times.
✗ Branch 1 not taken.
 20766176 static Scalar molarMass(int compIdx)
134 {
135 switch (compIdx)
136 {
137 case H2OIdx: return H2O::molarMass();
138 case CaIdx: return Ca::molarMass();
139 case NaIdx: return Na::molarMass();
140 case ClIdx: return Cl::molarMass();
141 };
142 ✗ DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
143 }
144
145 //initializing
146 static void init()
147 {
148 init(/*tempMin=*/273.15,
149 /*tempMax=*/623.15,
150 /*numTemp=*/100,
151 /*pMin=*/-10.,
152 /*pMax=*/20e6,
153 /*numP=*/200);
154 }
155 static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
156 Scalar pressMin, Scalar pressMax, unsigned nPress)
157 {
158
159 if (H2O::isTabulated)
160 {
161 std::cout << "Initializing tables for the H2O fluid properties ("
162 << nTemp*nPress
163 << " entries).\n";
164
165 H2O::init(tempMin, tempMax, nTemp, pressMin, pressMax, nPress);
166 }
167 }
168 // [[/exclude]]
169 //
170 // #### The Component-Phase Interactions
171 // In the following, we define the component-phase interactions such as // each component's fugacity coefficient in brine
172 // and each component's binary diffusion coefficient with water in brine
173 // [[codeblock]]
174 // The component fugacity coefficients
175 template <class FluidState>
176 static Scalar fugacityCoefficient(const FluidState &fluidState,
177 int phaseIdx,
178 int compIdx)
179 {
180 assert(0 <= phaseIdx && phaseIdx < numPhases);
181 assert(0 <= compIdx && compIdx < numComponents);
182
183 if (phaseIdx == compIdx)
184 // We could calculate the real fugacity coefficient of
185 // the component in the fluid. Probably that's not worth
186 // the effort, since the fugacity coefficient of the other
187 // component is infinite anyway...
188 return 1.0;
189 return std::numeric_limits<Scalar>::infinity();
190 }
191
192 // The binary diffusion coefficients of the components and the phases main component
193 template <class FluidState>
194 static Scalar binaryDiffusionCoefficient(const FluidState& fluidState,
195 int phaseIdx,
196 int compIIdx,
197 int compJIdx)
198 {
199 if (phaseIdx == liquidPhaseIdx)
200 {
201 if (compIIdx > compJIdx)
202 {
203 using std::swap;
204 swap(compIIdx, compJIdx);
205 }
206 if (compJIdx == NaIdx || compJIdx == ClIdx || compJIdx == CaIdx)
207 return 1.587e-9; //[m²/s] //J. Phys. D: Appl. Phys. 40 (2007) 2769-2776
208 else
209 DUNE_THROW(Dune::NotImplemented, "Binary diffusion coefficient of components "
210 << compIIdx << " and " << compJIdx
211 << " in phase " << phaseIdx);
212 }
213
214 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index: " << phaseIdx);
215 }
216 // [[/codeblock]]
217 //
218 // #### The Fluid (Brine) Properties
219 // In the following, all functions defining the phase properties are given:
220 // the density, viscosity, enthalpy, thermal conductivities, and heat capacities
221 // of each phase depending on temperature, pressure, and phase composition
222 // [[codeblock]]
223 // The phase density
224 template <class FluidState>
225 5191544 static Scalar density(const FluidState& fluidState, int phaseIdx = liquidPhaseIdx)
226 {
227
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5191544 times.
 5191544 assert(phaseIdx == liquidPhaseIdx);
228 5191544 const Scalar temperature = fluidState.temperature(phaseIdx);
229 5191544 const Scalar pressure = fluidState.pressure(phaseIdx);
230 5191544 const Scalar xNaCl = fluidState.massFraction(phaseIdx, NaIdx)
231 5191544 + fluidState.massFraction(phaseIdx, ClIdx)
232 5191544 + fluidState.massFraction(phaseIdx, CaIdx);
233
234 using std::max;
235 5191544 const Scalar TempC = temperature - 273.15;
236 5191544 const Scalar pMPa = pressure/1.0E6;
237
1/2
✓ Branch 0 taken 5191544 times.
✗ Branch 1 not taken.
 5191544 const Scalar salinity = max(0.0, xNaCl);
238
239 5191544 const Scalar rhow = H2O::liquidDensity(temperature, pressure);
240 15574632 const Scalar density = rhow + 1000*salinity*(0.668 +
241 10383088 0.44*salinity +
242 15574632 1.0E-6*(300*pMPa -
243 10383088 2400*pMPa*salinity +
244 10383088 TempC*(80.0 +
245 10383088 3*TempC -
246 10383088 3300*salinity -
247 10383088 13*pMPa +
248 5191544 47*pMPa*salinity)));
249
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5191544 times.
 5191544 assert(density > 0.0);
250 5191544 return density;
251 }
252
253
254
255 // The phase viscosity
256 template <class FluidState>
257 2595772 static Scalar viscosity(const FluidState& fluidState, int phaseIdx = liquidPhaseIdx)
258 {
259
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2595772 times.
 2595772 assert(phaseIdx == liquidPhaseIdx);
260 2595772 const Scalar temperature = fluidState.temperature(phaseIdx);
261 2595772 const Scalar xNaCl = fluidState.massFraction(phaseIdx, NaIdx)
262 2595772 + fluidState.massFraction(phaseIdx, ClIdx)
263 2595772 + fluidState.massFraction(phaseIdx, CaIdx);
264
265 using std::pow;
266 using Dune::power;
267 using std::exp;
268 using std::max;
269
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2595772 times.
 2595772 const Scalar T = max(temperature, 275.0);
270
1/2
✓ Branch 0 taken 2595772 times.
✗ Branch 1 not taken.
 2595772 const Scalar salinity = max(0.0, xNaCl);
271
272 2595772 const Scalar T_C = T - 273.15;
273 2595772 const Scalar A = ((0.42*power((pow(salinity, 0.8)-0.17), 2)) + 0.045)*pow(T_C, 0.8);
274 2595772 const Scalar mu_brine = 0.1 + (0.333*salinity) + (1.65+(91.9*salinity*salinity*salinity))*exp(-A); // [cP]
275
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2595772 times.
 2595772 assert(mu_brine > 0.0);
276 2595772 return mu_brine/1000.0; // [Pa·s]
277 }
278 // [[/codeblock]]
279 // [[exclude]]
280 // The vapor pressure
281 template <class FluidState>
282 static Scalar vaporPressure(const FluidState& fluidState, int compIdx)
283 {
284 if (compIdx == H2OIdx)
285 {
286 // simplified version of Eq 2.29 in Vishal Jambhekar's Promo
287 const Scalar temperature = fluidState.temperature(H2OIdx);
288 return H2O::vaporPressure(temperature)/fluidState.massFraction(phase0Idx, H2OIdx);
289 }
290 else if (compIdx == NaIdx)
291 DUNE_THROW(Dune::NotImplemented, "Na::vaporPressure(t)");
292 else if (compIdx == ClIdx)
293 DUNE_THROW(Dune::NotImplemented, "Cl::vaporPressure(t)");
294 else if (compIdx == CaIdx)
295 DUNE_THROW(Dune::NotImplemented, "Ca::vaporPressure(t)");
296 else
297 DUNE_THROW(Dune::NotImplemented, "Invalid component index " << compIdx);
298 }
299
300 // The phase enthalpies
301 template <class FluidState>
302 static Scalar enthalpy(const FluidState& fluidState, int phaseIdx)
303 {
304 /*Numerical coefficients from PALLISER*/
305 static constexpr Scalar f[] = { 2.63500E-1, 7.48368E-6, 1.44611E-6, -3.80860E-10 };
306
307 /*Numerical coefficients from MICHAELIDES for the enthalpy of brine*/
308 static constexpr Scalar a[4][3] = { { +9633.6, -4080.0, +286.49 },
309 { +166.58, +68.577, -4.6856 },
310 { -0.90963, -0.36524, +0.249667E-1 },
311 { +0.17965E-2, +0.71924E-3, -0.4900E-4 } };
312
313 const Scalar p = fluidState.pressure(phaseIdx);
314 const Scalar T = fluidState.temperature(phaseIdx);
315 const Scalar theta = T - 273.15;
316 const Scalar salSat = f[0] + f[1]*theta + f[2]*theta*theta + f[3]*theta*theta*theta;
317
318 /*Regularization*/
319 using std::min;
320 using std::max;
321 const Scalar xNaCl = fluidState.massFraction(phaseIdx, NaIdx)
322 + fluidState.massFraction(phaseIdx, ClIdx)
323 + fluidState.massFraction(phaseIdx, CaIdx);
324 const Scalar salinity = min(max(xNaCl,0.0), salSat);
325
326 const Scalar hw = H2O::liquidEnthalpy(T, p)/1E3; /* kJ/kg */
327
328 /*DAUBERT and DANNER*/
329 const Scalar h_NaCl = (3.6710E4*T + 0.5*(6.2770E1)*T*T - ((6.6670E-2)/3)*T*T*T
330 + ((2.8000E-5)/4)*(T*T*T*T))/(58.44E3)- 2.045698e+02; /* U [kJ/kg] */
331
332 const Scalar m = (1E3/58.44)*(salinity/(1-salinity));
333
334 using Dune::power;
335 Scalar d_h = 0;
336 for (int i = 0; i<=3; i++)
337 for (int j=0; j<=2; j++)
338 d_h = d_h + a[i][j] * power(theta, i) * power(m, j);
339
340 /* heat of dissolution for halite according to Michaelides 1971 */
341 const Scalar delta_h = (4.184/(1E3 + (58.44 * m)))*d_h;
342
343 /* Enthalpy of brine without any dissolved gas */
344 const Scalar h_ls1 =(1-salinity)*hw + salinity*h_NaCl + salinity*delta_h; /* kJ/kg */
345 return h_ls1*1E3; /*J/kg*/
346 }
347 // The molar density
348 template <class FluidState>
349 static Scalar molarDensity(const FluidState& fluidState, int phaseIdx = liquidPhaseIdx)
350 {
351 5191544 return density(fluidState, phaseIdx)/fluidState.averageMolarMass(phaseIdx);
352 }
353
354 template <class FluidState>
355 static Scalar diffusionCoefficient(const FluidState& fluidState, int phaseIdx, int compIdx)
356 {
357 DUNE_THROW(Dune::NotImplemented, "FluidSystems::ICPComplexSalinityBrine::diffusionCoefficient()");
358 }
359
360 // The thermal conductivity
361 template <class FluidState>
362 static Scalar thermalConductivity(const FluidState& fluidState, int phaseIdx)
363 {
364 if (phaseIdx == liquidPhaseIdx)
365 return H2O::liquidThermalConductivity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
366 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index: " << phaseIdx);
367 }
368
369 // The phase heat capacity
370 template <class FluidState>
371 static Scalar heatCapacity(const FluidState &fluidState, int phaseIdx)
372 {
373 if (phaseIdx == liquidPhaseIdx)
374 return H2O::liquidHeatCapacity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
375 DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
376 }
377 };
378
379 } // end namespace Dumux::FluidSystems
380 // [[/exclude]]
381 // [[/content]]
382 #endif
383