www.mooseframework.org
PorousFlowMassFractionAqueousEquilibriumChemistry.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 
11 
13 
14 template <>
15 InputParameters
17 {
18  InputParameters params = validParams<PorousFlowMaterialVectorBase>();
19  params.addRequiredCoupledVar(
20  "mass_fraction_vars",
21  "List of variables that represent the mass fractions. For the aqueous phase these are "
22  "concentrations of the primary species with units m^{3}(chemical)/m^{3}(fluid phase). For "
23  "the other phases (if any) these will typically be initialised to zero and will not change "
24  "throughout the simulation. Format is 'f_ph0^c0 f_ph0^c1 f_ph0^c2 ... f_ph0^c(N-1) f_ph1^c0 "
25  "f_ph1^c1 fph1^c2 ... fph1^c(N-1) ... fphP^c0 f_phP^c1 fphP^c2 ... fphP^c(N-1)' where "
26  "N=number of primary species and P=num_phases, and it is assumed that "
27  "f_ph^cN=1-sum(f_ph^c,{c,0,N-1}) so that f_ph^cN need not be given.");
28  params.addRequiredParam<unsigned>("num_reactions",
29  "Number of equations in the system of chemical reactions");
30  params.addParam<bool>("equilibrium_constants_as_log10",
31  false,
32  "If true, the equilibrium constants are written in their log10 form, eg, "
33  "-2. If false, the equilibrium constants are written in absolute terms, "
34  "eg, 0.01");
35  params.addRequiredCoupledVar("equilibrium_constants",
36  "Equilibrium constant for each equation (dimensionless). If these "
37  "are temperature dependent AuxVariables, the Jacobian will not be "
38  "exact");
39  params.addRequiredParam<std::vector<Real>>(
40  "primary_activity_coefficients",
41  "Activity coefficients for the primary species (dimensionless) (one for each)");
42  params.addRequiredParam<std::vector<Real>>(
43  "reactions",
44  "A matrix defining the aqueous reactions. The matrix is entered as a long vector: the first "
45  "row is "
46  "entered first, followed by the second row, etc. There should be num_reactions rows. All "
47  "primary species should appear only on the LHS of each reaction (and there should be just "
48  "one secondary species on the RHS, by definition) so they may have negative coefficients. "
49  "Each row should have number of primary_concentrations entries, which are the stoichiometric "
50  "coefficients. The first coefficient must always correspond to the first primary species, "
51  "etc");
52  params.addRequiredParam<std::vector<Real>>(
53  "secondary_activity_coefficients",
54  "Activity coefficients for the secondary species (dimensionless) (one for each reaction)");
55  params.addPrivateParam<std::string>("pf_material_type", "mass_fraction");
56  params.addClassDescription(
57  "This Material forms a std::vector<std::vector ...> of mass-fractions "
58  "(total concentrations of primary species (m^{3}(primary species)/m^{3}(solution)) and since "
59  "this is for an aqueous system only, mass-fraction equals volume-fraction) corresponding to "
60  "an "
61  "aqueous equilibrium chemistry system. The first mass fraction is the "
62  "concentration of the first primary species, etc, and the last mass "
63  "fraction is the concentration of H2O.");
64  return params;
65 }
66 
68  PorousFlowMassFractionAqueousEquilibriumChemistry(const InputParameters & parameters)
69  : PorousFlowMassFraction(parameters),
70  _sec_conc(_nodal_material
71  ? declareProperty<std::vector<Real>>("PorousFlow_secondary_concentration_nodal")
72  : declareProperty<std::vector<Real>>("PorousFlow_secondary_concentration_qp")),
73  _dsec_conc_dvar(_nodal_material ? declareProperty<std::vector<std::vector<Real>>>(
74  "dPorousFlow_secondary_concentration_nodal_dvar")
75  : declareProperty<std::vector<std::vector<Real>>>(
76  "dPorousFlow_secondary_concentration_qp_dvar")),
77 
78  _temperature(_nodal_material ? getMaterialProperty<Real>("PorousFlow_temperature_nodal")
79  : getMaterialProperty<Real>("PorousFlow_temperature_qp")),
80  _dtemperature_dvar(
81  _nodal_material
82  ? getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_nodal_dvar")
83  : getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")),
84 
85  _num_primary(_num_components - 1),
86  _aq_ph(_dictator.aqueousPhaseNumber()),
87  _aq_i(_aq_ph * _num_primary),
88  _num_reactions(getParam<unsigned>("num_reactions")),
89  _equilibrium_constants_as_log10(getParam<bool>("equilibrium_constants_as_log10")),
90  _num_equilibrium_constants(coupledComponents("equilibrium_constants")),
91  _equilibrium_constants(_num_equilibrium_constants),
92  _primary_activity_coefficients(getParam<std::vector<Real>>("primary_activity_coefficients")),
93  _reactions(getParam<std::vector<Real>>("reactions")),
94  _secondary_activity_coefficients(getParam<std::vector<Real>>("secondary_activity_coefficients"))
95 {
96  if (_dictator.numPhases() < 1)
97  mooseError("PorousFlowMassFractionAqueousEquilibriumChemistry: The number of fluid phases must "
98  "not be zero");
99 
100  // correct number of equilibrium constants
102  mooseError("PorousFlowMassFractionAqueousEquilibriumChemistry: The number of "
103  "equilibrium constants is ",
105  " which must be equal to the number of reactions (",
107  ")");
108 
109  // correct number of activity coefficients
111  mooseError("PorousFlowMassFractionAqueousEquilibriumChemistry: The number of primary activity "
112  "coefficients is ",
114  " which must be equal to the number of primary species (",
115  _num_primary,
116  ")");
117 
118  // correct number of stoichiometry coefficients
119  if (_reactions.size() != _num_reactions * _num_primary)
120  mooseError("PorousFlowMassFractionAqueousEquilibriumChemistry: The number of stoichiometric "
121  "coefficients specified in 'reactions' (",
122  _reactions.size(),
123  ") must be equal to the number of reactions (",
125  ") multiplied by the number of primary species (",
126  _num_primary,
127  ")");
128 
129  // correct number of secondary activity coefficients
131  mooseError(
132  "PorousFlowMassFractionAqueousEquilibriumChemistry: The number of secondary activity "
133  "coefficients is ",
135  " which must be equal to the number of secondary species (",
137  ")");
138 
139  // correct number of reactions
140  if (_num_reactions != _dictator.numAqueousEquilibrium())
141  mooseError("PorousFlowMassFractionAqueousEquilibriumChemistry: You have specified the number "
142  "of reactions to be ",
144  " but the Dictator knows that the number of aqueous equilibrium reactions is ",
145  _dictator.numAqueousEquilibrium());
146 
147  for (unsigned i = 0; i < _num_equilibrium_constants; ++i)
148  _equilibrium_constants[i] = (_nodal_material ? &coupledDofValues("equilibrium_constants", i)
149  : &coupledValue("equilibrium_constants", i));
150 }
151 
152 void
154 {
156 }
157 
158 void
160 {
161  // size all properties correctly and populate the non-aqueous phase info
163 
164  // size the secondary concentrations
165  _sec_conc[_qp].resize(_num_reactions);
166  _dsec_conc_dvar[_qp].resize(_num_reactions);
167  for (unsigned r = 0; r < _num_reactions; ++r)
168  _dsec_conc_dvar[_qp][r].assign(_num_var, 0.0);
169 
170  // Compute the secondary concentrations
171  if (_t_step == 0 && !_app.isRestarting())
173  else
175 
176  // compute _mass_frac[_qp][_aq_ph]
177  _mass_frac[_qp][_aq_ph][_num_components - 1] = 1.0; // the final component is H20
178  for (unsigned i = 0; i < _num_primary; ++i)
179  {
180  _mass_frac[_qp][_aq_ph][i] = (*_mf_vars[_aq_i + i])[_qp];
181  for (unsigned r = 0; r < _num_reactions; ++r)
182  _mass_frac[_qp][_aq_ph][i] += stoichiometry(r, i) * _sec_conc[_qp][r];
183 
184  // remove mass-fraction from the H20 component
185  _mass_frac[_qp][_aq_ph][_num_components - 1] -= _mass_frac[_qp][_aq_ph][i];
186  }
187 
188  // Compute the derivatives of the secondary concentrations
189  std::vector<std::vector<Real>> dsec(_num_reactions);
190  std::vector<Real> dsec_dT(_num_reactions);
191  for (unsigned r = 0; r < _num_reactions; ++r)
192  {
194  dsec_dT[r] = dQpSecondaryConcentration_dT(r);
195  }
196 
197  // Compute the derivatives of the mass_frac wrt the primary concentrations
198  // This is used in _dmass_frac_dvar as well as _grad_mass_frac
199  std::vector<std::vector<Real>> dmf(_num_components);
200  for (unsigned i = 0; i < _num_components; ++i)
201  dmf[i].assign(_num_primary, 0.0);
202  for (unsigned wrt = 0; wrt < _num_primary; ++wrt)
203  {
204  // run through the mass fractions (except the last one) adding to their derivatives
205  // The special case is:
206  dmf[wrt][wrt] = 1.0;
207  // The secondary-species contributions are:
208  for (unsigned i = 0; i < _num_primary; ++i)
209  for (unsigned r = 0; r < _num_reactions; ++r)
210  dmf[i][wrt] += stoichiometry(r, i) * dsec[r][wrt];
211 
212  // compute dmf[_num_components - 1]
213  for (unsigned i = 0; i < _num_primary; ++i)
214  dmf[_num_components - 1][wrt] -= dmf[i][wrt];
215  }
216 
217  // Compute the derivatives of the mass_frac wrt the temperature
218  // This is used in _dmass_frac_dvar
219  std::vector<Real> dmf_dT(_num_components, 0.0);
220  for (unsigned i = 0; i < _num_components - 1; ++i)
221  {
222  for (unsigned r = 0; r < _num_reactions; ++r)
223  dmf_dT[i] += stoichiometry(r, i) * dsec_dT[r];
224  dmf_dT[_num_components - 1] -= dmf_dT[i];
225  }
226 
227  // compute _dmass_frac_dvar[_qp][_aq_ph] and _dsec_conc_dvar[_qp]
228  for (unsigned wrt = 0; wrt < _num_primary; ++wrt)
229  {
230  // derivative with respect to the "wrt"^th primary species concentration
231  if (!_dictator.isPorousFlowVariable(_mf_vars_num[_aq_i + wrt]))
232  continue;
233  const unsigned pf_wrt = _dictator.porousFlowVariableNum(_mf_vars_num[_aq_i + wrt]);
234 
235  // run through the mass fractions, building the derivative using dmf
236  for (unsigned i = 0; i < _num_components; ++i)
237  _dmass_frac_dvar[_qp][_aq_ph][i][pf_wrt] = dmf[i][wrt];
238 
239  // run through the secondary concentrations, using dsec in the appropriate places
240  for (unsigned r = 0; r < _num_reactions; ++r)
241  _dsec_conc_dvar[_qp][r][pf_wrt] = dsec[r][wrt];
242  }
243 
244  // use the derivative wrt temperature
245  for (unsigned i = 0; i < _num_components; ++i)
246  for (unsigned v = 0; v < _num_var; ++v)
247  _dmass_frac_dvar[_qp][_aq_ph][i][v] += dmf_dT[i] * _dtemperature_dvar[_qp][v];
248  for (unsigned r = 0; r < _num_reactions; ++r)
249  for (unsigned v = 0; v < _num_var; ++v)
250  _dsec_conc_dvar[_qp][r][v] += dsec_dT[r] * _dtemperature_dvar[_qp][v];
251 
252  // compute the gradient, if needed
253  // NOTE: The derivative d(grad_mass_frac)/d(var) != d(mass_frac)/d(var) * grad_phi
254  // because mass fraction is a nonlinear function of the primary variables
255  // This means that the Jacobian in PorousFlowDispersiveFlux will be wrong
256  if (!_nodal_material)
257  {
258  (*_grad_mass_frac)[_qp][_aq_ph][_num_components - 1] = 0.0;
259  for (unsigned comp = 0; comp < _num_components - 1; ++comp)
260  {
261  (*_grad_mass_frac)[_qp][_aq_ph][comp] = 0.0;
262  for (unsigned wrt = 0; wrt < _num_primary; ++wrt)
263  (*_grad_mass_frac)[_qp][_aq_ph][comp] +=
264  dmf[comp][wrt] * (*_grad_mf_vars[_aq_i + wrt])[_qp];
265  (*_grad_mass_frac)[_qp][_aq_ph][_num_components - 1] -= (*_grad_mass_frac)[_qp][_aq_ph][comp];
266  }
267  }
268 }
269 
270 Real
272  unsigned primary_num) const
273 {
274  const unsigned index = reaction_num * _num_primary + primary_num;
275  return _reactions[index];
276 }
277 
278 void
280  unsigned & zero_conc_index, unsigned & zero_count) const
281 {
282  zero_count = 0;
283  for (unsigned i = 0; i < _num_primary; ++i)
284  {
285  if (_primary_activity_coefficients[i] * (*_mf_vars[_aq_i + i])[_qp] <= 0.0)
286  {
287  zero_count += 1;
288  zero_conc_index = i;
289  if (zero_count > 1)
290  return;
291  }
292  }
293  return;
294 }
295 
296 void
298 {
299  _sec_conc[_qp].assign(_num_reactions, 0.0);
300 }
301 
302 void
304 {
305  for (unsigned r = 0; r < _num_reactions; ++r)
306  {
307  _sec_conc[_qp][r] = 1.0;
308  for (unsigned i = 0; i < _num_primary; ++i)
309  {
310  const Real gamp = _primary_activity_coefficients[i] * (*_mf_vars[_aq_i + i])[_qp];
311  if (gamp <= 0.0)
312  {
313  if (stoichiometry(r, i) < 0.0)
314  _sec_conc[_qp][r] = std::numeric_limits<Real>::max();
315  else if (stoichiometry(r, i) == 0.0)
316  _sec_conc[_qp][r] *= 1.0;
317  else
318  {
319  _sec_conc[_qp][r] = 0.0;
320  break;
321  }
322  }
323  else
324  _sec_conc[_qp][r] *= std::pow(gamp, stoichiometry(r, i));
325  }
326  _sec_conc[_qp][r] *=
328  : (*_equilibrium_constants[r])[_qp]);
330  }
331 }
332 
333 void
335  unsigned reaction_num, std::vector<Real> & dsc) const
336 {
337  dsc.assign(_num_primary, 0.0);
338 
339  /*
340  * the derivatives are straightforward if all primary > 0.
341  *
342  * If more than one primary = 0 then I set the derivatives to zero, even though it could be argued
343  * that with certain stoichiometric coefficients you might have derivative = 0/0 and it might be
344  * appropriate to set this to a non-zero finite value.
345  *
346  * If exactly one primary = 0 and its stoichiometry = 1 then the derivative wrt this one is
347  * nonzero.
348  * If exactly one primary = 0 and its stoichiometry > 1 then all derivatives are zero.
349  * If exactly one primary = 0 and its stoichiometry < 1 then the derivative wrt this one is
350  * infinity
351  */
352 
353  unsigned zero_count = 0;
354  unsigned zero_conc_index = 0;
355  findZeroConcentration(zero_conc_index, zero_count);
356 
357  if (zero_count == 0)
358  {
359  for (unsigned i = 0; i < _num_primary; ++i)
360  dsc[i] = stoichiometry(reaction_num, i) * _sec_conc[_qp][reaction_num] /
361  (*_mf_vars[_aq_i + i])[_qp];
362  }
363  else
364  {
365  // count the number of primary <= 0, and record the one that's zero
366  if (zero_count == 1 and stoichiometry(reaction_num, zero_conc_index) == 1.0)
367  {
368  Real conc_without_zero = 1.0;
369  for (unsigned i = 0; i < _num_primary; ++i)
370  {
371  if (i == zero_conc_index)
372  conc_without_zero *= _primary_activity_coefficients[i];
373  else
374  conc_without_zero *=
376  stoichiometry(reaction_num, i));
377  }
378  conc_without_zero *= (_equilibrium_constants_as_log10
379  ? std::pow(10.0, (*_equilibrium_constants[reaction_num])[_qp])
380  : (*_equilibrium_constants[reaction_num])[_qp]);
381  conc_without_zero /= _secondary_activity_coefficients[reaction_num];
382  dsc[zero_conc_index] = conc_without_zero;
383  }
384  else if (zero_count == 1 && stoichiometry(reaction_num, zero_conc_index) < 1.0)
385  dsc[zero_conc_index] = std::numeric_limits<Real>::max();
386 
387  // all other cases have dsc = 0
388  }
389 }
390 
391 Real
393  unsigned /* reaction_num */) const
394 {
395  return 0.0;
396 }
const bool _equilibrium_constants_as_log10
Whether the equilibium constants are written in their log10 form, or in absolute terms.
InputParameters validParams< PorousFlowMassFractionAqueousEquilibriumChemistry >()
const unsigned int _num_reactions
Number of equations in the aqueous geochemistry system.
virtual void computeQpProperties() override
Material designed to form a std::vector<std::vector> of mass fractions from the individual mass fract...
MaterialProperty< std::vector< std::vector< RealGradient > > > *const _grad_mass_frac
Gradient of the mass fraction matrix at the quad points.
MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _dmass_frac_dvar
Derivative of the mass fraction matrix with respect to the porous flow variables. ...
const std::vector< Real > _reactions
Stoichiometry defining the aqeuous geochemistry equilibrium reactions.
MaterialProperty< std::vector< Real > > & _sec_conc
Secondary concentrations at quadpoint or nodes.
virtual void computeQpSecondaryConcentrations()
Compute the secondary-species concentration as defined by the chemistry Must be overridden by derived...
const std::vector< Real > _primary_activity_coefficients
Activity coefficients for the primary species (dimensionless)
virtual Real dQpSecondaryConcentration_dT(unsigned reaction_num) const
Computes derivative of the secondary concentration with respect to the temperature Must be overridden...
void findZeroConcentration(unsigned &zero_conc_index, unsigned &zero_count) const
Checks gamp[i] = _primary_activity_coefficients[i] * (*_primary[i])[qp].
Real stoichiometry(unsigned reaction_num, unsigned primary_num) const
The stoichiometric coefficient.
InputParameters validParams< PorousFlowMaterialVectorBase >()
std::vector< unsigned int > _mf_vars_num
The variable number of the mass-fraction variables.
const MaterialProperty< std::vector< Real > > & _dtemperature_dvar
d(temperature)/(d porflow variable)
const unsigned int _num_components
Number of fluid components.
virtual void dQpSecondaryConcentration_dprimary(unsigned reaction_num, std::vector< Real > &dsc) const
Computes derivative of the secondary concentration with respect to the primary concentrations Must be...
const unsigned int _num_var
Number of PorousFlow variables.
const unsigned int _aq_i
Index (into _mf_vars) of the first of the primary species.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
std::vector< const VariableValue * > _mf_vars
The mass-fraction variables.
virtual void initQpSecondaryConcentrations()
Initialises (at _t_step = 0) the secondary concentrations.
const unsigned _num_equilibrium_constants
Number of equilibrium_constants provided.
std::vector< const VariableGradient * > _grad_mf_vars
The gradient of the mass-fraction variables.
Material designed to form a std::vector<std::vector> of mass fractions from primary-species concentra...
std::vector< const VariableValue * > _equilibrium_constants
Equilibrium constants (dimensionless)
MaterialProperty< std::vector< std::vector< Real > > > & _mass_frac
Mass fraction matrix at quadpoint or nodes.
registerMooseObject("PorousFlowApp", PorousFlowMassFractionAqueousEquilibriumChemistry)
MaterialProperty< std::vector< std::vector< Real > > > & _dsec_conc_dvar
Derivative of the secondary concentrations with respect to the porous flow variables.
const std::vector< Real > _secondary_activity_coefficients
Activity coefficients for the secondary species.