www.mooseframework.org
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
PorousFlowAqueousPreDisChemistry Class Reference

Material designed to form a std::vector of mass fractions of mineral concentrations from primary-species concentrations for an equilibrium precipitation-dissolution chemistry reaction system. More...

#include <PorousFlowAqueousPreDisChemistry.h>

Inheritance diagram for PorousFlowAqueousPreDisChemistry:
[legend]

Public Member Functions

 PorousFlowAqueousPreDisChemistry (const InputParameters &parameters)
 

Protected Member Functions

void initQpStatefulProperties () override
 
void computeQpProperties () override
 
Real stoichiometry (unsigned reaction_num, unsigned primary_num) const
 The stoichiometric coefficient. More...
 
virtual void computeQpReactionRates ()
 Compute the secondary-species concentration as defined by the chemistry Must be overridden by derived classes. More...
 
void findZeroConcentration (unsigned &zero_conc_index, unsigned &zero_count) const
 Checks gamp[i] = _primary_activity_coefficients[i] * (*_primary[i])[qp]. More...
 
virtual void dQpReactionRate_dprimary (unsigned reaction_num, std::vector< Real > &drr) const
 Computes derivative of the reaction rate with respect to the primary concentrations. More...
 
virtual Real dQpReactionRate_dT (unsigned reaction_num) const
 Computes derivative of the reaction rate with respect to the temperature. More...
 
Real rateConstantQp (unsigned reaction_num) const
 

Protected Attributes

const MaterialProperty< Real > & _porosity_old
 Old values of the porosity. More...
 
const unsigned int _aq_ph
 Aqueous phase number. More...
 
const MaterialProperty< std::vector< Real > > & _saturation
 Saturation. More...
 
const MaterialProperty< Real > & _temperature
 Temperature. More...
 
const MaterialProperty< std::vector< Real > > & _dtemperature_dvar
 d(temperature)/(d porflow variable) More...
 
const unsigned int _num_primary
 Number of primary species. More...
 
const unsigned int _num_reactions
 Number of equations in the aqueous geochemistry system. More...
 
const bool _equilibrium_constants_as_log10
 Whether the equilibium constants are written in their log10 form, or in absolute terms. More...
 
const unsigned _num_equilibrium_constants
 Number of equilibrium_constants provided. More...
 
std::vector< const VariableValue * > _equilibrium_constants
 Equilibrium constants (dimensionless) More...
 
const std::vector< Real > _primary_activity_coefficients
 Activity coefficients for the primary species (dimensionless) More...
 
const std::vector< Real > _reactions
 Stoichiometry defining the aqeuous geochemistry equilibrium reactions. More...
 
std::vector< unsigned int > _primary_var_num
 The variable number of the primary variables. More...
 
std::vector< const VariableValue * > _primary
 Values of the primary species' concentrations (dimensionless) More...
 
const MaterialProperty< std::vector< Real > > & _sec_conc_old
 
std::vector< Real > _mineral_sat
 Mineral saturation ratio - a useful temporary variable during computeQpProperties. More...
 
std::vector< bool > _bounded_rate
 Whether the reaction rate has to be bounded in order that the precipitate stays inside [0, 1]. More...
 
MaterialProperty< std::vector< Real > > & _reaction_rate
 Reaction rate of mineralisation. More...
 
MaterialProperty< std::vector< std::vector< Real > > > & _dreaction_rate_dvar
 d(reaction rate of mineralisation)/d(porous flow var) More...
 
const std::vector< Real > _r_area
 Reactive surface area (m^2/L) for each reaction. More...
 
const std::vector< Real > _molar_volume
 Molar volume (L/mol) for each secondary species. More...
 
const std::vector< Real > _ref_kconst
 Rate constant (mol/(m^2 s)) at reference temperature for each reaction. More...
 
const std::vector< Real > _e_act
 Activation energy (J/mol) for each reaction. More...
 
const Real _gas_const
 Gas constant (J/(mol K)) More...
 
const Real _one_over_ref_temp
 1/reference_temperature (1/K) More...
 
const std::vector< Real > _theta_exponent
 Theta exponent for the precipitation-dissolution for each reaction. More...
 
const std::vector< Real > _eta_exponent
 Eta exponent for the precipitation-dissolution for each reaction. More...
 
std::vector< const VariableValue * > _initial_conc
 Initial values of the secondary species concentrations. More...
 
