[MPL] Added a new class of RelPermUdellNonwettingPhase, which is subtracted from RelPermUdell

MaterialLib/MPL/Properties/RelativePermeability/CreateRelPermUdellNonwettingPhase.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * Created on January 26, 2021, 11:11 AM
+ */
+
+#include "BaseLib/ConfigTree.h"
+#include "RelPermUdellNonwettingPhase.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<RelPermUdellNonwettingPhase> createRelPermUdellNonwettingPhase(
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{properties__property__type}
+    config.checkConfigParameter("type",
+                                "RelativePermeabilityUdellNonwettingPhase");
+
+    // Second access for storage.
+    //! \ogs_file_param{properties__property__name}
+    auto property_name = config.peekConfigParameter<std::string>("name");
+
+    DBUG("Create RelPermUdellNonwettingPhase medium property {:s}.",
+         property_name);
+
+    auto const residual_liquid_saturation =
+        //! \ogs_file_param{properties__property__RelativePermeabilityUdellNonwettingPhase__residual_liquid_saturation}
+        config.getConfigParameter<double>("residual_liquid_saturation");
+    auto const residual_gas_saturation =
+        //! \ogs_file_param{properties__property__RelativePermeabilityUdellNonwettingPhase__residual_gas_saturation}
+        config.getConfigParameter<double>("residual_gas_saturation");
+    auto const min_relative_permeability =
+        //! \ogs_file_param{properties__property__RelativePermeabilityUdellNonwettingPhase__min_relative_permeability}
+        config.getConfigParameter<double>("min_relative_permeability");
+
+    if (min_relative_permeability < 0)
+    {
+        OGS_FATAL("Minimal relative permeability must be non-negative.");
+    }
+
+    return std::make_unique<RelPermUdellNonwettingPhase>(
+        std::move(property_name), residual_liquid_saturation,
+        residual_gas_saturation, min_relative_permeability);
+}
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/RelativePermeability/CreateRelPermUdellNonwettingPhase.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * Created on January 26, 2021, 11:11 AM
+ */
+
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class RelPermUdellNonwettingPhase;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<RelPermUdellNonwettingPhase> createRelPermUdellNonwettingPhase(
+    BaseLib::ConfigTree const& config);
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/RelativePermeability/RelPermUdellNonwettingPhase.cpp

+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date   01.12.2020
+ *
+ * \copyright
+ * Copyright (c) 2012-2021, 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 "RelPermUdellNonwettingPhase.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "MaterialLib/MPL/Medium.h"
+
+namespace MaterialPropertyLib
+{
+RelPermUdellNonwettingPhase::RelPermUdellNonwettingPhase(
+    std::string name,
+    const double residual_liquid_saturation,
+    const double residual_gas_saturation,
+    const double min_relative_permeability)
+    : residual_liquid_saturation_(residual_liquid_saturation),
+      residual_gas_saturation_(residual_gas_saturation),
+      min_relative_permeability_(min_relative_permeability)
+{
+    name_ = std::move(name);
+}
+
+PropertyDataType RelPermUdellNonwettingPhase::value(
+    VariableArray const& variable_array,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+    double const /*dt*/) const
+{
+    const double S_L = std::get<double>(
+        variable_array[static_cast<int>(Variable::liquid_saturation)]);
+
+    if (std::isnan(S_L))
+    {
+        OGS_FATAL(
+            "In RelPermUdellNonwettingPhase::value, the liquid saturation is "
+            "NaN.");
+    }
+
+    auto const S_L_res = residual_liquid_saturation_;
+    auto const S_L_max = 1. - residual_gas_saturation_;
+
+    auto const S_e = (S_L - S_L_res) / (S_L_max - S_L_res);
+
+    if (S_e >= 1.0)
+    {
+        // fully saturated medium
+        return min_relative_permeability_;
+    }
+    if (S_e <= 0.0)
+    {
+        // dry medium
+        return 1.0;
+    }
+
+    auto const S_e_g = (1. - S_e);
+
+    return std::max(min_relative_permeability_, S_e_g * S_e_g * S_e_g);
+}
+PropertyDataType RelPermUdellNonwettingPhase::dValue(
+    VariableArray const& variable_array, Variable const primary_variable,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+    double const /*dt*/) const
+{
+    (void)primary_variable;
+    assert((primary_variable == Variable::liquid_saturation) &&
+           "RelPermUdellNonwettingPhase::dValue is implemented for "
+           " derivatives with respect to liquid saturation only.");
+
+    const double S_L = std::get<double>(
+        variable_array[static_cast<int>(Variable::liquid_saturation)]);
+
+    auto const S_L_res = residual_liquid_saturation_;
+    auto const S_L_max = 1. - residual_gas_saturation_;
+    auto const S_e = (S_L - S_L_res) / (S_L_max - S_L_res);
+
+    if ((S_e < 0.) || (S_e > 1.))
+    {
+        return 0.;
+    }
+
+    auto const dS_e_dS_L = 1. / (S_L_max - S_L_res);
+
+    auto const dk_rel_GR_dS_e = -3. * (1. - S_e) * (1. - S_e);
+    return dk_rel_GR_dS_e * dS_e_dS_L;
+}
+
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/RelativePermeability/RelPermUdellNonwettingPhase.h

+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date   01.12.2020
+ *
+ * \copyright
+ * Copyright (c) 2012-2021, 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 "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+class Medium;
+class Phase;
+class Component;
+/**
+ * A simple relative permeability function proposed by
+ * Kent S Udell \cite udell1985heat.
+ *
+ *  Definition:
+ *          \f[ k_{\text{rel}}^{\alpha}
+ * =\left(S^{\text{eff}}_{\alpha}\right)^{3}\f] where
+ *          - \f$k_{\text{rel}}^{\alpha}\f$ is relative permeability of phase
+ * \f$\alpha\f$
+ *          - \f$S^{\text{eff}}_{\alpha}\f$ is the effective saturation of
+ * phase \f$\alpha\f$
+ *
+ * This class handles the non-wetting (gas) phase portion of this relative
+ * permeability property, i,e with \f$\alpha=g\f$ for the  relative permeability
+ * function.
+ *
+ * \details This property must be a medium property, it
+ * computes the permeability reduction due to saturation as function of
+ * capillary pressure.
+ */
+class RelPermUdellNonwettingPhase final : public Property
+{
+private:
+    const double residual_liquid_saturation_;
+    const double residual_gas_saturation_;
+    const double min_relative_permeability_;
+
+public:
+    RelPermUdellNonwettingPhase(std::string name,
+                                const double residual_liquid_saturation,
+                                const double residual_gas_saturation,
+                                const double min_relative_permeability);
+
+    void checkScale() const override
+    {
+        if (!std::holds_alternative<Medium*>(scale_))
+        {
+            OGS_FATAL(
+                "The property 'RelativePermeabilityUdellNonwettingPhase' is "
+                "implemented on the 'media' scale only.");
+        }
+    }
+
+    PropertyDataType value(VariableArray const& variable_array,
+                           ParameterLib::SpatialPosition const& pos,
+                           double const t, double const dt) const override;
+    PropertyDataType dValue(VariableArray const& variable_array,
+                            Variable const primary_variable,
+                            ParameterLib::SpatialPosition const& pos,
+                            double const t, double const dt) const override;
+};
+
+}  // namespace MaterialPropertyLib