Commit de783df2 authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Melanie Lipp
Browse files

[stagged][fluxVars] Factor out the contribution of inflow/outflow BCs

* this should make it easier to properly reconstruct a stress at the boundary
parent 00518b70
......@@ -242,12 +242,6 @@ public:
// (pointing in opposite direction of the scvf's one).
frontalFlux += pressure * -1.0 * scvf.directionSign();
// Handle inflow or outflow conditions.
// Treat the staggered half-volume adjacent to the boundary as if it was on the opposite side of the boundary.
// The respective face's outer normal vector will point in the same direction as the scvf's one.
if(scvf.boundary() && problem.boundaryTypes(element, scvf).isDirichlet(Indices::pressureIdx))
frontalFlux += inflowOutflowBoundaryFlux_(problem, element, scvf, elemVolVars, elemFaceVars);
// Account for the staggered face's area. For rectangular elements, this equals the area of the scvf
// our velocity dof of interest lives on.
return frontalFlux * scvf.area() * insideVolVars.extrusionFactor();
......@@ -375,6 +369,53 @@ public:
return lateralFlux;
}
/*!
* \brief Returns the momentum flux over an inflow or outflow boundary face.
*
* \verbatim
* scvf //
* ---------=======// == and # staggered half-control-volume
* | || #// current scvf
* | || #// # staggered boundary face over which fluxes are calculated
* | || x~~~~> vel.Self
* | || #// x dof position
* scvf | || #//
* --------========// -- element
* scvf //
* // boundary
* \endverbatim
*/
FacePrimaryVariables inflowOutflowBoundaryFlux(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars) const
{
FacePrimaryVariables inOrOutflow(0.0);
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
// Advective momentum flux.
if (problem.enableInertiaTerms())
{
const Scalar velocitySelf = elemFaceVars[scvf].velocitySelf();
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
const auto& upVolVars = (scvf.directionSign() == sign(velocitySelf)) ?
insideVolVars : outsideVolVars;
inOrOutflow += velocitySelf * velocitySelf * upVolVars.density();
}
// Apply a pressure at the boundary.
const Scalar boundaryPressure = normalizePressure
? (problem.dirichlet(element, scvf)[Indices::pressureIdx] -
problem.initial(scvf)[Indices::pressureIdx])
: problem.dirichlet(element, scvf)[Indices::pressureIdx];
inOrOutflow += boundaryPressure;
// Account for the orientation of the face at the boundary,
return inOrOutflow * scvf.directionSign() * scvf.area() * insideVolVars.extrusionFactor();
}
private:
/*!
......@@ -577,54 +618,6 @@ private:
return lateralDiffusiveFlux * lateralFace.area() * 0.5 * extrusionFactor_(elemVolVars, lateralFace);
}
/*!
* \brief Returns the momentum flux over an inflow or outflow boundary face.
*
* \verbatim
* scvf //
* ---------=======// == and # staggered half-control-volume
* | || #// current scvf
* | || #// # staggered boundary face over which fluxes are calculated
* | || x~~~~> vel.Self
* | || #// x dof position
* scvf | || #//
* --------========// -- element
* scvf //
* // boundary
* \endverbatim
*/
FacePrimaryVariables inflowOutflowBoundaryFlux_(const Problem& problem,
const Element& element,
const SubControlVolumeFace& scvf,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars)
{
FacePrimaryVariables inOrOutflow(0.0);
// Advective momentum flux.
if (problem.enableInertiaTerms())
{
const Scalar velocitySelf = elemFaceVars[scvf].velocitySelf();
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
const auto& upVolVars = (scvf.directionSign() == sign(velocitySelf)) ?
insideVolVars : outsideVolVars;
inOrOutflow += velocitySelf * velocitySelf * upVolVars.density();
}
// Apply a pressure at the boundary.
const Scalar boundaryPressure = normalizePressure
? (problem.dirichlet(element, scvf)[Indices::pressureIdx] -
problem.initial(scvf)[Indices::pressureIdx])
: problem.dirichlet(element, scvf)[Indices::pressureIdx];
inOrOutflow += boundaryPressure;
// Account for the orientation of the face at the boundary,
return inOrOutflow * scvf.directionSign();
}
private:
/*!
* \brief Get the location of the lateral staggered face's center.
......
......@@ -279,6 +279,8 @@ public:
if (scvf.boundary())
{
FluxVariables fluxVars;
// handle the actual boundary conditions:
const auto bcTypes = problem.boundaryTypes(element, scvf);
if (bcTypes.isNeumann(Indices::velocity(scvf.directionIndex())))
......@@ -290,13 +292,16 @@ public:
* extrusionFactor * scvf.area();
// ... and treat the fluxes of the remaining (frontal and lateral) faces of the staggered control volume
result += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
result += fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache());
}
else if(bcTypes.isDirichlet(Indices::pressureIdx))
{
// if none of the above conditions apply, we are at an "fixed pressure" boundary for which the resdiual of the momentum balance needs to be assembled
// we are at an "fixed pressure" boundary for which the resdiual of the momentum balance needs to be assembled
// as if it where inside the domain and not on the boundary (source term has already been acounted for)
result = computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
result = fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache());
// incorporate the inflow or outflow contribution
result += fluxVars.inflowOutflowBoundaryFlux(problem, element, scvf, elemVolVars, elemFaceVars);
}
}
return result;
......
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