Commit 8c8a6c87 authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[freeflow] Correctly account for gravity when using (coupling) Neumann BCs

* add source (gravity) term to residual, also when face has Neumann BCs
* correct coupling condition
parent 79df24a9
......@@ -302,17 +302,18 @@ protected:
}
else if(bcTypes.isNeumann(Indices::velocity(scvf.directionIndex())))
{
// set a given Neumann flux for the face on the boundary itself
// the source term has already been accounted for, here we
// add a given Neumann flux for the face on the boundary itself ...
const auto extrusionFactor = elemVolVars[scvf.insideScvIdx()].extrusionFactor();
residual = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())]
residual += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())]
* extrusionFactor * scvf.area();
// treat the remaining (frontal and lateral) faces of the staggered control volume
// ... and treat the fluxes of the remaining (frontal and lateral) faces of the staggered control volume
residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
}
else if(bcTypes.isSymmetry())
{
// For symmetry boundary conditions, there is no flow accross the boundary and
// for symmetry boundary conditions, there is no flow accross the boundary and
// we therefore treat it like a Dirichlet boundary conditions with zero velocity
const Scalar velocity = elemFaceVars[scvf].velocitySelf();
const Scalar fixedValue = 0.0;
......@@ -320,7 +321,8 @@ protected:
}
else if(bcTypes.isDirichlet(Indices::pressureIdx))
{
// If none of the above conditions apply, we are at an "fixed pressure" boundary for which the velocity needs to be assembled
// 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
// as if it where inside the domain and not on the boundary (source term has already been acounted for)
residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
}
else
......
......@@ -294,9 +294,14 @@ public:
}
else // use pressure reconstruction for single phase models
{
// v = -K/mu * (gradP + rho*g)
const Scalar velocity = stokesElemFaceVars[scvf].velocitySelf();
const Scalar pressureInterFace = scvf.directionSign() * velocity * stokesContext.volVars.viscosity(darcyPhaseIdx)/darcyPermeability(scvf) * (stokesContext.element.geometry().center() - scvf.center()).two_norm() + darcyPressure;
momentumFlux = pressureInterFace;
const Scalar mu = stokesContext.volVars.viscosity(darcyPhaseIdx);
const Scalar rho = stokesContext.volVars.density(darcyPhaseIdx);
const Scalar distance = (stokesContext.element.geometry().center() - scvf.center()).two_norm();
const Scalar g = couplingManager_.problem(darcyIdx).gravity()[scvf.directionIndex()];
const Scalar interfacePressure = ((scvf.directionSign() * velocity * (mu/darcyPermeability(scvf))) + rho * g) * distance + darcyPressure;
momentumFlux = interfacePressure;
}
// normalize pressure
......
......@@ -30,6 +30,8 @@
#include <utility>
#include <memory>
#include <dune/common/float_cmp.hh>
#include <dune/common/exceptions.hh>
#include <dumux/common/properties.hh>
#include <dumux/multidomain/staggeredcouplingmanager.hh>
......@@ -119,7 +121,6 @@ public:
using ParentType::couplingStencil;
using ParentType::updateCouplingContext;
using DarcyGridVariables = GridVariables<2>;
using CouplingData = StokesDarcyCouplingData<MDTraits, StokesDarcyCouplingManager<MDTraits>>;
//! Constructor
......@@ -137,6 +138,9 @@ public:
std::shared_ptr<const Problem<darcyIdx>> darcyProblem,
const SolutionVector& curSol)
{
if(Dune::FloatCmp::ne(stokesProblem->gravity(), darcyProblem->gravity()))
DUNE_THROW(Dune::InvalidStateException, "Both models must use the same gravity vector");
ParentType::init(std::make_tuple(stokesProblem, stokesProblem, darcyProblem));
this->curSol() = curSol;
couplingData_ = std::make_shared<CouplingData>(*this);
......
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