[MPL] Added a class of CapillaryPressureRegularizedVanGenuchten

Documentation/bibliography.bib

@@ -284,3 +284,27 @@
   doi = "10.1007/978-3-319-68225-9",
   isbn = "978-3-319-68224-2"
 }
+
+@article{marchand2013fully,
+  title={{F}ully coupled generalized hybrid-mixed finite element approximation
+         of two-phase two-component flow in porous media. {Part I}:
+  formulation and properties of the mathematical model},
+  author={Marchand, E and M{\"u}ller, T and Knabner, P},
+  journal={Computational Geosciences},
+  volume={17},
+  number={2},
+  pages={431--442},
+  year={2013},
+  doi = "10.1007/s10596-013-9341-7",
+  publisher={Springer}
+}
+@article{huang2015extending,
+  title={Extending the persistent primary variable algorithm to simulate non-isothermal two-phase two-component flow with phase change phenomena},
+  author={Huang, Yonghui and Kolditz, Olaf and Shao, Haibing},
+  journal={Geothermal Energy},
+  volume={3},
+  number={1},
+  pages={13},
+  year={2015},
+  publisher={Springer}
+}

MaterialLib/MPL/Properties/CapillaryPressureSaturation/CapillaryPressureRegularizedVanGenuchten.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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 April 20, 2020, 9:30 AM
+ */
+
+#include "CapillaryPressureRegularizedVanGenuchten.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "BaseLib/Error.h"
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/VariableType.h"
+#include "MaterialLib/MPL/Utils/CheckVanGenuchtenExponentRange.h"
+
+namespace MaterialPropertyLib
+{
+void checkSaturationRange(const double Sl)
+{
+    if (Sl < 0 || Sl > 1)
+    {
+        OGS_FATAL("The saturation of {:e} is out of its range of [0, 1]", Sl);
+    }
+}
+
+CapillaryPressureRegularizedVanGenuchten::
+    CapillaryPressureRegularizedVanGenuchten(
+        double const residual_liquid_saturation,
+        double const maximum_liquid_saturation,
+        double const exponent,
+        double const p_b)
+    : Sg_r_(1.0 - maximum_liquid_saturation),
+      Sg_max_(1.0 - residual_liquid_saturation),
+      m_(exponent),
+      p_b_(p_b),
+      PcBarvGSg_Sg_max_(getPcBarvGSg(Sg_max_)),
+      dPcdSvGBarSg_max_(getdPcdSvGBar(Sg_max_))
+{
+    checkVanGenuchtenExponentRange(m_);
+}
+
+PropertyDataType CapillaryPressureRegularizedVanGenuchten::value(
+    VariableArray const& variable_array,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+    double const /*dt*/) const
+{
+    const double Sl = std::get<double>(
+        variable_array[static_cast<int>(Variable::liquid_saturation)]);
+
+    checkSaturationRange(Sl);
+
+    double const Sg = 1 - Sl;
+    if (!(Sg < Sg_r_ || Sg > Sg_max_))
+    {
+        return getPcBarvGSg(Sg);
+    }
+    if (Sg < Sg_r_)
+    {
+        return 0.0;
+    }
+
+    return PcBarvGSg_Sg_max_ + dPcdSvGBarSg_max_ * (Sg - Sg_max_);
+}
+
+PropertyDataType CapillaryPressureRegularizedVanGenuchten::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) &&
+        "CapillaryPressureRegularizedVanGenuchten::dValue is implemented for "
+        "derivatives with respect to liquid saturation only.");
+
+    const double Sl = std::get<double>(
+        variable_array[static_cast<int>(Variable::liquid_saturation)]);
+
+    checkSaturationRange(Sl);
+
+    double const Sg = 1 - Sl;
+    if (!(Sg < Sg_r_ || Sg > Sg_max_))
+    {
+        return -getdPcdSvGBar(Sg);
+    }
+    if (Sg < Sg_r_)
+    {
+        return 0.0;
+    }
+
+    return -dPcdSvGBarSg_max_;
+}
+
+double CapillaryPressureRegularizedVanGenuchten::getPcBarvGSg(
+    double const Sg) const
+{
+    double const S_bar = getSBar(Sg);
+    return getPcvGSg(S_bar) - getPcvGSg(Sg_r_ + (Sg_max_ - Sg_r_) * xi_ / 2);
+}
+
+double CapillaryPressureRegularizedVanGenuchten::getSBar(double const Sg) const
+{
+    return Sg_r_ + (1 - xi_) * (Sg - Sg_r_) + 0.5 * xi_ * (Sg_max_ - Sg_r_);
+}
+
+double CapillaryPressureRegularizedVanGenuchten::getPcvGSg(
+    double const Sg) const
+{
+    double const Se = (Sg_max_ - Sg) / (Sg_max_ - Sg_r_);
+    return p_b_ * std::pow(std::pow(Se, (-1.0 / m_)) - 1.0, 1.0 - m_);
+}
+
+double CapillaryPressureRegularizedVanGenuchten::getdPcdSvGBar(
+    double const Sg) const
+{
+    double S_bar = getSBar(Sg);
+    return getdPcdSvG(S_bar) * (1 - xi_);
+}
+
+double CapillaryPressureRegularizedVanGenuchten::getdPcdSvG(
+    const double Sg) const
+{
+    double const n = 1 / (1 - m_);
+    double const Se = (Sg_max_ - Sg) / (Sg_max_ - Sg_r_);
+    auto const temp = std::pow(Se, (-1 / m_));
+    return p_b_ * (1 / (m_ * n)) * (1 / (Sg_max_ - Sg_r_)) *
+           std::pow(temp - 1, (1 / n) - 1) * temp / Se;
+}
+
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CapillaryPressureSaturation/CapillaryPressureRegularizedVanGenuchten.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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 April 20, 2020, 9:30 AM
+ */
+
+#pragma once
+
+#include "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+class Medium;
+
+/**
+ *  \brief This class handles the computation of the capillary pressure,
+ *  \f$ p_c(S_l) \f$, with the
+ * capillary pressure regularization.
+ *
+ *   For the regularized van Genuchten model, please refer to
+ *  \cite huang2015extending.
+ *
+ *   For the van Genuchten capillary pressure mode,
+ *   \sa MaterialPropertyLib::CapillaryPressureVanGenuchten
+ */
+class CapillaryPressureRegularizedVanGenuchten final : public Property
+{
+public:
+    CapillaryPressureRegularizedVanGenuchten(
+        double const residual_liquid_saturation,
+        double const maximum_liquid_saturation,
+        double const exponent,
+        double const p_b);
+
+    void checkScale() const override
+    {
+        if (!std::holds_alternative<Medium*>(scale_))
+        {
+            OGS_FATAL(
+                "The property 'CapillaryPressureRegularizedVanGenuchten' is "
+                "implemented on the 'media' scale only.");
+        }
+    }
+
+    /// \return \f$ p_c(S_l) \f$.
+    PropertyDataType value(VariableArray const& variable_array,
+                           ParameterLib::SpatialPosition const& pos,
+                           double const t,
+                           double const dt) const override;
+
+    /// \return \f$ \frac{\partial p_c(S_l)}{\partial  S_l} \f$
+    PropertyDataType dValue(VariableArray const& variable_array,
+                            Variable const variable,
+                            ParameterLib::SpatialPosition const& pos,
+                            double const t,
+                            double const dt) const override;
+
+private:
+    double const Sg_r_;    ///< Residual saturation of gas phase
+    double const Sg_max_;  ///< Maximum saturation of gas phase
+    double const m_;       ///< Exponent.
+    /// Capillary pressure scaling factor. Sometimes, it is called apparent gas
+    /// entry pressure.
+    double const p_b_;
+    /// parameter in regularized van Genuchten model
+    static constexpr double xi_ = 1e-5;
+
+    double const PcBarvGSg_Sg_max_;
+    double const dPcdSvGBarSg_max_;
+
+    /// Gets regularized capillary pressure via the saturation of the
+    /// non-wetting phase.
+    double getPcBarvGSg(double const Sg) const;
+    /// Gets regularized saturation via the saturation of the non-wetting
+    /// phase.
+    double getSBar(double const Sg) const;
+    /// Gets capillary pressure via the non-wetting phase saturation with the
+    /// original van Genuchten model.
+    double getPcvGSg(double const Sg) const;
+    /// Gets \f$\frac{\partial p_c}{\partial {\bar S}_g }\f$.
+    double getdPcdSvGBar(double const Sg) const;
+    /// Gets \f$\frac{\partial p_c}{\partial {S}_g }\f$.
+    double getdPcdSvG(double const Sg) const;
+};
+
+}  // namespace MaterialPropertyLib