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