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 6756 : INSBase::validParams()
16 : {
17 6756 : InputParameters params = Kernel::validParams();
18 :
19 6756 : 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 13512 : params.addRequiredCoupledVar("u", "x-velocity");
24 13512 : params.addCoupledVar("v", 0, "y-velocity"); // only required in 2D and 3D
25 13512 : params.addCoupledVar("w", 0, "z-velocity"); // only required in 3D
26 6756 : params.addRequiredCoupledVar(NS::pressure, "pressure");
27 :
28 6756 : params.addParam<RealVectorValue>(
29 6756 : "gravity", RealVectorValue(0, 0, 0), "Direction of the gravity vector");
30 :
31 13512 : params.addParam<MaterialPropertyName>("mu_name", "mu", "The name of the dynamic viscosity");
32 13512 : params.addParam<MaterialPropertyName>("rho_name", "rho", "The name of the density");
33 :
34 13512 : params.addParam<Real>("alpha", 1., "Multiplicative factor on the stabilization parameter tau.");
35 13512 : params.addParam<bool>(
36 13512 : "laplace", true, "Whether the viscous term of the momentum equations is in laplace form.");
37 13512 : params.addParam<bool>("convective_term", true, "Whether to include the convective term.");
38 13512 : params.addParam<bool>("transient_term",
39 13512 : false,
40 : "Whether there should be a transient term in the momentum residuals.");
41 13512 : params.addCoupledVar("disp_x", "The x displacement");
42 13512 : params.addCoupledVar("disp_y", "The y displacement");
43 13512 : params.addCoupledVar("disp_z", "The z displacement");
44 13512 : params.addParam<bool>("picard",
45 13512 : false,
46 : "Whether we are applying a Picard strategy in which case we will linearize "
47 : "the nonlinear convective term.");
48 :
49 6756 : return params;
50 0 : }
51 :
52 3715 : INSBase::INSBase(const InputParameters & parameters)
53 : : Kernel(parameters),
54 3715 : _second_phi(_assembly.secondPhi()),
55 :
56 : // Coupled variables
57 3715 : _u_vel(coupledValue("u")),
58 3715 : _v_vel(coupledValue("v")),
59 3715 : _w_vel(coupledValue("w")),
60 3715 : _p(coupledValue(NS::pressure)),
61 :
62 7430 : _picard(getParam<bool>("picard")),
63 3715 : _u_vel_previous_nl(_picard ? &coupledValuePreviousNL("u") : nullptr),
64 3715 : _v_vel_previous_nl(_picard ? &coupledValuePreviousNL("v") : nullptr),
65 3715 : _w_vel_previous_nl(_picard ? &coupledValuePreviousNL("w") : nullptr),
66 :
67 : // Gradients
68 3715 : _grad_u_vel(coupledGradient("u")),
69 3715 : _grad_v_vel(coupledGradient("v")),
70 3715 : _grad_w_vel(coupledGradient("w")),
71 3715 : _grad_p(coupledGradient(NS::pressure)),
72 :
73 : // second derivative tensors
74 3715 : _second_u_vel(coupledSecond("u")),
75 3715 : _second_v_vel(coupledSecond("v")),
76 3715 : _second_w_vel(coupledSecond("w")),
77 :
78 : // time derivatives
79 3715 : _u_vel_dot(_is_transient ? coupledDot("u") : _zero),
80 3715 : _v_vel_dot(_is_transient ? coupledDot("v") : _zero),
81 3715 : _w_vel_dot(_is_transient ? coupledDot("w") : _zero),
82 :
83 : // derivatives of time derivatives
84 3715 : _d_u_vel_dot_du(_is_transient ? coupledDotDu("u") : _zero),
85 3715 : _d_v_vel_dot_dv(_is_transient ? coupledDotDu("v") : _zero),
86 3715 : _d_w_vel_dot_dw(_is_transient ? coupledDotDu("w") : _zero),
87 :
88 : // Variable numberings
89 3715 : _u_vel_var_number(coupled("u")),
90 3715 : _v_vel_var_number(coupled("v")),
91 3715 : _w_vel_var_number(coupled("w")),
92 3715 : _p_var_number(coupled(NS::pressure)),
93 :
94 7430 : _gravity(getParam<RealVectorValue>("gravity")),
95 :
96 : // Material properties
97 7430 : _mu(getMaterialProperty<Real>("mu_name")),
98 7430 : _rho(getMaterialProperty<Real>("rho_name")),
99 :
100 7430 : _alpha(getParam<Real>("alpha")),
101 7430 : _laplace(getParam<bool>("laplace")),
102 7430 : _convective_term(getParam<bool>("convective_term")),
103 7430 : _transient_term(getParam<bool>("transient_term")),
104 :
105 : // Displacements for mesh velocity for ALE simulations
106 7430 : _disps_provided(isParamValid("disp_x")),
107 7430 : _disp_x_dot(isParamValid("disp_x") ? coupledDot("disp_x") : _zero),
108 7430 : _disp_y_dot(isParamValid("disp_y") ? coupledDot("disp_y") : _zero),
109 7430 : _disp_z_dot(isParamValid("disp_z") ? coupledDot("disp_z") : _zero),
110 :
111 7430 : _rz_radial_coord(_mesh.getAxisymmetricRadialCoord())
112 : {
113 3715 : 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 3715 : }
118 :
119 : RealVectorValue
120 1295098124 : INSBase::relativeVelocity() const
121 : {
122 1295098124 : auto U = _picard ? RealVectorValue((*_u_vel_previous_nl)[_qp],
123 0 : (*_v_vel_previous_nl)[_qp],
124 0 : (*_w_vel_previous_nl)[_qp])
125 1295098124 : : RealVectorValue(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
126 1295098124 : if (_disps_provided)
127 0 : U -= RealVectorValue{_disp_x_dot[_qp], _disp_y_dot[_qp], _disp_z_dot[_qp]};
128 1295098124 : return U;
129 : }
130 :
131 : RealVectorValue
132 520348788 : INSBase::convectiveTerm()
133 : {
134 520348788 : const auto U = _picard ? RealVectorValue((*_u_vel_previous_nl)[_qp],
135 5650560 : (*_v_vel_previous_nl)[_qp],
136 5650560 : (*_w_vel_previous_nl)[_qp])
137 520348788 : : RealVectorValue(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
138 520348788 : return _rho[_qp] *
139 520348788 : RealVectorValue(U * _grad_u_vel[_qp], U * _grad_v_vel[_qp], U * _grad_w_vel[_qp]);
140 : }
141 :
142 : RealVectorValue
143 1635146406 : INSBase::dConvecDUComp(unsigned comp)
144 : {
145 1635146406 : if (_picard)
146 : {
147 : RealVectorValue U(
148 98910720 : (*_u_vel_previous_nl)[_qp], (*_v_vel_previous_nl)[_qp], (*_w_vel_previous_nl)[_qp]);
149 : RealVectorValue convective_term;
150 98910720 : convective_term(comp) = _rho[_qp] * U * _grad_phi[_j][_qp];
151 :
152 98910720 : return convective_term;
153 : }
154 : else
155 : {
156 1536235686 : RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
157 : RealVectorValue d_U_d_comp(0, 0, 0);
158 1536235686 : d_U_d_comp(comp) = _phi[_j][_qp];
159 :
160 1536235686 : RealVectorValue convective_term = _rho[_qp] * RealVectorValue(d_U_d_comp * _grad_u_vel[_qp],
161 1536235686 : d_U_d_comp * _grad_v_vel[_qp],
162 1536235686 : d_U_d_comp * _grad_w_vel[_qp]);
163 1536235686 : convective_term(comp) += _rho[_qp] * U * _grad_phi[_j][_qp];
164 :
165 1536235686 : return convective_term;
166 : }
167 : }
168 :
169 : RealVectorValue
170 326429666 : INSBase::strongViscousTermLaplace()
171 : {
172 326429666 : return -_mu[_qp] *
173 326429666 : RealVectorValue(_second_u_vel[_qp].tr(), _second_v_vel[_qp].tr(), _second_w_vel[_qp].tr());
174 : }
175 :
176 : RealVectorValue
177 37299204 : INSBase::strongViscousTermTraction()
178 : {
179 37299204 : return INSBase::strongViscousTermLaplace() -
180 37299204 : _mu[_qp] *
181 111897612 : (_second_u_vel[_qp].row(0) + _second_v_vel[_qp].row(1) + _second_w_vel[_qp].row(2));
182 : }
183 :
184 : RealVectorValue
185 256225032 : INSBase::dStrongViscDUCompLaplace(unsigned comp)
186 : {
187 : RealVectorValue viscous_term(0, 0, 0);
188 256225032 : viscous_term(comp) = -_mu[_qp] * _second_phi[_j][_qp].tr();
189 :
190 256225032 : return viscous_term;
191 : }
192 :
193 : RealVectorValue
194 30652668 : INSBase::dStrongViscDUCompTraction(unsigned comp)
195 : {
196 : RealVectorValue viscous_term(0, 0, 0);
197 30652668 : viscous_term(comp) = -_mu[_qp] * (_second_phi[_j][_qp](0, 0) + _second_phi[_j][_qp](1, 1) +
198 30652668 : _second_phi[_j][_qp](2, 2));
199 122610672 : for (unsigned i = 0; i < 3; i++)
200 91958004 : viscous_term(i) += -_mu[_qp] * _second_phi[_j][_qp](i, comp);
201 :
202 30652668 : 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 674609826 : INSBase::strongPressureTerm()
266 : {
267 674609826 : return _grad_p[_qp];
268 : }
269 :
270 : Real
271 85195500 : INSBase::weakPressureTerm()
272 : {
273 85195500 : return -_p[_qp];
274 : }
275 :
276 : RealVectorValue
277 175493519 : INSBase::dStrongPressureDPressure()
278 : {
279 175493519 : return _grad_phi[_j][_qp];
280 : }
281 :
282 : Real
283 248922631 : INSBase::dWeakPressureDPressure()
284 : {
285 248922631 : return -_phi[_j][_qp];
286 : }
287 :
288 : RealVectorValue
289 759805326 : INSBase::gravityTerm()
290 : {
291 759805326 : return -_rho[_qp] * _gravity;
292 : }
293 :
294 : RealVectorValue
295 12686664 : INSBase::timeDerivativeTerm()
296 : {
297 12686664 : return _rho[_qp] * RealVectorValue(_u_vel_dot[_qp], _v_vel_dot[_qp], _w_vel_dot[_qp]);
298 : }
299 :
300 : RealVectorValue
301 11707548 : INSBase::dTimeDerivativeDUComp(unsigned comp)
302 : {
303 11707548 : Real base = _rho[_qp] * _phi[_j][_qp];
304 11707548 : switch (comp)
305 : {
306 5853774 : case 0:
307 5853774 : return RealVectorValue(base * _d_u_vel_dot_du[_qp], 0, 0);
308 :
309 5853774 : case 1:
310 5853774 : 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 664086380 : INSBase::tau()
322 : {
323 664086380 : Real nu = _mu[_qp] / _rho[_qp];
324 664086380 : const auto U = relativeVelocity();
325 664086380 : Real h = _current_elem->hmax();
326 664086380 : Real transient_part = _transient_term ? 4. / (_dt * _dt) : 0.;
327 664086380 : return _alpha / std::sqrt(transient_part + (2. * U.norm() / h) * (2. * U.norm() / h) +
328 664086380 : 9. * (4. * nu / (h * h)) * (4. * nu / (h * h)));
329 : }
330 :
331 : Real
332 6540 : INSBase::tauNodal()
333 : {
334 6540 : Real nu = _mu[_qp] / _rho[_qp];
335 6540 : const auto U = relativeVelocity();
336 6540 : Real h = _current_elem->hmax();
337 : Real xi;
338 6540 : 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 6540 : return h / (2 * U.norm()) * xi;
346 : }
347 :
348 : Real
349 286877700 : INSBase::dTauDUComp(const unsigned int comp)
350 : {
351 286877700 : Real nu = _mu[_qp] / _rho[_qp];
352 286877700 : const auto U = relativeVelocity();
353 286877700 : Real h = _current_elem->hmax();
354 286877700 : Real transient_part = _transient_term ? 4. / (_dt * _dt) : 0.;
355 573755400 : return -_alpha / 2. *
356 286877700 : std::pow(transient_part + (2. * U.norm() / h) * (2. * U.norm() / h) +
357 286877700 : 9. * (4. * nu / (h * h)) * (4. * nu / (h * h)),
358 286877700 : -1.5) *
359 286877700 : 2. * (2. * U.norm() / h) * 2. / h * U(comp) * _phi[_j][_qp] /
360 286877700 : (U.norm() + std::numeric_limits<double>::epsilon());
361 : }
362 :
363 : RealVectorValue
364 61593786 : 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 61593786 : const auto r = _q_point[_qp](_rz_radial_coord);
374 : RealVectorValue rz_term;
375 61593786 : rz_term(0) = -_mu[_qp] * _grad_u_vel[_qp](_rz_radial_coord) / r;
376 61593786 : 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 61593786 : if (_rz_radial_coord == 0)
380 49497786 : rz_term(0) += _mu[_qp] * _u_vel[_qp] / (r * r);
381 : else
382 12096000 : rz_term(1) += _mu[_qp] * _v_vel[_qp] / (r * r);
383 :
384 61593786 : return rz_term;
385 : }
386 :
387 : RealVectorValue
388 49836060 : INSBase::dStrongViscDUCompLaplaceRZ(const unsigned int comp) const
389 : {
390 49836060 : const auto r = _q_point[_qp](_rz_radial_coord);
391 : RealVectorValue add_jac;
392 49836060 : add_jac(comp) = -_mu[_qp] * _grad_phi[_j][_qp](_rz_radial_coord) / r;
393 49836060 : if (comp == _rz_radial_coord)
394 24918030 : add_jac(comp) += _mu[_qp] * _phi[_j][_qp] / (r * r);
395 :
396 49836060 : return add_jac;
397 : }
398 :
399 : RealVectorValue
400 10241640 : INSBase::strongViscousTermTractionRZ() const
401 : {
402 10241640 : auto ret = strongViscousTermLaplaceRZ();
403 :
404 10241640 : const auto r = _q_point[_qp](_rz_radial_coord);
405 :
406 10241640 : const auto & grad_r_vel = (_rz_radial_coord == 0) ? _grad_u_vel[_qp] : _grad_v_vel[_qp];
407 10241640 : const auto & r_vel = (_rz_radial_coord == 0) ? _u_vel[_qp] : _v_vel[_qp];
408 10241640 : ret += -_mu[_qp] * grad_r_vel / r;
409 10241640 : ret(_rz_radial_coord) += _mu[_qp] * r_vel / (r * r);
410 :
411 10241640 : return ret;
412 : }
413 :
414 : RealVectorValue
415 6324480 : INSBase::dStrongViscDUCompTractionRZ(const unsigned int comp) const
416 : {
417 6324480 : auto ret = dStrongViscDUCompLaplaceRZ(comp);
418 6324480 : if (comp != _rz_radial_coord)
419 : return ret;
420 :
421 3162240 : const auto r = _q_point[_qp](_rz_radial_coord);
422 :
423 3162240 : ret += -_mu[_qp] * _grad_phi[_j][_qp] / r;
424 3162240 : ret(_rz_radial_coord) += _mu[_qp] * _phi[_j][_qp] / (r * r);
425 :
426 3162240 : return ret;
427 : }
|
