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 "CoupledConvectionReactionSub.h"
11 :
12 : using libMesh::RealGradient;
13 :
14 : registerMooseObject("ChemicalReactionsApp", CoupledConvectionReactionSub);
15 :
16 : InputParameters
17 3066 : CoupledConvectionReactionSub::validParams()
18 : {
19 3066 : InputParameters params = Kernel::validParams();
20 6132 : params.addParam<Real>("weight", 1.0, "Weight of the equilibrium species");
21 6132 : params.addCoupledVar("log_k", 0.0, "Equilibrium constant of dissociation equilibrium reaction");
22 6132 : params.addParam<Real>("sto_u",
23 6132 : 1.0,
24 : "Stoichiometric coef of the primary species the kernel "
25 : "operates on in the equilibrium reaction");
26 6132 : params.addCoupledVar(
27 : "gamma_u", 1.0, "Activity coefficient of primary species that this kernel operates on");
28 6132 : params.addParam<std::vector<Real>>(
29 : "sto_v",
30 : {},
31 : "The stoichiometric coefficients of coupled primary species in equilibrium reaction");
32 6132 : params.addRequiredCoupledVar("p", "Pressure");
33 6132 : params.addCoupledVar("v", "List of coupled primary species");
34 6132 : params.addCoupledVar("gamma_v", 1.0, "Activity coefficients of coupled primary species");
35 6132 : params.addCoupledVar("gamma_eq", 1.0, "Activity coefficient of this equilibrium species");
36 : RealVectorValue g(0, 0, 0);
37 6132 : params.addParam<RealVectorValue>("gravity", g, "Gravity vector (default is (0, 0, 0))");
38 3066 : params.addClassDescription("Convection of equilibrium species");
39 3066 : return params;
40 0 : }
41 :
42 1656 : CoupledConvectionReactionSub::CoupledConvectionReactionSub(const InputParameters & parameters)
43 : : DerivativeMaterialInterface<Kernel>(parameters),
44 1656 : _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 3312 : _cond(getMaterialProperty<Real>("conductivity")),
49 3312 : _gravity(getParam<RealVectorValue>("gravity")),
50 1656 : _density(getDefaultMaterialProperty<Real>("density")),
51 1656 : _grad_p(coupledGradient("p")),
52 1656 : _pvar(coupled("p")),
53 1656 : _vars(coupledIndices("v")),
54 1656 : _vals(coupledValues("v")),
55 1656 : _grad_vals(coupledGradients("v")),
56 1656 : _gamma_u(coupledValue("gamma_u")),
57 4968 : _gamma_v(isCoupled("gamma_v")
58 1656 : ? coupledValues("gamma_v") // have value
59 3312 : : std::vector<const VariableValue *>(coupledComponents("v"),
60 4968 : &coupledValue("gamma_v"))), // default
61 3312 : _gamma_eq(coupledValue("gamma_eq"))
62 : {
63 3312 : const unsigned int n = coupledComponents("v");
64 :
65 : // Check that the correct number of coupled values have been provided
66 1656 : if (_sto_v.size() != n)
67 0 : mooseError("The number of stoichiometric coefficients in sto_v is not equal to the number of "
68 : "coupled species in ",
69 0 : _name);
70 :
71 3312 : if (isCoupled("gamma_v"))
72 0 : if (coupledComponents("gamma_v") != n)
73 0 : mooseError("The number of activity coefficients in gamma_v is not equal to the number of "
74 : "coupled species in ",
75 0 : _name);
76 1656 : }
77 :
78 : Real
79 5002624 : CoupledConvectionReactionSub::computeQpResidual()
80 : {
81 5002624 : RealVectorValue darcy_vel = -_cond[_qp] * (_grad_p[_qp] - _density[_qp] * _gravity);
82 : RealGradient d_u =
83 5002624 : _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_u[_qp];
84 : RealGradient d_var_sum(0.0, 0.0, 0.0);
85 5002624 : const Real d_v_u = std::pow(_gamma_u[_qp] * _u[_qp], _sto_u);
86 :
87 8322432 : for (unsigned int i = 0; i < _vals.size(); ++i)
88 : {
89 3319808 : d_u *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
90 :
91 3319808 : RealGradient d_var = d_v_u * _sto_v[i] * (*_gamma_v[i])[_qp] *
92 3319808 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
93 3319808 : (*_grad_vals[i])[_qp];
94 :
95 6892288 : for (unsigned int j = 0; j < _vals.size(); ++j)
96 3572480 : if (j != i)
97 252672 : d_var *= std::pow((*_gamma_v[i])[_qp] * (*_vals[j])[_qp], _sto_v[j]);
98 :
99 : d_var_sum += d_var;
100 : }
101 :
102 : mooseAssert(_gamma_eq[_qp] > 0.0, "Activity coefficient must be greater than zero");
103 5002624 : return _weight * std::pow(10.0, _log_k[_qp]) * _test[_i][_qp] * darcy_vel * (d_u + d_var_sum) /
104 5002624 : _gamma_eq[_qp];
105 : }
106 :
107 : Real
108 3158912 : CoupledConvectionReactionSub::computeQpJacobian()
109 : {
110 3158912 : RealVectorValue darcy_vel = -_cond[_qp] * (_grad_p[_qp] - _density[_qp] * _gravity);
111 :
112 : RealGradient d_u_1 =
113 3158912 : _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_phi[_j][_qp];
114 3158912 : RealGradient d_u_2 = _phi[_j][_qp] * _sto_u * (_sto_u - 1.0) * _gamma_u[_qp] * _gamma_u[_qp] *
115 3158912 : std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 2.0) * _grad_u[_qp];
116 :
117 : RealGradient d_var_sum(0.0, 0.0, 0.0);
118 : const Real d_v_u =
119 3158912 : _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _phi[_j][_qp];
120 :
121 5286656 : for (unsigned int i = 0; i < _vals.size(); ++i)
122 : {
123 2127744 : d_u_1 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
124 2127744 : d_u_2 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
125 :
126 2127744 : RealGradient d_var = d_v_u * _sto_v[i] * (*_gamma_v[i])[_qp] *
127 2127744 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
128 2127744 : (*_grad_vals[i])[_qp];
129 4454784 : for (unsigned int j = 0; j < _vals.size(); ++j)
130 2327040 : if (j != i)
131 199296 : d_var *= std::pow((*_gamma_v[i])[_qp] * (*_vals[j])[_qp], _sto_v[j]);
132 :
133 : d_var_sum += d_var;
134 : }
135 :
136 : RealGradient d_u_j = d_u_1 + d_u_2;
137 3158912 : return _weight * std::pow(10.0, _log_k[_qp]) * _test[_i][_qp] * darcy_vel * (d_u_j + d_var_sum) /
138 3158912 : _gamma_eq[_qp];
139 : }
140 :
141 : Real
142 4043264 : CoupledConvectionReactionSub::computeQpOffDiagJacobian(unsigned int jvar)
143 : {
144 4043264 : if (jvar == _pvar)
145 : {
146 753280 : RealVectorValue ddarcy_vel_dp = -_cond[_qp] * _grad_phi[_j][_qp];
147 :
148 : RealGradient d_u =
149 753280 : _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_u[_qp];
150 : RealGradient d_var_sum(0.0, 0.0, 0.0);
151 753280 : const Real d_v_u = std::pow(_gamma_u[_qp] * _u[_qp], _sto_u);
152 :
153 1256320 : for (unsigned int i = 0; i < _vals.size(); ++i)
154 : {
155 503040 : d_u *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
156 :
157 503040 : RealGradient d_var = d_v_u * _sto_v[i] * (*_gamma_v[i])[_qp] *
158 503040 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
159 503040 : (*_grad_vals[i])[_qp];
160 1006080 : for (unsigned int j = 0; j < _vals.size(); ++j)
161 503040 : if (j != i)
162 0 : d_var *= std::pow((*_gamma_v[i])[_qp] * (*_vals[j])[_qp], _sto_v[j]);
163 :
164 : d_var_sum += d_var;
165 : }
166 753280 : return _weight * std::pow(10.0, _log_k[_qp]) * _test[_i][_qp] * ddarcy_vel_dp *
167 753280 : (d_u + d_var_sum) / _gamma_eq[_qp];
168 : }
169 :
170 3289984 : if (_vals.size() == 0)
171 : return 0.0;
172 :
173 2393472 : RealVectorValue darcy_vel = -_cond[_qp] * (_grad_p[_qp] - _density[_qp] * _gravity);
174 : RealGradient diff1 =
175 2393472 : _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_u[_qp];
176 4986240 : for (unsigned int i = 0; i < _vals.size(); ++i)
177 : {
178 2592768 : if (jvar == _vars[i])
179 2127744 : diff1 *= _sto_v[i] * (*_gamma_v[i])[_qp] *
180 2127744 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) * _phi[_j][_qp];
181 : else
182 465024 : diff1 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
183 : }
184 :
185 2393472 : Real val_u = std::pow(_gamma_u[_qp] * _u[_qp], _sto_u);
186 : RealGradient diff2_1(1.0, 1.0, 1.0);
187 : RealGradient diff2_2(1.0, 1.0, 1.0);
188 4986240 : for (unsigned int i = 0; i < _vals.size(); ++i)
189 2592768 : if (jvar == _vars[i])
190 : {
191 2127744 : diff2_1 = _sto_v[i] * (*_gamma_v[i])[_qp] * (*_gamma_v[i])[_qp] * (_sto_v[i] - 1.0) *
192 2127744 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 2.0) * _phi[_j][_qp] *
193 2127744 : (*_grad_vals[i])[_qp];
194 2127744 : diff2_2 = _sto_v[i] * (*_gamma_v[i])[_qp] *
195 2127744 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
196 2127744 : _grad_phi[_j][_qp];
197 : }
198 :
199 : RealGradient diff2 = val_u * (diff2_1 + diff2_2);
200 4986240 : for (unsigned int i = 0; i < _vals.size(); ++i)
201 2592768 : if (jvar != _vars[i])
202 : {
203 465024 : diff2 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
204 : }
205 :
206 : RealGradient diff3;
207 : RealGradient diff3_sum(0.0, 0.0, 0.0);
208 : Real val_jvar = 0.0;
209 : unsigned int var = 0;
210 :
211 4986240 : for (unsigned int i = 0; i < _vals.size(); ++i)
212 2592768 : if (jvar == _vars[i])
213 : {
214 : var = i;
215 2127744 : val_jvar = val_u * _sto_v[i] * (*_gamma_v[i])[_qp] *
216 2127744 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) * _phi[_j][_qp];
217 : }
218 :
219 4986240 : for (unsigned int i = 0; i < _vals.size(); ++i)
220 2592768 : if (i != var)
221 : {
222 199296 : diff3 = val_jvar * _sto_v[i] * (*_gamma_v[i])[_qp] *
223 199296 : std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
224 199296 : (*_grad_vals[i])[_qp];
225 :
226 597888 : for (unsigned int j = 0; j < _vals.size(); ++j)
227 398592 : if (j != var && j != i)
228 0 : diff3 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[j])[_qp], _sto_v[j]);
229 :
230 : diff3_sum += diff3;
231 : }
232 :
233 2393472 : return _weight * std::pow(10.0, _log_k[_qp]) * _test[_i][_qp] * darcy_vel *
234 2393472 : (diff1 + diff2 + diff3_sum) / _gamma_eq[_qp];
235 : }
|
