diff --git a/dumux/freeflow/navierstokes/momentum/fluxhelper.hh b/dumux/freeflow/navierstokes/momentum/fluxhelper.hh
index 651b4827ca32154c74d9885fbd855f891ecddaf9..c1cf53d8df1f6433f27e1a91c36b49a1c5a521c5 100644
--- a/dumux/freeflow/navierstokes/momentum/fluxhelper.hh
+++ b/dumux/freeflow/navierstokes/momentum/fluxhelper.hh
@@ -8,23 +8,45 @@
/*!
* \file
* \ingroup NavierStokesModel
- * \copydoc Dumux::NavierStokesMomentumBoundaryFluxHelper
+ * \copydoc Dumux::NavierStokesMomentumBoundaryFlux
*/
-#ifndef DUMUX_NAVIERSTOKES_MOMENTUM_BOUNDARY_FLUXHELPER_HH
-#define DUMUX_NAVIERSTOKES_MOMENTUM_BOUNDARY_FLUXHELPER_HH
+#ifndef DUMUX_FREEFLOW_NAVIERSTOKES_MOMENTUM_BOUNDARY_FLUXHELPER_HH
+#define DUMUX_FREEFLOW_NAVIERSTOKES_MOMENTUM_BOUNDARY_FLUXHELPER_HH
#include <dumux/common/math.hh>
#include <dumux/common/parameters.hh>
#include <dumux/discretization/method.hh>
#include <dumux/freeflow/navierstokes/momentum/velocitygradients.hh>
+#include <dumux/freeflow/navierstokes/slipcondition.hh>
namespace Dumux {
/*!
* \ingroup NavierStokesModel
- * \brief Struct containing flux helper functions to be used in the momentum problem's Neumann function.
+ * \brief Class to compute the boundary flux for the momentum balance of the Navier-Stokes model
+ * This helper class is typically used in the Neumann function of the momentum problem.
+ *
+ * \tparam DiscretizationMethod The discretization method used for the momentum problem
+ * \tparam SlipVelocityPolicy The policy to compute the slip velocity at the boundary
+ *
+ * The slip velocity policy is a class with a static member function `velocity`
+ * that computes the slip velocity at the boundary.
*/
-struct NavierStokesMomentumBoundaryFluxHelper
+template<class DiscretizationMethod, class SlipVelocityPolicy = void>
+struct NavierStokesMomentumBoundaryFlux;
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief Class to compute the boundary flux for the momentum balance of the Navier-Stokes model
+ * This helper class is typically used in the Neumann function of the momentum problem.
+ *
+ * \tparam SlipVelocityPolicy The policy to compute the slip velocity at the boundary
+ *
+ * The slip velocity policy is a class with a static member function `velocity`
+ * that computes the slip velocity at the boundary.
+ */
+template<class SlipVelocityPolicy>
+struct NavierStokesMomentumBoundaryFlux<DiscretizationMethods::FCStaggered, SlipVelocityPolicy>
{
/*!
* \brief Returns the momentum flux a fixed-pressure boundary.
@@ -50,7 +72,6 @@ struct NavierStokesMomentumBoundaryFluxHelper
const bool zeroNormalVelocityGradient = true)
{
// TODO density upwinding?
- static_assert(FVElementGeometry::GridGeometry::discMethod == DiscretizationMethods::fcstaggered);
using MomentumFlux = typename Problem::MomentumFluxType;
constexpr std::size_t dim = static_cast<std::size_t>(FVElementGeometry::GridGeometry::GridView::dimension);
static_assert(
@@ -203,7 +224,6 @@ struct NavierStokesMomentumBoundaryFluxHelper
const ElementFluxVariablesCache& elemFluxVarsCache,
const TangentialVelocityGradient& tangentialVelocityGradient = TangentialVelocityGradient(0.0))
{
- static_assert(FVElementGeometry::GridGeometry::discMethod == DiscretizationMethods::fcstaggered);
using MomentumFlux = typename Problem::MomentumFluxType;
constexpr std::size_t dim = static_cast<std::size_t>(FVElementGeometry::GridGeometry::GridView::dimension);
static_assert(
@@ -239,61 +259,67 @@ struct NavierStokesMomentumBoundaryFluxHelper
* :: staggered half-control-volume (neighbor element)
*
*/
- const Scalar velI = elemVolVars[scvf.insideScvIdx()].velocity();
- // viscous terms
- const Scalar mu = problem.effectiveViscosity(fvGeometry.element(), fvGeometry, scvf);
- const Scalar distance = (scv.dofPosition()- scvf.ipGlobal()).two_norm();
- const Scalar velGradJI = [&]
+ if constexpr (!std::is_void_v<SlipVelocityPolicy>)
{
- if (elemVolVars.hasVolVars(orthogonalScvf.outsideScvIdx()))
- return StaggeredVelocityGradients::velocityGradJI(fvGeometry, scvf, elemVolVars);
- else
- return tangentialVelocityGradient;
- }();
+ const Scalar velI = elemVolVars[scvf.insideScvIdx()].velocity();
- const Scalar slipVelocity = problem.beaversJosephVelocity(fvGeometry, scvf, elemVolVars, velGradJI)[scv.dofAxis()]; // TODO rename to slipVelocity
- const Scalar velGradIJ = (slipVelocity - velI) / distance * scvf.directionSign();
+ // viscous terms
+ const Scalar mu = problem.effectiveViscosity(fvGeometry.element(), fvGeometry, scvf);
+ const Scalar distance = (scv.dofPosition()- scvf.ipGlobal()).two_norm();
+ const Scalar velGradJI = [&]
+ {
+ if (elemVolVars.hasVolVars(orthogonalScvf.outsideScvIdx()))
+ return StaggeredVelocityGradients::velocityGradJI(fvGeometry, scvf, elemVolVars);
+ else
+ return tangentialVelocityGradient;
+ }();
- flux[scv.dofAxis()] -= (mu * (velGradIJ + velGradJI))*scvf.directionSign();
+ const Scalar slipVelocity = SlipVelocityPolicy::velocity(problem, fvGeometry, scvf, elemVolVars, velGradJI)[scv.dofAxis()];
+ const Scalar velGradIJ = (slipVelocity - velI) / distance * scvf.directionSign();
- // advective terms
- if (problem.enableInertiaTerms())
- {
- // transporting velocity corresponds to velJ
- const auto transportingVelocity = [&]()
+ flux[scv.dofAxis()] -= (mu * (velGradIJ + velGradJI))*scvf.directionSign();
+
+ // advective terms
+ if (problem.enableInertiaTerms())
{
- const auto innerTransportingVelocity = elemVolVars[orthogonalScvf.insideScvIdx()].velocity();
+ // transporting velocity corresponds to velJ
+ const auto transportingVelocity = [&]()
+ {
+ const auto innerTransportingVelocity = elemVolVars[orthogonalScvf.insideScvIdx()].velocity();
- if (!elemVolVars.hasVolVars(orthogonalScvf.outsideScvIdx()))
- return innerTransportingVelocity; // fallback
+ if (!elemVolVars.hasVolVars(orthogonalScvf.outsideScvIdx()))
+ return innerTransportingVelocity; // fallback
- const auto outerTransportingVelocity = elemVolVars[orthogonalScvf.outsideScvIdx()].velocity();
+ const auto outerTransportingVelocity = elemVolVars[orthogonalScvf.outsideScvIdx()].velocity();
- // if the orthogonal scvf lies on a boundary and if there are outside volvars, we assume that these come from a Dirichlet condition
- if (orthogonalScvf.boundary())
- return outerTransportingVelocity;
+ // if the orthogonal scvf lies on a boundary and if there are outside volvars, we assume that these come from a Dirichlet condition
+ if (orthogonalScvf.boundary())
+ return outerTransportingVelocity;
- static const bool useOldScheme = getParam<bool>("FreeFlow.UseOldTransportingVelocity", true); // TODO how to deprecate?
- if (useOldScheme)
- return innerTransportingVelocity;
- else
- {
- // average the transporting velocity by weighting with the scv volumes
- const auto insideVolume = fvGeometry.scv(orthogonalScvf.insideScvIdx()).volume();
- const auto outsideVolume = fvGeometry.scv(orthogonalScvf.outsideScvIdx()).volume();
- const auto outerTransportingVelocity = elemVolVars[orthogonalScvf.outsideScvIdx()].velocity();
- return (insideVolume*innerTransportingVelocity + outsideVolume*outerTransportingVelocity) / (insideVolume + outsideVolume);
- }
- }();
+ static const bool useOldScheme = getParam<bool>("FreeFlow.UseOldTransportingVelocity", true); // TODO how to deprecate?
+ if (useOldScheme)
+ return innerTransportingVelocity;
+ else
+ {
+ // average the transporting velocity by weighting with the scv volumes
+ const auto insideVolume = fvGeometry.scv(orthogonalScvf.insideScvIdx()).volume();
+ const auto outsideVolume = fvGeometry.scv(orthogonalScvf.outsideScvIdx()).volume();
+ const auto outerTransportingVelocity = elemVolVars[orthogonalScvf.outsideScvIdx()].velocity();
+ return (insideVolume*innerTransportingVelocity + outsideVolume*outerTransportingVelocity) / (insideVolume + outsideVolume);
+ }
+ }();
- // Do not use upwinding here but directly take the slip velocity located on the boundary. Upwinding with a weight of 0.5
- // would actually prevent second order grid convergence.
- const auto rho = problem.insideAndOutsideDensity(fvGeometry.element(), fvGeometry, scvf);
- const auto transportedMomentum = slipVelocity * rho.second;
+ // Do not use upwinding here but directly take the slip velocity located on the boundary. Upwinding with a weight of 0.5
+ // would actually prevent second order grid convergence.
+ const auto rho = problem.insideAndOutsideDensity(fvGeometry.element(), fvGeometry, scvf);
+ const auto transportedMomentum = slipVelocity * rho.second;
- flux[scv.dofAxis()] += transportingVelocity * transportedMomentum * scvf.directionSign();
+ flux[scv.dofAxis()] += transportingVelocity * transportedMomentum * scvf.directionSign();
+ }
}
+ else
+ DUNE_THROW(Dune::InvalidStateException, "SlipVelocityPolicy needs to be specified as template argument");
}
else if (scv.boundary() && problem.onSlipBoundary(fvGeometry, fvGeometry.frontalScvfOnBoundary(scv)))
{
@@ -317,67 +343,78 @@ struct NavierStokesMomentumBoundaryFluxHelper
* | | |
*
*/
- const Scalar velJ = elemVolVars[orthogonalScvf.insideScvIdx()].velocity();
-
- // viscous terms
- const Scalar mu = problem.effectiveViscosity(fvGeometry.element(), fvGeometry, scvf);
- const Scalar distance = (fvGeometry.scv(orthogonalScvf.insideScvIdx()).dofPosition()- scvf.ipGlobal()).two_norm();
- const Scalar velGradIJ = [&]
+ if constexpr (!std::is_void_v<SlipVelocityPolicy>)
{
- if (elemVolVars.hasVolVars(scvf.outsideScvIdx()))
- return StaggeredVelocityGradients::velocityGradIJ(fvGeometry, scvf, elemVolVars);
- else
- return tangentialVelocityGradient;
- }();
+ const Scalar velJ = elemVolVars[orthogonalScvf.insideScvIdx()].velocity();
- const Scalar slipVelocity = problem.beaversJosephVelocity(fvGeometry, orthogonalScvf, elemVolVars, velGradIJ)[scvf.normalAxis()]; // TODO rename to slipVelocity
- const Scalar velGradJI = (slipVelocity - velJ) / distance * orthogonalScvf.directionSign();
+ // viscous terms
+ const Scalar mu = problem.effectiveViscosity(fvGeometry.element(), fvGeometry, scvf);
+ const Scalar distance = (fvGeometry.scv(orthogonalScvf.insideScvIdx()).dofPosition()- scvf.ipGlobal()).two_norm();
- flux[scv.dofAxis()] -= (mu * (velGradIJ + velGradJI))*scvf.directionSign();
+ const Scalar velGradIJ = [&]
+ {
+ if (elemVolVars.hasVolVars(scvf.outsideScvIdx()))
+ return StaggeredVelocityGradients::velocityGradIJ(fvGeometry, scvf, elemVolVars);
+ else
+ return tangentialVelocityGradient;
+ }();
- // advective terms
- if (problem.enableInertiaTerms())
- {
- // transporting velocity corresponds to velJ
- const auto transportingVelocity = slipVelocity;
+ const Scalar slipVelocity = SlipVelocityPolicy::velocity(problem, fvGeometry, orthogonalScvf, elemVolVars, velGradIJ)[scvf.normalAxis()];
+ const Scalar velGradJI = (slipVelocity - velJ) / distance * orthogonalScvf.directionSign();
- // if the scvf lies on a boundary and if there are outside volvars, we assume that these come from a Dirichlet condition
- if (scvf.boundary() && elemVolVars.hasVolVars(scvf.outsideScvIdx()))
- {
- flux[scv.dofAxis()] += problem.density(fvGeometry.element(), scv)
- * elemVolVars[scvf.outsideScvIdx()].velocity() * transportingVelocity * scvf.directionSign(); // TODO revise density
- return flux;
- }
+ flux[scv.dofAxis()] -= (mu * (velGradIJ + velGradJI))*scvf.directionSign();
- const auto innerVelocity = elemVolVars[scvf.insideScvIdx()].velocity();
- const auto outerVelocity = [&]
+ // advective terms
+ if (problem.enableInertiaTerms())
{
- if (!elemVolVars.hasVolVars(scvf.outsideScvIdx()))
- return innerVelocity; // fallback
- else
- return elemVolVars[scvf.outsideScvIdx()].velocity();
- }();
+ // transporting velocity corresponds to velJ
+ const auto transportingVelocity = slipVelocity;
+
+ // if the scvf lies on a boundary and if there are outside volvars, we assume that these come from a Dirichlet condition
+ if (scvf.boundary() && elemVolVars.hasVolVars(scvf.outsideScvIdx()))
+ {
+ flux[scv.dofAxis()] += problem.density(fvGeometry.element(), scv)
+ * elemVolVars[scvf.outsideScvIdx()].velocity() * transportingVelocity * scvf.directionSign(); // TODO revise density
+ return flux;
+ }
- const auto rho = problem.insideAndOutsideDensity(fvGeometry.element(), fvGeometry, scvf);
- const bool selfIsUpstream = scvf.directionSign() == sign(transportingVelocity);
+ const auto innerVelocity = elemVolVars[scvf.insideScvIdx()].velocity();
+ const auto outerVelocity = [&]
+ {
+ if (!elemVolVars.hasVolVars(scvf.outsideScvIdx()))
+ return innerVelocity; // fallback
+ else
+ return elemVolVars[scvf.outsideScvIdx()].velocity();
+ }();
+
+ const auto rho = problem.insideAndOutsideDensity(fvGeometry.element(), fvGeometry, scvf);
+ const bool selfIsUpstream = scvf.directionSign() == sign(transportingVelocity);
- const auto insideMomentum = innerVelocity * rho.first;
- const auto outsideMomentum = outerVelocity * rho.second;
+ const auto insideMomentum = innerVelocity * rho.first;
+ const auto outsideMomentum = outerVelocity * rho.second;
- static const auto upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight");
+ static const auto upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight");
- const auto transportedMomentum = selfIsUpstream ? (upwindWeight * insideMomentum + (1.0 - upwindWeight) * outsideMomentum)
- : (upwindWeight * outsideMomentum + (1.0 - upwindWeight) * insideMomentum);
+ const auto transportedMomentum = selfIsUpstream ? (upwindWeight * insideMomentum + (1.0 - upwindWeight) * outsideMomentum)
+ : (upwindWeight * outsideMomentum + (1.0 - upwindWeight) * insideMomentum);
- flux[scv.dofAxis()] += transportingVelocity * transportedMomentum * scvf.directionSign();
+ flux[scv.dofAxis()] += transportingVelocity * transportedMomentum * scvf.directionSign();
+ }
}
+ else
+ DUNE_THROW(Dune::InvalidStateException, "SlipVelocityPolicy needs to be specified as template argument");
}
return flux;
}
};
+using NavierStokesMomentumBoundaryFluxHelper
+ [[deprecated("Replace with implementation class `NavierStokesMomentumBoundaryFlux`with template arguments `DiscretizationMethod` and `SlipVelocityPolicy`. This will be removed after 3.9.")]]
+ = NavierStokesMomentumBoundaryFlux<DiscretizationMethods::FCStaggered,
+ NavierStokesSlipVelocity<DiscretizationMethods::FCStaggered, NavierStokes::SlipConditions::BJ>>;
+
} // end namespace Dumux
#endif
diff --git a/dumux/freeflow/navierstokes/momentum/problem.hh b/dumux/freeflow/navierstokes/momentum/problem.hh
index 6a39312e6f6350afa238984990f2d9d565420698..3df164689440324574c54359478b6974014d0f25 100644
--- a/dumux/freeflow/navierstokes/momentum/problem.hh
+++ b/dumux/freeflow/navierstokes/momentum/problem.hh
@@ -427,6 +427,7 @@ public:
* \brief Returns the intrinsic permeability of 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.
*/
+ [[deprecated("Will be removed after release 3.9.")]]
Scalar permeability(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
{
DUNE_THROW(Dune::NotImplemented,
@@ -439,6 +440,7 @@ public:
* \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.
*/
+ [[deprecated("Will be removed after release 3.9. Implement betaBJ instead. ")]]
Scalar alphaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
{
DUNE_THROW(Dune::NotImplemented,
@@ -450,9 +452,11 @@ public:
/*!
* \brief Returns the beta value which is the alpha value divided by the square root of the (scalar-valued) interface permeability.
*/
+ [[deprecated("Needs to be implemented in test problem. Will be removed after release 3.9.")]]
Scalar betaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, const GlobalPosition& tangentialVector) const
{
- const Scalar interfacePermeability = interfacePermeability_(fvGeometry, scvf, tangentialVector);
+ const auto& K = asImp_().permeability(fvGeometry, scvf);
+ const auto interfacePermeability = vtmv(tangentialVector, K, tangentialVector);
using std::sqrt;
return asImp_().alphaBJ(fvGeometry, scvf) / sqrt(interfacePermeability);
}
@@ -460,6 +464,7 @@ public:
/*!
* \brief Returns the velocity in the porous medium (which is 0 by default according to Saffmann).
*/
+ [[deprecated("Needs to be implemented in test problem. Will be removed after release 3.9.")]]
VelocityVector porousMediumVelocity(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
{
return VelocityVector(0.0);
@@ -469,6 +474,7 @@ public:
* \brief Returns the slip velocity at a porous boundary based on the Beavers-Joseph(-Saffman) condition.
* \note This only returns a vector filled with one component of the slip velocity (corresponding to the dof axis of the scv the svf belongs to)
*/
+ [[deprecated("Will be removed after release 3.9.")]]
const VelocityVector beaversJosephVelocity(const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
@@ -511,19 +517,6 @@ public:
}
private:
- //! Returns a scalar permeability value at the coupling interface
- template<class Scvf>
- Scalar interfacePermeability_(const FVElementGeometry& fvGeometry, const Scvf& scvf, const GlobalPosition& tangentialVector) const
- {
- const auto& K = asImp_().permeability(fvGeometry, scvf);
-
- // use t*K*t for permeability tensors
- if constexpr (Dune::IsNumber<std::decay_t<decltype(K)>>::value)
- return K;
- else
- return vtmv(tangentialVector, K, tangentialVector);
- }
-
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation& asImp_()
{ return *static_cast<Implementation *>(this); }
diff --git a/dumux/freeflow/navierstokes/slipcondition.hh b/dumux/freeflow/navierstokes/slipcondition.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5888d263c84cd4ccc8d8d6bd969b1cd36cdc287e
--- /dev/null
+++ b/dumux/freeflow/navierstokes/slipcondition.hh
@@ -0,0 +1,126 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup NavierStokesModel
+ * \brief Navier Stokes slip condition
+ */
+#ifndef DUMUX_FREEFLOW_NAVIERSTOKES_SLIPCONDITION_HH
+#define DUMUX_FREEFLOW_NAVIERSTOKES_SLIPCONDITION_HH
+
+#include <dumux/common/tag.hh>
+#include <dumux/discretization/method.hh>
+
+namespace Dumux::NavierStokes::SlipConditions {
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief Tag for the Beavers-Joseph slip condition
+ */
+struct BJ : public Utility::Tag<BJ> {
+ static std::string name() { return "Beavers-Joseph"; }
+};
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief Tag for the Beavers-Joseph-Saffman slip condition
+ */
+struct BJS : public Utility::Tag<BJS> {
+ static std::string name() { return "Beavers-Joseph-Saffman"; }
+};
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief Tag for the Beavers-Joseph slip condition
+ */
+inline constexpr BJ bj{};
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief Tag for the Beavers-Joseph-Saffman slip condition
+ */
+inline constexpr BJS bjs{};
+
+} // end namespace Dumux::NavierStokes::SlipConditions
+
+namespace Dumux {
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief Navier Stokes slip velocity policy
+ */
+template<class DiscretizationMethod, class SlipCondition>
+class NavierStokesSlipVelocity;
+
+/*!
+ * \ingroup NavierStokesModel
+ * \brief Navier Stokes slip velocity helper for fcstaggered discretization
+ *
+ * For now, this class implements the Beavers-Joseph or the Beavers-Joseph-Saffman condition
+ * which models the slip velocity at a porous boundary. The condition is chosen by passing
+ * the corresponding SlipCondition tag to the class.
+ */
+template<class SlipCondition>
+struct NavierStokesSlipVelocity<DiscretizationMethods::FCStaggered, SlipCondition>
+{
+ static constexpr SlipCondition slipCondition{};
+
+ /*!
+ * \brief Returns the slip velocity at a porous boundary based on the Beavers-Joseph(-Saffman) condition.
+ * \note This only returns a vector filled with one component of the slip velocity (corresponding to the dof axis of the scv the svf belongs to)
+ * \param problem The problem
+ * \param fvGeometry The finite-volume geometry
+ * \param scvf The sub control volume face
+ * \param elemVolVars The volume variables for the element
+ * \param tangentialVelocityDeriv Pre-calculated velocity derivative
+ */
+ template<class Problem, class FVElementGeometry, class ElementVolumeVariables, class Scalar>
+ static auto velocity(const Problem& problem,
+ const FVElementGeometry& fvGeometry,
+ const typename FVElementGeometry::SubControlVolumeFace& scvf,
+ const ElementVolumeVariables& elemVolVars,
+ Scalar tangentialVelocityDeriv)
+ {
+ assert(scvf.isLateral());
+ assert(scvf.boundary());
+
+ static constexpr int dimWorld = FVElementGeometry::GridGeometry::GridView::dimensionworld;
+ using Vector = Dune::FieldVector<Scalar, dimWorld>;
+
+ Vector porousMediumVelocity(0.0);
+
+ if constexpr (slipCondition == NavierStokes::SlipConditions::bj)
+ porousMediumVelocity = problem.porousMediumVelocity(fvGeometry, scvf);
+ else if (!(slipCondition == NavierStokes::SlipConditions::bjs))
+ DUNE_THROW(Dune::NotImplemented, "Fcstaggered currently only implements BJ or BJS slip conditions");
+
+ const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
+
+ // create a unit normal vector oriented in positive coordinate direction
+ Vector tangent(0.0);
+ tangent[scv.dofAxis()] = 1.0;
+
+ // du/dy + dv/dx = beta * (u_boundary-uPM)
+ // beta = alpha/sqrt(K)
+ const Scalar betaBJ = problem.betaBJ(fvGeometry, scvf, tangent);
+ const Scalar distanceNormalToBoundary = (scvf.ipGlobal() - scv.dofPosition()).two_norm();
+
+ static const bool onlyNormalGradient = getParamFromGroup<bool>(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false);
+ if (onlyNormalGradient)
+ tangentialVelocityDeriv = 0.0;
+
+ const Scalar scalarSlipVelocity = (tangentialVelocityDeriv*distanceNormalToBoundary
+ + porousMediumVelocity * tangent * betaBJ * distanceNormalToBoundary
+ + elemVolVars[scv].velocity()) / (betaBJ*distanceNormalToBoundary + 1.0);
+
+ return scalarSlipVelocity*tangent;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/boundary/freeflowporenetwork/1p_1p/problem_freeflow.hh b/test/multidomain/boundary/freeflowporenetwork/1p_1p/problem_freeflow.hh
index 14005f17791440f4ad9e12c1e7112631a6d95f75..0c2cde830521373443d436259c25a29442c57ccc 100644
--- a/test/multidomain/boundary/freeflowporenetwork/1p_1p/problem_freeflow.hh
+++ b/test/multidomain/boundary/freeflowporenetwork/1p_1p/problem_freeflow.hh
@@ -166,7 +166,8 @@ public:
{
BoundaryFluxes values(0.0);
const auto& globalPos = scvf.ipGlobal();
- using FluxHelper = NavierStokesMomentumBoundaryFluxHelper;
+ using SlipVelocityPolicy = NavierStokesSlipVelocity<typename GridGeometry::DiscretizationMethod, NavierStokes::SlipConditions::BJ>;
+ using FluxHelper = NavierStokesMomentumBoundaryFlux<typename GridGeometry::DiscretizationMethod, SlipVelocityPolicy>;
if constexpr (ParentType::isMomentumProblem())
{
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_freeflow.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_freeflow.hh
index be65c269cf22a0d86921b831fa1fe32a4f87a0e1..3b3912ef1154cec99c744bc6374495e670d755af 100644
--- a/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_freeflow.hh
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/convergence/problem_freeflow.hh
@@ -164,7 +164,8 @@ public:
{
BoundaryFluxes values(0.0);
const auto& globalPos = scvf.ipGlobal();
- using FluxHelper = NavierStokesMomentumBoundaryFluxHelper;
+ using SlipVelocityPolicy = NavierStokesSlipVelocity<typename GridGeometry::DiscretizationMethod, NavierStokes::SlipConditions::BJS>;
+ using FluxHelper = NavierStokesMomentumBoundaryFlux<typename GridGeometry::DiscretizationMethod, SlipVelocityPolicy>;
// momentum boundary conditions
if constexpr (ParentType::isMomentumProblem())
@@ -283,18 +284,16 @@ public:
{ return InitialValues(0.0); }
/*!
- * \brief Returns the intrinsic permeability of required as input parameter
- for the Beavers-Joseph-Saffman boundary condition
+ * \brief Returns the beta value which is the alpha value divided by the square root of the (scalar-valued) interface permeability.
*/
- auto permeability(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
- { return couplingManager().darcyPermeability(fvGeometry, scvf); }
-
- /*!
- * \brief Returns the alpha value required as input parameter for the
- Beavers-Joseph-Saffman boundary condition.
- */
- Scalar alphaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
- { return couplingManager().problem(CouplingManager::porousMediumIndex).spatialParams().beaversJosephCoeffAtPos(scvf.ipGlobal()); }
+ Scalar betaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, const GlobalPosition& tangentialVector) const
+ {
+ const auto& K = couplingManager().darcyPermeability(fvGeometry, scvf);
+ const auto alphaBJ = couplingManager().problem(CouplingManager::porousMediumIndex).spatialParams().beaversJosephCoeffAtPos(scvf.ipGlobal());
+ const auto interfacePermeability = vtmv(tangentialVector, K, tangentialVector);
+ using std::sqrt;
+ return alphaBJ / sqrt(interfacePermeability);
+ }
/*!
* \brief Returns the analytical solution of the problem at a given position.
diff --git a/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_freeflow.hh b/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_freeflow.hh
index f1a8da197408ab30c22f2bd9c3f039f7d5d4368b..a11c9c939711d4344f928f861ef287d24f5a7463 100644
--- a/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_freeflow.hh
+++ b/test/multidomain/boundary/freeflowporousmedium/1p_1p/problem_freeflow.hh
@@ -180,7 +180,8 @@ public:
{
BoundaryFluxes values(0.0);
const auto& globalPos = scvf.ipGlobal();
- using FluxHelper = NavierStokesMomentumBoundaryFluxHelper;
+ using SlipVelocityPolicy = NavierStokesSlipVelocity<typename GridGeometry::DiscretizationMethod, NavierStokes::SlipConditions::BJS>;
+ using FluxHelper = NavierStokesMomentumBoundaryFlux<typename GridGeometry::DiscretizationMethod, SlipVelocityPolicy>;
if constexpr (ParentType::isMomentumProblem())
{
@@ -236,16 +237,16 @@ public:
}
/*!
- * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ * \brief Returns the beta value which is the alpha value divided by the square root of the (scalar-valued) interface permeability.
*/
- Scalar permeability(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
- { return couplingManager_->darcyPermeability(fvGeometry, scvf); }
-
- /*!
- * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
- */
- Scalar alphaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
- { return couplingManager_->problem(CouplingManager::porousMediumIndex).spatialParams().beaversJosephCoeffAtPos(scvf.ipGlobal()); }
+ Scalar betaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, const GlobalPosition& tangentialVector) const
+ {
+ const auto& K = couplingManager_->darcyPermeability(fvGeometry, scvf);
+ const auto alphaBJ = couplingManager_->problem(CouplingManager::porousMediumIndex).spatialParams().beaversJosephCoeffAtPos(scvf.ipGlobal());
+ const auto interfacePermeability = vtmv(tangentialVector, K, tangentialVector);
+ using std::sqrt;
+ return alphaBJ / sqrt(interfacePermeability);
+ }
/*!
* \brief Returns the pressure difference between inlet and outlet for the horizontal flow scenario