Merge branch 'peng_robinson' into 'master'

Peng Robinson Equation of State

See merge request ogs/ogs!5097

Documentation/ProjectFile/properties/property/PengRobinson/c_PengRobinson.md

@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::PengRobinson

Documentation/ProjectFile/properties/property/PengRobinson/t_acentric_factor.md

@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::PengRobinson::omega_

Documentation/ProjectFile/properties/property/PengRobinson/t_critical_pressure.md

@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::PengRobinson::pc_

Documentation/ProjectFile/properties/property/PengRobinson/t_critical_temperature.md

@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::PengRobinson::Tc_

MaterialLib/MPL/CreateProperty.cpp

@@ -98,6 +98,11 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
         return createIdealGasLaw(config);
     }
 
+    if (boost::iequals(property_type, "PengRobinson"))
+    {
+        return createPengRobinson(config);
+    }
+
     if (boost::iequals(property_type, "IdealGasLawBinaryMixture"))
     {
         return createIdealGasLawBinaryMixture(config);

MaterialLib/MPL/Properties/CreatePengRobinson.cpp

@@ -0,0 +1,38 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, 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 "CreatePengRobinson.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "ParameterLib/Utils.h"
+#include "PengRobinson.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<Property> createPengRobinson(BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{properties__property__type}
+    config.checkConfigParameter("type", "PengRobinson");
+    DBUG("Create PengRobinson EOS.");
+
+    auto const Tc =
+        //! \ogs_file_param{properties__property__PengRobinson__critical_temperature}
+        config.getConfigParameter<double>("critical_temperature");
+
+    auto const pc =
+        //! \ogs_file_param{properties__property__PengRobinson__critical_pressure}
+        config.getConfigParameter<double>("critical_pressure");
+
+    auto const omega =
+        //! \ogs_file_param{properties__property__PengRobinson__acentric_factor}
+        config.getConfigParameter<double>("acentric_factor");
+
+    return std::make_unique<PengRobinson>(Tc, pc, omega);
+}
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CreatePengRobinson.h

@@ -0,0 +1,27 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, 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 <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class Property;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<Property> createPengRobinson(BaseLib::ConfigTree const& config);
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CreateProperties.h

@@ -37,6 +37,7 @@
 #include "CreateLinear.h"
 #include "CreateOrthotropicEmbeddedFracturePermeability.h"
 #include "CreateParameter.h"
+#include "CreatePengRobinson.h"
 #include "CreatePermeabilityMohrCoulombFailureIndexModel.h"
 #include "CreatePermeabilityOrthotropicPowerLaw.h"
 #include "CreatePorosityFromMassBalance.h"

MaterialLib/MPL/Properties/PengRobinson.cpp

