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 "CoupledDiffusionReactionSub.h"
11 :
12 : using libMesh::RealGradient;
13 :
14 : registerMooseObject("ChemicalReactionsApp", CoupledDiffusionReactionSub);
15 :
16 : InputParameters
17 1676 : CoupledDiffusionReactionSub::validParams()
18 : {
19 1676 : InputParameters params = Kernel::validParams();
20 3352 : params.addParam<Real>(
21 : "weight",
22 3352 : 1.0,
23 : "Weight of equilibrium species concentration in the primary species concentration");
24 3352 : params.addCoupledVar(
25 : "log_k", 0.0, "Equilibrium constant of the equilibrium reaction in dissociation form");
26 3352 : params.addParam<Real>("sto_u",
27 3352 : 1.0,
28 : "Stoichiometric coef of the primary species this kernel "
29 : "operates on in the equilibrium reaction");
30 3352 : params.addCoupledVar(
31 : "gamma_u", 1.0, "Activity coefficient of primary species that this kernel operates on");
32 3352 : params.addParam<std::vector<Real>>(
33 : "sto_v", {}, "The stoichiometric coefficients of coupled primary species");
34 3352 : params.addCoupledVar("v", "List of coupled primary species in this equilibrium species");
35 3352 : params.addCoupledVar("gamma_v", 1.0, "Activity coefficients of coupled primary species");
36 3352 : params.addCoupledVar("gamma_eq", 1.0, "Activity coefficient of this equilibrium species");
37 1676 : params.addClassDescription("Diffusion of equilibrium species");
38 1676 : return params;
39 0 : }
40 :
41 884 : CoupledDiffusionReactionSub::CoupledDiffusionReactionSub(const InputParameters & parameters)
42 : : Kernel(parameters),
43 884 : _diffusivity(getMaterialProperty<Real>("diffusivity")),
44 1768 : _weight(getParam<Real>("weight")),
45 884 : _log_k(coupledValue("log_k")),
46 1768 : _sto_u(getParam<Real>("sto_u")),
47 1768 : _sto_v(getParam<std::vector<Real>>("sto_v")),
48 884 : _vars(coupledIndices("v")),
49 884 : _vals(coupledValues("v")),
50 884 : _grad_vals(coupledGradients("v")),
51 884 : _gamma_u(coupledValue("gamma_u")),
52 2652 : _gamma_v(isCoupled("gamma_v")
53 884 : ? coupledValues("gamma_v") // have value
54 1768 : : std::vector<const VariableValue *>(coupledComponents("v"),
55 2652 : &coupledValue("gamma_v"))), // default
56 1768 : _gamma_eq(coupledValue("gamma_eq"))
57 : {
58 1768 : const unsigned int n = coupledComponents("v");
59 :
60 : // Check that the correct number of coupled values have been provided
61 884 : if (_sto_v.size() != n)
62 0 : mooseError("The number of stoichiometric coefficients in sto_v is not equal to the number of "
63 : "coupled species in ",
64 0 : _name);
65 :
66 1768 : if (isCoupled("gamma_v"))
67 0 : if (coupledComponents("gamma_v") != n)
68 0 : mooseError("The number of activity coefficients in gamma_v is not equal to the number of "
69 : "coupled species in ",
70 0 : _name);
71 884 : }
72 :
73 : Real
74 4396336 : CoupledDiffusionReactionSub::computeQpResidual()
75 : {
76 : RealGradient diff1 =
77 4396336 : _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_u[_qp];
78 7105552 : for (unsigned int i = 0; i < _vals.size(); ++i)
79 2709216 : diff1 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
80 :
81 : RealGradient diff2_sum(0.0, 0.0, 0.0);
82 4396336 : const Real d_val = std::pow(_gamma_u[_qp] * _u[_qp], _sto_u);
83 7105552 : for (unsigned int i = 0; i < _vals.size(); ++i)
84 : {
85 2709216 : RealGradient diff2 = d_val * _sto_v[i] * (*_gamma_v[i])[_qp] *
86 2709216 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
87 2709216 : (*_grad_vals[i])[_qp];
88 :
89 5587872 : for (unsigned int j = 0; j < _vals.size(); ++j)
90 2878656 : if (j != i)
91 169440 : diff2 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[j])[_qp], _sto_v[j]);
92 :
93 : diff2_sum += diff2;
94 : }
95 :
96 : mooseAssert(_gamma_eq[_qp] > 0.0, "Activity coefficient must be greater than zero");
97 4396336 : return _weight * std::pow(10.0, _log_k[_qp]) * _diffusivity[_qp] * _grad_test[_i][_qp] *
98 4396336 : (diff1 + diff2_sum) / _gamma_eq[_qp];
99 : ;
100 : }
101 :
102 : Real
103 2829440 : CoupledDiffusionReactionSub::computeQpJacobian()
104 : {
105 : RealGradient diff1_1 =
106 2829440 : _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_phi[_j][_qp];
107 2829440 : RealGradient diff1_2 = _phi[_j][_qp] * _sto_u * (_sto_u - 1.0) * _gamma_u[_qp] * _gamma_u[_qp] *
108 2829440 : std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 2.0) * _grad_u[_qp];
109 4585088 : for (unsigned int i = 0; i < _vals.size(); ++i)
110 : {
111 1755648 : diff1_1 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
112 1755648 : diff1_2 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
113 : }
114 :
115 : RealGradient diff1 = diff1_1 + diff1_2;
116 : Real d_val =
117 2829440 : _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _phi[_j][_qp];
118 : RealGradient diff2_sum(0.0, 0.0, 0.0);
119 4585088 : for (unsigned int i = 0; i < _vals.size(); ++i)
120 : {
121 1755648 : RealGradient diff2 = d_val * _sto_v[i] * (*_gamma_v[i])[_qp] *
122 1755648 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
123 1755648 : (*_grad_vals[i])[_qp];
124 3644928 : for (unsigned int j = 0; j < _vals.size(); ++j)
125 1889280 : if (j != i)
126 133632 : diff2 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[j])[_qp], _sto_v[j]);
127 :
128 : diff2_sum += diff2;
129 : }
130 :
131 2829440 : return _weight * std::pow(10.0, _log_k[_qp]) * _diffusivity[_qp] * _grad_test[_i][_qp] *
132 2829440 : (diff1 + diff2_sum) / _gamma_eq[_qp];
133 : }
134 :
135 : Real
136 3698176 : CoupledDiffusionReactionSub::computeQpOffDiagJacobian(unsigned int jvar)
137 : {
138 : // If no coupled species, return 0
139 3698176 : if (_vals.size() == 0)
140 : return 0.0;
141 :
142 : // If jvar is not one of the coupled species, return 0
143 2600448 : if (std::find(_vars.begin(), _vars.end(), jvar) == _vars.end())
144 : return 0.0;
145 :
146 : RealGradient diff1 =
147 1755648 : _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_u[_qp];
148 3644928 : for (unsigned int i = 0; i < _vals.size(); ++i)
149 : {
150 1889280 : if (jvar == _vars[i])
151 1755648 : diff1 *= _sto_v[i] * (*_gamma_v[i])[_qp] *
152 1755648 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) * _phi[_j][_qp];
153 : else
154 133632 : diff1 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
155 : }
156 :
157 1755648 : Real val_u = std::pow(_gamma_u[_qp] * _u[_qp], _sto_u);
158 :
159 : RealGradient diff2_1(1.0, 1.0, 1.0);
160 : RealGradient diff2_2(1.0, 1.0, 1.0);
161 :
162 3644928 : for (unsigned int i = 0; i < _vals.size(); ++i)
163 1889280 : if (jvar == _vars[i])
164 : {
165 1755648 : diff2_1 = _sto_v[i] * (_sto_v[i] - 1.0) * (*_gamma_v[i])[_qp] * (*_gamma_v[i])[_qp] *
166 1755648 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 2.0) * _phi[_j][_qp] *
167 1755648 : (*_grad_vals[i])[_qp];
168 1755648 : diff2_2 = _sto_v[i] * (*_gamma_v[i])[_qp] *
169 1755648 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
170 1755648 : _grad_phi[_j][_qp];
171 : }
172 :
173 : RealGradient diff2 = val_u * (diff2_1 + diff2_2);
174 :
175 3644928 : for (unsigned int i = 0; i < _vals.size(); ++i)
176 1889280 : if (jvar != _vars[i])
177 133632 : diff2 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
178 :
179 : RealGradient diff3;
180 : RealGradient diff3_sum(0.0, 0.0, 0.0);
181 : Real val_jvar = 0.0;
182 : unsigned int var = 0;
183 :
184 3644928 : for (unsigned int i = 0; i < _vals.size(); ++i)
185 1889280 : if (jvar == _vars[i])
186 : {
187 : var = i;
188 1755648 : val_jvar = val_u * _sto_v[i] * (*_gamma_v[i])[_qp] *
189 1755648 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) * _phi[_j][_qp];
190 : }
191 :
192 3644928 : for (unsigned int i = 0; i < _vals.size(); ++i)
193 1889280 : if (i != var)
194 : {
195 133632 : diff3 = val_jvar * _sto_v[i] * (*_gamma_v[i])[_qp] *
196 133632 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
197 133632 : (*_grad_vals[i])[_qp];
198 :
199 400896 : for (unsigned int j = 0; j < _vals.size(); ++j)
200 267264 : if (j != var && j != i)
201 0 : diff3 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[j])[_qp], _sto_v[j]);
202 :
203 : diff3_sum += diff3;
204 : }
205 :
206 1755648 : return _weight * std::pow(10.0, _log_k[_qp]) * _diffusivity[_qp] * _grad_test[_i][_qp] *
207 1755648 : (diff1 + diff2 + diff3_sum) / _gamma_eq[_qp];
208 : }
|
