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 252 : PorousFlowWaterNCG::validParams()
19 : {
20 252 : InputParameters params = PorousFlowFluidStateMultiComponentBase::validParams();
21 504 : params.addRequiredParam<UserObjectName>("water_fp", "The name of the user object for water");
22 504 : params.addRequiredParam<UserObjectName>(
23 : "gas_fp", "The name of the user object for the non-condensable gas");
24 252 : params.addClassDescription("Fluid state class for water and non-condensable gas");
25 252 : return params;
26 0 : }
27 :
28 126 : PorousFlowWaterNCG::PorousFlowWaterNCG(const InputParameters & parameters)
29 : : PorousFlowFluidStateMultiComponentBase(parameters),
30 126 : _water_fp(getUserObject<SinglePhaseFluidProperties>("water_fp")),
31 126 : _water97_fp(getUserObject<Water97FluidProperties>("water_fp")),
32 126 : _ncg_fp(getUserObject<SinglePhaseFluidProperties>("gas_fp")),
33 126 : _Mh2o(_water_fp.molarMass()),
34 126 : _Mncg(_ncg_fp.molarMass()),
35 126 : _water_triple_temperature(_water_fp.triplePointTemperature()),
36 126 : _water_critical_temperature(_water_fp.criticalTemperature()),
37 252 : _ncg_henry(_ncg_fp.henryCoefficients())
38 : {
39 : // Check that the correct FluidProperties UserObjects have been provided
40 252 : 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 126 : _num_phases = 2;
45 126 : _num_components = 2;
46 126 : _gas_phase_number = 1 - _aqueous_phase_number;
47 126 : _gas_fluid_component = 1 - _aqueous_fluid_component;
48 :
49 : // Check that _aqueous_phase_number is <= total number of phases
50 126 : 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 126 : 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 126 : _empty_fsp = FluidStateProperties(_num_components);
62 126 : }
63 :
64 : std::string
65 1 : PorousFlowWaterNCG::fluidStateName() const
66 : {
67 1 : return "water-ncg";
68 : }
69 :
70 : void
71 613382 : 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 613382 : ADReal p = pressure;
80 613382 : Moose::derivInsert(p.derivatives(), _pidx, 1.0);
81 613382 : ADReal T = temperature;
82 613382 : Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
83 613382 : ADReal Zncg = Z;
84 613382 : Moose::derivInsert(Zncg.derivatives(), _Zidx, 1.0);
85 : // X is not used, but needed for consistency with PorousFlowFluidState interface
86 613382 : ADReal X = Xnacl;
87 613382 : Moose::derivInsert(X.derivatives(), _Xidx, 1.0);
88 :
89 613382 : thermophysicalProperties(p, T, X, Zncg, qp, fsp);
90 613375 : }
91 :
92 : void
93 613382 : 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 613382 : FluidStateProperties & liquid = fsp[_aqueous_phase_number];
101 613382 : FluidStateProperties & gas = fsp[_gas_phase_number];
102 :
103 : // Check whether the input temperature is within the region of validity
104 613382 : checkVariables(temperature.value());
105 :
106 : // Clear all of the FluidStateProperties data
107 613375 : clearFluidStateProperties(fsp);
108 :
109 : FluidStatePhaseEnum phase_state;
110 613375 : massFractions(pressure, temperature, Z, phase_state, fsp);
111 :
112 613375 : switch (phase_state)
113 : {
114 6544 : case FluidStatePhaseEnum::GAS:
115 : {
116 : // Set the gas saturations
117 6544 : gas.saturation = 1.0;
118 :
119 : // Calculate gas properties
120 6544 : gasProperties(pressure, temperature, fsp);
121 :
122 : break;
123 : }
124 :
125 569643 : case FluidStatePhaseEnum::LIQUID:
126 : {
127 : // Calculate the liquid properties
128 569643 : const ADReal liquid_pressure = pressure - _pc.capillaryPressure(1.0, qp);
129 569643 : liquidProperties(liquid_pressure, temperature, fsp);
130 :
131 : break;
132 : }
133 :
134 37188 : case FluidStatePhaseEnum::TWOPHASE:
135 : {
136 : // Calculate the gas and liquid properties in the two phase region
137 37188 : twoPhaseProperties(pressure, temperature, Z, qp, fsp);
138 :
139 : break;
140 : }
141 : }
142 :
143 : // Liquid saturations can now be set
144 1226750 : liquid.saturation = 1.0 - gas.saturation;
145 :
146 : // Save pressures to FluidStateProperties object
147 613375 : gas.pressure = pressure;
148 1226750 : liquid.pressure = pressure - _pc.capillaryPressure(liquid.saturation, qp);
149 613375 : }
150 :
151 : void
152 613411 : 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 613411 : FluidStateProperties & liquid = fsp[_aqueous_phase_number];
159 613411 : 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 613411 : equilibriumMassFractions(pressure, temperature, Xncg, Yh2o);
164 :
165 613411 : ADReal Yncg = 1.0 - Yh2o;
166 :
167 : // Determine which phases are present based on the value of Z
168 613411 : 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 613411 : ADReal Xh2o = 0.0;
174 :
175 613411 : switch (phase_state)
176 : {
177 569647 : case FluidStatePhaseEnum::LIQUID:
178 : {
179 569647 : Xncg = Z;
180 569647 : Yncg = 0.0;
181 1139294 : Xh2o = 1.0 - Z;
182 569647 : Yh2o = 0.0;
183 569647 : break;
184 : }
185 :
186 6549 : case FluidStatePhaseEnum::GAS:
187 : {
188 6549 : Xncg = 0.0;
189 6549 : Yncg = Z;
190 13098 : Yh2o = 1.0 - Z;
191 6549 : break;
192 : }
193 :
194 37215 : case FluidStatePhaseEnum::TWOPHASE:
195 : {
196 : // Keep equilibrium mass fractions
197 74430 : Xh2o = 1.0 - Xncg;
198 37215 : break;
199 : }
200 : }
201 :
202 : // Save the mass fractions in the FluidStateMassFractions object
203 613411 : liquid.mass_fraction[_aqueous_fluid_component] = Xh2o;
204 613411 : liquid.mass_fraction[_gas_fluid_component] = Xncg;
205 613411 : gas.mass_fraction[_aqueous_fluid_component] = Yh2o;
206 613411 : gas.mass_fraction[_gas_fluid_component] = Yncg;
207 613411 : }
208 :
209 : void
210 43749 : PorousFlowWaterNCG::gasProperties(const ADReal & pressure,
211 : const ADReal & temperature,
212 : std::vector<FluidStateProperties> & fsp) const
213 : {
214 43749 : FluidStateProperties & liquid = fsp[_aqueous_phase_number];
215 43749 : FluidStateProperties & gas = fsp[_gas_phase_number];
216 :
217 43749 : const ADReal psat = _water_fp.vaporPressure(temperature);
218 :
219 43749 : const ADReal Yncg = gas.mass_fraction[_gas_fluid_component];
220 43749 : 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 87498 : _ncg_fp.rho_mu_from_p_T(Yncg * pressure, temperature, ncg_density, ncg_viscosity);
226 87498 : 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 87498 : _water_fp.rho_mu_from_p_T((1.0 - Xncg) * psat, temperature, vapor_density, vapor_viscosity);
233 87498 : 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 43749 : gas.density = ncg_density + vapor_density;
237 :
238 : // Viscosity of the gas phase is a weighted sum of the individual viscosities
239 87498 : 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 87498 : 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 43749 : gas.internal_energy = gas.enthalpy - pressure / gas.density;
247 43749 : }
248 :
249 : void
250 606849 : PorousFlowWaterNCG::liquidProperties(const ADReal & pressure,
251 : const ADReal & temperature,
252 : std::vector<FluidStateProperties> & fsp) const
253 : {
254 606849 : 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 606849 : _water_fp.rho_mu_from_p_T(pressure, temperature, liquid_density, liquid_viscosity);
262 :
263 606849 : liquid.density = liquid_density;
264 606849 : liquid.viscosity = liquid_viscosity;
265 :
266 : // Enthalpy does include a contribution due to the enthalpy of dissolution
267 606849 : const ADReal hdis = enthalpyOfDissolution(temperature);
268 :
269 606849 : const ADReal water_enthalpy = _water_fp.h_from_p_T(pressure, temperature);
270 606849 : const ADReal ncg_enthalpy = _ncg_fp.h_from_p_T(pressure, temperature);
271 :
272 606849 : const ADReal Xncg = liquid.mass_fraction[_gas_fluid_component];
273 1213698 : 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 606849 : liquid.internal_energy = liquid.enthalpy - pressure / liquid.density;
278 606849 : }
279 :
280 : ADReal
281 37225 : PorousFlowWaterNCG::liquidDensity(const ADReal & pressure, const ADReal & temperature) const
282 : {
283 37225 : return _water_fp.rho_from_p_T(pressure, temperature);
284 : }
285 :
286 : ADReal
287 37225 : PorousFlowWaterNCG::gasDensity(const ADReal & pressure,
288 : const ADReal & temperature,
289 : std::vector<FluidStateProperties> & fsp) const
290 : {
291 37225 : auto & liquid = fsp[_aqueous_phase_number];
292 37225 : auto & gas = fsp[_gas_phase_number];
293 :
294 37225 : ADReal psat = _water_fp.vaporPressure(temperature);
295 :
296 37225 : const ADReal Yncg = gas.mass_fraction[_gas_fluid_component];
297 37225 : const ADReal Xncg = liquid.mass_fraction[_gas_fluid_component];
298 :
299 74450 : ADReal ncg_density = _ncg_fp.rho_from_p_T(Yncg * pressure, temperature);
300 74450 : 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 37225 : return ncg_density + vapor_density;
304 : }
305 :
306 : ADReal
307 37199 : PorousFlowWaterNCG::saturation(const ADReal & pressure,
308 : const ADReal & temperature,
309 : const ADReal & Z,
310 : std::vector<FluidStateProperties> & fsp) const
311 : {
312 37199 : auto & gas = fsp[_gas_phase_number];
313 37199 : 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 37199 : const ADReal gas_density = gasDensity(pressure, temperature, fsp);
323 37199 : const ADReal liquid_density = liquidDensity(pressure, temperature);
324 :
325 : // Set mass equilibrium constants used in the calculation of vapor mass fraction
326 37199 : const ADReal Xncg = liquid.mass_fraction[_gas_fluid_component];
327 37199 : const ADReal Yncg = gas.mass_fraction[_gas_fluid_component];
328 :
329 : const ADReal K0 = Yncg / Xncg;
330 74398 : const ADReal K1 = (1.0 - Yncg) / (1.0 - Xncg);
331 37199 : const ADReal vapor_mass_fraction = vaporMassFraction(Z, K0, K1);
332 :
333 : // The gas saturation in the two phase case
334 37199 : const ADReal saturation = vapor_mass_fraction * liquid_density /
335 37199 : (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
336 :
337 37199 : return saturation;
338 : }
339 :
340 : void
341 37189 : 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 37189 : auto & gas = fsp[_gas_phase_number];
348 :
349 : // Calculate all of the gas phase properties, as these don't depend on saturation
350 37189 : gasProperties(pressure, temperature, fsp);
351 :
352 : // The gas saturation in the two phase case
353 37189 : gas.saturation = saturation(pressure, temperature, Z, fsp);
354 :
355 : // The liquid pressure and properties can now be calculated
356 74378 : const ADReal liquid_pressure = pressure - _pc.capillaryPressure(1.0 - gas.saturation, qp);
357 37189 : liquidProperties(liquid_pressure, temperature, fsp);
358 37189 : }
359 :
360 : void
361 613442 : 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 613442 : const ADReal Kh = _water97_fp.henryConstant(temperature, _ncg_henry);
369 613442 : 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 1226884 : const ADReal xncg = (1.0 - Kh2o) / (Kncg - Kh2o);
377 : const ADReal yncg = Kncg * xncg;
378 :
379 : // Convert mole fractions to mass fractions
380 613442 : Xncg = moleFractionToMassFraction(xncg);
381 1226884 : Yh2o = 1.0 - moleFractionToMassFraction(yncg);
382 613442 : }
383 :
384 : ADReal
385 1226884 : PorousFlowWaterNCG::moleFractionToMassFraction(const ADReal & xmol) const
386 : {
387 3680652 : return xmol * _Mncg / (xmol * _Mncg + (1.0 - xmol) * _Mh2o);
388 : }
389 :
390 : void
391 613407 : 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 613407 : if (temperature < _water_triple_temperature || temperature > _water_critical_temperature)
396 21 : mooseException(name() + ": temperature " + Moose::stringify(temperature) +
397 : " is outside range 273.16 K <= T <= 647.096 K");
398 613400 : }
399 :
400 : ADReal
401 606853 : PorousFlowWaterNCG::enthalpyOfDissolution(const ADReal & temperature) const
402 : {
403 : // Henry's constant
404 606853 : const ADReal Kh = _water97_fp.henryConstant(temperature, _ncg_henry);
405 :
406 1213706 : 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 606853 : const Real dT = temperature.value() * 1.0e-8;
411 : const ADReal t2 = temperature + dT;
412 606853 : const ADReal Kh2 = _water97_fp.henryConstant(t2, _ncg_henry);
413 :
414 : const Real dhdis_dT =
415 1820559 : (-_R * t2 * t2 * Kh2.derivatives()[_Tidx] / Kh2 / _Mncg - hdis).value() / dT;
416 :
417 1213706 : hdis.derivatives() = temperature.derivatives() * dhdis_dT;
418 :
419 606853 : return hdis;
420 : }
421 :
422 : Real
423 25 : 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 25 : checkVariables(temperature);
428 :
429 : // As we do not require derivatives, we can simply ignore their initialisation
430 25 : const ADReal p = pressure;
431 25 : const ADReal T = temperature;
432 :
433 : // FluidStateProperties data structure
434 25 : std::vector<FluidStateProperties> fsp(_num_phases, FluidStateProperties(_num_components));
435 25 : auto & liquid = fsp[_aqueous_phase_number];
436 25 : auto & gas = fsp[_gas_phase_number];
437 :
438 : // Calculate equilibrium mass fractions in the two-phase state
439 : ADReal Xncg, Yh2o;
440 25 : equilibriumMassFractions(p, T, Xncg, Yh2o);
441 :
442 : // Save the mass fractions in the FluidStateMassFractions object to calculate gas density
443 25 : const ADReal Yncg = 1.0 - Yh2o;
444 50 : liquid.mass_fraction[_aqueous_fluid_component] = 1.0 - Xncg;
445 25 : liquid.mass_fraction[_gas_fluid_component] = Xncg;
446 25 : gas.mass_fraction[_aqueous_fluid_component] = Yh2o;
447 25 : gas.mass_fraction[_gas_fluid_component] = Yncg;
448 :
449 : // Gas density
450 25 : const Real gas_density = gasDensity(p, T, fsp).value();
451 :
452 : // Liquid density
453 25 : const ADReal liquid_pressure = p - _pc.capillaryPressure(1.0 - saturation, qp);
454 25 : const Real liquid_density = liquidDensity(liquid_pressure, T).value();
455 :
456 : // The total mass fraction of ncg (Z) can now be calculated
457 25 : const Real Z = (saturation * gas_density * Yncg.value() +
458 25 : (1.0 - saturation) * liquid_density * Xncg.value()) /
459 25 : (saturation * gas_density + (1.0 - saturation) * liquid_density);
460 :
461 25 : return Z;
462 25 : }
|
