[MPL] added OrthotropicEmbeddedFracturePermeability

Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/c_OrthotropicEmbeddedFracturePermeability.md

+\copydoc MaterialPropertyLib::OrthotropicEmbeddedFracturePermeability

Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_normals.md

+A set of the first two orthogonal fracture normal vectors. The third normal vector is calculated as the cross product.
\ No newline at end of file

Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_rotation_xy.md

+The angle by which the three fracture normals are rotated in the xy-plane.
\ No newline at end of file

Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_rotation_yz.md

+The angle by which the three fracture normals are rotated in the yz-plane.
\ No newline at end of file

Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_intrinsic_permeability.md

+The permeability of the undisturbed material.

Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_mean_frac_distances.md

+The mean distance between neighboring fractures for each fracture plane.

Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_threshold_strains.md

+Threshold strain for each fracture plane, which has to be exceeded to create additional permeability due to fracture opening.

MaterialLib/MPL/CreateProperty.cpp

@@ -104,6 +104,12 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
         return createEmbeddedFracturePermeability(geometry_dimension, config);
     }
 
+    if (boost::iequals(property_type, "OrthotropicEmbeddedFracturePermeability"))
+    {
+        return createOrthotropicEmbeddedFracturePermeability(
+            geometry_dimension, config, parameters);
+    }
+
     if (boost::iequals(property_type,
                        "PermeabilityMohrCoulombFailureIndexModel"))
     {

MaterialLib/MPL/Properties/CreateOrthotropicEmbeddedFracturePermeability.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
+ *
+ */
+
+#include "BaseLib/ConfigTree.h"
+#include "OrthotropicEmbeddedFracturePermeability.h"
+#include "ParameterLib/Utils.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<Property> createOrthotropicEmbeddedFracturePermeability(
+    int const geometry_dimension, BaseLib::ConfigTree const& config,
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters)
+{
+    if ((geometry_dimension != 2) && (geometry_dimension != 3))
+    {
+        OGS_FATAL(
+            "The OrthotropicEmbeddedFracturePermeability is implemented only "
+            "for 2D or 3D problems");
+    }
+
+    //! \ogs_file_param{properties__property__type}
+    config.checkConfigParameter("type",
+                                "OrthotropicEmbeddedFracturePermeability");
+
+    // Second access for storage.
+    //! \ogs_file_param{properties__property__name}
+    auto property_name = config.peekConfigParameter<std::string>("name");
+
+    DBUG("Create OrthotropicEmbeddedFracturePermeability medium property");
+
+    auto const a_i =
+        //! \ogs_file_param{properties__property__OrthotropicEmbeddedFracturePermeability__mean_frac_distances}
+        config.getConfigParameter<std::vector<double>>("mean_frac_distances");
+    if (a_i.size() != 3)
+    {
+        OGS_FATAL(
+            "The size of the mean fracture distances vector must be 3, but is "
+            "{}.",
+            a_i.size());
+    }
+
+    auto const e_i0 =
+        //! \ogs_file_param{properties__property__OrthotropicEmbeddedFracturePermeability__threshold_strains}
+        config.getConfigParameter<std::vector<double>>("threshold_strains");
+    if (e_i0.size() != 3)
+    {
+        OGS_FATAL(
+            "The size of the mean threshold strains vector must be 3, but is "
+            "{}.",
+            e_i0.size());
+    }
+
+    auto const n =
+        //! \ogs_file_param{properties__property__OrthotropicEmbeddedFracturePermeability__fracture_normals}
+        config.getConfigParameter<std::vector<double>>("fracture_normals");
+    if (n.size() != 6)
+    {
+        OGS_FATAL(
+            "The size of the fracture normals vector must be 6, but is {}.",
+            n.size());
+    }
+    Eigen::Vector3d const n1 = Eigen::Vector3d({n[0], n[1], n[2]}).normalized();
+    Eigen::Vector3d const n2 = Eigen::Vector3d({n[3], n[4], n[5]}).normalized();
+
+    if (n1.dot(n2) > std::numeric_limits<double>::epsilon())
+    {
+        OGS_FATAL(
+            "The given fracture normals are not orthogonal. Please provide two "
+            "orthogonal fracture normals");
+    }
+
+    Eigen::Matrix3d const n_i =
+        (Eigen::Matrix3d() << n1, n2, n1.cross(n2)).finished();
+
+    std::string const intrinsic_permeability_param_name =
+        //! \ogs_file_param{properties__property__OrthotropicEmbeddedFracturePermeability__intrinsic_permeability}
+        config.getConfigParameter<std::string>("intrinsic_permeability");
+
+    auto const& k = ParameterLib::findParameter<double>(
+        intrinsic_permeability_param_name, parameters, 0, nullptr);
+
+    std::string const fracture_rotation_xy_param_name =
+        //! \ogs_file_param{properties__property__OrthotropicEmbeddedFracturePermeability__fracture_rotation_xy}
+        config.getConfigParameter<std::string>("fracture_rotation_xy");
+
+    auto const& phi_xy = ParameterLib::findParameter<double>(
+        fracture_rotation_xy_param_name, parameters, 0, nullptr);
+
+    std::string const fracture_rotation_yz_param_name =
+        //! \ogs_file_param{properties__property__OrthotropicEmbeddedFracturePermeability__fracture_rotation_yz}
+        config.getConfigParameter<std::string>("fracture_rotation_yz");
+
+    auto const& phi_yz = ParameterLib::findParameter<double>(
+        fracture_rotation_yz_param_name, parameters, 0, nullptr);
+
+    if (geometry_dimension == 2)
+    {
+        return std::make_unique<OrthotropicEmbeddedFracturePermeability<2>>(
+            std::move(property_name), a_i, e_i0, n_i, k, phi_xy, phi_yz);
+    }
+    return std::make_unique<OrthotropicEmbeddedFracturePermeability<3>>(
+        std::move(property_name), a_i, e_i0, n_i, k, phi_xy, phi_yz);
+}
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CreateOrthotropicEmbeddedFracturePermeability.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
+ *
+ */
+
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<Property> createOrthotropicEmbeddedFracturePermeability(
+    int const geometry_dimension, BaseLib::ConfigTree const& config,
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
+        parameters);
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CreateProperties.h

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

