24#ifndef DUMUX_NAVIERSTOKES_MOMENTUM_FLUXVARIABLES_HH
25#define DUMUX_NAVIERSTOKES_MOMENTUM_FLUXVARIABLES_HH
47template<
class TypeTag>
52 using GridVolumeVariables =
typename GridVariables::GridVolumeVariables;
53 using ElementVolumeVariables =
typename GridVolumeVariables::LocalView;
54 using VolumeVariables =
typename GridVolumeVariables::VolumeVariables;
56 using GridFluxVariablesCache =
typename GridVariables::GridFluxVariablesCache;
57 using ElementFluxVariablesCache =
typename GridFluxVariablesCache::LocalView;
58 using FluxVariablesCache =
typename GridFluxVariablesCache::FluxVariablesCache;
60 using GridGeometry =
typename GridVariables::GridGeometry;
61 using FVElementGeometry =
typename GridGeometry::LocalView;
62 using SubControlVolume =
typename FVElementGeometry::SubControlVolume;
63 using SubControlVolumeFace =
typename FVElementGeometry::SubControlVolumeFace;
68 using GridView =
typename GridGeometry::GridView;
69 using Element =
typename GridView::template Codim<0>::Entity;
72 using Indices =
typename ModelTraits::Indices;
75 using GlobalPosition =
typename Element::Geometry::GlobalCoordinate;
85 const SubControlVolumeFace&
scvFace,
99 ==
static_cast<std::size_t
>(GridView::dimension),
100 "Expects problem.dirichlet to return an array with as many entries as dimensions."
105 {
return *problemPtr_; }
108 {
return *elementPtr_; }
111 {
return *scvFacePtr_; }
114 {
return *fvGeometryPtr_; }
117 {
return *elemVolVarsPtr_; }
120 {
return *elemFluxVarsCachePtr_; }
123 {
return *elemBcTypesPtr_; }
130 if (!this->
problem().enableInertiaTerms())
131 return NumEqVector(0.0);
133 if (this->
scvFace().isFrontal())
144 if (this->
scvFace().isFrontal())
170 const auto& scvf = this->
scvFace();
171 assert(scvf.isFrontal());
173 NumEqVector result(0.0);
183 GlobalPosition gradVn(0.0);
184 gradV.mv(scvf.unitOuterNormal(), gradVn);
185 const auto& scv =
fvGeometry.scv(scvf.insideScvIdx());
186 const Scalar velocityGrad_ii = gradVn[scv.dofAxis()];
188 static const bool enableUnsymmetrizedVelocityGradient
189 = getParamFromGroup<bool>(this->
problem().paramGroup(),
"FreeFlow.EnableUnsymmetrizedVelocityGradient",
false);
190 const Scalar factor = enableUnsymmetrizedVelocityGradient ? 1.0 : 2.0;
193 result -= factor * mu * velocityGrad_ii * Extrusion::area(scvf) *
elemVolVars[scvf.insideScvIdx()].extrusionFactor();
195 static const bool enableDilatationTerm = getParamFromGroup<bool>(this->
problem().paramGroup(),
"FreeFlow.EnableDilatationTerm",
false);
196 if (enableDilatationTerm)
198 Scalar divergence = 0.0;
201 const auto frontalScvf = *(scvfs(
fvGeometry, scv).begin());
202 assert(frontalScvf.isFrontal() && !frontalScvf.boundary());
207 result += 2.0/3.0 * mu * divergence * scvf.directionSign() * Extrusion::area(scvf) *
elemVolVars[scvf.insideScvIdx()].extrusionFactor();
238 const auto& scvf = this->
scvFace();
239 assert(scvf.isLateral());
241 NumEqVector result(0.0);
245 const auto& scv =
fvGeometry.scv(scvf.insideScvIdx());
247 static const bool enableUnsymmetrizedVelocityGradient
248 = getParamFromGroup<bool>(
problem.paramGroup(),
"FreeFlow.EnableUnsymmetrizedVelocityGradient",
false);
256 GlobalPosition gradVn(0.0);
257 gradV.mv(scvf.unitOuterNormal(), gradVn);
258 const Scalar velocityGrad_ij = gradVn[scv.dofAxis()];
259 result -= mu * velocityGrad_ij;
262 if (!enableUnsymmetrizedVelocityGradient)
264 GlobalPosition gradVTransposedN(0.0);
265 gradV.mtv(scvf.unitOuterNormal(), gradVTransposedN);
266 const Scalar velocityGrad_ji = gradVTransposedN[scv.dofAxis()];
267 result -= mu * velocityGrad_ji;
271 return result * Extrusion::area(scvf) *
elemVolVars[scvf.insideScvIdx()].extrusionFactor();
292 NumEqVector result(0.0);
293 const auto& scvf = this->
scvFace();
294 if (scvf.isLateral() || scvf.boundary())
299 result =
pressure*Extrusion::area(scvf)*this->
elemVolVars()[scvf.insideScvIdx()].extrusionFactor();
310 result -= referencePressure*scvf.area();
313 result *= scvf.directionSign();
337 NumEqVector flux(0.0);
338 const auto& scvf = this->
scvFace();
339 assert(scvf.isFrontal());
343 const auto velocitySelf =
elemVolVars[scvf.insideScvIdx()].velocity();
344 const auto velocityOpposite =
elemVolVars[scvf.outsideScvIdx()].velocity();
347 const Scalar transportingVelocity = (velocitySelf + velocityOpposite) * 0.5;
349 const bool selfIsUpstream = scvf.directionSign() ==
sign(transportingVelocity);
352 static const auto upwindWeight = getParamFromGroup<Scalar>(
problem.paramGroup(),
"Flux.UpwindWeight");
353 const Scalar transportedMomentum = selfIsUpstream ? (upwindWeight * velocitySelf + (1.0 - upwindWeight) * velocityOpposite) *
density
354 : (upwindWeight * velocityOpposite + (1.0 - upwindWeight) * velocitySelf) *
density;
356 return transportingVelocity * transportedMomentum * scvf.directionSign() * Extrusion::area(scvf) * extrusionFactor_(
elemVolVars, scvf);
383 NumEqVector flux(0.0);
384 const auto& scvf = this->
scvFace();
385 assert(scvf.isLateral());
392 const Scalar transportingVelocity = [&]()
394 const auto& orthogonalScvf =
fvGeometry.lateralOrthogonalScvf(scvf);
395 const Scalar innerTransportingVelocity =
elemVolVars[orthogonalScvf.insideScvIdx()].velocity();
397 static const bool useOldScheme = getParam<bool>(
"FreeFlow.UseOldTransportingVelocity",
true);
401 if (scvf.boundary() &&
fvGeometry.scv(scvf.insideScvIdx()).boundary())
403 if (this->
elemBcTypes()[scvf.localIndex()].isDirichlet(scvf.normalAxis()))
404 return problem.dirichlet(this->
element(), scvf)[scvf.normalAxis()];
407 return innerTransportingVelocity;
413 if (this->
elemBcTypes()[scvf.localIndex()].isDirichlet(scvf.normalAxis()))
414 return 0.5*(
problem.dirichlet(this->
element(), scvf)[scvf.normalAxis()] + innerTransportingVelocity);
417 if (orthogonalScvf.boundary())
419 if (this->
elemBcTypes()[orthogonalScvf.localIndex()].isDirichlet(scvf.normalAxis()))
420 return 0.5*(
problem.dirichlet(this->element(), scvf)[scvf.normalAxis()] + innerTransportingVelocity);
422 return innerTransportingVelocity;
426 const auto insideVolume =
fvGeometry.scv(orthogonalScvf.insideScvIdx()).
volume();
427 const auto outsideVolume =
fvGeometry.scv(orthogonalScvf.outsideScvIdx()).
volume();
428 const auto outerTransportingVelocity =
elemVolVars[orthogonalScvf.outsideScvIdx()].velocity();
429 return (insideVolume*innerTransportingVelocity + outsideVolume*outerTransportingVelocity) / (insideVolume + outsideVolume);
432 const Scalar transportedMomentum = [&]()
434 const auto& insideScv =
fvGeometry.scv(scvf.insideScvIdx());
436 auto getDirichletMomentumFlux = [&]()
444 if (!this->
elemBcTypes()[scvf.localIndex()].isDirichlet(insideScv.dofAxis()))
445 DUNE_THROW(Dune::InvalidStateException,
"Neither Dirichlet nor Neumann BC set at " << scvf.ipGlobal());
447 return getDirichletMomentumFlux();
451 if (
fvGeometry.scvfIntegrationPointInConcaveCorner(scvf))
454 const auto& outsideScvfWithSameIntegrationPoint =
fvGeometry.outsideScvfWithSameIntegrationPoint(scvf);
455 if (!this->
problem().boundaryTypes(this->
element(), outsideScvfWithSameIntegrationPoint).isDirichlet(insideScv.dofAxis()))
456 DUNE_THROW(Dune::InvalidStateException,
"Neither Dirichlet nor Neumann BC set at " << outsideScvfWithSameIntegrationPoint.ipGlobal());
458 return getDirichletMomentumFlux();
462 const bool selfIsUpstream = scvf.directionSign() ==
sign(transportingVelocity);
464 const auto innerVelocity =
elemVolVars[scvf.insideScvIdx()].velocity();
465 const auto outerVelocity =
elemVolVars[scvf.outsideScvIdx()].velocity();
469 const auto insideMomentum = innerVelocity * rho.first;
470 const auto outsideMomentum = outerVelocity * rho.second;
473 static const auto upwindWeight = getParamFromGroup<Scalar>(
problem.paramGroup(),
"Flux.UpwindWeight");
475 return selfIsUpstream ? (upwindWeight * insideMomentum + (1.0 - upwindWeight) * outsideMomentum)
476 : (upwindWeight * outsideMomentum + (1.0 - upwindWeight) * insideMomentum);
479 return transportingVelocity * transportedMomentum * scvf.directionSign() * Extrusion::area(scvf) * extrusionFactor_(
elemVolVars, scvf);
484 template<
class ElementVolumeVariables,
class SubControlVolumeFace>
485 Scalar extrusionFactor_(
const ElementVolumeVariables&
elemVolVars,
const SubControlVolumeFace& scvf)
const
487 const auto& insideVolVars =
elemVolVars[scvf.insideScvIdx()];
488 const auto& outsideVolVars =
elemVolVars[scvf.outsideScvIdx()];
489 return harmonicMean(insideVolVars.extrusionFactor(), outsideVolVars.extrusionFactor());
493 const Problem* problemPtr_;
494 const Element* elementPtr_;
495 const FVElementGeometry* fvGeometryPtr_;
496 const SubControlVolumeFace* scvFacePtr_;
497 const ElementVolumeVariables* elemVolVarsPtr_;
498 const ElementFluxVariablesCache* elemFluxVarsCachePtr_;
499 const ElementBoundaryTypes* elemBcTypesPtr_;
Some exceptions thrown in DuMux
A helper to deduce a vector with the same size as numbers of equations.
Define some often used mathematical functions.
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
Helper classes to compute the integration elements.
The available discretization methods in Dumux.
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:69
constexpr int sign(const ValueType &value) noexcept
Sign or signum function.
Definition: math.hh:641
typename NumEqVectorTraits< PrimaryVariables >::type NumEqVector
A vector with the same size as numbers of equations This is the default implementation and has to be ...
Definition: numeqvector.hh:46
typename Extrusion< T >::type Extrusion_t
Convenience alias for obtaining the extrusion type.
Definition: extrusion.hh:177
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:150
std::string pressure(int phaseIdx) noexcept
I/O name of pressure for multiphase systems.
Definition: name.hh:34
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Scalar volume(Shape shape, Scalar inscribedRadius)
Returns the volume of a given geometry based on the inscribed radius.
Definition: poreproperties.hh:73
The flux variables class for the Navier-Stokes model using the staggered grid discretization.
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:49
NumEqVector pressureContribution() const
Returns the frontal pressure contribution.
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:290
NumEqVector diffusiveMomentumFlux() const
Returns the diffusive momentum flux due to viscous forces.
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:142
NumEqVector frontalAdvectiveMomentumFlux() const
Returns the frontal part of the momentum flux. This treats the flux over the staggered face at the ce...
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:335
NumEqVector lateralAdvectiveMomentumFlux() const
Returns the advective momentum flux over the staggered face perpendicular to the scvf where the veloc...
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:381
NumEqVector frontalDiffusiveMomentumFlux() const
Returns the frontal part of the momentum flux. This treats the flux over the staggered face at the ce...
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:168
const ElementBoundaryTypes & elemBcTypes() const
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:122
const Problem & problem() const
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:104
const SubControlVolumeFace & scvFace() const
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:110
const FVElementGeometry & fvGeometry() const
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:113
NumEqVector advectiveMomentumFlux() const
Returns the diffusive momentum flux due to viscous forces.
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:128
NavierStokesMomentumFluxVariables(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const SubControlVolumeFace &scvFace, const ElementVolumeVariables &elemVolVars, const ElementFluxVariablesCache &elemFluxVarsCache, const ElementBoundaryTypes &elemBcTypes)
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:82
const ElementVolumeVariables & elemVolVars() const
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:116
const ElementFluxVariablesCache & elemFluxVarsCache() const
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:119
const Element & element() const
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:107
NumEqVector lateralDiffusiveMomentumFlux() const
Returns the diffusive momentum flux over the staggered face perpendicular to the scvf where the veloc...
Definition: freeflow/navierstokes/momentum/fluxvariables.hh:236
Helper class for calculating the velocity gradients for the Navier-Stokes model using the staggered g...
Definition: momentum/velocitygradients.hh:58
static auto velocityGradient(const FVElementGeometry &fvGeometry, const typename FVElementGeometry::SubControlVolumeFace &scvf, const ElemVolVars &elemVolVars, bool fullGradient=false)
Definition: momentum/velocitygradients.hh:62
static auto velocityGradII(const FVElementGeometry &fvGeometry, const typename FVElementGeometry::SubControlVolumeFace &scvf, const ElemVolVars &elemVolVars)
Returns the in-axis velocity gradient.
Definition: momentum/velocitygradients.hh:158
Declares all properties used in Dumux.