https://mooseframework.inl.gov
Loading...
Searching...
No Matches
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
16TEST_F(PorousFlowWaterNCGTest, name) { EXPECT_EQ("water-ncg", _fs->fluidStateName()); }
17
18/*
19 * Verify calculation of equilibrium mass fraction and derivatives
20 */
21TEST_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 */
57TEST_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 */
183TEST_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 */
212TEST_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 */
380TEST_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 */
400TEST_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 */
534TEST_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 */
613TEST_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 */
666TEST_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 */
695TEST_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}
ADReal
DualNumber< Real, DNDerivativeType, true > ADReal
p
const Real p
Definition GasLiquidMassTransferTest.C:18
T
const double T
Definition GasLiquidMassTransferTest.C:19
FluidStatePhaseEnum
FluidStatePhaseEnum
Phase state enum.
Definition PorousFlowFluidStateBase.h:17
TEST_F
TEST_F(PorousFlowWaterNCGTest, name)
Verify that the correct name is supplied.
Definition PorousFlowWaterNCGTest.C:16
name
const std::string name
Definition Setup.h:21
Moose::derivInsert
void derivInsert(SemiDynamicSparseNumberArray< Real, libMesh::dof_id_type, NWrapper< N > > &derivs, libMesh::dof_id_type index, Real value)
FluidStateProperties
AD data structure to pass calculated thermophysical properties.
Definition PorousFlowFluidStateBase.h:25
Generated on Fri Aug 21 2026 13:41:31 for https://mooseframework.inl.gov by doxygen 1.9.8