Commit 57427c3d authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Ned Coltman
Browse files

[freeflow][problem] Adapt signature of bjsVelocity

* pass only the face normal to our scvf of interest
* calculate distance within function
parent 9869aa19
......@@ -74,6 +74,7 @@ class NavierStokesProblem : public NavierStokesParentProblem<TypeTag>
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
......@@ -205,18 +206,18 @@ public:
//! helper function to evaluate the slip velocity on the boundary when the Beavers-Joseph-Saffman condition is used
const Scalar bjsVelocity(const Element& element,
const SubControlVolumeFace& scvf,
const SubControlVolumeFace& normalFace,
const Scalar& localSubFaceIdx,
const Scalar& velocitySelf) const
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, normalFace);
const Scalar alpha = asImp_().alphaBJ(normalFace);
return velocitySelf / (alpha / sqrt(K) * scvf.cellCenteredParallelDistance(localSubFaceIdx,0) + 1.0);
const Scalar K = asImp_().permeability(element, faceOnPorousBoundary);
const Scalar alpha = asImp_().alphaBJ(faceOnPorousBoundary);
const Scalar distance = (faceOnPorousBoundary.center() - scv.center()).two_norm();
return velocitySelf / (alpha / sqrt(K) * distance + 1.0);
}
private:
......
......@@ -493,7 +493,10 @@ private:
else if (lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex())))
return problem.dirichlet(element, makeParallelGhostFace_(scvf, localSubFaceIdx))[Indices::velocity(scvf.directionIndex())];
else
return problem.bjsVelocity(element, scvf, normalFace, localSubFaceIdx, innerParallelVelocity);
{
const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
return problem.bjsVelocity(element, scv, normalFace, innerParallelVelocity);
}
};
const Scalar outerParallelVelocity = getParallelVelocity();
......
......@@ -64,6 +64,7 @@ class StaggeredUpwindFluxVariables
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using Indices = typename ModelTraits::Indices;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>;
......@@ -383,9 +384,9 @@ private:
// thus we always use this value for the computation of the transported momentum.
if (!ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
{
const Scalar boundaryMomentum = getParallelVelocityFromBoundary_(problem, ownScvf, lateralFace,
const Scalar boundaryMomentum = getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf, lateralFace,
faceVars.velocitySelf(), localSubFaceIdx,
element, lateralFaceHasDirichletPressure,
lateralFaceHasDirichletPressure,
lateralFaceHasBJS) * insideVolVars.density();
return std::array<Scalar, 3>{boundaryMomentum, boundaryMomentum, boundaryMomentum};
......@@ -403,9 +404,9 @@ private:
if (ownScvf.hasParallelNeighbor(localSubFaceIdx, 0))
momenta[0] = faceVars.velocityParallel(localSubFaceIdx, 0) * insideVolVars.density();
else
momenta[0] = getParallelVelocityFromBoundary_(problem, ownScvf, lateralFace,
momenta[0] = getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf, lateralFace,
faceVars.velocitySelf(), localSubFaceIdx,
element, lateralFaceHasDirichletPressure,
lateralFaceHasDirichletPressure,
lateralFaceHasBJS) * insideVolVars.density();
// The local index of the faces that is opposite to localSubFaceIdx
......@@ -466,8 +467,8 @@ private:
const Scalar momentumParallel = ownScvf.hasParallelNeighbor(localSubFaceIdx, 0)
? faceVars.velocityParallel(localSubFaceIdx, 0) * outsideVolVars.density()
: (getParallelVelocityFromBoundary_(problem, ownScvf, lateralFace,
faceVars.velocitySelf(), localSubFaceIdx, element,
: (getParallelVelocityFromBoundary_(problem, element, fvGeometry, ownScvf, lateralFace,
faceVars.velocitySelf(), localSubFaceIdx,
lateralFaceHasDirichletPressure, lateralFaceHasBJS)
* insideVolVars.density());
......@@ -573,11 +574,12 @@ private:
* \param lateralFaceHasBJS @c true if there is a BJS condition fot the velocity on the boundary
*/
static Scalar getParallelVelocityFromBoundary_(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf,
const SubControlVolumeFace& normalFace,
const Scalar velocitySelf,
const int localSubFaceIdx,
const Element& element,
const bool lateralFaceHasDirichletPressure,
const bool lateralFaceHasBJS)
{
......@@ -587,7 +589,10 @@ private:
return velocitySelf;
if (lateralFaceHasBJS)
return problem.bjsVelocity(element, scvf, normalFace, localSubFaceIdx, velocitySelf);
{
const SubControlVolume& scv = fvGeometry.scv(scvf.insideScvIdx());
return problem.bjsVelocity(element, scv, normalFace, velocitySelf);
}
else
{
// ________________
......@@ -658,7 +663,10 @@ private:
else if (bcTypes.isSymmetry() || bcTypes.isDirichlet(Indices::pressureIdx))
return parallelVelocity;
else if (bcTypes.isBJS(Indices::velocity(scvf.directionIndex())))
return problem.bjsVelocity(boundaryElement, scvf, boundaryNormalFace, localIdx, parallelVelocity);
{
const SubControlVolume& scv = fvGeometry.scv(scvf.insideScvIdx());
return problem.bjsVelocity(boundaryElement, scv,boundaryNormalFace, parallelVelocity);
}
else
{
// Neumann conditions are not well implemented
......
......@@ -64,6 +64,7 @@ class RANSProblemBase : public NavierStokesProblem<TypeTag>
using FVElementGeometry = typename FVGridGeometry::LocalView;
using GridView = typename FVGridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
......@@ -332,21 +333,22 @@ public:
const int numSubFaces = scvf.pairData().size();
for(int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx)
{
const auto& normalFace = fvGeometry.scvf(scvf.insideScvIdx(), scvf.pairData()[localSubFaceIdx].localNormalFaceIdx);
const auto& lateralFace = fvGeometry.scvf(scvf.insideScvIdx(), scvf.pairData()[localSubFaceIdx].localLateralFaceIdx);
// 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 normalNormDim = lateralFace.directionIndex();
if (lateralFace.boundary() && (asImp_().boundaryTypes(element, lateralFace).isBJS(Indices::velocity(velIdx))))
{
unsigned int neighborIdx = neighborIdx_[elementIdx][normalNormDim][0];
if (normalFace.center()[normalNormDim] < cellCenter_[elementIdx][normalNormDim])
if (lateralFace.center()[normalNormDim] < cellCenter_[elementIdx][normalNormDim])
neighborIdx = neighborIdx_[elementIdx][normalNormDim][1];
bjsVelocityAverage[normalNormDim] += ParentType::bjsVelocity(element, scvf, normalFace, localSubFaceIdx, velocity_[elementIdx][velIdx]);
const SubControlVolume& scv = fvGeometry.scv(scvf.insideScvIdx());
bjsVelocityAverage[normalNormDim] += ParentType::bjsVelocity(element, scv, lateralFace, velocity_[elementIdx][velIdx]);
if (bjsNumFaces[normalNormDim] > 0 && neighborIdx != bjsNeighbor[normalNormDim])
DUNE_THROW(Dune::InvalidStateException, "Two different neighborIdx should not occur");
bjsNeighbor[normalNormDim] = neighborIdx;
normalNormCoordinate[normalNormDim] = normalFace.center()[normalNormDim];
normalNormCoordinate[normalNormDim] = lateralFace.center()[normalNormDim];
bjsNumFaces[normalNormDim]++;
}
}
......@@ -446,7 +448,7 @@ public:
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace& scvf,
const SubControlVolumeFace& localSubFace) const
const SubControlVolumeFace& lateralBoundaryFace) const
{ return FacePrimaryVariables(0.0); }
/*!
......
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