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