From 35eeebb33714efb8250e2a3429398f113c9bd94a Mon Sep 17 00:00:00 2001
From: Markus Wolff <markus.wolff@twt-gmbh.de>
Date: Mon, 3 Sep 2012 14:36:46 +0000
Subject: [PATCH] improved decoupled 2p model
- correct treatment of mixed boundary conditions with a total
velocity formulation if gravity and capillarity are enabled
- reviewed by Yufei
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@9028 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/decoupled/2p/fluxData2p.hh | 4 +-
.../2p/transport/fv/capillarydiffusion.hh | 18 +++++++-
.../decoupled/2p/transport/fv/gravitypart.hh | 41 +++++++++++++++----
dumux/decoupled/common/impetproblem.hh | 4 --
4 files changed, 51 insertions(+), 16 deletions(-)
diff --git a/dumux/decoupled/2p/fluxData2p.hh b/dumux/decoupled/2p/fluxData2p.hh
index 075402499b..f3dd3667e3 100644
--- a/dumux/decoupled/2p/fluxData2p.hh
+++ b/dumux/decoupled/2p/fluxData2p.hh
@@ -201,7 +201,7 @@ public:
*/
bool isUpwindCell(int phaseIdx, int indexInInside)
{
- return (potential_[indexInInside][phaseIdx] >= 0.);
+ return (potential_[indexInInside][phaseIdx] > 0.);
}
/*! \brief Checks for upwind direction
@@ -212,7 +212,7 @@ public:
*/
bool isUpwindCell(int phaseIdx, int indexInInside) const
{
- return (potential_[indexInInside][phaseIdx] >= 0.);
+ return (potential_[indexInInside][phaseIdx] > 0.);
}
/*! \brief Returns the phase potential at a cell-cell interface
diff --git a/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh b/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh
index 69304bd2a8..4f0da8c426 100644
--- a/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh
+++ b/dumux/decoupled/2p/transport/fv/capillarydiffusion.hh
@@ -58,6 +58,9 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
enum
{
@@ -65,7 +68,8 @@ private:
};
enum
{
- wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
+ pressEqIdx = Indices::pressEqIdx
};
typedef typename GridView::Traits::template Codim<0>::Entity Element;
@@ -176,6 +180,18 @@ public:
}//end intersection with neighbor
else
{
+ BoundaryTypes bcTypes;
+ problem_.boundaryTypes(bcTypes, intersection);
+ if (bcTypes.isNeumann(pressEqIdx))
+ {
+ PrimaryVariables priVars;
+ problem_.neumann(priVars, intersection);
+ if (priVars[wPhaseIdx] == 0)
+ {
+ flux = 0;
+ return;
+ }
+ }
// get permeability
problem_.spatialParams().meanK(meanPermeability,
problem_.spatialParams().intrinsicPermeability(*element));
diff --git a/dumux/decoupled/2p/transport/fv/gravitypart.hh b/dumux/decoupled/2p/transport/fv/gravitypart.hh
index e70d63c98c..33faa766ea 100644
--- a/dumux/decoupled/2p/transport/fv/gravitypart.hh
+++ b/dumux/decoupled/2p/transport/fv/gravitypart.hh
@@ -114,7 +114,13 @@ public:
lambdaNWI /= viscosity_[nPhaseIdx];
}
+ Scalar potentialW = cellDataI.fluxData().potential(wPhaseIdx, indexInInside);
+ Scalar potentialNW = cellDataI.fluxData().potential(nPhaseIdx, indexInInside);
+
DimMatrix meanPermeability(0);
+ GlobalPosition distVec(0);
+ Scalar lambdaW = 0;
+ Scalar lambdaNW = 0;
if (intersection.neighbor())
{
@@ -124,6 +130,8 @@ public:
int globalIdxJ = problem_.variables().index(*neighborPointer);
CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
+ distVec = neighborPointer->geometry().center() - element->geometry().center();
+
// get permeability
problem_.spatialParams().meanK(meanPermeability,
problem_.spatialParams().intrinsicPermeability(*element),
@@ -142,6 +150,11 @@ public:
lambdaNWJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*neighborPointer), satJ);
lambdaNWJ /= viscosity_[nPhaseIdx];
}
+
+ lambdaW = (potentialW >= 0) ? lambdaWI : lambdaWJ;
+ lambdaW = (potentialW == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
+ lambdaNW = (potentialNW >= 0) ? lambdaNWI : lambdaNWJ;
+ lambdaNW = (potentialNW == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNW;
}
else
{
@@ -149,27 +162,37 @@ public:
problem_.spatialParams().meanK(meanPermeability,
problem_.spatialParams().intrinsicPermeability(*element));
+ distVec = intersection.geometry().center() - element->geometry().center();
+
//calculate lambda_n*f_w at the boundary
lambdaWJ = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*element), satJ);
lambdaWJ /= viscosity_[wPhaseIdx];
lambdaNWJ = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*element), satJ);
lambdaNWJ /= viscosity_[nPhaseIdx];
+
+ //If potential is zero always take value from the boundary!
+ lambdaW = (potentialW > 0) ? lambdaWI : lambdaWJ;
+ lambdaNW = (potentialNW > 0) ? lambdaNWI : lambdaNWJ;
}
// set result to K*grad(pc)
- meanPermeability.mv(problem_.gravity(), flux);
+ const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
+ Scalar dist = distVec.two_norm();
+ //calculate unit distVec
+ distVec /= dist;
+ Scalar areaScaling = (unitOuterNormal * distVec);
- Scalar potentialW = cellDataI.fluxData().potential(wPhaseIdx, indexInInside);
- Scalar potentialNW = cellDataI.fluxData().potential(nPhaseIdx, indexInInside);
+ Dune::FieldVector<Scalar, dim> permeability(0);
+ meanPermeability.mv(unitOuterNormal, permeability);
+
+ Scalar scalarPerm = permeability.two_norm();
+
+ Scalar scalarGravity = problem_.gravity() * distVec;
- Scalar lambdaW = (potentialW >= 0) ? lambdaWI : lambdaWJ;
- lambdaW = (potentialW == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
- Scalar lambdaNW = (potentialNW >= 0) ? lambdaNWI : lambdaNWJ;
- lambdaNW = (potentialNW == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNW;
+ flux = unitOuterNormal;
// set result to f_w*lambda_n*K*grad(pc)
- flux *= lambdaW*lambdaNW/(lambdaW+lambdaNW);
- flux *= (density_[wPhaseIdx] - density_[nPhaseIdx]);
+ flux *= lambdaW*lambdaNW/(lambdaW+lambdaNW) * scalarPerm * (density_[wPhaseIdx] - density_[nPhaseIdx]) * scalarGravity * areaScaling;
}
/*! \brief Constructs a GravityPart object
*
diff --git a/dumux/decoupled/common/impetproblem.hh b/dumux/decoupled/common/impetproblem.hh
index bfd6b6c1a0..ea5886e435 100644
--- a/dumux/decoupled/common/impetproblem.hh
+++ b/dumux/decoupled/common/impetproblem.hh
@@ -102,7 +102,6 @@ public:
bboxMin_(std::numeric_limits<double>::max()),
bboxMax_(-std::numeric_limits<double>::max()),
timeManager_(&timeManager),
- deleteTimeManager_(false),
variables_(gridView),
outputInterval_(1),
outputTimeInterval_(0.0)
@@ -143,8 +142,6 @@ public:
delete transportModel_;
delete model_;
delete resultWriter_;
- if (deleteTimeManager_)
- delete timeManager_;
if (adaptiveGrid)
delete gridAdapt_;
}
@@ -813,7 +810,6 @@ private:
GlobalPosition bboxMax_;
TimeManager *timeManager_;
- bool deleteTimeManager_;
Variables variables_;
--
GitLab