diff --git a/dumux/decoupled/2p/fluxData2p.hh b/dumux/decoupled/2p/fluxData2p.hh index 075402499b9d92d61f0a3eee3918a381237a83a1..f3dd3667e332be0da1136267982615b47c9bc34a 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 69304bd2a8ffa7f3efacd5f8f650b373c1dc222f..4f0da8c4267411d40b3ff330755054e9a9feca41 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 e70d63c98cfcea7cbc2ab70c48445e7de149b7d1..33faa766ea2268ea664d0b062e61e40e679fa94f 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 bfd6b6c1a09f2f95e5b049debe7190812fb0915d..ea5886e435f0a1bd9076651d268dd5ab96945663 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_;