13#include <Eigen/Eigenvalues>
14#include <Eigen/Geometry>
21template <
int DisplacementDim>
24 Eigen::Matrix<double, 3, 1>
const fracture_normal,
25 bool const fracture_normal_is_constant,
26 double const intrinsic_permeability,
27 double const initial_aperture,
28 double const mean_fracture_distance,
29 double const threshold_strain,
32 double const jacobian_factor)
33 : _n(fracture_normal),
34 _n_const(fracture_normal_is_constant),
35 _k(intrinsic_permeability),
36 _b0(initial_aperture),
37 _a(mean_fracture_distance),
38 _e0(threshold_strain),
39 _phi_xy(fracture_rotation_xy),
40 _phi_yz(fracture_rotation_yz),
46template <
int DisplacementDim>
49 if (!std::holds_alternative<Medium*>(scale_))
52 "The property 'EmbeddedFracturePermeability' is "
53 "implemented on the 'medium' scale only.");
57template <
int DisplacementDim>
63 Eigen::Matrix<double, 3, 1>
const n = [&]
71 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(sigma);
72 return (Eigen::Matrix<double, 3, 1>)e_s.eigenvectors().col(2);
75 auto const rotMat_xy =
76 Eigen::AngleAxisd(_phi_xy(t, pos)[0], Eigen::Vector3d::UnitZ());
77 auto const rotMat_yz =
78 Eigen::AngleAxisd(_phi_yz(t, pos)[0], Eigen::Vector3d::UnitX());
80 Eigen::Matrix<double, 3, 1>
const n_r = rotMat_yz * (rotMat_xy * n);
85 double const e_n = (eps * n_r).dot(n_r.transpose());
86 double const H_de = (e_n > _e0) ? 1.0 : 0.0;
87 double const b_f = _b0 + H_de * _a * (e_n - _e0);
88 double const coeff = H_de * (b_f / _a) * ((b_f * b_f / 12.0) - _k);
90 Eigen::Matrix3d I = Eigen::Matrix3d::Identity();
91 return (coeff * (I - n_r * n_r.transpose()) + _k * I)
92 .template topLeftCorner<DisplacementDim, DisplacementDim>()
96template <
int DisplacementDim>
105 "EmbeddedFracturePermeability::dValue is implemented for "
106 "derivatives with respect to strain only.");
109 Eigen::Matrix<double, 3, 1>
const n = [&]
116 std::get<SymmetricTensor>(variable_array.
total_stress));
117 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(sigma);
118 return (Eigen::Matrix<double, 3, 1>)e_s.eigenvectors().col(2);
121 auto const rotMat_xy =
122 Eigen::AngleAxisd(_phi_xy(t, pos)[0], Eigen::Vector3d::UnitZ());
123 auto const rotMat_yz =
124 Eigen::AngleAxisd(_phi_yz(t, pos)[0], Eigen::Vector3d::UnitX());
126 Eigen::Matrix<double, 3, 1>
const n_r = rotMat_yz * (rotMat_xy * n);
131 double const e_n = (eps * n_r).dot(n_r.transpose());
132 double const H_de = (e_n > _e0) ? 1.0 : 0.0;
133 double const b_f = _b0 + H_de * _a * (e_n - _e0);
135 Eigen::Matrix3d
const M = n_r * n_r.transpose();
136 return Eigen::MatrixXd(
137 _jf * H_de * (b_f * b_f / 4 - _k) *
138 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
139 Eigen::Matrix3d::Identity() - M) *
140 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(M).transpose());
Permeability model as proposed by Olivella&Alonso.
EmbeddedFracturePermeability(std::string name, Eigen::Matrix< double, 3, 1 > const fracture_normal, bool const fracture_normal_is_constant, double const intrinsic_permeability, 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)
PropertyDataType dValue(VariableArray const &variable_array, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const override
void checkScale() const override
virtual PropertyDataType value() const
KelvinVector mechanical_strain
KelvinVector total_stress
template Eigen::Matrix< double, 3, 3 > formEigenTensor< 3 >(MaterialPropertyLib::PropertyDataType const &values)
std::variant< double, Eigen::Matrix< double, 2, 1 >, Eigen::Matrix< double, 3, 1 >, Eigen::Matrix< double, 2, 2 >, Eigen::Matrix< double, 3, 3 >, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 >, Eigen::MatrixXd > PropertyDataType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Eigen::Matrix< double, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)