www.mooseframework.org
CoupledConvectionReactionSub.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
8 
9 template <>
10 InputParameters
12 {
13  InputParameters params = validParams<Kernel>();
14  params.addParam<Real>("weight", 1.0, "Weight of the equilibrium species");
15  params.addParam<Real>("log_k", 0.0, "Equilibrium constant of dissociation equilibrium reaction");
16  params.addParam<Real>("sto_u",
17  1.0,
18  "Stoichiometric coef of the primary species the kernel "
19  "operates on in the equilibrium reaction");
20  params.addRequiredParam<std::vector<Real>>(
21  "sto_v",
22  "The stoichiometric coefficients of coupled primary species in equilibrium reaction");
23  params.addRequiredCoupledVar("p", "Pressure");
24  params.addCoupledVar("v", "List of coupled primary species");
25  RealVectorValue g(0, 0, 0);
26  params.addParam<RealVectorValue>("gravity", g, "Gravity vector (default is (0, 0, 0))");
27  params.addClassDescription("Convection of equilibrium species");
28  return params;
29 }
30 
32  : DerivativeMaterialInterface<Kernel>(parameters),
33  _weight(getParam<Real>("weight")),
34  _log_k(getParam<Real>("log_k")),
35  _sto_u(getParam<Real>("sto_u")),
36  _sto_v(getParam<std::vector<Real>>("sto_v")),
37  _cond(getMaterialProperty<Real>("conductivity")),
38  _gravity(getParam<RealVectorValue>("gravity")),
39  _density(getDefaultMaterialProperty<Real>("density")),
40  _grad_p(coupledGradient("p")),
41  _pvar(coupled("p"))
42 {
43  const unsigned int n = coupledComponents("v");
44  _vars.resize(n);
45  _vals.resize(n);
46  _grad_vals.resize(n);
47 
48  for (unsigned int i = 0; i < n; ++i)
49  {
50  _vars[i] = coupled("v", i);
51  _vals[i] = &coupledValue("v", i);
52  _grad_vals[i] = &coupledGradient("v", i);
53  }
54 }
55 
56 Real
58 {
59  RealVectorValue darcy_vel = -_cond[_qp] * (_grad_p[_qp] - _density[_qp] * _gravity);
60  RealGradient d_u = _sto_u * std::pow(_u[_qp], _sto_u - 1.0) * _grad_u[_qp];
61  RealGradient d_var_sum(0.0, 0.0, 0.0);
62  const Real d_v_u = std::pow(_u[_qp], _sto_u);
63 
64  for (unsigned int i = 0; i < _vals.size(); ++i)
65  {
66  d_u *= std::pow((*_vals[i])[_qp], _sto_v[i]);
67 
68  RealGradient d_var =
69  d_v_u * _sto_v[i] * std::pow((*_vals[i])[_qp], _sto_v[i] - 1.0) * (*_grad_vals[i])[_qp];
70  for (unsigned int j = 0; j < _vals.size(); ++j)
71  if (j != i)
72  d_var *= std::pow((*_vals[j])[_qp], _sto_v[j]);
73 
74  d_var_sum += d_var;
75  }
76 
77  return _weight * std::pow(10.0, _log_k) * _test[_i][_qp] * darcy_vel * (d_u + d_var_sum);
78 }
79 
80 Real
82 {
83  RealVectorValue darcy_vel = -_cond[_qp] * (_grad_p[_qp] - _density[_qp] * _gravity);
84 
85  RealGradient d_u_1 = _sto_u * std::pow(_u[_qp], _sto_u - 1.0) * _grad_phi[_j][_qp];
86  RealGradient d_u_2 =
87  _phi[_j][_qp] * _sto_u * (_sto_u - 1.0) * std::pow(_u[_qp], _sto_u - 2.0) * _grad_u[_qp];
88 
89  RealGradient d_var_sum(0.0, 0.0, 0.0);
90  const Real d_v_u = _sto_u * std::pow(_u[_qp], _sto_u - 1.0) * _phi[_j][_qp];
91 
92  for (unsigned int i = 0; i < _vals.size(); ++i)
93  {
94  d_u_1 *= std::pow((*_vals[i])[_qp], _sto_v[i]);
95  d_u_2 *= std::pow((*_vals[i])[_qp], _sto_v[i]);
96 
97  RealGradient d_var =
98  d_v_u * _sto_v[i] * std::pow((*_vals[i])[_qp], _sto_v[i] - 1.0) * (*_grad_vals[i])[_qp];
99  for (unsigned int j = 0; j < _vals.size(); ++j)
100  if (j != i)
101  d_var *= std::pow((*_vals[j])[_qp], _sto_v[j]);
102 
103  d_var_sum += d_var;
104  }
105 
106  RealGradient d_u_j = d_u_1 + d_u_2;
107  return _weight * std::pow(10.0, _log_k) * _test[_i][_qp] * darcy_vel * (d_u_j + d_var_sum);
108 }
109 
110 Real
112 {
113  if (jvar == _pvar)
114  {
115  RealVectorValue ddarcy_vel_dp = -_cond[_qp] * _grad_phi[_j][_qp];
116 
117  RealGradient d_u = _sto_u * std::pow(_u[_qp], _sto_u - 1.0) * _grad_u[_qp];
118  RealGradient d_var_sum(0.0, 0.0, 0.0);
119  const Real d_v_u = std::pow(_u[_qp], _sto_u);
120 
121  for (unsigned int i = 0; i < _vals.size(); ++i)
122  {
123  d_u *= std::pow((*_vals[i])[_qp], _sto_v[i]);
124 
125  RealGradient d_var =
126  d_v_u * _sto_v[i] * std::pow((*_vals[i])[_qp], _sto_v[i] - 1.0) * (*_grad_vals[i])[_qp];
127  for (unsigned int j = 0; j < _vals.size(); ++j)
128  if (j != i)
129  d_var *= std::pow((*_vals[j])[_qp], _sto_v[j]);
130 
131  d_var_sum += d_var;
132  }
133  return _weight * std::pow(10.0, _log_k) * _test[_i][_qp] * ddarcy_vel_dp * (d_u + d_var_sum);
134  }
135 
136  if (_vals.size() == 0)
137  return 0.0;
138 
139  RealVectorValue darcy_vel = -_cond[_qp] * (_grad_p[_qp] - _density[_qp] * _gravity);
140  RealGradient diff1 = _sto_u * std::pow(_u[_qp], _sto_u - 1.0) * _grad_u[_qp];
141  for (unsigned int i = 0; i < _vals.size(); ++i)
142  {
143  if (jvar == _vars[i])
144  diff1 *= _sto_v[i] * std::pow((*_vals[i])[_qp], _sto_v[i] - 1.0) * _phi[_j][_qp];
145  else
146  diff1 *= std::pow((*_vals[i])[_qp], _sto_v[i]);
147  }
148 
149  Real val_u = std::pow(_u[_qp], _sto_u);
150  RealGradient diff2_1(1.0, 1.0, 1.0);
151  RealGradient diff2_2(1.0, 1.0, 1.0);
152  for (unsigned int i = 0; i < _vals.size(); ++i)
153  if (jvar == _vars[i])
154  {
155  diff2_1 = _sto_v[i] * (_sto_v[i] - 1.0) * std::pow((*_vals[i])[_qp], _sto_v[i] - 2.0) *
156  _phi[_j][_qp] * (*_grad_vals[i])[_qp];
157  diff2_2 = _sto_v[i] * std::pow((*_vals[i])[_qp], _sto_v[i] - 1.0) * _grad_phi[_j][_qp];
158  }
159 
160  RealGradient diff2 = val_u * (diff2_1 + diff2_2);
161  for (unsigned int i = 0; i < _vals.size(); ++i)
162  if (jvar != _vars[i])
163  {
164  diff2 *= std::pow((*_vals[i])[_qp], _sto_v[i]);
165  }
166 
167  RealGradient diff3;
168  RealGradient diff3_sum(0.0, 0.0, 0.0);
169  Real val_jvar = 0.0;
170  unsigned int var = 0;
171 
172  for (unsigned int i = 0; i < _vals.size(); ++i)
173  if (jvar == _vars[i])
174  {
175  var = i;
176  val_jvar = val_u * _sto_v[i] * std::pow((*_vals[i])[_qp], _sto_v[i] - 1.0) * _phi[_j][_qp];
177  }
178 
179  for (unsigned int i = 0; i < _vals.size(); ++i)
180  if (i != var)
181  {
182  diff3 = val_jvar * _sto_v[i] * std::pow((*_vals[i])[_qp], _sto_v[i] - 1.0) *
183  (*_grad_vals[i])[_qp];
184 
185  for (unsigned int j = 0; j < _vals.size(); ++j)
186  if (j != var && j != i)
187  diff3 *= std::pow((*_vals[j])[_qp], _sto_v[j]);
188 
189  diff3_sum += diff3;
190  }
191 
192  return _weight * std::pow(10.0, _log_k) * _test[_i][_qp] * darcy_vel *
193  (diff1 + diff2 + diff3_sum);
194 }
InputParameters validParams< CoupledConvectionReactionSub >()
std::vector< const VariableGradient * > _grad_vals
Coupled gradients of primary species concentrations.
const MaterialProperty< Real > & _density
Fluid density.
const Real _log_k
Equilibrium constant for the equilibrium species in association form.
const std::vector< Real > _sto_v
Stoichiometric coefficients of the coupled primary species.
const VariableGradient & _grad_p
Pressure gradient.
virtual Real computeQpOffDiagJacobian(unsigned int jvar) override
const unsigned int _pvar
Pressure variable number.
const Real _weight
Weight of the equilibrium species concentration in the total primary species concentration.
const RealVectorValue _gravity
Gravity.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
std::vector< unsigned int > _vars
Coupled primary species variable numbers.
const Real _sto_u
Stoichiometric coefficient of the primary species.
const MaterialProperty< Real > & _cond
Hydraulic conductivity.
std::vector< const VariableValue * > _vals
Coupled primary species concentrations.
CoupledConvectionReactionSub(const InputParameters &parameters)