www.mooseframework.org
PorousFlowFluidStateWaterNCG.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
10 
11 template <>
12 InputParameters
14 {
15  InputParameters params = validParams<PorousFlowFluidStateFlashBase>();
16  params.addRequiredParam<UserObjectName>("water_fp", "The name of the user object for water");
17  params.addRequiredParam<UserObjectName>(
18  "gas_fp", "The name of the user object for the non-condensable gas");
19  params.addClassDescription("Fluid state class for water and non-condensable gas");
20  return params;
21 }
22 
24  : PorousFlowFluidStateFlashBase(parameters),
25 
26  _water_fp(getUserObject<Water97FluidProperties>("water_fp")),
27  _ncg_fp(getUserObject<SinglePhaseFluidPropertiesPT>("gas_fp")),
28  _Mh2o(_water_fp.molarMass()),
29  _Mncg(_ncg_fp.molarMass()),
30  _water_triple_temperature(_water_fp.triplePointTemperature()),
31  _water_critical_temperature(_water_fp.criticalTemperature())
32 {
33  // Check that the correct FluidProperties UserObjects have been provided
34  if (_water_fp.fluidName() != "water")
35  mooseError("Only a valid water FluidProperties UserObject can be provided in water_fp");
36 }
37 
38 void
40 {
41  // The FluidProperty objects use temperature in K
42  Real Tk = _temperature[_qp] + _T_c2k;
43 
44  // Check whether the input temperature is within the region of validity of this equation
45  // of state (T_triple <= T <= T_critical)
46  if (Tk < _water_triple_temperature || Tk > _water_critical_temperature)
47  mooseError("PorousFlowFluidStateWaterNCG: Temperature is outside range 273.16 K <= T "
48  "<= 647.096 K");
49 
50  // Equilibrium constants for each component (Henry's law for the NCG
51  // component, and Raoult's law for water).
52  // Note: these are in terms of mole fraction
53  Real psat, dpsat_dT;
54  _water_fp.vaporPressure_dT(Tk, psat, dpsat_dT);
55  Real K0 = _ncg_fp.henryConstant(Tk) / _gas_porepressure[_qp];
56  Real K1 = psat / _gas_porepressure[_qp];
57 
58  // The mole fractions for the NCG component in the two component
59  // case can be expressed in terms of the equilibrium constants only
60  Real x0 = (1.0 - K1) / (K0 - K1);
61  Real y0 = K0 * x0;
62 
63  // Convert mole fractions to mass fractions
64  Real X0 = moleFractionToMassFraction(x0);
65  Real Y0 = moleFractionToMassFraction(y0);
66 
67  // Determine which phases are present based on the value of z
68  bool is_liquid = false;
69  bool is_gas = false;
70  bool is_twophase = false;
71 
72  if ((*_z[0])[_qp] <= X0)
73  {
74  // In this case, there is not enough NCG to form a gas phase,
75  // so only a liquid phase is present
76  is_liquid = true;
77  }
78  else if ((*_z[0])[_qp] > X0 && (*_z[0])[_qp] < Y0)
79  {
80  // Two phases are present
81  is_twophase = true;
82  }
83  else // ((*_z[0])[_qp] >= Y0)
84  {
85  // In this case, there is not enough water to form a liquid
86  // phase, so only a gas phase is present
87  is_gas = true;
88  }
89 
90  // Material must provide the following properties
91  Real vapor_mass_fraction = 0.0;
92  Real gas_saturation, liquid_saturation;
93  Real gas_density, liquid_density;
94  Real gas_viscosity, liquid_viscosity;
95 
96  // And the following derivatives
97  Real dX0_dp = 0.0, dX0_dT = 0.0, dX0_dz = 0.0;
98  Real dY0_dp = 0.0, dY0_dT = 0.0, dY0_dz = 0.0;
99  Real dgas_density_dp = 0.0, dgas_density_dT = 0.0, dgas_density_dz = 0.0;
100  Real dliquid_density_dp = 0.0, dliquid_density_dT = 0.0, dliquid_density_dz = 0.0;
101  Real dgas_viscosity_dp = 0.0, dgas_viscosity_dT = 0.0, dgas_viscosity_dz = 0.0;
102  Real dliquid_viscosity_dp = 0.0, dliquid_viscosity_dT = 0.0, dliquid_viscosity_dz = 0.0;
103 
104  if (is_liquid)
105  {
106  X0 = (*_z[0])[_qp];
107  Y0 = 0.0;
108  gas_saturation = 0.0;
109  dX0_dz = 1.0;
110  gas_density = 0.0;
111  gas_viscosity = 1.0; // To guard against division by 0
112  }
113  else if (is_twophase)
114  {
115  // Two phase are present. Set mass equilibrium constants used in the
116  // calculation of vapor mass fraction
117  std::vector<Real> Ki(_num_components);
118  Ki[0] = Y0 / X0;
119  Ki[1] = (1.0 - Y0) / (1.0 - X0);
120  vapor_mass_fraction = vaporMassFraction(Ki);
121 
122  // Derivatives of mass fractions wrt PorousFlow variables
123  Real Kh, dKh_dT;
124  _ncg_fp.henryConstant_dT(Tk, Kh, dKh_dT);
125  Real dK0_dp = -Kh / _gas_porepressure[_qp] / _gas_porepressure[_qp];
126  Real dK0_dT = dKh_dT / _gas_porepressure[_qp];
127 
128  Real dK1_dp = -psat / _gas_porepressure[_qp] / _gas_porepressure[_qp];
129  Real dK1_dT = dpsat_dT / _gas_porepressure[_qp];
130 
131  Real dx0_dp = ((K1 - 1.0) * dK0_dp + (1 - K0) * dK1_dp) / (K0 - K1) / (K0 - K1);
132  Real dy0_dp = x0 * dK0_dp + K0 * dx0_dp;
133  Real dx0_dT = ((K1 - 1.0) * dK0_dT + (1 - K0) * dK1_dT) / (K0 - K1) / (K0 - K1);
134  Real dy0_dT = x0 * dK0_dT + K0 * dx0_dT;
135 
136  Real dX0_dx0 =
137  _Mncg * _Mh2o / (x0 * _Mncg + (1.0 - x0) * _Mh2o) / (x0 * _Mncg + (1.0 - x0) * _Mh2o);
138  Real dY0_dy0 =
139  _Mncg * _Mh2o / (y0 * _Mncg + (1.0 - y0) * _Mh2o) / (y0 * _Mncg + (1.0 - y0) * _Mh2o);
140 
141  dX0_dp = dX0_dx0 * dx0_dp;
142  dX0_dT = dX0_dx0 * dx0_dT;
143  dY0_dp = dY0_dy0 * dy0_dp;
144  dY0_dT = dY0_dy0 * dy0_dT;
145  }
146  else // if (is_gas)
147  {
148  X0 = 0.0;
149  Y0 = (*_z[0])[_qp];
150  vapor_mass_fraction = 1.0;
151  gas_saturation = 1.0;
152  dY0_dz = 1.0;
153  liquid_density = 0.0;
154  liquid_viscosity = 1.0; // To guard against division by 0
155  }
156 
157  // Calculate the gas density and viscosity in the single phase gas or two
158  // phase region
159  if (is_gas || is_twophase)
160  {
161  Real ncg_density, dncg_density_dp, dncg_density_dT;
162  Real vapor_density, dvapor_density_dp, dvapor_density_dT;
163  // NCG density calculated using partial pressure Y0 * gas_poreressure (Dalton's law)
164  _ncg_fp.rho_dpT(Y0 * _gas_porepressure[_qp], Tk, ncg_density, dncg_density_dp, dncg_density_dT);
165  // Vapor density calculated using partial pressure X1 * psat (Raoult's law)
166  _water_fp.rho_dpT((1.0 - X0) * psat, Tk, vapor_density, dvapor_density_dp, dvapor_density_dT);
167 
168  // The derivatives wrt pressure above must be multiplied by the derivative of the pressure
169  // variable using the chain rule
170  gas_density = ncg_density + vapor_density;
171  dgas_density_dp = (Y0 + dY0_dp * _gas_porepressure[_qp]) * dncg_density_dp -
172  dX0_dp * psat * dvapor_density_dp;
173  dgas_density_dT = dncg_density_dT + dvapor_density_dT;
174 
175  Real ncg_viscosity, dncg_viscosity_drho, dncg_viscosity_dT;
176  Real vapor_viscosity, dvapor_viscosity_drho, dvapor_viscosity_dT;
178  ncg_density, Tk, dncg_density_dT, ncg_viscosity, dncg_viscosity_drho, dncg_viscosity_dT);
179  _water_fp.mu_drhoT_from_rho_T(vapor_density,
180  Tk,
181  dvapor_density_dT,
182  vapor_viscosity,
183  dvapor_viscosity_drho,
184  dvapor_viscosity_dT);
185 
186  // Assume that the viscosity of the gas phase is a weighted sum of the
187  // individual viscosities
188  gas_viscosity = Y0 * ncg_viscosity + (1.0 - Y0) * vapor_viscosity;
189  dgas_viscosity_dp =
190  dY0_dp * (ncg_viscosity - vapor_viscosity) +
191  Y0 * (Y0 + dY0_dp * _gas_porepressure[_qp]) * dncg_viscosity_drho * dncg_density_dp -
192  dX0_dp * psat * (1.0 - Y0) * dvapor_viscosity_drho * dvapor_density_dp;
193  dgas_viscosity_dT = dY0_dT * (ncg_viscosity - vapor_viscosity) + Y0 * dncg_viscosity_dT +
194  (1.0 - Y0) * dvapor_density_dT;
195 
196  // Also calculate derivatives wrt z in the gas phase (these are 0 in the two phase region)
197  if (is_gas)
198  {
199  dgas_density_dz =
200  dY0_dz * _gas_porepressure[_qp] * dncg_density_dp - dX0_dz * psat * dvapor_density_dp;
201 
202  dgas_viscosity_dz =
203  dY0_dz * (ncg_viscosity - vapor_viscosity) +
204  Y0 * dncg_viscosity_drho * dncg_density_dp * dY0_dz * _gas_porepressure[_qp] -
205  dX0_dz * psat * (1.0 - Y0) * dvapor_viscosity_drho * dvapor_density_dp;
206  }
207  }
208 
209  // Calculate the saturation in the two phase case using the vapor mass fraction
210  if (is_twophase)
211  {
212  // Liquid density is approximated by the water density.
213  // Use old value of gas saturation to estimate liquid saturation
214  Real liqsat = 1.0;
215  if (!_is_initqp)
216  liqsat -= _saturation_old[_qp][_gas_phase_number];
217  liquid_density = _water_fp.rho(_gas_porepressure[_qp] + _pc_uo.capillaryPressure(liqsat), Tk);
218 
219  // The gas saturation in the two phase case
220  gas_saturation = vapor_mass_fraction * liquid_density /
221  (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
222  }
223 
224  // Calculate the saturations and pressures for each phase
225  liquid_saturation = 1.0 - gas_saturation;
226  Real liquid_porepressure = _gas_porepressure[_qp] - _pc_uo.capillaryPressure(liquid_saturation);
227 
228  // Calculate liquid density and viscosity if in the two phase or single phase
229  // liquid region, assuming they are not affected by the presence of dissolved
230  // NCG. Note: the (small) contribution due to derivative of capillary pressure
231  // wrt pressure (using the chain rule) is not implemented.
232  if (is_liquid || is_twophase)
233  {
234  Real dliquid_viscosity_drho;
236  liquid_porepressure, Tk, liquid_density, dliquid_density_dp, dliquid_density_dT);
237  _water_fp.mu_drhoT_from_rho_T(liquid_density,
238  Tk,
239  dliquid_density_dT,
240  liquid_viscosity,
241  dliquid_viscosity_drho,
242  dliquid_viscosity_dT);
243 
244  // The derivative of viscosity wrt pressure is given by the chain rule
245  dliquid_viscosity_dp = dliquid_viscosity_drho * dliquid_density_dp;
246  }
247 
248  // Save properties in FluidStateProperties vector
249  auto & liquid = _fsp[_aqueous_phase_number];
250  auto & gas = _fsp[_gas_phase_number];
251 
252  liquid.saturation = liquid_saturation;
253  gas.saturation = gas_saturation;
254  liquid.pressure = liquid_porepressure;
255  gas.pressure = _gas_porepressure[_qp];
256 
257  liquid.mass_fraction[_aqueous_fluid_component] = 1.0 - X0;
258  liquid.mass_fraction[_gas_fluid_component] = X0;
259  gas.mass_fraction[_aqueous_fluid_component] = 1.0 - Y0;
260  gas.mass_fraction[_gas_fluid_component] = Y0;
261 
262  liquid.fluid_density = liquid_density;
263  gas.fluid_density = gas_density;
264  liquid.fluid_viscosity = liquid_viscosity;
265  gas.fluid_viscosity = gas_viscosity;
266 
267  // Derivatives wrt PorousFlow variables are required by the kernels
268  // Note: these don't need to be stateful so don't calculate them in
269  // initQpStatefulProperties
270  if (!_is_initqp)
271  {
272  // Derivative of gas saturation wrt variables
273  Real ds_dp = 0.0, ds_dT = 0.0, ds_dz = 0.0;
274  if (is_twophase)
275  {
276  K0 = Y0 / X0;
277  K1 = (1.0 - Y0) / (1.0 - X0);
278  Real dv_dz = (K1 - K0) / ((K0 - 1.0) * (K1 - 1.0));
279  ds_dz = gas_density * liquid_density * dv_dz +
280  vapor_mass_fraction * (1.0 - vapor_mass_fraction) *
281  (gas_density * dliquid_density_dz - dgas_density_dz * liquid_density);
282  ds_dz /= (gas_density + vapor_mass_fraction * (liquid_density - gas_density)) *
283  (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
284 
285  Real Kh, dKh_dT;
286  _ncg_fp.henryConstant_dT(Tk, Kh, dKh_dT);
287  Real dK0_dp = -Kh / _gas_porepressure[_qp] / _gas_porepressure[_qp];
288  Real dK0_dT = dKh_dT / _gas_porepressure[_qp];
289 
290  Real dK1_dp = -psat / _gas_porepressure[_qp] / _gas_porepressure[_qp];
291  Real dK1_dT = dpsat_dT / _gas_porepressure[_qp];
292 
293  Real dv_dp = (*_z[0])[_qp] * dK1_dp / (K1 - 1.0) / (K1 - 1.0) +
294  (1.0 - (*_z[0])[_qp]) * dK0_dp / (K0 - 1.0) / (K0 - 1.0);
295  ds_dp = gas_density * liquid_density * dv_dp +
296  vapor_mass_fraction * (1.0 - vapor_mass_fraction) *
297  (gas_density * dliquid_density_dp - dgas_density_dp * liquid_density);
298 
299  ds_dp /= (gas_density + vapor_mass_fraction * (liquid_density - gas_density)) *
300  (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
301 
302  Real dv_dT = (*_z[0])[_qp] * dK1_dT / (K1 - 1.0) / (K1 - 1.0) +
303  (1.0 - (*_z[0])[_qp]) * dK0_dT / (K0 - 1.0) / (K0 - 1.0);
304  ds_dT = gas_density * liquid_density * dv_dT +
305  vapor_mass_fraction * (1.0 - vapor_mass_fraction) *
306  (gas_density * dliquid_density_dT - dgas_density_dT * liquid_density);
307 
308  ds_dT /= (gas_density + vapor_mass_fraction * (liquid_density - gas_density)) *
309  (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
310  }
311 
312  liquid.dsaturation_dp = -ds_dp;
313  liquid.dsaturation_dT = -ds_dT;
314  liquid.dsaturation_dz = -ds_dz;
315  gas.dsaturation_dp = ds_dp;
316  gas.dsaturation_dT = ds_dT;
317  gas.dsaturation_dz = ds_dz;
318 
319  liquid.dmass_fraction_dp[_aqueous_fluid_component] = -dX0_dp;
320  liquid.dmass_fraction_dp[_gas_fluid_component] = dX0_dp;
321  liquid.dmass_fraction_dT[_aqueous_fluid_component] = -dX0_dT;
322  liquid.dmass_fraction_dT[_gas_fluid_component] = dX0_dT;
323  liquid.dmass_fraction_dz[_aqueous_fluid_component] = -dX0_dz;
324  liquid.dmass_fraction_dz[_gas_fluid_component] = dX0_dz;
325 
326  gas.dmass_fraction_dp[_aqueous_fluid_component] = -dY0_dp;
327  gas.dmass_fraction_dp[_gas_fluid_component] = dY0_dp;
328  gas.dmass_fraction_dT[_aqueous_fluid_component] = -dY0_dT;
329  gas.dmass_fraction_dT[_gas_fluid_component] = dY0_dT;
330  gas.dmass_fraction_dz[_aqueous_fluid_component] = -dY0_dz;
331  gas.dmass_fraction_dz[_gas_fluid_component] = dY0_dz;
332 
333  liquid.dfluid_density_dp = dliquid_density_dp;
334  liquid.dfluid_density_dT = dliquid_density_dT;
335  liquid.dfluid_density_dz = dliquid_density_dz;
336  gas.dfluid_density_dp = dgas_density_dp;
337  gas.dfluid_density_dT = dgas_density_dT;
338  gas.dfluid_density_dz = dgas_density_dz;
339 
340  liquid.dfluid_viscosity_dp = dliquid_viscosity_dp;
341  liquid.dfluid_viscosity_dT = dliquid_viscosity_dT;
342  liquid.dfluid_viscosity_dz = dliquid_viscosity_dz;
343  gas.dfluid_viscosity_dp = dgas_viscosity_dp;
344  gas.dfluid_viscosity_dT = dgas_viscosity_dT;
345  gas.dfluid_viscosity_dz = dgas_viscosity_dz;
346  }
347 }
348 
349 Real
351 {
352  return -_R * temperature * temperature * _Mncg * dKh_dT / Kh;
353 }
354 
355 Real
357 {
358  return xmol * _Mncg / (xmol * _Mncg + (1.0 - xmol) * _Mh2o);
359 }
const Real _Mncg
Molar mass of non-condensable gas (kg/mol)
const Water97FluidProperties & _water_fp
Fluid properties UserObject for water.
const SinglePhaseFluidPropertiesPT & _ncg_fp
Fluid properties UserObject for the NCG.
const MaterialProperty< Real > & _temperature
Temperature.
virtual std::string fluidName() const override
Fluid name.
const Real _water_critical_temperature
Critical temperature of water (K)
void vaporPressure_dT(Real temperature, Real &psat, Real &dpsat_dT) const
Saturation pressure as a function of temperature and derivative wrt temperature.
virtual void rho_dpT(Real pressure, Real temperature, Real &rho, Real &drho_dp, Real &drho_dT) const =0
Density and its derivatives wrt pressure and temperature.
const Real _Mh2o
Molar mass of water (kg/mol)
Real moleFractionToMassFraction(Real xmol) const
Convert mole fraction to mass fraction.
const std::string temperature
Definition: NS.h:25
const VariableValue & _gas_porepressure
Porepressure.
virtual Real rho(Real pressure, Real temperature) const override
Density.
const MaterialProperty< std::vector< Real > > & _saturation_old
Old value of saturation.
Common class for single phase fluid properties using a pressure and temperature formulation.
virtual Real henryConstant(Real temperature) const =0
Henry&#39;s law constant for dissolution in water.
virtual void henryConstant_dT(Real temperature, Real &Kh, Real &dKh_dT) const =0
Henry&#39;s law constant for dissolution in water and derivative wrt temperature.
virtual Real capillaryPressure(Real saturation) const
Capillary pressure is calculated as a function of true saturation.
const unsigned int _gas_fluid_component
Fluid component number of the gas phase.
const PorousFlowCapillaryPressure & _pc_uo
Capillary pressure UserObject.
virtual void mu_drhoT_from_rho_T(Real density, Real temperature, Real ddensity_dT, Real &mu, Real &dmu_drho, Real &dmu_dT) const override
Dynamic viscosity and its derivatives wrt density and temperature.
std::vector< FluidStateProperties > _fsp
FluidStateProperties data structure.
Real enthalpyOfDissolution(Real temperature, Real Kh, Real dKh_dT) const
Enthalpy of dissolution of NCG in water calculated using Henry&#39;s constant From Himmelblau, Partial molal heats and entropies of solution for gases dissolved in water from the freezing to the near critical point, J.
const unsigned int _aqueous_phase_number
Phase number of the aqueous phase.
virtual void mu_drhoT_from_rho_T(Real density, Real temperature, Real ddensity_dT, Real &mu, Real &dmu_drho, Real &dmu_dT) const =0
Dynamic viscosity and its derivatives wrt density and temperature.
InputParameters validParams< PorousFlowFluidStateWaterNCG >()
const Real _R
Universal gas constant (J/mol/K)
const unsigned int _gas_phase_number
Phase number of the gas phase.
bool _is_initqp
Flag to indicate whether to calculate stateful properties.
virtual void rho_dpT(Real pressure, Real temperature, Real &rho, Real &drho_dp, Real &drho_dT) const override
Density and its derivatives wrt pressure and temperature.
const unsigned int _aqueous_fluid_component
Fluid component number of the aqueous component.
InputParameters validParams< PorousFlowFluidStateFlashBase >()
const Real _T_c2k
Conversion from degrees Celsius to degrees Kelvin.
const unsigned int _num_components
Number of components.
Water (H2O) fluid properties as a function of pressure (Pa) and temperature (K) from IAPWS-IF97: Revi...
virtual void thermophysicalProperties() override
Calculates all required thermophysical properties and derivatives for each phase and fluid component...
PorousFlowFluidStateWaterNCG(const InputParameters &parameters)
std::vector< const VariableValue * > _z
Total mass fraction(s) of the gas component(s) summed over all phases.
Base class for fluid states using a persistent set of primary variables for the mutliphase, multicomponent case.
Real vaporMassFraction(std::vector< Real > &Ki) const
Solves Rachford-Rice equation to provide vapor mass fraction.