3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
freeflow/navierstokes/staggered/fluxvariables.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH
25#define DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH
26
27#include <array>
28#include <dune/common/std/optional.hh>
29
30#include <dumux/common/math.hh>
34
37
38
40#include "velocitygradients.hh"
41
42namespace Dumux {
43
44// forward declaration
45template<class TypeTag, DiscretizationMethod discMethod>
46class NavierStokesFluxVariablesImpl;
47
52template<class TypeTag>
54: public FluxVariablesBase<GetPropType<TypeTag, Properties::Problem>,
55 typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView,
56 typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView,
57 typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView>
58{
60
61 using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
62 using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
63 using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
64
65 using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache;
66 using FluxVariablesCache = typename GridFluxVariablesCache::FluxVariablesCache;
67
68 using GridFaceVariables = typename GridVariables::GridFaceVariables;
69 using ElementFaceVariables = typename GridFaceVariables::LocalView;
70 using FaceVariables = typename GridFaceVariables::FaceVariables;
71
74 using FVElementGeometry = typename GridGeometry::LocalView;
75 using GridView = typename GridGeometry::GridView;
76 using Element = typename GridView::template Codim<0>::Entity;
79 using Indices = typename ModelTraits::Indices;
80 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
81 using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
85
86 static constexpr bool normalizePressure = getPropValue<TypeTag, Properties::NormalizePressure>();
87
88 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
89
90 static constexpr auto cellCenterIdx = GridGeometry::cellCenterIdx();
91 static constexpr auto faceIdx = GridGeometry::faceIdx();
92
93 static constexpr int upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
94 static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
95
96public:
97
99
108 template<class UpwindTerm>
109 static Scalar advectiveFluxForCellCenter(const Problem& problem,
110 const ElementVolumeVariables& elemVolVars,
111 const ElementFaceVariables& elemFaceVars,
112 const SubControlVolumeFace &scvf,
113 UpwindTerm upwindTerm)
114 {
115 const Scalar velocity = elemFaceVars[scvf].velocitySelf();
116 const bool insideIsUpstream = scvf.directionSign() == sign(velocity);
117 static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight");
118
119 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
120 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
121
122 const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
123 const auto& downstreamVolVars = insideIsUpstream ? outsideVolVars : insideVolVars;
124
125 const Scalar flux = (upwindWeight * upwindTerm(upstreamVolVars) +
126 (1.0 - upwindWeight) * upwindTerm(downstreamVolVars))
127 * velocity * scvf.area() * scvf.directionSign();
128
129 return flux * extrusionFactor_(elemVolVars, scvf);
130 }
131
147 CellCenterPrimaryVariables computeMassFlux(const Problem& problem,
148 const Element& element,
149 const FVElementGeometry& fvGeometry,
150 const ElementVolumeVariables& elemVolVars,
151 const ElementFaceVariables& elemFaceVars,
152 const SubControlVolumeFace& scvf,
153 const FluxVariablesCache& fluxVarsCache)
154 {
155 // The advectively transported quantity (i.e density for a single-phase system).
156 auto upwindTerm = [](const auto& volVars) { return volVars.density(); };
157
158 // Call the generic flux function.
159 CellCenterPrimaryVariables result(0.0);
160 result[Indices::conti0EqIdx - ModelTraits::dim()] = advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTerm);
161
162 return result;
163 }
164
168 FacePrimaryVariables computeMomentumFlux(const Problem& problem,
169 const Element& element,
170 const SubControlVolumeFace& scvf,
171 const FVElementGeometry& fvGeometry,
172 const ElementVolumeVariables& elemVolVars,
173 const ElementFaceVariables& elemFaceVars,
174 const GridFluxVariablesCache& gridFluxVarsCache)
175 {
176 return computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache) +
177 computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache);
178 }
179
197 FacePrimaryVariables computeFrontalMomentumFlux(const Problem& problem,
198 const Element& element,
199 const SubControlVolumeFace& scvf,
200 const FVElementGeometry& fvGeometry,
201 const ElementVolumeVariables& elemVolVars,
202 const ElementFaceVariables& elemFaceVars,
203 const GridFluxVariablesCache& gridFluxVarsCache)
204 {
205 FacePrimaryVariables frontalFlux(0.0);
206
207 // The velocities of the dof at interest and the one of the opposite scvf.
208 const Scalar velocitySelf = elemFaceVars[scvf].velocitySelf();
209 const Scalar velocityOpposite = elemFaceVars[scvf].velocityOpposite();
210
211 // Advective flux.
212 if (problem.enableInertiaTerms())
213 {
214 // Get the average velocity at the center of the element (i.e. the location of the staggered face).
215 const Scalar transportingVelocity = (velocitySelf + velocityOpposite) * 0.5;
216
217 frontalFlux += StaggeredUpwindFluxVariables<TypeTag, upwindSchemeOrder>::computeUpwindedFrontalMomentum(scvf, elemFaceVars, elemVolVars, gridFluxVarsCache, transportingVelocity)
218 * transportingVelocity * -1.0 * scvf.directionSign();
219 }
220
221 // The volume variables within the current element. We only require those (and none of neighboring elements)
222 // because the fluxes are calculated over the staggered face at the center of the element.
223 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
224
225 // Diffusive flux.
226 const Scalar velocityGrad_ii = VelocityGradients::velocityGradII(scvf, elemFaceVars[scvf]) * scvf.directionSign();
227
228 static const bool enableUnsymmetrizedVelocityGradient
229 = getParamFromGroup<bool>(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradient", false);
230 const Scalar factor = enableUnsymmetrizedVelocityGradient ? 1.0 : 2.0;
231 frontalFlux -= factor * insideVolVars.effectiveViscosity() * velocityGrad_ii;
232
233 // The pressure term.
234 // If specified, the pressure can be normalized using the initial value on the scfv of interest.
235 // The scvf is used to normalize by the same value from the left and right side.
236 // Can potentially help to improve the condition number of the system matrix.
237 const Scalar pressure = normalizePressure ?
238 insideVolVars.pressure() - problem.initial(scvf)[Indices::pressureIdx]
239 : insideVolVars.pressure();
240
241 // Account for the orientation of the staggered face's normal outer normal vector
242 // (pointing in opposite direction of the scvf's one).
243 frontalFlux += pressure * -1.0 * scvf.directionSign();
244
245 // Account for the staggered face's area. For rectangular elements, this equals the area of the scvf
246 // our velocity dof of interest lives on.
247 return frontalFlux * scvf.area() * insideVolVars.extrusionFactor();
248 }
249
267 FacePrimaryVariables computeLateralMomentumFlux(const Problem& problem,
268 const Element& element,
269 const SubControlVolumeFace& scvf,
270 const FVElementGeometry& fvGeometry,
271 const ElementVolumeVariables& elemVolVars,
272 const ElementFaceVariables& elemFaceVars,
273 const GridFluxVariablesCache& gridFluxVarsCache)
274 {
275 FacePrimaryVariables lateralFlux(0.0);
276 const auto& faceVars = elemFaceVars[scvf];
277 const std::size_t numSubFaces = scvf.pairData().size();
278
279 // If the current scvf is on a boundary, check if there is a Neumann BC for the stress in tangential direction.
280 // Create a boundaryTypes object (will be empty if not at a boundary).
281 Dune::Std::optional<BoundaryTypes> currentScvfBoundaryTypes;
282 if (scvf.boundary())
283 currentScvfBoundaryTypes.emplace(problem.boundaryTypes(element, scvf));
284
285 // Account for all sub faces.
286 for (int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx)
287 {
288 const auto eIdx = scvf.insideScvIdx();
289 // Get the face normal to the face the dof lives on. The staggered sub face conincides with half of this lateral face.
290 const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
291
292 // Create a boundaryTypes object (will be empty if not at a boundary).
293 Dune::Std::optional<BoundaryTypes> lateralFaceBoundaryTypes;
294
295 // Check if there is face/element parallel to our face of interest where the dof lives on. If there is no parallel neighbor,
296 // we are on a boundary where we have to check for boundary conditions.
297 if (lateralScvf.boundary())
298 {
299 // Retrieve the boundary types that correspond to the center of the lateral scvf. As a convention, we always query
300 // the type of BCs at the center of the element's "actual" lateral scvf (not the face of the staggered control volume).
301 // The value of the BC will be evaluated at the center of the staggered face.
302 // --------T######V || frontal face of staggered half-control-volume
303 // | || | current scvf # lateral staggered face of interest (may lie on a boundary)
304 // | || | x dof position
305 // | || x~~~~> vel.Self -- element boundaries
306 // | || | T position at which the type of boundary conditions will be evaluated
307 // | || | (center of lateral scvf)
308 // ---------------- V position at which the value of the boundary conditions will be evaluated
309 // (center of the staggered lateral face)
310 lateralFaceBoundaryTypes.emplace(problem.boundaryTypes(element, lateralScvf));
311 }
312
313 // If the current scvf is on a bounary and if there is a Neumann or Beavers-Joseph-(Saffmann) BC for the stress in tangential direction,
314 // assign this value for the lateral momentum flux. No further calculations are required. We assume that all lateral faces
315 // have the same type of BC (Neumann or Beavers-Joseph-(Saffmann)), but we sample the value at their actual positions.
316 if (currentScvfBoundaryTypes)
317 {
318 // Handle Neumann BCs.
319 if (currentScvfBoundaryTypes->isNeumann(Indices::velocity(lateralScvf.directionIndex())))
320 {
321 // Get the location of the lateral staggered face's center.
322 const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
323 lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars,
324 scvf.makeBoundaryFace(lateralStaggeredFaceCenter))[Indices::velocity(lateralScvf.directionIndex())]
325 * extrusionFactor_(elemVolVars, lateralScvf) * lateralScvf.area() * 0.5
326 * lateralScvf.directionSign();
327 continue;
328 }
329
330 // Treat the edge case when our current scvf has a BJ condition while the lateral face has a Neumann BC.
331 // It is not clear how to evaluate the BJ condition here.
332 // For symmetry reasons, our own scvf should then have the same Neumann flux as the lateral face.
333 // TODO: We should clarify if this is the correct approach.
334 if (currentScvfBoundaryTypes->isBJS(Indices::velocity(lateralScvf.directionIndex())) && lateralFaceBoundaryTypes &&
335 lateralFaceBoundaryTypes->isNeumann(Indices::velocity(scvf.directionIndex())))
336 {
337 const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
338 lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars,
339 lateralScvf.makeBoundaryFace(lateralStaggeredFaceCenter))[Indices::velocity(scvf.directionIndex())]
340 * extrusionFactor_(elemVolVars, lateralScvf) * lateralScvf.area() * 0.5
341 * lateralScvf.directionSign();
342 continue;
343 }
344 }
345
346 // Check if the lateral face (perpendicular to our current scvf) lies on a boundary. If yes, boundary conditions might need to be treated
347 // and further calculations can be skipped.
348 if (lateralScvf.boundary())
349 {
350 // Check if we have a symmetry boundary condition. If yes, the tangential part of the momentum flux can be neglected
351 // and we may skip any further calculations for the given sub face.
352 if (lateralFaceBoundaryTypes->isSymmetry())
353 continue;
354
355 // Handle Neumann boundary conditions. No further calculations are then required for the given sub face.
356 if (lateralFaceBoundaryTypes->isNeumann(Indices::velocity(scvf.directionIndex())))
357 {
358 // Construct a temporary scvf which corresponds to the staggered sub face, featuring the location
359 // the staggered faces's center.
360 const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
361 const auto lateralBoundaryFace = lateralScvf.makeBoundaryFace(lateralStaggeredFaceCenter);
362 lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())]
363 * elemVolVars[lateralScvf.insideScvIdx()].extrusionFactor() * lateralScvf.area() * 0.5;
364 continue;
365 }
366
367 // Handle wall-function fluxes (only required for RANS models)
368 if (incorporateWallFunction_(lateralFlux, problem, element, fvGeometry, scvf, elemVolVars, elemFaceVars, localSubFaceIdx))
369 continue;
370 }
371
372 // Check the consistency of the boundary conditions, exactly one of the following must be set
373 if (lateralFaceBoundaryTypes)
374 {
375 std::bitset<3> admittableBcTypes;
376 admittableBcTypes.set(0, lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx));
377 admittableBcTypes.set(1, lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex())));
378 admittableBcTypes.set(2, lateralFaceBoundaryTypes->isBJS(Indices::velocity(scvf.directionIndex())));
379 if (admittableBcTypes.count() != 1)
380 {
381 DUNE_THROW(Dune::InvalidStateException, "Invalid boundary conditions for lateral scvf "
382 "for the momentum equations at global position " << lateralStaggeredFaceCenter_(scvf, localSubFaceIdx)
383 << ", current scvf global position " << scvf.center());
384 }
385 }
386
387 // If none of the above boundary conditions apply for the given sub face, proceed to calculate the tangential momentum flux.
388 if (problem.enableInertiaTerms())
389 lateralFlux += computeAdvectivePartOfLateralMomentumFlux_(problem, fvGeometry, element,
390 scvf, elemVolVars, faceVars,
391 gridFluxVarsCache,
392 currentScvfBoundaryTypes, lateralFaceBoundaryTypes,
393 localSubFaceIdx);
394
395 lateralFlux += computeDiffusivePartOfLateralMomentumFlux_(problem, fvGeometry, element,
396 scvf, elemVolVars, faceVars,
397 currentScvfBoundaryTypes, lateralFaceBoundaryTypes,
398 localSubFaceIdx);
399 }
400 return lateralFlux;
401 }
402
419 FacePrimaryVariables inflowOutflowBoundaryFlux(const Problem& problem,
420 const Element& element,
421 const SubControlVolumeFace& scvf,
422 const ElementVolumeVariables& elemVolVars,
423 const ElementFaceVariables& elemFaceVars) const
424 {
425 FacePrimaryVariables inOrOutflow(0.0);
426 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
427
428 // Advective momentum flux.
429 if (problem.enableInertiaTerms())
430 {
431 const Scalar velocitySelf = elemFaceVars[scvf].velocitySelf();
432 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
433 const auto& upVolVars = (scvf.directionSign() == sign(velocitySelf)) ?
434 insideVolVars : outsideVolVars;
435
436 inOrOutflow += velocitySelf * velocitySelf * upVolVars.density();
437 }
438
439 // Apply a pressure at the boundary.
440 const Scalar boundaryPressure = normalizePressure
441 ? (problem.dirichlet(element, scvf)[Indices::pressureIdx] -
442 problem.initial(scvf)[Indices::pressureIdx])
443 : problem.dirichlet(element, scvf)[Indices::pressureIdx];
444 inOrOutflow += boundaryPressure;
445
446 // Account for the orientation of the face at the boundary,
447 return inOrOutflow * scvf.directionSign() * scvf.area() * insideVolVars.extrusionFactor();
448 }
449
450private:
451
473 FacePrimaryVariables computeAdvectivePartOfLateralMomentumFlux_(const Problem& problem,
474 const FVElementGeometry& fvGeometry,
475 const Element& element,
476 const SubControlVolumeFace& scvf,
477 const ElementVolumeVariables& elemVolVars,
478 const FaceVariables& faceVars,
479 const GridFluxVariablesCache& gridFluxVarsCache,
480 const Dune::Std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
481 const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
482 const int localSubFaceIdx)
483 {
484 const auto eIdx = scvf.insideScvIdx();
485 const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
486
487 // Get the transporting velocity, located at the scvf perpendicular to the current scvf where the dof
488 // of interest is located.
489 const Scalar transportingVelocity = [&]()
490 {
491 if (!scvf.boundary())
492 return faceVars.velocityLateralInside(localSubFaceIdx);
493 else
494 {
495 // Create a boundaryTypes object. Get the boundary conditions. We sample the type of BC at the center of the current scvf.
496 const BoundaryTypes bcTypes = problem.boundaryTypes(element, scvf);
497
498 if (bcTypes.isDirichlet(Indices::velocity(lateralFace.directionIndex())))
499 {
500 // Construct a temporary scvf which corresponds to the staggered sub face, featuring the location
501 // the staggered faces's center.
502 const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
503 return problem.dirichlet(element, scvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(lateralFace.directionIndex())];
504 }
505 else if (bcTypes.isBJS(Indices::velocity(lateralFace.directionIndex())))
506 {
507 return VelocityGradients::beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars,
508 currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
509 }
510 else
511 return faceVars.velocityLateralInside(localSubFaceIdx);
512 }
513 }();
514
515 return StaggeredUpwindFluxVariables<TypeTag, upwindSchemeOrder>::computeUpwindedLateralMomentum(problem, fvGeometry, element, scvf, elemVolVars, faceVars,
516 gridFluxVarsCache, localSubFaceIdx, currentScvfBoundaryTypes, lateralFaceBoundaryTypes)
517 * transportingVelocity * lateralFace.directionSign() * lateralFace.area() * 0.5 * extrusionFactor_(elemVolVars, lateralFace);
518 }
519
542 FacePrimaryVariables computeDiffusivePartOfLateralMomentumFlux_(const Problem& problem,
543 const FVElementGeometry& fvGeometry,
544 const Element& element,
545 const SubControlVolumeFace& scvf,
546 const ElementVolumeVariables& elemVolVars,
547 const FaceVariables& faceVars,
548 const Dune::Std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
549 const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
550 const int localSubFaceIdx)
551 {
552 const auto eIdx = scvf.insideScvIdx();
553 const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
554
555 FacePrimaryVariables lateralDiffusiveFlux(0.0);
556
557 static const bool enableUnsymmetrizedVelocityGradient
558 = getParamFromGroup<bool>(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradient", false);
559
560 // Get the volume variables of the own and the neighboring element. The neighboring
561 // element is adjacent to the staggered face normal to the current scvf
562 // where the dof of interest is located.
563 const auto& insideVolVars = elemVolVars[lateralFace.insideScvIdx()];
564 const auto& outsideVolVars = elemVolVars[lateralFace.outsideScvIdx()];
565
566 // Get the averaged viscosity at the staggered face normal to the current scvf.
567 const Scalar muAvg = lateralFace.boundary()
568 ? insideVolVars.effectiveViscosity()
569 : (insideVolVars.effectiveViscosity() + outsideVolVars.effectiveViscosity()) * 0.5;
570
571 // Consider the shear stress caused by the gradient of the velocities normal to our face of interest.
572 if (!enableUnsymmetrizedVelocityGradient)
573 {
574 if (!scvf.boundary() ||
575 currentScvfBoundaryTypes->isDirichlet(Indices::velocity(lateralFace.directionIndex())) ||
576 currentScvfBoundaryTypes->isBJS(Indices::velocity(lateralFace.directionIndex())))
577 {
578 const Scalar velocityGrad_ji = VelocityGradients::velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
579 // Account for the orientation of the staggered normal face's outer normal vector.
580 lateralDiffusiveFlux -= muAvg * velocityGrad_ji * lateralFace.directionSign();
581 }
582 }
583
584 // Consider the shear stress caused by the gradient of the velocities parallel to our face of interest.
585 // If we have a Dirichlet condition for the pressure at the lateral face we assume to have a zero velocityGrad_ij velocity gradient
586 // so we can skip the computation.
587 if (!lateralFace.boundary() || !lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx))
588 {
589 const Scalar velocityGrad_ij = VelocityGradients::velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
590 lateralDiffusiveFlux -= muAvg * velocityGrad_ij * lateralFace.directionSign();
591 }
592
593 // Account for the area of the staggered lateral face (0.5 of the coinciding scfv).
594 return lateralDiffusiveFlux * lateralFace.area() * 0.5 * extrusionFactor_(elemVolVars, lateralFace);
595 }
596
611 const GlobalPosition& lateralStaggeredFaceCenter_(const SubControlVolumeFace& scvf, const int localSubFaceIdx) const
612 {
613 return scvf.pairData(localSubFaceIdx).lateralStaggeredFaceCenter;
614 };
615
617 static Scalar extrusionFactor_(const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf)
618 {
619 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
620 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
621 return harmonicMean(insideVolVars.extrusionFactor(), outsideVolVars.extrusionFactor());
622 }
623
625 template<class ...Args, bool turbulenceModel = ModelTraits::usesTurbulenceModel(), std::enable_if_t<!turbulenceModel, int> = 0>
626 bool incorporateWallFunction_(Args&&... args) const
627 { return false; }
628
630 template<bool turbulenceModel = ModelTraits::usesTurbulenceModel(), std::enable_if_t<turbulenceModel, int> = 0>
631 bool incorporateWallFunction_(FacePrimaryVariables& lateralFlux,
632 const Problem& problem,
633 const Element& element,
634 const FVElementGeometry& fvGeometry,
635 const SubControlVolumeFace& scvf,
636 const ElementVolumeVariables& elemVolVars,
637 const ElementFaceVariables& elemFaceVars,
638 const std::size_t localSubFaceIdx) const
639 {
640 const auto eIdx = scvf.insideScvIdx();
641 const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
642
643 if (problem.useWallFunction(element, lateralScvf, Indices::velocity(scvf.directionIndex())))
644 {
645 const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
646 const auto lateralBoundaryFace = lateralScvf.makeBoundaryFace(lateralStaggeredFaceCenter);
647 lateralFlux += problem.wallFunction(element, fvGeometry, elemVolVars, elemFaceVars, scvf, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())]
648 * extrusionFactor_(elemVolVars, lateralScvf) * lateralScvf.area() * 0.5;
649 return true;
650 }
651 else
652 return false;
653 }
654};
655
656} // end namespace Dumux
657
658#endif // DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH
Some exceptions thrown in DuMux
Define some often used mathematical functions.
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
The available discretization methods in Dumux.
Base class for the flux variables living on a sub control volume face.
DiscretizationMethod
The available discretization methods in Dumux.
Definition: method.hh:37
constexpr Scalar harmonicMean(Scalar x, Scalar y, Scalar wx=1.0, Scalar wy=1.0) noexcept
Calculate the (weighted) harmonic mean of two scalar values.
Definition: math.hh:68
constexpr int sign(const ValueType &value) noexcept
Sign or signum function.
Definition: math.hh:618
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Property tag Indices
Definition: porousmediumflow/sequential/properties.hh:59
std::string pressure(int phaseIdx) noexcept
I/O name of pressure for multiphase systems.
Definition: name.hh:34
Class to specify the type of a boundary.
Definition: common/boundarytypes.hh:38
bool isDirichlet(unsigned eqIdx) const
Returns true if an equation is used to specify a Dirichlet condition.
Definition: common/boundarytypes.hh:236
Base class for the flux variables living on a sub control volume face.
Definition: fluxvariablesbase.hh:45
Definition: freeflow/navierstokes/fluxvariables.hh:35
FacePrimaryVariables computeLateralMomentumFlux(const Problem &problem, const Element &element, const SubControlVolumeFace &scvf, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const ElementFaceVariables &elemFaceVars, const GridFluxVariablesCache &gridFluxVarsCache)
Returns the momentum flux over the staggered faces perpendicular to the scvf where the velocity dof o...
Definition: freeflow/navierstokes/staggered/fluxvariables.hh:267
FacePrimaryVariables computeFrontalMomentumFlux(const Problem &problem, const Element &element, const SubControlVolumeFace &scvf, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const ElementFaceVariables &elemFaceVars, const GridFluxVariablesCache &gridFluxVarsCache)
Returns the frontal part of the momentum flux. This treats the flux over the staggered face at the ce...
Definition: freeflow/navierstokes/staggered/fluxvariables.hh:197
static Scalar advectiveFluxForCellCenter(const Problem &problem, const ElementVolumeVariables &elemVolVars, const ElementFaceVariables &elemFaceVars, const SubControlVolumeFace &scvf, UpwindTerm upwindTerm)
Returns the advective flux over a sub control volume face.
Definition: freeflow/navierstokes/staggered/fluxvariables.hh:109
FacePrimaryVariables computeMomentumFlux(const Problem &problem, const Element &element, const SubControlVolumeFace &scvf, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const ElementFaceVariables &elemFaceVars, const GridFluxVariablesCache &gridFluxVarsCache)
Returns the momentum flux over all staggered faces.
Definition: freeflow/navierstokes/staggered/fluxvariables.hh:168
GetPropType< TypeTag, Properties::HeatConductionType > HeatConductionType
Definition: freeflow/navierstokes/staggered/fluxvariables.hh:98
FacePrimaryVariables inflowOutflowBoundaryFlux(const Problem &problem, const Element &element, const SubControlVolumeFace &scvf, const ElementVolumeVariables &elemVolVars, const ElementFaceVariables &elemFaceVars) const
Returns the momentum flux over an inflow or outflow boundary face.
Definition: freeflow/navierstokes/staggered/fluxvariables.hh:419
CellCenterPrimaryVariables computeMassFlux(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const ElementFaceVariables &elemFaceVars, const SubControlVolumeFace &scvf, const FluxVariablesCache &fluxVarsCache)
Computes the flux for the cell center residual (mass balance).
Definition: freeflow/navierstokes/staggered/fluxvariables.hh:147
static FacePrimaryVariables computeUpwindedFrontalMomentum(const SubControlVolumeFace &scvf, const ElementFaceVariables &elemFaceVars, const ElementVolumeVariables &elemVolVars, const GridFluxVariablesCache &gridFluxVarsCache, const Scalar transportingVelocity)
Returns the momentum in the frontal directon.
Definition: staggeredupwindfluxvariables.hh:87
static FacePrimaryVariables computeUpwindedLateralMomentum(const Problem &problem, const FVElementGeometry &fvGeometry, const Element &element, const SubControlVolumeFace &scvf, const ElementVolumeVariables &elemVolVars, const FaceVariables &faceVars, const GridFluxVariablesCache &gridFluxVarsCache, const int localSubFaceIdx, const Dune::Std::optional< BoundaryTypes > &currentScvfBoundaryTypes, const Dune::Std::optional< BoundaryTypes > &lateralFaceBoundaryTypes)
Returns the momentum in the lateral directon.
Definition: staggeredupwindfluxvariables.hh:118
Helper class for calculating the velocity gradients for the Navier-Stokes model using the staggered g...
Definition: velocitygradients.hh:39
Declares all properties used in Dumux.