www.mooseframework.org
PorousFlowAqueousPreDisChemistry.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  "primary_concentrations",
21  "List of MOOSE Variables that represent the concentrations of the primary species");
22  params.addRequiredParam<unsigned>("num_reactions",
23  "Number of equations in the system of chemical reactions");
24  params.addParam<bool>("equilibrium_constants_as_log10",
25  false,
26  "If true, the equilibrium constants are written in their log10 form, eg, "
27  "-2. If false, the equilibrium constants are written in absolute terms, "
28  "eg, 0.01");
29  params.addRequiredCoupledVar("equilibrium_constants",
30  "Equilibrium constant for each equation (dimensionless). If these "
31  "are temperature dependent AuxVariables, the Jacobian will not be "
32  "exact");
33  params.addRequiredParam<std::vector<Real>>(
34  "primary_activity_coefficients",
35  "Activity coefficients for the primary species (dimensionless) (one for each)");
36  params.addRequiredParam<std::vector<Real>>(
37  "reactions",
38  "A matrix defining the aqueous reactions. The matrix is entered as a long vector: the first "
39  "row is "
40  "entered first, followed by the second row, etc. There should be num_reactions rows. All "
41  "primary species should appear only on the LHS of each reaction (and there should be just "
42  "one secondary species on the RHS, by definition) so they may have negative coefficients. "
43  "Each row should have number of primary_concentrations entries, which are the stoichiometric "
44  "coefficients. The first coefficient must always correspond to the first primary species, "
45  "etc");
46  params.addRequiredParam<std::vector<Real>>("specific_reactive_surface_area",
47  "Specific reactive surface area in m^2/(L solution).");
48  params.addRequiredParam<std::vector<Real>>(
49  "kinetic_rate_constant",
50  "Kinetic rate constant in mol/(m^2 s), at the reference temperature (one for each reaction)");
51  params.addRequiredParam<std::vector<Real>>("molar_volume",
52  "Volume occupied by one mole of the secondary species "
53  "(L(solution)/mol) (one for each reaction)");
54  params.addRequiredParam<std::vector<Real>>("activation_energy",
55  "Activation energy, J/mol (one for each reaction)");
56  params.addParam<Real>("gas_constant", 8.31434, "Gas constant, in J/(mol K)");
57  params.addParam<Real>("reference_temperature", 298.15, "Reference temperature, K");
58  params.addParam<std::vector<Real>>("theta_exponent",
59  "Theta exponent. Defaults to 1. (one for each reaction)");
60  params.addParam<std::vector<Real>>("eta_exponent",
61  "Eta exponent. Defaults to 1. (one for each reaction)");
62  params.addPrivateParam<std::string>("pf_material_type", "chemistry");
63  params.addClassDescription("This Material forms a std::vector of mineralisation reaction rates "
64  "(L(precipitate)/L(solution)/s) appropriate to the aqueous "
65  "precipitation-dissolution system provided. Note: the "
66  "PorousFlowTemperature must be measured in Kelvin.");
67  return params;
68 }
69 
71  const InputParameters & parameters)
72  : PorousFlowMaterialVectorBase(parameters),
73  _porosity_old(_nodal_material ? getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")
74  : getMaterialPropertyOld<Real>("PorousFlow_porosity_qp")),
75  _aq_ph(_dictator.aqueousPhaseNumber()),
76  _saturation(_nodal_material
77  ? getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")
78  : getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_qp")),
79 
80  _temperature(_nodal_material ? getMaterialProperty<Real>("PorousFlow_temperature_nodal")
81  : getMaterialProperty<Real>("PorousFlow_temperature_qp")),
82  _dtemperature_dvar(
83  _nodal_material
84  ? getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_nodal_dvar")
85  : getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")),
86 
87  _num_primary(coupledComponents("primary_concentrations")),
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 
95  _sec_conc_old(
96  _nodal_material
97  ? getMaterialPropertyOld<std::vector<Real>>("PorousFlow_mineral_concentration_nodal")
98  : getMaterialPropertyOld<std::vector<Real>>("PorousFlow_mineral_concentration_qp")),
99 
100  _mineral_sat(_num_reactions),
101  _bounded_rate(_num_reactions),
102  _reaction_rate(
103  _nodal_material
104  ? declareProperty<std::vector<Real>>("PorousFlow_mineral_reaction_rate_nodal")
105  : declareProperty<std::vector<Real>>("PorousFlow_mineral_reaction_rate_qp")),
106  _dreaction_rate_dvar(_nodal_material ? declareProperty<std::vector<std::vector<Real>>>(
107  "dPorousFlow_mineral_reaction_rate_nodal_dvar")
108  : declareProperty<std::vector<std::vector<Real>>>(
109  "dPorousFlow_mineral_reaction_rate_qp_dvar")),
110 
111  _r_area(getParam<std::vector<Real>>("specific_reactive_surface_area")),
112  _molar_volume(getParam<std::vector<Real>>("molar_volume")),
113  _ref_kconst(getParam<std::vector<Real>>("kinetic_rate_constant")),
114  _e_act(getParam<std::vector<Real>>("activation_energy")),
115  _gas_const(getParam<Real>("gas_constant")),
116  _one_over_ref_temp(1.0 / getParam<Real>("reference_temperature")),
117  _theta_exponent(isParamValid("theta_exponent") ? getParam<std::vector<Real>>("theta_exponent")
118  : std::vector<Real>(_num_reactions, 1.0)),
119  _eta_exponent(isParamValid("eta_exponent") ? getParam<std::vector<Real>>("eta_exponent")
120  : std::vector<Real>(_num_reactions, 1.0))
121 {
122  if (_dictator.numPhases() < 1)
123  mooseError("PorousFlowAqueousPreDisChemistry: The number of fluid phases must not be zero");
124 
125  if (_num_primary != _num_components - 1)
126  mooseError("PorousFlowAqueousPreDisChemistry: The number of mass_fraction_vars is ",
128  " which must be one greater than the number of primary concentrations (which is ",
129  _num_primary,
130  ")");
131 
132  // correct number of equilibrium constants
134  mooseError("PorousFlowAqueousPreDisChemistry: The number of equilibrium constants is ",
136  " which must be equal to the number of reactions (",
138  ")");
139 
140  // correct number of activity coefficients
142  mooseError("PorousFlowAqueousPreDisChemistry: The number of primary activity "
143  "coefficients is ",
145  " which must be equal to the number of primary species (",
146  _num_primary,
147  ")");
148 
149  // correct number of stoichiometry coefficients
150  if (_reactions.size() != _num_reactions * _num_primary)
151  mooseError("PorousFlowAqueousPreDisChemistry: The number of stoichiometric "
152  "coefficients specified in 'reactions' (",
153  _reactions.size(),
154  ") must be equal to the number of reactions (",
156  ") multiplied by the number of primary species (",
157  _num_primary,
158  ")");
159 
160  if (_r_area.size() != _num_reactions)
161  mooseError("PorousFlowAqueousPreDisChemistry: The number of specific reactive "
162  "surface areas provided is ",
163  _r_area.size(),
164  " which must be equal to the number of reactions (",
166  ")");
167 
168  if (_ref_kconst.size() != _num_reactions)
169  mooseError("PorousFlowAqueousPreDisChemistry: The number of kinetic rate constants is ",
170  _ref_kconst.size(),
171  " which must be equal to the number of reactions (",
173  ")");
174 
175  if (_e_act.size() != _num_reactions)
176  mooseError("PorousFlowAqueousPreDisChemistry: The number of activation energies is ",
177  _e_act.size(),
178  " which must be equal to the number of reactions (",
180  ")");
181 
182  if (_molar_volume.size() != _num_reactions)
183  mooseError("PorousFlowAqueousPreDisChemistry: The number of molar volumes is ",
184  _molar_volume.size(),
185  " which must be equal to the number of reactions (",
187  ")");
188 
189  if (_theta_exponent.size() != _num_reactions)
190  mooseError("PorousFlowAqueousPreDisChemistry: The number of theta exponents is ",
191  _theta_exponent.size(),
192  " which must be equal to the number of reactions (",
194  ")");
195 
196  if (_eta_exponent.size() != _num_reactions)
197  mooseError("PorousFlowAqueousPreDisChemistry: The number of eta exponents is ",
198  _eta_exponent.size(),
199  " which must be equal to the number of reactions (",
201  ")");
202 
203  if (_num_reactions != _dictator.numAqueousKinetic())
204  mooseError("PorousFlowAqueousPreDisChemistry: You have specified the number of "
205  "reactions to be ",
207  " but the Dictator knows that the number of aqueous kinetic "
208  "(precipitation-dissolution) reactions is ",
209  _dictator.numAqueousKinetic());
210 
212  _primary.resize(_num_primary);
213  for (unsigned i = 0; i < _num_primary; ++i)
214  {
215  _primary_var_num[i] = coupled("primary_concentrations", i);
216  _primary[i] = (_nodal_material ? &coupledDofValues("primary_concentrations", i)
217  : &coupledValue("primary_concentrations", i));
218  }
219 
220  for (unsigned i = 0; i < _num_equilibrium_constants; ++i)
221  _equilibrium_constants[i] = (_nodal_material ? &coupledDofValues("equilibrium_constants", i)
222  : &coupledValue("equilibrium_constants", i));
223 }
224 
225 void
227 {
228  _reaction_rate[_qp].resize(_num_reactions);
230  for (unsigned r = 0; r < _num_reactions; ++r)
231  _dreaction_rate_dvar[_qp][r].assign(_num_var, 0.0);
232 }
233 
234 void
236 {
237  _reaction_rate[_qp].resize(_num_reactions);
239  for (unsigned r = 0; r < _num_reactions; ++r)
240  _dreaction_rate_dvar[_qp][r].assign(_num_var, 0.0);
241 
242  // Compute the reaction rates
244 
245  // Compute the derivatives of the reaction rates
246  std::vector<std::vector<Real>> drr(_num_reactions);
247  std::vector<Real> drr_dT(_num_reactions);
248  for (unsigned r = 0; r < _num_reactions; ++r)
249  {
250  dQpReactionRate_dprimary(r, drr[r]);
251  drr_dT[r] = dQpReactionRate_dT(r);
252  }
253 
254  // compute _dreaction_rate_dvar[_qp]
255  for (unsigned wrt = 0; wrt < _num_primary; ++wrt)
256  {
257  // derivative with respect to the "wrt"^th primary species concentration
258  if (!_dictator.isPorousFlowVariable(_primary_var_num[wrt]))
259  continue;
260  const unsigned pf_wrt = _dictator.porousFlowVariableNum(_primary_var_num[wrt]);
261 
262  // run through the reactions, using drr in the appropriate places
263  for (unsigned r = 0; r < _num_reactions; ++r)
264  _dreaction_rate_dvar[_qp][r][pf_wrt] = drr[r][wrt];
265  }
266 
267  // use the derivative wrt temperature
268  for (unsigned r = 0; r < _num_reactions; ++r)
269  for (unsigned v = 0; v < _num_var; ++v)
270  _dreaction_rate_dvar[_qp][r][v] += drr_dT[r] * _dtemperature_dvar[_qp][v];
271 }
272 
273 Real
274 PorousFlowAqueousPreDisChemistry::stoichiometry(unsigned reaction_num, unsigned primary_num) const
275 {
276  const unsigned index = reaction_num * _num_primary + primary_num;
277  return _reactions[index];
278 }
279 
280 void
282  unsigned & zero_count) const
283 {
284  zero_count = 0;
285  for (unsigned i = 0; i < _num_primary; ++i)
286  {
287  if (_primary_activity_coefficients[i] * (*_primary[i])[_qp] <= 0.0)
288  {
289  zero_count += 1;
290  zero_conc_index = i;
291  if (zero_count > 1)
292  return;
293  }
294  }
295  return;
296 }
297 
298 void
300 {
301  for (unsigned r = 0; r < _num_reactions; ++r)
302  {
303  _mineral_sat[r] =
305  : 1.0 / (*_equilibrium_constants[r])[_qp]);
306  for (unsigned j = 0; j < _num_primary; ++j)
307  {
308  const Real gamp = _primary_activity_coefficients[j] * (*_primary[j])[_qp];
309  if (gamp <= 0.0)
310  {
311  if (stoichiometry(r, j) < 0.0)
312  _mineral_sat[r] = std::numeric_limits<Real>::max();
313  else if (stoichiometry(r, j) == 0.0)
314  _mineral_sat[r] *= 1.0;
315  else
316  {
317  _mineral_sat[r] = 0.0;
318  break;
319  }
320  }
321  else
322  _mineral_sat[r] *= std::pow(gamp, stoichiometry(r, j));
323  }
324  const Real fac = 1.0 - std::pow(_mineral_sat[r], _theta_exponent[r]);
325  // if fac > 0 then dissolution occurs; if fac < 0 then precipitation occurs.
326  const Real sgn = (fac < 0 ? -1.0 : 1.0);
327  const Real unbounded_rr = -sgn * rateConstantQp(r) * _r_area[r] * _molar_volume[r] *
328  std::pow(std::abs(fac), _eta_exponent[r]);
329 
330  /*
331  *
332  * Note the use of the OLD value of porosity here.
333  * This strategy, which breaks the cyclic dependency between porosity
334  * and mineral concentration, is used in
335  * Kernel: PorousFlowPreDis
336  * Material: PorousFlowPorosity
337  * Material: PorousFlowAqueousPreDisChemistry
338  * Material: PorousFlowAqueousPreDisMineral
339  *
340  */
341 
342  // bound the reaction so _sec_conc lies between zero and unity
343  const Real por_times_rr_dt = _porosity_old[_qp] * _saturation[_qp][_aq_ph] * unbounded_rr * _dt;
344  if (_sec_conc_old[_qp][r] + por_times_rr_dt > 1.0)
345  {
346  _bounded_rate[r] = true;
347  _reaction_rate[_qp][r] =
348  (1.0 - _sec_conc_old[_qp][r]) / _porosity_old[_qp] / _saturation[_qp][_aq_ph] / _dt;
349  }
350  else if (_sec_conc_old[_qp][r] + por_times_rr_dt < 0.0)
351  {
352  _bounded_rate[r] = true;
353  _reaction_rate[_qp][r] =
354  -_sec_conc_old[_qp][r] / _porosity_old[_qp] / _saturation[_qp][_aq_ph] / _dt;
355  }
356  else
357  {
358  _bounded_rate[r] = false;
359  _reaction_rate[_qp][r] = unbounded_rr;
360  }
361  }
362 }
363 
364 void
366  std::vector<Real> & drr) const
367 {
368  drr.assign(_num_primary, 0.0);
369 
370  // handle corner case
371  if (_bounded_rate[reaction_num])
372  return;
373 
374  /*
375  * Form the derivative of _mineral_sat, and store it in drr for now.
376  * The derivatives are straightforward if all primary > 0.
377  *
378  * If more than one primary = 0 then I set the derivatives to zero, even though it could be
379  * argued that with certain stoichiometric coefficients you might have derivative = 0/0 and it
380  * might be appropriate to set this to a non-zero finite value.
381  *
382  * If exactly one primary = 0 and its stoichiometry = 1 then the derivative wrt this one is
383  * nonzero.
384  * If exactly one primary = 0 and its stoichiometry > 1 then all derivatives are zero.
385  * If exactly one primary = 0 and its stoichiometry < 1 then the derivative wrt this one is
386  * infinity
387  */
388 
389  unsigned zero_count = 0;
390  unsigned zero_conc_index = 0;
391  findZeroConcentration(zero_conc_index, zero_count);
392  if (zero_count == 0)
393  {
394  for (unsigned i = 0; i < _num_primary; ++i)
395  drr[i] = stoichiometry(reaction_num, i) * _mineral_sat[reaction_num] /
396  std::max((*_primary[i])[_qp], 0.0);
397  }
398  else
399  {
400  if (_theta_exponent[reaction_num] < 1.0)
401  // dfac = infinity (see below) so the derivative may be inf, inf * 0, or inf/inf. I simply
402  // return with drr = 0
403  return;
404 
405  // count the number of primary <= 0, and record the one that's zero
406  if (zero_count == 1 && stoichiometry(reaction_num, zero_conc_index) == 1.0)
407  {
408  Real conc_without_zero = (_equilibrium_constants_as_log10
409  ? std::pow(10.0, -(*_equilibrium_constants[reaction_num])[_qp])
410  : 1.0 / (*_equilibrium_constants[reaction_num])[_qp]);
411  for (unsigned i = 0; i < _num_primary; ++i)
412  {
413  if (i == zero_conc_index)
414  conc_without_zero *= _primary_activity_coefficients[i];
415  else
416  conc_without_zero *=
417  std::pow(_primary_activity_coefficients[i] * std::max((*_primary[i])[_qp], 0.0),
418  stoichiometry(reaction_num, i));
419  }
420  drr[zero_conc_index] = conc_without_zero;
421  }
422  else if (zero_count == 0 and stoichiometry(reaction_num, zero_conc_index) < 1.0)
423  drr[zero_conc_index] = std::numeric_limits<Real>::max();
424  else
425  // all other cases have drr = 0, so return without performing any other calculations
426  return;
427  }
428 
429  // At the moment _drr = d(mineral_sat)/d(primary)
430  // Multiply by appropriate term to give _drr = d(reaction_rate)/d(primary)
431  const Real fac = 1.0 - std::pow(_mineral_sat[reaction_num], _theta_exponent[reaction_num]);
432  const Real dfac = -_theta_exponent[reaction_num] *
433  std::pow(_mineral_sat[reaction_num], _theta_exponent[reaction_num] - 1.0);
434  const Real multiplier = -rateConstantQp(reaction_num) * _r_area[reaction_num] *
435  _molar_volume[reaction_num] *
436  std::pow(std::abs(fac), _eta_exponent[reaction_num] - 1.0) *
437  _eta_exponent[reaction_num] * dfac;
438  for (unsigned i = 0; i < _num_primary; ++i)
439  drr[i] *= multiplier;
440 }
441 
442 Real
444 {
445  return _ref_kconst[reaction_num] * std::exp(_e_act[reaction_num] / _gas_const *
446  (_one_over_ref_temp - 1.0 / _temperature[_qp]));
447 }
448 
449 Real
451 {
452  // handle corner case
453  if (_bounded_rate[reaction_num])
454  return 0.0;
455 
456  const Real drate_const = rateConstantQp(reaction_num) * _e_act[reaction_num] / _gas_const *
457  std::pow(_temperature[_qp], -2.0);
458  const Real fac = 1.0 - std::pow(_mineral_sat[reaction_num], _theta_exponent[reaction_num]);
459  const Real sgn = (fac < 0 ? -1.0 : 1.0);
460  const Real dkinetic_rate = -sgn * drate_const * _r_area[reaction_num] *
461  _molar_volume[reaction_num] *
462  std::pow(std::abs(fac), _eta_exponent[reaction_num]);
463 
464  return dkinetic_rate;
465 }
const MaterialProperty< std::vector< Real > > & _sec_conc_old
InputParameters validParams< PorousFlowAqueousPreDisChemistry >()
const unsigned _num_equilibrium_constants
Number of equilibrium_constants provided.
std::vector< Real > _mineral_sat
Mineral saturation ratio - a useful temporary variable during computeQpProperties.
virtual void computeQpReactionRates()
Compute the secondary-species concentration as defined by the chemistry Must be overridden by derived...
const std::vector< Real > _reactions
Stoichiometry defining the aqeuous geochemistry equilibrium reactions.
const bool _equilibrium_constants_as_log10
Whether the equilibium constants are written in their log10 form, or in absolute terms.
MaterialProperty< std::vector< Real > > & _reaction_rate
Reaction rate of mineralisation.
Material designed to form a std::vector of mass fractions of mineral concentrations from primary-spec...
const std::vector< Real > _molar_volume
Molar volume (L/mol) for each secondary species.
PorousFlowAqueousPreDisChemistry(const InputParameters &parameters)
registerMooseObject("PorousFlowApp", PorousFlowAqueousPreDisChemistry)
MaterialProperty< std::vector< std::vector< Real > > > & _dreaction_rate_dvar
d(reaction rate of mineralisation)/d(porous flow var)
const std::vector< Real > _r_area
Reactive surface area (m^2/L) for each reaction.
const Real _gas_const
Gas constant (J/(mol K))
InputParameters validParams< PorousFlowMaterialVectorBase >()
const MaterialProperty< Real > & _temperature
Temperature.
Real rateConstantQp(unsigned reaction_num) const
Base class for all PorousFlow vector materials.
const unsigned int _num_components
Number of fluid components.
std::vector< const VariableValue * > _primary
Values of the primary species&#39; concentrations (dimensionless)
virtual Real dQpReactionRate_dT(unsigned reaction_num) const
Computes derivative of the reaction rate with respect to the temperature.
Real stoichiometry(unsigned reaction_num, unsigned primary_num) const
The stoichiometric coefficient.
const std::vector< Real > _eta_exponent
Eta exponent for the precipitation-dissolution for each reaction.
const std::vector< Real > _theta_exponent
Theta exponent for the precipitation-dissolution for each reaction.
const unsigned int _num_var
Number of PorousFlow variables.
const MaterialProperty< std::vector< Real > > & _saturation
Saturation.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
const unsigned int _num_primary
Number of primary species.
std::vector< bool > _bounded_rate
Whether the reaction rate has to be bounded in order that the precipitate stays inside [0...
void findZeroConcentration(unsigned &zero_conc_index, unsigned &zero_count) const
Checks gamp[i] = _primary_activity_coefficients[i] * (*_primary[i])[qp].
const std::vector< Real > _ref_kconst
Rate constant (mol/(m^2 s)) at reference temperature for each reaction.
std::vector< unsigned int > _primary_var_num
The variable number of the primary variables.
const std::vector< Real > _primary_activity_coefficients
Activity coefficients for the primary species (dimensionless)
virtual void dQpReactionRate_dprimary(unsigned reaction_num, std::vector< Real > &drr) const
Computes derivative of the reaction rate with respect to the primary concentrations.
const unsigned int _num_reactions
Number of equations in the aqueous geochemistry system.
const MaterialProperty< Real > & _porosity_old
Old values of the porosity.
const unsigned int _aq_ph
Aqueous phase number.
const std::vector< Real > _e_act
Activation energy (J/mol) for each reaction.
const Real _one_over_ref_temp
1/reference_temperature (1/K)
const MaterialProperty< std::vector< Real > > & _dtemperature_dvar
d(temperature)/(d porflow variable)
std::vector< const VariableValue * > _equilibrium_constants
Equilibrium constants (dimensionless)