OGS
Coulomb.cpp
Go to the documentation of this file.
1
10#include "Coulomb.h"
11
12#include <Eigen/LU>
13#include <boost/math/constants/constants.hpp>
14
15#include "BaseLib/Error.h"
16#include "LogPenalty.h"
17#include "MathLib/MathTools.h"
18#include "NumLib/Exceptions.h"
19
20namespace MaterialLib
21{
22namespace Fracture
23{
24namespace Coulomb
25{
27{
28 double Kn = 0.0;
29 double Ks = 0.0;
30 double phi = 0.0; // friction angle
31 double psi = 0.0; // dilation angle
32 double c = 0.0;
33
34 template <typename MaterialProperties>
35 MaterialPropertyValues(MaterialProperties const& mp,
36 double const t,
38 {
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>(); // pi/180
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];
46 }
47};
48
49template <int DisplacementDim>
51 double const t,
53 double const aperture0,
54 Eigen::Ref<Eigen::VectorXd const>
55 /*sigma0*/,
56 Eigen::Ref<Eigen::VectorXd const>
57 w_prev,
58 Eigen::Ref<Eigen::VectorXd const>
59 w,
60 Eigen::Ref<Eigen::VectorXd const>
61 sigma_prev,
62 Eigen::Ref<Eigen::VectorXd>
63 sigma,
64 Eigen::Ref<Eigen::MatrixXd>
65 Kep,
67 material_state_variables)
68{
69 assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
70 &material_state_variables) != nullptr);
71
73 static_cast<StateVariables<DisplacementDim>&>(material_state_variables);
74
75 MaterialPropertyValues const mat(_mp, t, x);
76
77 const int index_ns = DisplacementDim - 1;
78 double const aperture = w[index_ns] + aperture0;
79
80 Eigen::MatrixXd Ke;
81 { // Elastic tangent stiffness
82 Ke = Eigen::MatrixXd::Zero(DisplacementDim, DisplacementDim);
83 for (int i = 0; i < index_ns; i++)
84 {
85 Ke(i, i) = mat.Ks;
86 }
87
88 Ke(index_ns, index_ns) = mat.Kn;
89 }
90
91 // Total plastic aperture compression
92 // NOTE: Initial condition sigma0 seems to be associated with an initial
93 // condition of the w0 = 0. Therefore the initial state is not associated
94 // with a plastic aperture change.
95 { // Exact elastic predictor
96 sigma.noalias() = Ke * (w - w_prev);
97
98 sigma.coeffRef(index_ns) *=
99 logPenaltyDerivative(aperture0, aperture, _penalty_aperture_cutoff);
100 sigma.noalias() += sigma_prev;
101 }
102
103 // correction for an opening fracture
104 if (_tension_cutoff && sigma[DisplacementDim - 1] >= 0)
105 {
106 Kep.setZero();
107 sigma.setZero();
108 state.w_p = w;
109 material_state_variables.setTensileStress(true);
110 return;
111 }
112
113 auto yield_function = [&mat](Eigen::VectorXd const& s)
114 {
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(); // magnitude
118 return mag_tau + sigma_n * std::tan(mat.phi) - mat.c;
119 };
120
121 { // Exit if still in elastic range by checking the shear yield function.
122 double const Fs = yield_function(sigma);
123 material_state_variables.setShearYieldFunctionValue(Fs);
124 if (Fs < .0)
125 {
126 Kep = Ke;
127 Kep(index_ns, index_ns) *= logPenaltyDerivative(
128 aperture0, aperture, _penalty_aperture_cutoff);
129 return;
130 }
131 }
132
133 auto yield_function_derivative = [&mat](Eigen::VectorXd const& s)
134 {
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);
139 return dFs_dS;
140 };
141
142 // plastic potential function: Qs = |tau| + Sn * tan da
143 auto plastic_potential_derivative = [&mat](Eigen::VectorXd const& s)
144 {
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);
149 return dQs_dS;
150 };
151
152 { // Newton
153
154 Eigen::FullPivLU<Eigen::Matrix<double, 1, 1, Eigen::RowMajor>>
155 linear_solver;
156 using ResidualVectorType = Eigen::Matrix<double, 1, 1, Eigen::RowMajor>;
157 using JacobianMatrix = Eigen::Matrix<double, 1, 1, Eigen::RowMajor>;
158
159 JacobianMatrix jacobian;
160 ResidualVectorType solution;
161 solution << 0;
162
163 auto const update_residual = [&](ResidualVectorType& residual)
164 { residual[0] = yield_function(sigma); };
165
166 auto const update_jacobian = [&](JacobianMatrix& jacobian)
167 {
168 jacobian(0, 0) = -yield_function_derivative(sigma).transpose() *
169 Ke * plastic_potential_derivative(sigma);
170 };
171
172 auto const update_solution = [&](ResidualVectorType const& increment)
173 {
174 solution += increment;
175 /*DBUG("analytical = {:g}",
176 Fs / (mat.Ks + mat.Kn * std::tan(mat.psi) * std::tan(mat.phi)))
177 */
178 state.w_p = state.w_p_prev +
179 solution[0] * plastic_potential_derivative(sigma);
180
181 sigma.noalias() = Ke * (w - w_prev - state.w_p + state.w_p_prev);
182
183 sigma.coeffRef(index_ns) *= logPenaltyDerivative(
184 aperture0, aperture, _penalty_aperture_cutoff);
185 sigma.noalias() += sigma_prev;
186 };
187
188 auto newton_solver = NumLib::NewtonRaphson(
189 linear_solver, update_jacobian, update_residual, update_solution,
190 _nonlinear_solver_parameters);
191
192 auto const success_iterations = newton_solver.solve(jacobian);
193
194 if (!success_iterations)
195 {
197 "FractureModel/Coulomb local nonlinear solver didn't "
198 "converge.");
199 }
200
201 // Solution containing lambda is not needed; w_p and sigma already
202 // up to date.
203 }
204
205 { // Update material state shear yield function value.
206 double const Fs = yield_function(sigma);
207 material_state_variables.setShearYieldFunctionValue(Fs);
208 }
209
210 Ke(index_ns, index_ns) *=
211 logPenaltyDerivative(aperture0, aperture, _penalty_aperture_cutoff);
212 Eigen::RowVectorXd const A = yield_function_derivative(sigma).transpose() *
213 Ke /
214 (yield_function_derivative(sigma).transpose() *
215 Ke * plastic_potential_derivative(sigma));
216 Kep = Ke - Ke * plastic_potential_derivative(sigma) * A;
217}
218
219template class Coulomb<2>;
220template class Coulomb<3>;
221
222} // namespace Coulomb
223} // namespace Fracture
224} // namespace MaterialLib
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
Definition Coulomb.cpp:50
double logPenaltyDerivative(double const aperture0, double const aperture, double const aperture_cutoff)
Definition LogPenalty.h:23
MaterialPropertyValues(MaterialProperties const &mp, double const t, ParameterLib::SpatialPosition const &x)
Definition Coulomb.cpp:35
Eigen::Matrix< double, DisplacementDim, 1 > w_p
Definition Coulomb.h:35
Eigen::Matrix< double, DisplacementDim, 1 > w_p_prev
Definition Coulomb.h:39