diff --git a/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/c_RelPermBrooksCorey.md b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/c_RelPermBrooksCorey.md
new file mode 100644
index 0000000000000000000000000000000000000000..5637495c26abe04ac29b497ddced0659d4dba443
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/c_RelPermBrooksCorey.md
@@ -0,0 +1 @@
+The relative permeability function proposed by Brooks&Corey (1964) \cite BrooksCorey1964.
diff --git a/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_k_rel_min_gas.md b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_k_rel_min_gas.md
new file mode 100644
index 0000000000000000000000000000000000000000..47e4b2a9521528a4d09b572ba23e376238ff6763
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_k_rel_min_gas.md
@@ -0,0 +1 @@
+The minimal relative permeability of the gas (or non-wetting) phase.
diff --git a/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_k_rel_min_liquid.md b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_k_rel_min_liquid.md
new file mode 100644
index 0000000000000000000000000000000000000000..9cb5411f5500a5b03715a932065b976795833924
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_k_rel_min_liquid.md
@@ -0,0 +1 @@
+The minimal relative permeability of the liquid (or wetting) phase.
diff --git a/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_lambda.md b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_lambda.md
new file mode 100644
index 0000000000000000000000000000000000000000..5015ac5ec5313f4071a60cceb056c6d8b7a6a9d2
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_lambda.md
@@ -0,0 +1 @@
+Exponent of the Brooks&Corey relative permeability function.
diff --git a/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_residual_gas_saturation.md b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_residual_gas_saturation.md
new file mode 100644
index 0000000000000000000000000000000000000000..d8631eb27f4627345443d3cb426330a84c325db1
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_residual_gas_saturation.md
@@ -0,0 +1 @@
+The smallest degree of saturation of the gas (or non-wetting) phase in the medium.
diff --git a/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_residual_liquid_saturation.md b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_residual_liquid_saturation.md
new file mode 100644
index 0000000000000000000000000000000000000000..95ce0ed4b319db3ecc0cfc87b42c611dcd572210
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_residual_liquid_saturation.md
@@ -0,0 +1 @@
+The smallest degree of saturation of the liquid (or wetting) phase in the medium.
diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp
index 02646eddf901dc4a1e81856cb37c3287cdee0f23..0c165ece8880ee4aa122b6e723e378813d8cb5de 100644
--- a/MaterialLib/MPL/CreateProperty.cpp
+++ b/MaterialLib/MPL/CreateProperty.cpp
@@ -65,6 +65,11 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
return createSaturationBrooksCorey(config);
}
+ if (boost::iequals(property_type, "RelPermBrooksCorey"))
+ {
+ return createRelPermBrooksCorey(config);
+ }
+
// If none of the above property types are found, OGS throws an error.
OGS_FATAL("The specified component property type '%s' was not recognized",
property_type.c_str());
diff --git a/MaterialLib/MPL/Properties/CreateProperties.h b/MaterialLib/MPL/Properties/CreateProperties.h
index 646a6e2213ab21729599674913def917a1b292cb..386d7be0313304e13b9f9f3e53e82fd6dab41e97 100644
--- a/MaterialLib/MPL/Properties/CreateProperties.h
+++ b/MaterialLib/MPL/Properties/CreateProperties.h
@@ -17,4 +17,5 @@
#include "CreateIdealGasLaw.h"
#include "CreateLinearProperty.h"
#include "CreateParameterProperty.h"
+#include "CreateRelPermBrooksCorey.h"
#include "CreateSaturationBrooksCorey.h"
\ No newline at end of file
diff --git a/MaterialLib/MPL/Properties/CreateRelPermBrooksCorey.cpp b/MaterialLib/MPL/Properties/CreateRelPermBrooksCorey.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b9b0f89539e2e6e0946e223883e7dc65eb60c98a
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreateRelPermBrooksCorey.cpp
@@ -0,0 +1,51 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "BaseLib/ConfigTree.h"
+#include "RelPermBrooksCorey.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<RelPermBrooksCorey> createRelPermBrooksCorey(
+ BaseLib::ConfigTree const& config)
+{
+ // check is reading the parameter, not peeking it...
+ //! \ogs_file_param{prj__media__medium__properties__property__RelPermBrooksCorey}
+ config.checkConfigParameter("type", "RelPermBrooksCorey");
+ DBUG("Create RelPermBrooksCorey medium property");
+
+ auto const residual_liquid_saturation =
+ //! \ogs_file_param{prj__media__medium__properties__property__RelPermBrooksCorey__residual_liquid_saturation}
+ config.getConfigParameter<double>("residual_liquid_saturation");
+ auto const residual_gas_saturation =
+ //! \ogs_file_param{prj__media__medium__properties__property__RelPermBrooksCorey__residual_gas_saturation}
+ config.getConfigParameter<double>("residual_gas_saturation");
+ auto const min_relative_permeability_liquid =
+ //! \ogs_file_param{prj__media__medium__properties__property__RelPermBrooksCorey__k_rel_min_liquid}
+ config.getConfigParameter<double>("min_relative_permeability_liquid");
+ auto const min_relative_permeability_gas =
+ //! \ogs_file_param{prj__media__medium__properties__property__RelPermBrooksCorey__k_rel_min_gas}
+ config.getConfigParameter<double>("min_relative_permeability_gas");
+ auto const exponent =
+ //! \ogs_file_param{prj__media__medium__properties__property__RelPermBrooksCorey__lambda}
+ config.getConfigParameter<double>("lambda");
+ if (exponent <= 0.)
+ {
+ OGS_FATAL("Exponent 'lambda' must be positive.");
+ }
+
+ return std::make_unique<RelPermBrooksCorey>(
+ residual_liquid_saturation,
+ residual_gas_saturation,
+ min_relative_permeability_liquid,
+ min_relative_permeability_gas,
+ exponent);
+}
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateRelPermBrooksCorey.h b/MaterialLib/MPL/Properties/CreateRelPermBrooksCorey.h
new file mode 100644
index 0000000000000000000000000000000000000000..691059e25b16e409b8cea4d0ef41fa41e52be81e
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreateRelPermBrooksCorey.h
@@ -0,0 +1,27 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class RelPermBrooksCorey;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<RelPermBrooksCorey> createRelPermBrooksCorey(
+ BaseLib::ConfigTree const& config);
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/Properties.h b/MaterialLib/MPL/Properties/Properties.h
index bd4b40595de10470c7f9ee174806a5f92839190d..24797c117090acba70042a4219a1865f09c7cc51 100644
--- a/MaterialLib/MPL/Properties/Properties.h
+++ b/MaterialLib/MPL/Properties/Properties.h
@@ -17,4 +17,5 @@
#include "IdealGasLaw.h"
#include "LinearProperty.h"
#include "ParameterProperty.h"
+#include "RelPermBrooksCorey.h"
#include "SaturationBrooksCorey.h"
diff --git a/MaterialLib/MPL/Properties/RelPermBrooksCorey.cpp b/MaterialLib/MPL/Properties/RelPermBrooksCorey.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6ff47331871b9777a3ca78f7023e4973622fef5b
--- /dev/null
+++ b/MaterialLib/MPL/Properties/RelPermBrooksCorey.cpp
@@ -0,0 +1,111 @@
+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date 02.07.2018
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "MaterialLib/MPL/Properties/RelPermBrooksCorey.h"
+#include "MaterialLib/MPL/Medium.h"
+
+#include <algorithm>
+#include <cmath>
+
+namespace MaterialPropertyLib
+{
+RelPermBrooksCorey::RelPermBrooksCorey(
+ const double residual_liquid_saturation,
+ const double residual_gas_saturation,
+ const double min_relative_permeability_liquid,
+ const double min_relative_permeability_gas,
+ const double exponent)
+ : _residual_liquid_saturation(residual_liquid_saturation),
+ _residual_gas_saturation(residual_gas_saturation),
+ _min_relative_permeability_liquid(min_relative_permeability_liquid),
+ _min_relative_permeability_gas(min_relative_permeability_gas),
+ _exponent(exponent){};
+
+PropertyDataType RelPermBrooksCorey::value(
+ VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& pos,
+ double const t) const
+{
+ /// here, an extra computation of saturation is forced, guaranteeing a
+ /// correct value. In order to speed up the computing time, saturation could
+ /// be insertred into the primary variable array after it is computed in the
+ /// FEM assembly.
+ auto const s_L = _medium->property(PropertyType::saturation)
+ .template value<double>(variable_array, pos, t);
+
+ auto const s_L_res = _residual_liquid_saturation;
+ auto const s_L_max = 1. - _residual_gas_saturation;
+ auto const k_rel_min_LR = _min_relative_permeability_liquid;
+ auto const k_rel_min_GR = _min_relative_permeability_gas;
+
+ auto const lambda = _exponent;
+
+ auto const s_eff = (s_L - s_L_res) / (s_L_max - s_L_res);
+
+ if (s_eff >= 1.0)
+ {
+ // fully saturated medium
+ return Pair{{1.0, k_rel_min_GR}};
+ }
+ if (s_eff <= 0.0)
+ {
+ // dry medium
+ return Pair{{k_rel_min_LR, 1.0}};
+ }
+
+ auto const k_rel_LR = std::pow(s_eff, (2. + 3. * lambda) / lambda);
+ auto const k_rel_GR = (1. - s_eff) * (1. - s_eff) *
+ (1. - std::pow(s_eff, (2. + lambda) / lambda));
+
+ const Pair kRel = {std::max(k_rel_LR, k_rel_min_LR),
+ std::max(k_rel_GR, k_rel_min_GR)};
+ return kRel;
+}
+PropertyDataType RelPermBrooksCorey::dValue(
+ VariableArray const& variable_array, Variable const primary_variable,
+ ParameterLib::SpatialPosition const& pos, double const t) const
+{
+ (void)primary_variable;
+ assert((primary_variable == Variable::liquid_saturation) &&
+ "RelPermBrooksCorey::dValue is implemented for "
+ " derivatives with respect to liquid saturation only.");
+ auto const s_L = _medium->property(PropertyType::saturation)
+ .template value<double>(variable_array, pos, t);
+
+ auto const s_L_res = _residual_liquid_saturation;
+ auto const s_L_max = 1. - _residual_gas_saturation;
+ auto const lambda = _exponent;
+
+ auto const s_eff = (s_L - s_L_res) / (s_L_max - s_L_res);
+ if ((s_eff < 0.) || (s_eff > 1.))
+ return Pair{0., 0.};
+
+ auto const d_se_d_sL = 1. / (s_L_max - s_L_res);
+ auto const dk_rel_LRdse =
+ (3 * lambda + 2.) / lambda * std::pow(s_eff, 2. / lambda + 2.);
+
+ auto const dk_rel_LRdsL = dk_rel_LRdse * d_se_d_sL;
+
+ auto const _2L_L = (2. + lambda) / lambda;
+ auto const dk_rel_GRdse =
+ -2. * (1 - s_eff) * (1. - std::pow(s_eff, _2L_L)) -
+ _2L_L * std::pow(s_eff, _2L_L - 1.) * (1. - s_eff) * (1. - s_eff);
+
+ auto const dk_rel_GRdsL = dk_rel_GRdse * d_se_d_sL;
+ const Pair dkReldsL = {{dk_rel_LRdsL, dk_rel_GRdsL}};
+
+ return dkReldsL;
+}
+
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/RelPermBrooksCorey.h b/MaterialLib/MPL/Properties/RelPermBrooksCorey.h
new file mode 100644
index 0000000000000000000000000000000000000000..a387ac129a5c78a0dbc15c0fc21f807d3d0d57d3
--- /dev/null
+++ b/MaterialLib/MPL/Properties/RelPermBrooksCorey.h
@@ -0,0 +1,76 @@
+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date 27.06.2018
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+#pragma once
+
+#include "BaseLib/ConfigTree.h"
+#include "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+class Medium;
+class Phase;
+class Component;
+/**
+ * \class RelPermBrooksCorey
+ * \brief Relative permeability function proposed by Brooks&Corey
+ * \details This property must be a medium property, it
+ * computes the permeability reduction due to saturation as function
+ * of capillary pressure.
+ */
+class RelPermBrooksCorey final : public Property
+{
+private:
+ Medium* _medium;
+ const double _residual_liquid_saturation;
+ const double _residual_gas_saturation;
+ const double _min_relative_permeability_liquid;
+ const double _min_relative_permeability_gas;
+ const double _exponent;
+
+public:
+ RelPermBrooksCorey(const double /*residual_liquid_saturation*/,
+ const double /*residual_gas_saturation*/,
+ const double /*_min_relative_permeability_liquid*/,
+ const double /*_min_relative_permeability_gas*/,
+ const double /*exponent*/
+ );
+ /// This method assigns a pointer to the material object that is the owner
+ /// of this property
+ void setScale(
+ std::variant<Medium*, Phase*, Component*> scale_pointer) override
+ {
+ if (std::holds_alternative<Medium*>(scale_pointer))
+ {
+ _medium = std::get<Medium*>(scale_pointer);
+ }
+ else
+ {
+ OGS_FATAL(
+ "The property 'RelPermBrooksCorey' is implemented on the "
+ "'media' scale only.");
+ }
+ }
+
+ /// Those methods override the base class implementations and
+ /// actually compute and set the property _values and _dValues.
+ PropertyDataType value(VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& pos,
+ double const t) const override;
+ PropertyDataType dValue(VariableArray const& variable_array,
+ Variable const variable,
+ ParameterLib::SpatialPosition const& pos,
+ double const t) const override;
+};
+
+} // namespace MaterialPropertyLib
diff --git a/Tests/MaterialLib/TestMPLRelPermBrooksCorey.cpp b/Tests/MaterialLib/TestMPLRelPermBrooksCorey.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a20582c2fc3f01ca1ef925cc79e7466cfe1808f2
--- /dev/null
+++ b/Tests/MaterialLib/TestMPLRelPermBrooksCorey.cpp
@@ -0,0 +1,125 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+#include <gtest/gtest.h>
+#include <sstream>
+
+#include "TestMPL.h"
+#include "Tests/TestTools.h"
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Properties/RelPermBrooksCorey.h"
+
+TEST(MaterialPropertyLib, RelPermBrooksCorey)
+{
+ const double ref_lambda = 2.5;
+ const double ref_residual_liquid_saturation = 0.01;
+ const double ref_residual_gas_saturation = 0.01;
+ const double ref_k_rel_L_min = 1.e-9;
+ const double ref_k_rel_G_min = 1.e-7;
+
+ std::stringstream m_beg;
+ std::stringstream m_end;
+
+ m_beg << "<medium>\n";
+ m_beg << "<phases></phases>\n";
+ m_beg << "<properties>\n";
+ m_beg << " <property>\n";
+ m_beg << " <name>saturation</name>\n";
+ m_beg << " <type>Constant</type>\n";
+ // constant saturation value will be inserted here...
+ m_end << " </property>\n";
+ m_end << " <property>\n";
+ m_end << " <name>relative_permeability</name>\n";
+ m_end << " <type>RelPermBrooksCorey</type>\n";
+ m_end << " <residual_liquid_saturation>"
+ << ref_residual_liquid_saturation
+ << "</residual_liquid_saturation>\n";
+ m_end << " <residual_gas_saturation>" << ref_residual_gas_saturation
+ << "</residual_gas_saturation>\n";
+ m_end << " <lambda>" << ref_lambda << "</lambda>\n";
+ m_end << " <min_relative_permeability_liquid>" << ref_k_rel_L_min
+ << "</k_rel_min_liquid>\n";
+ m_end << " <min_relative_permeability_gas>" << ref_k_rel_G_min
+ << "</k_rel_min_gas>\n";
+ m_end << " </property>\n";
+ m_end << "</properties>\n";
+ m_end << "</medium>\n";
+
+ std::array<double, 16> ref_saturation = {
+ -0.01, 0.00, 0.01, 0.02, 0.04, 0.10, 0.20, 0.40,
+ 0.60, 0.80, 0.90, 0.96, 0.98, 0.99, 1.00, 1.01};
+
+ std::array<double, 16> ref_k_rel_L = {
+ 1.0000000000E-09, 1.0000000000E-09, 1.0000000000E-09, 2.7123199126E-08,
+ 1.7636064574E-06, 1.1467333635E-04, 1.9615735463E-03, 3.0156880778E-02,
+ 1.4540473747E-01, 4.4088351375E-01, 6.9346243084E-01, 8.8856788446E-01,
+ 9.6177504114E-01, 1.0000000000E+00, 1.0000000000E+00, 1.0000000000E+00};
+
+ std::array<double, 16> ref_k_rel_G = {
+ 1.0000000000E+00, 1.0000000000E+00, 1.0000000000E+00, 9.7944075784E-01,
+ 9.3794411438E-01, 8.1354659673E-01, 6.1592154540E-01, 2.9343532599E-01,
+ 9.4837870477E-02, 1.2086349932E-02, 1.3426518034E-03, 5.1003060118E-05,
+ 1.9046571241E-06, 1.0000000000E-07, 1.0000000000E-07, 1.0000000000E-07};
+
+ std::array<double, 16> ref_dk_rel_L_ds_L = {
+ 0.0000000000E+00, 0.0000000000E+00, 0.0000000000E+00, 1.0306815668E-05,
+ 2.2339015126E-04, 4.8417630905E-03, 3.9231470925E-02, 2.9383627425E-01,
+ 9.3650508877E-01, 2.1207055092E+00, 2.9608508283E+00, 3.5542715378E+00,
+ 3.7677785117E+00, 3.8775510204E+00, 0.0000000000E+00, 0.0000000000E+00};
+
+ std::array<double, 16> ref_dk_rel_G_ds_L = {
+ 0.0000000000E+00, 0.0000000000E+00, -2.0408163265E+00,
+ -2.0654018726E+00, -2.0807295100E+00, -2.0524729938E+00,
+ -1.8805652791E+00, -1.3132397982E+00, -6.8017950205E-01,
+ -1.8533091175E-01, -4.4178730635E-02, -5.0791425970E-03,
+ -5.7061547092E-04, -0.0000000000E+00, 0.0000000000E+00,
+ 0.0000000000E+00};
+
+ for (size_t idx = 0; idx < ref_saturation.size(); idx++)
+ {
+ std::stringstream m_sat;
+ m_sat << " <value>" << ref_saturation[idx] << "</value>\n";
+ std::stringstream m;
+ m << m_beg.str() << m_sat.str() << m_end.str();
+
+ auto const& medium = createTestMaterial(m.str());
+ MaterialPropertyLib::VariableArray variable_array;
+ ParameterLib::SpatialPosition const pos;
+ double const time = std::numeric_limits<double>::quiet_NaN();
+
+ auto k_rel =
+ medium
+ ->property(
+ MaterialPropertyLib::PropertyType::relative_permeability)
+ .template value<MaterialPropertyLib::Pair>(variable_array, pos,
+ time);
+
+ auto dk_rel_ds_L =
+ medium
+ ->property(
+ MaterialPropertyLib::PropertyType::relative_permeability)
+ .template dValue<MaterialPropertyLib::Pair>(
+ variable_array,
+ MaterialPropertyLib::Variable::liquid_saturation, pos,
+ time);
+
+ auto k_rel_L = k_rel[0];
+ auto k_rel_G = k_rel[1];
+
+ auto dk_rel_L_ds_L = dk_rel_ds_L[0];
+ auto dk_rel_G_ds_L = dk_rel_ds_L[1];
+
+ ASSERT_NEAR(k_rel_L, ref_k_rel_L[idx], 1.0e-10);
+ ASSERT_NEAR(k_rel_G, ref_k_rel_G[idx], 1.0e-10);
+ ASSERT_NEAR(dk_rel_L_ds_L, ref_dk_rel_L_ds_L[idx], 1.0e-10);
+ ASSERT_NEAR(dk_rel_G_ds_L, ref_dk_rel_G_ds_L[idx], 1.0e-10);
+ }
+}
\ No newline at end of file