14 #include <boost/math/constants/constants.hpp>
29 template <
int DisplacementDim>
33 double const kr,
double const b,
double const c,
double const phi,
34 double const k_max,
double const t_sigma_max,
40 phi_(boost::math::constants::degree<double>() * phi),
42 t_sigma_max_(t_sigma_max),
43 local_coordinate_system_(local_coordinate_system)
45 const double t_sigma_upper =
c_ / std::tan(
phi_);
46 if (t_sigma_max_ <= 0.0 || t_sigma_max_ > t_sigma_upper ||
48 std::numeric_limits<double>::epsilon())
51 "Tensile strength parameter of {:e} is out of the range (0, "
52 "c/tan(phi)) = (0, {:e})",
59 template <
int DisplacementDim>
63 if (!std::holds_alternative<Medium*>(scale_))
66 "The property 'PermeabilityMohrCoulombFailureIndexModel' is "
67 "implemented on the 'medium' scale only.");
71 template <
int DisplacementDim>
78 auto const& stress_vector = std::get<SymmetricTensor<DisplacementDim>>(
81 auto const& stress_tensor =
84 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>>
85 eigenvalue_solver(stress_tensor);
88 auto const sigma = eigenvalue_solver.eigenvalues();
90 auto k_data = k0_(t, pos);
92 double const max_sigma = std::max(std::fabs(sigma[0]), std::fabs(sigma[2]));
94 if (max_sigma < std::numeric_limits<double>::epsilon())
99 double const sigma_m = 0.5 * (sigma[2] + sigma[0]);
101 double const tau_m = 0.5 * std::fabs(sigma[2] - sigma[0]);
103 if (sigma_m > t_sigma_max_)
106 f = sigma_m / t_sigma_max_;
108 double const tau_tt =
109 c_ * std::cos(phi_) - t_sigma_max_ * std::sin(phi_);
111 f = std::max(f, tau_m / tau_tt);
116 f = tau_m / (c_ * std::cos(phi_) - sigma_m * std::sin(phi_));
121 const double exp_value = std::exp(b_ * f);
122 for (
auto& k_i : k_data)
124 k_i = std::min(k_i + kr_ * exp_value, k_max_);
130 if (local_coordinate_system_ && (k_data.size() == DisplacementDim))
132 Eigen::Matrix<double, DisplacementDim, DisplacementDim>
const e =
133 local_coordinate_system_->transformation<DisplacementDim>(pos);
134 Eigen::Matrix<double, DisplacementDim, DisplacementDim> k =
135 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
137 for (
int i = 0; i < DisplacementDim; ++i)
139 Eigen::Matrix<double, DisplacementDim, DisplacementDim>
const
140 ei_otimes_ei = e.col(i) * e.col(i).transpose();
142 k += k_data[i] * ei_otimes_ei;
150 template <
int DisplacementDim>
163 "The derivative of the intrinsic permeability k(sigma, ...) with "
164 "respect to stress tensor (sigma) is not implemented because that "
165 "dk/du is normally omitted.");
A failure index dependent permeability model .
PermeabilityMohrCoulombFailureIndexModel(std::string name, ParameterLib::Parameter< double > const &k0, double const kr, double const b, double const c, double const phi, double const k_max, double const t_sigma_max, ParameterLib::CoordinateSystem const *const local_coordinate_system)
void checkScale() const override
PropertyDataType dValue(VariableArray const &variable_array, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const override
double const phi_
Angle of internal friction.
double const t_sigma_max_
Tensile strength parameter.
virtual PropertyDataType value() const
PropertyDataType fromVector(std::vector< double > 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 > > PropertyDataType
template Eigen::Matrix< double, 3, 3 > formEigenTensor< 3 >(MaterialPropertyLib::PropertyDataType const &values)
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray