Line data Source code
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 "PorousFlowWaterNCG.h"
11 : #include "SinglePhaseFluidProperties.h"
12 : #include "Water97FluidProperties.h"
13 : #include "Conversion.h"
14 :
15 : registerMooseObject("PorousFlowApp", PorousFlowWaterNCG);
16 :
17 : InputParameters
18 468 : PorousFlowWaterNCG::validParams()
19 : {
20 468 : InputParameters params = PorousFlowFluidStateMultiComponentBase::validParams();
21 936 : params.addRequiredParam<UserObjectName>("water_fp", "The name of the user object for water");
22 936 : params.addRequiredParam<UserObjectName>(
23 : "gas_fp", "The name of the user object for the non-condensable gas");
24 468 : params.addClassDescription("Fluid state class for water and non-condensable gas");
25 468 : return params;
26 0 : }
27 :
28 234 : PorousFlowWaterNCG::PorousFlowWaterNCG(const InputParameters & parameters)
29 : : PorousFlowFluidStateMultiComponentBase(parameters),
30 234 : _water_fp(getUserObject<SinglePhaseFluidProperties>("water_fp")),
31 234 : _water97_fp(getUserObject<Water97FluidProperties>("water_fp")),
32 234 : _ncg_fp(getUserObject<SinglePhaseFluidProperties>("gas_fp")),
33 234 : _Mh2o(_water_fp.molarMass()),
34 234 : _Mncg(_ncg_fp.molarMass()),
35 234 : _water_triple_temperature(_water_fp.triplePointTemperature()),
36 234 : _water_critical_temperature(_water_fp.criticalTemperature()),
37 468 : _ncg_henry(_ncg_fp.henryCoefficients())
38 : {
39 : // Check that the correct FluidProperties UserObjects have been provided
40 468 : if (_water_fp.fluidName() != "water")
41 0 : paramError("water_fp", "A valid water FluidProperties UserObject must be provided in water_fp");
42 :
43 : // Set the number of phases and components, and their indexes
44 234 : _num_phases = 2;
45 234 : _num_components = 2;
46 234 : _gas_phase_number = 1 - _aqueous_phase_number;
47 234 : _gas_fluid_component = 1 - _aqueous_fluid_component;
48 :
49 : // Check that _aqueous_phase_number is <= total number of phases
50 234 : if (_aqueous_phase_number >= _num_phases)
51 0 : paramError("liquid_phase_number",
52 : "This value is larger than the possible number of phases ",
53 : _num_phases);
54 :
55 : // Check that _aqueous_fluid_component is <= total number of fluid components
56 234 : if (_aqueous_fluid_component >= _num_components)
57 0 : paramError("liquid_fluid_component",
58 : "This value is larger than the possible number of fluid components",
59 : _num_components);
60 :
61 234 : _empty_fsp = FluidStateProperties(_num_components);
62 234 : }
63 :
64 : std::string
65 1 : PorousFlowWaterNCG::fluidStateName() const
66 : {
67 1 : return "water-ncg";
68 : }
69 :
70 : void
71 912456 : PorousFlowWaterNCG::thermophysicalProperties(Real pressure,
72 : Real temperature,
73 : Real Xnacl,
74 : Real Z,
75 : unsigned int qp,
76 : std::vector<FluidStateProperties> & fsp) const
77 : {
78 : // Make AD versions of primary variables then call AD thermophysicalProperties()
79 912456 : ADReal p = pressure;
80 912456 : Moose::derivInsert(p.derivatives(), _pidx, 1.0);
81 912456 : ADReal T = temperature;
82 912456 : Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
83 912456 : ADReal Zncg = Z;
84 912456 : Moose::derivInsert(Zncg.derivatives(), _Zidx, 1.0);
85 : // X is not used, but needed for consistency with PorousFlowFluidState interface
86 912456 : ADReal X = Xnacl;
87 912456 : Moose::derivInsert(X.derivatives(), _Xidx, 1.0);
88 :
89 912456 : thermophysicalProperties(p, T, X, Zncg, qp, fsp);
90 912445 : }
91 :
92 : void
93 912456 : PorousFlowWaterNCG::thermophysicalProperties(const ADReal & pressure,
94 : const ADReal & temperature,
95 : const ADReal & /* Xnacl */,
96 : const ADReal & Z,
97 : unsigned int qp,
98 : std::vector<FluidStateProperties> & fsp) const
99 : {
100 912456 : FluidStateProperties & liquid = fsp[_aqueous_phase_number];
101 912456 : FluidStateProperties & gas = fsp[_gas_phase_number];
102 :
103 : // Check whether the input temperature is within the region of validity
104 912456 : checkVariables(temperature.value());
105 :
106 : // Clear all of the FluidStateProperties data
107 912445 : clearFluidStateProperties(fsp);
108 :
109 : FluidStatePhaseEnum phase_state;
110 912445 : massFractions(pressure, temperature, Z, phase_state, fsp);
111 :
112 912445 : switch (phase_state)
113 : {
114 8248 : case FluidStatePhaseEnum::GAS:
115 : {
116 : // Set the gas saturations
117 8248 : gas.saturation = 1.0;
118 :
119 : // Calculate gas properties
120 8248 : gasProperties(pressure, temperature, fsp);
121 :
122 : break;
123 : }
124 :
125 852234 : case FluidStatePhaseEnum::LIQUID:
126 : {
127 : // Calculate the liquid properties
128 852234 : const ADReal liquid_pressure = pressure - _pc.capillaryPressure(1.0, qp);
129 852234 : liquidProperties(liquid_pressure, temperature, fsp);
130 :
131 : break;
132 : }
133 :
134 51963 : case FluidStatePhaseEnum::TWOPHASE:
135 : {
136 : // Calculate the gas and liquid properties in the two phase region
137 51963 : twoPhaseProperties(pressure, temperature, Z, qp, fsp);
138 :
139 : break;
140 : }
141 : }
142 :
143 : // Liquid saturations can now be set
144 1824890 : liquid.saturation = 1.0 - gas.saturation;
145 :
146 : // Save pressures to FluidStateProperties object
147 912445 : gas.pressure = pressure;
148 1824890 : liquid.pressure = pressure - _pc.capillaryPressure(liquid.saturation, qp);
149 912445 : }
150 :
151 : void
152 912481 : PorousFlowWaterNCG::massFractions(const ADReal & pressure,
153 : const ADReal & temperature,
154 : const ADReal & Z,
155 : FluidStatePhaseEnum & phase_state,
156 : std::vector<FluidStateProperties> & fsp) const
157 : {
158 912481 : FluidStateProperties & liquid = fsp[_aqueous_phase_number];
159 912481 : FluidStateProperties & gas = fsp[_gas_phase_number];
160 :
161 : // Equilibrium mass fraction of NCG in liquid and H2O in gas phases
162 : ADReal Xncg, Yh2o;
163 912481 : equilibriumMassFractions(pressure, temperature, Xncg, Yh2o);
164 :
165 912481 : ADReal Yncg = 1.0 - Yh2o;
166 :
167 : // Determine which phases are present based on the value of Z
168 912481 : phaseState(Z.value(), Xncg.value(), Yncg.value(), phase_state);
169 :
170 : // The equilibrium mass fractions calculated above are only correct in the two phase
171 : // state. If only liquid or gas phases are present, the mass fractions are given by
172 : // the total mass fraction Z.
173 912481 : ADReal Xh2o = 0.0;
174 :
175 912481 : switch (phase_state)
176 : {
177 852238 : case FluidStatePhaseEnum::LIQUID:
178 : {
179 852238 : Xncg = Z;
180 852238 : Yncg = 0.0;
181 1704476 : Xh2o = 1.0 - Z;
182 852238 : Yh2o = 0.0;
183 852238 : break;
184 : }
185 :
186 8253 : case FluidStatePhaseEnum::GAS:
187 : {
188 8253 : Xncg = 0.0;
189 8253 : Yncg = Z;
190 16506 : Yh2o = 1.0 - Z;
191 8253 : break;
192 : }
193 :
194 51990 : case FluidStatePhaseEnum::TWOPHASE:
195 : {
196 : // Keep equilibrium mass fractions
197 103980 : Xh2o = 1.0 - Xncg;
198 51990 : break;
199 : }
200 : }
201 :
202 : // Save the mass fractions in the FluidStateMassFractions object
203 912481 : liquid.mass_fraction[_aqueous_fluid_component] = Xh2o;
204 912481 : liquid.mass_fraction[_gas_fluid_component] = Xncg;
205 912481 : gas.mass_fraction[_aqueous_fluid_component] = Yh2o;
206 912481 : gas.mass_fraction[_gas_fluid_component] = Yncg;
207 912481 : }
208 :
209 : void
210 60228 : PorousFlowWaterNCG::gasProperties(const ADReal & pressure,
211 : const ADReal & temperature,
212 : std::vector<FluidStateProperties> & fsp) const
213 : {
214 60228 : FluidStateProperties & liquid = fsp[_aqueous_phase_number];
215 60228 : FluidStateProperties & gas = fsp[_gas_phase_number];
216 :
217 60228 : const ADReal psat = _water_fp.vaporPressure(temperature);
218 :
219 60228 : const ADReal Yncg = gas.mass_fraction[_gas_fluid_component];
220 60228 : const ADReal Xncg = liquid.mass_fraction[_gas_fluid_component];
221 :
222 : // NCG density, viscosity and enthalpy calculated using partial pressure
223 : // Yncg * gas_poreressure (Dalton's law)
224 : ADReal ncg_density, ncg_viscosity;
225 120456 : _ncg_fp.rho_mu_from_p_T(Yncg * pressure, temperature, ncg_density, ncg_viscosity);
226 120456 : ADReal ncg_enthalpy = _ncg_fp.h_from_p_T(Yncg * pressure, temperature);
227 :
228 : // Vapor density, viscosity and enthalpy calculated using partial pressure
229 : // X1 * psat (Raoult's law)
230 : ADReal vapor_density, vapor_viscosity;
231 :
232 120456 : _water_fp.rho_mu_from_p_T((1.0 - Xncg) * psat, temperature, vapor_density, vapor_viscosity);
233 120456 : ADReal vapor_enthalpy = _water_fp.h_from_p_T((1.0 - Xncg) * psat, temperature);
234 :
235 : // Density is just the sum of individual component densities
236 60228 : gas.density = ncg_density + vapor_density;
237 :
238 : // Viscosity of the gas phase is a weighted sum of the individual viscosities
239 120456 : gas.viscosity = Yncg * ncg_viscosity + (1.0 - Yncg) * vapor_viscosity;
240 :
241 : // Enthalpy of the gas phase is a weighted sum of the individual enthalpies
242 120456 : gas.enthalpy = Yncg * ncg_enthalpy + (1.0 - Yncg) * vapor_enthalpy;
243 :
244 : // Internal energy of the gas phase (e = h - pv)
245 : mooseAssert(gas.density.value() > 0.0, "Gas density must be greater than zero");
246 60228 : gas.internal_energy = gas.enthalpy - pressure / gas.density;
247 60228 : }
248 :
249 : void
250 904215 : PorousFlowWaterNCG::liquidProperties(const ADReal & pressure,
251 : const ADReal & temperature,
252 : std::vector<FluidStateProperties> & fsp) const
253 : {
254 904215 : FluidStateProperties & liquid = fsp[_aqueous_phase_number];
255 :
256 : // Calculate liquid density and viscosity if in the two phase or single phase
257 : // liquid region, assuming they are not affected by the presence of dissolved
258 : // NCG. Note: the (small) contribution due to derivative of capillary pressure
259 : // wrt pressure (using the chain rule) is not implemented.
260 : ADReal liquid_density, liquid_viscosity;
261 904215 : _water_fp.rho_mu_from_p_T(pressure, temperature, liquid_density, liquid_viscosity);
262 :
263 904215 : liquid.density = liquid_density;
264 904215 : liquid.viscosity = liquid_viscosity;
265 :
266 : // Enthalpy does include a contribution due to the enthalpy of dissolution
267 904215 : const ADReal hdis = enthalpyOfDissolution(temperature);
268 :
269 904215 : const ADReal water_enthalpy = _water_fp.h_from_p_T(pressure, temperature);
270 904215 : const ADReal ncg_enthalpy = _ncg_fp.h_from_p_T(pressure, temperature);
271 :
272 904215 : const ADReal Xncg = liquid.mass_fraction[_gas_fluid_component];
273 1808430 : liquid.enthalpy = (1.0 - Xncg) * water_enthalpy + Xncg * (ncg_enthalpy + hdis);
274 :
275 : // Internal energy of the liquid phase (e = h - pv)
276 : mooseAssert(liquid.density.value() > 0.0, "Liquid density must be greater than zero");
277 904215 : liquid.internal_energy = liquid.enthalpy - pressure / liquid.density;
278 904215 : }
279 :
280 : ADReal
281 52012 : PorousFlowWaterNCG::liquidDensity(const ADReal & pressure, const ADReal & temperature) const
282 : {
283 52012 : return _water_fp.rho_from_p_T(pressure, temperature);
284 : }
285 :
286 : ADReal
287 52012 : PorousFlowWaterNCG::gasDensity(const ADReal & pressure,
288 : const ADReal & temperature,
289 : std::vector<FluidStateProperties> & fsp) const
290 : {
291 52012 : auto & liquid = fsp[_aqueous_phase_number];
292 52012 : auto & gas = fsp[_gas_phase_number];
293 :
294 52012 : ADReal psat = _water_fp.vaporPressure(temperature);
295 :
296 52012 : const ADReal Yncg = gas.mass_fraction[_gas_fluid_component];
297 52012 : const ADReal Xncg = liquid.mass_fraction[_gas_fluid_component];
298 :
299 104024 : ADReal ncg_density = _ncg_fp.rho_from_p_T(Yncg * pressure, temperature);
300 104024 : ADReal vapor_density = _water_fp.rho_from_p_T((1.0 - Xncg) * psat, temperature);
301 :
302 : // Density is just the sum of individual component densities
303 52012 : return ncg_density + vapor_density;
304 : }
305 :
306 : ADReal
307 51974 : PorousFlowWaterNCG::saturation(const ADReal & pressure,
308 : const ADReal & temperature,
309 : const ADReal & Z,
310 : std::vector<FluidStateProperties> & fsp) const
311 : {
312 51974 : auto & gas = fsp[_gas_phase_number];
313 51974 : auto & liquid = fsp[_aqueous_fluid_component];
314 :
315 : // Approximate liquid density as saturation isn't known yet, by using the gas
316 : // pressure rather than the liquid pressure. This does result in a small error
317 : // in the calculated saturation, but this is below the error associated with
318 : // the correlations. A more accurate saturation could be found iteraviely,
319 : // at the cost of increased computational expense
320 :
321 : // The gas and liquid densities
322 51974 : const ADReal gas_density = gasDensity(pressure, temperature, fsp);
323 51974 : const ADReal liquid_density = liquidDensity(pressure, temperature);
324 :
325 : // Set mass equilibrium constants used in the calculation of vapor mass fraction
326 51974 : const ADReal Xncg = liquid.mass_fraction[_gas_fluid_component];
327 51974 : const ADReal Yncg = gas.mass_fraction[_gas_fluid_component];
328 :
329 : const ADReal K0 = Yncg / Xncg;
330 103948 : const ADReal K1 = (1.0 - Yncg) / (1.0 - Xncg);
331 51974 : const ADReal vapor_mass_fraction = vaporMassFraction(Z, K0, K1);
332 :
333 : // The gas saturation in the two phase case
334 51974 : const ADReal saturation = vapor_mass_fraction * liquid_density /
335 51974 : (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
336 :
337 51974 : return saturation;
338 : }
339 :
340 : void
341 51964 : PorousFlowWaterNCG::twoPhaseProperties(const ADReal & pressure,
342 : const ADReal & temperature,
343 : const ADReal & Z,
344 : unsigned int qp,
345 : std::vector<FluidStateProperties> & fsp) const
346 : {
347 51964 : auto & gas = fsp[_gas_phase_number];
348 :
349 : // Calculate all of the gas phase properties, as these don't depend on saturation
350 51964 : gasProperties(pressure, temperature, fsp);
351 :
352 : // The gas saturation in the two phase case
353 51964 : gas.saturation = saturation(pressure, temperature, Z, fsp);
354 :
355 : // The liquid pressure and properties can now be calculated
356 103928 : const ADReal liquid_pressure = pressure - _pc.capillaryPressure(1.0 - gas.saturation, qp);
357 51964 : liquidProperties(liquid_pressure, temperature, fsp);
358 51964 : }
359 :
360 : void
361 912524 : PorousFlowWaterNCG::equilibriumMassFractions(const ADReal & pressure,
362 : const ADReal & temperature,
363 : ADReal & Xncg,
364 : ADReal & Yh2o) const
365 : {
366 : // Equilibrium constants for each component (Henry's law for the NCG
367 : // component, and Raoult's law for water).
368 912524 : const ADReal Kh = _water97_fp.henryConstant(temperature, _ncg_henry);
369 912524 : const ADReal psat = _water_fp.vaporPressure(temperature);
370 :
371 : const ADReal Kncg = Kh / pressure;
372 : const ADReal Kh2o = psat / pressure;
373 :
374 : // The mole fractions for the NCG component in the two component
375 : // case can be expressed in terms of the equilibrium constants only
376 1825048 : const ADReal xncg = (1.0 - Kh2o) / (Kncg - Kh2o);
377 : const ADReal yncg = Kncg * xncg;
378 :
379 : // Convert mole fractions to mass fractions
380 912524 : Xncg = moleFractionToMassFraction(xncg);
381 1825048 : Yh2o = 1.0 - moleFractionToMassFraction(yncg);
382 912524 : }
383 :
384 : ADReal
385 1825048 : PorousFlowWaterNCG::moleFractionToMassFraction(const ADReal & xmol) const
386 : {
387 5475144 : return xmol * _Mncg / (xmol * _Mncg + (1.0 - xmol) * _Mh2o);
388 : }
389 :
390 : void
391 912493 : PorousFlowWaterNCG::checkVariables(Real temperature) const
392 : {
393 : // Check whether the input temperature is within the region of validity of this equation
394 : // of state (T_triple <= T <= T_critical)
395 912493 : if (temperature < _water_triple_temperature || temperature > _water_critical_temperature)
396 33 : mooseException(name() + ": temperature " + Moose::stringify(temperature) +
397 : " is outside range 273.16 K <= T <= 647.096 K");
398 912482 : }
399 :
400 : ADReal
401 904219 : PorousFlowWaterNCG::enthalpyOfDissolution(const ADReal & temperature) const
402 : {
403 : // Henry's constant
404 904219 : const ADReal Kh = _water97_fp.henryConstant(temperature, _ncg_henry);
405 :
406 1808438 : ADReal hdis = -_R * temperature * temperature * Kh.derivatives()[_Tidx] / Kh / _Mncg;
407 :
408 : // Derivative of enthalpy of dissolution wrt temperature requires the second derivative of
409 : // Henry's constant wrt temperature. For simplicity, approximate this numerically
410 904219 : const Real dT = temperature.value() * 1.0e-8;
411 : const ADReal t2 = temperature + dT;
412 904219 : const ADReal Kh2 = _water97_fp.henryConstant(t2, _ncg_henry);
413 :
414 : const Real dhdis_dT =
415 2712657 : (-_R * t2 * t2 * Kh2.derivatives()[_Tidx] / Kh2 / _Mncg - hdis).value() / dT;
416 :
417 1808438 : hdis.derivatives() = temperature.derivatives() * dhdis_dT;
418 :
419 904219 : return hdis;
420 : }
421 :
422 : Real
423 37 : PorousFlowWaterNCG::totalMassFraction(
424 : Real pressure, Real temperature, Real /* Xnacl */, Real saturation, unsigned int qp) const
425 : {
426 : // Check whether the input temperature is within the region of validity
427 37 : checkVariables(temperature);
428 :
429 : // As we do not require derivatives, we can simply ignore their initialisation
430 37 : const ADReal p = pressure;
431 37 : const ADReal T = temperature;
432 :
433 : // FluidStateProperties data structure
434 37 : std::vector<FluidStateProperties> fsp(_num_phases, FluidStateProperties(_num_components));
435 37 : auto & liquid = fsp[_aqueous_phase_number];
436 37 : auto & gas = fsp[_gas_phase_number];
437 :
438 : // Calculate equilibrium mass fractions in the two-phase state
439 : ADReal Xncg, Yh2o;
440 37 : equilibriumMassFractions(p, T, Xncg, Yh2o);
441 :
442 : // Save the mass fractions in the FluidStateMassFractions object to calculate gas density
443 37 : const ADReal Yncg = 1.0 - Yh2o;
444 74 : liquid.mass_fraction[_aqueous_fluid_component] = 1.0 - Xncg;
445 37 : liquid.mass_fraction[_gas_fluid_component] = Xncg;
446 37 : gas.mass_fraction[_aqueous_fluid_component] = Yh2o;
447 37 : gas.mass_fraction[_gas_fluid_component] = Yncg;
448 :
449 : // Gas density
450 37 : const Real gas_density = gasDensity(p, T, fsp).value();
451 :
452 : // Liquid density
453 37 : const ADReal liquid_pressure = p - _pc.capillaryPressure(1.0 - saturation, qp);
454 37 : const Real liquid_density = liquidDensity(liquid_pressure, T).value();
455 :
456 : // The total mass fraction of ncg (Z) can now be calculated
457 37 : const Real Z = (saturation * gas_density * Yncg.value() +
458 37 : (1.0 - saturation) * liquid_density * Xncg.value()) /
459 37 : (saturation * gas_density + (1.0 - saturation) * liquid_density);
460 :
461 37 : return Z;
462 37 : }
|
