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 3310 : LinearWCNSFVMomentumFlux::validParams()
22 : {
23 3310 : InputParameters params = LinearFVFluxKernel::validParams();
24 3310 : params.addClassDescription("Represents the matrix and right hand side contributions of the "
25 : "stress and advection terms of the momentum equation.");
26 6620 : params.addRequiredParam<SolverVariableName>("u", "The velocity in the x direction.");
27 6620 : params.addParam<SolverVariableName>("v", "The velocity in the y direction.");
28 6620 : params.addParam<SolverVariableName>("w", "The velocity in the z direction.");
29 6620 : params.addRequiredParam<UserObjectName>(
30 : "rhie_chow_user_object",
31 : "The rhie-chow user-object which is used to determine the face velocity.");
32 3310 : params.addRequiredParam<MooseFunctorName>(NS::mu, "The diffusion coefficient.");
33 6620 : MooseEnum momentum_component("x=0 y=1 z=2");
34 6620 : params.addRequiredParam<MooseEnum>(
35 : "momentum_component",
36 : momentum_component,
37 : "The component of the momentum equation that this kernel applies to.");
38 6620 : params.addParam<bool>(
39 : "use_nonorthogonal_correction",
40 6620 : true,
41 : "If the nonorthogonal correction should be used when computing the normal gradient.");
42 6620 : params.addParam<bool>(
43 6620 : "use_deviatoric_terms", false, "If deviatoric terms in the stress terms need to be used.");
44 :
45 3310 : params += Moose::FV::advectedInterpolationParameter();
46 3310 : return params;
47 3310 : }
48 :
49 1745 : LinearWCNSFVMomentumFlux::LinearWCNSFVMomentumFlux(const InputParameters & params)
50 : : LinearFVFluxKernel(params),
51 1745 : _dim(_subproblem.mesh().dimension()),
52 1745 : _mass_flux_provider(getUserObject<RhieChowMassFlux>("rhie_chow_user_object")),
53 1745 : _mu(getFunctor<Real>(getParam<MooseFunctorName>(NS::mu))),
54 3490 : _use_nonorthogonal_correction(getParam<bool>("use_nonorthogonal_correction")),
55 3490 : _use_deviatoric_terms(getParam<bool>("use_deviatoric_terms")),
56 1745 : _advected_interp_coeffs(std::make_pair<Real, Real>(0, 0)),
57 1745 : _face_mass_flux(0.0),
58 1745 : _boundary_normal_factor(1.0),
59 1745 : _stress_matrix_contribution(0.0),
60 1745 : _stress_rhs_contribution(0.0),
61 3490 : _index(getParam<MooseEnum>("momentum_component")),
62 1745 : _u_var(dynamic_cast<const MooseLinearVariableFVReal *>(
63 3490 : &_fe_problem.getVariable(_tid, getParam<SolverVariableName>("u")))),
64 1745 : _v_var(params.isParamValid("v")
65 3424 : ? dynamic_cast<const MooseLinearVariableFVReal *>(
66 6782 : &_fe_problem.getVariable(_tid, getParam<SolverVariableName>("v")))
67 : : nullptr),
68 1745 : _w_var(params.isParamValid("w")
69 1826 : ? dynamic_cast<const MooseLinearVariableFVReal *>(
70 1988 : &_fe_problem.getVariable(_tid, getParam<SolverVariableName>("w")))
71 1745 : : 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 1745 : if (_use_nonorthogonal_correction || _use_deviatoric_terms)
76 440 : _var.computeCellGradients();
77 :
78 1745 : Moose::FV::setInterpolationMethod(*this, _advected_interp_method, "advected_interp_method");
79 :
80 1745 : if (!_u_var)
81 0 : paramError("u", "the u velocity must be a MooseLinearVariableFVReal.");
82 :
83 1745 : 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 1745 : 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 1745 : }
93 :
94 : Real
95 51611281 : LinearWCNSFVMomentumFlux::computeElemMatrixContribution()
96 : {
97 51611281 : return (computeInternalAdvectionElemMatrixContribution() +
98 51611281 : computeInternalStressMatrixContribution()) *
99 51611281 : _current_face_area;
100 : }
101 :
102 : Real
103 51611281 : LinearWCNSFVMomentumFlux::computeNeighborMatrixContribution()
104 : {
105 51611281 : return (computeInternalAdvectionNeighborMatrixContribution() -
106 51611281 : computeInternalStressMatrixContribution()) *
107 51611281 : _current_face_area;
108 : }
109 :
110 : Real
111 51611281 : LinearWCNSFVMomentumFlux::computeElemRightHandSideContribution()
112 : {
113 51611281 : return computeInternalStressRHSContribution() * _current_face_area;
114 : }
115 :
116 : Real
117 51611281 : LinearWCNSFVMomentumFlux::computeNeighborRightHandSideContribution()
118 : {
119 51611281 : return -computeInternalStressRHSContribution() * _current_face_area;
120 : }
121 :
122 : Real
123 6357018 : 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 6357018 : return (computeStressBoundaryMatrixContribution(adv_diff_bc) +
128 6357018 : computeAdvectionBoundaryMatrixContribution(adv_diff_bc)) *
129 6357018 : _current_face_area;
130 : }
131 :
132 : Real
133 6357018 : 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 6357018 : return (computeStressBoundaryRHSContribution(adv_diff_bc) +
139 6357018 : computeAdvectionBoundaryRHSContribution(adv_diff_bc)) *
140 6357018 : _current_face_area;
141 : }
142 :
143 : Real
144 51611281 : LinearWCNSFVMomentumFlux::computeInternalAdvectionElemMatrixContribution()
145 : {
146 51611281 : return _advected_interp_coeffs.first * _face_mass_flux;
147 : }
148 :
149 : Real
150 51611281 : LinearWCNSFVMomentumFlux::computeInternalAdvectionNeighborMatrixContribution()
151 : {
152 51611281 : return _advected_interp_coeffs.second * _face_mass_flux;
153 : }
154 :
155 : Real
156 103222562 : LinearWCNSFVMomentumFlux::computeInternalStressMatrixContribution()
157 : {
158 : // If we don't have the value yet, we compute it
159 103222562 : if (!_cached_matrix_contribution)
160 : {
161 51611281 : 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 51611281 : const auto d = _use_nonorthogonal_correction
166 51611281 : ? std::abs(_current_face_info->dCN() * _current_face_info->normal())
167 48147457 : : _current_face_info->dCNMag();
168 :
169 : // Cache the matrix contribution
170 51611281 : _stress_matrix_contribution = _mu(face_arg, determineState()) / d;
171 51611281 : _cached_matrix_contribution = true;
172 : }
173 :
174 103222562 : return _stress_matrix_contribution;
175 : }
176 :
177 : Real
178 103222562 : 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 103222562 : 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 51611281 : if (_dim > 1 && _use_nonorthogonal_correction)
189 : {
190 3463824 : const auto face_arg = makeCDFace(*_current_face_info);
191 3463824 : const auto state_arg = determineState();
192 :
193 : // Get the gradients from the adjacent cells
194 3463824 : const auto grad_elem = _var.gradSln(*_current_face_info->elemInfo());
195 3463824 : 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 3463824 : interpCoeffs(Moose::FV::InterpMethod::Average, *_current_face_info, true);
200 :
201 : const auto correction_vector =
202 : _current_face_info->normal() -
203 3463824 : 1 / (_current_face_info->normal() * _current_face_info->eCN()) *
204 : _current_face_info->eCN();
205 :
206 : // Cache the matrix contribution
207 3463824 : _stress_rhs_contribution +=
208 3463824 : _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 51611281 : if (_use_deviatoric_terms)
214 : {
215 : // Interpolate the two gradients to the face
216 : const auto interp_coeffs =
217 8581032 : interpCoeffs(Moose::FV::InterpMethod::Average, *_current_face_info, true);
218 :
219 8581032 : const auto u_grad_elem = _u_var->gradSln(*_current_face_info->elemInfo());
220 8581032 : 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 8581032 : deviatoric_vector_elem(0) = u_grad_elem(_index);
228 8581032 : deviatoric_vector_neighbor(0) = u_grad_neighbor(_index);
229 8581032 : trace_elem += u_grad_elem(0);
230 8581032 : trace_neighbor += u_grad_neighbor(0);
231 :
232 8581032 : const auto face_arg = makeCDFace(*_current_face_info);
233 8581032 : const auto state_arg = determineState();
234 :
235 8581032 : if (_dim > 1)
236 : {
237 8581032 : const auto v_grad_elem = _v_var->gradSln(*_current_face_info->elemInfo());
238 8581032 : const auto v_grad_neighbor = _v_var->gradSln(*_current_face_info->neighborInfo());
239 :
240 8581032 : deviatoric_vector_elem(1) = v_grad_elem(_index);
241 8581032 : deviatoric_vector_neighbor(1) = v_grad_neighbor(_index);
242 8581032 : trace_elem += v_grad_elem(1);
243 8581032 : trace_neighbor += v_grad_neighbor(1);
244 8581032 : 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 8581032 : deviatoric_vector_elem(_index) -= 2 / 3 * trace_elem;
256 8581032 : deviatoric_vector_neighbor(_index) -= 2 / 3 * trace_neighbor;
257 :
258 8581032 : _stress_rhs_contribution += _mu(face_arg, state_arg) *
259 : (interp_coeffs.first * deviatoric_vector_elem +
260 : interp_coeffs.second * deviatoric_vector_neighbor) *
261 8581032 : _current_face_info->normal();
262 : }
263 51611281 : _cached_rhs_contribution = true;
264 : }
265 :
266 103222562 : return _stress_rhs_contribution;
267 : }
268 :
269 : Real
270 6357018 : LinearWCNSFVMomentumFlux::computeStressBoundaryMatrixContribution(
271 : const LinearFVAdvectionDiffusionBC * bc)
272 : {
273 6357018 : auto grad_contrib = bc->computeBoundaryGradientMatrixContribution();
274 : // If the boundary condition does not include the diffusivity contribution then
275 : // add it here.
276 6357018 : if (!bc->includesMaterialPropertyMultiplier())
277 : {
278 6357018 : const auto face_arg = singleSidedFaceArg(_current_face_info);
279 6357018 : grad_contrib *= _mu(face_arg, determineState());
280 : }
281 :
282 6357018 : return grad_contrib;
283 : }
284 :
285 : Real
286 6357018 : LinearWCNSFVMomentumFlux::computeStressBoundaryRHSContribution(
287 : const LinearFVAdvectionDiffusionBC * bc)
288 : {
289 6357018 : const auto face_arg = singleSidedFaceArg(_current_face_info);
290 6357018 : auto grad_contrib = bc->computeBoundaryGradientRHSContribution();
291 :
292 : // If the boundary condition does not include the diffusivity contribution then
293 : // add it here.
294 6357018 : if (!bc->includesMaterialPropertyMultiplier())
295 6357018 : 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 6357018 : 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 547056 : const auto elem_info = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
304 547056 : ? _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 547056 : const auto e_Cf = _current_face_info->faceCentroid() - elem_info->centroid();
310 : const auto correction_vector =
311 547056 : _current_face_info->normal() - 1 / (_current_face_info->normal() * e_Cf) * e_Cf;
312 :
313 547056 : grad_contrib += _mu(face_arg, determineState()) * _var.gradSln(*elem_info) *
314 : _boundary_normal_factor * correction_vector;
315 : }
316 :
317 6357018 : 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 1204240 : const auto elem_info = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
322 1204240 : ? _current_face_info->elemInfo()
323 16072 : : _current_face_info->neighborInfo();
324 :
325 1204240 : const auto u_grad_elem = _u_var->gradSln(*elem_info);
326 :
327 : Real trace_elem = 0;
328 : RealVectorValue frace_grad_approx;
329 :
330 1204240 : trace_elem += u_grad_elem(0);
331 1204240 : frace_grad_approx(0) = u_grad_elem(_index);
332 :
333 1204240 : const auto state_arg = determineState();
334 :
335 1204240 : if (_dim > 1)
336 : {
337 1204240 : const auto v_grad_elem = _v_var->gradSln(*elem_info);
338 1204240 : trace_elem += v_grad_elem(1);
339 1204240 : frace_grad_approx(1) = v_grad_elem(_index);
340 :
341 1204240 : 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 1204240 : 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 1204240 : grad_contrib += _mu(face_arg, state_arg) * frace_grad_approx * _boundary_normal_factor *
353 1204240 : _current_face_info->normal();
354 : }
355 :
356 6357018 : return grad_contrib;
357 : }
358 :
359 : Real
360 6357018 : LinearWCNSFVMomentumFlux::computeAdvectionBoundaryMatrixContribution(
361 : const LinearFVAdvectionDiffusionBC * bc)
362 : {
363 6357018 : const auto boundary_value_matrix_contrib = bc->computeBoundaryValueMatrixContribution();
364 6357018 : return boundary_value_matrix_contrib * _face_mass_flux;
365 : }
366 :
367 : Real
368 6357018 : LinearWCNSFVMomentumFlux::computeAdvectionBoundaryRHSContribution(
369 : const LinearFVAdvectionDiffusionBC * bc)
370 : {
371 6357018 : const auto boundary_value_rhs_contrib = bc->computeBoundaryValueRHSContribution();
372 6357018 : return -boundary_value_rhs_contrib * _face_mass_flux;
373 : }
374 :
375 : void
376 57968299 : LinearWCNSFVMomentumFlux::setupFaceData(const FaceInfo * face_info)
377 : {
378 57968299 : 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 57968299 : _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 57968299 : _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 57968299 : 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 57968299 : _stress_rhs_contribution = 0;
396 57968299 : }
|