const unsigned int _num_phases
 Number of phases. More...
 
const unsigned int _num_components
 Number of fluid components. More...
 
const unsigned int _num_var
 Number of PorousFlow variables. More...
 

Detailed Description

Material designed to form a std::vector of mass fractions of mineral concentrations from primary-species concentrations for an equilibrium precipitation-dissolution chemistry reaction system.

Definition at line 25 of file PorousFlowAqueousPreDisChemistry.h.

Constructor & Destructor Documentation

◆ PorousFlowAqueousPreDisChemistry()

PorousFlowAqueousPreDisChemistry::PorousFlowAqueousPreDisChemistry ( const InputParameters &  parameters)

Definition at line 70 of file PorousFlowAqueousPreDisChemistry.C.

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")),
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")),
92  _primary_activity_coefficients(getParam<std::vector<Real>>("primary_activity_coefficients")),
93  _reactions(getParam<std::vector<Real>>("reactions")),
94 
96  _nodal_material
97  ? getMaterialPropertyOld<std::vector<Real>>("PorousFlow_mineral_concentration_nodal")
98  : getMaterialPropertyOld<std::vector<Real>>("PorousFlow_mineral_concentration_qp")),
99 
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 ? &coupledNodalValue("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 ? &coupledNodalValue("equilibrium_constants", i)
222  : &coupledValue("equilibrium_constants", i));
223 }
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.
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.
const std::vector< Real > _molar_volume
Molar volume (L/mol) for each secondary species.
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))
const MaterialProperty< Real > & _temperature
Temperature.
const unsigned int _num_components
Number of fluid components.
std::vector< const VariableValue * > _primary
Values of the primary species&#39; concentrations (dimensionless)
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 MaterialProperty< std::vector< Real > > & _saturation
Saturation.
PorousFlowMaterialVectorBase(const InputParameters &parameters)
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...
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)
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)

Member Function Documentation

◆ computeQpProperties()

void PorousFlowAqueousPreDisChemistry::computeQpProperties ( )
overrideprotected

Definition at line 235 of file PorousFlowAqueousPreDisChemistry.C.

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 }
virtual void computeQpReactionRates()
Compute the secondary-species concentration as defined by the chemistry Must be overridden by derived...
MaterialProperty< std::vector< Real > > & _reaction_rate
Reaction rate of mineralisation.
MaterialProperty< std::vector< std::vector< Real > > > & _dreaction_rate_dvar
d(reaction rate of mineralisation)/d(porous flow var)
virtual Real dQpReactionRate_dT(unsigned reaction_num) const
Computes derivative of the reaction rate with respect to the temperature.
const unsigned int _num_var
Number of PorousFlow variables.
const unsigned int _num_primary
Number of primary species.
std::vector< unsigned int > _primary_var_num
The variable number of the primary variables.
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< std::vector< Real > > & _dtemperature_dvar
d(temperature)/(d porflow variable)

◆ computeQpReactionRates()

void PorousFlowAqueousPreDisChemistry::computeQpReactionRates ( )
protectedvirtual

Compute the secondary-species concentration as defined by the chemistry Must be overridden by derived classes.

Definition at line 299 of file PorousFlowAqueousPreDisChemistry.C.

Referenced by computeQpProperties().

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 }
const MaterialProperty< std::vector< Real > > & _sec_conc_old
std::vector< Real > _mineral_sat
Mineral saturation ratio - a useful temporary variable during computeQpProperties.
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.
const std::vector< Real > _molar_volume
Molar volume (L/mol) for each secondary species.
const std::vector< Real > _r_area
Reactive surface area (m^2/L) for each reaction.
Real rateConstantQp(unsigned reaction_num) const
std::vector< const VariableValue * > _primary
Values of the primary species&#39; concentrations (dimensionless)
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 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...
const std::vector< Real > _primary_activity_coefficients
Activity coefficients for the primary species (dimensionless)
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.
std::vector< const VariableValue * > _equilibrium_constants
Equilibrium constants (dimensionless)

◆ dQpReactionRate_dprimary()

void PorousFlowAqueousPreDisChemistry::dQpReactionRate_dprimary ( unsigned  reaction_num,
std::vector< Real > &  drr 
) const
protectedvirtual

Computes derivative of the reaction rate with respect to the primary concentrations.

Parameters
reaction_numThe reaction number corresponding to the secondary-species concentration
drrdrr[i] = d(reactionRate[reaction_num])/d(primary_species[i])

Definition at line 365 of file PorousFlowAqueousPreDisChemistry.C.

Referenced by computeQpProperties().

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 }
std::vector< Real > _mineral_sat
Mineral saturation ratio - a useful temporary variable during computeQpProperties.
const bool _equilibrium_constants_as_log10
Whether the equilibium constants are written in their log10 form, or in absolute terms.
const std::vector< Real > _molar_volume
Molar volume (L/mol) for each secondary species.
const std::vector< Real > _r_area
Reactive surface area (m^2/L) for each reaction.
Real rateConstantQp(unsigned reaction_num) const
std::vector< const VariableValue * > _primary
Values of the primary species&#39; concentrations (dimensionless)
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.
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 > _primary_activity_coefficients
Activity coefficients for the primary species (dimensionless)
std::vector< const VariableValue * > _equilibrium_constants
Equilibrium constants (dimensionless)

◆ dQpReactionRate_dT()

Real PorousFlowAqueousPreDisChemistry::dQpReactionRate_dT ( unsigned  reaction_num) const
protectedvirtual

Computes derivative of the reaction rate with respect to the temperature.

Parameters
reaction_numThe reaction number corresponding to the secondary-species concentration

Definition at line 450 of file PorousFlowAqueousPreDisChemistry.C.

Referenced by computeQpProperties().

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 }
std::vector< Real > _mineral_sat
Mineral saturation ratio - a useful temporary variable during computeQpProperties.
const std::vector< Real > _molar_volume
Molar volume (L/mol) for each secondary species.
const std::vector< Real > _r_area
Reactive surface area (m^2/L) for each reaction.
const Real _gas_const
Gas constant (J/(mol K))
const MaterialProperty< Real > & _temperature
Temperature.
Real rateConstantQp(unsigned reaction_num) const
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.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
std::vector< bool > _bounded_rate
Whether the reaction rate has to be bounded in order that the precipitate stays inside [0...
const std::vector< Real > _e_act
Activation energy (J/mol) for each reaction.

◆ findZeroConcentration()

void PorousFlowAqueousPreDisChemistry::findZeroConcentration ( unsigned &  zero_conc_index,
unsigned &  zero_count 
) const
protected

Checks gamp[i] = _primary_activity_coefficients[i] * (*_primary[i])[qp].

Returns: if all of these are positive, then zero_count = 0, zero_conc_index = 0 if one of these is zero, then zero_count = 1, zero_conc_index = the index of the zero gamp if more than one is zero, then zero_count = 2, and zero_conc_index is the index of the 2nd zero

Definition at line 281 of file PorousFlowAqueousPreDisChemistry.C.

Referenced by dQpReactionRate_dprimary().

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 }
std::vector< const VariableValue * > _primary
Values of the primary species&#39; concentrations (dimensionless)
const unsigned int _num_primary
Number of primary species.
const std::vector< Real > _primary_activity_coefficients
Activity coefficients for the primary species (dimensionless)

◆ initQpStatefulProperties()

void PorousFlowAqueousPreDisChemistry::initQpStatefulProperties ( )
overrideprotected

Definition at line 226 of file PorousFlowAqueousPreDisChemistry.C.

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 }
MaterialProperty< std::vector< Real > > & _reaction_rate
Reaction rate of mineralisation.
MaterialProperty< std::vector< std::vector< Real > > > & _dreaction_rate_dvar
d(reaction rate of mineralisation)/d(porous flow var)
const unsigned int _num_var
Number of PorousFlow variables.
const unsigned int _num_reactions
Number of equations in the aqueous geochemistry system.

◆ rateConstantQp()

Real PorousFlowAqueousPreDisChemistry::rateConstantQp ( unsigned  reaction_num) const
protected

Definition at line 443 of file PorousFlowAqueousPreDisChemistry.C.

Referenced by computeQpReactionRates(), dQpReactionRate_dprimary(), and dQpReactionRate_dT().

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 }
const Real _gas_const
Gas constant (J/(mol K))
const MaterialProperty< Real > & _temperature
Temperature.
const std::vector< Real > _ref_kconst
Rate constant (mol/(m^2 s)) at reference temperature for each reaction.
const std::vector< Real > _e_act
Activation energy (J/mol) for each reaction.
const Real _one_over_ref_temp
1/reference_temperature (1/K)

◆ stoichiometry()