MaterialLib/MPL/Properties/OrthotropicEmbeddedFracturePermeability.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
+ *
+ */
+
+#include "MaterialLib/MPL/Properties/OrthotropicEmbeddedFracturePermeability.h"
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
+
+namespace MaterialPropertyLib
+{
+template <int DisplacementDim>
+OrthotropicEmbeddedFracturePermeability<DisplacementDim>::
+    OrthotropicEmbeddedFracturePermeability(
+        std::string name,
+        std::vector<double> const& mean_fracture_distances,
+        std::vector<double> const& threshold_strains,
+        Eigen::Matrix<double, 3, 3> const fracture_normals,
+        ParameterLib::Parameter<double> const& intrinsic_permeability,
+        ParameterLib::Parameter<double> const& fracture_rotation_xy,
+        ParameterLib::Parameter<double> const& fracture_rotation_yz)
+    : _a(mean_fracture_distances),
+      _e0(threshold_strains),
+      _n(fracture_normals),
+      _k(intrinsic_permeability),
+      _phi_xy(fracture_rotation_xy),
+      _phi_yz(fracture_rotation_yz)
+{
+    name_ = std::move(name);
+}
+
+template <int DisplacementDim>
+void OrthotropicEmbeddedFracturePermeability<DisplacementDim>::checkScale()
+    const
+{
+    if (!std::holds_alternative<Medium*>(scale_))
+    {
+        OGS_FATAL(
+            "The property 'OrthotropicEmbeddedFracturePermeability' is "
+            "implemented on the 'media' scale only.");
+    }
+}
+
+template <int DisplacementDim>
+PropertyDataType
+OrthotropicEmbeddedFracturePermeability<DisplacementDim>::value(
+    VariableArray const& variable_array,
+    ParameterLib::SpatialPosition const& pos, double const t,
+    double const /*dt*/) const
+{
+    auto const eps = formEigenTensor<3>(std::get<SymmetricTensor>(
+        variable_array[static_cast<int>(Variable::mechanical_strain)]));
+    double const k = std::get<double>(fromVector(_k(t, pos)));
+    double const _b0 = std::sqrt(12.0 * k);
+
+    double const phi_xy = std::get<double>(fromVector(_phi_xy(t, pos)));
+    double const phi_yz = std::get<double>(fromVector(_phi_yz(t, pos)));
+    auto const rotMat_xy = Eigen::AngleAxisd(phi_xy, Eigen::Vector3d::UnitZ());
+    auto const rotMat_yz = Eigen::AngleAxisd(phi_yz, Eigen::Vector3d::UnitX());
+
+    Eigen::Matrix3d result = Eigen::Vector3d::Constant(k).asDiagonal();
+
+    for (int i = 0; i < 3; i++)
+    {
+        Eigen::Vector3d const n_i = rotMat_yz * (rotMat_xy * _n.col(i));
+        double const e_n = (eps * n_i).dot(n_i.transpose());
+        double const H_de = (e_n > _e0[i]) ? 1.0 : 0.0;
+        double const b_f = _b0 + H_de * _a[i] * (e_n - _e0[i]);
+
+        result += H_de * (b_f / _a[i]) * ((b_f * b_f / 12.0) - k) *
+                  (Eigen::Matrix3d::Identity() - n_i * n_i.transpose());
+    }
+    return result.template topLeftCorner<DisplacementDim, DisplacementDim>()
+        .eval();
+}
+template <int DisplacementDim>
+PropertyDataType
+OrthotropicEmbeddedFracturePermeability<DisplacementDim>::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::mechanical_strain) &&
+           "OrthotropicEmbeddedFracturePermeability::dValue is implemented for "
+           " derivatives with respect to strain only.");
+
+    auto const eps = formEigenTensor<3>(std::get<SymmetricTensor>(
+        variable_array[static_cast<int>(Variable::mechanical_strain)]));
+    double const k = std::get<double>(fromVector(_k(t, pos)));
+    double const _b0 = std::sqrt(12.0 * k);
+
+    double const phi_xy = std::get<double>(fromVector(_phi_xy(t, pos)));
+    double const phi_yz = std::get<double>(fromVector(_phi_yz(t, pos)));
+    auto const rotMat_xy = Eigen::AngleAxisd(phi_xy, Eigen::Vector3d::UnitZ());
+    auto const rotMat_yz = Eigen::AngleAxisd(phi_yz, Eigen::Vector3d::UnitX());
+
+    Eigen::Matrix3d result = Eigen::Matrix3d::Zero();
+
+    for (int i = 0; i < 3; i++)
+    {
+        Eigen::Vector3d const n_i = rotMat_yz * (rotMat_xy * _n.col(i));
+        Eigen::Matrix3d const M = n_i * n_i.transpose();
+        double const e_n = (eps * n_i).dot(n_i.transpose());
+        double const H_de = (e_n > _e0[i]) ? 1.0 : 0.0;
+        double const b_f = _b0 + H_de * _a[i] * (e_n - _e0[i]);
+
+        result += H_de * (b_f * b_f / 4.0 - k) *
+                  (Eigen::Matrix3d::Identity() - M) * M;
+    }
+    return result.template topLeftCorner<DisplacementDim, DisplacementDim>()
+        .eval();
+}
+
+template class OrthotropicEmbeddedFracturePermeability<2>;
+template class OrthotropicEmbeddedFracturePermeability<3>;
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/OrthotropicEmbeddedFracturePermeability.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
+ *
+ */
+#pragma once
+
+#include "MaterialLib/MPL/Property.h"
+#include "MathLib/KelvinVector.h"
+#include "ParameterLib/Parameter.h"
+
+namespace MaterialPropertyLib
+{
+class Medium;
+class Phase;
+class Component;
+/**
+ * \class OrthotropicEmbeddedFracturePermeability
+ * \brief Extended Permeability model based on Olivella&Alonso
+ * \details This property must be a medium property, it
+ * computes the permeability in dependence of the strain
+ *
+ * The base model was proposed
+ * in \cite alonso2006mechanisms and it was further investigated
+ * in \cite olivella2008gas . This extended Version features three
+ * orthotropic fracture planes.
+ *
+ * The model takes the form of
+ * \f[ \mathbf{k} = k_\text{m} \mathbf{I} + \sum \limits_{i=1}^3
+ * \frac{b_i}{a_i} \left( \frac{b_i^2}{12} - k_\text{m} \right) \left(
+ * \mathbf{I} - \mathbf{M}_i \right) \f]
+ * with
+ * \f[ \mathbf{M}_i = \vec{n}_i \otimes \vec{n}_i \f]
+ * and
+ * \f[ b_i = b_{i0} + \Delta b_i \\
+ * \Delta b_i = a_i \langle \mathbf{\epsilon} : \mathbf{M}_i -
+ * \varepsilon_{0i} \rangle
+ * \f]
+ * where
+ * <table>
+ * <tr><td>\f$ k_\text{m} \f$  <td> permeability of undisturbed material
+ * <tr><td>\f$ b_i \f$  <td> fracture aperture of each fracture plane
+ * <tr><td>\f$ b_{i0} \f$  <td> initial aperture of each fracture plane
+ * <tr><td>\f$ a_i \f$  <td> mean fracture distance of each fracture plane
+ * <tr><td>\f$ \vec{n}_i \f$  <td> fracture normal vector of each fracture
+ * plane
+ * <tr><td>\f$ \varepsilon_{i0} \f$  <td> threshold strain of each fracture
+ * plane
+ * </table>
+ */
+
+template <int DisplacementDim>
+class OrthotropicEmbeddedFracturePermeability final : public Property
+{
+private:
+    Medium* _medium = nullptr;
+    std::vector<double> const _a;
+    std::vector<double> const _e0;
+    Eigen::Matrix<double, 3, 3> const _n;
+    ParameterLib::Parameter<double> const& _k;
+    ParameterLib::Parameter<double> const& _phi_xy;
+    ParameterLib::Parameter<double> const& _phi_yz;
+
+public:
+    OrthotropicEmbeddedFracturePermeability(
+        std::string name,
+        std::vector<double> const& mean_fracture_distances,
+        std::vector<double> const& threshold_strains,
+        Eigen::Matrix<double, 3, 3> const fracture_normals,
+        ParameterLib::Parameter<double> const& intrinsic_permeability,
+        ParameterLib::Parameter<double> const& fracture_rotation_xy,
+        ParameterLib::Parameter<double> const& fracture_rotation_yz);
+
+    using SymmetricTensor = Eigen::Matrix<
+        double,
+        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim), 1>;
+
+    void checkScale() const override;
+
+    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

MaterialLib/MPL/Properties/Properties.h

@@ -32,6 +32,7 @@
 #include "IdealGasLaw.h"
 #include "Linear.h"
 #include "Parameter.h"
+#include "OrthotropicEmbeddedFracturePermeability.h"
 #include "PorosityFromMassBalance.h"
 #include "RelativePermeability/RelPermBrooksCorey.h"
 #include "RelativePermeability/RelPermBrooksCoreyNonwettingPhase.h"