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 1680 : CoupledBEEquilibriumSub::validParams()
16 : {
17 1680 : InputParameters params = TimeDerivative::validParams();
18 3360 : params.addParam<Real>(
19 3360 : "weight", 1.0, "The weight of the equilibrium species in total concentration");
20 3360 : params.addCoupledVar("log_k", 0.0, "The equilibrium constant of this equilibrium species");
21 3360 : params.addParam<Real>(
22 : "sto_u",
23 3360 : 1.0,
24 : "The stoichiometric coefficient of the primary species this kernel operates on");
25 3360 : params.addCoupledVar(
26 : "gamma_u", 1.0, "Activity coefficient of primary species that this kernel operates on");
27 3360 : params.addParam<std::vector<Real>>(
28 : "sto_v", {}, "The stoichiometric coefficients of coupled primary species");
29 3360 : params.addCoupledVar("gamma_v", 1.0, "Activity coefficients of coupled primary species");
30 3360 : params.addCoupledVar("gamma_eq", 1.0, "Activity coefficient of this equilibrium species");
31 3360 : params.addCoupledVar("v", "Coupled primary species constituting the equilibrium species");
32 1680 : params.addClassDescription("Derivative of equilibrium species concentration wrt time");
33 1680 : return params;
34 0 : }
35 :
36 888 : CoupledBEEquilibriumSub::CoupledBEEquilibriumSub(const InputParameters & parameters)
37 : : TimeDerivative(parameters),
38 888 : _weight(getParam<Real>("weight")),
39 888 : _log_k(coupledValue("log_k")),
40 1776 : _sto_u(getParam<Real>("sto_u")),
41 1776 : _sto_v(getParam<std::vector<Real>>("sto_v")),
42 888 : _gamma_u(coupledValue("gamma_u")),
43 888 : _gamma_u_old(coupledValueOld("gamma_u")),
44 2664 : _gamma_v(isCoupled("gamma_v")
45 888 : ? coupledValues("gamma_v") // have value
46 1774 : : std::vector<const VariableValue *>(coupledComponents("v"),
47 2660 : &coupledValue("gamma_v"))), // default
48 2664 : _gamma_v_old(isCoupled("gamma_v")
49 888 : ? coupledValuesOld("gamma_v") // have value
50 1774 : : std::vector<const VariableValue *>(coupledComponents("v"),
51 2660 : &coupledValue("gamma_v"))), // default
52 888 : _gamma_eq(coupledValue("gamma_eq")),
53 888 : _gamma_eq_old(coupledValueOld("gamma_eq")),
54 1776 : _porosity(getMaterialProperty<Real>("porosity")),
55 888 : _vars(coupledIndices("v")),
56 888 : _v_vals(coupledValues("v")),
57 888 : _v_vals_old(coupledValuesOld("v")),
58 1776 : _u_old(valueOld())
59 : {
60 1776 : const unsigned int n = coupledComponents("v");
61 :
62 : // Check that the correct number of coupled values have been provided
63 888 : 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 1772 : 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 884 : }
74 :
75 : Real
76 4396336 : 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 4396336 : std::pow(10.0, _log_k[_qp]) * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u) / _gamma_eq[_qp];
83 4396336 : Real val_old = std::pow(10.0, _log_k[_qp]) * std::pow(_gamma_u_old[_qp] * _u_old[_qp], _sto_u) /
84 4396336 : _gamma_eq_old[_qp];
85 :
86 : // Contribution due to coupled primary species
87 7105552 : for (unsigned int i = 0; i < _vars.size(); ++i)
88 : {
89 2709216 : val_new *= std::pow((*_gamma_v[i])[_qp] * (*_v_vals[i])[_qp], _sto_v[i]);
90 2709216 : val_old *= std::pow((*_gamma_v_old[i])[_qp] * (*_v_vals_old[i])[_qp], _sto_v[i]);
91 : }
92 :
93 4396336 : return _porosity[_qp] * _weight * _test[_i][_qp] * (val_new - val_old) / _dt;
94 : }
95 :
96 : Real
97 2829440 : CoupledBEEquilibriumSub::computeQpJacobian()
98 : {
99 2829440 : Real val_new = std::pow(10.0, _log_k[_qp]) * _sto_u * _gamma_u[_qp] *
100 2829440 : std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _phi[_j][_qp] / _gamma_eq[_qp];
101 :
102 4585088 : for (unsigned int i = 0; i < _vars.size(); ++i)
103 1755648 : val_new *= std::pow((*_gamma_v[i])[_qp] * (*_v_vals[i])[_qp], _sto_v[i]);
104 :
105 2829440 : return _porosity[_qp] * _test[_i][_qp] * _weight * val_new / _dt;
106 : }
107 :
108 : Real
109 3698176 : CoupledBEEquilibriumSub::computeQpOffDiagJacobian(unsigned int jvar)
110 : {
111 : // If no coupled species, return 0
112 3698176 : if (_v_vals.size() == 0)
113 : return 0.0;
114 :
115 : // If jvar is not one of the coupled species, return 0
116 2600448 : if (std::find(_vars.begin(), _vars.end(), jvar) == _vars.end())
117 : return 0.0;
118 :
119 : Real val_new =
120 1755648 : std::pow(10.0, _log_k[_qp]) * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u) / _gamma_eq[_qp];
121 :
122 3644928 : for (unsigned int i = 0; i < _vars.size(); ++i)
123 : {
124 1889280 : if (jvar == _vars[i])
125 1755648 : val_new *= _sto_v[i] * (*_gamma_v[i])[_qp] *
126 1755648 : std::pow((*_gamma_v[i])[_qp] * (*_v_vals[i])[_qp], _sto_v[i] - 1.0) *
127 1755648 : _phi[_j][_qp];
128 : else
129 133632 : val_new *= std::pow((*_gamma_v[i])[_qp] * (*_v_vals[i])[_qp], _sto_v[i]);
130 : }
131 :
132 1755648 : return _porosity[_qp] * _test[_i][_qp] * _weight * val_new / _dt;
133 : }
|
