From 4dea8c06423d8a8fcbbaffd0b4aafbad3ebdbb4b Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 28 Jul 2015 17:24:14 +0200 Subject: [PATCH] add invertion methods to some of the two-phase fluid-matrix interactions these are required for hysteresis --- .../fluidmatrixinteractions/BrooksCorey.hpp | 49 +++++++++- .../PiecewiseLinearTwoPhaseMaterial.hpp | 91 ++++++++++++++++++- .../RegularizedBrooksCorey.hpp | 26 ++++++ 3 files changed, 160 insertions(+), 6 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/BrooksCorey.hpp b/opm/material/fluidmatrixinteractions/BrooksCorey.hpp index 8482912b0..0aec7ba2c 100644 --- a/opm/material/fluidmatrixinteractions/BrooksCorey.hpp +++ b/opm/material/fluidmatrixinteractions/BrooksCorey.hpp @@ -30,6 +30,8 @@ #include "BrooksCoreyParams.hpp" #include +#include +#include #include #include @@ -167,7 +169,17 @@ public: assert(0 <= Sw && Sw <= 1); - return params.entryPressure()*Toolbox::pow(Sw, -1.0/params.lambda()); + return params.entryPressure()*Toolbox::pow(Sw, -1/params.lambda()); + } + + template + static Evaluation twoPhaseSatPcnwInv(const Params ¶ms, const Evaluation& pcnw) + { + typedef MathToolbox Toolbox; + + assert(pcnw > 0.0); + + return Toolbox::pow(params.entryPressure()/pcnw, -params.lambda()); } /*! @@ -243,6 +255,14 @@ public: return Toolbox::pow(Sw, 2.0/params.lambda() + 3); } + template + static Evaluation twoPhaseSatKrwInv(const Params ¶ms, const Evaluation& krw) + { + typedef MathToolbox Toolbox; + + return Toolbox::pow(krw, 1.0/(2.0/params.lambda() + 3)); + } + /*! * \brief The relative permeability for the non-wetting phase of * the medium as implied by the Brooks-Corey @@ -273,6 +293,33 @@ public: const Evaluation Sn = 1. - Sw; return Sn*Sn*(1. - Toolbox::pow(Sw, exponent)); } + + template + static Evaluation twoPhaseSatKrnInv(const Params ¶ms, const Evaluation& krn) + { + typedef MathToolbox Toolbox; + + // since inverting the formula for krn is hard to do analytically, we use the + // Newton-Raphson method + Evaluation Sw = Toolbox::createConstant(0.5); + Scalar eps = 1e-10; + for (int i = 0; i < 20; ++i) { + Evaluation f = krn - twoPhaseSatKrn(params, Sw); + Evaluation fStar = krn - twoPhaseSatKrn(params, Sw + eps); + Evaluation fPrime = (fStar - f)/eps; + + Evaluation delta = f/fPrime; + Sw -= delta; + if (Sw < 0) + Sw = Toolbox::createConstant(0.0); + if (Toolbox::abs(delta) < 1e-10) + return Sw; + } + + OPM_THROW(NumericalIssue, + "Couldn't invert the Brooks-Corey non-wetting phase" + " relperm within 20 iterations"); + } }; } // namespace Opm diff --git a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp index b019aaeba..cc1b9ebae 100644 --- a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp @@ -143,6 +143,14 @@ public: static Evaluation twoPhaseSatPcnw(const Params ¶ms, const Evaluation& Sw) { return eval_(params.SwSamples(), params.pcnwSamples(), Sw); } + template + static Evaluation twoPhaseSatPcnwInv(const Params ¶ms, const Evaluation& pcnw) + { + return eval_(params.pcnwSamples(), + params.SwSamples(), + pcnw); + } + /*! * \brief The saturation-capillary pressure curve */ @@ -181,7 +189,7 @@ public: } template - static Evaluation twoPhaseSatKrw(const Params ¶ms, const Evaluation Sw) + static Evaluation twoPhaseSatKrw(const Params ¶ms, const Evaluation& Sw) { typedef MathToolbox Toolbox; @@ -189,6 +197,14 @@ public: return Toolbox::max(0.0, Toolbox::min(1.0, res)); } + template + static Evaluation twoPhaseSatKrwInv(const Params ¶ms, const Evaluation& krw) + { + return eval_(params.krwSamples(), + params.SwSamples(), + krw); + } + /*! * \brief The relative permeability for the non-wetting phase * of the porous medium @@ -213,18 +229,36 @@ public: Sw))); } + template + static Evaluation twoPhaseSatKrnInv(const Params ¶ms, const Evaluation& krn) + { + return eval_(params.krnSamples(), + params.SwSamples(), + krn); + } + private: template static Evaluation eval_(const ValueVector &xValues, const ValueVector &yValues, const Evaluation& x) + { + if (xValues.front() < xValues.back()) + return evalAscending_(xValues, yValues, x); + return evalDescending_(xValues, yValues, x); + } + + template + static Evaluation evalAscending_(const ValueVector &xValues, + const ValueVector &yValues, + const Evaluation& x) { typedef MathToolbox Toolbox; - if (Toolbox::value(x) < xValues.front()) - return Toolbox::createConstant(yValues.front()); - if (Toolbox::value(x) > xValues.back()) - return Toolbox::createConstant(yValues.back()); + if (x <= xValues.front()) + return yValues.front(); + if (x >= xValues.back()) + return yValues.back(); int segIdx = findSegmentIndex_(xValues, Toolbox::value(x)); @@ -239,6 +273,31 @@ private: return y0 + (x - x0)*m; } + template + static Evaluation evalDescending_(const ValueVector &xValues, + const ValueVector &yValues, + const Evaluation& x) + { + typedef MathToolbox Toolbox; + + if (x >= xValues.front()) + return yValues.front(); + if (x <= xValues.back()) + return yValues.back(); + + int segIdx = findSegmentIndexDescending_(xValues, Toolbox::value(x)); + + Scalar x0 = xValues[segIdx]; + Scalar x1 = xValues[segIdx + 1]; + + Scalar y0 = yValues[segIdx]; + Scalar y1 = yValues[segIdx + 1]; + + Scalar m = (y1 - y0)/(x1 - x0); + + return y0 + (x - x0)*m; + } + template static Evaluation evalDeriv_(const ValueVector &xValues, const ValueVector &yValues, @@ -283,6 +342,28 @@ private: return lowIdx; } + + static int findSegmentIndexDescending_(const ValueVector &xValues, Scalar x) + { + int n = xValues.size() - 1; + assert(n >= 1); // we need at least two sampling points! + if (xValues[n] >= x) + return n - 1; + else if (xValues[0] <= x) + return 0; + + // bisection + int lowIdx = 0, highIdx = n; + while (lowIdx + 1 < highIdx) { + int curIdx = (lowIdx + highIdx)/2; + if (xValues[curIdx] >= x) + lowIdx = curIdx; + else + highIdx = curIdx; + } + + return lowIdx; + } }; } // namespace Opm diff --git a/opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp b/opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp index 13525730d..4400e4f50 100644 --- a/opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp +++ b/opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp @@ -299,6 +299,19 @@ public: return BrooksCorey::twoPhaseSatKrw(params, Sw); } + template + static Evaluation twoPhaseSatKrwInv(const Params ¶ms, const Evaluation& krw) + { + typedef MathToolbox Toolbox; + + if (krw <= 0.0) + return Toolbox::createConstant(0.0); + else if (krw >= 1.0) + return Toolbox::createConstant(1.0); + + return BrooksCorey::twoPhaseSatKrwInv(params, krw); + } + /*! * \brief Regularized version of the relative permeability of the * non-wetting phase of the Brooks-Corey curves. @@ -335,6 +348,19 @@ public: return BrooksCorey::twoPhaseSatKrn(params, Sw); } + + template + static Evaluation twoPhaseSatKrnInv(const Params ¶ms, const Evaluation& krn) + { + typedef MathToolbox Toolbox; + + if (krn <= 0.0) + return Toolbox::createConstant(1.0); + else if (krn >= 1.0) + return Toolbox::createConstant(0.0); + + return BrooksCorey::twoPhaseSatKrnInv(params, krn); + } }; } // namespace Opm