3.5-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
freeflow/navierstokes/momentum/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_MOMENTUM_FLUXVARIABLES_HH
25#define DUMUX_NAVIERSTOKES_MOMENTUM_FLUXVARIABLES_HH
26
27#include <array>
28
30#include <dumux/common/math.hh>
34
37
38// #include "staggeredupwindhelper.hh"
39#include "velocitygradients.hh"
40
41namespace Dumux {
42
47template<class TypeTag>
49{
51
52 using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
53 using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
54 using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
55
56 using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache;
57 using ElementFluxVariablesCache = typename GridFluxVariablesCache::LocalView;
58 using FluxVariablesCache = typename GridFluxVariablesCache::FluxVariablesCache;
59
60 using GridGeometry = typename GridVariables::GridGeometry;
61 using FVElementGeometry = typename GridGeometry::LocalView;
62 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
63 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
64
66
68 using GridView = typename GridGeometry::GridView;
69 using Element = typename GridView::template Codim<0>::Entity;
72 using Indices = typename ModelTraits::Indices;
74
75 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
76 using Extrusion = Extrusion_t<GridGeometry>;
77
79
80public:
81
83 const Element& element,
84 const FVElementGeometry& fvGeometry,
85 const SubControlVolumeFace& scvFace,
86 const ElementVolumeVariables& elemVolVars,
87 const ElementFluxVariablesCache& elemFluxVarsCache,
88 const ElementBoundaryTypes& elemBcTypes)
89 : problemPtr_(&problem)
90 , elementPtr_(&element)
91 , fvGeometryPtr_(&fvGeometry)
92 , scvFacePtr_(&scvFace)
93 , elemVolVarsPtr_(&elemVolVars)
94 , elemFluxVarsCachePtr_(&elemFluxVarsCache)
95 , elemBcTypesPtr_(&elemBcTypes)
96 {
97 static_assert(
98 std::decay_t<decltype(problem.dirichlet(element, scvFace))>::size()
99 == static_cast<std::size_t>(GridView::dimension),
100 "Expects problem.dirichlet to return an array with as many entries as dimensions."
101 );
102 }
103
104 const Problem& problem() const
105 { return *problemPtr_; }
106
107 const Element& element() const
108 { return *elementPtr_; }
109
110 const SubControlVolumeFace& scvFace() const
111 { return *scvFacePtr_; }
112
113 const FVElementGeometry& fvGeometry() const
114 { return *fvGeometryPtr_; }
115
116 const ElementVolumeVariables& elemVolVars() const
117 { return *elemVolVarsPtr_; }
118
119 const ElementFluxVariablesCache& elemFluxVarsCache() const
120 { return *elemFluxVarsCachePtr_; }
121
122 const ElementBoundaryTypes& elemBcTypes() const
123 { return *elemBcTypesPtr_; }
124
128 NumEqVector advectiveMomentumFlux() const
129 {
130 if (!this->problem().enableInertiaTerms())
131 return NumEqVector(0.0);
132
133 if (this->scvFace().isFrontal())
135 else
137 }
138
142 NumEqVector diffusiveMomentumFlux() const
143 {
144 if (this->scvFace().isFrontal())
146 else
148 }
149
169 {
170 const auto& scvf = this->scvFace();
171 assert(scvf.isFrontal());
172
173 NumEqVector result(0.0);
174 const auto& fvGeometry = this->fvGeometry();
175 const auto& elemVolVars = this->elemVolVars();
176
177 if (scvf.boundary())
178 return result;
179
180 // get the velocity gradient at the normal face's integration point
181 const auto gradV = VelocityGradients::velocityGradient(fvGeometry, scvf, elemVolVars, this->elemBcTypes(), false);
182
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()];
187
188 static const bool enableUnsymmetrizedVelocityGradient
189 = getParamFromGroup<bool>(this->problem().paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradient", false);
190 const Scalar factor = enableUnsymmetrizedVelocityGradient ? 1.0 : 2.0;
191
192 const auto mu = this->problem().effectiveViscosity(this->element(), this->fvGeometry(), this->scvFace());
193 result -= factor * mu * velocityGrad_ii * Extrusion::area(scvf) * elemVolVars[scvf.insideScvIdx()].extrusionFactor();
194
195 static const bool enableDilatationTerm = getParamFromGroup<bool>(this->problem().paramGroup(), "FreeFlow.EnableDilatationTerm", false);
196 if (enableDilatationTerm)
197 {
198 Scalar divergence = 0.0;
199 for (const auto& scv : scvs(fvGeometry))
200 {
201 const auto frontalScvf = *(scvfs(fvGeometry, scv).begin());
202 assert(frontalScvf.isFrontal() && !frontalScvf.boundary());
203 divergence += VelocityGradients::velocityGradII(fvGeometry, frontalScvf, elemVolVars);
204 }
205 // std::cout << "divergence at " << scvf.center() << " is " << divergence << std::endl;
206 // std::cout << std::setprecision(15) << "old term " << factor * mu * velocityGrad_ii * scvf.area() * elemVolVars[scvf.insideScvIdx()].extrusionFactor() << ", div term " << 2.0/3.0 * mu * divergence * scvf.directionSign() * scvf.area() * elemVolVars[scvf.insideScvIdx()].extrusionFactor() << std::endl;
207 result += 2.0/3.0 * mu * divergence * scvf.directionSign() * Extrusion::area(scvf) * elemVolVars[scvf.insideScvIdx()].extrusionFactor();
208 }
209
210
211 return result;
212 }
213
237 {
238 const auto& scvf = this->scvFace();
239 assert(scvf.isLateral());
240
241 NumEqVector result(0.0);
242 const auto& fvGeometry = this->fvGeometry();
243 const auto& elemVolVars = this->elemVolVars();
244 const auto& problem = this->problem();
245 const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
246
247 static const bool enableUnsymmetrizedVelocityGradient
248 = getParamFromGroup<bool>(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradient", false);
249
250 const auto mu = this->problem().effectiveViscosity(this->element(), this->fvGeometry(), this->scvFace());
251
252 // get the velocity gradient at the lateral face's integration point
253 const auto gradV = VelocityGradients::velocityGradient(fvGeometry, scvf, elemVolVars, this->elemBcTypes(), false);
254
255 // Consider the shear stress caused by the gradient of the velocities parallel to our face of interest.
256 GlobalPosition gradVn(0.0);
257 gradV.mv(scvf.unitOuterNormal(), gradVn);
258 const Scalar velocityGrad_ij = gradVn[scv.dofAxis()];
259 result -= mu * velocityGrad_ij;
260
261 // Consider the shear stress caused by the gradient of the velocities normal to our face of interest.
262 if (!enableUnsymmetrizedVelocityGradient)
263 {
264 GlobalPosition gradVTransposedN(0.0);
265 gradV.mtv(scvf.unitOuterNormal(), gradVTransposedN);
266 const Scalar velocityGrad_ji = gradVTransposedN[scv.dofAxis()];
267 result -= mu * velocityGrad_ji;
268 }
269
270 // Account for the area of the staggered lateral face.
271 return result * Extrusion::area(scvf) * elemVolVars[scvf.insideScvIdx()].extrusionFactor();
272 }
273
290 NumEqVector pressureContribution() const
291 {
292 NumEqVector result(0.0);
293 const auto& scvf = this->scvFace();
294 if (scvf.isLateral() || scvf.boundary())
295 return result;
296
297 // The pressure force needs to take the extruded scvf area into account.
298 const auto pressure = this->problem().pressure(this->element(), this->fvGeometry(), scvf);
299 result = pressure*Extrusion::area(scvf)*this->elemVolVars()[scvf.insideScvIdx()].extrusionFactor();
300
301 // The pressure contribution calculated above might have a much larger numerical value compared to the viscous or inertial forces.
302 // This may lead to numerical inaccuracies due to loss of significance (cancellantion) for the final residual value.
303 // In the end, we are only interested in a pressure difference between the two relevant faces so we can
304 // substract a reference value from the actual pressure contribution. Assuming an axisparallel cartesian grid,
305 // scvf.area() will have the same value at both opposing faces such that the reference pressure contribution
306 // cancels out in the final residual which combines the pressure contribution of two adjacent elements
307 // We explicitly do extrude the area here because that might yield different results in both elements.
308 // The multiplication by scvf.area() aims at having a reference value of the same order of magnitude as the actual pressure contribution.
309 const auto referencePressure = this->problem().referencePressure(this->element(), this->fvGeometry(), scvf);
310 result -= referencePressure*scvf.area();
311
312 // Account for the orientation of the face.
313 result *= scvf.directionSign();
314 return result;
315 }
316
336 {
337 NumEqVector flux(0.0);
338 const auto& scvf = this->scvFace();
339 assert(scvf.isFrontal());
340
341 const auto& problem = this->problem();
342 const auto& elemVolVars = this->elemVolVars();
343 const auto velocitySelf = elemVolVars[scvf.insideScvIdx()].velocity();
344 const auto velocityOpposite = elemVolVars[scvf.outsideScvIdx()].velocity();
345
346 // Get the average velocity at the center of the element (i.e. the location of the staggered face).
347 const Scalar transportingVelocity = (velocitySelf + velocityOpposite) * 0.5;
348 const Scalar density = this->problem().density(this->element(), this->fvGeometry(), scvf);
349 const bool selfIsUpstream = scvf.directionSign() == sign(transportingVelocity);
350
351 // TODO use higher order helper
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;
355
356 return transportingVelocity * transportedMomentum * scvf.directionSign() * Extrusion::area(scvf) * extrusionFactor_(elemVolVars, scvf);
357 }
358
382 {
383 NumEqVector flux(0.0);
384 const auto& scvf = this->scvFace();
385 assert(scvf.isLateral());
386
387 const auto& problem = this->problem();
388 const auto& elemVolVars = this->elemVolVars();
389 const auto fvGeometry = this->fvGeometry();
390
391 // get the transporting velocity which is perpendicular to our own (inner) velocity
392 const Scalar transportingVelocity = [&]()
393 {
394 const auto& orthogonalScvf = fvGeometry.lateralOrthogonalScvf(scvf);
395 const Scalar innerTransportingVelocity = elemVolVars[orthogonalScvf.insideScvIdx()].velocity();
396
397 static const bool useOldScheme = getParam<bool>("FreeFlow.UseOldTransportingVelocity", true); // TODO how to deprecate?
398
399 if (useOldScheme)
400 {
401 if (scvf.boundary() && fvGeometry.scv(scvf.insideScvIdx()).boundary())
402 {
403 if (this->elemBcTypes()[scvf.localIndex()].isDirichlet(scvf.normalAxis()))
404 return problem.dirichlet(this->element(), scvf)[scvf.normalAxis()];
405 }
406 else
407 return innerTransportingVelocity;
408 }
409
410 // use the Dirichlet velocity as transporting velocity if the lateral face is on a Dirichlet boundary
411 if (scvf.boundary())
412 {
413 if (this->elemBcTypes()[scvf.localIndex()].isDirichlet(scvf.normalAxis()))
414 return 0.5*(problem.dirichlet(this->element(), scvf)[scvf.normalAxis()] + innerTransportingVelocity);
415 }
416
417 if (orthogonalScvf.boundary())
418 {
419 if (this->elemBcTypes()[orthogonalScvf.localIndex()].isDirichlet(scvf.normalAxis()))
420 return 0.5*(problem.dirichlet(this->element(), scvf)[scvf.normalAxis()] + innerTransportingVelocity);
421 else
422 return innerTransportingVelocity; // fallback value, should actually never be called
423 }
424
425 // average the transporting velocity by weighting with the scv volumes
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);
430 }();
431
432 const Scalar transportedMomentum = [&]()
433 {
434 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
435
436 auto getDirichletMomentumFlux = [&]()
437 {
438 return problem.dirichlet(this->element(), scvf)[insideScv.dofAxis()] * this->problem().density(this->element(), insideScv);
439 };
440
441 // use the Dirichlet velocity as for transported momentum if the lateral face is on a Dirichlet boundary
442 if (scvf.boundary())
443 {
444 if (!this->elemBcTypes()[scvf.localIndex()].isDirichlet(insideScv.dofAxis()))
445 DUNE_THROW(Dune::InvalidStateException, "Neither Dirichlet nor Neumann BC set at " << scvf.ipGlobal());
446
447 return getDirichletMomentumFlux();
448 }
449 else
450 {
451 if (fvGeometry.scvfIntegrationPointInConcaveCorner(scvf))
452 {
453 // TODO: we could put this into an assert, as the construction of outsideScvfWithSameIntegrationPoint is quite expensive
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());
457
458 return getDirichletMomentumFlux();
459 }
460 }
461
462 const bool selfIsUpstream = scvf.directionSign() == sign(transportingVelocity);
463
464 const auto innerVelocity = elemVolVars[scvf.insideScvIdx()].velocity();
465 const auto outerVelocity = elemVolVars[scvf.outsideScvIdx()].velocity();
466
467 const auto rho = this->problem().insideAndOutsideDensity(this->element(), fvGeometry, scvf);
468
469 const auto insideMomentum = innerVelocity * rho.first;
470 const auto outsideMomentum = outerVelocity * rho.second;
471
472 // TODO use higher order helper
473 static const auto upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight");
474
475 return selfIsUpstream ? (upwindWeight * insideMomentum + (1.0 - upwindWeight) * outsideMomentum)
476 : (upwindWeight * outsideMomentum + (1.0 - upwindWeight) * insideMomentum);
477 }();
478
479 return transportingVelocity * transportedMomentum * scvf.directionSign() * Extrusion::area(scvf) * extrusionFactor_(elemVolVars, scvf);
480 }
481
482private:
483
484 template<class ElementVolumeVariables, class SubControlVolumeFace>
485 Scalar extrusionFactor_(const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf) const
486 {
487 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
488 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
489 return harmonicMean(insideVolVars.extrusionFactor(), outsideVolVars.extrusionFactor());
490 }
491
492
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_;
500};
501
502} // end namespace Dumux
503
504#endif // DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH
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
Definition: adapt.hh:29
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.