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