From a649f629d5ca9211f8dbcdf2e0f50d6a3b41d30f Mon Sep 17 00:00:00 2001
From: Norbert Grunwald <Norbert.Grunwald@ufz.de>
Date: Thu, 15 Aug 2019 08:27:35 +0200
Subject: [PATCH] Brooks-Corey saturation model added to MPL properties
---
MaterialLib/MPL/CreateProperty.cpp | 6 +
MaterialLib/MPL/Properties/Properties.h | 1 +
.../MPL/Properties/SaturationBrooksCorey.cpp | 110 ++++++++++++++++++
.../MPL/Properties/SaturationBrooksCorey.h | 100 ++++++++++++++++
4 files changed, 217 insertions(+)
create mode 100644 MaterialLib/MPL/Properties/SaturationBrooksCorey.cpp
create mode 100644 MaterialLib/MPL/Properties/SaturationBrooksCorey.h
diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp
index 4d6eb09f375..7f754fe95ba 100644
--- a/MaterialLib/MPL/CreateProperty.cpp
+++ b/MaterialLib/MPL/CreateProperty.cpp
@@ -58,6 +58,12 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
{
return createIdealGasLaw(config);
}
+
+ if (boost::iequals(property_type, "SaturationBrooksCorey"))
+ {
+ return createSaturationBrooksCorey(config);
+ }
+
/* TODO Additional properties go here, for example:
if (boost::iequals(property_type, "BilinearTemperaturePressure"))
{
diff --git a/MaterialLib/MPL/Properties/Properties.h b/MaterialLib/MPL/Properties/Properties.h
index f8ca051af51..bd4b40595de 100644
--- a/MaterialLib/MPL/Properties/Properties.h
+++ b/MaterialLib/MPL/Properties/Properties.h
@@ -17,3 +17,4 @@
#include "IdealGasLaw.h"
#include "LinearProperty.h"
#include "ParameterProperty.h"
+#include "SaturationBrooksCorey.h"
diff --git a/MaterialLib/MPL/Properties/SaturationBrooksCorey.cpp b/MaterialLib/MPL/Properties/SaturationBrooksCorey.cpp
new file mode 100644
index 00000000000..9f6d0121bcb
--- /dev/null
+++ b/MaterialLib/MPL/Properties/SaturationBrooksCorey.cpp
@@ -0,0 +1,110 @@
+/**
+ * \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
+ *
+ */
+
+#include "SaturationBrooksCorey.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MathLib/MathTools.h"
+
+namespace MaterialPropertyLib
+{
+SaturationBrooksCorey::SaturationBrooksCorey(
+ const double residual_liquid_saturation,
+ const double residual_gas_saturation,
+ const double exponent,
+ const double entry_pressure)
+ : _residual_liquid_saturation(residual_liquid_saturation),
+ _residual_gas_saturation(residual_gas_saturation),
+ _exponent(exponent),
+ _entry_pressure(entry_pressure){};
+
+PropertyDataType SaturationBrooksCorey::value(
+ VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/) const
+{
+ const double p_cap = std::max(
+ std::numeric_limits<double>::min(),
+ std::get<double>(
+ variable_array[static_cast<int>(Variable::capillary_pressure)]));
+
+ const double s_L_res = _residual_liquid_saturation;
+ const double s_L_max = 1.0 - _residual_gas_saturation;
+ const double lambda = _exponent;
+ const double p_b = _entry_pressure;
+
+ if (p_cap <= p_b)
+ return s_L_max;
+
+ const double s_eff = std::pow(p_b / p_cap, lambda);
+ const double s = s_eff * (s_L_max - s_L_res) + s_L_res;
+
+ return MathLib::limitValueInInterval(s, s_L_res + _minor_offset,
+ s_L_max - _minor_offset);
+}
+
+PropertyDataType SaturationBrooksCorey::dValue(
+ VariableArray const& variable_array, Variable const primary_variable,
+ ParameterLib::SpatialPosition const& pos, double const t) const
+{
+ (void)primary_variable;
+ assert((primary_variable == Variable::capillary_pressure) &&
+ "SaturationBrooksCorey::dValue is implemented for "
+ " derivatives with respect to capillary pressure only.");
+
+ const double p_b = _entry_pressure;
+ const double p_cap = std::max(
+ p_b,
+ std::get<double>(
+ variable_array[static_cast<int>(Variable::capillary_pressure)]));
+
+ auto const s_L = _medium->property(PropertyType::saturation)
+ .template value<double>(variable_array, pos, t);
+
+ const double lambda = _exponent;
+
+ return -lambda / p_cap * s_L;
+}
+
+PropertyDataType SaturationBrooksCorey::d2Value(
+ VariableArray const& variable_array, Variable const primary_variable1,
+ Variable const primary_variable2, ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/) const
+{
+ (void)primary_variable1;
+ (void)primary_variable2;
+ assert((primary_variable1 == Variable::capillary_pressure) &&
+ (primary_variable2 == Variable::capillary_pressure) &&
+ "SaturationBrooksCorey::d2Value is implemented for "
+ " derivatives with respect to capillary pressure only.");
+
+ const double p_b = _entry_pressure;
+ const double p_cap = std::max(
+ p_b,
+ std::get<double>(
+ variable_array[static_cast<int>(Variable::capillary_pressure)]));
+
+ const double s_L_res = _residual_liquid_saturation;
+ const double s_L_max = 1.0 - _residual_gas_saturation;
+
+ const double lambda = _exponent;
+
+ return lambda * (lambda + 1) * std::pow(p_b / p_cap, lambda) /
+ (p_cap * p_cap) * (s_L_max - s_L_res);
+}
+
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/SaturationBrooksCorey.h b/MaterialLib/MPL/Properties/SaturationBrooksCorey.h
new file mode 100644
index 00000000000..2a23e1c77a6
--- /dev/null
+++ b/MaterialLib/MPL/Properties/SaturationBrooksCorey.h
@@ -0,0 +1,100 @@
+/**
+ * \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 <limits>
+#include "BaseLib/ConfigTree.h"
+#include "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+class Medium;
+class Phase;
+class Component;
+/**
+ * \brief A well known soil characteristics function
+ * \details This property must be a medium property, it
+ * computes the saturation of the wetting phase as function
+ * of capillary pressure.
+ */
+class SaturationBrooksCorey final : public Property
+{
+private:
+ Medium* _medium = nullptr;
+ const double _residual_liquid_saturation;
+ const double _residual_gas_saturation;
+ const double _exponent;
+ const double _entry_pressure;
+
+public:
+ SaturationBrooksCorey(const double residual_liquid_saturation,
+ const double residual_gas_saturation,
+ const double exponent,
+ const double entry_pressure);
+
+ void setScale(
+ std::variant<Medium*, Phase*, Component*> scale_pointer) override
+ {
+ if (!std::holds_alternative<Medium*>(scale_pointer))
+ {
+ OGS_FATAL(
+ "The property 'SaturationBrooksCorey' is implemented on the "
+ "'media' scale only.");
+ }
+ _medium = std::get<Medium*>(scale_pointer);
+ }
+
+ /// 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;
+ PropertyDataType d2Value(VariableArray const& variable_array,
+ Variable const variable1,
+ Variable const variable2,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/) const override;
+};
+
+inline std::unique_ptr<SaturationBrooksCorey> createSaturationBrooksCorey(
+ BaseLib::ConfigTree const& config)
+{
+ // check is reading the parameter, not peeking it...
+ //! \ogs_file_param{prj__media__medium__properties__property__SaturationBrooksCorey}
+ // config.checkConfigParameter("type", "SaturationBrooksCorey");
+ DBUG("Create SaturationBrooksCorey medium property");
+
+ auto const residual_liquid_saturation =
+ //! \ogs_file_param{prj__media__medium__properties__property__SaturationBrooksCorey__residual_liquid_saturation}
+ config.getConfigParameter<double>("residual_liquid_saturation");
+ auto const residual_gas_saturation =
+ //! \ogs_file_param{prj__media__medium__properties__property__SaturationBrooksCorey__residual_gas_saturation}
+ config.getConfigParameter<double>("residual_gas_saturation");
+ auto const exponent =
+ //! \ogs_file_param{prj__media__medium__properties__property__SaturationBrooksCorey__lambda}
+ config.getConfigParameter<double>("lambda");
+ auto const entry_pressure =
+ //! \ogs_file_param{prj__media__medium__properties__property__SaturationBrooksCorey__entry_pressure}
+ config.getConfigParameter<double>("entry_pressure");
+
+ return std::make_unique<SaturationBrooksCorey>(residual_liquid_saturation,
+ residual_gas_saturation,
+ exponent, entry_pressure);
+}
+
+} // namespace MaterialPropertyLib
--
GitLab