13#include <boost/math/constants/constants.hpp>
34 template <
typename MaterialProperties>
39 Kn = mp.normal_stiffness(t, x)[0];
40 Ks = mp.shear_stiffness(t, x)[0];
41 auto constexpr degree =
42 boost::math::constants::degree<double>();
43 phi = mp.friction_angle(t, x)[0] * degree;
44 psi = mp.dilatancy_angle(t, x)[0] * degree;
45 c = mp.cohesion(t, x)[0];
49template <
int DisplacementDim>
53 double const aperture0,
54 Eigen::Ref<Eigen::VectorXd const>
56 Eigen::Ref<Eigen::VectorXd const>
58 Eigen::Ref<Eigen::VectorXd const>
60 Eigen::Ref<Eigen::VectorXd const>
62 Eigen::Ref<Eigen::VectorXd>
64 Eigen::Ref<Eigen::MatrixXd>
67 material_state_variables)
70 &material_state_variables) !=
nullptr);
77 const int index_ns = DisplacementDim - 1;
78 double const aperture = w[index_ns] + aperture0;
82 Ke = Eigen::MatrixXd::Zero(DisplacementDim, DisplacementDim);
83 for (
int i = 0; i < index_ns; i++)
88 Ke(index_ns, index_ns) = mat.Kn;
96 sigma.noalias() = Ke * (w - w_prev);
98 sigma.coeffRef(index_ns) *=
100 sigma.noalias() += sigma_prev;
104 if (_tension_cutoff && sigma[DisplacementDim - 1] >= 0)
113 auto yield_function = [&mat](Eigen::VectorXd
const& s)
115 double const sigma_n = s[DisplacementDim - 1];
116 Eigen::VectorXd
const sigma_s = s.head(DisplacementDim - 1);
117 double const mag_tau = sigma_s.norm();
118 return mag_tau + sigma_n * std::tan(mat.phi) - mat.c;
122 double const Fs = yield_function(sigma);
128 aperture0, aperture, _penalty_aperture_cutoff);
133 auto yield_function_derivative = [&mat](Eigen::VectorXd
const& s)
135 Eigen::Matrix<double, DisplacementDim, 1> dFs_dS;
136 dFs_dS.template head<DisplacementDim - 1>().noalias() =
137 s.template head<DisplacementDim - 1>().normalized();
138 dFs_dS.coeffRef(DisplacementDim - 1) = std::tan(mat.phi);
143 auto plastic_potential_derivative = [&mat](Eigen::VectorXd
const& s)
145 Eigen::Matrix<double, DisplacementDim, 1> dQs_dS;
146 dQs_dS.template head<DisplacementDim - 1>().noalias() =
147 s.template head<DisplacementDim - 1>().normalized();
148 dQs_dS.coeffRef(DisplacementDim - 1) = std::tan(mat.psi);
154 Eigen::FullPivLU<Eigen::Matrix<double, 1, 1, Eigen::RowMajor>>
156 using ResidualVectorType = Eigen::Matrix<double, 1, 1, Eigen::RowMajor>;
157 using JacobianMatrix = Eigen::Matrix<double, 1, 1, Eigen::RowMajor>;
159 JacobianMatrix jacobian;
160 ResidualVectorType solution;
163 auto const update_residual = [&](ResidualVectorType& residual)
164 { residual[0] = yield_function(sigma); };
166 auto const update_jacobian = [&](JacobianMatrix& jacobian)
168 jacobian(0, 0) = -yield_function_derivative(sigma).transpose() *
169 Ke * plastic_potential_derivative(sigma);
172 auto const update_solution = [&](ResidualVectorType
const& increment)
174 solution += increment;
179 solution[0] * plastic_potential_derivative(sigma);
181 sigma.noalias() = Ke * (w - w_prev - state.
w_p + state.
w_p_prev);
184 aperture0, aperture, _penalty_aperture_cutoff);
185 sigma.noalias() += sigma_prev;
190 decltype(update_jacobian), ResidualVectorType,
191 decltype(update_residual),
192 decltype(update_solution)>(
193 linear_solver, update_jacobian, update_residual,
194 update_solution, _nonlinear_solver_parameters);
196 auto const success_iterations = newton_solver.
solve(jacobian);
198 if (!success_iterations)
201 "FractureModel/Coulomb local nonlinear solver didn't "
210 double const Fs = yield_function(sigma);
214 Ke(index_ns, index_ns) *=
216 Eigen::RowVectorXd
const A = yield_function_derivative(sigma).transpose() *
218 (yield_function_derivative(sigma).transpose() *
219 Ke * plastic_potential_derivative(sigma));
220 Kep = Ke - Ke * plastic_potential_derivative(sigma) * A;
void computeConstitutiveRelation(double const t, ParameterLib::SpatialPosition const &x, double const aperture0, Eigen::Ref< Eigen::VectorXd const > sigma0, Eigen::Ref< Eigen::VectorXd const > w_prev, Eigen::Ref< Eigen::VectorXd const > w, Eigen::Ref< Eigen::VectorXd const > sigma_prev, Eigen::Ref< Eigen::VectorXd > sigma, Eigen::Ref< Eigen::MatrixXd > Kep, typename FractureModelBase< DisplacementDim >::MaterialStateVariables &material_state_variables) override
std::optional< int > solve(JacobianMatrix &jacobian) const
double logPenaltyDerivative(double const aperture0, double const aperture, double const aperture_cutoff)
MaterialPropertyValues(MaterialProperties const &mp, double const t, ParameterLib::SpatialPosition const &x)
Eigen::Matrix< double, DisplacementDim, 1 > w_p
Eigen::Matrix< double, DisplacementDim, 1 > w_p_prev
void setShearYieldFunctionValue(double Fs)
void setTensileStress(bool flag)