From b39b80bd847a8f56c58ea53d6f8126f19f5660a3 Mon Sep 17 00:00:00 2001
From: Nicolas Schwenck <nicolas.schwenck@iws.uni-stuttgart.de>
Date: Tue, 2 Dec 2014 09:28:07 +0000
Subject: [PATCH] fixed bug in potential gradient calculation with gravity in
fractures in 2pdfm.
reviewed by Alex K.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@13832 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/implicit/2pdfm/2pdfmfluxvariables.hh | 32 ++++++++++++++++++++++
1 file changed, 32 insertions(+)
diff --git a/dumux/implicit/2pdfm/2pdfmfluxvariables.hh b/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
index e5476d9b03..9bd7cbcf48 100644
--- a/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
+++ b/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
@@ -214,6 +214,38 @@ private:
potentialGradFracture_[phaseIdx] =
(elemVolVars[j].pressure(phaseIdx) - elemVolVars[i].pressure(phaseIdx))
/ diff_ij.two_norm();
+
+ // correct the pressure gradient by the gravitational acceleration
+ if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
+ {
+ // ask for the gravitational acceleration at the given SCV face
+ GlobalPosition g(problem.gravityAtPos(this->face().ipGlobal));
+
+ // calculate the phase density at the integration point. we
+ // only do this if the wetting phase is present in both cells
+ Scalar SI = elemVolVars[this->face().i].fluidState().saturation(phaseIdx);
+ Scalar SJ = elemVolVars[this->face().j].fluidState().saturation(phaseIdx);
+ Scalar rhoI = elemVolVars[this->face().i].fluidState().density(phaseIdx);
+ Scalar rhoJ = elemVolVars[this->face().j].fluidState().density(phaseIdx);
+ Scalar fI = std::max(0.0, std::min(SI/1e-5, 0.5));
+ Scalar fJ = std::max(0.0, std::min(SJ/1e-5, 0.5));
+ if (fI + fJ == 0)
+ // doesn't matter because no wetting phase is present in
+ // both cells!
+ fI = fJ = 0.5;
+ const Scalar density = (fI*rhoI + fJ*rhoJ)/(fI + fJ);
+
+ // make gravity acceleration a force
+ GlobalPosition f(g);
+ f *= density;
+
+ //transform into fracture local coordinates
+ diff_ij/=diff_ij.two_norm();
+ Scalar fractureLocalGravityForce=diff_ij*f;
+
+ // calculate the final potential gradient
+ potentialGradFracture_[phaseIdx] -= fractureLocalGravityForce;
+ }// gravity
}
else
{
--
GitLab