From 3be867a18c380f756e16f2c22660a794d2d3cc6b Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Mon, 23 Mar 2020 15:37:23 +0100
Subject: [PATCH] [MPL] Added a class of
CapillaryPressureRegularizedVanGenuchten
---
Documentation/bibliography.bib | 24 ++++
...pillaryPressureRegularizedVanGenuchten.cpp | 136 ++++++++++++++++++
...CapillaryPressureRegularizedVanGenuchten.h | 91 ++++++++++++
3 files changed, 251 insertions(+)
create mode 100644 MaterialLib/MPL/Properties/CapillaryPressureSaturation/CapillaryPressureRegularizedVanGenuchten.cpp
create mode 100644 MaterialLib/MPL/Properties/CapillaryPressureSaturation/CapillaryPressureRegularizedVanGenuchten.h
diff --git a/Documentation/bibliography.bib b/Documentation/bibliography.bib
index b55ba6ca395..0d1386d341c 100644
--- a/Documentation/bibliography.bib
+++ b/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}
+}
diff --git a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CapillaryPressureRegularizedVanGenuchten.cpp b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CapillaryPressureRegularizedVanGenuchten.cpp
new file mode 100644
index 00000000000..c0c042bb450
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CapillaryPressureRegularizedVanGenuchten.cpp
@@ -0,0 +1,136 @@
+/**
+ * \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
diff --git a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CapillaryPressureRegularizedVanGenuchten.h b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CapillaryPressureRegularizedVanGenuchten.h
new file mode 100644
index 00000000000..36ed1a0f6b0
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CapillaryPressureRegularizedVanGenuchten.h
@@ -0,0 +1,91 @@
+/**
+ * \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
--
GitLab