https://mooseframework.inl.gov
PorousFlowWaterNCGTest.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 "PorousFlowWaterNCGTest.h"
11 #include "FluidPropertiesTestUtils.h"
12 
16 TEST_F(PorousFlowWaterNCGTest, name) { EXPECT_EQ("water-ncg", _fs->fluidStateName()); }
17 
18 /*
19  * Verify calculation of equilibrium mass fraction and derivatives
20  */
21 TEST_F(PorousFlowWaterNCGTest, equilibriumMassFraction)
22 {
23  ADReal p = 1.0e6;
24  Moose::derivInsert(p.derivatives(), _pidx, 1.0);
25 
26  ADReal T = 350.0;
27  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
28 
29  const Real dp = 1.0e-1;
30  const Real dT = 1.0e-6;
31 
32  ADReal Xncg, Yh2o, Xncg1, Yh2o1, Xncg2, Yh2o2;
33  _fs->equilibriumMassFractions(p, T, Xncg, Yh2o);
34  _fs->equilibriumMassFractions(p - dp, T, Xncg1, Yh2o1);
35  _fs->equilibriumMassFractions(p + dp, T, Xncg2, Yh2o2);
36 
37  Real dXncg_dp_fd = (Xncg2.value() - Xncg1.value()) / (2.0 * dp);
38  Real dYh2o_dp_fd = (Yh2o2.value() - Yh2o1.value()) / (2.0 * dp);
39 
40  REL_TEST(Xncg.derivatives()[_pidx], dXncg_dp_fd, 1.0e-8);
41  REL_TEST(Yh2o.derivatives()[_pidx], dYh2o_dp_fd, 1.0e-7);
42 
43  _fs->equilibriumMassFractions(p, T - dT, Xncg1, Yh2o1);
44  _fs->equilibriumMassFractions(p, T + dT, Xncg2, Yh2o2);
45 
46  Real dXncg_dT_fd = (Xncg2.value() - Xncg1.value()) / (2.0 * dT);
47  Real dYh2o_dT_fd = (Yh2o2.value() - Yh2o1.value()) / (2.0 * dT);
48 
49  REL_TEST(Xncg.derivatives()[_Tidx], dXncg_dT_fd, 1.0e-7);
50  REL_TEST(Yh2o.derivatives()[_Tidx], dYh2o_dT_fd, 1.0e-7);
51 }
52 
53 /*
54  * Verify calculation of actual mass fraction and derivatives depending on value of
55  * total mass fraction
56  */
57 TEST_F(PorousFlowWaterNCGTest, MassFraction)
58 {
59  ADReal p = 1.0e6;
60  Moose::derivInsert(p.derivatives(), _pidx, 1.0);
61 
62  ADReal T = 350.0;
63  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
64 
65  FluidStatePhaseEnum phase_state;
66  std::vector<FluidStateProperties> fsp(2, FluidStateProperties(2));
67 
68  // Liquid region
69  ADReal Z = 0.0001;
70  Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
71 
72  _fs->massFractions(p, T, Z, phase_state, fsp);
73  EXPECT_EQ(phase_state, FluidStatePhaseEnum::LIQUID);
74 
75  // Verfify mass fraction values
76  ADReal Xncg = fsp[0].mass_fraction[1];
77  ADReal Yncg = fsp[1].mass_fraction[1];
78  ADReal Xh2o = fsp[0].mass_fraction[0];
79  ADReal Yh2o = fsp[1].mass_fraction[0];
80  ABS_TEST(Xncg.value(), Z.value(), 1.0e-12);
81  ABS_TEST(Yncg.value(), 0.0, 1.0e-12);
82  ABS_TEST(Xh2o.value(), 1.0 - Z.value(), 1.0e-12);
83  ABS_TEST(Yh2o.value(), 0.0, 1.0e-12);
84 
85  // Verify derivatives
86  Real dXncg_dp = Xncg.derivatives()[_pidx];
87  Real dXncg_dT = Xncg.derivatives()[_Tidx];
88  Real dXncg_dZ = Xncg.derivatives()[_Zidx];
89  Real dYncg_dp = Yncg.derivatives()[_pidx];
90  Real dYncg_dT = Yncg.derivatives()[_Tidx];
91  Real dYncg_dZ = Yncg.derivatives()[_Zidx];
92  ABS_TEST(dXncg_dp, 0.0, 1.0e-12);
93  ABS_TEST(dXncg_dT, 0.0, 1.0e-12);
94  ABS_TEST(dXncg_dZ, 1.0, 1.0e-12);
95  ABS_TEST(dYncg_dp, 0.0, 1.0e-12);
96  ABS_TEST(dYncg_dT, 0.0, 1.0e-12);
97  ABS_TEST(dYncg_dZ, 0.0, 1.0e-12);
98 
99  // Gas region
100  Z = 0.995;
101  Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
102 
103  _fs->massFractions(p, T, Z, phase_state, fsp);
104  EXPECT_EQ(phase_state, FluidStatePhaseEnum::GAS);
105 
106  // Verfify mass fraction values
107  Xncg = fsp[0].mass_fraction[1];
108  Yncg = fsp[1].mass_fraction[1];
109  Xh2o = fsp[0].mass_fraction[0];
110  Yh2o = fsp[1].mass_fraction[0];
111  ABS_TEST(Xncg.value(), 0.0, 1.0e-12);
112  ABS_TEST(Yncg.value(), Z.value(), 1.0e-12);
113  ABS_TEST(Xh2o.value(), 0.0, 1.0e-12);
114  ABS_TEST(Yh2o.value(), 1.0 - Z.value(), 1.0e-12);
115 
116  // Verify derivatives
117  dXncg_dp = Xncg.derivatives()[_pidx];
118  dXncg_dT = Xncg.derivatives()[_Tidx];
119  dXncg_dZ = Xncg.derivatives()[_Zidx];
120  dYncg_dp = Yncg.derivatives()[_pidx];
121  dYncg_dT = Yncg.derivatives()[_Tidx];
122  dYncg_dZ = Yncg.derivatives()[_Zidx];
123  ABS_TEST(dXncg_dp, 0.0, 1.0e-12);
124  ABS_TEST(dXncg_dT, 0.0, 1.0e-12);
125  ABS_TEST(dXncg_dZ, 0.0, 1.0e-12);
126  ABS_TEST(dYncg_dp, 0.0, 1.0e-12);
127  ABS_TEST(dYncg_dT, 0.0, 1.0e-12);
128  ABS_TEST(dYncg_dZ, 1.0, 1.0e-12);
129 
130  // Two phase region. In this region, the mass fractions and derivatives can
131  // be verified using the equilibrium mass fraction derivatives that have
132  // been verified above
133  Z = 0.45;
134  Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
135 
136  _fs->massFractions(p, T, Z, phase_state, fsp);
137  EXPECT_EQ(phase_state, FluidStatePhaseEnum::TWOPHASE);
138 
139  // Equilibrium mass fractions and derivatives
140  ADReal Xncg_eq, Yh2o_eq;
141  _fs->equilibriumMassFractions(p, T, Xncg_eq, Yh2o_eq);
142 
143  // Verfify mass fraction values
144  Xncg = fsp[0].mass_fraction[1];
145  Yncg = fsp[1].mass_fraction[1];
146  Xh2o = fsp[0].mass_fraction[0];
147  Yh2o = fsp[1].mass_fraction[0];
148  ABS_TEST(Xncg, Xncg_eq, 1.0e-12);
149  ABS_TEST(Yncg, 1.0 - Yh2o_eq, 1.0e-12);
150  ABS_TEST(Xh2o, 1.0 - Xncg_eq, 1.0e-12);
151  ABS_TEST(Yh2o, Yh2o_eq, 1.0e-12);
152 
153  // Verify derivatives wrt p and T
154  dXncg_dp = Xncg.derivatives()[_pidx];
155  dXncg_dT = Xncg.derivatives()[_Tidx];
156  dXncg_dZ = Xncg.derivatives()[_Zidx];
157  dYncg_dp = Yncg.derivatives()[_pidx];
158  dYncg_dT = Yncg.derivatives()[_Tidx];
159  dYncg_dZ = Yncg.derivatives()[_Zidx];
160  ABS_TEST(dXncg_dp, Xncg_eq.derivatives()[_pidx], 1.0e-8);
161  ABS_TEST(dXncg_dT, Xncg_eq.derivatives()[_Tidx], 1.0e-8);
162  ABS_TEST(dXncg_dZ, 0.0, 1.0e-8);
163  ABS_TEST(dYncg_dp, -Yh2o_eq.derivatives()[_pidx], 1.0e-8);
164  ABS_TEST(dYncg_dT, -Yh2o_eq.derivatives()[_Tidx], 1.0e-8);
165  ABS_TEST(dYncg_dZ, 0.0, 1.0e-8);
166 
167  // Use finite differences to verify derivative wrt Z is unaffected by Z
168  const Real dZ = 1.0e-8;
169  _fs->massFractions(p, T, Z + dZ, phase_state, fsp);
170  ADReal Xncg1 = fsp[0].mass_fraction[1];
171  ADReal Yncg1 = fsp[1].mass_fraction[1];
172  _fs->massFractions(p, T, Z - dZ, phase_state, fsp);
173  ADReal Xncg2 = fsp[0].mass_fraction[1];
174  ADReal Yncg2 = fsp[1].mass_fraction[1];
175 
176  ABS_TEST(dXncg_dZ, (Xncg1 - Xncg2).value() / (2.0 * dZ), 1.0e-8);
177  ABS_TEST(dYncg_dZ, (Yncg1 - Yncg2).value() / (2.0 * dZ), 1.0e-8);
178 }
179 
180 /*
181  * Verify calculation of gas density and derivatives
182  */
183 TEST_F(PorousFlowWaterNCGTest, gasDensity)
184 {
185  ADReal p = 1.0e6;
186  Moose::derivInsert(p.derivatives(), _pidx, 1.0);
187 
188  ADReal T = 350.0;
189  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
190 
191  FluidStatePhaseEnum phase_state;
192  std::vector<FluidStateProperties> fsp(2, FluidStateProperties(2));
193 
194  // Gas region
195  ADReal Z = 0.995;
196  Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
197 
198  _fs->massFractions(p, T, Z, phase_state, fsp);
199  EXPECT_EQ(phase_state, FluidStatePhaseEnum::GAS);
200 
201  const ADReal gas_density = _fs->gasDensity(p, T, fsp);
202 
203  const ADReal density =
204  _ncg_fp->rho_from_p_T(Z * p, T) + _water_fp->rho_from_p_T(_water_fp->vaporPressure(T), T);
205 
206  ABS_TEST(gas_density, density, 1.0e-12);
207 }
208 
209 /*
210  * Verify calculation of gas density, viscosity, enthalpy and derivatives
211  */
212 TEST_F(PorousFlowWaterNCGTest, gasProperties)
213 {
214  ADReal p = 1.0e6;
215  Moose::derivInsert(p.derivatives(), _pidx, 1.0);
216 
217  ADReal T = 350.0;
218  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
219 
220  FluidStatePhaseEnum phase_state;
221  std::vector<FluidStateProperties> fsp(2, FluidStateProperties(2));
222 
223  // Gas region
224  ADReal Z = 0.995;
225  Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
226 
227  _fs->massFractions(p, T, Z, phase_state, fsp);
228  EXPECT_EQ(phase_state, FluidStatePhaseEnum::GAS);
229 
230  // Verify fluid density, viscosity and enthalpy
231  _fs->gasProperties(p, T, fsp);
232  ADReal gas_density = fsp[1].density;
233  ADReal gas_viscosity = fsp[1].viscosity;
234  ADReal gas_enthalpy = fsp[1].enthalpy;
235 
236  ADReal density =
237  _ncg_fp->rho_from_p_T(Z * p, T) + _water_fp->rho_from_p_T(_water_fp->vaporPressure(T), T);
238  ADReal viscosity = Z * _ncg_fp->mu_from_p_T(Z * p, T) +
239  (1.0 - Z) * _water_fp->mu_from_p_T(_water_fp->vaporPressure(T), T);
240  ADReal enthalpy = Z * _ncg_fp->h_from_p_T(Z * p, T) +
241  (1.0 - Z) * _water_fp->h_from_p_T(_water_fp->vaporPressure(T), T);
242 
243  ABS_TEST(gas_density, density, 1.0e-10);
244  ABS_TEST(gas_viscosity, viscosity, 1.0e-10);
245  ABS_TEST(gas_enthalpy, enthalpy, 1.0e-10);
246 
247  // Verify derivatives
248  Real ddensity_dp = gas_density.derivatives()[_pidx];
249  Real ddensity_dT = gas_density.derivatives()[_Tidx];
250  Real ddensity_dZ = gas_density.derivatives()[_Zidx];
251  Real dviscosity_dp = gas_viscosity.derivatives()[_pidx];
252  Real dviscosity_dT = gas_viscosity.derivatives()[_Tidx];
253  Real dviscosity_dZ = gas_viscosity.derivatives()[_Zidx];
254  Real denthalpy_dp = gas_enthalpy.derivatives()[_pidx];
255  Real denthalpy_dT = gas_enthalpy.derivatives()[_Tidx];
256  Real denthalpy_dZ = gas_enthalpy.derivatives()[_Zidx];
257 
258  const Real dp = 1.0e-1;
259  _fs->gasProperties(p + dp, T, fsp);
260  ADReal rho1 = fsp[1].density;
261  ADReal mu1 = fsp[1].viscosity;
262  ADReal h1 = fsp[1].enthalpy;
263 
264  _fs->gasProperties(p - dp, T, fsp);
265  ADReal rho2 = fsp[1].density;
266  ADReal mu2 = fsp[1].viscosity;
267  ADReal h2 = fsp[1].enthalpy;
268 
269  REL_TEST(ddensity_dp, (rho1 - rho2).value() / (2.0 * dp), 1.0e-7);
270  REL_TEST(dviscosity_dp, (mu1 - mu2).value() / (2.0 * dp), 1.0e-7);
271  REL_TEST(denthalpy_dp, (h1 - h2).value() / (2.0 * dp), 5.0e-7);
272 
273  const Real dT = 1.0e-3;
274  _fs->gasProperties(p, T + dT, fsp);
275  rho1 = fsp[1].density;
276  mu1 = fsp[1].viscosity;
277  h1 = fsp[1].enthalpy;
278 
279  _fs->gasProperties(p, T - dT, fsp);
280  rho2 = fsp[1].density;
281  mu2 = fsp[1].viscosity;
282  h2 = fsp[1].enthalpy;
283 
284  REL_TEST(ddensity_dT, (rho1 - rho2).value() / (2.0 * dT), 1.0e-7);
285  REL_TEST(dviscosity_dT, (mu1 - mu2).value() / (2.0 * dT), 1.0e-8);
286  REL_TEST(denthalpy_dT, (h1 - h2).value() / (2.0 * dT), 1.0e-8);
287 
288  // Note: mass fraction changes with Z
289  const Real dZ = 1.0e-8;
290  _fs->massFractions(p, T, Z + dZ, phase_state, fsp);
291  _fs->gasProperties(p, T, fsp);
292  rho1 = fsp[1].density;
293  mu1 = fsp[1].viscosity;
294  h1 = fsp[1].enthalpy;
295 
296  _fs->massFractions(p, T, Z - dZ, phase_state, fsp);
297  _fs->gasProperties(p, T, fsp);
298  rho2 = fsp[1].density;
299  mu2 = fsp[1].viscosity;
300  h2 = fsp[1].enthalpy;
301 
302  REL_TEST(ddensity_dZ, (rho1 - rho2).value() / (2.0 * dZ), 5.0e-8);
303  REL_TEST(dviscosity_dZ, (mu1 - mu2).value() / (2.0 * dZ), 5.0e-8);
304  REL_TEST(denthalpy_dZ, (h1 - h2).value() / (2.0 * dZ), 1.0e-8);
305 
306  // Check derivatives in the two phase region as well. Note that the mass fractions
307  // vary with pressure and temperature in this region
308  Z = 0.45;
309  Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
310 
311  _fs->massFractions(p, T, Z, phase_state, fsp);
312  EXPECT_EQ(phase_state, FluidStatePhaseEnum::TWOPHASE);
313 
314  _fs->gasProperties(p, T, fsp);
315  gas_density = fsp[1].density;
316  gas_viscosity = fsp[1].viscosity;
317  gas_enthalpy = fsp[1].enthalpy;
318  ddensity_dp = gas_density.derivatives()[_pidx];
319  ddensity_dT = gas_density.derivatives()[_Tidx];
320  ddensity_dZ = gas_density.derivatives()[_Zidx];
321  dviscosity_dp = gas_viscosity.derivatives()[_pidx];
322  dviscosity_dT = gas_viscosity.derivatives()[_Tidx];
323  dviscosity_dZ = gas_viscosity.derivatives()[_Zidx];
324  denthalpy_dp = gas_enthalpy.derivatives()[_pidx];
325  denthalpy_dT = gas_enthalpy.derivatives()[_Tidx];
326  denthalpy_dZ = gas_enthalpy.derivatives()[_Zidx];
327 
328  _fs->massFractions(p + dp, T, Z, phase_state, fsp);
329  _fs->gasProperties(p + dp, T, fsp);
330  rho1 = fsp[1].density;
331  mu1 = fsp[1].viscosity;
332  h1 = fsp[1].enthalpy;
333 
334  _fs->massFractions(p - dp, T, Z, phase_state, fsp);
335  _fs->gasProperties(p - dp, T, fsp);
336  rho2 = fsp[1].density;
337  mu2 = fsp[1].viscosity;
338  h2 = fsp[1].enthalpy;
339 
340  REL_TEST(ddensity_dp, (rho1 - rho2).value() / (2.0 * dp), 1.0e-7);
341  REL_TEST(dviscosity_dp, (mu1 - mu2).value() / (2.0 * dp), 1.0e-7);
342  REL_TEST(denthalpy_dp, (h1 - h2).value() / (2.0 * dp), 1.0e-7);
343 
344  _fs->massFractions(p, T + dT, Z, phase_state, fsp);
345  _fs->gasProperties(p, T + dT, fsp);
346  rho1 = fsp[1].density;
347  mu1 = fsp[1].viscosity;
348  h1 = fsp[1].enthalpy;
349 
350  _fs->massFractions(p, T - dT, Z, phase_state, fsp);
351  _fs->gasProperties(p, T - dT, fsp);
352  rho2 = fsp[1].density;
353  mu2 = fsp[1].viscosity;
354  h2 = fsp[1].enthalpy;
355 
356  REL_TEST(ddensity_dT, (rho1 - rho2).value() / (2.0 * dT), 1.0e-7);
357  REL_TEST(dviscosity_dT, (mu1 - mu2).value() / (2.0 * dT), 1.0e-8);
358  REL_TEST(denthalpy_dT, (h1 - h2).value() / (2.0 * dT), 1.0e-8);
359 
360  _fs->massFractions(p, T, Z + dZ, phase_state, fsp);
361  _fs->gasProperties(p, T, fsp);
362  rho1 = fsp[1].density;
363  mu1 = fsp[1].viscosity;
364  h1 = fsp[1].enthalpy;
365 
366  _fs->massFractions(p, T, Z - dZ, phase_state, fsp);
367  _fs->gasProperties(p, T, fsp);
368  rho2 = fsp[1].density;
369  mu2 = fsp[1].viscosity;
370  h2 = fsp[1].enthalpy;
371 
372  ABS_TEST(ddensity_dZ, (rho1 - rho2).value() / (2.0 * dZ), 1.0e-8);
373  ABS_TEST(dviscosity_dT, (mu1 - mu2).value() / (2.0 * dZ), 1.0e-7);
374  ABS_TEST(denthalpy_dZ, (h1 - h2).value() / (2.0 * dZ), 1.0e-8);
375 }
376 
377 /*
378  * Verify calculation of liquid density and derivatives
379  */
380 TEST_F(PorousFlowWaterNCGTest, liquidDensity)
381 {
382  ADReal p = 1.0e6;
383  Moose::derivInsert(p.derivatives(), _pidx, 1.0);
384 
385  ADReal T = 350.0;
386  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
387 
388  const ADReal liquid_density = _fs->liquidDensity(p, T);
389 
390  const ADReal density = _water_fp->rho_from_p_T(p, T);
391  ABS_TEST(liquid_density, density, 1.0e-12);
392 }
393 
394 /*
395  * Verify calculation of liquid density, viscosity, enthalpy and derivatives. Note that as
396  * density and viscosity don't depend on mass fraction, only the liquid region needs to be
397  * tested (the calculations are identical in the two phase region). The enthalpy does depend
398  * on mass fraction, so should be tested in the two phase region as well as the liquid region
399  */
400 TEST_F(PorousFlowWaterNCGTest, liquidProperties)
401 {
402  ADReal p = 1.0e6;
403  Moose::derivInsert(p.derivatives(), _pidx, 1.0);
404 
405  ADReal T = 350.0;
406  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
407 
408  FluidStatePhaseEnum phase_state;
409  std::vector<FluidStateProperties> fsp(2, FluidStateProperties(2));
410 
411  // Verify fluid density and viscosity
412  _fs->liquidProperties(p, T, fsp);
413  ADReal liquid_density = fsp[0].density;
414  ADReal liquid_viscosity = fsp[0].viscosity;
415  ADReal liquid_enthalpy = fsp[0].enthalpy;
416 
417  ADReal density = _water_fp->rho_from_p_T(p, T);
418  ADReal viscosity = _water_fp->mu_from_p_T(p, T);
419  ADReal enthalpy = _water_fp->h_from_p_T(p, T);
420 
421  ABS_TEST(liquid_density, density, 1.0e-12);
422  ABS_TEST(liquid_viscosity, viscosity, 1.0e-12);
423  ABS_TEST(liquid_enthalpy, enthalpy, 1.0e-12);
424 
425  // Verify derivatives
426  Real ddensity_dp = liquid_density.derivatives()[_pidx];
427  Real ddensity_dT = liquid_density.derivatives()[_Tidx];
428  Real ddensity_dZ = liquid_density.derivatives()[_Zidx];
429  Real dviscosity_dp = liquid_viscosity.derivatives()[_pidx];
430  Real dviscosity_dT = liquid_viscosity.derivatives()[_Tidx];
431  Real dviscosity_dZ = liquid_viscosity.derivatives()[_Zidx];
432  Real denthalpy_dp = liquid_enthalpy.derivatives()[_pidx];
433  Real denthalpy_dT = liquid_enthalpy.derivatives()[_Tidx];
434  Real denthalpy_dZ = liquid_enthalpy.derivatives()[_Zidx];
435 
436  const Real dp = 10.0;
437  _fs->liquidProperties(p + dp, T, fsp);
438  ADReal rho1 = fsp[0].density;
439  ADReal mu1 = fsp[0].viscosity;
440  ADReal h1 = fsp[0].enthalpy;
441 
442  _fs->liquidProperties(p - dp, T, fsp);
443  ADReal rho2 = fsp[0].density;
444  ADReal mu2 = fsp[0].viscosity;
445  ADReal h2 = fsp[0].enthalpy;
446 
447  REL_TEST(ddensity_dp, (rho1 - rho2) / (2.0 * dp), 1.0e-6);
448  REL_TEST(dviscosity_dp, (mu1 - mu2) / (2.0 * dp), 1.0e-6);
449  REL_TEST(denthalpy_dp, (h1 - h2) / (2.0 * dp), 1.0e-6);
450 
451  const Real dT = 1.0e-4;
452  _fs->liquidProperties(p, T + dT, fsp);
453  rho1 = fsp[0].density;
454  mu1 = fsp[0].viscosity;
455  h1 = fsp[0].enthalpy;
456 
457  _fs->liquidProperties(p, T - dT, fsp);
458  rho2 = fsp[0].density;
459  mu2 = fsp[0].viscosity;
460  h2 = fsp[0].enthalpy;
461 
462  REL_TEST(ddensity_dT, (rho1 - rho2) / (2.0 * dT), 1.0e-8);
463  REL_TEST(dviscosity_dT, (mu1 - mu2) / (2.0 * dT), 1.0e-8);
464  REL_TEST(denthalpy_dT, (h1 - h2) / (2.0 * dT), 1.0e-8);
465 
466  ADReal Z = 0.0001;
467  Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
468 
469  const Real dZ = 1.0e-8;
470 
471  _fs->massFractions(p, T, Z, phase_state, fsp);
472  _fs->liquidProperties(p, T, fsp);
473  denthalpy_dZ = fsp[0].enthalpy.derivatives()[_Zidx];
474 
475  _fs->massFractions(p, T, Z + dZ, phase_state, fsp);
476  _fs->liquidProperties(p, T, fsp);
477  h1 = fsp[0].enthalpy;
478 
479  _fs->massFractions(p, T, Z - dZ, phase_state, fsp);
480  _fs->liquidProperties(p, T, fsp);
481  h2 = fsp[0].enthalpy;
482 
483  REL_TEST(denthalpy_dZ, (h1 - h2) / (2.0 * dZ), 1.0e-8);
484 
485  // Density and viscosity don't depend on Z, so derivatives should be 0
486  ABS_TEST(ddensity_dZ, 0.0, 1.0e-12);
487  ABS_TEST(dviscosity_dZ, 0.0, 1.0e-12);
488 
489  // Check enthalpy calculations in the two phase region as well. Note that the mass fractions
490  // vary with pressure and temperature in this region
491  Z = 0.45;
492  Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
493 
494  _fs->massFractions(p, T, Z, phase_state, fsp);
495  _fs->liquidProperties(p, T, fsp);
496  denthalpy_dp = fsp[0].enthalpy.derivatives()[_pidx];
497  denthalpy_dT = fsp[0].enthalpy.derivatives()[_Tidx];
498  denthalpy_dZ = fsp[0].enthalpy.derivatives()[_Zidx];
499 
500  _fs->massFractions(p + dp, T, Z, phase_state, fsp);
501  _fs->liquidProperties(p + dp, T, fsp);
502  h1 = fsp[0].enthalpy;
503 
504  _fs->massFractions(p - dp, T, Z, phase_state, fsp);
505  _fs->liquidProperties(p - dp, T, fsp);
506  h2 = fsp[0].enthalpy;
507 
508  REL_TEST(denthalpy_dp, (h1 - h2) / (2.0 * dp), 1.0e-8);
509 
510  _fs->massFractions(p, T + dT, Z, phase_state, fsp);
511  _fs->liquidProperties(p, T + dT, fsp);
512  h1 = fsp[0].enthalpy;
513 
514  _fs->massFractions(p, T - dT, Z, phase_state, fsp);
515  _fs->liquidProperties(p, T - dT, fsp);
516  h2 = fsp[0].enthalpy;
517 
518  REL_TEST(denthalpy_dT, (h1 - h2) / (2.0 * dT), 1.0e-8);
519 
520  _fs->massFractions(p, T, Z + dZ, phase_state, fsp);
521  _fs->liquidProperties(p, T, fsp);
522  h1 = fsp[0].enthalpy;
523 
524  _fs->massFractions(p, T, Z - dZ, phase_state, fsp);
525  _fs->liquidProperties(p, T, fsp);
526  h2 = fsp[0].enthalpy;
527 
528  ABS_TEST(denthalpy_dZ, (h1 - h2) / (2.0 * dZ), 1.0e-8);
529 }
530 
531 /*
532  * Verify calculation of gas saturation and derivatives in the two-phase region
533  */
534 TEST_F(PorousFlowWaterNCGTest, saturation)
535 {
536  ADReal p = 1.0e6;
537  Moose::derivInsert(p.derivatives(), _pidx, 1.0);
538 
539  ADReal T = 350.0;
540  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
541 
542  FluidStatePhaseEnum phase_state;
543  std::vector<FluidStateProperties> fsp(2, FluidStateProperties(2));
544 
545  // In the two-phase region, the mass fractions are the equilibrium values, so
546  // a temporary value of Z can be used (as long as it corresponds to the two-phase
547  // region)
548  ADReal Z = 0.45;
549  Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
550 
551  _fs->massFractions(p, T, Z, phase_state, fsp);
552  EXPECT_EQ(phase_state, FluidStatePhaseEnum::TWOPHASE);
553 
554  // Calculate Z that gives a saturation of 0.25
555  ADReal target_gas_saturation = 0.25;
556  fsp[1].saturation = target_gas_saturation;
557 
558  // Calculate gas density and liquid density
559  _fs->gasProperties(p, T, fsp);
560  _fs->liquidProperties(p, T, fsp);
561 
562  // The mass fraction that corresponds to a gas_saturation = 0.25
563  ADReal Zc =
564  (target_gas_saturation * fsp[1].density * fsp[1].mass_fraction[1] +
565  (1.0 - target_gas_saturation) * fsp[0].density * fsp[0].mass_fraction[1]) /
566  (target_gas_saturation * fsp[1].density + (1.0 - target_gas_saturation) * fsp[0].density);
567 
568  // Calculate the gas saturation and derivatives
569  ADReal gas_saturation = _fs->saturation(p, T, Zc, fsp);
570 
571  REL_TEST(gas_saturation, target_gas_saturation, 1.0e-6);
572 
573  // Test the derivatives of gas saturation
574  gas_saturation = _fs->saturation(p, T, Z, fsp);
575  const Real dp = 1.0e-1;
576 
577  Real dgas_saturation_dp = gas_saturation.derivatives()[_pidx];
578  Real dgas_saturation_dT = gas_saturation.derivatives()[_Tidx];
579  Real dgas_saturation_dZ = gas_saturation.derivatives()[_Zidx];
580 
581  _fs->massFractions(p + dp, T, Z, phase_state, fsp);
582  ADReal gsat1 = _fs->saturation(p + dp, T, Z, fsp);
583 
584  _fs->massFractions(p - dp, T, Z, phase_state, fsp);
585  ADReal gsat2 = _fs->saturation(p - dp, T, Z, fsp);
586 
587  REL_TEST(dgas_saturation_dp, (gsat1 - gsat2).value() / (2.0 * dp), 1.0e-7);
588 
589  const Real dT = 1.0e-4;
590 
591  _fs->massFractions(p, T + dT, Z, phase_state, fsp);
592  gsat1 = _fs->saturation(p, T + dT, Z, fsp);
593 
594  _fs->massFractions(p, T - dT, Z, phase_state, fsp);
595  gsat2 = _fs->saturation(p, T - dT, Z, fsp);
596 
597  REL_TEST(dgas_saturation_dT, (gsat1 - gsat2).value() / (2.0 * dT), 1.0e-7);
598 
599  const Real dZ = 1.0e-8;
600 
601  _fs->massFractions(p, T, Z + dZ, phase_state, fsp);
602  gsat1 = _fs->saturation(p, T, Z + dZ, fsp);
603 
604  _fs->massFractions(p, T, Z - dZ, phase_state, fsp);
605  gsat2 = _fs->saturation(p, T, Z - dZ, fsp);
606 
607  REL_TEST(dgas_saturation_dZ, (gsat1 - gsat2).value() / (2.0 * dZ), 1.0e-7);
608 }
609 
610 /*
611  * Verify calculation of gas saturation and derivatives in the two-phase region
612  */
613 TEST_F(PorousFlowWaterNCGTest, twoPhase)
614 {
615  ADReal p = 1.0e6;
616  Moose::derivInsert(p.derivatives(), _pidx, 1.0);
617 
618  ADReal T = 350.0;
619  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
620 
621  FluidStatePhaseEnum phase_state;
622  std::vector<FluidStateProperties> fsp(2, FluidStateProperties(2));
623 
624  ADReal Z = 0.45;
625  Moose::derivInsert(Z.derivatives(), _Zidx, 1.0);
626 
627  // Dummy qp value that isn't needed to test this
628  const unsigned int qp = 0;
629 
630  _fs->massFractions(p, T, Z, phase_state, fsp);
631  EXPECT_EQ(phase_state, FluidStatePhaseEnum::TWOPHASE);
632 
633  _fs->twoPhaseProperties(p, T, Z, qp, fsp);
634  ADReal liquid_density = fsp[0].density;
635  ADReal liquid_viscosity = fsp[0].viscosity;
636  ADReal liquid_enthalpy = fsp[0].enthalpy;
637  ADReal gas_saturation = fsp[1].saturation;
638  ADReal gas_density = fsp[1].density;
639  ADReal gas_viscosity = fsp[1].viscosity;
640  ADReal gas_enthalpy = fsp[1].enthalpy;
641 
642  // As all the values are provided by the gasProperties(), liquidProperties() and
643  // saturation(), we can simply compare the ADReal results are the same (really
644  // only need to make sure that liquid properties are calculated correctly given the saturation)
645  std::vector<FluidStateProperties> fsp2(2, FluidStateProperties(2));
646 
647  _fs->massFractions(p, T, Z, phase_state, fsp2);
648  _fs->gasProperties(p, T, fsp2);
649  fsp2[1].saturation = _fs->saturation(p, T, Z, fsp2);
650 
651  const ADReal pliq = p - _pc->capillaryPressure(1.0 - fsp2[1].saturation, qp);
652  _fs->liquidProperties(pliq, T, fsp2);
653 
654  ABS_TEST(gas_density, fsp2[1].density, 1.0e-12);
655  ABS_TEST(gas_viscosity, fsp2[1].viscosity, 1.0e-12);
656  ABS_TEST(gas_enthalpy, fsp2[1].enthalpy, 1.0e-12);
657 
658  ABS_TEST(liquid_density, fsp2[0].density, 1.0e-12);
659  ABS_TEST(liquid_viscosity, fsp2[0].viscosity, 1.0e-12);
660  ABS_TEST(liquid_enthalpy, fsp2[0].enthalpy, 1.0e-12);
661 }
662 
663 /*
664  * Verify calculation of total mass fraction given a gas saturation
665  */
666 TEST_F(PorousFlowWaterNCGTest, totalMassFraction)
667 {
668  const Real p = 1.0e6;
669  const Real T = 350.0;
670  const Real s = 0.2;
671  const Real Xnacl = 0.1;
672  const unsigned qp = 0;
673 
674  Real Z = _fs->totalMassFraction(p, T, Xnacl, s, qp);
675 
676  // Test that the saturation calculated in this fluid state using Z is equal to s
677  FluidStatePhaseEnum phase_state;
678  std::vector<FluidStateProperties> fsp(2, FluidStateProperties(2));
679 
680  ADReal pressure = p;
681  ADReal temperature = T;
682  ADReal z = Z;
683  _fs->massFractions(pressure, temperature, z, phase_state, fsp);
684  EXPECT_EQ(phase_state, FluidStatePhaseEnum::TWOPHASE);
685 
686  ADReal gas_saturation = _fs->saturation(pressure, temperature, z, fsp);
687  REL_TEST(gas_saturation, s, 1.0e-5);
688 }
689 
690 /*
691  * Verify calculation of enthalpy of dissolution. Note: the values calculated compare
692  * well by eye to the values presented in Figure 4 of Battistelli et al, "A fluid property
693  * module for the TOUGH2 simulator for saline brines with non-condensible gas"
694  */
695 TEST_F(PorousFlowWaterNCGTest, enthalpyOfDissolution)
696 {
697  // T = 50C
698  ADReal T = 323.15;
699  Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
700 
701  // Enthalpy of dissolution of NCG in water
702  ADReal hdis = _fs->enthalpyOfDissolution(T);
703  REL_TEST(hdis.value(), -3.45731e5, 1.0e-5);
704 
705  // T = 350C
706  T = 623.15;
707  Moose::derivInsert(T.derivatives(), 2, 1.0);
708 
709  // Enthalpy of dissolution of NCG in water
710  hdis = _fs->enthalpyOfDissolution(T);
711  REL_TEST(hdis.value(), 1.23423e+06, 1.0e-5);
712 
713  // Test the derivative wrt temperature
714  const Real dT = 1.0e-4;
715  ADReal hdis2 = _fs->enthalpyOfDissolution(T + dT);
716  ADReal hdis3 = _fs->enthalpyOfDissolution(T - dT);
717 
718  REL_TEST(hdis.derivatives()[_Tidx], (hdis2 - hdis3) / (2.0 * dT), 1.0e-6);
719 }
NS::density
static const std::string density
Definition: NS.h:33
NS::temperature
static const std::string temperature
Definition: NS.h:59
ADReal
DualNumber< Real, DNDerivativeType, true > ADReal
TEST_F
TEST_F(PorousFlowWaterNCGTest, name)
Verify that the correct name is supplied.
Definition: PorousFlowWaterNCGTest.C:16
value
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
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 Wed Jul 16 2025 14:09:40 for https://mooseframework.inl.gov by   doxygen 1.8.14