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