https://mooseframework.inl.gov
Loading...
Searching...
No Matches
PorousFlowBrineCO2Test.C
Go to the documentation of this file.
1//* This file is part of the MOOSE framework
2//* https://mooseframework.inl.gov
3//*
4//* All rights reserved, see COPYRIGHT for full restrictions
5//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6//*
7//* Licensed under LGPL 2.1, please see LICENSE for details
8//* https://www.gnu.org/licenses/lgpl-2.1.html
9
10#include "PorousFlowBrineCO2Test.h"
11#include "FluidPropertiesTestUtils.h"
12
16TEST_F(PorousFlowBrineCO2Test, name) { EXPECT_EQ("brine-co2", _fs->fluidStateName()); }
17
21TEST_F(PorousFlowBrineCO2Test, indices)
22{
23 EXPECT_EQ((unsigned int)2, _fs->numPhases());
24 EXPECT_EQ((unsigned int)3, _fs->numComponents());
25 EXPECT_EQ((unsigned int)0, _fs->aqueousPhaseIndex());
26 EXPECT_EQ((unsigned int)1, _fs->gasPhaseIndex());
27 EXPECT_EQ((unsigned int)0, _fs->aqueousComponentIndex());
28 EXPECT_EQ((unsigned int)1, _fs->gasComponentIndex());
29 EXPECT_EQ((unsigned int)2, _fs->saltComponentIndex());
30}
31
32/*
33 * Verify calculation of the equilibrium constants and their derivatives wrt temperature
34 */
35TEST_F(PorousFlowBrineCO2Test, equilibriumConstants)
36{
37 ADReal T = 350.0;
38 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
39
40 const Real dT = 1.0e-6;
41
42 ADReal K0H2O = _fs->equilibriumConstantH2O(T);
43 ADReal K0CO2 = _fs->equilibriumConstantCO2(T);
44
45 ABS_TEST(K0H2O.value(), 0.412597711705, 1.0e-10);
46 ABS_TEST(K0CO2.value(), 74.0435888596, 1.0e-10);
47
48 Real K0H2O_2 = _fs->equilibriumConstantH2O(T + dT).value();
49 Real K0CO2_2 = _fs->equilibriumConstantCO2(T + dT).value();
50
51 Real dK0H2O_dT_fd = (K0H2O_2 - K0H2O.value()) / dT;
52 Real dK0CO2_dT_fd = (K0CO2_2 - K0CO2.value()) / dT;
53 REL_TEST(K0H2O.derivatives()[_Tidx], dK0H2O_dT_fd, 1.0e-6);
54 REL_TEST(K0CO2.derivatives()[_Tidx], dK0CO2_dT_fd, 1.0e-6);
55
56 // Test the high temperature formulation
57 T = 450.0;
58 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
59
60 K0H2O = _fs->equilibriumConstantH2O(T);
61 K0CO2 = _fs->equilibriumConstantCO2(T);
62
63 ABS_TEST(K0H2O.value(), 8.75944517916, 1.0e-10);
64 ABS_TEST(K0CO2.value(), 105.013867434, 1.0e-10);
65
66 K0H2O_2 = _fs->equilibriumConstantH2O(T + dT).value();
67 K0CO2_2 = _fs->equilibriumConstantCO2(T + dT).value();
68
69 dK0H2O_dT_fd = (K0H2O_2 - K0H2O.value()) / dT;
70 dK0CO2_dT_fd = (K0CO2_2 - K0CO2.value()) / dT;
71 REL_TEST(K0H2O.derivatives()[_Tidx], dK0H2O_dT_fd, 1.0e-6);
72 REL_TEST(K0CO2.derivatives()[_Tidx], dK0CO2_dT_fd, 1.0e-6);
73}
74
75/*
76 * Verify calculation of the fugacity coefficients
77 */
78TEST_F(PorousFlowBrineCO2Test, fugacityCoefficients)
79{
80 // Test the low temperature formulation
81 ADReal p = 40.0e6;
82 Moose::derivInsert(p.derivatives(), _pidx, 1.0);
83
84 ADReal T = 350.0;
85 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
86
87 ADReal Xnacl = 0.0;
88 Moose::derivInsert(Xnacl.derivatives(), _Xidx, 1.0);
89
90 ADReal co2_density = _co2_fp->rho_from_p_T(p, T);
91
92 ADReal phiH2O, phiCO2;
93 _fs->fugacityCoefficientsLowTemp(p, T, co2_density, phiCO2, phiH2O);
94 ABS_TEST(phiCO2.value(), 0.401939386415, 1.0e-8);
95 ABS_TEST(phiH2O.value(), 0.0898968578757, 1.0e-8);
96
97 // Test the high temperature formulation
98 T = 423.15;
99
100 ADReal xco2, yh2o;
101 co2_density = _co2_fp->rho_from_p_T(p, T);
102
103 _fs->solveEquilibriumMoleFractionHighTemp(
104 p.value(), T.value(), Xnacl.value(), co2_density.value(), xco2.value(), yh2o.value());
105 phiH2O = _fs->fugacityCoefficientH2OHighTemp(p, T, co2_density, xco2, yh2o);
106 phiCO2 = _fs->fugacityCoefficientCO2HighTemp(p, T, co2_density, xco2, yh2o);
107
108 ABS_TEST(phiH2O.value(), 0.156303437579, 1.0e-8);
109 ABS_TEST(phiCO2.value(), 0.641936639599, 1.0e-8);
110
111 // Test that the same results are returned in fugacityCoefficientsHighTemp()
112 ADReal phiCO2_2, phiH2O_2;
113 _fs->fugacityCoefficientsHighTemp(p, T, co2_density, xco2, yh2o, phiCO2_2, phiH2O_2);
114
115 ABS_TEST(phiH2O, phiH2O_2, 1.0e-12);
116 ABS_TEST(phiCO2, phiCO2_2, 1.0e-12);
117}
118
119/*
120 * Verify calculation of the H2O and CO2 activity coefficients
121 */
122TEST_F(PorousFlowBrineCO2Test, activityCoefficients)
123{
124 ADReal xco2 = 0.01;
125 ADReal T = 350.0;
126
127 ADReal gammaH2O = _fs->activityCoefficientH2O(T, xco2);
128 ADReal gammaCO2 = _fs->activityCoefficientCO2(T, xco2);
129
130 ABS_TEST(gammaH2O.value(), 1.0, 1.0e-10);
131 ABS_TEST(gammaCO2.value(), 1.0, 1.0e-10);
132
133 T = 450.0;
134
135 gammaH2O = _fs->activityCoefficientH2O(T, xco2);
136 gammaCO2 = _fs->activityCoefficientCO2(T, xco2);
137
138 ABS_TEST(gammaH2O.value(), 1.00022113664, 1.0e-10);
139 ABS_TEST(gammaCO2.value(), 0.956736800255, 1.0e-10);
140}
141
142/*
143 * Verify calculation of the Duan and Sun activity coefficient and its derivatives wrt
144 * pressure, temperature and salt mass fraction
145 */
146TEST_F(PorousFlowBrineCO2Test, activityCoefficientCO2Brine)
147{
148 ADReal p = 10.0e6;
149 Moose::derivInsert(p.derivatives(), _pidx, 1.0);
150
151 ADReal T = 350.0;
152 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
153
154 ADReal Xnacl = 0.1;
155 Moose::derivInsert(Xnacl.derivatives(), _Xidx, 1.0);
156
157 const Real dp = 1.0e-1;
158 const Real dT = 1.0e-6;
159 const Real dx = 1.0e-8;
160
161 // Low temperature regime
162 ADReal gamma = _fs->activityCoefficient(p, T, Xnacl);
163 ABS_TEST(gamma.value(), 1.43276649338, 1.0e-8);
164
165 ADReal gamma_2 = _fs->activityCoefficient(p + dp, T, Xnacl);
166
167 Real dgamma_dp_fd = (gamma_2.value() - gamma.value()) / dp;
168 REL_TEST(gamma.derivatives()[_pidx], dgamma_dp_fd, 1.0e-6);
169
170 gamma_2 = _fs->activityCoefficient(p, T + dT, Xnacl);
171
172 Real dgamma_dT_fd = (gamma_2.value() - gamma.value()) / dT;
173 REL_TEST(gamma.derivatives()[_Tidx], dgamma_dT_fd, 1.0e-6);
174
175 gamma_2 = _fs->activityCoefficient(p, T, Xnacl + dx);
176
177 Real dgamma_dX_fd = (gamma_2.value() - gamma.value()) / dx;
178 REL_TEST(gamma.derivatives()[_Xidx], dgamma_dX_fd, 1.0e-6);
179
180 // High temperature regime
181 T = 523.15;
182 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
183
184 gamma = _fs->activityCoefficientHighTemp(T, Xnacl);
185 ABS_TEST(gamma.value(), 1.50047006243, 1.0e-8);
186
187 gamma_2 = _fs->activityCoefficientHighTemp(T + dT, Xnacl);
188 dgamma_dT_fd = (gamma_2.value() - gamma.value()) / dT;
189 REL_TEST(gamma.derivatives()[_Tidx], dgamma_dT_fd, 1.0e-6);
190
191 gamma_2 = _fs->activityCoefficientHighTemp(T, Xnacl + dx);
192 dgamma_dX_fd = (gamma_2.value() - gamma.value()) / dx;
193 REL_TEST(gamma.derivatives()[_Xidx], dgamma_dX_fd, 1.0e-6);
194
195 // Check that both formulations return gamma = 1 for Xnacl = 0
196 Xnacl = 0.0;
197 T = 350.0;
198 gamma = _fs->activityCoefficient(p, T, Xnacl);
199 ABS_TEST(gamma.value(), 1.0, 1.0e-12);
200
201 T = 523.15;
202 gamma = _fs->activityCoefficientHighTemp(T, Xnacl);
203 ABS_TEST(gamma.value(), 1.0, 1.0e-12);
204}
205
206/*
207 * Verify calculation of the partial density of CO2 and its derivative wrt temperature
208 */
209TEST_F(PorousFlowBrineCO2Test, partialDensity)
210{
211 ADReal T = 473.15;
212 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
213
214 const Real dT = 1.0e-6;
215
216 ADReal partial_density = _fs->partialDensityCO2(T);
217 ABS_TEST(partial_density.value(), 893.332, 1.0e-3);
218
219 ADReal partial_density_2 = _fs->partialDensityCO2(T + dT);
220
221 Real dpartial_density_dT_fd = (partial_density_2.value() - partial_density.value()) / dT;
222 REL_TEST(partial_density.derivatives()[_Tidx], dpartial_density_dT_fd, 1.0e-6);
223}
224
225/*
226 * Verify calculation of equilibrium mass fraction and derivatives
227 */
228TEST_F(PorousFlowBrineCO2Test, equilibriumMassFraction)
229{
230 // Low temperature regime
231 ADReal p = 1.0e6;
232 Moose::derivInsert(p.derivatives(), _pidx, 1.0);
233
234 ADReal T = 350.0;
235 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
236
237 ADReal Xnacl = 0.1;
238 Moose::derivInsert(Xnacl.derivatives(), _Xidx, 1.0);
239
240 const Real dp = 1.0e-2;
241 const Real dT = 1.0e-6;
242 const Real dx = 1.0e-8;
243
244 ADReal X, Y;
245 _fs->equilibriumMassFractions(p, T, Xnacl, X, Y);
246
247 ABS_TEST(X.value(), 0.0035573020148, 1.0e-10);
248 ABS_TEST(Y.value(), 0.0171977397214, 1.0e-10);
249
250 // Derivative wrt pressure
251 ADReal X1, Y1, X2, Y2;
252 _fs->equilibriumMassFractions(p - dp, T, Xnacl, X1, Y1);
253 _fs->equilibriumMassFractions(p + dp, T, Xnacl, X2, Y2);
254
255 Real dX_dp_fd = (X2.value() - X1.value()) / (2.0 * dp);
256 Real dY_dp_fd = (Y2.value() - Y1.value()) / (2.0 * dp);
257
258 REL_TEST(X.derivatives()[_pidx], dX_dp_fd, 1.0e-6);
259 REL_TEST(Y.derivatives()[_pidx], dY_dp_fd, 1.0e-6);
260
261 // Derivative wrt temperature
262 _fs->equilibriumMassFractions(p, T - dT, Xnacl, X1, Y1);
263 _fs->equilibriumMassFractions(p, T + dT, Xnacl, X2, Y2);
264
265 Real dX_dT_fd = (X2.value() - X1.value()) / (2.0 * dT);
266 Real dY_dT_fd = (Y2.value() - Y1.value()) / (2.0 * dT);
267
268 REL_TEST(X.derivatives()[_Tidx], dX_dT_fd, 1.0e-6);
269 REL_TEST(Y.derivatives()[_Tidx], dY_dT_fd, 1.0e-6);
270
271 // Derivative wrt salt mass fraction
272 _fs->equilibriumMassFractions(p, T, Xnacl - dx, X1, Y1);
273 _fs->equilibriumMassFractions(p, T, Xnacl + dx, X2, Y2);
274
275 Real dX_dX_fd = (X2.value() - X1.value()) / (2.0 * dx);
276 Real dY_dX_fd = (Y2.value() - Y1).value() / (2.0 * dx);
277
278 REL_TEST(X.derivatives()[_Xidx], dX_dX_fd, 1.0e-6);
279 REL_TEST(Y.derivatives()[_Xidx], dY_dX_fd, 1.0e-6);
280
281 // High temperature regime
282 p = 10.0e6;
283 Moose::derivInsert(p.derivatives(), _pidx, 1.0);
284
285 T = 525.15;
286 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
287
288 _fs->equilibriumMassFractions(p, T, Xnacl, X, Y);
289
290 ABS_TEST(X.value(), 0.016299479086, 1.0e-10);
291 ABS_TEST(Y.value(), 0.249471400766, 1.0e-10);
292}
293
294/*
295 * Verify calculation of actual mass fraction and derivatives depending on value of
296 * total mass fraction
297 */
298TEST_F(PorousFlowBrineCO2Test, MassFraction)
299{
300 ADReal p = 1.0e6;
301 Moose::derivInsert(p.derivatives(), _pidx, 1.0);
302
303 ADReal T = 350.0;
304 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
305
306 ADReal Xnacl = 0.1;
307 Moose::derivInsert(Xnacl.derivatives(), _Xidx, 1.0);
308
309 FluidStatePhaseEnum phase_state;
310 const unsigned int np = _fs->numPhases();
311 const unsigned int nc = _fs->numComponents();
312 std::vector<FluidStateProperties> fsp(np, FluidStateProperties(nc));
313
314 // Liquid region
315 ADReal Z = 0.0001;
316 Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
317
318 _fs->massFractions(p, T, Xnacl, Z, phase_state, fsp);
319 EXPECT_EQ(phase_state, FluidStatePhaseEnum::LIQUID);
320
321 // Verfify mass fraction values
322 ADReal Xco2 = fsp[0].mass_fraction[1];
323 ADReal Yco2 = fsp[1].mass_fraction[1];
324 ADReal Xh2o = fsp[0].mass_fraction[0];
325 ADReal Yh2o = fsp[1].mass_fraction[0];
326 ADReal Xnacl2 = fsp[0].mass_fraction[2];
327 ABS_TEST(Xco2.value(), Z.value(), 1.0e-8);
328 ABS_TEST(Yco2.value(), 0.0, 1.0e-8);
329 ABS_TEST(Xh2o.value(), 1.0 - Z.value(), 1.0e-8);
330 ABS_TEST(Yh2o.value(), 0.0, 1.0e-8);
331 ABS_TEST(Xnacl2.value(), Xnacl.value(), 1.0e-8);
332
333 // Verify derivatives
334 ABS_TEST(Xco2.derivatives()[_pidx], 0.0, 1.0e-8);
335 ABS_TEST(Xco2.derivatives()[_Tidx], 0.0, 1.0e-8);
336 ABS_TEST(Xco2.derivatives()[_Xidx], 0.0, 1.0e-8);
337 ABS_TEST(Xco2.derivatives()[_Zidx], 1.0, 1.0e-8);
338 ABS_TEST(Yco2.derivatives()[_pidx], 0.0, 1.0e-8);
339 ABS_TEST(Yco2.derivatives()[_Tidx], 0.0, 1.0e-8);
340 ABS_TEST(Yco2.derivatives()[_Xidx], 0.0, 1.0e-8);
341 ABS_TEST(Yco2.derivatives()[_Zidx], 0.0, 1.0e-8);
342 ABS_TEST(Xnacl.derivatives()[_Xidx], 1.0, 1.0e-8);
343
344 // Gas region
345 Z = 0.995;
346 Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
347
348 _fs->massFractions(p, T, Xnacl, Z, phase_state, fsp);
349 EXPECT_EQ(phase_state, FluidStatePhaseEnum::GAS);
350
351 // Verfify mass fraction values
352 Xco2 = fsp[0].mass_fraction[1];
353 Yco2 = fsp[1].mass_fraction[1];
354 Xh2o = fsp[0].mass_fraction[0];
355 Yh2o = fsp[1].mass_fraction[0];
356 ADReal Ynacl = fsp[1].mass_fraction[2];
357 ABS_TEST(Xco2.value(), 0.0, 1.0e-8);
358 ABS_TEST(Yco2.value(), Z.value(), 1.0e-8);
359 ABS_TEST(Xh2o.value(), 0.0, 1.0e-8);
360 ABS_TEST(Yh2o.value(), 1.0 - Z.value(), 1.0e-8);
361 ABS_TEST(Ynacl.value(), 0.0, 1.0e-8);
362
363 // Verify derivatives
364 ABS_TEST(Xco2.derivatives()[_pidx], 0.0, 1.0e-8);
365 ABS_TEST(Xco2.derivatives()[_Tidx], 0.0, 1.0e-8);
366 ABS_TEST(Xco2.derivatives()[_Xidx], 0.0, 1.0e-8);
367 ABS_TEST(Xco2.derivatives()[_Zidx], 0.0, 1.0e-8);
368 ABS_TEST(Yco2.derivatives()[_pidx], 0.0, 1.0e-8);
369 ABS_TEST(Yco2.derivatives()[_Tidx], 0.0, 1.0e-8);
370 ABS_TEST(Yco2.derivatives()[_Xidx], 0.0, 1.0e-8);
371 ABS_TEST(Yco2.derivatives()[_Zidx], 1.0, 1.0e-8);
372 ABS_TEST(Ynacl.derivatives()[_Xidx], 0.0, 1.0e-8);
373
374 // Two phase region. In this region, the mass fractions and derivatives can
375 // be verified using the equilibrium mass fraction derivatives that have
376 // been verified above
377 Z = 0.45;
378 Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
379
380 _fs->massFractions(p, T, Xnacl, Z, phase_state, fsp);
381 EXPECT_EQ(phase_state, FluidStatePhaseEnum::TWOPHASE);
382
383 // Equilibrium mass fractions and derivatives
384 ADReal Xco2_eq, Yh2o_eq;
385 _fs->equilibriumMassFractions(p, T, Xnacl, Xco2_eq, Yh2o_eq);
386
387 // Verfify mass fraction values - comparing ADReals verifies that derivatives
388 // are also identical
389 Xco2 = fsp[0].mass_fraction[1];
390 Yco2 = fsp[1].mass_fraction[1];
391 Xh2o = fsp[0].mass_fraction[0];
392 Yh2o = fsp[1].mass_fraction[0];
393 ABS_TEST(Xco2, Xco2_eq, 1.0e-8);
394 ABS_TEST(Yco2, 1.0 - Yh2o_eq, 1.0e-8);
395 ABS_TEST(Xh2o, 1.0 - Xco2_eq, 1.0e-8);
396 ABS_TEST(Yh2o, Yh2o_eq, 1.0e-8);
397
398 // Use finite differences to verify derivative wrt Z is unaffected by Z
399 const Real dZ = 1.0e-8;
400 _fs->massFractions(p, T, Xnacl, Z + dZ, phase_state, fsp);
401 ADReal Xco21 = fsp[0].mass_fraction[1];
402 ADReal Yco21 = fsp[1].mass_fraction[1];
403 _fs->massFractions(p, T, Xnacl, Z - dZ, phase_state, fsp);
404 ADReal Xco22 = fsp[0].mass_fraction[1];
405 ADReal Yco22 = fsp[1].mass_fraction[1];
406
407 ABS_TEST(Xco2.derivatives()[_Zidx], (Xco21.value() - Xco22.value()) / (2.0 * dZ), 1.0e-8);
408 ABS_TEST(Yco2.derivatives()[_Zidx], (Yco21.value() - Yco22.value()) / (2.0 * dZ), 1.0e-8);
409}
410
411/*
412 * Verify calculation of gas density, viscosity enthalpy, and derivatives. Note that as
413 * these properties don't depend on mass fraction, only the gas region needs to be
414 * tested (the calculations are identical in the two phase region)
415 */
416TEST_F(PorousFlowBrineCO2Test, gasProperties)
417{
418 ADReal p = 1.0e6;
419 Moose::derivInsert(p.derivatives(), _pidx, 1.0);
420
421 ADReal T = 350.0;
422 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
423
424 ADReal Xnacl = 0.1;
425 Moose::derivInsert(Xnacl.derivatives(), _Xidx, 1.0);
426
427 FluidStatePhaseEnum phase_state;
428 const unsigned int np = _fs->numPhases();
429 const unsigned int nc = _fs->numComponents();
430 std::vector<FluidStateProperties> fsp(np, FluidStateProperties(nc));
431
432 // Gas region
433 ADReal Z = 0.995;
434 Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
435
436 _fs->massFractions(p, T, Xnacl, Z, phase_state, fsp);
437 EXPECT_EQ(phase_state, FluidStatePhaseEnum::GAS);
438
439 // Verify fluid density, viscosity and enthalpy
440 _fs->gasProperties(p, T, fsp);
441 ADReal gas_density = fsp[1].density;
442 ADReal gas_viscosity = fsp[1].viscosity;
443 ADReal gas_enthalpy = fsp[1].enthalpy;
444
445 Real density = _co2_fp->rho_from_p_T(p.value(), T.value());
446 Real viscosity = _co2_fp->mu_from_p_T(p.value(), T.value());
447 Real enthalpy = _co2_fp->h_from_p_T(p.value(), T.value());
448
449 ABS_TEST(gas_density.value(), density, 1.0e-8);
450 ABS_TEST(gas_viscosity.value(), viscosity, 1.0e-8);
451 ABS_TEST(gas_enthalpy.value(), enthalpy, 1.0e-8);
452
453 // Verify derivatives
454 const Real dp = 1.0e-1;
455 _fs->gasProperties(p + dp, T, fsp);
456 Real rho1 = fsp[1].density.value();
457 Real mu1 = fsp[1].viscosity.value();
458 Real h1 = fsp[1].enthalpy.value();
459
460 _fs->gasProperties(p - dp, T, fsp);
461 Real rho2 = fsp[1].density.value();
462 Real mu2 = fsp[1].viscosity.value();
463 Real h2 = fsp[1].enthalpy.value();
464
465 REL_TEST(gas_density.derivatives()[_pidx], (rho1 - rho2) / (2.0 * dp), 1.0e-6);
466 REL_TEST(gas_viscosity.derivatives()[_pidx], (mu1 - mu2) / (2.0 * dp), 1.0e-6);
467 REL_TEST(gas_enthalpy.derivatives()[_pidx], (h1 - h2) / (2.0 * dp), 1.0e-6);
468
469 const Real dT = 1.0e-3;
470 _fs->gasProperties(p, T + dT, fsp);
471 rho1 = fsp[1].density.value();
472 mu1 = fsp[1].viscosity.value();
473 h1 = fsp[1].enthalpy.value();
474
475 _fs->gasProperties(p, T - dT, fsp);
476 rho2 = fsp[1].density.value();
477 mu2 = fsp[1].viscosity.value();
478 h2 = fsp[1].enthalpy.value();
479
480 REL_TEST(gas_density.derivatives()[_Tidx], (rho1 - rho2) / (2.0 * dT), 1.0e-6);
481 REL_TEST(gas_viscosity.derivatives()[_Tidx], (mu1 - mu2) / (2.0 * dT), 1.0e-6);
482 REL_TEST(gas_enthalpy.derivatives()[_Tidx], (h1 - h2) / (2.0 * dT), 1.0e-6);
483
484 // Note: mass fraction changes with Z
485 const Real dZ = 1.0e-8;
486 _fs->massFractions(p, T, Xnacl, Z + dZ, phase_state, fsp);
487 _fs->gasProperties(p, T, fsp);
488 rho1 = fsp[1].density.value();
489 mu1 = fsp[1].viscosity.value();
490 h1 = fsp[1].enthalpy.value();
491
492 _fs->massFractions(p, T, Xnacl, Z - dZ, phase_state, fsp);
493 _fs->gasProperties(p, T, fsp);
494 rho2 = fsp[1].density.value();
495 mu2 = fsp[1].viscosity.value();
496 h2 = fsp[1].enthalpy.value();
497
498 ABS_TEST(gas_density.derivatives()[_Zidx], (rho1 - rho2) / (2.0 * dZ), 1.0e-8);
499 ABS_TEST(gas_viscosity.derivatives()[_Zidx], (mu1 - mu2) / (2.0 * dZ), 1.0e-8);
500 ABS_TEST(gas_enthalpy.derivatives()[_Zidx], (h1 - h2) / (2.0 * dZ), 1.0e-6);
501}
502
503/*
504 * Verify calculation of liquid density, viscosity, enthalpy, and derivatives
505 */
506TEST_F(PorousFlowBrineCO2Test, liquidProperties)
507{
508 ADReal p = 1.0e6;
509 Moose::derivInsert(p.derivatives(), _pidx, 1.0);
510
511 ADReal T = 350.0;
512 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
513
514 ADReal Xnacl = 0.1;
515 Moose::derivInsert(Xnacl.derivatives(), _Xidx, 1.0);
516
517 FluidStatePhaseEnum phase_state;
518 const unsigned int np = _fs->numPhases();
519 const unsigned int nc = _fs->numComponents();
520 std::vector<FluidStateProperties> fsp(np, FluidStateProperties(nc));
521
522 // Liquid region
523 ADReal Z = 0.0001;
524 Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
525
526 _fs->massFractions(p, T, Xnacl, Z, phase_state, fsp);
527 EXPECT_EQ(phase_state, FluidStatePhaseEnum::LIQUID);
528
529 // Verify fluid density and viscosity
530 _fs->liquidProperties(p, T, Xnacl, fsp);
531
532 ADReal liquid_density = fsp[0].density;
533 ADReal liquid_viscosity = fsp[0].viscosity;
534 ADReal liquid_enthalpy = fsp[0].enthalpy;
535
536 Real co2_partial_density = _fs->partialDensityCO2(T).value();
537 Real brine_density = _brine_fp->rho_from_p_T_X(p.value(), T.value(), Xnacl.value());
538
539 Real density = 1.0 / (Z.value() / co2_partial_density + (1.0 - Z.value()) / brine_density);
540
541 Real viscosity = _brine_fp->mu_from_p_T_X(p.value(), T.value(), Xnacl.value());
542
543 Real brine_enthalpy = _brine_fp->h_from_p_T_X(p.value(), T.value(), Xnacl.value());
544 Real hdis = _fs->enthalpyOfDissolution(T).value();
545 Real co2_enthalpy = _co2_fp->h_from_p_T(p.value(), T.value());
546 Real enthalpy = (1.0 - Z.value()) * brine_enthalpy + Z.value() * (co2_enthalpy + hdis);
547
548 ABS_TEST(liquid_density.value(), density, 1.0e-12);
549 ABS_TEST(liquid_viscosity.value(), viscosity, 1.0e-12);
550 // Tolerance loosened to pass on Apple M4
551 ABS_TEST(liquid_enthalpy.value(), enthalpy, 3.e-10);
552
553 // Verify derivatives
554 // Derivatives wrt pressure
555 const Real dp = 1.0;
556 _fs->liquidProperties(p + dp, T, Xnacl, fsp);
557 Real rho1 = fsp[0].density.value();
558 Real mu1 = fsp[0].viscosity.value();
559 Real h1 = fsp[0].enthalpy.value();
560
561 _fs->liquidProperties(p - dp, T, Xnacl, fsp);
562 Real rho2 = fsp[0].density.value();
563 Real mu2 = fsp[0].viscosity.value();
564 Real h2 = fsp[0].enthalpy.value();
565
566 REL_TEST(liquid_density.derivatives()[_pidx], (rho1 - rho2) / (2.0 * dp), 1.0e-6);
567 REL_TEST(liquid_viscosity.derivatives()[_pidx], (mu1 - mu2) / (2.0 * dp), 2.0e-6);
568 REL_TEST(liquid_enthalpy.derivatives()[_pidx], (h1 - h2) / (2.0 * dp), 1.0e-6);
569
570 // Derivatives wrt temperature
571 const Real dT = 1.0e-4;
572 _fs->liquidProperties(p, T + dT, Xnacl, fsp);
573 rho1 = fsp[0].density.value();
574 mu1 = fsp[0].viscosity.value();
575 h1 = fsp[0].enthalpy.value();
576
577 _fs->liquidProperties(p, T - dT, Xnacl, fsp);
578 rho2 = fsp[0].density.value();
579 mu2 = fsp[0].viscosity.value();
580 h2 = fsp[0].enthalpy.value();
581
582 REL_TEST(liquid_density.derivatives()[_Tidx], (rho1 - rho2) / (2.0 * dT), 1.0e-6);
583 REL_TEST(liquid_viscosity.derivatives()[_Tidx], (mu1 - mu2) / (2.0 * dT), 1.0e-6);
584 REL_TEST(liquid_enthalpy.derivatives()[_Tidx], (h1 - h2) / (2.0 * dT), 1.0e-6);
585
586 // Derivatives wrt Xnacl
587 const Real dx = 1.0e-8;
588 _fs->liquidProperties(p, T, Xnacl + dx, fsp);
589 rho1 = fsp[0].density.value();
590 mu1 = fsp[0].viscosity.value();
591 h1 = fsp[0].enthalpy.value();
592
593 _fs->liquidProperties(p, T, Xnacl - dx, fsp);
594 rho2 = fsp[0].density.value();
595 mu2 = fsp[0].viscosity.value();
596 h2 = fsp[0].enthalpy.value();
597
598 REL_TEST(liquid_density.derivatives()[_Xidx], (rho1 - rho2) / (2.0 * dx), 1.0e-6);
599 REL_TEST(liquid_viscosity.derivatives()[_Xidx], (mu1 - mu2) / (2.0 * dx), 1.0e-6);
600 REL_TEST(liquid_enthalpy.derivatives()[_Xidx], (h1 - h2) / (2.0 * dx), 1.0e-6);
601
602 // Derivatives wrt Z
603 const Real dZ = 1.0e-8;
604 _fs->massFractions(p, T, Xnacl, Z + dZ, phase_state, fsp);
605 _fs->liquidProperties(p, T, Xnacl, fsp);
606 rho1 = fsp[0].density.value();
607 mu1 = fsp[0].viscosity.value();
608 h1 = fsp[0].enthalpy.value();
609
610 _fs->massFractions(p, T, Xnacl, Z - dZ, phase_state, fsp);
611 _fs->liquidProperties(p, T, Xnacl, fsp);
612 rho2 = fsp[0].density.value();
613 mu2 = fsp[0].viscosity.value();
614 h2 = fsp[0].enthalpy.value();
615
616 REL_TEST(liquid_density.derivatives()[_Zidx], (rho1 - rho2) / (2.0 * dZ), 1.0e-6);
617 ABS_TEST(liquid_viscosity.derivatives()[_Zidx], (mu1 - mu2) / (2.0 * dZ), 1.0e-6);
618 REL_TEST(liquid_enthalpy.derivatives()[_Zidx], (h1 - h2) / (2.0 * dZ), 1.0e-6);
619
620 // Two-phase region
621 Z = 0.045;
622 Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
623
624 _fs->massFractions(p, T, Xnacl, Z, phase_state, fsp);
625 EXPECT_EQ(phase_state, FluidStatePhaseEnum::TWOPHASE);
626
627 // Verify fluid density and viscosity derivatives
628 _fs->liquidProperties(p, T, Xnacl, fsp);
629
630 liquid_density = fsp[0].density;
631 liquid_viscosity = fsp[0].viscosity;
632 liquid_enthalpy = fsp[0].enthalpy;
633
634 // Derivatives wrt pressure
635 _fs->massFractions(p + dp, T, Xnacl, Z, phase_state, fsp);
636 _fs->liquidProperties(p + dp, T, Xnacl, fsp);
637 rho1 = fsp[0].density.value();
638 mu1 = fsp[0].viscosity.value();
639 h1 = fsp[0].enthalpy.value();
640
641 _fs->massFractions(p - dp, T, Xnacl, Z, phase_state, fsp);
642 _fs->liquidProperties(p - dp, T, Xnacl, fsp);
643 rho2 = fsp[0].density.value();
644 mu2 = fsp[0].viscosity.value();
645 h2 = fsp[0].enthalpy.value();
646
647 REL_TEST(liquid_density.derivatives()[_pidx], (rho1 - rho2) / (2.0 * dp), 1.0e-6);
648 REL_TEST(liquid_viscosity.derivatives()[_pidx], (mu1 - mu2) / (2.0 * dp), 2.0e-6);
649 REL_TEST(liquid_enthalpy.derivatives()[_pidx], (h1 - h2) / (2.0 * dp), 1.0e-6);
650
651 // Derivatives wrt temperature
652 _fs->massFractions(p, T + dT, Xnacl, Z, phase_state, fsp);
653 _fs->liquidProperties(p, T + dT, Xnacl, fsp);
654 rho1 = fsp[0].density.value();
655 mu1 = fsp[0].viscosity.value();
656 h1 = fsp[0].enthalpy.value();
657
658 _fs->massFractions(p, T - dT, Xnacl, Z, phase_state, fsp);
659 _fs->liquidProperties(p, T - dT, Xnacl, fsp);
660 rho2 = fsp[0].density.value();
661 mu2 = fsp[0].viscosity.value();
662 h2 = fsp[0].enthalpy.value();
663
664 REL_TEST(liquid_density.derivatives()[_Tidx], (rho1 - rho2) / (2.0 * dT), 1.0e-6);
665 REL_TEST(liquid_viscosity.derivatives()[_Tidx], (mu1 - mu2) / (2.0 * dT), 1.0e-6);
666 REL_TEST(liquid_enthalpy.derivatives()[_Tidx], (h1 - h2) / (2.0 * dT), 1.0e-6);
667
668 // Derivatives wrt Xnacl
669 _fs->massFractions(p, T, Xnacl + dx, Z, phase_state, fsp);
670 _fs->liquidProperties(p, T, Xnacl + dx, fsp);
671 rho1 = fsp[0].density.value();
672 mu1 = fsp[0].viscosity.value();
673 h1 = fsp[0].enthalpy.value();
674
675 _fs->massFractions(p, T, Xnacl - dx, Z, phase_state, fsp);
676 _fs->liquidProperties(p, T, Xnacl - dx, fsp);
677 rho2 = fsp[0].density.value();
678 mu2 = fsp[0].viscosity.value();
679 h2 = fsp[0].enthalpy.value();
680
681 REL_TEST(liquid_density.derivatives()[_Xidx], (rho1 - rho2) / (2.0 * dx), 1.0e-6);
682 REL_TEST(liquid_viscosity.derivatives()[_Xidx], (mu1 - mu2) / (2.0 * dx), 1.0e-6);
683 REL_TEST(liquid_enthalpy.derivatives()[_Xidx], (h1 - h2) / (2.0 * dx), 1.0e-6);
684
685 // Derivatives wrt Z
686 _fs->massFractions(p, T, Xnacl, Z + dZ, phase_state, fsp);
687 _fs->liquidProperties(p, T, Xnacl, fsp);
688 rho1 = fsp[0].density.value();
689 mu1 = fsp[0].viscosity.value();
690 h1 = fsp[0].enthalpy.value();
691
692 _fs->massFractions(p, T, Xnacl, Z - dZ, phase_state, fsp);
693 _fs->liquidProperties(p, T, Xnacl, fsp);
694 rho2 = fsp[0].density.value();
695 mu2 = fsp[0].viscosity.value();
696 h2 = fsp[0].enthalpy.value();
697
698 ABS_TEST(liquid_density.derivatives()[_Zidx], (rho1 - rho2) / (2.0 * dZ), 1.0e-6);
699 ABS_TEST(liquid_viscosity.derivatives()[_Zidx], (mu1 - mu2) / (2.0 * dZ), 1.0e-6);
700 ABS_TEST(liquid_enthalpy.derivatives()[_Zidx], (h1 - h2) / (2.0 * dZ), 1.0e-6);
701}
702
703/*
704 * Verify calculation of gas saturation and derivatives in the two-phase region
705 */
706TEST_F(PorousFlowBrineCO2Test, saturation)
707{
708 ADReal p = 1.0e6;
709 Moose::derivInsert(p.derivatives(), _pidx, 1.0);
710
711 ADReal T = 350.0;
712 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
713
714 ADReal Xnacl = 0.1;
715 Moose::derivInsert(Xnacl.derivatives(), _Xidx, 1.0);
716
717 FluidStatePhaseEnum phase_state;
718 const unsigned int np = _fs->numPhases();
719 const unsigned int nc = _fs->numComponents();
720 std::vector<FluidStateProperties> fsp(np, FluidStateProperties(nc));
721
722 // In the two-phase region, the mass fractions are the equilibrium values, so
723 // a temporary value of Z can be used (as long as it corresponds to the two-phase
724 // region)
725 ADReal Z = 0.45;
726 Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
727
728 _fs->massFractions(p, T, Xnacl, Z, phase_state, fsp);
729 EXPECT_EQ(phase_state, FluidStatePhaseEnum::TWOPHASE);
730
731 // Calculate Z that gives a saturation of 0.25
732 ADReal gas_saturation = 0.25;
733 ADReal liquid_pressure = p - _pc->capillaryPressure(1.0 - gas_saturation);
734
735 // Calculate gas density and liquid density
736 _fs->gasProperties(p, T, fsp);
737 _fs->liquidProperties(liquid_pressure, T, Xnacl, fsp);
738
739 // The mass fraction that corresponds to a gas_saturation = 0.25
740 ADReal Zc = (gas_saturation * fsp[1].density * fsp[1].mass_fraction[1] +
741 (1.0 - gas_saturation) * fsp[0].density * fsp[0].mass_fraction[1]) /
742 (gas_saturation * fsp[1].density + (1.0 - gas_saturation) * fsp[0].density);
743
744 // Calculate the gas saturation and derivatives
745 ADReal saturation = _fs->saturation(p, T, Xnacl, Zc, fsp);
746 ABS_TEST(saturation.value(), gas_saturation.value(), 1.0e-6);
747
748 // Test the derivatives of gas saturation
749 gas_saturation = _fs->saturation(p, T, Xnacl, Z, fsp);
750
751 const Real dp = 1.0e-1;
752
753 // Derivative wrt pressure
754 _fs->massFractions(p + dp, T, Xnacl, Z, phase_state, fsp);
755 Real gsat1 = _fs->saturation(p + dp, T, Xnacl, Z, fsp).value();
756
757 _fs->massFractions(p - dp, T, Xnacl, Z, phase_state, fsp);
758 Real gsat2 = _fs->saturation(p - dp, T, Xnacl, Z, fsp).value();
759
760 REL_TEST(gas_saturation.derivatives()[_pidx], (gsat1 - gsat2) / (2.0 * dp), 1.0e-6);
761
762 // Derivative wrt temperature
763 const Real dT = 1.0e-4;
764 _fs->massFractions(p, T + dT, Xnacl, Z, phase_state, fsp);
765 gsat1 = _fs->saturation(p, T + dT, Xnacl, Z, fsp).value();
766
767 _fs->massFractions(p, T - dT, Xnacl, Z, phase_state, fsp);
768 gsat2 = _fs->saturation(p, T - dT, Xnacl, Z, fsp).value();
769
770 REL_TEST(gas_saturation.derivatives()[_Tidx], (gsat1 - gsat2) / (2.0 * dT), 1.0e-6);
771
772 // Derivative wrt Xnacl
773 const Real dx = 1.0e-8;
774 _fs->massFractions(p, T, Xnacl + dx, Z, phase_state, fsp);
775 gsat1 = _fs->saturation(p, T, Xnacl + dx, Z, fsp).value();
776
777 _fs->massFractions(p, T, Xnacl - dx, Z, phase_state, fsp);
778 gsat2 = _fs->saturation(p, T, Xnacl - dx, Z, fsp).value();
779
780 REL_TEST(gas_saturation.derivatives()[_Xidx], (gsat1 - gsat2) / (2.0 * dx), 1.0e-6);
781
782 // Derivative wrt Z
783 const Real dZ = 1.0e-8;
784
785 _fs->massFractions(p, T, Xnacl, Z, phase_state, fsp);
786 gsat1 = _fs->saturation(p, T, Xnacl, Z + dZ, fsp).value();
787
788 gsat2 = _fs->saturation(p, T, Xnacl, Z - dZ, fsp).value();
789
790 REL_TEST(gas_saturation.derivatives()[_Zidx], (gsat1 - gsat2) / (2.0 * dZ), 1.0e-6);
791}
792
793/*
794 * Verify calculation of total mass fraction given a gas saturation
795 */
796TEST_F(PorousFlowBrineCO2Test, totalMassFraction)
797{
798 const Real p = 1.0e6;
799 const Real T = 350.0;
800 const Real Xnacl = 0.1;
801 const Real s = 0.2;
802 const unsigned qp = 0;
803
804 Real Z = _fs->totalMassFraction(p, T, Xnacl, s, qp);
805
806 // Test that the saturation calculated in this fluid state using Z is equal to s
807 FluidStatePhaseEnum phase_state;
808 const unsigned int np = _fs->numPhases();
809 const unsigned int nc = _fs->numComponents();
810 std::vector<FluidStateProperties> fsp(np, FluidStateProperties(nc));
811
812 ADReal pressure = p;
813 ADReal temperature = T;
814 ADReal z = Z;
815 ADReal x = Xnacl;
816 _fs->massFractions(pressure, temperature, x, z, phase_state, fsp);
817 EXPECT_EQ(phase_state, FluidStatePhaseEnum::TWOPHASE);
818
819 ADReal gas_saturation = _fs->saturation(pressure, temperature, x, z, fsp);
820 ABS_TEST(gas_saturation, s, 1.0e-6);
821}
822
823/*
824 * Verify calculation of Henry constant with brine correction. Note: the values
825 * calculated compare well by eye to the values presented in Figure 4 of
826 * Battistelli et al, "A fluid property module for the TOUGH2 simulator for saline
827 * brines with non-condensible gas"
828 */
829TEST_F(PorousFlowBrineCO2Test, henryConstant)
830{
831 ADReal T = 373.15;
832 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
833
834 ADReal Xnacl = 0.1;
835 Moose::derivInsert(Xnacl.derivatives(), _Xidx, 1.0);
836
837 ADReal Kh = _fs->henryConstant(T, Xnacl);
838 REL_TEST(Kh.value(), 7.46559e+08, 1.0e-3);
839
840 T = 473.15;
841 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
842
843 Xnacl = 0.2;
844 Moose::derivInsert(Xnacl.derivatives(), _Xidx, 1.0);
845
846 Kh = _fs->henryConstant(T, Xnacl);
847 REL_TEST(Kh.value(), 1.66069e+09, 1.0e-3);
848
849 // Test the derivative wrt temperature
850 const Real dT = 1.0e-4;
851 ADReal Kh1 = _fs->henryConstant(T + dT, Xnacl);
852 ADReal Kh2 = _fs->henryConstant(T - dT, Xnacl);
853
854 REL_TEST(Kh.derivatives()[_Tidx], (Kh1.value() - Kh2.value()) / (2.0 * dT), 1.0e-6);
855
856 // Test the derivative wrt Xnacl
857 const Real dx = 1.0e-8;
858 Kh1 = _fs->henryConstant(T, Xnacl + dx);
859 Kh2 = _fs->henryConstant(T, Xnacl - dx);
860
861 REL_TEST(Kh.derivatives()[_Xidx], (Kh1.value() - Kh2.value()) / (2.0 * dx), 1.0e-6);
862}
863
864/*
865 * Verify calculation of enthalpy of dissolution. Note: the values calculated compare
866 * well by eye to the values presented in Figure 4 of Battistelli et al, "A fluid property
867 * module for the TOUGH2 simulator for saline brines with non-condensible gas"
868 */
869TEST_F(PorousFlowBrineCO2Test, enthalpyOfDissolutionGas)
870{
871 // T = 50C
872 ADReal T = 323.15;
873 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
874
875 ADReal Xnacl = 0.1;
876 Moose::derivInsert(Xnacl.derivatives(), _Xidx, 1.0);
877
878 // Enthalpy of dissolution of CO2 in brine
879 ADReal hdis = _fs->enthalpyOfDissolutionGas(T, Xnacl);
880 REL_TEST(hdis, -3.20130e5, 1.0e-3);
881
882 // T = 350C
883 T = 623.15;
884 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
885
886 // Enthalpy of dissolution of CO2 in brine
887 hdis = _fs->enthalpyOfDissolutionGas(T, Xnacl);
888 REL_TEST(hdis, 9.83813e+05, 1.0e-3);
889
890 // Test the derivative wrt temperature
891 const Real dT = 1.0e-4;
892 Real hdis1 = _fs->enthalpyOfDissolutionGas(T + dT, Xnacl).value();
893 Real hdis2 = _fs->enthalpyOfDissolutionGas(T - dT, Xnacl).value();
894
895 REL_TEST(hdis.derivatives()[_Tidx], (hdis1 - hdis2) / (2.0 * dT), 1.0e-5);
896
897 // Test the derivative wrt salt mass fraction
898 const Real dx = 1.0e-8;
899 hdis1 = _fs->enthalpyOfDissolutionGas(T, Xnacl + dx).value();
900 hdis2 = _fs->enthalpyOfDissolutionGas(T, Xnacl - dx).value();
901
902 REL_TEST(hdis.derivatives()[_Xidx], (hdis1 - hdis2) / (2.0 * dx), 1.0e-5);
903}
904
905/*
906 * Verify calculation of enthalpy of dissolution of CO2. Note: the values calculated compare
907 * well by eye to the values presented in Figure 6 of Duan and Sun, An improved model
908 * calculating CO2 solubility in pure water and aqueous NaCl solutions from 273 to 533 K
909 * and from 0 to 2000 bar, Chemical geology, 193, 257--271 (2003)
910 */
911TEST_F(PorousFlowBrineCO2Test, enthalpyOfDissolution)
912{
913 // T = 50C
914 ADReal T = 323.15;
915 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
916
917 // Enthalpy of dissolution of CO2 in water
918 ADReal hdis = _fs->enthalpyOfDissolution(T);
919 REL_TEST(hdis.value(), -3.38185e5, 1.0e-3);
920
921 // T = 350C
922 T = 623.15;
923 Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
924
925 // Enthalpy of dissolution of CO2 in water
926 hdis = _fs->enthalpyOfDissolution(T);
927 REL_TEST(hdis.value(), 5.78787e5, 1.0e-3);
928
929 // Test the derivative wrt temperature
930 const Real dT = 1.0e-4;
931 Real hdis1 = _fs->enthalpyOfDissolution(T + dT).value();
932 Real hdis2 = _fs->enthalpyOfDissolution(T - dT).value();
933
934 REL_TEST(hdis.derivatives()[_Tidx], (hdis1 - hdis2) / (2.0 * dT), 1.0e-5);
935}
936
941TEST_F(PorousFlowBrineCO2Test, solveEquilibriumMoleFractionHighTemp)
942{
943 Real p = 20.0e6;
944 Real T = 473.15;
945 Real Xnacl = 0.0;
946
947 Real yh2o, xco2;
948 Real co2_density = _co2_fp->rho_from_p_T(p, T);
949 _fs->solveEquilibriumMoleFractionHighTemp(p, T, Xnacl, co2_density, xco2, yh2o);
950 ABS_TEST(yh2o, 0.161429743509, 1.0e-10);
951 ABS_TEST(xco2, 0.0236966821858, 1.0e-10);
952
953 // Mass fraction equivalent to molality of 2
954 p = 30.0e6;
955 T = 423.15;
956 Xnacl = 0.105;
957
958 _fs->solveEquilibriumMoleFractionHighTemp(p, T, Xnacl, co2_density, xco2, yh2o);
959 ABS_TEST(yh2o, 0.0468197608955, 1.0e-10);
960 ABS_TEST(xco2, 0.0236644599437, 1.0e-10);
961
962 // Mass fraction equivalent to molality of 4
963 T = 523.15;
964 Xnacl = 0.189;
965
966 co2_density = _co2_fp->rho_from_p_T(p, T);
967 _fs->solveEquilibriumMoleFractionHighTemp(p, T, Xnacl, co2_density, xco2, yh2o);
968 ABS_TEST(yh2o, 0.253292728991, 1.0e-10);
969 ABS_TEST(xco2, 0.0168344513321, 1.0e-10);
970}
971
975TEST_F(PorousFlowBrineCO2Test, equilibriumMoleFractions)
976{
977 // Test pure water (Xnacl = 0)
978 // Low temperature regime
979 ADReal p = 20.0e6;
980 ADReal T = 323.15;
981 ADReal Xnacl = 0.0;
982
983 ADReal x, y;
984 _fs->equilibriumMoleFractions(p, T, Xnacl, x, y);
985 ABS_TEST(y.value(), 0.00696393845155, 1.0e-8);
986 ABS_TEST(x.value(), 0.0236554537395, 1.0e-8);
987
988 // Intermediate temperature regime
989 T = 373.15;
990
991 _fs->equilibriumMoleFractions(p, T, Xnacl, x, y);
992 ABS_TEST(y.value(), 0.0194394631944, 1.0e-8);
993 ABS_TEST(x.value(), 0.020195776139, 1.0e-8);
994
995 // High temperature regime
996 p = 30.0e6;
997 T = 523.15;
998
999 _fs->equilibriumMoleFractions(p, T, Xnacl, x, y);
1000 ABS_TEST(y.value(), 0.286117195565, 1.0e-8);
1001 ABS_TEST(x.value(), 0.0409622051253, 1.0e-8);
1002
1003 // High temperature regime with low pressure
1004 p = 1.0e6;
1005 T = 523.15;
1006
1007 _fs->equilibriumMoleFractions(p, T, Xnacl, x, y);
1008 ABS_TEST(y.value(), 1.0, 1.0e-8);
1009 ABS_TEST(x.value(), 0.0, 1.0e-8);
1010
1011 // Test brine (Xnacl = 0.1)
1012 // Low temperature regime
1013 p = 20.0e6;
1014 T = 323.15;
1015 Xnacl = 0.1;
1016
1017 _fs->equilibriumMoleFractions(p, T, Xnacl, x, y);
1018 ABS_TEST(y.value(), 0.00657335846826, 1.0e-8);
1019 ABS_TEST(x.value(), 0.0152863933134, 1.0e-8);
1020
1021 // Intermediate temperature regime
1022 T = 373.15;
1023
1024 _fs->equilibriumMoleFractions(p, T, Xnacl, x, y);
1025 ABS_TEST(y.value(), 0.01831360857, 1.0e-8);
1026 ABS_TEST(x.value(), 0.0132653916293, 1.0e-8);
1027
1028 // High temperature regime
1029 p = 30.0e6;
1030 T = 523.15;
1031
1032 _fs->equilibriumMoleFractions(p, T, Xnacl, x, y);
1033 ABS_TEST(y.value(), 0.270258924237, 1.0e-8);
1034 ABS_TEST(x.value(), 0.0246589135776, 1.0e-8);
1035}
ADReal
DualNumber< Real, DNDerivativeType, true > ADReal
y
const std::vector< double > y
Definition EquilibriumConstantFitTest.C:18
x
const std::vector< double > x
Definition EquilibriumConstantFitTest.C:17
p
const Real p
Definition GasLiquidMassTransferTest.C:18
T
const double T
Definition GasLiquidMassTransferTest.C:19
TEST_F
TEST_F(PorousFlowBrineCO2Test, name)
Verify that the correct name is supplied.
Definition PorousFlowBrineCO2Test.C:16
FluidStatePhaseEnum
FluidStatePhaseEnum
Phase state enum.
Definition PorousFlowFluidStateBase.h:17
name
const std::string name
Definition Setup.h:21
Moose::derivInsert
void derivInsert(SemiDynamicSparseNumberArray< Real, libMesh::dof_id_type, NWrapper< N > > &derivs, libMesh::dof_id_type index, Real value)
FluidStateProperties
AD data structure to pass calculated thermophysical properties.
Definition PorousFlowFluidStateBase.h:25
Generated on Fri Aug 21 2026 13:41:31 for https://mooseframework.inl.gov by doxygen 1.9.8