From bdb7b01cd8f395d0f0a94d12633fcc5374eacb58 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 15 Mar 2016 02:21:41 +0100
Subject: [PATCH] [fluxes] Use average density as face density approximation
---
.../implicit/cellcentered/tpfa/darcyslaw.hh | 20 +++++++++++--------
.../implicit/cellcentered/tpfa/fickslaw.hh | 6 +++---
2 files changed, 15 insertions(+), 11 deletions(-)
diff --git a/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh b/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh
index 2667980165..75f7a891a9 100644
--- a/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh
+++ b/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh
@@ -124,20 +124,24 @@ public:
if (enableGravity_)
{
+ // do averaging for the density
+ const auto rhoInside = insideVolVars->density(phaseIdx);
+ const auto rhoOutide = outsideVolVars->density(phaseIdx);
+ const auto rho = (rhoInside + rhoOutide)*0.5;
+
+
// ask for the gravitational acceleration in the inside neighbor
const auto xInside = insideScv.center();
const auto gInside = problem_().gravityAtPos(xInside);
- const auto rhoInside = insideVolVars->density(phaseIdx);
- hInside -= rhoInside*(gInside*xInside);
+ hInside -= rho*(gInside*xInside);
// and the outside neighbor
- auto rhoOutside = outsideVolVars->density(phaseIdx);
if (scvFace_().boundary())
{
const auto xOutside = scvFace_().center();
const auto gOutside = problem_().gravityAtPos(xOutside);
- hOutside -= rhoOutside*(gOutside*xOutside);
+ hOutside -= rho*(gOutside*xOutside);
}
else
{
@@ -145,7 +149,7 @@ public:
const auto& outsideScv = problem_().model().fvGeometries().subControlVolume(outsideScvIdx);
const auto xOutside = outsideScv.center();
const auto gOutside = problem_().gravityAtPos(xOutside);
- hOutside -= rhoOutside*(gOutside*xOutside);
+ hOutside -= rho*(gOutside*xOutside);
}
}
@@ -180,7 +184,7 @@ public:
return kGradPNormal_[phaseIdx]*upwindFunction(upVolVars(phaseIdx), dnVolVars(phaseIdx));
}
- std::set<IndexType> stencil() const
+ const std::set<IndexType>& stencil() const
{
return stencil_;
}
@@ -234,10 +238,10 @@ private:
{
// fill the stencil
if (!scvFace_().boundary())
- stencil_= {scvFace_().insideScvIdx(), scvFace_().outsideScvIdx()};
+ stencil_.insert({scvFace_().insideScvIdx(), scvFace_().outsideScvIdx()});
else
// fill the stencil
- stencil_ = {scvFace_().insideScvIdx()};
+ stencil_.insert(scvFace_().insideScvIdx());
}
Scalar calculateOmega_(const DimWorldMatrix &K, const SubControlVolume &scv) const
diff --git a/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh b/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh
index 646c7b6000..be8a44f8a1 100644
--- a/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh
+++ b/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh
@@ -119,9 +119,9 @@ public:
}
// compute the diffusive flux
- auto xInside = insideVolVars->moleFraction(phaseIdx_, compIdx_);
- auto xOutside = outsideVolVars->moleFraction(phaseIdx_, compIdx_);
- auto rho = 0.5*(insideVolVars->molarDensity(phaseIdx_) + outsideVolVars->molarDensity(phaseIdx_));
+ const auto xInside = insideVolVars->moleFraction(phaseIdx_, compIdx_);
+ const auto xOutside = outsideVolVars->moleFraction(phaseIdx_, compIdx_);
+ const auto rho = 0.5*(insideVolVars->molarDensity(phaseIdx_) + outsideVolVars->molarDensity(phaseIdx_));
rhoDGradXNormal_ = rho*tij_*(xInside - xOutside);
}
--
GitLab