www.mooseframework.org
INSSplitMomentum.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 "INSSplitMomentum.h"
9 #include "MooseMesh.h"
10 
11 template <>
12 InputParameters
14 {
15  InputParameters params = validParams<Kernel>();
16 
17  params.addClassDescription("This class computes the 'split' momentum equation residual.");
18  // Coupled variables
19  params.addRequiredCoupledVar("u", "x-velocity");
20  params.addCoupledVar("v", "y-velocity"); // only required in 2D and 3D
21  params.addCoupledVar("w", "z-velocity"); // only required in 3D
22 
23  params.addRequiredCoupledVar("a1", "x-acceleration");
24  params.addCoupledVar("a2", "y-acceleration"); // only required in 2D and 3D
25  params.addCoupledVar("a3", "z-acceleration"); // only required in 3D
26 
27  // Required parameters
28  params.addRequiredParam<RealVectorValue>("gravity", "Direction of the gravity vector");
29  params.addRequiredParam<unsigned>(
30  "component",
31  "0,1,2 depending on if we are solving the x,y,z component of the momentum equation");
32 
33  // Optional parameters
34  params.addParam<MaterialPropertyName>("mu_name", "mu", "The name of the dynamic viscosity");
35  params.addParam<MaterialPropertyName>("rho_name", "rho", "The name of the density");
36 
37  return params;
38 }
39 
40 INSSplitMomentum::INSSplitMomentum(const InputParameters & parameters)
41  : Kernel(parameters),
42 
43  // Coupled variables
44  _u_vel(coupledValue("u")),
45  _v_vel(_mesh.dimension() >= 2 ? coupledValue("v") : _zero),
46  _w_vel(_mesh.dimension() == 3 ? coupledValue("w") : _zero),
47 
48  _a1(coupledValue("a1")),
49  _a2(_mesh.dimension() >= 2 ? coupledValue("a2") : _zero),
50  _a3(_mesh.dimension() == 3 ? coupledValue("a3") : _zero),
51 
52  // Gradients
53  _grad_u_vel(coupledGradient("u")),
54  _grad_v_vel(_mesh.dimension() >= 2 ? coupledGradient("v") : _grad_zero),
55  _grad_w_vel(_mesh.dimension() == 3 ? coupledGradient("w") : _grad_zero),
56 
57  // Variable numberings
58  _u_vel_var_number(coupled("u")),
59  _v_vel_var_number(_mesh.dimension() >= 2 ? coupled("v") : libMesh::invalid_uint),
60  _w_vel_var_number(_mesh.dimension() == 3 ? coupled("w") : libMesh::invalid_uint),
61 
62  _a1_var_number(coupled("a1")),
63  _a2_var_number(_mesh.dimension() >= 2 ? coupled("a2") : libMesh::invalid_uint),
64  _a3_var_number(_mesh.dimension() == 3 ? coupled("a3") : libMesh::invalid_uint),
65 
66  // Required parameters
67  _gravity(getParam<RealVectorValue>("gravity")),
68  _component(getParam<unsigned>("component")),
69 
70  // Material properties
71  _mu(getMaterialProperty<Real>("mu_name")),
72  _rho(getMaterialProperty<Real>("rho_name"))
73 {
74 }
75 
76 Real
78 {
79  // Vector object for a
80  RealVectorValue a(_a1[_qp], _a2[_qp], _a3[_qp]);
81 
82  // Vector object for test function
83  RealVectorValue test;
84  test(_component) = _test[_i][_qp];
85 
86  // Tensor object for "grad U" = d(u_i)/d(x_j)
87  RealTensorValue grad_U(_grad_u_vel[_qp], _grad_v_vel[_qp], _grad_w_vel[_qp]);
88 
89  // Vector object for U
90  RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
91 
92  // Viscous stress tensor object
93  RealTensorValue tau(
94  // Row 1
95  2. * _grad_u_vel[_qp](0),
96  _grad_u_vel[_qp](1) + _grad_v_vel[_qp](0),
97  _grad_u_vel[_qp](2) + _grad_w_vel[_qp](0),
98  // Row 2
99  _grad_v_vel[_qp](0) + _grad_u_vel[_qp](1),
100  2. * _grad_v_vel[_qp](1),
101  _grad_v_vel[_qp](2) + _grad_w_vel[_qp](1),
102  // Row 3
103  _grad_w_vel[_qp](0) + _grad_u_vel[_qp](2),
104  _grad_w_vel[_qp](1) + _grad_v_vel[_qp](2),
105  2. * _grad_w_vel[_qp](2));
106 
107  // Tensor object for test function gradient
108  RealTensorValue grad_test;
109  for (unsigned k = 0; k < 3; ++k)
110  grad_test(_component, k) = _grad_test[_i][_qp](k);
111 
112  //
113  // Compute the different pieces...
114  //
115 
116  // "Symmetric" part, a.test
117  Real symmetric_part = a(_component) * _test[_i][_qp];
118 
119  // The convection part, (u.grad) * u_component * test_scalar. Which can also be
120  // written as (grad_U * U) * test_vec
121  Real convective_part = (grad_U * U) * test;
122 
123  // The viscous part, tau : grad(v)
124  Real viscous_part = (_mu[_qp] / _rho[_qp]) * tau.contract(grad_test);
125 
126  return symmetric_part + convective_part + viscous_part;
127 }
128 
129 Real
131 {
132  // The on-diagonal Jacobian contribution (derivative of a.test wrt a)
133  // is just the mass matrix entry.
134  return _phi[_j][_qp] * _test[_i][_qp];
135 }
136 
137 Real
139 {
140  if ((jvar == _u_vel_var_number) || (jvar == _v_vel_var_number) || (jvar == _w_vel_var_number))
141  {
142  // Derivative of viscous stress tensor
143  RealTensorValue dtau;
144 
145  // Initialize to invalid value, then determine correct value.
146  unsigned vel_index = 99;
147 
148  // Set index and build dtau for that index
149  if (jvar == _u_vel_var_number)
150  {
151  vel_index = 0;
152  dtau(0, 0) = 2. * _grad_phi[_j][_qp](0);
153  dtau(0, 1) = _grad_phi[_j][_qp](1);
154  dtau(0, 2) = _grad_phi[_j][_qp](2);
155  dtau(1, 0) = _grad_phi[_j][_qp](1);
156  dtau(2, 0) = _grad_phi[_j][_qp](2);
157  }
158  else if (jvar == _v_vel_var_number)
159  {
160  vel_index = 1;
161  /* */ dtau(0, 1) = _grad_phi[_j][_qp](0);
162  dtau(1, 0) = _grad_phi[_j][_qp](0);
163  dtau(1, 1) = 2. * _grad_phi[_j][_qp](1);
164  dtau(1, 2) = _grad_phi[_j][_qp](2);
165  /* */ dtau(2, 1) = _grad_phi[_j][_qp](2);
166  }
167  else if (jvar == _w_vel_var_number)
168  {
169  vel_index = 2;
170  /* */ dtau(0, 2) =
171  _grad_phi[_j][_qp](0);
172  /* */ dtau(1, 2) =
173  _grad_phi[_j][_qp](1);
174  dtau(2, 0) = _grad_phi[_j][_qp](0);
175  dtau(2, 1) = _grad_phi[_j][_qp](1);
176  dtau(2, 2) = 2. * _grad_phi[_j][_qp](2);
177  }
178 
179  // Vector object for test function
180  RealVectorValue test;
181  test(_component) = _test[_i][_qp];
182 
183  // Vector object for U
184  RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
185 
186  // Tensor object for test function gradient
187  RealTensorValue grad_test;
188  for (unsigned k = 0; k < 3; ++k)
189  grad_test(_component, k) = _grad_test[_i][_qp](k);
190 
191  // Compute the convective part
192  RealVectorValue convective_jac = _phi[_j][_qp] * RealVectorValue(_grad_u_vel[_qp](vel_index),
193  _grad_v_vel[_qp](vel_index),
194  _grad_w_vel[_qp](vel_index));
195 
196  // Extra contribution in vel_index component
197  convective_jac(vel_index) += U * _grad_phi[_j][_qp];
198  Real convective_part = convective_jac * test;
199 
200  // Compute the viscous part
201  Real viscous_part = (_mu[_qp] / _rho[_qp]) * dtau.contract(grad_test);
202 
203  // Return the result
204  return convective_part + viscous_part;
205  }
206 
207  else
208  return 0;
209 }
InputParameters validParams< INSSplitMomentum >()
const VariableValue & _a2
unsigned _u_vel_var_number
const MaterialProperty< Real > & _mu
virtual Real computeQpResidual()
const VariableValue & _a1
const MaterialProperty< Real > & _rho
virtual Real computeQpOffDiagJacobian(unsigned jvar)
const VariableGradient & _grad_u_vel
virtual Real computeQpJacobian()
const VariableGradient & _grad_v_vel
const VariableValue & _u_vel
const VariableValue & _w_vel
const VariableGradient & _grad_w_vel
unsigned _v_vel_var_number
unsigned _w_vel_var_number
const VariableValue & _v_vel
const VariableValue & _a3
INSSplitMomentum(const InputParameters &parameters)