From 9e9e2facb0c8d9e1078c59d90a62ba6cb3db7a9f Mon Sep 17 00:00:00 2001 From: Maziar Veyskarami Date: Thu, 3 Feb 2022 14:02:48 +0100 Subject: [PATCH 1/5] [material][fluidmatrixinteraction] modify the pcSnapOff for porenetwork model to acount the corner half angle of the throat --- .../throat/thresholdcapillarypressures.hh | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/dumux/material/fluidmatrixinteractions/porenetwork/throat/thresholdcapillarypressures.hh b/dumux/material/fluidmatrixinteractions/porenetwork/throat/thresholdcapillarypressures.hh index a79caf04c9..17ca9e6888 100644 --- a/dumux/material/fluidmatrixinteractions/porenetwork/throat/thresholdcapillarypressures.hh +++ b/dumux/material/fluidmatrixinteractions/porenetwork/throat/thresholdcapillarypressures.hh @@ -35,14 +35,18 @@ struct ThresholdCapillaryPressures template static constexpr Scalar pcSnapoff(const Scalar surfaceTension, const Scalar contactAngle, - const Scalar inscribedRadius) noexcept + const Scalar inscribedRadius, + const Scalar cornerHalfAngle) noexcept { using std::sin; using std::cos; const Scalar theta = contactAngle; - const Scalar cosTheta = std::cos(theta); - const Scalar sinTheta = std::sin(theta); - return surfaceTension / inscribedRadius * (cosTheta - sinTheta); + const Scalar cosTheta = cos(theta); + const Scalar sinTheta = sin(theta); + const Scalar Beta = cornerHalfAngle; + const Scalar cosBeta = cos(Beta); + const Scalar sinBeta = sin(Beta); + return surfaceTension / inscribedRadius * (cosTheta * cosBeta/sinBeta - sinTheta); } /*! \brief The entry capillary pressure of a pore throat. -- GitLab From 3c9b20b0958154721060b5a77e4ac6ba9ae58cac Mon Sep 17 00:00:00 2001 From: Maziar Veyskarami Date: Thu, 3 Feb 2022 14:04:33 +0100 Subject: [PATCH 2/5] [material][spatialparams] modify the 2p spatialparams for pnm to include the cornerhalfangle in pcSnapOff] --- .../spatialparams/porenetwork/porenetwork2p.hh | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/dumux/material/spatialparams/porenetwork/porenetwork2p.hh b/dumux/material/spatialparams/porenetwork/porenetwork2p.hh index dfa9788623..b9ab69ff2d 100644 --- a/dumux/material/spatialparams/porenetwork/porenetwork2p.hh +++ b/dumux/material/spatialparams/porenetwork/porenetwork2p.hh @@ -149,11 +149,19 @@ public: template const Scalar pcSnapoff(const Element& element, const ElementVolumeVariables& elemVolVars) const { - // take the average of both adjacent pores TODO: is this correct? - const Scalar surfaceTension = 0.5*(elemVolVars[0].surfaceTension() + elemVolVars[1].surfaceTension()); - return ThresholdCapillaryPressures::pcSnapoff(surfaceTension, - this->asImp_().contactAngle(element, elemVolVars), - this->asImp_().throatInscribedRadius(element, elemVolVars)); + const auto eIdx = this->gridGeometry().elementMapper().index(element); + const auto shape = this->gridGeometry().throatCrossSectionShape(eIdx); + if(shape == Throat::Shape::circle) + return -1.0; // Snap-off doesn't happen in a circular throat + else + { + // take the average of both adjacent pores TODO: is this correct? + const Scalar surfaceTension = 0.5*(elemVolVars[0].surfaceTension() + elemVolVars[1].surfaceTension()); + return ThresholdCapillaryPressures::pcSnapoff(surfaceTension, + this->asImp_().contactAngle(element, elemVolVars), + this->asImp_().throatInscribedRadius(element, elemVolVars), + cornerHalfAngles(element)[0]); + } } /*! -- GitLab From 40549ec3a92e12d34d5f8be0985d2c4ce41a7f74 Mon Sep 17 00:00:00 2001 From: Maziar Veyskarami Date: Tue, 8 Feb 2022 10:42:07 +0100 Subject: [PATCH 3/5] [material][spatialparams] set the lowest neagtive value and comment on it for circular throat --- .../material/spatialparams/porenetwork/porenetwork2p.hh | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/dumux/material/spatialparams/porenetwork/porenetwork2p.hh b/dumux/material/spatialparams/porenetwork/porenetwork2p.hh index b9ab69ff2d..655e2e0141 100644 --- a/dumux/material/spatialparams/porenetwork/porenetwork2p.hh +++ b/dumux/material/spatialparams/porenetwork/porenetwork2p.hh @@ -151,8 +151,13 @@ public: { const auto eIdx = this->gridGeometry().elementMapper().index(element); const auto shape = this->gridGeometry().throatCrossSectionShape(eIdx); - if(shape == Throat::Shape::circle) - return -1.0; // Snap-off doesn't happen in a circular throat + + // Snap-off doesn't happen in a circular throat. + // To make sure that snap-off never happens in case of circular throat, + // i.e. capillary pressure in the pore always be greater than pc snap-off, + // it returns the lowest negative value. + if (shape == Throat::Shape::circle) + return std::numeric_limits::lowest(); else { // take the average of both adjacent pores TODO: is this correct? -- GitLab From a94bb54e342f6ca7d1931ca6937fbb85d3a7f4e0 Mon Sep 17 00:00:00 2001 From: Maziar Veyskarami Date: Tue, 8 Feb 2022 10:43:55 +0100 Subject: [PATCH 4/5] [material][fluidmatrixinteraction] set corner half angle of square for pc snap-off in a throat and clean up --- .../porenetwork/throat/thresholdcapillarypressures.hh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dumux/material/fluidmatrixinteractions/porenetwork/throat/thresholdcapillarypressures.hh b/dumux/material/fluidmatrixinteractions/porenetwork/throat/thresholdcapillarypressures.hh index 17ca9e6888..5e01f27107 100644 --- a/dumux/material/fluidmatrixinteractions/porenetwork/throat/thresholdcapillarypressures.hh +++ b/dumux/material/fluidmatrixinteractions/porenetwork/throat/thresholdcapillarypressures.hh @@ -36,16 +36,16 @@ struct ThresholdCapillaryPressures static constexpr Scalar pcSnapoff(const Scalar surfaceTension, const Scalar contactAngle, const Scalar inscribedRadius, - const Scalar cornerHalfAngle) noexcept + const Scalar cornerHalfAngle = 0.25*M_PI /*corner half angle of square*/) noexcept { using std::sin; using std::cos; const Scalar theta = contactAngle; const Scalar cosTheta = cos(theta); const Scalar sinTheta = sin(theta); - const Scalar Beta = cornerHalfAngle; - const Scalar cosBeta = cos(Beta); - const Scalar sinBeta = sin(Beta); + const Scalar beta = cornerHalfAngle; + const Scalar cosBeta = cos(beta); + const Scalar sinBeta = sin(beta); return surfaceTension / inscribedRadius * (cosTheta * cosBeta/sinBeta - sinTheta); } -- GitLab From 38bd239bebdc2e900420037a80cc6c73eb928ff2 Mon Sep 17 00:00:00 2001 From: Maziar Veyskarami Date: Wed, 9 Feb 2022 08:48:23 +0100 Subject: [PATCH 5/5] [material][fluidmatrixinteraction] fix the pc snap-off relation and add a reference --- .../porenetwork/throat/thresholdcapillarypressures.hh | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/dumux/material/fluidmatrixinteractions/porenetwork/throat/thresholdcapillarypressures.hh b/dumux/material/fluidmatrixinteractions/porenetwork/throat/thresholdcapillarypressures.hh index 5e01f27107..012b315ee1 100644 --- a/dumux/material/fluidmatrixinteractions/porenetwork/throat/thresholdcapillarypressures.hh +++ b/dumux/material/fluidmatrixinteractions/porenetwork/throat/thresholdcapillarypressures.hh @@ -31,7 +31,10 @@ namespace Dumux::PoreNetwork { struct ThresholdCapillaryPressures { - //! The snap-off capillary pressure of a pore throat + /*! \brief The snap-off capillary pressure of a pore throat. + * For details, see Eq.4.8 in + * Blunt, Martin J. "Multiphase Flow in Permeable Media: A Pore-scale Perspective." + */ template static constexpr Scalar pcSnapoff(const Scalar surfaceTension, const Scalar contactAngle, @@ -46,7 +49,7 @@ struct ThresholdCapillaryPressures const Scalar beta = cornerHalfAngle; const Scalar cosBeta = cos(beta); const Scalar sinBeta = sin(beta); - return surfaceTension / inscribedRadius * (cosTheta * cosBeta/sinBeta - sinTheta); + return surfaceTension / inscribedRadius * (cosTheta - sinTheta * sinBeta/cosBeta); } /*! \brief The entry capillary pressure of a pore throat. -- GitLab