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 "INSBase.h"
11 : #include "Assembly.h"
12 : #include "NS.h"
13 :
14 : InputParameters
15 3248 : INSBase::validParams()
16 : {
17 3248 : InputParameters params = Kernel::validParams();
18 :
19 3248 : params.addClassDescription("This class computes various strong and weak components of the "
20 : "incompressible navier stokes equations which can then be assembled "
21 : "together in child classes.");
22 : // Coupled variables
23 6496 : params.addRequiredCoupledVar("u", "x-velocity");
24 6496 : params.addCoupledVar("v", 0, "y-velocity"); // only required in 2D and 3D
25 6496 : params.addCoupledVar("w", 0, "z-velocity"); // only required in 3D
26 3248 : params.addRequiredCoupledVar(NS::pressure, "pressure");
27 :
28 3248 : params.addParam<RealVectorValue>(
29 3248 : "gravity", RealVectorValue(0, 0, 0), "Direction of the gravity vector");
30 :
31 6496 : params.addParam<MaterialPropertyName>("mu_name", "mu", "The name of the dynamic viscosity");
32 6496 : params.addParam<MaterialPropertyName>("rho_name", "rho", "The name of the density");
33 :
34 6496 : params.addParam<Real>("alpha", 1., "Multiplicative factor on the stabilization parameter tau.");
35 6496 : params.addParam<bool>(
36 6496 : "laplace", true, "Whether the viscous term of the momentum equations is in laplace form.");
37 6496 : params.addParam<bool>("convective_term", true, "Whether to include the convective term.");
38 6496 : params.addParam<bool>("transient_term",
39 6496 : false,
40 : "Whether there should be a transient term in the momentum residuals.");
41 6496 : params.addCoupledVar("disp_x", "The x displacement");
42 6496 : params.addCoupledVar("disp_y", "The y displacement");
43 6496 : params.addCoupledVar("disp_z", "The z displacement");
44 6496 : params.addParam<bool>("picard",
45 6496 : false,
46 : "Whether we are applying a Picard strategy in which case we will linearize "
47 : "the nonlinear convective term.");
48 :
49 3248 : return params;
50 0 : }
51 :
52 1751 : INSBase::INSBase(const InputParameters & parameters)
53 : : Kernel(parameters),
54 1751 : _second_phi(_assembly.secondPhi()),
55 :
56 : // Coupled variables
57 1751 : _u_vel(coupledValue("u")),
58 1751 : _v_vel(coupledValue("v")),
59 1751 : _w_vel(coupledValue("w")),
60 1751 : _p(coupledValue(NS::pressure)),
61 :
62 3502 : _picard(getParam<bool>("picard")),
63 1751 : _u_vel_previous_nl(_picard ? &coupledValuePreviousNL("u") : nullptr),
64 1751 : _v_vel_previous_nl(_picard ? &coupledValuePreviousNL("v") : nullptr),
65 1751 : _w_vel_previous_nl(_picard ? &coupledValuePreviousNL("w") : nullptr),
66 :
67 : // Gradients
68 1751 : _grad_u_vel(coupledGradient("u")),
69 1751 : _grad_v_vel(coupledGradient("v")),
70 1751 : _grad_w_vel(coupledGradient("w")),
71 1751 : _grad_p(coupledGradient(NS::pressure)),
72 :
73 : // second derivative tensors
74 1751 : _second_u_vel(coupledSecond("u")),
75 1751 : _second_v_vel(coupledSecond("v")),
76 1751 : _second_w_vel(coupledSecond("w")),
77 :
78 : // time derivatives
79 1751 : _u_vel_dot(_is_transient ? coupledDot("u") : _zero),
80 1751 : _v_vel_dot(_is_transient ? coupledDot("v") : _zero),
81 1751 : _w_vel_dot(_is_transient ? coupledDot("w") : _zero),
82 :
83 : // derivatives of time derivatives
84 1751 : _d_u_vel_dot_du(_is_transient ? coupledDotDu("u") : _zero),
85 1751 : _d_v_vel_dot_dv(_is_transient ? coupledDotDu("v") : _zero),
86 1751 : _d_w_vel_dot_dw(_is_transient ? coupledDotDu("w") : _zero),
87 :
88 : // Variable numberings
89 1751 : _u_vel_var_number(coupled("u")),
90 1751 : _v_vel_var_number(coupled("v")),
91 1751 : _w_vel_var_number(coupled("w")),
92 1751 : _p_var_number(coupled(NS::pressure)),
93 :
94 3502 : _gravity(getParam<RealVectorValue>("gravity")),
95 :
96 : // Material properties
97 3502 : _mu(getMaterialProperty<Real>("mu_name")),
98 3502 : _rho(getMaterialProperty<Real>("rho_name")),
99 :
100 3502 : _alpha(getParam<Real>("alpha")),
101 3502 : _laplace(getParam<bool>("laplace")),
102 3502 : _convective_term(getParam<bool>("convective_term")),
103 3502 : _transient_term(getParam<bool>("transient_term")),
104 :
105 : // Displacements for mesh velocity for ALE simulations
106 3502 : _disps_provided(isParamValid("disp_x")),
107 3502 : _disp_x_dot(isParamValid("disp_x") ? coupledDot("disp_x") : _zero),
108 3502 : _disp_y_dot(isParamValid("disp_y") ? coupledDot("disp_y") : _zero),
109 3502 : _disp_z_dot(isParamValid("disp_z") ? coupledDot("disp_z") : _zero),
110 :
111 3502 : _rz_radial_coord(_mesh.getAxisymmetricRadialCoord())
112 : {
113 1751 : if (_picard && _disps_provided)
114 0 : paramError("picard",
115 : "Picard is not currently supported for ALE-type simulations in which we subtract "
116 : "the mesh velocity from the velocity variables");
117 1751 : }
118 :
119 : RealVectorValue
120 894082430 : INSBase::relativeVelocity() const
121 : {
122 894082430 : auto U = _picard ? RealVectorValue((*_u_vel_previous_nl)[_qp],
123 0 : (*_v_vel_previous_nl)[_qp],
124 0 : (*_w_vel_previous_nl)[_qp])
125 894082430 : : RealVectorValue(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
126 894082430 : if (_disps_provided)
127 0 : U -= RealVectorValue{_disp_x_dot[_qp], _disp_y_dot[_qp], _disp_z_dot[_qp]};
128 894082430 : return U;
129 : }
130 :
131 : RealVectorValue
132 346332822 : INSBase::convectiveTerm()
133 : {
134 346332822 : const auto U = _picard ? RealVectorValue((*_u_vel_previous_nl)[_qp],
135 4520448 : (*_v_vel_previous_nl)[_qp],
136 4520448 : (*_w_vel_previous_nl)[_qp])
137 346332822 : : RealVectorValue(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
138 346332822 : return _rho[_qp] *
139 346332822 : RealVectorValue(U * _grad_u_vel[_qp], U * _grad_v_vel[_qp], U * _grad_w_vel[_qp]);
140 : }
141 :
142 : RealVectorValue
143 1140163506 : INSBase::dConvecDUComp(unsigned comp)
144 : {
145 1140163506 : if (_picard)
146 : {
147 : RealVectorValue U(
148 79128576 : (*_u_vel_previous_nl)[_qp], (*_v_vel_previous_nl)[_qp], (*_w_vel_previous_nl)[_qp]);
149 : RealVectorValue convective_term;
150 79128576 : convective_term(comp) = _rho[_qp] * U * _grad_phi[_j][_qp];
151 :
152 79128576 : return convective_term;
153 : }
154 : else
155 : {
156 1061034930 : RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
157 : RealVectorValue d_U_d_comp(0, 0, 0);
158 1061034930 : d_U_d_comp(comp) = _phi[_j][_qp];
159 :
160 1061034930 : RealVectorValue convective_term = _rho[_qp] * RealVectorValue(d_U_d_comp * _grad_u_vel[_qp],
161 1061034930 : d_U_d_comp * _grad_v_vel[_qp],
162 1061034930 : d_U_d_comp * _grad_w_vel[_qp]);
163 1061034930 : convective_term(comp) += _rho[_qp] * U * _grad_phi[_j][_qp];
164 :
165 1061034930 : return convective_term;
166 : }
167 : }
168 :
169 : RealVectorValue
170 224944424 : INSBase::strongViscousTermLaplace()
171 : {
172 224944424 : return -_mu[_qp] *
173 224944424 : RealVectorValue(_second_u_vel[_qp].tr(), _second_v_vel[_qp].tr(), _second_w_vel[_qp].tr());
174 : }
175 :
176 : RealVectorValue
177 23842008 : INSBase::strongViscousTermTraction()
178 : {
179 23842008 : return INSBase::strongViscousTermLaplace() -
180 23842008 : _mu[_qp] *
181 71526024 : (_second_u_vel[_qp].row(0) + _second_v_vel[_qp].row(1) + _second_w_vel[_qp].row(2));
182 : }
183 :
184 : RealVectorValue
185 178286004 : INSBase::dStrongViscDUCompLaplace(unsigned comp)
186 : {
187 : RealVectorValue viscous_term(0, 0, 0);
188 178286004 : viscous_term(comp) = -_mu[_qp] * _second_phi[_j][_qp].tr();
189 :
190 178286004 : return viscous_term;
191 : }
192 :
193 : RealVectorValue
194 20151396 : INSBase::dStrongViscDUCompTraction(unsigned comp)
195 : {
196 : RealVectorValue viscous_term(0, 0, 0);
197 20151396 : viscous_term(comp) = -_mu[_qp] * (_second_phi[_j][_qp](0, 0) + _second_phi[_j][_qp](1, 1) +
198 20151396 : _second_phi[_j][_qp](2, 2));
199 80605584 : for (unsigned i = 0; i < 3; i++)
200 60454188 : viscous_term(i) += -_mu[_qp] * _second_phi[_j][_qp](i, comp);
201 :
202 20151396 : return viscous_term;
203 : }
204 :
205 : RealVectorValue
206 0 : INSBase::weakViscousTermLaplace(unsigned comp)
207 : {
208 0 : switch (comp)
209 : {
210 0 : case 0:
211 0 : return _mu[_qp] * _grad_u_vel[_qp];
212 :
213 0 : case 1:
214 0 : return _mu[_qp] * _grad_v_vel[_qp];
215 :
216 0 : case 2:
217 0 : return _mu[_qp] * _grad_w_vel[_qp];
218 :
219 0 : default:
220 0 : return _zero[_qp];
221 : }
222 : }
223 :
224 : RealVectorValue
225 0 : INSBase::weakViscousTermTraction(unsigned comp)
226 : {
227 0 : switch (comp)
228 : {
229 0 : case 0:
230 : {
231 0 : RealVectorValue transpose(_grad_u_vel[_qp](0), _grad_v_vel[_qp](0), _grad_w_vel[_qp](0));
232 0 : return _mu[_qp] * _grad_u_vel[_qp] + _mu[_qp] * transpose;
233 : }
234 :
235 0 : case 1:
236 : {
237 0 : RealVectorValue transpose(_grad_u_vel[_qp](1), _grad_v_vel[_qp](1), _grad_w_vel[_qp](1));
238 0 : return _mu[_qp] * _grad_v_vel[_qp] + _mu[_qp] * transpose;
239 : }
240 :
241 0 : case 2:
242 : {
243 0 : RealVectorValue transpose(_grad_u_vel[_qp](2), _grad_v_vel[_qp](2), _grad_w_vel[_qp](2));
244 0 : return _mu[_qp] * _grad_w_vel[_qp] + _mu[_qp] * transpose;
245 : }
246 :
247 0 : default:
248 0 : return _zero[_qp];
249 : }
250 : }
251 :
252 : RealVectorValue
253 0 : INSBase::dWeakViscDUCompLaplace()
254 : {
255 0 : return _mu[_qp] * _grad_phi[_j][_qp];
256 : }
257 :
258 : RealVectorValue
259 0 : INSBase::dWeakViscDUCompTraction()
260 : {
261 0 : return _mu[_qp] * _grad_phi[_j][_qp];
262 : }
263 :
264 : RealVectorValue
265 451771788 : INSBase::strongPressureTerm()
266 : {
267 451771788 : return _grad_p[_qp];
268 : }
269 :
270 : Real
271 59530728 : INSBase::weakPressureTerm()
272 : {
273 59530728 : return -_p[_qp];
274 : }
275 :
276 : RealVectorValue
277 120785137 : INSBase::dStrongPressureDPressure()
278 : {
279 120785137 : return _grad_phi[_j][_qp];
280 : }
281 :
282 : Real
283 175929725 : INSBase::dWeakPressureDPressure()
284 : {
285 175929725 : return -_phi[_j][_qp];
286 : }
287 :
288 : RealVectorValue
289 511302516 : INSBase::gravityTerm()
290 : {
291 511302516 : return -_rho[_qp] * _gravity;
292 : }
293 :
294 : RealVectorValue
295 10608768 : INSBase::timeDerivativeTerm()
296 : {
297 10608768 : return _rho[_qp] * RealVectorValue(_u_vel_dot[_qp], _v_vel_dot[_qp], _w_vel_dot[_qp]);
298 : }
299 :
300 : RealVectorValue
301 9790896 : INSBase::dTimeDerivativeDUComp(unsigned comp)
302 : {
303 9790896 : Real base = _rho[_qp] * _phi[_j][_qp];
304 9790896 : switch (comp)
305 : {
306 4895448 : case 0:
307 4895448 : return RealVectorValue(base * _d_u_vel_dot_du[_qp], 0, 0);
308 :
309 4895448 : case 1:
310 4895448 : return RealVectorValue(0, base * _d_v_vel_dot_dv[_qp], 0);
311 :
312 0 : case 2:
313 0 : return RealVectorValue(0, 0, base * _d_w_vel_dot_dw[_qp]);
314 :
315 0 : default:
316 0 : mooseError("comp must be 0, 1, or 2");
317 : }
318 : }
319 :
320 : Real
321 458570030 : INSBase::tau()
322 : {
323 458570030 : Real nu = _mu[_qp] / _rho[_qp];
324 458570030 : const auto U = relativeVelocity();
325 458570030 : Real h = _current_elem->hmax();
326 458570030 : Real transient_part = _transient_term ? 4. / (_dt * _dt) : 0.;
327 458570030 : return _alpha / std::sqrt(transient_part + (2. * U.norm() / h) * (2. * U.norm() / h) +
328 458570030 : 9. * (4. * nu / (h * h)) * (4. * nu / (h * h)));
329 : }
330 :
331 : Real
332 4200 : INSBase::tauNodal()
333 : {
334 4200 : Real nu = _mu[_qp] / _rho[_qp];
335 4200 : const auto U = relativeVelocity();
336 4200 : Real h = _current_elem->hmax();
337 : Real xi;
338 4200 : if (nu < std::numeric_limits<Real>::epsilon())
339 : xi = 1;
340 : else
341 : {
342 0 : Real alpha = U.norm() * h / (2 * nu);
343 0 : xi = 1. / std::tanh(alpha) - 1. / alpha;
344 : }
345 4200 : return h / (2 * U.norm()) * xi;
346 : }
347 :
348 : Real
349 198437400 : INSBase::dTauDUComp(const unsigned int comp)
350 : {
351 198437400 : Real nu = _mu[_qp] / _rho[_qp];
352 198437400 : const auto U = relativeVelocity();
353 198437400 : Real h = _current_elem->hmax();
354 198437400 : Real transient_part = _transient_term ? 4. / (_dt * _dt) : 0.;
355 396874800 : return -_alpha / 2. *
356 198437400 : std::pow(transient_part + (2. * U.norm() / h) * (2. * U.norm() / h) +
357 198437400 : 9. * (4. * nu / (h * h)) * (4. * nu / (h * h)),
358 198437400 : -1.5) *
359 198437400 : 2. * (2. * U.norm() / h) * 2. / h * U(comp) * _phi[_j][_qp] /
360 198437400 : (U.norm() + std::numeric_limits<double>::epsilon());
361 : }
362 :
363 : RealVectorValue
364 40551552 : INSBase::strongViscousTermLaplaceRZ() const
365 : {
366 : // To understand the code below, visit
367 : // https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates.
368 : // The u_r / r^2 term comes from the vector Laplacian. The -du_r/dr * 1/r term comes from
369 : // the scalar Laplacian. The scalar Laplacian in axisymmetric cylindrical coordinates is
370 : // equivalent to the Cartesian Laplacian plus a 1/r * df/dr term. And of course we are
371 : // applying a minus sign here because the strong form is -\nabala^2 * \vec{u}
372 :
373 40551552 : const auto r = _q_point[_qp](_rz_radial_coord);
374 : RealVectorValue rz_term;
375 40551552 : rz_term(0) = -_mu[_qp] * _grad_u_vel[_qp](_rz_radial_coord) / r;
376 40551552 : rz_term(1) = -_mu[_qp] * _grad_v_vel[_qp](_rz_radial_coord) / r;
377 : mooseAssert((_rz_radial_coord == 0) || (_rz_radial_coord == 1),
378 : "We expect X or Y as the possible radial coordinate");
379 40551552 : if (_rz_radial_coord == 0)
380 32487552 : rz_term(0) += _mu[_qp] * _u_vel[_qp] / (r * r);
381 : else
382 8064000 : rz_term(1) += _mu[_qp] * _v_vel[_qp] / (r * r);
383 :
384 40551552 : return rz_term;
385 : }
386 :
387 : RealVectorValue
388 33131376 : INSBase::dStrongViscDUCompLaplaceRZ(const unsigned int comp) const
389 : {
390 33131376 : const auto r = _q_point[_qp](_rz_radial_coord);
391 : RealVectorValue add_jac;
392 33131376 : add_jac(comp) = -_mu[_qp] * _grad_phi[_j][_qp](_rz_radial_coord) / r;
393 33131376 : if (comp == _rz_radial_coord)
394 16565688 : add_jac(comp) += _mu[_qp] * _phi[_j][_qp] / (r * r);
395 :
396 33131376 : return add_jac;
397 : }
398 :
399 : RealVectorValue
400 6366420 : INSBase::strongViscousTermTractionRZ() const
401 : {
402 6366420 : auto ret = strongViscousTermLaplaceRZ();
403 :
404 6366420 : const auto r = _q_point[_qp](_rz_radial_coord);
405 :
406 6366420 : const auto & grad_r_vel = (_rz_radial_coord == 0) ? _grad_u_vel[_qp] : _grad_v_vel[_qp];
407 6366420 : const auto & r_vel = (_rz_radial_coord == 0) ? _u_vel[_qp] : _v_vel[_qp];
408 6366420 : ret += -_mu[_qp] * grad_r_vel / r;
409 6366420 : ret(_rz_radial_coord) += _mu[_qp] * r_vel / (r * r);
410 :
411 6366420 : return ret;
412 : }
413 :
414 : RealVectorValue
415 4168800 : INSBase::dStrongViscDUCompTractionRZ(const unsigned int comp) const
416 : {
417 4168800 : auto ret = dStrongViscDUCompLaplaceRZ(comp);
418 4168800 : if (comp != _rz_radial_coord)
419 : return ret;
420 :
421 2084400 : const auto r = _q_point[_qp](_rz_radial_coord);
422 :
423 2084400 : ret += -_mu[_qp] * _grad_phi[_j][_qp] / r;
424 2084400 : ret(_rz_radial_coord) += _mu[_qp] * _phi[_j][_qp] / (r * r);
425 :
426 2084400 : return ret;
427 : }
|
