diff --git a/dumux/discretization/staggered/freeflow/boundarytypes.hh b/dumux/discretization/staggered/freeflow/boundarytypes.hh index c06008b763185e152a6a2a3bc8aec1a88ce10d1b..0be8e3490b534f70c61a241f590240b2835f71d5 100644 --- a/dumux/discretization/staggered/freeflow/boundarytypes.hh +++ b/dumux/discretization/staggered/freeflow/boundarytypes.hh @@ -57,6 +57,7 @@ public: boundaryInfo_[eqIdx].visited = false; boundaryInfo_[eqIdx].isDirichletCell = false; boundaryInfo_[eqIdx].isSymmetry = false; + boundaryInfo_[eqIdx].isBJS = false; } @@ -120,11 +121,46 @@ public: static_assert(AlwaysFalse::value, "Setting all boundary types to Neumann not permitted!"); } + /*! + * \brief Set a boundary condition for a single equation to + * Beavers-Joseph-Saffman (special case of Dirichlet b.c.). + */ + void setBJS(int eqIdx) + { + resetEq(eqIdx); + boundaryInfo_[eqIdx].visited = true; + boundaryInfo_[eqIdx].isBJS = true; + } + + /*! + * \brief Returns true if an equation is used to specify a + * Beavers-Joseph-Saffman boundary condition. + * + * \param eqIdx The index of the equation + */ + bool isBJS(unsigned eqIdx) const + { + return boundaryInfo_[eqIdx].isBJS; + } + + /*! + * \brief Returns true if some equation is used to specify a + * Beavers-Joseph-Saffman boundary condition. + */ + bool hasBJS() const + { + for (int i = 0; i < numEq; ++i) + if (boundaryInfo_[i].isBJS) + return true; + return false; + } + protected: struct StaggeredFreeFlowBoundaryInfo { bool visited; bool isDirichletCell; bool isSymmetry; + bool isBJS; }; std::array boundaryInfo_; diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh index 18f14136f434a88cdfd71b4a8aadc6adfe0d213c..11c85417f2ce2ba1ca59eda1dca2876831686e65 100644 --- a/dumux/freeflow/navierstokes/problem.hh +++ b/dumux/freeflow/navierstokes/problem.hh @@ -24,6 +24,7 @@ #ifndef DUMUX_NAVIERSTOKES_PROBLEM_HH #define DUMUX_NAVIERSTOKES_PROBLEM_HH +#include #include #include #include @@ -199,6 +200,26 @@ public: const SubControlVolumeFace& scvf) const { return CellCenterPrimaryVariables(0.0); } + /*! + * \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. + */ + Scalar permeability(const SubControlVolumeFace& scvf) const + { + DUNE_THROW(Dune::NotImplemented, "When using the Beavers-Joseph-Saffman boundary condition, the permeability 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 alphaBJ(const SubControlVolumeFace& scvf) const + { + DUNE_THROW(Dune::NotImplemented, "When using the Beavers-Joseph-Saffman boundary condition, the alpha value must be returned in the acutal problem"); + } + private: //! Returns the implementation of the problem (i.e. static polymorphism) diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh index 381d4f202c2b0d88506d7d8494f667d4355cdea0..0b57648768aa3733dd39f9c974b0c359e12afced 100644 --- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh +++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh @@ -291,6 +291,8 @@ public: // Get the face normal to the face the dof lives on. The staggered sub face conincides with half of this normal face. const auto& normalFace = fvGeometry.scvf(eIdx, scvf.pairData()[localSubFaceIdx].localNormalFaceIdx); + bool isBJS = false; // check for Beavers-Joseph-Saffman boundary condition + // Check if there is face/element parallel to our face of interest where the dof lives on. If there is no parallel neighbor, // we are on a boundary where we have to check for boundary conditions. if(!scvf.hasParallelNeighbor(localSubFaceIdx)) @@ -325,13 +327,18 @@ public: * elemVolVars[normalFace.insideScvIdx()].extrusionFactor() * normalFace.area() * 0.5; continue; } + + // Check if we have a Beavers-Joseph-Saffman condition. If yes, the parallel velocity at the boundary is calculated + // accordingly for the advective part and the diffusive part of the normal momentum flux + if(bcTypes.isBJS(Indices::velocity(scvf.directionIndex()))) + isBJS = true; } // If there is no symmetry or Neumann boundary condition for the given sub face, proceed to calculate the tangential momentum flux. if(enableInertiaTerms) - normalFlux += computeAdvectivePartOfLateralMomentumFlux_(problem, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx); + normalFlux += computeAdvectivePartOfLateralMomentumFlux_(problem, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx, isBJS); - normalFlux += computeDiffusivePartOfLateralMomentumFlux_(problem, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx); + normalFlux += computeDiffusivePartOfLateralMomentumFlux_(problem, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx, isBJS); } return normalFlux; } @@ -365,7 +372,8 @@ private: const SubControlVolumeFace& normalFace, const ElementVolumeVariables& elemVolVars, const FaceVariables& faceVars, - const int localSubFaceIdx) + const int localSubFaceIdx, + const bool isBJS) { // Get the transporting velocity, located at the scvf perpendicular to the current scvf where the dof // of interest is located. @@ -382,6 +390,8 @@ private: auto getParallelVelocityFromBoundary = [&]() { const auto ghostFace = makeParallelGhostFace_(scvf, localSubFaceIdx); + if (isBJS) + return bjsVelocity_(problem, scvf, normalFace, localSubFaceIdx, velocitySelf); return problem.dirichlet(element, ghostFace)[Indices::velocity(scvf.directionIndex())]; }; @@ -443,7 +453,8 @@ private: const SubControlVolumeFace& normalFace, const ElementVolumeVariables& elemVolVars, const FaceVariables& faceVars, - const int localSubFaceIdx) + const int localSubFaceIdx, + const bool isBJS) { FacePrimaryVariables normalDiffusiveFlux(0.0); @@ -495,6 +506,8 @@ private: auto getParallelVelocityFromBoundary = [&]() { const auto ghostFace = makeParallelGhostFace_(scvf, localSubFaceIdx); + if (isBJS) + return bjsVelocity_(problem, scvf, normalFace, localSubFaceIdx, innerParallelVelocity); return problem.dirichlet(element, ghostFace)[Indices::velocity(scvf.directionIndex())]; }; @@ -581,6 +594,22 @@ private: const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; return harmonicMean(insideVolVars.extrusionFactor(), outsideVolVars.extrusionFactor()); } + + //! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph-Saffman condition is used + Scalar bjsVelocity_(const Problem& problem, + const SubControlVolumeFace& scvf, + const SubControlVolumeFace& normalFace, + const Scalar& localSubFaceIdx, + const Scalar& velocitySelf) + { + // 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 = problem.permeability(normalFace); + const Scalar alpha = problem.alphaBJ(normalFace); + return velocitySelf / (alpha / sqrt(K) * scvf.pairData(localSubFaceIdx).parallelDistance + 1.0); + } }; } // end namespace