Real PorousFlowAqueousPreDisChemistry::stoichiometry ( unsigned  reaction_num,
unsigned  primary_num 
) const
protected

The stoichiometric coefficient.

Parameters
reaction_numReaction number (0, ..., _num_reactions - 1)
primary_numThe number of the primary species (0, ..., _num_primary - 1)

Definition at line 274 of file PorousFlowAqueousPreDisChemistry.C.

Referenced by computeQpReactionRates(), and dQpReactionRate_dprimary().

275 {
276  const unsigned index = reaction_num * _num_primary + primary_num;
277  return _reactions[index];
278 }
const std::vector< Real > _reactions
Stoichiometry defining the aqeuous geochemistry equilibrium reactions.
const unsigned int _num_primary
Number of primary species.

Member Data Documentation

◆ _aq_ph

const unsigned int PorousFlowAqueousPreDisChemistry::_aq_ph
protected

Aqueous phase number.

Definition at line 75 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpReactionRates().

◆ _bounded_rate

std::vector<bool> PorousFlowAqueousPreDisChemistry::_bounded_rate
protected

Whether the reaction rate has to be bounded in order that the precipitate stays inside [0, 1].

Definition at line 120 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpReactionRates(), dQpReactionRate_dprimary(), and dQpReactionRate_dT().

◆ _dreaction_rate_dvar

MaterialProperty<std::vector<std::vector<Real> > >& PorousFlowAqueousPreDisChemistry::_dreaction_rate_dvar
protected

d(reaction rate of mineralisation)/d(porous flow var)

Definition at line 126 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpProperties(), and initQpStatefulProperties().

◆ _dtemperature_dvar

const MaterialProperty<std::vector<Real> >& PorousFlowAqueousPreDisChemistry::_dtemperature_dvar
protected

d(temperature)/(d porflow variable)

Definition at line 84 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpProperties().

◆ _e_act

const std::vector<Real> PorousFlowAqueousPreDisChemistry::_e_act
protected

Activation energy (J/mol) for each reaction.

Definition at line 138 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by dQpReactionRate_dT(), PorousFlowAqueousPreDisChemistry(), and rateConstantQp().

◆ _equilibrium_constants

std::vector<const VariableValue *> PorousFlowAqueousPreDisChemistry::_equilibrium_constants
protected

Equilibrium constants (dimensionless)

Definition at line 99 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpReactionRates(), dQpReactionRate_dprimary(), and PorousFlowAqueousPreDisChemistry().

◆ _equilibrium_constants_as_log10

const bool PorousFlowAqueousPreDisChemistry::_equilibrium_constants_as_log10
protected

Whether the equilibium constants are written in their log10 form, or in absolute terms.

Definition at line 93 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpReactionRates(), and dQpReactionRate_dprimary().

◆ _eta_exponent

const std::vector<Real> PorousFlowAqueousPreDisChemistry::_eta_exponent
protected

Eta exponent for the precipitation-dissolution for each reaction.

Definition at line 150 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpReactionRates(), dQpReactionRate_dprimary(), dQpReactionRate_dT(), and PorousFlowAqueousPreDisChemistry().

◆ _gas_const

const Real PorousFlowAqueousPreDisChemistry::_gas_const
protected

Gas constant (J/(mol K))

Definition at line 141 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by dQpReactionRate_dT(), and rateConstantQp().

◆ _initial_conc

std::vector<const VariableValue *> PorousFlowAqueousPreDisChemistry::_initial_conc
protected

Initial values of the secondary species concentrations.

Definition at line 153 of file PorousFlowAqueousPreDisChemistry.h.

◆ _mineral_sat

std::vector<Real> PorousFlowAqueousPreDisChemistry::_mineral_sat
protected

Mineral saturation ratio - a useful temporary variable during computeQpProperties.

Definition at line 117 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpReactionRates(), dQpReactionRate_dprimary(), and dQpReactionRate_dT().

◆ _molar_volume

const std::vector<Real> PorousFlowAqueousPreDisChemistry::_molar_volume
protected

Molar volume (L/mol) for each secondary species.

Definition at line 132 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpReactionRates(), dQpReactionRate_dprimary(), dQpReactionRate_dT(), and PorousFlowAqueousPreDisChemistry().

◆ _num_components

const unsigned int PorousFlowMaterialVectorBase::_num_components
protectedinherited

◆ _num_equilibrium_constants

const unsigned PorousFlowAqueousPreDisChemistry::_num_equilibrium_constants
protected

Number of equilibrium_constants provided.

Definition at line 96 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by PorousFlowAqueousPreDisChemistry().

◆ _num_phases

const unsigned int PorousFlowMaterialVectorBase::_num_phases
protectedinherited

◆ _num_primary

const unsigned int PorousFlowAqueousPreDisChemistry::_num_primary
protected

◆ _num_reactions

const unsigned int PorousFlowAqueousPreDisChemistry::_num_reactions
protected

Number of equations in the aqueous geochemistry system.

Definition at line 90 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpProperties(), computeQpReactionRates(), initQpStatefulProperties(), and PorousFlowAqueousPreDisChemistry().

◆ _num_var

const unsigned int PorousFlowMaterialVectorBase::_num_var
protectedinherited

◆ _one_over_ref_temp

const Real PorousFlowAqueousPreDisChemistry::_one_over_ref_temp
protected

1/reference_temperature (1/K)

Definition at line 144 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by rateConstantQp().

◆ _porosity_old

const MaterialProperty<Real>& PorousFlowAqueousPreDisChemistry::_porosity_old
protected

Old values of the porosity.

Definition at line 72 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpReactionRates().

◆ _primary

std::vector<const VariableValue *> PorousFlowAqueousPreDisChemistry::_primary
protected

Values of the primary species' concentrations (dimensionless)

Definition at line 111 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpReactionRates(), dQpReactionRate_dprimary(), findZeroConcentration(), and PorousFlowAqueousPreDisChemistry().

◆ _primary_activity_coefficients

const std::vector<Real> PorousFlowAqueousPreDisChemistry::_primary_activity_coefficients
protected

Activity coefficients for the primary species (dimensionless)

Definition at line 102 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpReactionRates(), dQpReactionRate_dprimary(), findZeroConcentration(), and PorousFlowAqueousPreDisChemistry().

◆ _primary_var_num

std::vector<unsigned int> PorousFlowAqueousPreDisChemistry::_primary_var_num
protected

The variable number of the primary variables.

Definition at line 108 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpProperties(), and PorousFlowAqueousPreDisChemistry().

◆ _r_area

const std::vector<Real> PorousFlowAqueousPreDisChemistry::_r_area
protected

Reactive surface area (m^2/L) for each reaction.

Definition at line 129 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpReactionRates(), dQpReactionRate_dprimary(), dQpReactionRate_dT(), and PorousFlowAqueousPreDisChemistry().

◆ _reaction_rate

MaterialProperty<std::vector<Real> >& PorousFlowAqueousPreDisChemistry::_reaction_rate
protected

Reaction rate of mineralisation.

Definition at line 123 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpProperties(), computeQpReactionRates(), and initQpStatefulProperties().

◆ _reactions

const std::vector<Real> PorousFlowAqueousPreDisChemistry::_reactions
protected

Stoichiometry defining the aqeuous geochemistry equilibrium reactions.

Definition at line 105 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by PorousFlowAqueousPreDisChemistry(), and stoichiometry().

◆ _ref_kconst

const std::vector<Real> PorousFlowAqueousPreDisChemistry::_ref_kconst
protected

Rate constant (mol/(m^2 s)) at reference temperature for each reaction.

Definition at line 135 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by PorousFlowAqueousPreDisChemistry(), and rateConstantQp().

◆ _saturation

const MaterialProperty<std::vector<Real> >& PorousFlowAqueousPreDisChemistry::_saturation
protected

Saturation.

Definition at line 78 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpReactionRates().

◆ _sec_conc_old

const MaterialProperty<std::vector<Real> >& PorousFlowAqueousPreDisChemistry::_sec_conc_old
protected

Definition at line 114 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpReactionRates().

◆ _temperature

const MaterialProperty<Real>& PorousFlowAqueousPreDisChemistry::_temperature
protected

Temperature.

Definition at line 81 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by dQpReactionRate_dT(), and rateConstantQp().

◆ _theta_exponent

const std::vector<Real> PorousFlowAqueousPreDisChemistry::_theta_exponent
protected

Theta exponent for the precipitation-dissolution for each reaction.

Definition at line 147 of file PorousFlowAqueousPreDisChemistry.h.

Referenced by computeQpReactionRates(), dQpReactionRate_dprimary(), dQpReactionRate_dT(), and PorousFlowAqueousPreDisChemistry().


The documentation for this class was generated from the following files: