diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseBrookCoreyOilGas.cpp b/MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseBrookCoreyOilGas.cpp
index 5a6d09da7976996570438a115f7105538b8a74e4..f81ba6d8b7b6a464f152141fb9f77c0e69762cbd 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 32a91922a796eba81e641bf057ba042f3f837a19..9b34f2e54f9c36e6139d154d57f1ce09a377e0d2 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 435c52dc3872e47de4db2a03ceef3a5f953e9775..f3d397ad26abfde45c492d07ebd29c2c63f59df9 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 f24ecd4d22c73687d64ba8ac71b012fca6a30a87..a039616fdb6ce4d127769b0405462373debd0c1f 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 0b62d8ac50ba0e1e20cc6f4bb4616b5901bf8f36..da69bb2a25be16b6b1ce68af25e3ae135342b6e9 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 8232932a5ef3a91af62e2961860f709e6f906228..b89b4c7524eb1b61fb97e39e65c5106d12a0113e 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 938bf57d2505b8e040b25ee297bddc473b595a57..cb8cebc589a2cc4156acb64c15824c5c61d83972 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 28b095d9982e1670e7bd985352e0807f768df88a..93283367763b4b5062f4bd9e5d5dbfa1dc6a1aaa 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 fdcb70908ac17d658756c9425a3c01c44312f99c..e1429bce38b1d9fbd72f3d1535a9b36f9771b9a0 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 8fd580f378e0bfa0c03dd284f5aa630de67e7b65..c69b0850b8d6d7521918b81d312fc84e5c9c7b4a 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 ae340266fab13e01d2aa5928b0cfcbfba6b15a56..7612ed4d8e010f28b8d1f32e67599413cdd2969a 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);
}
}