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 "LinearWCNSFVMomentumFlux.h"
11 : #include "MooseLinearVariableFV.h"
12 : #include "NS.h"
13 : #include "RhieChowMassFlux.h"
14 : #include "LinearFVBoundaryCondition.h"
15 : #include "LinearFVAdvectionDiffusionBC.h"
16 :
17 : registerMooseObject("NavierStokesApp", LinearWCNSFVMomentumFlux);
18 :
19 : InputParameters
20 1997 : LinearWCNSFVMomentumFlux::validParams()
21 : {
22 1997 : InputParameters params = LinearFVFluxKernel::validParams();
23 1997 : params.addClassDescription("Represents the matrix and right hand side contributions of the "
24 : "stress and advection terms of the momentum equation.");
25 3994 : params.addRequiredParam<SolverVariableName>("u", "The velocity in the x direction.");
26 3994 : params.addParam<SolverVariableName>("v", "The velocity in the y direction.");
27 3994 : params.addParam<SolverVariableName>("w", "The velocity in the z direction.");
28 3994 : params.addRequiredParam<UserObjectName>(
29 : "rhie_chow_user_object",
30 : "The rhie-chow user-object which is used to determine the face velocity.");
31 1997 : params.addRequiredParam<MooseFunctorName>(NS::mu, "The diffusion coefficient.");
32 3994 : MooseEnum momentum_component("x=0 y=1 z=2");
33 3994 : params.addRequiredParam<MooseEnum>(
34 : "momentum_component",
35 : momentum_component,
36 : "The component of the momentum equation that this kernel applies to.");
37 3994 : params.addParam<bool>(
38 : "use_nonorthogonal_correction",
39 3994 : true,
40 : "If the nonorthogonal correction should be used when computing the normal gradient.");
41 3994 : params.addParam<bool>(
42 3994 : "use_deviatoric_terms", false, "If deviatoric terms in the stress terms need to be used.");
43 :
44 3994 : params.addRequiredParam<InterpolationMethodName>(
45 : "advected_interp_method_name",
46 : "Name of the FVInterpolationMethod to use for the advected velocity.");
47 1997 : return params;
48 1997 : }
49 :
50 1070 : LinearWCNSFVMomentumFlux::LinearWCNSFVMomentumFlux(const InputParameters & params)
51 : : LinearFVFluxKernel(params),
52 : FVInterpolationMethodInterface(this),
53 1070 : _dim(_subproblem.mesh().dimension()),
54 1070 : _mass_flux_provider(getUserObject<RhieChowMassFlux>("rhie_chow_user_object")),
55 1070 : _mu(getFunctor<Real>(getParam<MooseFunctorName>(NS::mu))),
56 2140 : _use_nonorthogonal_correction(getParam<bool>("use_nonorthogonal_correction")),
57 2140 : _use_deviatoric_terms(getParam<bool>("use_deviatoric_terms")),
58 2140 : _adv_interp_method(getFVAdvectedInterpolationMethod(
59 : getParam<InterpolationMethodName>("advected_interp_method_name"))),
60 1070 : _face_mass_flux(0.0),
61 1070 : _boundary_normal_factor(1.0),
62 1070 : _stress_matrix_contribution(0.0),
63 1070 : _stress_rhs_contribution(0.0),
64 2140 : _index(getParam<MooseEnum>("momentum_component")),
65 1070 : _velocity_vars{nullptr, nullptr, nullptr},
66 1070 : _coord_type(getBlockCoordSystem()),
67 2140 : _rz_radial_coord(_fe_problem.mesh().getAxisymmetricRadialCoord())
68 : {
69 : // We only need gradients if the nonorthogonal correction is enabled or when we request the
70 : // computation of the deviatoric parts of the stress tensor.
71 1070 : if (_use_nonorthogonal_correction || _use_deviatoric_terms)
72 298 : _var.computeCellGradients();
73 :
74 1070 : if (_adv_interp_method.needsGradients())
75 0 : _var.computeCellGradients(_adv_interp_method.gradientLimiter());
76 :
77 2130 : auto get_velocity_var = [&](const std::string & param_name)
78 : {
79 2130 : return dynamic_cast<const MooseLinearVariableFVReal *>(
80 2130 : &_fe_problem.getVariable(_tid, getParam<SolverVariableName>(param_name)));
81 1070 : };
82 :
83 1070 : _velocity_vars[0] = get_velocity_var("u");
84 1070 : if (!_velocity_vars[0])
85 0 : paramError("u", "the u velocity must be a MooseLinearVariableFVReal.");
86 :
87 1070 : if (_dim >= 2)
88 : {
89 2012 : if (!params.isParamValid("v"))
90 0 : paramError("v", "In two or more dimensions, the v velocity must be supplied.");
91 1006 : _velocity_vars[1] = get_velocity_var("v");
92 1006 : if (!_velocity_vars[1])
93 0 : paramError("v",
94 : "In two or more dimensions, the v velocity must be supplied and it must be a "
95 : "MooseLinearVariableFVReal.");
96 : }
97 :
98 1070 : if (_dim >= 3)
99 : {
100 108 : if (!params.isParamValid("w"))
101 0 : paramError("w", "In three-dimensions, the w velocity must be supplied.");
102 54 : _velocity_vars[2] = get_velocity_var("w");
103 54 : if (!_velocity_vars[2])
104 0 : paramError("w",
105 : "In three-dimensions, the w velocity must be supplied and it must be a "
106 : "MooseLinearVariableFVReal.");
107 : }
108 1070 : }
109 :
110 : Real
111 34568947 : LinearWCNSFVMomentumFlux::computeElemMatrixContribution()
112 : {
113 34568947 : return (computeInternalAdvectionElemMatrixContribution() +
114 34568947 : computeInternalStressMatrixContribution()) *
115 34568947 : _current_face_area;
116 : }
117 :
118 : Real
119 34568947 : LinearWCNSFVMomentumFlux::computeNeighborMatrixContribution()
120 : {
121 34568947 : return (computeInternalAdvectionNeighborMatrixContribution() -
122 34568947 : computeInternalStressMatrixContribution()) *
123 34568947 : _current_face_area;
124 : }
125 :
126 : Real
127 34568947 : LinearWCNSFVMomentumFlux::computeElemRightHandSideContribution()
128 : {
129 34568947 : return (computeInternalStressRHSContribution() +
130 34568947 : _adv_interp_result.rhs_face_value * _face_mass_flux) *
131 34568947 : _current_face_area;
132 : }
133 :
134 : Real
135 34568947 : LinearWCNSFVMomentumFlux::computeNeighborRightHandSideContribution()
136 : {
137 34568947 : return -(computeInternalStressRHSContribution() +
138 34568947 : _adv_interp_result.rhs_face_value * _face_mass_flux) *
139 34568947 : _current_face_area;
140 : }
141 :
142 : Real
143 4904954 : LinearWCNSFVMomentumFlux::computeBoundaryMatrixContribution(const LinearFVBoundaryCondition & bc)
144 : {
145 : const auto * const adv_diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
146 :
147 : mooseAssert(adv_diff_bc, "This should be a valid BC!");
148 4904954 : return (computeStressBoundaryMatrixContribution(adv_diff_bc) +
149 4904954 : computeAdvectionBoundaryMatrixContribution(adv_diff_bc)) *
150 4904954 : _current_face_area;
151 : }
152 :
153 : Real
154 4904954 : LinearWCNSFVMomentumFlux::computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc)
155 : {
156 : const auto * const adv_diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
157 : mooseAssert(adv_diff_bc, "This should be a valid BC!");
158 4904954 : return (computeStressBoundaryRHSContribution(adv_diff_bc) +
159 4904954 : computeAdvectionBoundaryRHSContribution(adv_diff_bc)) *
160 4904954 : _current_face_area;
161 : }
162 :
163 : Real
164 34568947 : LinearWCNSFVMomentumFlux::computeInternalAdvectionElemMatrixContribution()
165 : {
166 34568947 : return _adv_interp_result.weights_matrix.first * _face_mass_flux;
167 : }
168 :
169 : Real
170 34568947 : LinearWCNSFVMomentumFlux::computeInternalAdvectionNeighborMatrixContribution()
171 : {
172 34568947 : return _adv_interp_result.weights_matrix.second * _face_mass_flux;
173 : }
174 :
175 : Real
176 69137894 : LinearWCNSFVMomentumFlux::computeInternalStressMatrixContribution()
177 : {
178 : // If we don't have the value yet, we compute it
179 69137894 : if (!_cached_matrix_contribution)
180 : {
181 34568947 : const auto face_arg = makeCDFace(*_current_face_info);
182 :
183 : // If we requested nonorthogonal correction, we use the normal component of the
184 : // cell to face vector.
185 34568947 : const auto d = _use_nonorthogonal_correction
186 34568947 : ? std::abs(_current_face_info->dCN() * _current_face_info->normal())
187 32378899 : : _current_face_info->dCNMag();
188 :
189 : // Cache the matrix contribution
190 34568947 : _stress_matrix_contribution = _mu(face_arg, determineState()) / d;
191 34568947 : _cached_matrix_contribution = true;
192 : }
193 :
194 69137894 : return _stress_matrix_contribution;
195 : }
196 :
197 : Real
198 69137894 : LinearWCNSFVMomentumFlux::computeInternalStressRHSContribution()
199 : {
200 : // We can have contributions to the right hand side in two occasions:
201 : // (1) when we use nonorthogonal correction for the normal gradients
202 : // (2) when we request the deviatoric parts of the stress tensor. (needed for space-dependent
203 : // viscosities for example)
204 69137894 : if (!_cached_rhs_contribution)
205 : {
206 : // scenario (1), we need to add the nonorthogonal correction. In 1D, we don't have
207 : // any correction so we just skip this part
208 34568947 : if (_dim > 1 && _use_nonorthogonal_correction)
209 : {
210 2190048 : const auto face_arg = makeCDFace(*_current_face_info);
211 2190048 : const auto state_arg = determineState();
212 :
213 : // Get the gradients from the adjacent cells
214 2190048 : const auto grad_elem = _var.gradSln(*_current_face_info->elemInfo(), state_arg);
215 2190048 : const auto & grad_neighbor = _var.gradSln(*_current_face_info->neighborInfo(), state_arg);
216 :
217 : // Interpolate the two gradients to the face
218 : const auto interp_coeffs =
219 2190048 : interpCoeffs(Moose::FV::InterpMethod::Average, *_current_face_info, true);
220 :
221 : const auto correction_vector =
222 : _current_face_info->normal() -
223 2190048 : 1 / (_current_face_info->normal() * _current_face_info->eCN()) *
224 : _current_face_info->eCN();
225 :
226 : // Cache the matrix contribution
227 2190048 : _stress_rhs_contribution +=
228 2190048 : _mu(face_arg, state_arg) *
229 : (interp_coeffs.first * grad_elem + interp_coeffs.second * grad_neighbor) *
230 : correction_vector;
231 : }
232 : // scenario (2), we will have to account for the deviatoric parts of the stress tensor.
233 34568947 : if (_use_deviatoric_terms)
234 : {
235 2614032 : const auto state_arg = determineState();
236 :
237 : // Interpolate the two gradients to the face
238 : const auto interp_coeffs =
239 2614032 : interpCoeffs(Moose::FV::InterpMethod::Average, *_current_face_info, true);
240 :
241 10456128 : RealGradient grad_elem[3];
242 10456128 : RealGradient grad_neighbor[3];
243 : Real trace_elem = 0;
244 : Real trace_neighbor = 0;
245 : RealVectorValue deviatoric_vector_elem;
246 : RealVectorValue deviatoric_vector_neighbor;
247 :
248 : // Loop over every velocity component so we can form the symmetric gradient pieces
249 7842096 : for (const auto dir : make_range(_dim))
250 : {
251 5228064 : grad_elem[dir] = velocityVar(dir).gradSln(*_current_face_info->elemInfo(), state_arg);
252 : grad_neighbor[dir] =
253 5228064 : velocityVar(dir).gradSln(*_current_face_info->neighborInfo(), state_arg);
254 5228064 : trace_elem += grad_elem[dir](dir);
255 5228064 : trace_neighbor += grad_neighbor[dir](dir);
256 : }
257 :
258 2614032 : const auto face_arg = makeCDFace(*_current_face_info);
259 :
260 2614032 : if (_coord_type == Moose::CoordinateSystemType::COORD_RZ)
261 : {
262 : Real elem_value = 0.0;
263 : Real neighbor_value = 0.0;
264 235848 : const auto & radial_var = velocityVar(_rz_radial_coord);
265 235848 : elem_value = radial_var.getElemValue(*_current_face_info->elemInfo(), state_arg) /
266 235848 : _current_face_info->elemInfo()->centroid()(_rz_radial_coord);
267 235848 : neighbor_value = radial_var.getElemValue(*_current_face_info->neighborInfo(), state_arg) /
268 235848 : _current_face_info->neighborInfo()->centroid()(_rz_radial_coord);
269 :
270 235848 : trace_elem += elem_value;
271 235848 : trace_neighbor += neighbor_value;
272 : }
273 :
274 : // Assemble the explicit transpose/trace contribution component by component
275 7842096 : for (const auto dir : make_range(_dim))
276 : {
277 5228064 : grad_elem[dir](dir) -= 2. / 3 * trace_elem;
278 5228064 : grad_neighbor[dir](dir) -= 2. / 3 * trace_neighbor;
279 :
280 5228064 : deviatoric_vector_elem(dir) = grad_elem[dir](_index);
281 5228064 : deviatoric_vector_neighbor(dir) = grad_neighbor[dir](_index);
282 : }
283 :
284 2614032 : _stress_rhs_contribution += _mu(face_arg, state_arg) *
285 : (interp_coeffs.first * deviatoric_vector_elem +
286 : interp_coeffs.second * deviatoric_vector_neighbor) *
287 2614032 : _current_face_info->normal();
288 : }
289 34568947 : _cached_rhs_contribution = true;
290 : }
291 :
292 69137894 : return _stress_rhs_contribution;
293 : }
294 :
295 : Real
296 4904954 : LinearWCNSFVMomentumFlux::computeStressBoundaryMatrixContribution(
297 : const LinearFVAdvectionDiffusionBC * bc)
298 : {
299 4904954 : auto grad_contrib = bc->computeBoundaryGradientMatrixContribution();
300 : // If the boundary condition does not include the diffusivity contribution then
301 : // add it here.
302 4904954 : if (!bc->includesMaterialPropertyMultiplier())
303 : {
304 3905729 : const auto face_arg = singleSidedFaceArg(_current_face_info);
305 3905729 : grad_contrib *= _mu(face_arg, determineState());
306 : }
307 :
308 4904954 : return grad_contrib;
309 : }
310 :
311 : Real
312 4904954 : LinearWCNSFVMomentumFlux::computeStressBoundaryRHSContribution(
313 : const LinearFVAdvectionDiffusionBC * bc)
314 : {
315 4904954 : const auto face_arg = singleSidedFaceArg(_current_face_info);
316 4904954 : auto grad_contrib = bc->computeBoundaryGradientRHSContribution();
317 : // If the boundary condition does not include the diffusivity contribution then
318 : // add it here.
319 4904954 : if (!bc->includesMaterialPropertyMultiplier())
320 3905729 : grad_contrib *= _mu(face_arg, determineState());
321 :
322 : // We add the nonorthogonal corrector for the face here. Potential idea: we could do
323 : // this in the boundary condition too. For now, however, we keep it like this.
324 4904954 : if (_use_nonorthogonal_correction && bc->useBoundaryGradientExtrapolation())
325 : {
326 : // We support internal boundaries as well. In that case we have to decide on which side
327 : // of the boundary we are on.
328 274148 : const auto elem_info = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
329 274148 : ? _current_face_info->elemInfo()
330 17952 : : _current_face_info->neighborInfo();
331 :
332 : // Unit vector to the boundary. Unfortunately, we have to recompute it because the value
333 : // stored in the face info is only correct for external boundaries
334 274148 : const auto e_Cf = _current_face_info->faceCentroid() - elem_info->centroid();
335 : const auto correction_vector =
336 274148 : _current_face_info->normal() - 1 / (_current_face_info->normal() * e_Cf) * e_Cf;
337 :
338 274148 : const auto state_arg = determineState();
339 274148 : grad_contrib += _mu(face_arg, state_arg) * _var.gradSln(*elem_info, state_arg) *
340 : _boundary_normal_factor * correction_vector;
341 : }
342 :
343 4904954 : if (_use_deviatoric_terms && bc->useBoundaryGradientExtrapolation())
344 : {
345 : // We might be on a face which is an internal boundary so we want to make sure we
346 : // get the gradient from the right side.
347 698460 : const auto elem_info = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
348 698460 : ? _current_face_info->elemInfo()
349 10664 : : _current_face_info->neighborInfo();
350 :
351 698460 : const auto state_arg = determineState();
352 :
353 2793840 : RealGradient grad_elem[3];
354 : Real trace_elem = 0;
355 : RealVectorValue deviatoric_vector_elem;
356 :
357 2095380 : for (const auto dir : make_range(_dim))
358 : {
359 1396920 : grad_elem[dir] = velocityVar(dir).gradSln(*elem_info, state_arg);
360 1396920 : trace_elem += grad_elem[dir](dir);
361 : }
362 :
363 698460 : if (_coord_type == Moose::CoordinateSystemType::COORD_RZ)
364 : {
365 29052 : const auto & radial_var = velocityVar(_rz_radial_coord);
366 : const Real elem_value =
367 29052 : radial_var.getElemValue(*elem_info, state_arg) / elem_info->centroid()(_rz_radial_coord);
368 29052 : trace_elem += elem_value;
369 : }
370 :
371 2095380 : for (const auto dir : make_range(_dim))
372 : {
373 1396920 : grad_elem[dir](dir) -= 2. / 3 * trace_elem;
374 1396920 : deviatoric_vector_elem(dir) = grad_elem[dir](_index);
375 : }
376 :
377 : // We support internal boundaries too so we have to make sure the normal points always outward
378 698460 : grad_contrib += _mu(face_arg, state_arg) * deviatoric_vector_elem * _boundary_normal_factor *
379 698460 : _current_face_info->normal();
380 : }
381 :
382 4904954 : return grad_contrib;
383 : }
384 :
385 : Real
386 4904954 : LinearWCNSFVMomentumFlux::computeAdvectionBoundaryMatrixContribution(
387 : const LinearFVAdvectionDiffusionBC * bc)
388 : {
389 4904954 : const auto boundary_value_matrix_contrib = bc->computeBoundaryValueMatrixContribution();
390 4904954 : return boundary_value_matrix_contrib * _face_mass_flux;
391 : }
392 :
393 : Real
394 4904954 : LinearWCNSFVMomentumFlux::computeAdvectionBoundaryRHSContribution(
395 : const LinearFVAdvectionDiffusionBC * bc)
396 : {
397 4904954 : const auto boundary_value_rhs_contrib = bc->computeBoundaryValueRHSContribution();
398 4904954 : return -boundary_value_rhs_contrib * _face_mass_flux;
399 : }
400 :
401 : void
402 39473901 : LinearWCNSFVMomentumFlux::setupFaceData(const FaceInfo * face_info)
403 : {
404 39473901 : LinearFVFluxKernel::setupFaceData(face_info);
405 :
406 : // Multiplier that ensures the normal of the boundary always points outwards, even in cases
407 : // when the boundary is within the mesh.
408 39473901 : _boundary_normal_factor = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) ? 1.0 : -1.0;
409 :
410 : // Caching the mass flux on the face which will be reused in the advection term's matrix and
411 : // right hand side contributions
412 39473901 : _face_mass_flux = _mass_flux_provider.getMassFlux(*face_info);
413 :
414 39473901 : if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
415 : {
416 34568947 : const auto state = determineState();
417 34568947 : const auto & elem_info = *_current_face_info->elemInfo();
418 : const auto & neighbor_info = *_current_face_info->neighborInfo();
419 :
420 34568947 : const Real elem_value = _var.getElemValue(elem_info, state);
421 34568947 : const Real neighbor_value = _var.getElemValue(neighbor_info, state);
422 :
423 34568947 : if (_adv_interp_method.needsGradients())
424 : {
425 0 : const auto limiter_type = _adv_interp_method.gradientLimiter();
426 0 : _elem_grad_storage = _var.gradSln(elem_info, state, limiter_type);
427 0 : _neighbor_grad_storage = _var.gradSln(neighbor_info, state, limiter_type);
428 : }
429 :
430 34568947 : _adv_interp_result = _adv_interp_method.advectedInterpolate(*_current_face_info,
431 : elem_value,
432 : neighbor_value,
433 34568947 : &_elem_grad_storage,
434 34568947 : &_neighbor_grad_storage,
435 : _face_mass_flux);
436 : }
437 :
438 : // We'll have to set this to zero to make sure that we don't accumulate values over multiple
439 : // faces. The matrix contribution should be fine.
440 39473901 : _stress_rhs_contribution = 0;
441 39473901 : }
442 :
443 : const MooseLinearVariableFVReal &
444 12117948 : LinearWCNSFVMomentumFlux::velocityVar(unsigned int dir) const
445 : {
446 : mooseAssert(dir < _velocity_vars.size() && _velocity_vars[dir],
447 : "Velocity variable for requested direction is not available.");
448 12117948 : return *_velocity_vars[dir];
449 : }
|
