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  ABS_TEST(liquid_enthalpy.value(), enthalpy, 1.0e-12);
551 
552  // Verify derivatives
553  // Derivatives wrt pressure
554  const Real dp = 1.0;
555  _fs->liquidProperties(p + dp, T, Xnacl, fsp);
556  Real rho1 = fsp[0].density.value();
557  Real mu1 = fsp[0].viscosity.value();
558  Real h1 = fsp[0].enthalpy.value();
559 
560  _fs->liquidProperties(p - dp, T, Xnacl, fsp);
561  Real rho2 = fsp[0].density.value();
562  Real mu2 = fsp[0].viscosity.value();
563  Real h2 = fsp[0].enthalpy.value();
564 
565  REL_TEST(liquid_density.derivatives()[_pidx], (rho1 - rho2) / (2.0 * dp), 1.0e-6);
566  REL_TEST(liquid_viscosity.derivatives()[_pidx], (mu1 - mu2) / (2.0 * dp), 2.0e-6);
567  REL_TEST(liquid_enthalpy.derivatives()[_pidx], (h1 - h2) / (2.0 * dp), 1.0e-6);
568 
569  // Derivatives wrt temperature
570  const Real dT = 1.0e-4;
571  _fs->liquidProperties(p, T + dT, Xnacl, fsp);
572  rho1 = fsp[0].density.value();
573  mu1 = fsp[0].viscosity.value();
574  h1 = fsp[0].enthalpy.value();
575 
576  _fs->liquidProperties(p, T - dT, Xnacl, fsp);
577  rho2 = fsp[0].density.value();
578  mu2 = fsp[0].viscosity.value();
579  h2 = fsp[0].enthalpy.value();
580 
581  REL_TEST(liquid_density.derivatives()[_Tidx], (rho1 - rho2) / (2.0 * dT), 1.0e-6);
582  REL_TEST(liquid_viscosity.derivatives()[_Tidx], (mu1 - mu2) / (2.0 * dT), 1.0e-6);
583  REL_TEST(liquid_enthalpy.derivatives()[_Tidx], (h1 - h2) / (2.0 * dT), 1.0e-6);
584 
585  // Derivatives wrt Xnacl
586  const Real dx = 1.0e-8;
587  _fs->liquidProperties(p, T, Xnacl + dx, fsp);
588  rho1 = fsp[0].density.value();
589  mu1 = fsp[0].viscosity.value();
590  h1 = fsp[0].enthalpy.value();
591 
592  _fs->liquidProperties(p, T, Xnacl - dx, fsp);
593  rho2 = fsp[0].density.value();
594  mu2 = fsp[0].viscosity.value();
595  h2 = fsp[0].enthalpy.value();
596 
597  REL_TEST(liquid_density.derivatives()[_Xidx], (rho1 - rho2) / (2.0 * dx), 1.0e-6);
598  REL_TEST(liquid_viscosity.derivatives()[_Xidx], (mu1 - mu2) / (2.0 * dx), 1.0e-6);
599  REL_TEST(liquid_enthalpy.derivatives()[_Xidx], (h1 - h2) / (2.0 * dx), 1.0e-6);
600 
601  // Derivatives wrt Z
602  const Real dZ = 1.0e-8;
603  _fs->massFractions(p, T, Xnacl, Z + dZ, phase_state, fsp);
604  _fs->liquidProperties(p, T, Xnacl, fsp);
605  rho1 = fsp[0].density.value();
606  mu1 = fsp[0].viscosity.value();
607  h1 = fsp[0].enthalpy.value();
608 
609  _fs->massFractions(p, T, Xnacl, Z - dZ, phase_state, fsp);
610  _fs->liquidProperties(p, T, Xnacl, fsp);
611  rho2 = fsp[0].density.value();
612  mu2 = fsp[0].viscosity.value();
613  h2 = fsp[0].enthalpy.value();
614 
615  REL_TEST(liquid_density.derivatives()[_Zidx], (rho1 - rho2) / (2.0 * dZ), 1.0e-6);
616  ABS_TEST(liquid_viscosity.derivatives()[_Zidx], (mu1 - mu2) / (2.0 * dZ), 1.0e-6);
617  REL_TEST(liquid_enthalpy.derivatives()[_Zidx], (h1 - h2) / (2.0 * dZ), 1.0e-6);
618 
619  // Two-phase region
620  Z = 0.045;
621  Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
622 
623  _fs->massFractions(p, T, Xnacl, Z, phase_state, fsp);
624  EXPECT_EQ(phase_state, FluidStatePhaseEnum::TWOPHASE);
625 
626  // Verify fluid density and viscosity derivatives
627  _fs->liquidProperties(p, T, Xnacl, fsp);
628 
629  liquid_density = fsp[0].density;
630  liquid_viscosity = fsp[0].viscosity;
631  liquid_enthalpy = fsp[0].enthalpy;
632 
633  // Derivatives wrt pressure
634  _fs->massFractions(p + dp, T, Xnacl, Z, phase_state, fsp);
635  _fs->liquidProperties(p + dp, T, Xnacl, fsp);
636  rho1 = fsp[0].density.value();
637  mu1 = fsp[0].viscosity.value();
638  h1 = fsp[0].enthalpy.value();
639 
640  _fs->massFractions(p - dp, T, Xnacl, Z, phase_state, fsp);
641  _fs->liquidProperties(p - dp, T, Xnacl, fsp);
642  rho2 = fsp[0].density.value();
643  mu2 = fsp[0].viscosity.value();
644  h2 = fsp[0].enthalpy.value();
645 
646  REL_TEST(liquid_density.derivatives()[_pidx], (rho1 - rho2) / (2.0 * dp), 1.0e-6);
647  REL_TEST(liquid_viscosity.derivatives()[_pidx], (mu1 - mu2) / (2.0 * dp), 2.0e-6);
648  REL_TEST(liquid_enthalpy.derivatives()[_pidx], (h1 - h2) / (2.0 * dp), 1.0e-6);
649 
650  // Derivatives wrt temperature
651  _fs->massFractions(p, T + dT, Xnacl, Z, phase_state, fsp);
652  _fs->liquidProperties(p, T + dT, Xnacl, fsp);
653  rho1 = fsp[0].density.value();
654  mu1 = fsp[0].viscosity.value();
655  h1 = fsp[0].enthalpy.value();
656 
657  _fs->massFractions(p, T - dT, Xnacl, Z, phase_state, fsp);
658  _fs->liquidProperties(p, T - dT, Xnacl, fsp);
659  rho2 = fsp[0].density.value();
660  mu2 = fsp[0].viscosity.value();
661  h2 = fsp[0].enthalpy.value();
662 
663  REL_TEST(liquid_density.derivatives()[_Tidx], (rho1 - rho2) / (2.0 * dT), 1.0e-6);
664  REL_TEST(liquid_viscosity.derivatives()[_Tidx], (mu1 - mu2) / (2.0 * dT), 1.0e-6);
665  REL_TEST(liquid_enthalpy.derivatives()[_Tidx], (h1 - h2) / (2.0 * dT), 1.0e-6);
666 
667  // Derivatives wrt Xnacl
668  _fs->massFractions(p, T, Xnacl + dx, Z, phase_state, fsp);
669  _fs->liquidProperties(p, T, Xnacl + dx, fsp);
670  rho1 = fsp[0].density.value();
671  mu1 = fsp[0].viscosity.value();
672  h1 = fsp[0].enthalpy.value();
673 
674  _fs->massFractions(p, T, Xnacl - dx, Z, phase_state, fsp);
675  _fs->liquidProperties(p, T, Xnacl - dx, fsp);
676  rho2 = fsp[0].density.value();
677  mu2 = fsp[0].viscosity.value();
678  h2 = fsp[0].enthalpy.value();
679 
680  REL_TEST(liquid_density.derivatives()[_Xidx], (rho1 - rho2) / (2.0 * dx), 1.0e-6);
681  REL_TEST(liquid_viscosity.derivatives()[_Xidx], (mu1 - mu2) / (2.0 * dx), 1.0e-6);
682  REL_TEST(liquid_enthalpy.derivatives()[_Xidx], (h1 - h2) / (2.0 * dx), 1.0e-6);
683 
684  // Derivatives wrt Z
685  _fs->massFractions(p, T, Xnacl, Z + dZ, phase_state, fsp);
686  _fs->liquidProperties(p, T, Xnacl, fsp);
687  rho1 = fsp[0].density.value();
688  mu1 = fsp[0].viscosity.value();
689  h1 = fsp[0].enthalpy.value();
690 
691  _fs->massFractions(p, T, Xnacl, Z - dZ, phase_state, fsp);
692  _fs->liquidProperties(p, T, Xnacl, fsp);
693  rho2 = fsp[0].density.value();
694  mu2 = fsp[0].viscosity.value();
695  h2 = fsp[0].enthalpy.value();
696 
697  ABS_TEST(liquid_density.derivatives()[_Zidx], (rho1 - rho2) / (2.0 * dZ), 1.0e-6);
698  ABS_TEST(liquid_viscosity.derivatives()[_Zidx], (mu1 - mu2) / (2.0 * dZ), 1.0e-6);
699  ABS_TEST(liquid_enthalpy.derivatives()[_Zidx], (h1 - h2) / (2.0 * dZ), 1.0e-6);
700 }
701 
702 /*
703  * Verify calculation of gas saturation and derivatives in the two-phase region
704  */
705 TEST_F(PorousFlowBrineCO2Test, saturation)
706 {
707  ADReal p = 1.0e6;
708  Moose::derivInsert(p.derivatives(), _pidx, 1.0);
709 
710  ADReal T = 350.0;
711  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
712 
713  ADReal Xnacl = 0.1;
714  Moose::derivInsert(Xnacl.derivatives(), _Xidx, 1.0);
715 
716  FluidStatePhaseEnum phase_state;
717  const unsigned int np = _fs->numPhases();
718  const unsigned int nc = _fs->numComponents();
719  std::vector<FluidStateProperties> fsp(np, FluidStateProperties(nc));
720 
721  // In the two-phase region, the mass fractions are the equilibrium values, so
722  // a temporary value of Z can be used (as long as it corresponds to the two-phase
723  // region)
724  ADReal Z = 0.45;
725  Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
726 
727  _fs->massFractions(p, T, Xnacl, Z, phase_state, fsp);
728  EXPECT_EQ(phase_state, FluidStatePhaseEnum::TWOPHASE);
729 
730  // Calculate Z that gives a saturation of 0.25
731  ADReal gas_saturation = 0.25;
732  ADReal liquid_pressure = p - _pc->capillaryPressure(1.0 - gas_saturation);
733 
734  // Calculate gas density and liquid density
735  _fs->gasProperties(p, T, fsp);
736  _fs->liquidProperties(liquid_pressure, T, Xnacl, fsp);
737 
738  // The mass fraction that corresponds to a gas_saturation = 0.25
739  ADReal Zc = (gas_saturation * fsp[1].density * fsp[1].mass_fraction[1] +
740  (1.0 - gas_saturation) * fsp[0].density * fsp[0].mass_fraction[1]) /
741  (gas_saturation * fsp[1].density + (1.0 - gas_saturation) * fsp[0].density);
742 
743  // Calculate the gas saturation and derivatives
744  ADReal saturation = _fs->saturation(p, T, Xnacl, Zc, fsp);
745  ABS_TEST(saturation.value(), gas_saturation.value(), 1.0e-6);
746 
747  // Test the derivatives of gas saturation
748  gas_saturation = _fs->saturation(p, T, Xnacl, Z, fsp);
749 
750  const Real dp = 1.0e-1;
751 
752  // Derivative wrt pressure
753  _fs->massFractions(p + dp, T, Xnacl, Z, phase_state, fsp);
754  Real gsat1 = _fs->saturation(p + dp, T, Xnacl, Z, fsp).value();
755 
756  _fs->massFractions(p - dp, T, Xnacl, Z, phase_state, fsp);
757  Real gsat2 = _fs->saturation(p - dp, T, Xnacl, Z, fsp).value();
758 
759  REL_TEST(gas_saturation.derivatives()[_pidx], (gsat1 - gsat2) / (2.0 * dp), 1.0e-6);
760 
761  // Derivative wrt temperature
762  const Real dT = 1.0e-4;
763  _fs->massFractions(p, T + dT, Xnacl, Z, phase_state, fsp);
764  gsat1 = _fs->saturation(p, T + dT, Xnacl, Z, fsp).value();
765 
766  _fs->massFractions(p, T - dT, Xnacl, Z, phase_state, fsp);
767  gsat2 = _fs->saturation(p, T - dT, Xnacl, Z, fsp).value();
768 
769  REL_TEST(gas_saturation.derivatives()[_Tidx], (gsat1 - gsat2) / (2.0 * dT), 1.0e-6);
770 
771  // Derivative wrt Xnacl
772  const Real dx = 1.0e-8;
773  _fs->massFractions(p, T, Xnacl + dx, Z, phase_state, fsp);
774  gsat1 = _fs->saturation(p, T, Xnacl + dx, Z, fsp).value();
775 
776  _fs->massFractions(p, T, Xnacl - dx, Z, phase_state, fsp);
777  gsat2 = _fs->saturation(p, T, Xnacl - dx, Z, fsp).value();
778 
779  REL_TEST(gas_saturation.derivatives()[_Xidx], (gsat1 - gsat2) / (2.0 * dx), 1.0e-6);
780 
781  // Derivative wrt Z
782  const Real dZ = 1.0e-8;
783 
784  _fs->massFractions(p, T, Xnacl, Z, phase_state, fsp);
785  gsat1 = _fs->saturation(p, T, Xnacl, Z + dZ, fsp).value();
786 
787  gsat2 = _fs->saturation(p, T, Xnacl, Z - dZ, fsp).value();
788 
789  REL_TEST(gas_saturation.derivatives()[_Zidx], (gsat1 - gsat2) / (2.0 * dZ), 1.0e-6);
790 }
791 
792 /*
793  * Verify calculation of total mass fraction given a gas saturation
794  */
795 TEST_F(PorousFlowBrineCO2Test, totalMassFraction)
796 {
797  const Real p = 1.0e6;
798  const Real T = 350.0;
799  const Real Xnacl = 0.1;
800  const Real s = 0.2;
801  const unsigned qp = 0;
802 
803  Real Z = _fs->totalMassFraction(p, T, Xnacl, s, qp);
804 
805  // Test that the saturation calculated in this fluid state using Z is equal to s
806  FluidStatePhaseEnum phase_state;
807  const unsigned int np = _fs->numPhases();
808  const unsigned int nc = _fs->numComponents();
809  std::vector<FluidStateProperties> fsp(np, FluidStateProperties(nc));
810 
811  ADReal pressure = p;
812  ADReal temperature = T;
813  ADReal z = Z;
814  ADReal x = Xnacl;
815  _fs->massFractions(pressure, temperature, x, z, phase_state, fsp);
816  EXPECT_EQ(phase_state, FluidStatePhaseEnum::TWOPHASE);
817 
818  ADReal gas_saturation = _fs->saturation(pressure, temperature, x, z, fsp);
819  ABS_TEST(gas_saturation, s, 1.0e-6);
820 }
821 
822 /*
823  * Verify calculation of Henry constant with brine correction. Note: the values
824  * calculated compare well by eye to the values presented in Figure 4 of
825  * Battistelli et al, "A fluid property module for the TOUGH2 simulator for saline
826  * brines with non-condensible gas"
827  */
828 TEST_F(PorousFlowBrineCO2Test, henryConstant)
829 {
830  ADReal T = 373.15;
831  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
832 
833  ADReal Xnacl = 0.1;
834  Moose::derivInsert(Xnacl.derivatives(), _Xidx, 1.0);
835 
836  ADReal Kh = _fs->henryConstant(T, Xnacl);
837  REL_TEST(Kh.value(), 7.46559e+08, 1.0e-3);
838 
839  T = 473.15;
840  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
841 
842  Xnacl = 0.2;
843  Moose::derivInsert(Xnacl.derivatives(), _Xidx, 1.0);
844 
845  Kh = _fs->henryConstant(T, Xnacl);
846  REL_TEST(Kh.value(), 1.66069e+09, 1.0e-3);
847 
848  // Test the derivative wrt temperature
849  const Real dT = 1.0e-4;
850  ADReal Kh1 = _fs->henryConstant(T + dT, Xnacl);
851  ADReal Kh2 = _fs->henryConstant(T - dT, Xnacl);
852 
853  REL_TEST(Kh.derivatives()[_Tidx], (Kh1.value() - Kh2.value()) / (2.0 * dT), 1.0e-6);
854 
855  // Test the derivative wrt Xnacl
856  const Real dx = 1.0e-8;
857  Kh1 = _fs->henryConstant(T, Xnacl + dx);
858  Kh2 = _fs->henryConstant(T, Xnacl - dx);
859 
860  REL_TEST(Kh.derivatives()[_Xidx], (Kh1.value() - Kh2.value()) / (2.0 * dx), 1.0e-6);
861 }
862 
863 /*
864  * Verify calculation of enthalpy of dissolution. Note: the values calculated compare
865  * well by eye to the values presented in Figure 4 of Battistelli et al, "A fluid property
866  * module for the TOUGH2 simulator for saline brines with non-condensible gas"
867  */
868 TEST_F(PorousFlowBrineCO2Test, enthalpyOfDissolutionGas)
869 {
870  // T = 50C
871  ADReal T = 323.15;
872  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
873 
874  ADReal Xnacl = 0.1;
875  Moose::derivInsert(Xnacl.derivatives(), _Xidx, 1.0);
876 
877  // Enthalpy of dissolution of CO2 in brine
878  ADReal hdis = _fs->enthalpyOfDissolutionGas(T, Xnacl);
879  REL_TEST(hdis, -3.20130e5, 1.0e-3);
880 
881  // T = 350C
882  T = 623.15;
883  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
884 
885  // Enthalpy of dissolution of CO2 in brine
886  hdis = _fs->enthalpyOfDissolutionGas(T, Xnacl);
887  REL_TEST(hdis, 9.83813e+05, 1.0e-3);
888 
889  // Test the derivative wrt temperature
890  const Real dT = 1.0e-4;
891  Real hdis1 = _fs->enthalpyOfDissolutionGas(T + dT, Xnacl).value();
892  Real hdis2 = _fs->enthalpyOfDissolutionGas(T - dT, Xnacl).value();
893 
894  REL_TEST(hdis.derivatives()[_Tidx], (hdis1 - hdis2) / (2.0 * dT), 1.0e-5);
895 
896  // Test the derivative wrt salt mass fraction
897  const Real dx = 1.0e-8;
898  hdis1 = _fs->enthalpyOfDissolutionGas(T, Xnacl + dx).value();
899  hdis2 = _fs->enthalpyOfDissolutionGas(T, Xnacl - dx).value();
900 
901  REL_TEST(hdis.derivatives()[_Xidx], (hdis1 - hdis2) / (2.0 * dx), 1.0e-5);
902 }
903 
904 /*
905  * Verify calculation of enthalpy of dissolution of CO2. Note: the values calculated compare
906  * well by eye to the values presented in Figure 6 of Duan and Sun, An improved model
907  * calculating CO2 solubility in pure water and aqueous NaCl solutions from 273 to 533 K
908  * and from 0 to 2000 bar, Chemical geology, 193, 257--271 (2003)
909  */
910 TEST_F(PorousFlowBrineCO2Test, enthalpyOfDissolution)
911 {
912  // T = 50C
913  ADReal T = 323.15;
914  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
915 
916  // Enthalpy of dissolution of CO2 in water
917  ADReal hdis = _fs->enthalpyOfDissolution(T);
918  REL_TEST(hdis.value(), -3.38185e5, 1.0e-3);
919 
920  // T = 350C
921  T = 623.15;
922  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
923 
924  // Enthalpy of dissolution of CO2 in water
925  hdis = _fs->enthalpyOfDissolution(T);
926  REL_TEST(hdis.value(), 5.78787e5, 1.0e-3);
927 
928  // Test the derivative wrt temperature
929  const Real dT = 1.0e-4;
930  Real hdis1 = _fs->enthalpyOfDissolution(T + dT).value();
931  Real hdis2 = _fs->enthalpyOfDissolution(T - dT).value();
932 
933  REL_TEST(hdis.derivatives()[_Tidx], (hdis1 - hdis2) / (2.0 * dT), 1.0e-5);
934 }
935 
940 TEST_F(PorousFlowBrineCO2Test, solveEquilibriumMoleFractionHighTemp)
941 {
942  Real p = 20.0e6;
943  Real T = 473.15;
944  Real Xnacl = 0.0;
945 
946  Real yh2o, xco2;
947  Real co2_density = _co2_fp->rho_from_p_T(p, T);
948  _fs->solveEquilibriumMoleFractionHighTemp(p, T, Xnacl, co2_density, xco2, yh2o);
949  ABS_TEST(yh2o, 0.161429743509, 1.0e-10);
950  ABS_TEST(xco2, 0.0236966821858, 1.0e-10);
951 
952  // Mass fraction equivalent to molality of 2
953  p = 30.0e6;
954  T = 423.15;
955  Xnacl = 0.105;
956 
957  _fs->solveEquilibriumMoleFractionHighTemp(p, T, Xnacl, co2_density, xco2, yh2o);
958  ABS_TEST(yh2o, 0.0468197608955, 1.0e-10);
959  ABS_TEST(xco2, 0.0236644599437, 1.0e-10);
960 
961  // Mass fraction equivalent to molality of 4
962  T = 523.15;
963  Xnacl = 0.189;
964 
965  co2_density = _co2_fp->rho_from_p_T(p, T);
966  _fs->solveEquilibriumMoleFractionHighTemp(p, T, Xnacl, co2_density, xco2, yh2o);
967  ABS_TEST(yh2o, 0.253292728991, 1.0e-10);
968  ABS_TEST(xco2, 0.0168344513321, 1.0e-10);
969 }
970 
974 TEST_F(PorousFlowBrineCO2Test, equilibriumMoleFractions)
975 {
976  // Test pure water (Xnacl = 0)
977  // Low temperature regime
978  ADReal p = 20.0e6;
979  ADReal T = 323.15;
980  ADReal Xnacl = 0.0;
981 
982  ADReal x, y;
983  _fs->equilibriumMoleFractions(p, T, Xnacl, x, y);
984  ABS_TEST(y.value(), 0.00696393845155, 1.0e-8);
985  ABS_TEST(x.value(), 0.0236554537395, 1.0e-8);
986 
987  // Intermediate temperature regime
988  T = 373.15;
989 
990  _fs->equilibriumMoleFractions(p, T, Xnacl, x, y);
991  ABS_TEST(y.value(), 0.0194394631944, 1.0e-8);
992  ABS_TEST(x.value(), 0.020195776139, 1.0e-8);
993 
994  // High temperature regime
995  p = 30.0e6;
996  T = 523.15;
997 
998  _fs->equilibriumMoleFractions(p, T, Xnacl, x, y);
999  ABS_TEST(y.value(), 0.286117195565, 1.0e-8);
1000  ABS_TEST(x.value(), 0.0409622051253, 1.0e-8);
1001 
1002  // High temperature regime with low pressure
1003  p = 1.0e6;
1004  T = 523.15;
1005 
1006  _fs->equilibriumMoleFractions(p, T, Xnacl, x, y);
1007  ABS_TEST(y.value(), 1.0, 1.0e-8);
1008  ABS_TEST(x.value(), 0.0, 1.0e-8);
1009 
1010  // Test brine (Xnacl = 0.1)
1011  // Low temperature regime
1012  p = 20.0e6;
1013  T = 323.15;
1014  Xnacl = 0.1;
1015 
1016  _fs->equilibriumMoleFractions(p, T, Xnacl, x, y);
1017  ABS_TEST(y.value(), 0.00657335846826, 1.0e-8);
1018  ABS_TEST(x.value(), 0.0152863933134, 1.0e-8);
1019 
1020  // Intermediate temperature regime
1021  T = 373.15;
1022 
1023  _fs->equilibriumMoleFractions(p, T, Xnacl, x, y);
1024  ABS_TEST(y.value(), 0.01831360857, 1.0e-8);
1025  ABS_TEST(x.value(), 0.0132653916293, 1.0e-8);
1026 
1027  // High temperature regime
1028  p = 30.0e6;
1029  T = 523.15;
1030 
1031  _fs->equilibriumMoleFractions(p, T, Xnacl, x, y);
1032  ABS_TEST(y.value(), 0.270258924237, 1.0e-8);
1033  ABS_TEST(x.value(), 0.0246589135776, 1.0e-8);
1034 }
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 Thu Jul 17 2025 01:44:23 for https://mooseframework.inl.gov by   doxygen 1.8.14