https://mooseframework.inl.gov
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 
16 TEST_F(PorousFlowBrineCO2Test, name) { EXPECT_EQ("brine-co2", _fs->fluidStateName()); }
17 
21 TEST_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  */
35 TEST_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  */
78 TEST_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  */
122 TEST_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  */
146 TEST_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  */
209 TEST_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  */
228 TEST_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  */
298 TEST_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  */
416 TEST_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  */
506 TEST_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  */
706 TEST_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  */
796 TEST_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  */
829 TEST_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  */
869 TEST_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  */
911 TEST_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 
941 TEST_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 
975 TEST_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 }
TEST_F
TEST_F(PorousFlowBrineCO2Test, name)
Verify that the correct name is supplied.
Definition: PorousFlowBrineCO2Test.C:16
NS::density
static const std::string density
Definition: NS.h:33
y
const std::vector< double > y
Definition: EquilibriumConstantFitTest.C:18
NS::temperature
static const std::string temperature
Definition: NS.h:59
ADReal
DualNumber< Real, DNDerivativeType, true > ADReal
value
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
x
const std::vector< double > x
Definition: EquilibriumConstantFitTest.C:17
name
const std::string name
Definition: Setup.h:20
FluidStateProperties
AD data structure to pass calculated thermophysical properties.
Definition: PorousFlowFluidStateBase.h:24
FluidStatePhaseEnum
FluidStatePhaseEnum
Phase state enum.
Definition: PorousFlowFluidStateBase.h:16
NS::Z
static const std::string Z
Definition: NS.h:169
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
NS::pressure
static const std::string pressure
Definition: NS.h:56
Moose::derivInsert
void derivInsert(SemiDynamicSparseNumberArray< Real, libMesh::dof_id_type, NWrapper< N >> &derivs, libMesh::dof_id_type index, Real value)
Generated on Fri Aug 29 2025 20:18:28 for https://mooseframework.inl.gov by   doxygen 1.8.14