From 19e6e1c99f82fa55a3c55e7b1bfc55a8da0e9a2a Mon Sep 17 00:00:00 2001 From: FZill <florian.zill@ufz.de> Date: Wed, 21 Apr 2021 14:30:11 +0200 Subject: [PATCH] [MPL] added OrthotropicEmbeddedFracturePermeability --- ...OrthotropicEmbeddedFracturePermeability.md | 1 + .../t_fracture_normals.md | 1 + .../t_fracture_rotation_xy.md | 1 + .../t_fracture_rotation_yz.md | 1 + .../t_intrinsic_permeability.md | 1 + .../t_mean_frac_distances.md | 1 + .../t_threshold_strains.md | 1 + MaterialLib/MPL/CreateProperty.cpp | 6 + ...rthotropicEmbeddedFracturePermeability.cpp | 111 ++++++++++++++++ ...eOrthotropicEmbeddedFracturePermeability.h | 26 ++++ MaterialLib/MPL/Properties/CreateProperties.h | 1 + ...rthotropicEmbeddedFracturePermeability.cpp | 123 ++++++++++++++++++ .../OrthotropicEmbeddedFracturePermeability.h | 93 +++++++++++++ MaterialLib/MPL/Properties/Properties.h | 1 + 14 files changed, 368 insertions(+) create mode 100644 Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/c_OrthotropicEmbeddedFracturePermeability.md create mode 100644 Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_normals.md create mode 100644 Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_rotation_xy.md create mode 100644 Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_rotation_yz.md create mode 100644 Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_intrinsic_permeability.md create mode 100644 Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_mean_frac_distances.md create mode 100644 Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_threshold_strains.md create mode 100644 MaterialLib/MPL/Properties/CreateOrthotropicEmbeddedFracturePermeability.cpp create mode 100644 MaterialLib/MPL/Properties/CreateOrthotropicEmbeddedFracturePermeability.h create mode 100644 MaterialLib/MPL/Properties/OrthotropicEmbeddedFracturePermeability.cpp create mode 100644 MaterialLib/MPL/Properties/OrthotropicEmbeddedFracturePermeability.h diff --git a/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/c_OrthotropicEmbeddedFracturePermeability.md b/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/c_OrthotropicEmbeddedFracturePermeability.md new file mode 100644 index 00000000000..cee45d4fc97 --- /dev/null +++ b/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/c_OrthotropicEmbeddedFracturePermeability.md @@ -0,0 +1 @@ +\copydoc MaterialPropertyLib::OrthotropicEmbeddedFracturePermeability diff --git a/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_normals.md b/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_normals.md new file mode 100644 index 00000000000..2676df188ba --- /dev/null +++ b/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_normals.md @@ -0,0 +1 @@ +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 diff --git a/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_rotation_xy.md b/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_rotation_xy.md new file mode 100644 index 00000000000..3bfb6ff085c --- /dev/null +++ b/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_rotation_xy.md @@ -0,0 +1 @@ +The angle by which the three fracture normals are rotated in the xy-plane. \ No newline at end of file diff --git a/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_rotation_yz.md b/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_rotation_yz.md new file mode 100644 index 00000000000..67c3459530a --- /dev/null +++ b/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_rotation_yz.md @@ -0,0 +1 @@ +The angle by which the three fracture normals are rotated in the yz-plane. \ No newline at end of file diff --git a/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_intrinsic_permeability.md b/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_intrinsic_permeability.md new file mode 100644 index 00000000000..1c710341f15 --- /dev/null +++ b/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_intrinsic_permeability.md @@ -0,0 +1 @@ +The permeability of the undisturbed material. diff --git a/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_mean_frac_distances.md b/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_mean_frac_distances.md new file mode 100644 index 00000000000..0622e8b190c --- /dev/null +++ b/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_mean_frac_distances.md @@ -0,0 +1 @@ +The mean distance between neighboring fractures for each fracture plane. diff --git a/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_threshold_strains.md b/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_threshold_strains.md new file mode 100644 index 00000000000..8bf090a1869 --- /dev/null +++ b/Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_threshold_strains.md @@ -0,0 +1 @@ +Threshold strain for each fracture plane, which has to be exceeded to create additional permeability due to fracture opening. diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp index 76975329ab4..3f845abd5fb 100644 --- a/MaterialLib/MPL/CreateProperty.cpp +++ b/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")) { diff --git a/MaterialLib/MPL/Properties/CreateOrthotropicEmbeddedFracturePermeability.cpp b/MaterialLib/MPL/Properties/CreateOrthotropicEmbeddedFracturePermeability.cpp new file mode 100644 index 00000000000..abaae34dccc --- /dev/null +++ b/MaterialLib/MPL/Properties/CreateOrthotropicEmbeddedFracturePermeability.cpp @@ -0,0 +1,111 @@ +/** + * \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 diff --git a/MaterialLib/MPL/Properties/CreateOrthotropicEmbeddedFracturePermeability.h b/MaterialLib/MPL/Properties/CreateOrthotropicEmbeddedFracturePermeability.h new file mode 100644 index 00000000000..a5d4b9ee10a --- /dev/null +++ b/MaterialLib/MPL/Properties/CreateOrthotropicEmbeddedFracturePermeability.h @@ -0,0 +1,26 @@ +/** + * \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 diff --git a/MaterialLib/MPL/Properties/CreateProperties.h b/MaterialLib/MPL/Properties/CreateProperties.h index 5cf7b9628bc..7b45b0a11a4 100644 --- a/MaterialLib/MPL/Properties/CreateProperties.h +++ b/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" diff --git a/MaterialLib/MPL/Properties/OrthotropicEmbeddedFracturePermeability.cpp b/MaterialLib/MPL/Properties/OrthotropicEmbeddedFracturePermeability.cpp new file mode 100644 index 00000000000..b1470477bd7 --- /dev/null +++ b/MaterialLib/MPL/Properties/OrthotropicEmbeddedFracturePermeability.cpp @@ -0,0 +1,123 @@ +/** + * \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 diff --git a/MaterialLib/MPL/Properties/OrthotropicEmbeddedFracturePermeability.h b/MaterialLib/MPL/Properties/OrthotropicEmbeddedFracturePermeability.h new file mode 100644 index 00000000000..14aa85c204c --- /dev/null +++ b/MaterialLib/MPL/Properties/OrthotropicEmbeddedFracturePermeability.h @@ -0,0 +1,93 @@ +/** + * \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 diff --git a/MaterialLib/MPL/Properties/Properties.h b/MaterialLib/MPL/Properties/Properties.h index 4af58ab9e47..f22a6061d3c 100644 --- a/MaterialLib/MPL/Properties/Properties.h +++ b/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" -- GitLab