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 3066 : CoupledDiffusionReactionSub::validParams()
18 : {
19 3066 : InputParameters params = Kernel::validParams();
20 6132 : params.addParam<Real>(
21 : "weight",
22 6132 : 1.0,
23 : "Weight of equilibrium species concentration in the primary species concentration");
24 6132 : params.addCoupledVar(
25 : "log_k", 0.0, "Equilibrium constant of the equilbrium reaction in dissociation form");
26 6132 : params.addParam<Real>("sto_u",
27 6132 : 1.0,
28 : "Stoichiometric coef of the primary species this kernel "
29 : "operates on in the equilibrium reaction");
30 6132 : params.addCoupledVar(
31 : "gamma_u", 1.0, "Activity coefficient of primary species that this kernel operates on");
32 6132 : params.addParam<std::vector<Real>>(
33 : "sto_v", {}, "The stoichiometric coefficients of coupled primary species");
34 6132 : params.addCoupledVar("v", "List of coupled primary species in this equilibrium species");
35 6132 : params.addCoupledVar("gamma_v", 1.0, "Activity coefficients of coupled primary species");
36 6132 : params.addCoupledVar("gamma_eq", 1.0, "Activity coefficient of this equilibrium species");
37 3066 : params.addClassDescription("Diffusion of equilibrium species");
38 3066 : return params;
39 0 : }
40 :
41 1656 : CoupledDiffusionReactionSub::CoupledDiffusionReactionSub(const InputParameters & parameters)
42 : : Kernel(parameters),
43 1656 : _diffusivity(getMaterialProperty<Real>("diffusivity")),
44 3312 : _weight(getParam<Real>("weight")),
45 1656 : _log_k(coupledValue("log_k")),
46 3312 : _sto_u(getParam<Real>("sto_u")),
47 3312 : _sto_v(getParam<std::vector<Real>>("sto_v")),
48 1656 : _vars(coupledIndices("v")),
49 1656 : _vals(coupledValues("v")),
50 1656 : _grad_vals(coupledGradients("v")),
51 1656 : _gamma_u(coupledValue("gamma_u")),
52 4968 : _gamma_v(isCoupled("gamma_v")
53 1656 : ? coupledValues("gamma_v") // have value
54 3312 : : std::vector<const VariableValue *>(coupledComponents("v"),
55 4968 : &coupledValue("gamma_v"))), // default
56 3312 : _gamma_eq(coupledValue("gamma_eq"))
57 : {
58 3312 : const unsigned int n = coupledComponents("v");
59 :
60 : // Check that the correct number of coupled values have been provided
61 1656 : 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 3312 : 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 1656 : }
72 :
73 : Real
74 5002624 : CoupledDiffusionReactionSub::computeQpResidual()
75 : {
76 : RealGradient diff1 =
77 5002624 : _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_u[_qp];
78 8322432 : for (unsigned int i = 0; i < _vals.size(); ++i)
79 3319808 : 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 5002624 : const Real d_val = std::pow(_gamma_u[_qp] * _u[_qp], _sto_u);
83 8322432 : for (unsigned int i = 0; i < _vals.size(); ++i)
84 : {
85 3319808 : RealGradient diff2 = d_val * _sto_v[i] * (*_gamma_v[i])[_qp] *
86 3319808 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
87 3319808 : (*_grad_vals[i])[_qp];
88 :
89 6892288 : for (unsigned int j = 0; j < _vals.size(); ++j)
90 3572480 : if (j != i)
91 252672 : 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 5002624 : return _weight * std::pow(10.0, _log_k[_qp]) * _diffusivity[_qp] * _grad_test[_i][_qp] *
98 5002624 : (diff1 + diff2_sum) / _gamma_eq[_qp];
99 : ;
100 : }
101 :
102 : Real
103 3158912 : CoupledDiffusionReactionSub::computeQpJacobian()
104 : {
105 : RealGradient diff1_1 =
106 3158912 : _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_phi[_j][_qp];
107 3158912 : RealGradient diff1_2 = _phi[_j][_qp] * _sto_u * (_sto_u - 1.0) * _gamma_u[_qp] * _gamma_u[_qp] *
108 3158912 : std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 2.0) * _grad_u[_qp];
109 5286656 : for (unsigned int i = 0; i < _vals.size(); ++i)
110 : {
111 2127744 : diff1_1 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
112 2127744 : 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 3158912 : _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 5286656 : for (unsigned int i = 0; i < _vals.size(); ++i)
120 : {
121 2127744 : RealGradient diff2 = d_val * _sto_v[i] * (*_gamma_v[i])[_qp] *
122 2127744 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
123 2127744 : (*_grad_vals[i])[_qp];
124 4454784 : for (unsigned int j = 0; j < _vals.size(); ++j)
125 2327040 : if (j != i)
126 199296 : diff2 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[j])[_qp], _sto_v[j]);
127 :
128 : diff2_sum += diff2;
129 : }
130 :
131 3158912 : return _weight * std::pow(10.0, _log_k[_qp]) * _diffusivity[_qp] * _grad_test[_i][_qp] *
132 3158912 : (diff1 + diff2_sum) / _gamma_eq[_qp];
133 : }
134 :
135 : Real
136 4043264 : CoupledDiffusionReactionSub::computeQpOffDiagJacobian(unsigned int jvar)
137 : {
138 : // If no coupled species, return 0
139 4043264 : if (_vals.size() == 0)
140 : return 0.0;
141 :
142 : // If jvar is not one of the coupled species, return 0
143 2896512 : if (std::find(_vars.begin(), _vars.end(), jvar) == _vars.end())
144 : return 0.0;
145 :
146 : RealGradient diff1 =
147 2127744 : _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_u[_qp];
148 4454784 : for (unsigned int i = 0; i < _vals.size(); ++i)
149 : {
150 2327040 : if (jvar == _vars[i])
151 2127744 : diff1 *= _sto_v[i] * (*_gamma_v[i])[_qp] *
152 2127744 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) * _phi[_j][_qp];
153 : else
154 199296 : diff1 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
155 : }
156 :
157 2127744 : 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 4454784 : for (unsigned int i = 0; i < _vals.size(); ++i)
163 2327040 : if (jvar == _vars[i])
164 : {
165 2127744 : diff2_1 = _sto_v[i] * (_sto_v[i] - 1.0) * (*_gamma_v[i])[_qp] * (*_gamma_v[i])[_qp] *
166 2127744 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 2.0) * _phi[_j][_qp] *
167 2127744 : (*_grad_vals[i])[_qp];
168 2127744 : diff2_2 = _sto_v[i] * (*_gamma_v[i])[_qp] *
169 2127744 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
170 2127744 : _grad_phi[_j][_qp];
171 : }
172 :
173 : RealGradient diff2 = val_u * (diff2_1 + diff2_2);
174 :
175 4454784 : for (unsigned int i = 0; i < _vals.size(); ++i)
176 2327040 : if (jvar != _vars[i])
177 199296 : 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 4454784 : for (unsigned int i = 0; i < _vals.size(); ++i)
185 2327040 : if (jvar == _vars[i])
186 : {
187 : var = i;
188 2127744 : val_jvar = val_u * _sto_v[i] * (*_gamma_v[i])[_qp] *
189 2127744 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) * _phi[_j][_qp];
190 : }
191 :
192 4454784 : for (unsigned int i = 0; i < _vals.size(); ++i)
193 2327040 : if (i != var)
194 : {
195 199296 : diff3 = val_jvar * _sto_v[i] * (*_gamma_v[i])[_qp] *
196 199296 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
197 199296 : (*_grad_vals[i])[_qp];
198 :
199 597888 : for (unsigned int j = 0; j < _vals.size(); ++j)
200 398592 : 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 2127744 : return _weight * std::pow(10.0, _log_k[_qp]) * _diffusivity[_qp] * _grad_test[_i][_qp] *
207 2127744 : (diff1 + diff2 + diff3_sum) / _gamma_eq[_qp];
208 : }
|