@@ -0,0 +1,105 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2024, 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 "PengRobinson.h"
+
+#include <algorithm>
+#include <boost/math/special_functions/pow.hpp>
+#include <cmath>
+
+#include "MathLib/Nonlinear/CubicRoots.h"
+
+namespace MaterialPropertyLib
+{
+
+PropertyDataType PengRobinson::value(
+    MaterialPropertyLib::VariableArray const& variable_array,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+    double const /*dt*/) const
+{
+    const double gas_constant = MaterialLib::PhysicalConstant::IdealGasConstant;
+    const double pressure = variable_array.gas_phase_pressure;
+    const double temperature = variable_array.temperature;
+    const double molar_mass = variable_array.molar_mass;
+
+    const double Tr = temperature / Tc_;  // reduced Temperature
+
+    const double kappa = 0.37464 + 1.5422 * omega_ - 0.26992 * omega_ * omega_;
+
+    const double alpha = boost::math::pow<2>(1 + kappa * (1 - std::sqrt(Tr)));
+
+    // EOS in the form: 0 = rho^3 + z1*rho^2 + z2*rho + z3
+    const double denominator =
+        b_ *
+        (pressure * b_ * b_ + b_ * gas_constant * temperature - a_ * alpha);
+
+    const double z1 =
+        (molar_mass * a_ * alpha - 3 * molar_mass * b_ * b_ * pressure -
+         2 * molar_mass * gas_constant * temperature * b_) /
+        denominator;
+    const auto z2 = (molar_mass * molar_mass *
+                     (b_ * pressure - gas_constant * temperature)) /
+                    denominator;
+    const auto z3 =
+        (molar_mass * molar_mass * molar_mass * pressure) / denominator;
+
+    MathLib::CubicSolver cubic_solver_(1., z1, z2, z3);
+    return cubic_solver_.smallestPositiveRealRoot();
+}
+
+PropertyDataType PengRobinson::dValue(
+    MaterialPropertyLib::VariableArray const& variable_array,
+    Variable const variable, ParameterLib::SpatialPosition const& pos,
+    double const t, double const dt) const
+{
+    if (variable == Variable::temperature)
+    {
+        const double temperature = variable_array.temperature;
+        const double epsilon = 1.e-6;
+
+        MaterialPropertyLib::VariableArray perturbed = variable_array;
+        // Increase temperature by +epsilon
+        perturbed.temperature = temperature + epsilon;
+        const double rho_plus = std::get<double>(value(perturbed, pos, t, dt));
+
+        // Decrease temperature by -epsilon
+        perturbed.temperature = temperature - epsilon;
+        const double rho_minus = std::get<double>(value(perturbed, pos, t, dt));
+
+        // Calculate the central difference quotient
+        return (rho_plus - rho_minus) / (2 * epsilon);
+    }
+
+    if (variable == Variable::gas_phase_pressure)
+    {
+        const double pressure = variable_array.gas_phase_pressure;
+        const double epsilon = 1.e-3;
+
+        MaterialPropertyLib::VariableArray perturbed = variable_array;
+        // Increase pressure by +epsilon
+        perturbed.gas_phase_pressure = pressure + epsilon;
+        const double rho_plus = std::get<double>(value(perturbed, pos, t, dt));
+
+        // Decrease pressure by -epsilon
+        perturbed.gas_phase_pressure = pressure - epsilon;
+        const double rho_minus = std::get<double>(value(perturbed, pos, t, dt));
+
+        // Calculate the central difference quotient
+        return (rho_plus - rho_minus) / (2 * epsilon);
+    }
+
+    OGS_FATAL(
+        "PengRobinson::dValue is implemented for derivatives with respect to "
+        "gas phase pressure or temperature only.");
+
+    return 0.;
+}
+
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/PengRobinson.h

@@ -0,0 +1,91 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2024, 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"
+#include "MaterialLib/MPL/VariableType.h"
+#include "MaterialLib/PhysicalConstant.h"
+#include "ParameterLib/Parameter.h"
+
+namespace MaterialPropertyLib
+{
+/**
+ * \brief Peng-Robinson equation of state
+ *
+ * This class implements the Peng-Robinson equation of state (PR-EOS),
+ * a widely used cubic equation of state for describing the behaviour
+ * of real gases, particularly hydrocarbons. It accounts for non-ideal
+ * behaviour of fluids over a range of temperatures and pressures.
+ *
+ * \details
+ * The equation is given in terms of the molar density \f$\rho\f$ as:
+ * \f[
+ * P = \frac{R T \rho}{1 - b \rho} - \frac{a \rho^2}{1 + 2b \rho - b^2 \rho^2}
+ * \f]
+ * where \f$P\f$ is the pressure, \f$T\f$ the temperature, \f$\rho\f$ the molar
+ * density,
+ * \f$R\f$ the universal gas constant, and \f$a\f$, \f$b\f$ are
+ * substance-specific parameters.
+ *
+ * The parameters \f$a\f$ and \f$b\f$ are computed from the critical temperature
+ * \f$T_c\f$ in Kelvin, critical pressure \f$p_c\f$ in Pascal, and
+ * (dimensionless) acentric factor \f$\omega\f$ as follows: \f[ a = 0.457235
+ * \frac{R^2 T_c^2}{p_c} \f] \f[ b = 0.077796 \frac{R T_c}{p_c} \f]
+ *
+ * The Peng-Robinson equation is applicable for a wide range of
+ * substances (gases and liquids), particularly hydrocarbons, at conditions
+ * ranging from subcritical to supercritical. The EOS is not suitable for solid
+ * phases or very low-temperature applications where real gases behave ideally.
+ *
+ * All input parameters (temperature, pressure, density) are assumed to
+ * be in SI units:
+ * - Temperature \f$T\f$ in Kelvin [K]
+ * - Pressure \f$P\f$ in Pascal [Pa]
+ * - Mass density \f$\rho\f$ in \f$[kg/m^3]\f$
+ * - The resulting properties will also follow SI units.
+ *
+ * Original source: D.-Y. Peng and D.B. Robinson, "A New Two-Constant
+ * Equation of State," Industrial & Engineering Chemistry Fundamentals, vol. 15,
+ * pp. 59-64, 1976.
+ */
+class PengRobinson final : public Property
+{
+public:
+    explicit PengRobinson(const double Tc, const double pc, const double omega)
+        : Tc_(Tc), pc_(pc), omega_(omega)
+    {
+        const double gas_constant =
+            MaterialLib::PhysicalConstant::IdealGasConstant;
+        a_ = 0.457235 * gas_constant * gas_constant * Tc_ * Tc_ / pc_;
+        b_ = 0.077796 * gas_constant * Tc_ / pc_;
+    }
+
+    PropertyDataType value(
+        MaterialPropertyLib::VariableArray const& variable_array,
+        ParameterLib::SpatialPosition const& pos, double const t,
+        double const dt) const override;
+    PropertyDataType dValue(
+        MaterialPropertyLib::VariableArray const& variable_array,
+        Variable const variable, ParameterLib::SpatialPosition const& pos,
+        double const t, double const dt) const override;
+
+private:
+    /// Critical temperature
+    const double Tc_;
+    /// Critical pressure
+    const double pc_;
+    /// Acentric factor
+    const double omega_;
+    /// Parameter a (cohesion pressure)
+    double a_;
+    /// Parameter b (covolume)
+    double b_;
+};
+}  // namespace MaterialPropertyLib

MathLib/CMakeLists.txt

@@ -34,6 +34,7 @@ target_link_libraries(
         $<$<BOOL:${OGS_USE_CVODE}>:CVODE::CVODE>
         $<$<BOOL:${OGS_USE_PETSC}>:PkgConfig::PETSC>
         Eigen3::Eigen
+        Boost::math
         $<$<AND:$<TARGET_EXISTS:OpenMP::OpenMP_CXX>,$<BOOL:$<STREQUAL:${OGS_EIGEN_PARALLEL_BACKEND},OpenMP>>>:OpenMP::OpenMP_CXX>
     PRIVATE $<$<BOOL:${OGS_USE_MKL}>:MKL::MKL>
 )

MathLib/Nonlinear/CubicRoots.h

@@ -0,0 +1,75 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, 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 <algorithm>
+#include <limits>
+#include <range/v3/algorithm/min.hpp>
+#include <range/v3/algorithm/sort.hpp>
+#include <range/v3/range/conversion.hpp>
+#include <range/v3/view/filter.hpp>
+#include <vector>
+
+#include "BaseLib/Error.h"
+#include "cubic_roots.hpp"
+namespace MathLib
+{
+
+class CubicSolver
+{
+public:
+    CubicSolver(const double a, const double b, const double c, const double d)
+        : a_(a), b_(b), c_(c), d_(d)
+    {
+        if (std::abs(a_) < 1e-9)
+        {
+            OGS_FATAL("'a' must be non-zero for a cubic equation.");
+        }
+    }
+
+    // Method to solve the cubic equation using Boost
+    std::vector<double> solve()
+    {
+        // Solve using Boost's cubic_roots
+        std::array<double, 3> const roots =
+            boost::math::tools::cubic_roots<double>(a_, b_, c_, d_);
+
+        std::vector<double> filtered_roots =
+            ranges::views::ref(roots) |
+            ranges::views::filter([](double const root)
+                                  { return !std::isnan(root); }) |
+            ranges::to<std::vector>();
+        ranges::sort(filtered_roots);
+
+        return filtered_roots;
+    }
+
+    // Method that returns the smallest positive real root
+    double smallestPositiveRealRoot()
+    {
+        std::vector<double> const roots = solve();
+
+        auto positive_roots =
+            roots | ranges::views::filter([](double root) { return root > 0; });
+
+        // If no positive root exists, return NaN
+        if (ranges::empty(positive_roots))
+        {
+            return std::numeric_limits<double>::quiet_NaN();
+        }
+
+        return ranges::min(positive_roots);
+    }
+
+private:
+    double a_, b_, c_, d_;  // Coefficients of the cubic equation
+};
+
+}  // namespace MathLib