Merge pull request #2656 from Scinopode/mpl_brooks_corey_relative_permeability

New MPL property: Brooks&Corey relative permeability

Documentation/ProjectFile/properties/property/RelPermBrooksCorey/c_RelPermBrooksCorey.md

+The relative permeability function proposed by Brooks&Corey (1964) \cite BrooksCorey1964.

Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_k_rel_min_gas.md

+The minimal relative permeability of the gas (or non-wetting) phase.

Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_k_rel_min_liquid.md

+The minimal relative permeability of the liquid (or wetting) phase.

Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_lambda.md

+Exponent of the Brooks&Corey relative permeability function.

Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_residual_gas_saturation.md

+The smallest degree of saturation of the gas (or non-wetting) phase in the medium.

Documentation/ProjectFile/properties/property/RelPermBrooksCorey/t_residual_liquid_saturation.md

+The smallest degree of saturation of the liquid (or wetting) phase in the medium.

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());

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

MaterialLib/MPL/Properties/CreateRelPermBrooksCorey.cpp

+/**
+ * \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

MaterialLib/MPL/Properties/CreateRelPermBrooksCorey.h

+/**
+ * \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

MaterialLib/MPL/Properties/Properties.h

@@ -17,4 +17,5 @@
 #include "IdealGasLaw.h"
 #include "LinearProperty.h"
 #include "ParameterProperty.h"
+#include "RelPermBrooksCorey.h"
 #include "SaturationBrooksCorey.h"

MaterialLib/MPL/Properties/RelPermBrooksCorey.cpp

+/**
+ * \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

MaterialLib/MPL/Properties/RelPermBrooksCorey.h

+/**
+ * \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

Tests/MaterialLib/TestMPLRelPermBrooksCorey.cpp

+/**
+ * \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