From 4d122ec94acdb8aeaea36673a286eca548f60d82 Mon Sep 17 00:00:00 2001
From: Norbert Grunwald <Norbert.Grunwald@ufz.de>
Date: Tue, 20 Aug 2019 17:54:08 +0200
Subject: [PATCH] introduce Brooks&Corey relative permeability to MPL
---
MaterialLib/MPL/CreateProperty.cpp | 5 +
.../MPL/Properties/RelPermBrooksCorey.cpp | 108 ++++++++++++++++++
.../MPL/Properties/RelPermBrooksCorey.h | 107 +++++++++++++++++
3 files changed, 220 insertions(+)
create mode 100644 MaterialLib/MPL/Properties/RelPermBrooksCorey.cpp
create mode 100644 MaterialLib/MPL/Properties/RelPermBrooksCorey.h
diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp
index 02646eddf90..0c165ece888 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/RelPermBrooksCorey.cpp b/MaterialLib/MPL/Properties/RelPermBrooksCorey.cpp
new file mode 100644
index 00000000000..eb760f3a5b7
--- /dev/null
+++ b/MaterialLib/MPL/Properties/RelPermBrooksCorey.cpp
@@ -0,0 +1,108 @@
+/**
+ * \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
+{
+ 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);
+ 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 00000000000..010b5985e00
--- /dev/null
+++ b/MaterialLib/MPL/Properties/RelPermBrooksCorey.h
@@ -0,0 +1,107 @@
+/**
+ * \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 meterial 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;
+};
+
+inline 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");
+
+ return std::make_unique<RelPermBrooksCorey>(
+ residual_liquid_saturation,
+ residual_gas_saturation,
+ min_relative_permeability_liquid,
+ min_relative_permeability_gas,
+ exponent);
+}
+
+} // namespace MaterialPropertyLib
--
GitLab