Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "SolidMechanicsPlasticMohrCoulomb.h"
11 : #include "RankFourTensor.h"
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("SolidMechanicsApp", SolidMechanicsPlasticMohrCoulomb);
15 : registerMooseObjectRenamed("SolidMechanicsApp",
16 : TensorMechanicsPlasticMohrCoulomb,
17 : "01/01/2025 00:00",
18 : SolidMechanicsPlasticMohrCoulomb);
19 :
20 : InputParameters
21 152 : SolidMechanicsPlasticMohrCoulomb::validParams()
22 : {
23 152 : InputParameters params = SolidMechanicsPlasticModel::validParams();
24 304 : params.addRequiredParam<UserObjectName>(
25 : "cohesion",
26 : "A SolidMechanicsHardening UserObject that defines hardening of the cohesion. "
27 : "Physically the cohesion should not be negative.");
28 304 : params.addRequiredParam<UserObjectName>(
29 : "friction_angle",
30 : "A SolidMechanicsHardening UserObject that defines hardening of the "
31 : "friction angle (in radians). Physically the friction angle should be "
32 : "between 0 and 90deg.");
33 304 : params.addRequiredParam<UserObjectName>(
34 : "dilation_angle",
35 : "A SolidMechanicsHardening UserObject that defines hardening of the "
36 : "dilation angle (in radians). Usually the dilation angle is not greater "
37 : "than the friction angle, and it is between 0 and 90deg.");
38 456 : params.addRangeCheckedParam<Real>(
39 : "mc_edge_smoother",
40 304 : 25.0,
41 : "mc_edge_smoother>=0 & mc_edge_smoother<=30",
42 : "Smoothing parameter: the edges of the cone are smoothed by the given amount.");
43 304 : MooseEnum tip_scheme("hyperbolic cap", "hyperbolic");
44 304 : params.addParam<MooseEnum>(
45 : "tip_scheme", tip_scheme, "Scheme by which the pyramid's tip will be smoothed.");
46 304 : params.addRequiredRangeCheckedParam<Real>("mc_tip_smoother",
47 : "mc_tip_smoother>=0",
48 : "Smoothing parameter: the cone vertex at mean = "
49 : "cohesion*cot(friction_angle), will be smoothed by "
50 : "the given amount. Typical value is 0.1*cohesion");
51 304 : params.addParam<Real>(
52 : "cap_start",
53 304 : 0.0,
54 : "For the 'cap' tip_scheme, smoothing is performed in the stress_mean > cap_start region");
55 456 : params.addRangeCheckedParam<Real>("cap_rate",
56 304 : 0.0,
57 : "cap_rate>=0",
58 : "For the 'cap' tip_scheme, this controls how quickly the cap "
59 : "degenerates to a hemisphere: small values mean a slow "
60 : "degeneration to a hemisphere (and zero means the 'cap' will "
61 : "be totally inactive). Typical value is 1/tensile_strength");
62 304 : params.addParam<Real>("mc_lode_cutoff",
63 : "If the second invariant of stress is less than this "
64 : "amount, the Lode angle is assumed to be zero. This is "
65 : "to gaurd against precision-loss problems, and this "
66 : "parameter should be set small. Default = "
67 : "0.00001*((yield_Function_tolerance)^2)");
68 152 : params.addClassDescription("Non-associative Mohr-Coulomb plasticity with hardening/softening");
69 :
70 152 : return params;
71 152 : }
72 :
73 78 : SolidMechanicsPlasticMohrCoulomb::SolidMechanicsPlasticMohrCoulomb(
74 78 : const InputParameters & parameters)
75 : : SolidMechanicsPlasticModel(parameters),
76 78 : _cohesion(getUserObject<SolidMechanicsHardeningModel>("cohesion")),
77 78 : _phi(getUserObject<SolidMechanicsHardeningModel>("friction_angle")),
78 78 : _psi(getUserObject<SolidMechanicsHardeningModel>("dilation_angle")),
79 156 : _tip_scheme(getParam<MooseEnum>("tip_scheme")),
80 156 : _small_smoother2(Utility::pow<2>(getParam<Real>("mc_tip_smoother"))),
81 156 : _cap_start(getParam<Real>("cap_start")),
82 156 : _cap_rate(getParam<Real>("cap_rate")),
83 156 : _tt(getParam<Real>("mc_edge_smoother") * libMesh::pi / 180.0),
84 78 : _costt(std::cos(_tt)),
85 78 : _sintt(std::sin(_tt)),
86 78 : _cos3tt(std::cos(3 * _tt)),
87 78 : _sin3tt(std::sin(3 * _tt)),
88 78 : _cos6tt(std::cos(6 * _tt)),
89 78 : _sin6tt(std::sin(6 * _tt)),
90 80 : _lode_cutoff(parameters.isParamValid("mc_lode_cutoff") ? getParam<Real>("mc_lode_cutoff")
91 78 : : 1.0E-5 * Utility::pow<2>(_f_tol))
92 :
93 : {
94 78 : if (_lode_cutoff < 0)
95 1 : mooseError("mc_lode_cutoff must not be negative");
96 :
97 : // With arbitary UserObjects, it is impossible to check everything, and
98 : // I think this is the best I can do
99 77 : if (phi(0) < 0 || psi(0) < 0 || phi(0) > libMesh::pi / 2.0 || psi(0) > libMesh::pi / 2.0)
100 0 : mooseError("Mohr-Coulomb friction and dilation angles must lie in [0, Pi/2]");
101 77 : if (phi(0) < psi(0))
102 1 : mooseError("Mohr-Coulomb friction angle must not be less than Mohr-Coulomb dilation angle");
103 76 : if (cohesion(0) < 0)
104 0 : mooseError("Mohr-Coulomb cohesion must not be negative");
105 :
106 : // check Abbo et al's convexity constraint (Eqn c.18 in their paper)
107 : // With an arbitrary UserObject, it is impossible to check for all angles
108 : // I think the following is the best we can do
109 76 : Real sin_angle = std::sin(std::max(phi(0), psi(0)));
110 76 : sin_angle = std::max(sin_angle, std::sin(std::max(phi(1E6), psi(1E6))));
111 76 : Real rhs = std::sqrt(3) * (35 * std::sin(_tt) + 14 * std::sin(5 * _tt) - 5 * std::sin(7 * _tt)) /
112 76 : 16 / Utility::pow<5>(std::cos(_tt)) / (11 - 10 * std::cos(2 * _tt));
113 76 : if (rhs <= sin_angle)
114 2 : mooseError("Mohr-Coulomb edge smoothing angle is too small and a non-convex yield surface will "
115 : "result. Please choose a larger value");
116 74 : }
117 :
118 : Real
119 791688 : SolidMechanicsPlasticMohrCoulomb::yieldFunction(const RankTwoTensor & stress, Real intnl) const
120 : {
121 791688 : Real mean_stress = stress.trace() / 3.0;
122 791688 : Real sinphi = std::sin(phi(intnl));
123 791688 : Real cosphi = std::cos(phi(intnl));
124 791688 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
125 791688 : if (std::abs(sin3Lode) <= _sin3tt)
126 : {
127 : // the non-edge-smoothed version
128 : std::vector<Real> eigvals;
129 211800 : stress.symmetricEigenvalues(eigvals);
130 423600 : return mean_stress * sinphi +
131 211800 : std::sqrt(smooth(stress) +
132 211800 : 0.25 * Utility::pow<2>(eigvals[2] - eigvals[0] +
133 211800 : (eigvals[2] + eigvals[0] - 2 * mean_stress) * sinphi)) -
134 211800 : cohesion(intnl) * cosphi;
135 211800 : }
136 : else
137 : {
138 : // the edge-smoothed version
139 : Real aaa, bbb, ccc;
140 579888 : abbo(sin3Lode, sinphi, aaa, bbb, ccc);
141 579888 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
142 579888 : Real sibar2 = stress.secondInvariant();
143 579888 : return mean_stress * sinphi + std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk)) -
144 579888 : cohesion(intnl) * cosphi;
145 : }
146 : }
147 :
148 : RankTwoTensor
149 1308816 : SolidMechanicsPlasticMohrCoulomb::df_dsig(const RankTwoTensor & stress, const Real sin_angle) const
150 : {
151 1308816 : Real mean_stress = stress.trace() / 3.0;
152 1308816 : RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
153 1308816 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
154 1308816 : if (std::abs(sin3Lode) <= _sin3tt)
155 : {
156 : // the non-edge-smoothed version
157 : std::vector<Real> eigvals;
158 : std::vector<RankTwoTensor> deigvals;
159 315956 : stress.dsymmetricEigenvalues(eigvals, deigvals);
160 315956 : Real tmp = eigvals[2] - eigvals[0] + (eigvals[2] + eigvals[0] - 2.0 * mean_stress) * sin_angle;
161 : RankTwoTensor dtmp =
162 315956 : deigvals[2] - deigvals[0] + (deigvals[2] + deigvals[0] - 2.0 * dmean_stress) * sin_angle;
163 315956 : Real denom = std::sqrt(smooth(stress) + 0.25 * Utility::pow<2>(tmp));
164 315956 : return dmean_stress * sin_angle +
165 315956 : (0.5 * dsmooth(stress) * dmean_stress + 0.25 * tmp * dtmp) / denom;
166 315956 : }
167 : else
168 : {
169 : // the edge-smoothed version
170 : Real aaa, bbb, ccc;
171 992860 : abbo(sin3Lode, sin_angle, aaa, bbb, ccc);
172 992860 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
173 992860 : RankTwoTensor dkk = (bbb + 2 * ccc * sin3Lode) * stress.dsin3Lode(_lode_cutoff);
174 992860 : Real sibar2 = stress.secondInvariant();
175 992860 : RankTwoTensor dsibar2 = stress.dsecondInvariant();
176 992860 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
177 992860 : return dmean_stress * sin_angle + (0.5 * dsmooth(stress) * dmean_stress +
178 992860 : 0.5 * dsibar2 * Utility::pow<2>(kk) + sibar2 * kk * dkk) /
179 992860 : denom;
180 : }
181 : }
182 :
183 : RankTwoTensor
184 306792 : SolidMechanicsPlasticMohrCoulomb::dyieldFunction_dstress(const RankTwoTensor & stress,
185 : Real intnl) const
186 : {
187 306792 : Real sinphi = std::sin(phi(intnl));
188 306792 : return df_dsig(stress, sinphi);
189 : }
190 :
191 : Real
192 306792 : SolidMechanicsPlasticMohrCoulomb::dyieldFunction_dintnl(const RankTwoTensor & stress,
193 : Real intnl) const
194 : {
195 306792 : Real sin_angle = std::sin(phi(intnl));
196 306792 : Real cos_angle = std::cos(phi(intnl));
197 306792 : Real dsin_angle = cos_angle * dphi(intnl);
198 306792 : Real dcos_angle = -sin_angle * dphi(intnl);
199 :
200 306792 : Real mean_stress = stress.trace() / 3.0;
201 306792 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
202 306792 : if (std::abs(sin3Lode) <= _sin3tt)
203 : {
204 : // the non-edge-smoothed version
205 : std::vector<Real> eigvals;
206 72812 : stress.symmetricEigenvalues(eigvals);
207 72812 : Real tmp = eigvals[2] - eigvals[0] + (eigvals[2] + eigvals[0] - 2 * mean_stress) * sin_angle;
208 72812 : Real dtmp = (eigvals[2] + eigvals[0] - 2 * mean_stress) * dsin_angle;
209 72812 : Real denom = std::sqrt(smooth(stress) + 0.25 * Utility::pow<2>(tmp));
210 72812 : return mean_stress * dsin_angle + 0.25 * tmp * dtmp / denom - dcohesion(intnl) * cos_angle -
211 72812 : cohesion(intnl) * dcos_angle;
212 72812 : }
213 : else
214 : {
215 : // the edge-smoothed version
216 : Real aaa, bbb, ccc;
217 233980 : abbo(sin3Lode, sin_angle, aaa, bbb, ccc);
218 : Real daaa, dbbb, dccc;
219 233980 : dabbo(sin3Lode, sin_angle, daaa, dbbb, dccc);
220 233980 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
221 233980 : Real dkk = (daaa + dbbb * sin3Lode + dccc * Utility::pow<2>(sin3Lode)) * dsin_angle;
222 233980 : Real sibar2 = stress.secondInvariant();
223 233980 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
224 233980 : return mean_stress * dsin_angle + sibar2 * kk * dkk / denom - dcohesion(intnl) * cos_angle -
225 233980 : cohesion(intnl) * dcos_angle;
226 : }
227 : }
228 :
229 : RankTwoTensor
230 1002024 : SolidMechanicsPlasticMohrCoulomb::flowPotential(const RankTwoTensor & stress, Real intnl) const
231 : {
232 1002024 : Real sinpsi = std::sin(psi(intnl));
233 1002024 : return df_dsig(stress, sinpsi);
234 : }
235 :
236 : RankFourTensor
237 306792 : SolidMechanicsPlasticMohrCoulomb::dflowPotential_dstress(const RankTwoTensor & stress,
238 : Real intnl) const
239 : {
240 306792 : RankFourTensor dr_dstress;
241 306792 : Real sin_angle = std::sin(psi(intnl));
242 306792 : Real mean_stress = stress.trace() / 3.0;
243 306792 : RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
244 306792 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
245 306792 : if (std::abs(sin3Lode) <= _sin3tt)
246 : {
247 : // the non-edge-smoothed version
248 : std::vector<Real> eigvals;
249 : std::vector<RankTwoTensor> deigvals;
250 : std::vector<RankFourTensor> d2eigvals;
251 72812 : stress.dsymmetricEigenvalues(eigvals, deigvals);
252 72812 : stress.d2symmetricEigenvalues(d2eigvals);
253 :
254 72812 : Real tmp = eigvals[2] - eigvals[0] + (eigvals[2] + eigvals[0] - 2.0 * mean_stress) * sin_angle;
255 : RankTwoTensor dtmp =
256 72812 : deigvals[2] - deigvals[0] + (deigvals[2] + deigvals[0] - 2.0 * dmean_stress) * sin_angle;
257 72812 : Real denom = std::sqrt(smooth(stress) + 0.25 * Utility::pow<2>(tmp));
258 : Real denom3 = Utility::pow<3>(denom);
259 72812 : Real d2smooth_over_denom = d2smooth(stress) / denom;
260 72812 : RankTwoTensor numer = dsmooth(stress) * dmean_stress + 0.5 * tmp * dtmp;
261 :
262 72812 : dr_dstress = 0.25 * tmp *
263 145624 : (d2eigvals[2] - d2eigvals[0] + (d2eigvals[2] + d2eigvals[0]) * sin_angle) / denom;
264 :
265 291248 : for (unsigned i = 0; i < 3; ++i)
266 873744 : for (unsigned j = 0; j < 3; ++j)
267 2621232 : for (unsigned k = 0; k < 3; ++k)
268 7863696 : for (unsigned l = 0; l < 3; ++l)
269 : {
270 5897772 : dr_dstress(i, j, k, l) +=
271 5897772 : 0.5 * d2smooth_over_denom * dmean_stress(i, j) * dmean_stress(k, l);
272 5897772 : dr_dstress(i, j, k, l) += 0.25 * dtmp(i, j) * dtmp(k, l) / denom;
273 5897772 : dr_dstress(i, j, k, l) -= 0.25 * numer(i, j) * numer(k, l) / denom3;
274 : }
275 72812 : }
276 : else
277 : {
278 : // the edge-smoothed version
279 : Real aaa, bbb, ccc;
280 233980 : abbo(sin3Lode, sin_angle, aaa, bbb, ccc);
281 233980 : RankTwoTensor dsin3Lode = stress.dsin3Lode(_lode_cutoff);
282 233980 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
283 233980 : RankTwoTensor dkk = (bbb + 2 * ccc * sin3Lode) * dsin3Lode;
284 233980 : RankFourTensor d2kk = (bbb + 2 * ccc * sin3Lode) * stress.d2sin3Lode(_lode_cutoff);
285 935920 : for (unsigned i = 0; i < 3; ++i)
286 2807760 : for (unsigned j = 0; j < 3; ++j)
287 8423280 : for (unsigned k = 0; k < 3; ++k)
288 25269840 : for (unsigned l = 0; l < 3; ++l)
289 18952380 : d2kk(i, j, k, l) += 2 * ccc * dsin3Lode(i, j) * dsin3Lode(k, l);
290 :
291 233980 : Real sibar2 = stress.secondInvariant();
292 233980 : RankTwoTensor dsibar2 = stress.dsecondInvariant();
293 233980 : RankFourTensor d2sibar2 = stress.d2secondInvariant();
294 :
295 233980 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
296 : Real denom3 = Utility::pow<3>(denom);
297 233980 : Real d2smooth_over_denom = d2smooth(stress) / denom;
298 : RankTwoTensor numer_full =
299 233980 : 0.5 * dsmooth(stress) * dmean_stress + 0.5 * dsibar2 * kk * kk + sibar2 * kk * dkk;
300 :
301 233980 : dr_dstress = (0.5 * d2sibar2 * Utility::pow<2>(kk) + sibar2 * kk * d2kk) / denom;
302 935920 : for (unsigned i = 0; i < 3; ++i)
303 2807760 : for (unsigned j = 0; j < 3; ++j)
304 8423280 : for (unsigned k = 0; k < 3; ++k)
305 25269840 : for (unsigned l = 0; l < 3; ++l)
306 : {
307 18952380 : dr_dstress(i, j, k, l) +=
308 18952380 : 0.5 * d2smooth_over_denom * dmean_stress(i, j) * dmean_stress(k, l);
309 18952380 : dr_dstress(i, j, k, l) +=
310 18952380 : (dsibar2(i, j) * dkk(k, l) * kk + dkk(i, j) * dsibar2(k, l) * kk +
311 18952380 : sibar2 * dkk(i, j) * dkk(k, l)) /
312 : denom;
313 18952380 : dr_dstress(i, j, k, l) -= numer_full(i, j) * numer_full(k, l) / denom3;
314 : }
315 : }
316 306792 : return dr_dstress;
317 : }
318 :
319 : RankTwoTensor
320 306792 : SolidMechanicsPlasticMohrCoulomb::dflowPotential_dintnl(const RankTwoTensor & stress,
321 : Real intnl) const
322 : {
323 306792 : Real sin_angle = std::sin(psi(intnl));
324 306792 : Real dsin_angle = std::cos(psi(intnl)) * dpsi(intnl);
325 :
326 306792 : Real mean_stress = stress.trace() / 3.0;
327 306792 : RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
328 306792 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
329 :
330 306792 : if (std::abs(sin3Lode) <= _sin3tt)
331 : {
332 : // the non-edge-smoothed version
333 : std::vector<Real> eigvals;
334 : std::vector<RankTwoTensor> deigvals;
335 72812 : stress.dsymmetricEigenvalues(eigvals, deigvals);
336 72812 : Real tmp = eigvals[2] - eigvals[0] + (eigvals[2] + eigvals[0] - 2.0 * mean_stress) * sin_angle;
337 72812 : Real dtmp_dintnl = (eigvals[2] + eigvals[0] - 2 * mean_stress) * dsin_angle;
338 : RankTwoTensor dtmp_dstress =
339 72812 : deigvals[2] - deigvals[0] + (deigvals[2] + deigvals[0] - 2.0 * dmean_stress) * sin_angle;
340 : RankTwoTensor d2tmp_dstress_dintnl =
341 72812 : (deigvals[2] + deigvals[0] - 2.0 * dmean_stress) * dsin_angle;
342 72812 : Real denom = std::sqrt(smooth(stress) + 0.25 * Utility::pow<2>(tmp));
343 72812 : return dmean_stress * dsin_angle + 0.25 * dtmp_dintnl * dtmp_dstress / denom +
344 72812 : 0.25 * tmp * d2tmp_dstress_dintnl / denom -
345 72812 : 0.5 * (dsmooth(stress) * dmean_stress + 0.5 * tmp * dtmp_dstress) * 0.25 * tmp *
346 72812 : dtmp_dintnl / Utility::pow<3>(denom);
347 72812 : }
348 : else
349 : {
350 : // the edge-smoothed version
351 : Real aaa, bbb, ccc;
352 233980 : abbo(sin3Lode, sin_angle, aaa, bbb, ccc);
353 233980 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
354 :
355 : Real daaa, dbbb, dccc;
356 233980 : dabbo(sin3Lode, sin_angle, daaa, dbbb, dccc);
357 233980 : Real dkk_dintnl = (daaa + dbbb * sin3Lode + dccc * Utility::pow<2>(sin3Lode)) * dsin_angle;
358 233980 : RankTwoTensor dkk_dstress = (bbb + 2 * ccc * sin3Lode) * stress.dsin3Lode(_lode_cutoff);
359 : RankTwoTensor d2kk_dstress_dintnl =
360 233980 : (dbbb + 2 * dccc * sin3Lode) * stress.dsin3Lode(_lode_cutoff) * dsin_angle;
361 :
362 233980 : Real sibar2 = stress.secondInvariant();
363 233980 : RankTwoTensor dsibar2 = stress.dsecondInvariant();
364 233980 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
365 :
366 233980 : return dmean_stress * dsin_angle +
367 233980 : (dsibar2 * kk * dkk_dintnl + sibar2 * dkk_dintnl * dkk_dstress +
368 233980 : sibar2 * kk * d2kk_dstress_dintnl) /
369 : denom -
370 233980 : (0.5 * dsmooth(stress) * dmean_stress + 0.5 * dsibar2 * Utility::pow<2>(kk) +
371 467960 : sibar2 * kk * dkk_dstress) *
372 233980 : sibar2 * kk * dkk_dintnl / Utility::pow<3>(denom);
373 : }
374 : }
375 :
376 : Real
377 1098556 : SolidMechanicsPlasticMohrCoulomb::cohesion(const Real internal_param) const
378 : {
379 1098556 : return _cohesion.value(internal_param);
380 : }
381 :
382 : Real
383 306792 : SolidMechanicsPlasticMohrCoulomb::dcohesion(const Real internal_param) const
384 : {
385 306792 : return _cohesion.derivative(internal_param);
386 : }
387 :
388 : Real
389 2504135 : SolidMechanicsPlasticMohrCoulomb::phi(const Real internal_param) const
390 : {
391 2504135 : return _phi.value(internal_param);
392 : }
393 :
394 : Real
395 613584 : SolidMechanicsPlasticMohrCoulomb::dphi(const Real internal_param) const
396 : {
397 613584 : return _phi.derivative(internal_param);
398 : }
399 :
400 : Real
401 1922783 : SolidMechanicsPlasticMohrCoulomb::psi(const Real internal_param) const
402 : {
403 1922783 : return _psi.value(internal_param);
404 : }
405 :
406 : Real
407 306792 : SolidMechanicsPlasticMohrCoulomb::dpsi(const Real internal_param) const
408 : {
409 306792 : return _psi.derivative(internal_param);
410 : }
411 :
412 : void
413 2274688 : SolidMechanicsPlasticMohrCoulomb::abbo(
414 : const Real sin3lode, const Real sin_angle, Real & aaa, Real & bbb, Real & ccc) const
415 : {
416 2274688 : Real tmp1 = (sin3lode >= 0 ? _costt - sin_angle * _sintt / std::sqrt(3.0)
417 891620 : : _costt + sin_angle * _sintt / std::sqrt(3.0));
418 2274688 : Real tmp2 = (sin3lode >= 0 ? _sintt + sin_angle * _costt / std::sqrt(3.0)
419 891620 : : -_sintt + sin_angle * _costt / std::sqrt(3.0));
420 :
421 2274688 : ccc = -_cos3tt * tmp1;
422 2274688 : ccc += (sin3lode >= 0 ? -3 * _sin3tt * tmp2 : 3 * _sin3tt * tmp2);
423 2274688 : ccc /= 18 * Utility::pow<3>(_cos3tt);
424 :
425 2274688 : bbb = (sin3lode >= 0 ? _sin6tt * tmp1 : -_sin6tt * tmp1);
426 2274688 : bbb -= 6 * _cos6tt * tmp2;
427 2274688 : bbb /= 18 * Utility::pow<3>(_cos3tt);
428 :
429 2274688 : aaa = (sin3lode >= 0 ? -sin_angle * _sintt / std::sqrt(3.0) - bbb * _sin3tt
430 891620 : : sin_angle * _sintt / std::sqrt(3.0) + bbb * _sin3tt);
431 2274688 : aaa += -ccc * Utility::pow<2>(_sin3tt) + _costt;
432 2274688 : }
433 :
434 : void
435 467960 : SolidMechanicsPlasticMohrCoulomb::dabbo(
436 : const Real sin3lode, const Real /*sin_angle*/, Real & daaa, Real & dbbb, Real & dccc) const
437 : {
438 467960 : Real dtmp1 = (sin3lode >= 0 ? -_sintt / std::sqrt(3.0) : _sintt / std::sqrt(3.0));
439 467960 : Real dtmp2 = _costt / std::sqrt(3.0);
440 :
441 467960 : dccc = -_cos3tt * dtmp1;
442 467960 : dccc += (sin3lode >= 0 ? -3 * _sin3tt * dtmp2 : 3 * _sin3tt * dtmp2);
443 467960 : dccc /= 18 * Utility::pow<3>(_cos3tt);
444 :
445 467960 : dbbb = (sin3lode >= 0 ? _sin6tt * dtmp1 : -_sin6tt * dtmp1);
446 467960 : dbbb -= 6 * _cos6tt * dtmp2;
447 467960 : dbbb /= 18 * Utility::pow<3>(_cos3tt);
448 :
449 467960 : daaa = (sin3lode >= 0 ? -_sintt / std::sqrt(3.0) - dbbb * _sin3tt
450 180344 : : _sintt / std::sqrt(3.0) + dbbb * _sin3tt);
451 467960 : daaa += -dccc * Utility::pow<2>(_sin3tt);
452 467960 : }
453 :
454 : Real
455 3020880 : SolidMechanicsPlasticMohrCoulomb::smooth(const RankTwoTensor & stress) const
456 : {
457 3020880 : Real smoother2 = _small_smoother2;
458 3020880 : if (_tip_scheme == "cap")
459 : {
460 368160 : Real x = stress.trace() / 3.0 - _cap_start;
461 : Real p = 0;
462 368160 : if (x > 0)
463 330288 : p = x * (1 - std::exp(-_cap_rate * x));
464 368160 : smoother2 += Utility::pow<2>(p);
465 : }
466 3020880 : return smoother2;
467 : }
468 :
469 : Real
470 1922400 : SolidMechanicsPlasticMohrCoulomb::dsmooth(const RankTwoTensor & stress) const
471 : {
472 : Real dsmoother2 = 0;
473 1922400 : if (_tip_scheme == "cap")
474 : {
475 236088 : Real x = stress.trace() / 3.0 - _cap_start;
476 : Real p = 0;
477 : Real dp_dx = 0;
478 236088 : if (x > 0)
479 : {
480 211392 : p = x * (1 - std::exp(-_cap_rate * x));
481 211392 : dp_dx = (1 - std::exp(-_cap_rate * x)) + x * _cap_rate * std::exp(-_cap_rate * x);
482 : }
483 236088 : dsmoother2 += 2 * p * dp_dx;
484 : }
485 1922400 : return dsmoother2;
486 : }
487 :
488 : Real
489 306792 : SolidMechanicsPlasticMohrCoulomb::d2smooth(const RankTwoTensor & stress) const
490 : {
491 : Real d2smoother2 = 0;
492 306792 : if (_tip_scheme == "cap")
493 : {
494 37776 : Real x = stress.trace() / 3.0 - _cap_start;
495 : Real p = 0;
496 : Real dp_dx = 0;
497 : Real d2p_dx2 = 0;
498 37776 : if (x > 0)
499 : {
500 33360 : p = x * (1 - std::exp(-_cap_rate * x));
501 33360 : dp_dx = (1 - std::exp(-_cap_rate * x)) + x * _cap_rate * std::exp(-_cap_rate * x);
502 33360 : d2p_dx2 = 2 * _cap_rate * std::exp(-_cap_rate * x) -
503 33360 : x * Utility::pow<2>(_cap_rate) * std::exp(-_cap_rate * x);
504 : }
505 37776 : d2smoother2 += 2 * Utility::pow<2>(dp_dx) + 2 * p * d2p_dx2;
506 : }
507 306792 : return d2smoother2;
508 : }
509 :
510 : std::string
511 0 : SolidMechanicsPlasticMohrCoulomb::modelName() const
512 : {
513 0 : return "MohrCoulomb";
514 : }
|
