From 2c4ea3769d48f9af8c7a4fb39ef7430c8ed6638a Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Tue, 1 Nov 2016 16:31:19 +0100
Subject: [PATCH] [Unsat] Added Brook-Corey capillary pressure model
---
.../BrookCoreyCapillaryPressure.cpp | 52 ++++++++++
.../BrookCoreyCapillaryPressure.h | 97 +++++++++++++++++++
MathLib/MathTools.h | 10 ++
3 files changed, 159 insertions(+)
create mode 100644 MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressure.cpp
create mode 100644 MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressure.h
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressure.cpp b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressure.cpp
new file mode 100644
index 00000000000..d76a57c5c5b
--- /dev/null
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressure.cpp
@@ -0,0 +1,52 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file BrookCoreyCapillaryPressure.cpp
+ *
+ * Created on November 1, 2016, 9:45 AM
+ */
+
+#include "BrookCoreyCapillaryPressure.h"
+
+#include <cmath>
+
+#include "MathLib/MathTools.h"
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+double BrookCoreyCapillaryPressure::getCapillaryPressure(
+ const double saturation) const
+{
+ const double S = MathLib::limitValueInInterval(
+ saturation, _Sr + _perturbation, _Smax - _perturbation);
+ const double Se = (S - _Sr) / (_Smax - _Sr);
+ const double pc = _pb * std::pow(Se, -1.0 / _mm);
+ return MathLib::limitValueInInterval(pc, _perturbation, _Pc_max);
+}
+
+double BrookCoreyCapillaryPressure::getSturation(
+ const double capillary_pressure) const
+{
+ const double pc = (capillary_pressure < 0.0) ? 0.0 : capillary_pressure;
+ const double Se = std::pow(pc / _pb, -_mm);
+ const double S = Se * (_Smax - _Sr) + _Sr;
+ return MathLib::limitValueInInterval(S, _Sr + _perturbation,
+ _Smax - _perturbation);
+}
+
+double BrookCoreyCapillaryPressure::getdPcdS(const double saturation) const
+{
+ const double S = MathLib::limitValueInInterval(
+ saturation, _Sr + _perturbation, _Smax - _perturbation);
+ const double val = std::pow(((S - _Sr) / (_Smax - _Sr)), -1.0 / _mm);
+ return (_pb * val) / (_mm * (_Sr - S));
+}
+
+} // end namespace
+} // end namespace
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressure.h b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressure.h
new file mode 100644
index 00000000000..5aabe718117
--- /dev/null
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressure.h
@@ -0,0 +1,97 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file BrookCoreyCapillaryPressure.h
+ *
+ * Created on November 1, 2016, 9:45 AM
+ */
+
+#ifndef OGS_BROOK_COREY_CAPILLARY_PRESSURE_H
+#define OGS_BROOK_COREY_CAPILLARY_PRESSURE_H
+
+#include <array>
+#include <limits> // std::numeric_limits
+
+#include "CapillaryPressureSaturation.h"
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+/**
+ * \brief Brook-Corey capillary pressure saturation model
+ *
+ * \f[p_c=p_b S_e^{-1/m}\f]
+ * with
+ * \f[S_e=\dfrac{S-S_r}{S_{\mbox{max}}-S_r}\f]
+ * where
+ * \f{eqnarray*}{
+ * &p_b& \mbox{ entry pressure,}\\
+ * &S_r& \mbox{ residual saturation,}\\
+ * &S_{\mbox{max}}& \mbox{ maximum saturation,}\\
+ * &m(>=1) & \mbox{ exponent.}\\
+ * \f}
+ *
+ * Note:
+ * \f[m=1/(1-n)\f]
+ */
+class BrookCoreyCapillaryPressure final : public CapillaryPressureSaturation
+{
+public:
+ /** \param parameters An array contains the five parameters:
+ * [0] $f p_b $f
+ * [1] $f S_r $f
+ * [2] $f S_{\mbox{max}} $f
+ * [3] $f m $f
+ * [4] $f P_c^{\mbox{max}}$f
+ */
+ BrookCoreyCapillaryPressure(std::array<double, 5> const& parameters)
+ : _pb(parameters[0]),
+ _Sr(parameters[1]),
+ _Smax(parameters[2]),
+ _mm(parameters[3]),
+ _Pc_max(parameters[4])
+ {
+ }
+
+ BrookCoreyCapillaryPressure(const BrookCoreyCapillaryPressure& orig)
+ : _pb(orig._pb),
+ _Sr(orig._Sr),
+ _Smax(orig._Smax),
+ _mm(orig._mm),
+ _Pc_max(orig._Pc_max)
+ {
+ }
+
+ /// Get model name.
+ std::string getName() const override
+ {
+ return "Brook-Corey capillary pressure saturation model.";
+ }
+
+ /// Get capillary pressure.
+ double getCapillaryPressure(const double saturation) const override;
+
+ /// Get capillary pressure.
+ double getSturation(const double capillary_pressure) const override;
+
+ /// Get the derivative of the capillary pressure with respect to saturation
+ double getdPcdS(const double saturation) const override;
+
+private:
+ const double _pb; ///< Entry pressure.
+ const double _Sr; ///< Residual saturation.
+ const double _Smax; ///< Maximum saturation.
+ const double _mm; ///< Exponent (<=1.0), n=1/(1-mm).
+ const double _Pc_max; ///< Maximum capillaray pressure
+
+ const double _perturbation = std::numeric_limits<double>::epsilon();
+};
+
+} // end namespace
+} // end namespace
+#endif /* OGS_BROOK_COREY_CAPILLARY_PRESSURE_H */
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index e578be54451..0ee2bd42ae9 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -120,6 +120,16 @@ inline double to_radians(double degrees) {
return degrees*boost::math::constants::pi<double>()/180.;
}
+template<typename Type> Type limitValueInInterval(const Type variable,
+ const Type lower_bound,
+ const Type upper_bound)
+{
+ if (variable < lower_bound)
+ return lower_bound;
+ if (variable > upper_bound)
+ return upper_bound;
+ return variable;
+}
} // namespace
#endif /* MATHTOOLS_H_ */
--
GitLab