Commit 1dda7c37 authored by Ned Coltman's avatar Ned Coltman
Browse files

Merge branch 'feature/improve-beavers-joseph' into 'master'

[staggered] Improve Beavers-Joseph slip condition

Closes #713

See merge request !1639
parents a82f5aff 98c6e38f
......@@ -368,7 +368,6 @@ private:
std::vector<GridIndexType> scvIndices_;
bool boundary_;
int dofIdx_;
Scalar selfToOppositeDistance_;
AxisData axisData_;
std::array<PairData, numPairs> pairData_;
......
......@@ -24,6 +24,8 @@
#ifndef DUMUX_NAVIERSTOKES_PROBLEM_HH
#define DUMUX_NAVIERSTOKES_PROBLEM_HH
#include <dune/common/deprecated.hh>
#include <dune/common/exceptions.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/staggeredfvproblem.hh>
......@@ -69,6 +71,7 @@ class NavierStokesProblem : public NavierStokesParentProblem<TypeTag>
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using GridFaceVariables = typename GridVariables::GridFaceVariables;
using ElementFaceVariables = typename GridFaceVariables::LocalView;
using FaceVariables = typename GridFaceVariables::FaceVariables;
using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
......@@ -204,20 +207,46 @@ public:
DUNE_THROW(Dune::NotImplemented, "When using the Beavers-Joseph-Saffman boundary condition, the alpha value must be returned in the acutal problem");
}
/*!
* \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
*
* This member function must be overloaded in the problem implementation, if the BJS boundary condition is used.
*/
Scalar betaBJ(const Element& element, const SubControlVolumeFace& scvf) const
{
using std::sqrt;
return asImp_().alphaBJ(scvf) / sqrt(asImp_().permeability(element, scvf));
}
/*!
* \brief Returns the velocity in the porous medium (which is 0 by default according to Saffmann).
*/
Scalar velocityPorousMedium(const Element& element, const SubControlVolumeFace& scvf) const
{ return 0.0; }
//! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph-Saffman condition is used
DUNE_DEPRECATED_MSG("Use beaversJosephVelocity(element, scv, faceOnPorousBoundary, velocitySelf, tangentialVelocityGradient) instead")
const Scalar bjsVelocity(const Element& element,
const SubControlVolume& scv,
const SubControlVolumeFace& faceOnPorousBoundary,
const Scalar velocitySelf) const
{
// du/dy = alpha/sqrt(K) * u_boundary
// du/dy = (u_center - u_boundary) / deltaY
// u_boundary = u_center / (alpha/sqrt(K)*deltaY + 1)
using std::sqrt;
const Scalar K = asImp_().permeability(element, faceOnPorousBoundary);
const Scalar alpha = asImp_().alphaBJ(faceOnPorousBoundary);
// assume tangential velocity gradient of zero and
return beaversJosephVelocity(element, scv, faceOnPorousBoundary, velocitySelf, 0.0);
}
//! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph condition is used
const Scalar beaversJosephVelocity(const Element& element,
const SubControlVolume& scv,
const SubControlVolumeFace& faceOnPorousBoundary,
const Scalar velocitySelf,
const Scalar tangentialVelocityGradient) const
{
// du/dy + dv/dx = alpha/sqrt(K) * (u_boundary-uPM)
// beta = alpha/sqrt(K)
const Scalar betaBJ = asImp_().betaBJ(element, faceOnPorousBoundary);
const Scalar distance = (faceOnPorousBoundary.center() - scv.center()).two_norm();
return velocitySelf / (alpha / sqrt(K) * distance + 1.0);
return (tangentialVelocityGradient*distance + asImp_().velocityPorousMedium(element,faceOnPorousBoundary)*betaBJ*distance + velocitySelf) / (betaBJ*distance + 1.0);
}
private:
......
......@@ -33,6 +33,7 @@
#include <dumux/discretization/method.hh>
#include <dumux/freeflow/staggeredupwindmethods.hh>
#include "velocitygradients.hh"
namespace Dumux {
......@@ -69,6 +70,7 @@ class StaggeredUpwindFluxVariables
using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using VelocityGradients = StaggeredVelocityGradients<Scalar, FVGridGeometry, BoundaryTypes, Indices>;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
......@@ -121,6 +123,7 @@ public:
const FaceVariables& faceVars,
const GridFluxVariablesCache& gridFluxVarsCache,
const int localSubFaceIdx,
const Dune::Std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes)
{
const auto eIdx = scvf.insideScvIdx();
......@@ -138,15 +141,15 @@ public:
if (canHigherOrder)
{
const auto parallelUpwindingMomenta = getLateralUpwindingMomenta_(problem, fvGeometry, element, scvf, elemVolVars, faceVars,
transportingVelocity, localSubFaceIdx, lateralFaceBoundaryTypes,
std::integral_constant<bool, useHigherOrder>{});
transportingVelocity, localSubFaceIdx, currentScvfBoundaryTypes,
lateralFaceBoundaryTypes, std::integral_constant<bool, useHigherOrder>{});
return doLateralMomentumUpwinding_(fvGeometry, scvf, parallelUpwindingMomenta, transportingVelocity, localSubFaceIdx, gridFluxVarsCache);
}
else
{
const auto parallelUpwindingMomenta = getLateralUpwindingMomenta_(problem, fvGeometry, element, scvf, elemVolVars, faceVars,
transportingVelocity, localSubFaceIdx, lateralFaceBoundaryTypes,
std::integral_constant<bool, false>{});
transportingVelocity, localSubFaceIdx, currentScvfBoundaryTypes,
lateralFaceBoundaryTypes, std::integral_constant<bool, false>{});
return doLateralMomentumUpwinding_(fvGeometry, scvf, parallelUpwindingMomenta, transportingVelocity, localSubFaceIdx, gridFluxVarsCache);
}
}
......@@ -186,11 +189,7 @@ private:
{
// Depending on selfIsUpstream I have to check if I have a forward or a backward neighbor to retrieve
const bool selfIsUpstream = ownScvf.directionSign() != sign(transportingVelocity);
if (selfIsUpstream)
return ownScvf.hasForwardNeighbor(0);
else
return ownScvf.hasBackwardNeighbor(0);
return selfIsUpstream ? ownScvf.hasForwardNeighbor(0) : ownScvf.hasBackwardNeighbor(0);
}
/*!
......@@ -369,6 +368,7 @@ private:
const FaceVariables& faceVars,
const Scalar transportingVelocity,
const int localSubFaceIdx,
const Dune::Std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
std::true_type)
{
......@@ -383,8 +383,8 @@ private:
if (!ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
{
const Scalar boundaryMomentum = getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf,
faceVars.velocitySelf(), localSubFaceIdx,
lateralFaceBoundaryTypes) * insideVolVars.density();
faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes,
localSubFaceIdx) * insideVolVars.density();
return std::array<Scalar, 3>{boundaryMomentum, boundaryMomentum, boundaryMomentum};
}
......@@ -402,8 +402,8 @@ private:
momenta[0] = faceVars.velocityParallel(localSubFaceIdx, 0) * insideVolVars.density();
else
momenta[0] = getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf,
faceVars.velocitySelf(), localSubFaceIdx,
lateralFaceBoundaryTypes) * insideVolVars.density();
faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes,
localSubFaceIdx) * insideVolVars.density();
// The local index of the faces that is opposite to localSubFaceIdx
const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1;
......@@ -450,6 +450,7 @@ private:
const FaceVariables& faceVars,
const Scalar transportingVelocity,
const int localSubFaceIdx,
const Dune::Std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
std::false_type)
{
......@@ -463,8 +464,8 @@ private:
const Scalar momentumParallel = ownScvf.hasParallelNeighbor(localSubFaceIdx, 0)
? faceVars.velocityParallel(localSubFaceIdx, 0) * outsideVolVars.density()
: (getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf,
faceVars.velocitySelf(), localSubFaceIdx,
lateralFaceBoundaryTypes)
faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes,
localSubFaceIdx)
* insideVolVars.density());
// If the lateral face lies on a boundary, we assume that the parallel velocity on the boundary is actually known,
......@@ -572,13 +573,15 @@ private:
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const Scalar velocitySelf,
const int localSubFaceIdx,
const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes)
const FaceVariables& faceVars,
const Dune::Std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
const int localSubFaceIdx)
{
// Find out what boundary type is set on the lateral face
const bool lateralFaceHasDirichletPressure = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx);
const bool lateralFaceHasBJS = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isBJS(Indices::velocity(scvf.directionIndex()));
const Scalar velocitySelf = faceVars.velocitySelf();
// If there is a Dirichlet condition for the pressure we assume zero gradient for the velocity,
// so the velocity at the boundary equal to that on the scvf.
......@@ -589,9 +592,11 @@ private:
{
const auto eIdx = scvf.insideScvIdx();
const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
const Scalar velocityGrad_ji = VelocityGradients::velocityGradJI(problem, element, fvGeometry, scvf, faceVars,
currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
const SubControlVolume& scv = fvGeometry.scv(scvf.insideScvIdx());
return problem.bjsVelocity(element, scv, lateralFace, velocitySelf);
return problem.beaversJosephVelocity(element, scv, lateralFace, velocitySelf, velocityGrad_ji);
}
else
{
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup NavierStokesModel
* \copydoc Dumux::StaggeredVelocityGradients
*/
#ifndef DUMUX_NAVIERSTOKES_STAGGERED_VELOCITYGRADIENTS_HH
#define DUMUX_NAVIERSTOKES_STAGGERED_VELOCITYGRADIENTS_HH
#include <dune/common/std/optional.hh>
#include <dumux/common/exceptions.hh>
#include <dumux/common/parameters.hh>
namespace Dumux {
/*!
* \ingroup NavierStokesModel
* \brief Helper class for calculating the velocity gradients for the Navier-Stokes model using the staggered grid discretization.
*/
template<class Scalar, class FVGridGeometry, class BoundaryTypes, class Indices>
class StaggeredVelocityGradients
{
using FVElementGeometry = typename FVGridGeometry::LocalView;
using GridView = typename FVGridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
/*!
* \brief Returns the in-axis velocity gradient.
*
* \verbatim
* ---------======= == and # staggered half-control-volume
* | # | current scvf
* | # | # staggered face over which fluxes are calculated
* vel.Opp <~~| O~~> x~~~~> vel.Self
* | # | x dof position
* scvf | # |
* --------======== -- element
*
* O position at which gradient is evaluated
* \endverbatim
*/
template<class FaceVariables>
static Scalar velocityGradII(const SubControlVolumeFace& scvf,
const FaceVariables& faceVars)
{
// The velocities of the dof at interest and the one of the opposite scvf.
const Scalar velocitySelf = faceVars.velocitySelf();
const Scalar velocityOpposite = faceVars.velocityOpposite();
return ((velocityOpposite - velocitySelf) / scvf.selfToOppositeDistance()) * scvf.directionSign();
}
/*!
* \brief Returns the velocity gradient perpendicular to the orientation of our current scvf.
*
* \verbatim
* ----------------
* | |vel.
* | |Parallel
* | |~~~~> ------->
* | | ------> * gradient
* | | ----->
* scvf ---------######O::::::::: ----> || and # staggered half-control-volume (own element)
* | || | curr. :: --->
* | || | scvf :: --> :: staggered half-control-volume (neighbor element)
* | || x~~~~> :: ->
* | || | vel. :: # lateral staggered faces over which fluxes are calculated
* scvf | || | Self ::
* ---------#######::::::::: x dof position
* scvf
* -- elements
*
* O position at which gradient is evaluated
* \endverbatim
*/
template<class Problem, class FaceVariables>
static Scalar velocityGradIJ(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const FaceVariables& faceVars,
const Dune::Std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
const std::size_t localSubFaceIdx)
{
const auto eIdx = scvf.insideScvIdx();
const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
// For the velocityGrad_ij derivative, get the velocities at the current (own) scvf
// and at the parallel one at the neighboring scvf.
const Scalar innerParallelVelocity = faceVars.velocitySelf();
const auto outerParallelVelocity = [&]()
{
if (!lateralScvf.boundary())
return faceVars.velocityParallel(localSubFaceIdx, 0);
else if (lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex())))
{
// Sample the value of the Dirichlet BC at the center of the staggered lateral face.
const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
return problem.dirichlet(element, lateralScvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(scvf.directionIndex())];
}
else if (lateralFaceBoundaryTypes->isBJS(Indices::velocity(scvf.directionIndex())))
{
const auto tangentialVelocityGradient = [&]()
{
// If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a BJ condition for
// the slip velocity is set there, assume a tangential velocity gradient of zero along the lateral face
// (towards the current scvf).
static const bool unsymmetrizedGradientForBJ = getParamFromGroup<bool>(problem.paramGroup(),
"FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false);
if (unsymmetrizedGradientForBJ)
return 0.0;
if (scvf.boundary())
{
if (currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
currentScvfBoundaryTypes->isBJS(Indices::velocity(lateralScvf.directionIndex())))
return 0.0;
}
return velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
}();
return problem.beaversJosephVelocity(element, fvGeometry.scv(scvf.insideScvIdx()),
lateralScvf, innerParallelVelocity,
tangentialVelocityGradient);
}
else
DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary type at " << lateralScvf.center());
}();
// The velocity gradient already accounts for the orientation
// of the staggered face's outer normal vector. This also correctly accounts for the reduced
// distance used in the gradient if the lateral scvf lies on a boundary.
return (outerParallelVelocity - innerParallelVelocity)
/ scvf.parallelDofsDistance(localSubFaceIdx, 0) * lateralScvf.directionSign();
}
/*!
* \brief Returns the velocity gradient in line with our current scvf.
*
* \verbatim
* ^ gradient
* | ^
* | | ^
* | | | ^
* | | | | ^
* | | | | | ^
* | | | | | |
*
* ----------------
* | |
* | in.norm. |
* | vel. |
* | ^ | ^ out.norm.vel.
* | | | |
* scvf ---------######O::::::::: || and # staggered half-control-volume (own element)
* | || | curr. ::
* | || | scvf :: :: staggered half-control-volume (neighbor element)
* | || x~~~~> ::
* | || | vel. :: # lateral staggered faces over which fluxes are calculated
* scvf | || | Self ::
* ---------#######::::::::: x dof position
* scvf
* -- elements
*
* O position at which gradient is evaluated
* \endverbatim
*/
template<class Problem, class FaceVariables>
static Scalar velocityGradJI(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const FaceVariables& faceVars,
const Dune::Std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
const Dune::Std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
const std::size_t localSubFaceIdx)
{
const auto eIdx = scvf.insideScvIdx();
const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
// Assume a zero velocity gradient for pressure boundary conditions.
if (currentScvfBoundaryTypes && currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx))
return 0.0;
// For the velocityGrad_ji gradient, get the velocities perpendicular to the velocity at the current scvf.
// The inner one is located at staggered face within the own element,
// the outer one at the respective staggered face of the element on the other side of the
// current scvf.
const Scalar innerLateralVelocity = faceVars.velocityLateralInside(localSubFaceIdx);
const Scalar outerLateralVelocity = [&]()
{
if (!scvf.boundary())
return faceVars.velocityLateralOutside(localSubFaceIdx);
else if (currentScvfBoundaryTypes->isDirichlet(Indices::velocity(lateralScvf.directionIndex())))
{
// Sample the value of the Dirichlet BC at the center of the lateral face intersecting with the boundary.
const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx);
return problem.dirichlet(element, scvf.makeBoundaryFace(lateralBoundaryFacePos))[Indices::velocity(lateralScvf.directionIndex())];
}
else if (currentScvfBoundaryTypes->isBJS(Indices::velocity(lateralScvf.directionIndex())))
{
const auto tangentialVelocityGradient = [&]()
{
// If the current scvf is on a boundary and if a Dirichlet BC for the pressure or a BJ condition for
// the slip velocity is set there, assume a tangential velocity gradient of zero along the lateral face
// (towards the current scvf).
static const bool unsymmetrizedGradientForBJ = getParamFromGroup<bool>(problem.paramGroup(),
"FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false);
if (unsymmetrizedGradientForBJ)
return 0.0;
if (lateralScvf.boundary())
{
if (lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
lateralFaceBoundaryTypes->isBJS(Indices::velocity(scvf.directionIndex())))
return 0.0;
}
return velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
}();
return problem.beaversJosephVelocity(element, fvGeometry.scv(scvf.insideScvIdx()),
scvf, innerLateralVelocity,
tangentialVelocityGradient);
}
else
DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary types at " << lateralScvf.center());
}();
// Calculate the velocity gradient in positive coordinate direction.
const Scalar lateralDeltaV = scvf.normalInPosCoordDir()
? (outerLateralVelocity - innerLateralVelocity)
: (innerLateralVelocity - outerLateralVelocity);
return lateralDeltaV / scvf.pairData(localSubFaceIdx).lateralDistance;
}
private:
/*!
* \brief Get the location of the lateral staggered face's center.
* Only needed for boundary conditions if the current scvf or the lateral one is on a bounary.
*
* \verbatim
* --------#######o || frontal face of staggered half-control-volume
* | || | current scvf # lateral staggered face of interest (may lie on a boundary)
* | || | x dof position
* | || x~~~~> vel.Self -- element boundaries, current scvf may lie on a boundary
* | || | o position at which the boundary conditions will be evaluated
* | || | (lateralStaggeredFaceCenter)
* ----------------
* \endverbatim
*/
static const GlobalPosition& lateralStaggeredFaceCenter_(const SubControlVolumeFace& scvf, const int localSubFaceIdx)
{
return scvf.pairData(localSubFaceIdx).lateralStaggeredFaceCenter;
};
};
} // end namespace Dumux
#endif // DUMUX_NAVIERSTOKES_STAGGERED_VELOCITYGRADIENTS_HH
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment