www.mooseframework.org
PorousFlowFluidStateFlashBase.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<PorousFlowVariableBase>();
16  params.addRequiredCoupledVar("gas_porepressure",
17  "Variable that is the porepressure of the gas phase");
18  params.addRequiredCoupledVar("z", "Total mass fraction of component i summed over all phases");
19  params.addParam<unsigned int>("liquid_phase_number", 0, "The phase number of the liquid phase");
20  params.addParam<unsigned int>(
21  "liquid_fluid_component", 0, "The fluid component number of the liquid phase");
22  MooseEnum unit_choice("Kelvin=0 Celsius=1", "Kelvin");
23  params.addParam<MooseEnum>(
24  "temperature_unit", unit_choice, "The unit of the temperature variable");
25  params.addRequiredParam<UserObjectName>("capillary_pressure",
26  "Name of the UserObject defining the capillary pressure");
27  params.addClassDescription("Base class for fluid state calculations using persistent primary "
28  "variables and a vapor-liquid flash");
29  return params;
30 }
31 
33  : PorousFlowVariableBase(parameters),
34 
35  _gas_porepressure(_nodal_material ? coupledNodalValue("gas_porepressure")
36  : coupledValue("gas_porepressure")),
37  _gas_gradp_qp(coupledGradient("gas_porepressure")),
38  _gas_porepressure_varnum(coupled("gas_porepressure")),
39  _pvar(_dictator.isPorousFlowVariable(_gas_porepressure_varnum)
40  ? _dictator.porousFlowVariableNum(_gas_porepressure_varnum)
41  : 0),
42 
43  _num_z_vars(coupledComponents("z")),
44 
45  _aqueous_phase_number(getParam<unsigned int>("liquid_phase_number")),
46  _gas_phase_number(1 - _aqueous_phase_number),
47  _aqueous_fluid_component(getParam<unsigned int>("liquid_fluid_component")),
48  _gas_fluid_component(1 - _aqueous_fluid_component),
49 
50  _temperature(_nodal_material ? getMaterialProperty<Real>("PorousFlow_temperature_nodal")
51  : getMaterialProperty<Real>("PorousFlow_temperature_qp")),
52  _gradT_qp(getMaterialProperty<RealGradient>("PorousFlow_grad_temperature_qp")),
53  _dtemperature_dvar(
54  _nodal_material
55  ? getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_nodal_dvar")
56  : getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")),
57 
58  _mass_frac(_nodal_material
59  ? declareProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")
60  : declareProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_qp")),
61  _grad_mass_frac_qp(_nodal_material ? nullptr
62  : &declareProperty<std::vector<std::vector<RealGradient>>>(
63  "PorousFlow_grad_mass_frac_qp")),
64  _dmass_frac_dvar(_nodal_material ? declareProperty<std::vector<std::vector<std::vector<Real>>>>(
65  "dPorousFlow_mass_frac_nodal_dvar")
66  : declareProperty<std::vector<std::vector<std::vector<Real>>>>(
67  "dPorousFlow_mass_frac_qp_dvar")),
68  _saturation_old(_nodal_material
69  ? getMaterialPropertyOld<std::vector<Real>>("PorousFlow_saturation_nodal")
70  : getMaterialPropertyOld<std::vector<Real>>("PorousFlow_saturation_qp")),
71 
72  _fluid_density(_nodal_material
73  ? declareProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")
74  : declareProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
75  _dfluid_density_dvar(_nodal_material ? declareProperty<std::vector<std::vector<Real>>>(
76  "dPorousFlow_fluid_phase_density_nodal_dvar")
77  : declareProperty<std::vector<std::vector<Real>>>(
78  "dPorousFlow_fluid_phase_density_qp_dvar")),
79  _fluid_viscosity(_nodal_material
80  ? declareProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")
81  : declareProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
82  _dfluid_viscosity_dvar(
83  _nodal_material
84  ? declareProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")
85  : declareProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_qp_dvar")),
86 
87  _T_c2k(getParam<MooseEnum>("temperature_unit") == 0 ? 0.0 : 273.15),
88  _R(8.3144598),
89  _nr_max_its(42),
90  _nr_tol(1.0e-12),
91  _is_initqp(false),
92  _pc_uo(getUserObject<PorousFlowCapillaryPressure>("capillary_pressure"))
93 {
94  // Check that the number of total mass fractions provided as primary variables is correct
95  if (_num_z_vars != _num_components - 1)
96  mooseError("The number of supplied mass fraction variables should be ",
97  _num_components - 1,
98  " in ",
99  _name,
100  " but ",
101  _num_z_vars,
102  " are supplied");
103 
104  // Store all total mass fractions and associated variable numbers
105  _z.resize(_num_z_vars);
106  _gradz_qp.resize(_num_z_vars);
107  _z_varnum.resize(_num_z_vars);
108  _zvar.resize(_num_z_vars);
109 
110  for (unsigned int i = 0; i < _num_z_vars; ++i)
111  {
112  _z[i] = (_nodal_material ? &coupledNodalValue("z", i) : &coupledValue("z", i));
113  _gradz_qp[i] = &coupledGradient("z", i);
114  _z_varnum[i] = coupled("z", i);
115  _zvar[i] = (_dictator.isPorousFlowVariable(_z_varnum[i])
116  ? _dictator.porousFlowVariableNum(_z_varnum[i])
117  : 0);
118  }
119 
120  // Set the size of the FluidStateProperties vector
121  _fsp.resize(_num_phases);
122 
123  // Set the size of the mass fraction vectors for each phase
124  for (unsigned int ph = 0; ph < _num_phases; ++ph)
125  {
126  _fsp[ph].mass_fraction.resize(_num_components);
127  _fsp[ph].dmass_fraction_dp.resize(_num_components);
128  _fsp[ph].dmass_fraction_dT.resize(_num_components);
129  _fsp[ph].dmass_fraction_dz.resize(_num_components);
130  }
131 }
132 
133 void
135 {
136  _is_initqp = true;
137  // Set the size of pressure and saturation vectors
139 
140  // Set the size of all other vectors
142 
143  // Set the initial values of the properties at the nodes.
144  // Note: not required for qp materials as no old values at the qps are requested
145  if (_nodal_material)
146  {
148 
149  for (unsigned int ph = 0; ph < _num_phases; ++ph)
150  {
151  _saturation[_qp][ph] = _fsp[ph].saturation;
152  _porepressure[_qp][ph] = _fsp[ph].pressure;
153  _fluid_density[_qp][ph] = _fsp[ph].fluid_density;
154  _fluid_viscosity[_qp][ph] = _fsp[ph].fluid_viscosity;
155  _mass_frac[_qp][ph] = _fsp[ph].mass_fraction;
156  }
157  }
158 }
159 
160 void
162 {
163  _is_initqp = false;
164  // Prepare the derivative vectors
166 
167  // Set the size of all other vectors
169 
170  // Calculate all required thermophysical properties
172 
173  for (unsigned int ph = 0; ph < _num_phases; ++ph)
174  {
175  _saturation[_qp][ph] = _fsp[ph].saturation;
176  _porepressure[_qp][ph] = _fsp[ph].pressure;
177  _fluid_density[_qp][ph] = _fsp[ph].fluid_density;
178  _fluid_viscosity[_qp][ph] = _fsp[ph].fluid_viscosity;
179  _mass_frac[_qp][ph] = _fsp[ph].mass_fraction;
180  }
181 
182  // Derivative of saturation wrt variables
183  for (unsigned int ph = 0; ph < _num_phases; ++ph)
184  {
185  _dsaturation_dvar[_qp][ph][_zvar[0]] = _fsp[ph].dsaturation_dz;
186  _dsaturation_dvar[_qp][ph][_pvar] = _fsp[ph].dsaturation_dp;
187  }
188  // Derivative of capillary pressure
190 
191  // Derivative of porepressure wrt variables
192  if (_dictator.isPorousFlowVariable(_gas_porepressure_varnum))
193  {
194  for (unsigned int ph = 0; ph < _num_phases; ++ph)
195  {
196  _dporepressure_dvar[_qp][ph][_pvar] = 1.0;
197  if (!_nodal_material)
198  (*_dgradp_qp_dgradv)[_qp][ph][_pvar] = 1.0;
199  }
200 
201  if (!_nodal_material)
202  (*_dgradp_qp_dgradv)[_qp][_aqueous_phase_number][_zvar[0]] =
203  -dpc * _fsp[_aqueous_phase_number].dsaturation_dp -
204  dpc * _fsp[_aqueous_phase_number].dsaturation_dz;
205 
206  // The aqueous phase porepressure is also a function of liquid saturation,
207  // which depends on both gas porepressure and z
211  -dpc * _dsaturation_dvar[_qp][_aqueous_phase_number][_zvar[0]];
212  }
213 
214  // Calculate derivatives of material properties wrt primary variables
215  // Derivative of z wrt variables
216  std::vector<Real> dz_dvar;
217  dz_dvar.assign(_num_pf_vars, 0.0);
218  if (_dictator.isPorousFlowVariable(_z_varnum[0]))
219  dz_dvar[_zvar[0]] = 1.0;
220 
221  // Derivatives of properties wrt primary variables
222  for (unsigned int v = 0; v < _num_pf_vars; ++v)
223  {
224  for (unsigned int ph = 0; ph < _num_phases; ++ph)
225  {
226  // Derivative of density in each phase
227  _dfluid_density_dvar[_qp][ph][v] =
228  _fsp[ph].dfluid_density_dp * _dporepressure_dvar[_qp][ph][v];
229  _dfluid_density_dvar[_qp][ph][v] += _fsp[ph].dfluid_density_dT * _dtemperature_dvar[_qp][v];
230  _dfluid_density_dvar[_qp][ph][v] += _fsp[ph].dfluid_density_dz * dz_dvar[v];
231 
232  // Derivative of viscosity in each phase
233  _dfluid_viscosity_dvar[_qp][ph][v] =
234  _fsp[ph].dfluid_viscosity_dp * _dporepressure_dvar[_qp][ph][v];
235  _dfluid_viscosity_dvar[_qp][ph][v] +=
236  _fsp[ph].dfluid_viscosity_dT * _dtemperature_dvar[_qp][v];
237  _dfluid_viscosity_dvar[_qp][ph][v] += _fsp[ph].dfluid_viscosity_dz * dz_dvar[v];
238 
239  // The derivative of the mass fractions for each fluid component in each phase
240  for (unsigned int comp = 0; comp < _num_components; ++comp)
241  {
242  _dmass_frac_dvar[_qp][ph][comp][v] =
243  _fsp[ph].dmass_fraction_dp[comp] * _dporepressure_dvar[_qp][ph][v];
244  _dmass_frac_dvar[_qp][ph][comp][v] +=
245  _fsp[ph].dmass_fraction_dT[comp] * _dtemperature_dvar[_qp][v];
246  _dmass_frac_dvar[_qp][ph][comp][v] += _fsp[ph].dmass_fraction_dz[comp] * dz_dvar[v];
247  }
248  }
249  }
250 
251  // If the material properties are being evaluated at the qps, calculate the
252  // gradients as well. Note: only nodal properties are evaluated in
253  // initQpStatefulProperties(), so no need to check _is_initqp flag for qp
254  // properties
255  if (!_nodal_material)
256  {
257  // Second derivative of capillary pressure
258  Real d2pc = _pc_uo.d2CapillaryPressure(_fsp[_aqueous_phase_number].saturation);
259 
260  (*_grads_qp)[_qp][_gas_phase_number] =
262  _dsaturation_dvar[_qp][_gas_phase_number][_zvar[0]] * (*_gradz_qp[0])[_qp];
263  (*_grads_qp)[_qp][_aqueous_phase_number] = -(*_grads_qp)[_qp][_gas_phase_number];
264 
265  (*_gradp_qp)[_qp][_gas_phase_number] = _gas_gradp_qp[_qp];
266  (*_gradp_qp)[_qp][_aqueous_phase_number] =
267  _gas_gradp_qp[_qp] - dpc * (*_grads_qp)[_qp][_aqueous_phase_number];
268 
269  (*_dgradp_qp_dv)[_qp][_aqueous_phase_number][_zvar[0]] =
270  -d2pc * (*_grads_qp)[_qp][_aqueous_phase_number];
271 
272  (*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_aqueous_fluid_component] =
274  _gas_gradp_qp[_qp] +
276  (*_gradz_qp[0])[_qp];
277  (*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_gas_fluid_component] =
278  -(*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_aqueous_fluid_component];
279  (*_grad_mass_frac_qp)[_qp][_gas_phase_number][_aqueous_fluid_component] =
280  _fsp[_gas_phase_number].dmass_fraction_dp[_aqueous_fluid_component] * _gas_gradp_qp[_qp] +
281  _fsp[_gas_phase_number].dmass_fraction_dz[_aqueous_fluid_component] * (*_gradz_qp[0])[_qp];
282  (*_grad_mass_frac_qp)[_qp][_gas_phase_number][_gas_fluid_component] =
283  -(*_grad_mass_frac_qp)[_qp][_gas_phase_number][_aqueous_fluid_component];
284  }
285 }
286 
287 void
289 {
290  _fluid_density[_qp].assign(_num_phases, 0.0);
291  _fluid_viscosity[_qp].assign(_num_phases, 0.0);
292  _mass_frac[_qp].resize(_num_phases);
293 
294  // Derivatives and gradients are not required in initQpStatefulProperties
295  if (!_is_initqp)
296  {
297  _dfluid_density_dvar[_qp].resize(_num_phases);
298  _dfluid_viscosity_dvar[_qp].resize(_num_phases);
299  _dmass_frac_dvar[_qp].resize(_num_phases);
300 
301  if (!_nodal_material)
302  (*_grad_mass_frac_qp)[_qp].resize(_num_phases);
303 
304  for (unsigned int ph = 0; ph < _num_phases; ++ph)
305  {
306  _dfluid_density_dvar[_qp][ph].assign(_num_pf_vars, 0.0);
307  _dfluid_viscosity_dvar[_qp][ph].assign(_num_pf_vars, 0.0);
308  _dmass_frac_dvar[_qp][ph].resize(_num_components);
309 
310  for (unsigned int comp = 0; comp < _num_components; ++comp)
311  _dmass_frac_dvar[_qp][ph][comp].assign(_num_pf_vars, 0.0);
312 
313  if (!_nodal_material)
314  (*_grad_mass_frac_qp)[_qp][ph].assign(_num_components, RealGradient());
315  }
316  }
317 }
318 
319 Real
320 PorousFlowFluidStateFlashBase::rachfordRice(Real x, std::vector<Real> & Ki) const
321 {
322  // Check that the size of the equilibrium constant vector is correct
323  if (Ki.size() != _num_components)
324  mooseError("The number of equilibrium components passed to rachfordRice is not correct");
325 
326  Real f = 0.0;
327  Real z_total = 0.0;
328 
329  for (unsigned int i = 0; i < _num_z_vars; ++i)
330  {
331  f += (*_z[i])[_qp] * (Ki[i] - 1.0) / (1.0 + x * (Ki[i] - 1.0));
332  z_total += (*_z[i])[_qp];
333  }
334 
335  // Add the last component (with total mass fraction = 1 - z_total)
336  f += (1.0 - z_total) * (Ki[_num_z_vars] - 1.0) / (1.0 + x * (Ki[_num_z_vars] - 1.0));
337 
338  return f;
339 }
340 
341 Real
342 PorousFlowFluidStateFlashBase::rachfordRiceDeriv(Real x, std::vector<Real> & Ki) const
343 {
344  // Check that the size of the equilibrium constant vector is correct
345  if (Ki.size() != _num_components)
346  mooseError("The number of equilibrium components passed to rachfordRiceDeriv is not correct");
347 
348  Real df = 0.0;
349  Real z_total = 0.0;
350 
351  for (unsigned int i = 0; i < _num_z_vars; ++i)
352  {
353  df -= (*_z[i])[_qp] * (Ki[i] - 1.0) * (Ki[i] - 1.0) / (1.0 + x * (Ki[i] - 1.0)) /
354  (1.0 + x * (Ki[i] - 1.0));
355  z_total += (*_z[i])[_qp];
356  }
357 
358  // Add the last component (with total mass fraction = 1 - z_total)
359  df -= (1.0 - z_total) * (Ki[_num_z_vars] - 1.0) * (Ki[_num_z_vars] - 1.0) /
360  (1.0 + x * (Ki[_num_z_vars] - 1.0)) / (1.0 + x * (Ki[_num_z_vars] - 1.0));
361 
362  return df;
363 }
364 
365 Real
367 {
368  // Check that the size of the equilibrium constant vector is correct
369  if (Ki.size() != _num_components)
370  mooseError("The number of equilibrium components passed to vaporMassFraction is not correct");
371 
372  Real v;
373 
374  // If there are only two components, an analytical solution is possible
375  if (_num_components == 2)
376  v = ((*_z[0])[_qp] * (Ki[1] - Ki[0]) - (Ki[1] - 1.0)) / ((Ki[0] - 1.0) * (Ki[1] - 1.0));
377  else
378  {
379  // More than two components - solve the Rachford-Rice equation using
380  // Newton-Raphson method.
381  // Initial guess for vapor mass fraction
382  Real v0 = 0.5;
383  unsigned int iter = 0;
384 
385  while (std::abs(rachfordRice(v0, Ki)) > _nr_tol)
386  {
387  v0 = v0 - rachfordRice(v0, Ki) / rachfordRiceDeriv(v0, Ki);
388  iter++;
389 
390  if (iter > _nr_max_its)
391  break;
392  }
393  v = v0;
394  }
395  return v;
396 }
InputParameters validParams< PorousFlowFluidStateFlashBase >()
std::vector< const VariableGradient * > _gradz_qp
Gradient(s) of total mass fraction(s) of the gas component(s) (only defined at the qps) ...
void setMaterialVectorSize() const
Size material property vectors and initialise with zeros.
Real rachfordRiceDeriv(Real x, std::vector< Real > &Ki) const
Derivative of Rachford-Rice equation wrt vapor fraction.
MaterialProperty< std::vector< std::vector< Real > > > & _dporepressure_dvar
d(porepressure)/d(PorousFlow variable)
PorousFlowFluidStateFlashBase(const InputParameters &parameters)
const MaterialProperty< std::vector< Real > > & _dtemperature_dvar
Derivative of temperature wrt PorousFlow variables.
const VariableGradient & _gas_gradp_qp
Gradient of porepressure (only defined at the qps)
MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_viscosity_dvar
Derivative of the fluid viscosity for each phase wrt PorousFlow variables.
const Real _nr_tol
Tolerance for Newton-Raphson iterations.
const unsigned int _num_pf_vars
Number of PorousFlow variables.
const unsigned int _gas_porepressure_varnum
Moose variable number of the gas porepressure.
std::vector< unsigned int > _z_varnum
Moose variable number of z.
const unsigned int _pvar
PorousFlow variable number of the gas porepressure.
virtual void initQpStatefulProperties() override
MaterialProperty< std::vector< Real > > & _saturation
Computed nodal or qp saturation of the phases.
MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_density_dvar
Derivative of the fluid density for each phase wrt PorousFlow variables.
Base class for capillary pressure for multiphase flow in porous media.
MaterialProperty< std::vector< std::vector< Real > > > & _dsaturation_dvar
d(saturation)/d(PorousFlow variable)
virtual Real dCapillaryPressure(Real saturation) const
Derivative of capillary pressure wrt true saturation.
MaterialProperty< std::vector< Real > > & _fluid_density
Fluid density of each phase.
virtual Real d2CapillaryPressure(Real saturation) const
Second derivative of capillary pressure wrt true saturation.
InputParameters validParams< PorousFlowVariableBase >()
MaterialProperty< std::vector< std::vector< Real > > > & _mass_frac
Mass fraction matrix.
const unsigned int _gas_fluid_component
Fluid component number of the gas phase.
const PorousFlowCapillaryPressure & _pc_uo
Capillary pressure UserObject.
virtual void initQpStatefulProperties() override
std::vector< FluidStateProperties > _fsp
FluidStateProperties data structure.
std::vector< unsigned int > _zvar
PorousFlow variable number of z.
const unsigned int _aqueous_phase_number
Phase number of the aqueous phase.
Real rachfordRice(Real x, std::vector< Real > &Ki) const
Rachford-Rice equation for vapor fraction.
virtual void computeQpProperties() override
const unsigned int _num_z_vars
Number of coupled total mass fractions. Should be _num_phases - 1.
const unsigned int _gas_phase_number
Phase number of the gas phase.
bool _is_initqp
Flag to indicate whether to calculate stateful properties.
const unsigned int _aqueous_fluid_component
Fluid component number of the aqueous component.
const unsigned int _num_components
Number of components.
const unsigned int _num_phases
Number of phases.
const Real _nr_max_its
Maximum number of iterations for the Newton-Raphson iterations.
void FORTRAN_CALL() saturation(double &P, double &T, int &N, int &nerr)
Base class for thermophysical variable materials, which assemble materials for primary variables such...
MaterialProperty< std::vector< Real > > & _porepressure
Computed nodal or quadpoint values of porepressure of the phases.
MaterialProperty< std::vector< Real > > & _fluid_viscosity
Viscosity of each phase.
std::vector< const VariableValue * > _z
Total mass fraction(s) of the gas component(s) summed over all phases.
MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _dmass_frac_dvar
Derivative of the mass fraction matrix with respect to the Porous Flow variables. ...
Real vaporMassFraction(std::vector< Real > &Ki) const
Solves Rachford-Rice equation to provide vapor mass fraction.
virtual void thermophysicalProperties()=0
Calculates all required thermophysical properties and derivatives for each phase and fluid component...