From 2d9cff730988e52d483b9bd983de97e952fbf816 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Thu, 10 Nov 2016 15:48:13 +0100
Subject: [PATCH] [Unsat] Added the calculation of the derivative of relative
permeability
with respect to saturation.
---
.../NonWettingPhaseBrookCoreyOilGas.cpp | 14 ++++++-
.../NonWettingPhaseBrookCoreyOilGas.h | 4 ++
.../NonWettingPhaseVanGenuchten.cpp | 15 ++++++-
.../NonWettingPhaseVanGenuchten.h | 4 ++
.../RelativePermeability.h | 4 ++
.../RelativePermeabilityCurve.h | 10 +++++
.../WettingPhaseBrookCoreyOilGas.cpp | 9 ++++-
.../WettingPhaseBrookCoreyOilGas.h | 4 ++
.../WettingPhaseVanGenuchten.cpp | 16 +++++++-
.../WettingPhaseVanGenuchten.h | 4 ++
.../TestRelativePermeabilityModel.cpp | 39 +++++++++++++++++++
11 files changed, 119 insertions(+), 4 deletions(-)
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseBrookCoreyOilGas.cpp b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseBrookCoreyOilGas.cpp
index 5a6d09da797..f81ba6d8b7b 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseBrookCoreyOilGas.cpp
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseBrookCoreyOilGas.cpp
@@ -28,7 +28,19 @@ double NonWettingPhaseBrookCoreyOilGas::getValue(
const double Se = (S - _Sr) / (_Smax - _Sr);
const double krel =
(1.0 - Se) * (1.0 - Se) * (1.0 - std::pow(Se, 1.0 + 2.0 / _mm));
- return krel < _krel_min ? _krel_min : krel;
+ return std::max(_krel_min, krel);
+}
+
+double NonWettingPhaseBrookCoreyOilGas::getdValue(
+ const double saturation_w) const
+{
+ const double S = MathLib::limitValueInInterval(
+ saturation_w, _Sr + _minor_offset, _Smax - _minor_offset);
+ const double Se = (S - _Sr) / (_Smax - _Sr);
+ return (-2. * (1.0 - Se) * (1.0 - std::pow(Se, 1.0 + 2.0 / _mm)) -
+ (1.0 + 2.0 / _mm) * (1.0 - Se) * (1.0 - Se) *
+ std::pow(Se, 2.0 / _mm)) /
+ (_Smax - _Sr);
}
} // end namespace
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseBrookCoreyOilGas.h b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseBrookCoreyOilGas.h
index 32a91922a79..9b34f2e54f9 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseBrookCoreyOilGas.h
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseBrookCoreyOilGas.h
@@ -60,6 +60,10 @@ public:
/// \param saturation_w Wetting phase saturation
double getValue(const double saturation_w) const override;
+ /// Get the derivative of relative permeability with respect to saturation.
+ /// \param saturation_w Wetting phase saturation
+ double getdValue(const double saturation_w) const override;
+
private:
const double _Sr; ///< Residual saturation of wetting phase, 1-Snmax.
const double _Smax; ///< Maximum saturation of wetting phase, 1-Snr.
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseVanGenuchten.cpp b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseVanGenuchten.cpp
index 435c52dc387..f3d397ad26a 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseVanGenuchten.cpp
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseVanGenuchten.cpp
@@ -27,7 +27,20 @@ double NonWettingPhaseVanGenuchten::getValue(const double saturation_w) const
const double Se = (S - _Sr) / (_Smax - _Sr);
const double krel = std::cbrt(1.0 - Se) *
std::pow(1.0 - std::pow(Se, 1.0 / _mm), 2.0 * _mm);
- return krel < _krel_min ? _krel_min : krel;
+ return std::max(_krel_min, krel);
+}
+
+double NonWettingPhaseVanGenuchten::getdValue(const double saturation_w) const
+{
+ const double S = MathLib::limitValueInInterval(
+ saturation_w, _Sr + _minor_offset, _Smax - _minor_offset);
+ const double Se = (S - _Sr) / (_Smax - _Sr);
+ const double cbrt1_Se = std::cbrt(1.0 - Se);
+ const double temp_val = 1.0 - std::pow(Se, 1.0 / _mm);
+ return (-std::pow(temp_val, 2. * _mm) / (3. * cbrt1_Se * cbrt1_Se) -
+ 2. * cbrt1_Se * std::pow(temp_val, 2. * _mm - 1.) *
+ std::pow(Se, (1. - _mm) / _mm)) /
+ (_Smax - _Sr);
}
} // end namespace
} // end namespace
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseVanGenuchten.h b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseVanGenuchten.h
index f24ecd4d22c..a039616fdb6 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseVanGenuchten.h
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseVanGenuchten.h
@@ -62,6 +62,10 @@ public:
/// \param saturation_w Wetting phase saturation
double getValue(const double saturation_w) const override;
+ /// Get the derivative of relative permeability with respect to saturation.
+ /// \param saturation_w Wetting phase saturation
+ double getdValue(const double saturation_w) const override;
+
private:
const double _Sr; ///< Residual saturation of wetting phase, 1-Snmax.
const double _Smax; ///< Maximum saturation of wetting phase, 1-Snr.
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h
index 0b62d8ac50b..da69bb2a25b 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h
@@ -32,6 +32,10 @@ public:
/// \param saturation Wetting phase saturation
virtual double getValue(const double saturation) const = 0;
+ /// Get the derivative of relative permeability with respect to saturation.
+ /// \param saturation Wetting phase saturation
+ virtual double getdValue(const double saturation) const = 0;
+
protected:
/** A small number for an offset:
* 1. to set the bound of S, the saturation, such that
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeabilityCurve.h b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeabilityCurve.h
index 8232932a5ef..b89b4c7524e 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeabilityCurve.h
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeabilityCurve.h
@@ -49,6 +49,16 @@ public:
return _curve_data->getValue(S);
}
+ /// Get the derivative of relative permeability with respect to saturation.
+ /// \param saturation Wetting phase saturation
+ double getdValue(const double saturation) const override
+ {
+ const double S = MathLib::limitValueInInterval(
+ saturation, _Sr + _minor_offset, _Smax - _minor_offset);
+
+ return _curve_data->getDerivative(S);
+ }
+
private:
const double _Sr; ///< Residual saturation.
const double _Smax; ///< Maximum saturation.
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseBrookCoreyOilGas.cpp b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseBrookCoreyOilGas.cpp
index 938bf57d250..cb8cebc589a 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseBrookCoreyOilGas.cpp
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseBrookCoreyOilGas.cpp
@@ -26,8 +26,15 @@ double WettingPhaseBrookCoreyOilGas::getValue(const double saturation) const
saturation, _Sr + _minor_offset, _Smax - _minor_offset);
const double Se = (S - _Sr) / (_Smax - _Sr);
const double krel = std::pow(Se, 3.0 + 2.0 / _mm);
- return krel < _krel_min ? _krel_min : krel;
+ return std::max(_krel_min, krel);
}
+double WettingPhaseBrookCoreyOilGas::getdValue(const double saturation) const
+{
+ const double S = MathLib::limitValueInInterval(
+ saturation, _Sr + _minor_offset, _Smax - _minor_offset);
+ const double Se = (S - _Sr) / (_Smax - _Sr);
+ return ((3.0 + 2.0 / _mm) * std::pow(Se, 2.0 + 2.0 / _mm)) / (_Smax - _Sr);
+}
} // end namespace
} // end namespace
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseBrookCoreyOilGas.h b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseBrookCoreyOilGas.h
index 28b095d9982..93283367763 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseBrookCoreyOilGas.h
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseBrookCoreyOilGas.h
@@ -56,6 +56,10 @@ public:
/// Get relative permeability value.
double getValue(const double saturation) const override;
+ /// Get the derivative of relative permeability with respect to saturation.
+ /// \param saturation Wetting phase saturation
+ double getdValue(const double saturation) const override;
+
private:
const double _Sr; ///< Residual saturation.
const double _Smax; ///< Maximum saturation.
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseVanGenuchten.cpp b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseVanGenuchten.cpp
index fdcb70908ac..e1429bce38b 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseVanGenuchten.cpp
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseVanGenuchten.cpp
@@ -27,7 +27,21 @@ double WettingPhaseVanGenuchten::getValue(const double saturation) const
const double Se = (S - _Sr) / (_Smax - _Sr);
const double val = 1.0 - std::pow(1.0 - std::pow(Se, 1.0 / _mm), _mm);
const double krel = std::sqrt(Se) * val * val;
- return krel < _krel_min ? _krel_min : krel;
+ return std::max(_krel_min, krel);
+}
+
+double WettingPhaseVanGenuchten::getdValue(const double saturation) const
+{
+ const double S = MathLib::limitValueInInterval(
+ saturation, _Sr + _minor_offset, _Smax - _minor_offset);
+ const double Se = (S - _Sr) / (_Smax - _Sr);
+ const double sqrtSe = std::sqrt(Se);
+ const double temp_val = 1.0 - std::pow(1.0 - std::pow(Se, 1.0 / _mm), _mm);
+ return (0.5 * temp_val * temp_val / sqrtSe +
+ 2. * sqrtSe * temp_val *
+ std::pow(1.0 - std::pow(Se, 1.0 / _mm), _mm - 1.) *
+ std::pow(Se, (1.0 - _mm) / _mm)) /
+ (_Smax - _Sr);
}
} // end namespace
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseVanGenuchten.h b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseVanGenuchten.h
index 8fd580f378e..c69b0850b8d 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseVanGenuchten.h
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseVanGenuchten.h
@@ -59,6 +59,10 @@ public:
/// Get relative permeability value.
double getValue(const double saturation) const override;
+ /// Get the derivative of relative permeability with respect to saturation.
+ /// \param saturation Wetting phase saturation
+ double getdValue(const double saturation) const override;
+
private:
const double _Sr; ///< Residual saturation.
const double _Smax; ///< Maximum saturation.
diff --git a/Tests/MaterialLib/TestRelativePermeabilityModel.cpp b/Tests/MaterialLib/TestRelativePermeabilityModel.cpp
index ae340266fab..7612ed4d8e0 100644
--- a/Tests/MaterialLib/TestRelativePermeabilityModel.cpp
+++ b/Tests/MaterialLib/TestRelativePermeabilityModel.cpp
@@ -58,9 +58,16 @@ TEST(MaterialPorousMedium, checkWettingPhaseVanGenuchten)
0.0074049488610241927,
0.13546615958442240};
+ const double purterbation = 1.e-9;
for (std::size_t i = 0; i < S.size(); i++)
{
ASSERT_NEAR(krel[i], perm_model->getValue(S[i]), 1.e-9);
+
+ // Compare the derivative with numerical one.
+ ASSERT_NEAR((perm_model->getValue(S[i] + purterbation) -
+ perm_model->getValue(S[i])) /
+ purterbation,
+ perm_model->getdValue(S[i]), 1.e-6);
}
}
@@ -81,9 +88,16 @@ TEST(MaterialPorousMedium, checkNonWettingPhaseVanGenuchten)
0.59527539448807487, 0.49976666464188485,
0.38520070797257489, 0.041219134248319585};
+ const double purterbation = 1.e-9;
for (std::size_t i = 0; i < S.size(); i++)
{
ASSERT_NEAR(krel[i], perm_model->getValue(S[i]), 1.e-9);
+
+ // Compare the derivative with numerical one.
+ ASSERT_NEAR((perm_model->getValue(S[i] + purterbation) -
+ perm_model->getValue(S[i])) /
+ purterbation,
+ perm_model->getdValue(S[i]), 1.e-6);
}
}
@@ -107,9 +121,20 @@ TEST(MaterialPorousMedium, checkWettingPhaseBrookCoreyOilGas)
0.19753086419753069,
1.0};
+ const double purterbation = 1.e-9;
for (std::size_t i = 0; i < S.size(); i++)
{
ASSERT_NEAR(krel[i], perm_model->getValue(S[i]), 1.e-9);
+
+ // i == S.size() -1, S(=0.85) is limited to 0.8, and the numerical
+ // derivative of krel is unavailable.
+ if (i == S.size() - 1)
+ continue;
+ // Compare the derivative with numerical one.
+ ASSERT_NEAR((perm_model->getValue(S[i] + purterbation) -
+ perm_model->getValue(S[i])) /
+ purterbation,
+ perm_model->getdValue(S[i]), 1.e-6);
}
}
@@ -133,9 +158,16 @@ TEST(MaterialPorousMedium, checkNonWettingPhaseBrookCoreyOilGas)
0.061728395061728412,
.0};
+ const double purterbation = 1.e-9;
for (std::size_t i = 0; i < S.size(); i++)
{
ASSERT_NEAR(krel[i], perm_model->getValue(S[i]), 1.e-9);
+
+ // Compare the derivative with numerical one.
+ ASSERT_NEAR((perm_model->getValue(S[i] + purterbation) -
+ perm_model->getValue(S[i])) /
+ purterbation,
+ perm_model->getdValue(S[i]), 1.e-6);
}
}
@@ -155,8 +187,15 @@ TEST(MaterialPorousMedium, checkReletivePermeabilityCurve)
std::vector<double> S = {0.2, 0.33, 0.45, 0.52, 0.6, 0.85};
std::vector<double> krel = {0.7, 0.57, 0.451, 0.3824, 0.304, 0.059};
+ const double purterbation = 1.e-9;
for (std::size_t i = 0; i < S.size(); i++)
{
ASSERT_NEAR(krel[i], perm_model->getValue(S[i]), 1.e-9);
+
+ // Compare the derivative with numerical one.
+ ASSERT_NEAR((perm_model->getValue(S[i] + purterbation) -
+ perm_model->getValue(S[i])) /
+ purterbation,
+ perm_model->getdValue(S[i]), 1.e-6);
}
}
--
GitLab