www.mooseframework.org
INSBase.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 /****************************************************************/
7 
8 #include "INSBase.h"
9 #include "Assembly.h"
10 
11 template <>
12 InputParameters
14 {
15  InputParameters params = validParams<Kernel>();
16 
17  params.addClassDescription("This class computes various strong and weak components of the "
18  "incompressible navier stokes equations which can then be assembled "
19  "together in child classes.");
20  // Coupled variables
21  params.addRequiredCoupledVar("u", "x-velocity");
22  params.addCoupledVar("v", 0, "y-velocity"); // only required in 2D and 3D
23  params.addCoupledVar("w", 0, "z-velocity"); // only required in 3D
24  params.addRequiredCoupledVar("p", "pressure");
25 
26  params.addParam<RealVectorValue>(
27  "gravity", RealVectorValue(0, 0, 0), "Direction of the gravity vector");
28 
29  params.addParam<MaterialPropertyName>("mu_name", "mu", "The name of the dynamic viscosity");
30  params.addParam<MaterialPropertyName>("rho_name", "rho", "The name of the density");
31 
32  params.addParam<Real>("alpha", 1., "Multiplicative factor on the stabilization parameter tau.");
33  params.addParam<bool>(
34  "laplace", true, "Whether the viscous term of the momentum equations is in laplace form.");
35  params.addParam<bool>("convective_term", true, "Whether to include the convective term.");
36  params.addParam<bool>("transient_term",
37  false,
38  "Whether there should be a transient term in the momentum residuals.");
39 
40  return params;
41 }
42 
43 INSBase::INSBase(const InputParameters & parameters)
44  : Kernel(parameters),
45  _second_phi(_assembly.secondPhi()),
46 
47  // Coupled variables
48  _u_vel(coupledValue("u")),
49  _v_vel(coupledValue("v")),
50  _w_vel(coupledValue("w")),
51  _p(coupledValue("p")),
52 
53  // Gradients
54  _grad_u_vel(coupledGradient("u")),
55  _grad_v_vel(coupledGradient("v")),
56  _grad_w_vel(coupledGradient("w")),
57  _grad_p(coupledGradient("p")),
58 
59  // second derivative tensors
60  _second_u_vel(coupledSecond("u")),
61  _second_v_vel(coupledSecond("v")),
62  _second_w_vel(coupledSecond("w")),
63 
64  // time derivatives
65  _u_vel_dot(coupledDot("u")),
66  _v_vel_dot(coupledDot("v")),
67  _w_vel_dot(coupledDot("w")),
68 
69  // derivatives of time derivatives
70  _d_u_vel_dot_du(coupledDotDu("u")),
71  _d_v_vel_dot_dv(coupledDotDu("v")),
72  _d_w_vel_dot_dw(coupledDotDu("w")),
73 
74  // Variable numberings
75  _u_vel_var_number(coupled("u")),
76  _v_vel_var_number(coupled("v")),
77  _w_vel_var_number(coupled("w")),
78  _p_var_number(coupled("p")),
79 
80  _gravity(getParam<RealVectorValue>("gravity")),
81 
82  // Material properties
83  _mu(getMaterialProperty<Real>("mu_name")),
84  _rho(getMaterialProperty<Real>("rho_name")),
85 
86  _alpha(getParam<Real>("alpha")),
87  _laplace(getParam<bool>("laplace")),
88  _convective_term(getParam<bool>("convective_term")),
89  _transient_term(getParam<bool>("transient_term"))
90 {
91 }
92 
93 RealVectorValue
95 {
96  RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
97  return _rho[_qp] *
98  RealVectorValue(U * _grad_u_vel[_qp], U * _grad_v_vel[_qp], U * _grad_w_vel[_qp]);
99 }
100 
101 RealVectorValue
103 {
104  RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
105  RealVectorValue d_U_d_comp(0, 0, 0);
106  d_U_d_comp(comp) = _phi[_j][_qp];
107 
108  RealVectorValue convective_term = _rho[_qp] * RealVectorValue(d_U_d_comp * _grad_u_vel[_qp],
109  d_U_d_comp * _grad_v_vel[_qp],
110  d_U_d_comp * _grad_w_vel[_qp]);
111  convective_term(comp) += _rho[_qp] * U * _grad_phi[_j][_qp];
112 
113  return convective_term;
114 }
115 
116 RealVectorValue
118 {
119  return -_mu[_qp] *
120  RealVectorValue(_second_u_vel[_qp].tr(), _second_v_vel[_qp].tr(), _second_w_vel[_qp].tr());
121 }
122 
123 RealVectorValue
125 {
126  return strongViscousTermLaplace() -
127  _mu[_qp] *
128  (_second_u_vel[_qp].row(0) + _second_v_vel[_qp].row(1) + _second_w_vel[_qp].row(2));
129 }
130 
131 RealVectorValue
133 {
134  RealVectorValue viscous_term(0, 0, 0);
135  viscous_term(comp) = -_mu[_qp] * _second_phi[_j][_qp].tr();
136 
137  return viscous_term;
138 }
139 
140 RealVectorValue
142 {
143  RealVectorValue viscous_term(0, 0, 0);
144  viscous_term(comp) = -_mu[_qp] * (_second_phi[_j][_qp](0, 0) + _second_phi[_j][_qp](1, 1) +
145  _second_phi[_j][_qp](2, 2));
146  for (unsigned i = 0; i < 3; i++)
147  viscous_term(i) += -_mu[_qp] * _second_phi[_j][_qp](i, comp);
148 
149  return viscous_term;
150 }
151 
152 RealVectorValue
154 {
155  switch (comp)
156  {
157  case 0:
158  return _mu[_qp] * _grad_u_vel[_qp];
159 
160  case 1:
161  return _mu[_qp] * _grad_v_vel[_qp];
162 
163  case 2:
164  return _mu[_qp] * _grad_w_vel[_qp];
165 
166  default:
167  return _zero[_qp];
168  }
169 }
170 
171 RealVectorValue
173 {
174  switch (comp)
175  {
176  case 0:
177  {
178  RealVectorValue transpose(_grad_u_vel[_qp](0), _grad_v_vel[_qp](0), _grad_w_vel[_qp](0));
179  return _mu[_qp] * _grad_u_vel[_qp] + _mu[_qp] * transpose;
180  }
181 
182  case 1:
183  {
184  RealVectorValue transpose(_grad_u_vel[_qp](1), _grad_v_vel[_qp](1), _grad_w_vel[_qp](1));
185  return _mu[_qp] * _grad_v_vel[_qp] + _mu[_qp] * transpose;
186  }
187 
188  case 2:
189  {
190  RealVectorValue transpose(_grad_u_vel[_qp](2), _grad_v_vel[_qp](2), _grad_w_vel[_qp](2));
191  return _mu[_qp] * _grad_w_vel[_qp] + _mu[_qp] * transpose;
192  }
193 
194  default:
195  return _zero[_qp];
196  }
197 }
198 
199 RealVectorValue
201 {
202  return _mu[_qp] * _grad_phi[_j][_qp];
203 }
204 
205 RealVectorValue
207 {
208  return _mu[_qp] * _grad_phi[_j][_qp];
209 }
210 
211 RealVectorValue
213 {
214  return _grad_p[_qp];
215 }
216 
217 Real
219 {
220  return -_p[_qp];
221 }
222 
223 RealVectorValue
225 {
226  return _grad_phi[_j][_qp];
227 }
228 
229 Real
231 {
232  return -_phi[_j][_qp];
233 }
234 
235 RealVectorValue
237 {
238  return -_rho[_qp] * _gravity;
239 }
240 
241 RealVectorValue
243 {
244  return _rho[_qp] * RealVectorValue(_u_vel_dot[_qp], _v_vel_dot[_qp], _w_vel_dot[_qp]);
245 }
246 
247 RealVectorValue
249 {
250  Real base = _rho[_qp] * _phi[_j][_qp];
251  switch (comp)
252  {
253  case 0:
254  return RealVectorValue(base * _d_u_vel_dot_du[_qp], 0, 0);
255 
256  case 1:
257  return RealVectorValue(0, base * _d_v_vel_dot_dv[_qp], 0);
258 
259  case 2:
260  return RealVectorValue(0, 0, base * _d_w_vel_dot_dw[_qp]);
261 
262  default:
263  mooseError("comp must be 0, 1, or 2");
264  }
265 }
266 
267 Real
269 {
270  Real nu = _mu[_qp] / _rho[_qp];
271  RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
272  Real h = _current_elem->hmax();
273  Real transient_part = _transient_term ? 4. / (_dt * _dt) : 0.;
274  return _alpha / std::sqrt(transient_part + (2. * U.norm() / h) * (2. * U.norm() / h) +
275  9. * (4. * nu / (h * h)) * (4. * nu / (h * h)));
276 }
277 
278 Real
280 {
281  Real nu = _mu[_qp] / _rho[_qp];
282  RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
283  Real h = _current_elem->hmax();
284  Real xi;
285  if (nu < std::numeric_limits<Real>::epsilon())
286  xi = 1;
287  else
288  {
289  Real alpha = U.norm() * h / (2 * nu);
290  xi = 1. / std::tanh(alpha) - 1. / alpha;
291  }
292  return h / (2 * U.norm()) * xi;
293 }
294 
295 Real
296 INSBase::dTauDUComp(unsigned comp)
297 {
298  Real nu = _mu[_qp] / _rho[_qp];
299  RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
300  Real h = _current_elem->hmax();
301  Real transient_part = _transient_term ? 4. / (_dt * _dt) : 0.;
302  return -_alpha / 2. * std::pow(transient_part + (2. * U.norm() / h) * (2. * U.norm() / h) +
303  9. * (4. * nu / (h * h)) * (4. * nu / (h * h)),
304  -1.5) *
305  2. * (2. * U.norm() / h) * 2. / h * U(comp) * _phi[_j][_qp] /
306  (U.norm() + std::numeric_limits<double>::epsilon());
307 }
virtual RealVectorValue strongViscousTermLaplace()
Definition: INSBase.C:117
const VariableValue & _p
Definition: INSBase.h:70
const VariablePhiSecond & _second_phi
second derivatives of the shape function
Definition: INSBase.h:64
virtual RealVectorValue dStrongViscDUCompTraction(unsigned comp)
Definition: INSBase.C:141
virtual Real tau()
Definition: INSBase.C:268
const MaterialProperty< Real > & _rho
Definition: INSBase.h:103
virtual Real tauNodal()
Provides tau which yields superconvergence for 1D advection-diffusion.
Definition: INSBase.C:279
const VariableValue & _u_vel
Definition: INSBase.h:67
InputParameters validParams< INSBase >()
Definition: INSBase.C:13
const MaterialProperty< Real > & _mu
Definition: INSBase.h:102
virtual RealVectorValue dWeakViscDUCompTraction()
Definition: INSBase.C:206
virtual Real dTauDUComp(unsigned comp)
Definition: INSBase.C:296
const VariableGradient & _grad_v_vel
Definition: INSBase.h:74
const VariableSecond & _second_u_vel
Definition: INSBase.h:79
const VariableValue & _w_vel_dot
Definition: INSBase.h:86
const VariableValue & _w_vel
Definition: INSBase.h:69
virtual RealVectorValue dStrongPressureDPressure()
Definition: INSBase.C:224
virtual RealVectorValue timeDerivativeTerm()
Definition: INSBase.C:242
virtual RealVectorValue convectiveTerm()
Definition: INSBase.C:94
const VariableValue & _d_v_vel_dot_dv
Definition: INSBase.h:90
virtual RealVectorValue weakViscousTermTraction(unsigned comp)
Definition: INSBase.C:172
const VariableGradient & _grad_u_vel
Definition: INSBase.h:73
const VariableSecond & _second_v_vel
Definition: INSBase.h:80
virtual RealVectorValue strongPressureTerm()
Definition: INSBase.C:212
virtual Real dWeakPressureDPressure()
Definition: INSBase.C:230
bool _transient_term
Definition: INSBase.h:108
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
const VariableValue & _u_vel_dot
Definition: INSBase.h:84
virtual RealVectorValue dConvecDUComp(unsigned comp)
Definition: INSBase.C:102
virtual RealVectorValue dStrongViscDUCompLaplace(unsigned comp)
Definition: INSBase.C:132
const VariableValue & _d_w_vel_dot_dw
Definition: INSBase.h:91
const VariableGradient & _grad_w_vel
Definition: INSBase.h:75
const VariableSecond & _second_w_vel
Definition: INSBase.h:81
virtual RealVectorValue gravityTerm()
Definition: INSBase.C:236
INSBase(const InputParameters &parameters)
Definition: INSBase.C:43
virtual RealVectorValue strongViscousTermTraction()
Definition: INSBase.C:124
const VariableGradient & _grad_p
Definition: INSBase.h:76
virtual RealVectorValue dWeakViscDUCompLaplace()
Definition: INSBase.C:200
const VariableValue & _v_vel
Definition: INSBase.h:68
virtual RealVectorValue weakViscousTermLaplace(unsigned comp)
Definition: INSBase.C:153
RealVectorValue _gravity
Definition: INSBase.h:99
virtual Real weakPressureTerm()
Definition: INSBase.C:218
virtual RealVectorValue dTimeDerivativeDUComp(unsigned comp)
Definition: INSBase.C:248
const VariableValue & _d_u_vel_dot_du
Definition: INSBase.h:89
const VariableValue & _v_vel_dot
Definition: INSBase.h:85
const Real & _alpha
Definition: INSBase.h:105