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