diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh
index 11c85417f2ce2ba1ca59eda1dca2876831686e65..75572cfa64a75140ffd7d0cfc3e11363305366bd 100644
--- a/dumux/freeflow/navierstokes/problem.hh
+++ b/dumux/freeflow/navierstokes/problem.hh
@@ -220,6 +220,21 @@ public:
DUNE_THROW(Dune::NotImplemented, "When using the Beavers-Joseph-Saffman boundary condition, the alpha value must be returned in the acutal problem");
}
+ //! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph-Saffman condition is used
+ const Scalar bjsVelocity(const SubControlVolumeFace& scvf,
+ const SubControlVolumeFace& normalFace,
+ const Scalar& localSubFaceIdx,
+ 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(normalFace);
+ const Scalar alpha = asImp_().alphaBJ(normalFace);
+ return velocitySelf / (alpha / sqrt(K) * scvf.pairData(localSubFaceIdx).parallelDistance + 1.0);
+ }
+
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 2396290a5c15e2bbb0886ae536d787f149394c8c..b3161ea85b9312b83022f989bd724f133a9c3754 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
@@ -391,7 +391,7 @@ private:
{
const auto ghostFace = makeParallelGhostFace_(scvf, localSubFaceIdx);
if (isBJS)
- return bjsVelocity_(problem, scvf, normalFace, localSubFaceIdx, velocitySelf);
+ return problem.bjsVelocity(scvf, normalFace, localSubFaceIdx, velocitySelf);
return problem.dirichlet(element, ghostFace)[Indices::velocity(scvf.directionIndex())];
};
@@ -507,7 +507,7 @@ private:
{
const auto ghostFace = makeParallelGhostFace_(scvf, localSubFaceIdx);
if (isBJS)
- return bjsVelocity_(problem, scvf, normalFace, localSubFaceIdx, innerParallelVelocity);
+ return problem.bjsVelocity(scvf, normalFace, localSubFaceIdx, innerParallelVelocity);
return problem.dirichlet(element, ghostFace)[Indices::velocity(scvf.directionIndex())];
};
@@ -594,22 +594,6 @@ 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
diff --git a/dumux/freeflow/rans/problem.hh b/dumux/freeflow/rans/problem.hh
index 88289aaa4a556896f5843e6943f8e3c4ed6fb8c6..255a3163ef6169c12c47ae29c4b02240abb4113e 100644
--- a/dumux/freeflow/rans/problem.hh
+++ b/dumux/freeflow/rans/problem.hh
@@ -291,23 +291,67 @@ public:
fvGeometry.bindElement(element);
for (auto&& scvf : scvfs(fvGeometry))
{
- unsigned int normDim = scvf.directionIndex();
- if (scvf.boundary() && asImp_().boundaryTypes(element, scvf).isDirichlet(Indices::velocity(normDim)))
+ // adapt calculations for Dirichlet condition
+ unsigned int scvfNormDim = scvf.directionIndex();
+ if (scvf.boundary())
{
for (unsigned int velIdx = 0; velIdx < dim; ++velIdx)
{
- // face Value
+ if (!asImp_().boundaryTypes(element, scvf).isDirichlet(Indices::velocity(velIdx)))
+ continue;
+
Scalar dirichletVelocity = asImp_().dirichlet(element, scvf)[Indices::velocity(velIdx)];
- unsigned int neighborID = neighborID_[elementID][normDim][0];
- if (scvf.center()[normDim] < cellCenter_[elementID][normDim])
- neighborID = neighborID_[elementID][normDim][1];
+ unsigned int neighborID = neighborID_[elementID][scvfNormDim][0];
+ if (scvf.center()[scvfNormDim] < cellCenter_[elementID][scvfNormDim])
+ neighborID = neighborID_[elementID][scvfNormDim][1];
- velocityGradients_[elementID][velIdx][normDim]
+ velocityGradients_[elementID][velIdx][scvfNormDim]
= (velocity_[neighborID][velIdx] - dirichletVelocity)
- / (cellCenter_[neighborID][normDim] - scvf.center()[normDim]);
+ / (cellCenter_[neighborID][scvfNormDim] - scvf.center()[scvfNormDim]);
+ }
+ }
+
+ // Calculate the BJS-velocity by accounting for all sub faces.
+ std::vector<int> bjsNumFaces(dim, 0);
+ std::vector<unsigned int> bjsNeighbor(dim, 0);
+ DimVector bjsVelocityAverage(0.0);
+ DimVector normalNormCoordinate(0.0);
+ unsigned int velIdx = Indices::velocity(scvfNormDim);
+ const int numSubFaces = scvf.pairData().size();
+ for(int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx)
+ {
+ const auto& normalFace = fvGeometry.scvf(scvf.insideScvIdx(), scvf.pairData()[localSubFaceIdx].localNormalFaceIdx);
+
+ // adapt calculations for Beavers-Joseph-Saffman condition
+ unsigned int normalNormDim = normalFace.directionIndex();
+ if (normalFace.boundary() && (asImp_().boundaryTypes(element, normalFace).isBJS(Indices::velocity(velIdx))))
+ {
+ unsigned int neighborID = neighborID_[elementID][normalNormDim][0];
+ if (normalFace.center()[normalNormDim] < cellCenter_[elementID][normalNormDim])
+ neighborID = neighborID_[elementID][normalNormDim][1];
+
+ bjsVelocityAverage[normalNormDim] += ParentType::bjsVelocity(scvf, normalFace, localSubFaceIdx, velocity_[elementID][velIdx]);
+ if (bjsNumFaces[normalNormDim] > 0 && neighborID != bjsNeighbor[normalNormDim])
+ DUNE_THROW(Dune::InvalidStateException, "Two different neighborID should not occur");
+ bjsNeighbor[normalNormDim] = neighborID;
+ normalNormCoordinate[normalNormDim] = normalFace.center()[normalNormDim];
+ bjsNumFaces[normalNormDim]++;
}
}
+ for (unsigned dirIdx = 0; dirIdx < dim; ++dirIdx)
+ {
+ if (bjsNumFaces[dirIdx] == 0)
+ continue;
+
+ unsigned int neighborID = bjsNeighbor[dirIdx];
+ bjsVelocityAverage[dirIdx] /= bjsNumFaces[dirIdx];
+
+ velocityGradients_[elementID][velIdx][dirIdx]
+ = (velocity_[neighborID][velIdx] - bjsVelocityAverage[dirIdx])
+ / (cellCenter_[neighborID][dirIdx] - normalNormCoordinate[dirIdx]);
+
+ }
}
}