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 "CoupledBEEquilibriumSub.h"
11 :
12 : registerMooseObject("ChemicalReactionsApp", CoupledBEEquilibriumSub);
13 :
14 : InputParameters
15 3070 : CoupledBEEquilibriumSub::validParams()
16 : {
17 3070 : InputParameters params = TimeDerivative::validParams();
18 6140 : params.addParam<Real>(
19 6140 : "weight", 1.0, "The weight of the equilibrium species in total concentration");
20 6140 : params.addCoupledVar("log_k", 0.0, "The equilibrium constant of this equilibrium species");
21 6140 : params.addParam<Real>(
22 : "sto_u",
23 6140 : 1.0,
24 : "The stoichiometric coefficient of the primary species this kernel operates on");
25 6140 : params.addCoupledVar(
26 : "gamma_u", 1.0, "Activity coefficient of primary species that this kernel operates on");
27 6140 : params.addParam<std::vector<Real>>(
28 : "sto_v", {}, "The stoichiometric coefficients of coupled primary species");
29 6140 : params.addCoupledVar("gamma_v", 1.0, "Activity coefficients of coupled primary species");
30 6140 : params.addCoupledVar("gamma_eq", 1.0, "Activity coefficient of this equilibrium species");
31 6140 : params.addCoupledVar("v", "Coupled primary species constituting the equilibrium species");
32 3070 : params.addClassDescription("Derivative of equilibrium species concentration wrt time");
33 3070 : return params;
34 0 : }
35 :
36 1660 : CoupledBEEquilibriumSub::CoupledBEEquilibriumSub(const InputParameters & parameters)
37 : : TimeDerivative(parameters),
38 1660 : _weight(getParam<Real>("weight")),
39 1660 : _log_k(coupledValue("log_k")),
40 3320 : _sto_u(getParam<Real>("sto_u")),
41 3320 : _sto_v(getParam<std::vector<Real>>("sto_v")),
42 1660 : _gamma_u(coupledValue("gamma_u")),
43 1660 : _gamma_u_old(coupledValueOld("gamma_u")),
44 4980 : _gamma_v(isCoupled("gamma_v")
45 1660 : ? coupledValues("gamma_v") // have value
46 3318 : : std::vector<const VariableValue *>(coupledComponents("v"),
47 4976 : &coupledValue("gamma_v"))), // default
48 4980 : _gamma_v_old(isCoupled("gamma_v")
49 1660 : ? coupledValuesOld("gamma_v") // have value
50 3318 : : std::vector<const VariableValue *>(coupledComponents("v"),
51 4976 : &coupledValue("gamma_v"))), // default
52 1660 : _gamma_eq(coupledValue("gamma_eq")),
53 1660 : _gamma_eq_old(coupledValueOld("gamma_eq")),
54 3320 : _porosity(getMaterialProperty<Real>("porosity")),
55 1660 : _vars(coupledIndices("v")),
56 1660 : _v_vals(coupledValues("v")),
57 1660 : _v_vals_old(coupledValuesOld("v")),
58 3320 : _u_old(valueOld())
59 : {
60 3320 : const unsigned int n = coupledComponents("v");
61 :
62 : // Check that the correct number of coupled values have been provided
63 1660 : if (_sto_v.size() != n)
64 2 : mooseError("The number of stoichiometric coefficients in sto_v is not equal to the number of "
65 : "coupled species in ",
66 2 : _name);
67 :
68 3316 : if (isCoupled("gamma_v"))
69 4 : if (coupledComponents("gamma_v") != n)
70 2 : mooseError("The number of activity coefficients in gamma_v is not equal to the number of "
71 : "coupled species in ",
72 2 : _name);
73 1656 : }
74 :
75 : Real
76 5002624 : CoupledBEEquilibriumSub::computeQpResidual()
77 : {
78 : mooseAssert(_gamma_eq[_qp] > 0.0, "Activity coefficient must be greater than zero");
79 :
80 : // Contribution due to primary species that this kernel acts on
81 : Real val_new =
82 5002624 : std::pow(10.0, _log_k[_qp]) * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u) / _gamma_eq[_qp];
83 5002624 : Real val_old = std::pow(10.0, _log_k[_qp]) * std::pow(_gamma_u_old[_qp] * _u_old[_qp], _sto_u) /
84 5002624 : _gamma_eq_old[_qp];
85 :
86 : // Contribution due to coupled primary species
87 8322432 : for (unsigned int i = 0; i < _vars.size(); ++i)
88 : {
89 3319808 : val_new *= std::pow((*_gamma_v[i])[_qp] * (*_v_vals[i])[_qp], _sto_v[i]);
90 3319808 : val_old *= std::pow((*_gamma_v_old[i])[_qp] * (*_v_vals_old[i])[_qp], _sto_v[i]);
91 : }
92 :
93 5002624 : return _porosity[_qp] * _weight * _test[_i][_qp] * (val_new - val_old) / _dt;
94 : }
95 :
96 : Real
97 3158912 : CoupledBEEquilibriumSub::computeQpJacobian()
98 : {
99 3158912 : Real val_new = std::pow(10.0, _log_k[_qp]) * _sto_u * _gamma_u[_qp] *
100 3158912 : std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _phi[_j][_qp] / _gamma_eq[_qp];
101 :
102 5286656 : for (unsigned int i = 0; i < _vars.size(); ++i)
103 2127744 : val_new *= std::pow((*_gamma_v[i])[_qp] * (*_v_vals[i])[_qp], _sto_v[i]);
104 :
105 3158912 : return _porosity[_qp] * _test[_i][_qp] * _weight * val_new / _dt;
106 : }
107 :
108 : Real
109 4043264 : CoupledBEEquilibriumSub::computeQpOffDiagJacobian(unsigned int jvar)
110 : {
111 : // If no coupled species, return 0
112 4043264 : if (_v_vals.size() == 0)
113 : return 0.0;
114 :
115 : // If jvar is not one of the coupled species, return 0
116 2896512 : if (std::find(_vars.begin(), _vars.end(), jvar) == _vars.end())
117 : return 0.0;
118 :
119 : Real val_new =
120 2127744 : std::pow(10.0, _log_k[_qp]) * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u) / _gamma_eq[_qp];
121 :
122 4454784 : for (unsigned int i = 0; i < _vars.size(); ++i)
123 : {
124 2327040 : if (jvar == _vars[i])
125 2127744 : val_new *= _sto_v[i] * (*_gamma_v[i])[_qp] *
126 2127744 : std::pow((*_gamma_v[i])[_qp] * (*_v_vals[i])[_qp], _sto_v[i] - 1.0) *
127 2127744 : _phi[_j][_qp];
128 : else
129 199296 : val_new *= std::pow((*_gamma_v[i])[_qp] * (*_v_vals[i])[_qp], _sto_v[i]);
130 : }
131 :
132 2127744 : return _porosity[_qp] * _test[_i][_qp] * _weight * val_new / _dt;
133 : }
|
