Skip to content
GitLab
Commit 1f4213f7 authored by Florian Zill's avatar Florian Zill Committed by Dmitry Yu. Naumov
Browse files

[MPL] enable fracture normal rotation with EmbeddedFracturePermeability 

the same as in OrthotropicEmbeddedFracturePermeability
parent e72a8da2
Hide whitespace changes
Inline Side-by-side
Documentation/ProjectFile/properties/property/EmbeddedFracturePermeability/t_fracture_rotation_xy.md 0 → 100644
View file @ 1f4213f7
The angle in radians by which the fracture normal is rotated in the xy-plane.
\ No newline at end of file
Documentation/ProjectFile/properties/property/EmbeddedFracturePermeability/t_fracture_rotation_yz.md 0 → 100644
View file @ 1f4213f7
The angle in radians by which the fracture normal is rotated in the yz-plane.
\ No newline at end of file
Documentation/ProjectFile/properties/property/OrthotropicEmbeddedFracturePermeability/t_fracture_rotation_xy.md
View file @ 1f4213f7
The angle by which the three fracture normals are rotated in the xy-plane.
\ No newline at end of file
The angle in radians 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
View file @ 1f4213f7
The angle by which the three fracture normals are rotated in the yz-plane.
\ No newline at end of file
The angle in radians by which the three fracture normals are rotated in the yz-plane.
\ No newline at end of file
MaterialLib/MPL/CreateProperty.cpp
View file @ 1f4213f7
......@@ -112,7 +112,8 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
if (boost::iequals(property_type, "EmbeddedFracturePermeability"))
{
return createEmbeddedFracturePermeability(geometry_dimension, config);
return createEmbeddedFracturePermeability(geometry_dimension, config,
parameters);
}
if (boost::iequals(property_type,
......
MaterialLib/MPL/Properties/CreateEmbeddedFracturePermeability.cpp
View file @ 1f4213f7
......@@ -10,12 +10,13 @@
#include "BaseLib/ConfigTree.h"
#include "EmbeddedFracturePermeability.h"
#include "ParameterLib/Parameter.h"
#include "ParameterLib/Utils.h"
namespace MaterialPropertyLib
{
std::unique_ptr<Property> createEmbeddedFracturePermeability(
int const geometry_dimension, BaseLib::ConfigTree const& config)
int const geometry_dimension, BaseLib::ConfigTree const& config,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters)
{
if ((geometry_dimension != 2) && (geometry_dimension != 3))
{
......@@ -74,6 +75,20 @@ std::unique_ptr<Property> createEmbeddedFracturePermeability(
"determined as the third principal stress vector.");
}
std::string const fracture_rotation_xy_param_name =
//! \ogs_file_param{properties__property__EmbeddedFracturePermeability__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__EmbeddedFracturePermeability__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);
auto const jf =
//! \ogs_file_param{properties__property__EmbeddedFracturePermeability__jacobian_factor}
config.getConfigParameter<double>("jacobian_factor", 0.);
......@@ -81,9 +96,10 @@ std::unique_ptr<Property> createEmbeddedFracturePermeability(
if (geometry_dimension == 2)
{
return std::make_unique<EmbeddedFracturePermeability<2>>(
std::move(property_name), n, n_const, k, b0, a, e0, jf);
std::move(property_name), n, n_const, k, b0, a, e0, phi_xy, phi_yz,
jf);
}
return std::make_unique<EmbeddedFracturePermeability<3>>(
std::move(property_name), n, n_const, k, b0, a, e0, jf);
std::move(property_name), n, n_const, k, b0, a, e0, phi_xy, phi_yz, jf);
}
} // namespace MaterialPropertyLib
MaterialLib/MPL/Properties/CreateEmbeddedFracturePermeability.h
View file @ 1f4213f7
......@@ -22,5 +22,7 @@ namespace MaterialPropertyLib
class Property;
std::unique_ptr<Property> createEmbeddedFracturePermeability(
int const geometry_dimension, BaseLib::ConfigTree const& config);
int const geometry_dimension, BaseLib::ConfigTree const& config,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
parameters);
} // namespace MaterialPropertyLib
MaterialLib/MPL/Properties/EmbeddedFracturePermeability.cpp
View file @ 1f4213f7
......@@ -12,7 +12,6 @@
#include "MaterialLib/MPL/Medium.h"
#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
#include "MathLib/KelvinVector.h"
namespace MaterialPropertyLib
{
......@@ -25,6 +24,8 @@ EmbeddedFracturePermeability<DisplacementDim>::EmbeddedFracturePermeability(
double const initial_aperture,
double const mean_fracture_distance,
double const threshold_strain,
ParameterLib::Parameter<double> const& fracture_rotation_xy,
ParameterLib::Parameter<double> const& fracture_rotation_yz,
double const jacobian_factor)
: _n(fracture_normal),
_n_const(fracture_normal_is_constant),
......@@ -32,6 +33,8 @@ EmbeddedFracturePermeability<DisplacementDim>::EmbeddedFracturePermeability(
_b0(initial_aperture),
_a(mean_fracture_distance),
_e0(threshold_strain),
_phi_xy(fracture_rotation_xy),
_phi_yz(fracture_rotation_yz),
_jf(jacobian_factor)
{
name_ = std::move(name);
......@@ -51,7 +54,7 @@ void EmbeddedFracturePermeability<DisplacementDim>::checkScale() const
template <int DisplacementDim>
PropertyDataType EmbeddedFracturePermeability<DisplacementDim>::value(
VariableArray const& variable_array,
ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
ParameterLib::SpatialPosition const& pos, double const t,
double const /*dt*/) const
{
Eigen::Matrix<double, 3, 1> const n = [&]
......@@ -66,16 +69,23 @@ PropertyDataType EmbeddedFracturePermeability<DisplacementDim>::value(
return (Eigen::Matrix<double, 3, 1>)e_s.eigenvectors().col(2);
}();
auto const rotMat_xy =
Eigen::AngleAxisd(_phi_xy(t, pos)[0], Eigen::Vector3d::UnitZ());
auto const rotMat_yz =
Eigen::AngleAxisd(_phi_yz(t, pos)[0], Eigen::Vector3d::UnitX());
Eigen::Matrix<double, 3, 1> const n_r = rotMat_yz * (rotMat_xy * n);
auto const eps = MathLib::KelvinVector::kelvinVectorToTensor(
std::get<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
variable_array[static_cast<int>(Variable::mechanical_strain)]));
double const e_n = (eps * n).dot(n.transpose());
double const e_n = (eps * n_r).dot(n_r.transpose());
double const H_de = (e_n > _e0) ? 1.0 : 0.0;
double const b_f = _b0 + H_de * _a * (e_n - _e0);
double const coeff = H_de * (b_f / _a) * ((b_f * b_f / 12.0) - _k);
Eigen::Matrix3d I = Eigen::Matrix3d::Identity();
return (coeff * (I - n * n.transpose()) + _k * I)
return (coeff * (I - n_r * n_r.transpose()) + _k * I)
.template topLeftCorner<DisplacementDim, DisplacementDim>()
.eval();
}
......@@ -83,7 +93,7 @@ PropertyDataType EmbeddedFracturePermeability<DisplacementDim>::value(
template <int DisplacementDim>
PropertyDataType EmbeddedFracturePermeability<DisplacementDim>::dValue(
VariableArray const& variable_array, Variable const primary_variable,
ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
ParameterLib::SpatialPosition const& pos, double const t,
double const /*dt*/) const
{
if (primary_variable != Variable::mechanical_strain)
......@@ -105,14 +115,21 @@ PropertyDataType EmbeddedFracturePermeability<DisplacementDim>::dValue(
return (Eigen::Matrix<double, 3, 1>)e_s.eigenvectors().col(2);
}();
auto const rotMat_xy =
Eigen::AngleAxisd(_phi_xy(t, pos)[0], Eigen::Vector3d::UnitZ());
auto const rotMat_yz =
Eigen::AngleAxisd(_phi_yz(t, pos)[0], Eigen::Vector3d::UnitX());
Eigen::Matrix<double, 3, 1> const n_r = rotMat_yz * (rotMat_xy * n);
auto const eps = MathLib::KelvinVector::kelvinVectorToTensor(
std::get<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
variable_array[static_cast<int>(Variable::mechanical_strain)]));
double const e_n = (eps * n).dot(n.transpose());
double const e_n = (eps * n_r).dot(n_r.transpose());
double const H_de = (e_n > _e0) ? 1.0 : 0.0;
double const b_f = _b0 + H_de * _a * (e_n - _e0);
Eigen::Matrix3d const M = n * n.transpose();
Eigen::Matrix3d const M = n_r * n_r.transpose();
return Eigen::MatrixXd(
_jf * H_de * (b_f * b_f / 4 - _k) *
MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
......
MaterialLib/MPL/Properties/EmbeddedFracturePermeability.h
View file @ 1f4213f7
......@@ -11,6 +11,7 @@
#include "MaterialLib/MPL/Property.h"
#include "MathLib/KelvinVector.h"
#include "ParameterLib/Parameter.h"
namespace MaterialPropertyLib
{
......@@ -56,6 +57,8 @@ private:
double const _b0;
double const _a;
double const _e0;
ParameterLib::Parameter<double> const& _phi_xy;
ParameterLib::Parameter<double> const& _phi_yz;
double const _jf;
public:
......@@ -67,6 +70,8 @@ public:
double const initial_aperture,
double const mean_fracture_distance,
double const threshold_strain,
ParameterLib::Parameter<double> const& fracture_rotation_xy,
ParameterLib::Parameter<double> const& fracture_rotation_yz,
double const jacobian_factor);
static int const KelvinVectorSize =
......
Tests/Data/HydroMechanics/EmbeddedFracturePermeability/cube.prj
View file @ 1f4213f7
......@@ -77,7 +77,9 @@
<initial_aperture>0</initial_aperture>
<mean_frac_distance>0.01</mean_frac_distance>
<threshold_strain>1e-5</threshold_strain>
<fracture_normal>0 0 1</fracture_normal>
<fracture_normal>1 0 1</fracture_normal>
<fracture_rotation_xy>pi_over_2</fracture_rotation_xy>
<fracture_rotation_yz>pi_over_4</fracture_rotation_yz>
</property>
</properties>
</medium>
......@@ -163,6 +165,21 @@
<type>Constant</type>
<value>1e-4</value>
</parameter>
<parameter>
<name>zero</name>
<type>Constant</type>
<value>0</value>
</parameter>
<parameter>
<name>pi_over_4</name>
<type>Constant</type>
<value>0.78539816339744831</value>
</parameter>
<parameter>
<name>pi_over_2</name>
<type>Constant</type>
<value>1.5707963267948966</value>
</parameter>
</parameters>
<process_variables>
<process_variable>
......
Tests/Data/HydroMechanics/EmbeddedFracturePermeability/square.prj
View file @ 1f4213f7
......@@ -79,6 +79,8 @@
<mean_frac_distance>0.01</mean_frac_distance>
<threshold_strain>1e-5</threshold_strain>
<fracture_normal>1 0 0</fracture_normal>
<fracture_rotation_xy>zero</fracture_rotation_xy>
<fracture_rotation_yz>zero</fracture_rotation_yz>
</property>
</properties>
</medium>
......@@ -165,6 +167,11 @@
<type>Constant</type>
<value>1e-4</value>
</parameter>
<parameter>
<name>zero</name>
<type>Constant</type>
<value>0</value>
</parameter>
</parameters>
<process_variables>
<process_variable>
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment