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 240 : SolidMechanicsPlasticMohrCoulomb::validParams()
22 : {
23 240 : InputParameters params = SolidMechanicsPlasticModel::validParams();
24 480 : params.addRequiredParam<UserObjectName>(
25 : "cohesion",
26 : "A SolidMechanicsHardening UserObject that defines hardening of the cohesion. "
27 : "Physically the cohesion should not be negative.");
28 480 : 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 480 : 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 720 : params.addRangeCheckedParam<Real>(
39 : "mc_edge_smoother",
40 480 : 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 480 : MooseEnum tip_scheme("hyperbolic cap", "hyperbolic");
44 480 : params.addParam<MooseEnum>(
45 : "tip_scheme", tip_scheme, "Scheme by which the pyramid's tip will be smoothed.");
46 480 : 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 480 : params.addParam<Real>(
52 : "cap_start",
53 480 : 0.0,
54 : "For the 'cap' tip_scheme, smoothing is performed in the stress_mean > cap_start region");
55 720 : params.addRangeCheckedParam<Real>("cap_rate",
56 480 : 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 480 : 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 240 : params.addClassDescription("Non-associative Mohr-Coulomb plasticity with hardening/softening");
69 :
70 240 : return params;
71 240 : }
72 :
73 124 : SolidMechanicsPlasticMohrCoulomb::SolidMechanicsPlasticMohrCoulomb(
74 124 : const InputParameters & parameters)
75 : : SolidMechanicsPlasticModel(parameters),
76 124 : _cohesion(getUserObject<SolidMechanicsHardeningModel>("cohesion")),
77 124 : _phi(getUserObject<SolidMechanicsHardeningModel>("friction_angle")),
78 124 : _psi(getUserObject<SolidMechanicsHardeningModel>("dilation_angle")),
79 248 : _tip_scheme(getParam<MooseEnum>("tip_scheme")),
80 248 : _small_smoother2(Utility::pow<2>(getParam<Real>("mc_tip_smoother"))),
81 248 : _cap_start(getParam<Real>("cap_start")),
82 248 : _cap_rate(getParam<Real>("cap_rate")),
83 248 : _tt(getParam<Real>("mc_edge_smoother") * libMesh::pi / 180.0),
84 124 : _costt(std::cos(_tt)),
85 124 : _sintt(std::sin(_tt)),
86 124 : _cos3tt(std::cos(3 * _tt)),
87 124 : _sin3tt(std::sin(3 * _tt)),
88 124 : _cos6tt(std::cos(6 * _tt)),
89 124 : _sin6tt(std::sin(6 * _tt)),
90 128 : _lode_cutoff(parameters.isParamValid("mc_lode_cutoff") ? getParam<Real>("mc_lode_cutoff")
91 124 : : 1.0E-5 * Utility::pow<2>(_f_tol))
92 :
93 : {
94 124 : if (_lode_cutoff < 0)
95 2 : 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 122 : 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 122 : if (phi(0) < psi(0))
102 2 : mooseError("Mohr-Coulomb friction angle must not be less than Mohr-Coulomb dilation angle");
103 120 : 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 120 : Real sin_angle = std::sin(std::max(phi(0), psi(0)));
110 120 : sin_angle = std::max(sin_angle, std::sin(std::max(phi(1E6), psi(1E6))));
111 120 : Real rhs = std::sqrt(3) * (35 * std::sin(_tt) + 14 * std::sin(5 * _tt) - 5 * std::sin(7 * _tt)) /
112 120 : 16 / Utility::pow<5>(std::cos(_tt)) / (11 - 10 * std::cos(2 * _tt));
113 120 : if (rhs <= sin_angle)
114 4 : mooseError("Mohr-Coulomb edge smoothing angle is too small and a non-convex yield surface will "
115 : "result. Please choose a larger value");
116 116 : }
117 :
118 : Real
119 991968 : SolidMechanicsPlasticMohrCoulomb::yieldFunction(const RankTwoTensor & stress, Real intnl) const
120 : {
121 991968 : Real mean_stress = stress.trace() / 3.0;
122 991968 : Real sinphi = std::sin(phi(intnl));
123 991968 : Real cosphi = std::cos(phi(intnl));
124 991968 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
125 991968 : if (std::abs(sin3Lode) <= _sin3tt)
126 : {
127 : // the non-edge-smoothed version
128 : std::vector<Real> eigvals;
129 258416 : stress.symmetricEigenvalues(eigvals);
130 516832 : return mean_stress * sinphi +
131 258416 : std::sqrt(smooth(stress) +
132 258416 : 0.25 * Utility::pow<2>(eigvals[2] - eigvals[0] +
133 258416 : (eigvals[2] + eigvals[0] - 2 * mean_stress) * sinphi)) -
134 258416 : cohesion(intnl) * cosphi;
135 : }
136 : else
137 : {
138 : // the edge-smoothed version
139 : Real aaa, bbb, ccc;
140 733552 : abbo(sin3Lode, sinphi, aaa, bbb, ccc);
141 733552 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
142 733552 : Real sibar2 = stress.secondInvariant();
143 733552 : return mean_stress * sinphi + std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk)) -
144 733552 : cohesion(intnl) * cosphi;
145 : }
146 : }
147 :
148 : RankTwoTensor
149 1652528 : SolidMechanicsPlasticMohrCoulomb::df_dsig(const RankTwoTensor & stress, const Real sin_angle) const
150 : {
151 1652528 : Real mean_stress = stress.trace() / 3.0;
152 1652528 : RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
153 1652528 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
154 1652528 : if (std::abs(sin3Lode) <= _sin3tt)
155 : {
156 : // the non-edge-smoothed version
157 : std::vector<Real> eigvals;
158 : std::vector<RankTwoTensor> deigvals;
159 385648 : stress.dsymmetricEigenvalues(eigvals, deigvals);
160 385648 : Real tmp = eigvals[2] - eigvals[0] + (eigvals[2] + eigvals[0] - 2.0 * mean_stress) * sin_angle;
161 : RankTwoTensor dtmp =
162 385648 : deigvals[2] - deigvals[0] + (deigvals[2] + deigvals[0] - 2.0 * dmean_stress) * sin_angle;
163 385648 : Real denom = std::sqrt(smooth(stress) + 0.25 * Utility::pow<2>(tmp));
164 385648 : return dmean_stress * sin_angle +
165 385648 : (0.5 * dsmooth(stress) * dmean_stress + 0.25 * tmp * dtmp) / denom;
166 : }
167 : else
168 : {
169 : // the edge-smoothed version
170 : Real aaa, bbb, ccc;
171 1266880 : abbo(sin3Lode, sin_angle, aaa, bbb, ccc);
172 1266880 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
173 1266880 : RankTwoTensor dkk = (bbb + 2 * ccc * sin3Lode) * stress.dsin3Lode(_lode_cutoff);
174 1266880 : Real sibar2 = stress.secondInvariant();
175 1266880 : RankTwoTensor dsibar2 = stress.dsecondInvariant();
176 1266880 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
177 1266880 : return dmean_stress * sin_angle + (0.5 * dsmooth(stress) * dmean_stress +
178 1266880 : 0.5 * dsibar2 * Utility::pow<2>(kk) + sibar2 * kk * dkk) /
179 1266880 : denom;
180 : }
181 : }
182 :
183 : RankTwoTensor
184 393216 : SolidMechanicsPlasticMohrCoulomb::dyieldFunction_dstress(const RankTwoTensor & stress,
185 : Real intnl) const
186 : {
187 393216 : Real sinphi = std::sin(phi(intnl));
188 393216 : return df_dsig(stress, sinphi);
189 : }
190 :
191 : Real
192 393216 : SolidMechanicsPlasticMohrCoulomb::dyieldFunction_dintnl(const RankTwoTensor & stress,
193 : Real intnl) const
194 : {
195 393216 : Real sin_angle = std::sin(phi(intnl));
196 393216 : Real cos_angle = std::cos(phi(intnl));
197 393216 : Real dsin_angle = cos_angle * dphi(intnl);
198 393216 : Real dcos_angle = -sin_angle * dphi(intnl);
199 :
200 393216 : Real mean_stress = stress.trace() / 3.0;
201 393216 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
202 393216 : if (std::abs(sin3Lode) <= _sin3tt)
203 : {
204 : // the non-edge-smoothed version
205 : std::vector<Real> eigvals;
206 90320 : stress.symmetricEigenvalues(eigvals);
207 90320 : Real tmp = eigvals[2] - eigvals[0] + (eigvals[2] + eigvals[0] - 2 * mean_stress) * sin_angle;
208 90320 : Real dtmp = (eigvals[2] + eigvals[0] - 2 * mean_stress) * dsin_angle;
209 90320 : Real denom = std::sqrt(smooth(stress) + 0.25 * Utility::pow<2>(tmp));
210 90320 : return mean_stress * dsin_angle + 0.25 * tmp * dtmp / denom - dcohesion(intnl) * cos_angle -
211 90320 : cohesion(intnl) * dcos_angle;
212 : }
213 : else
214 : {
215 : // the edge-smoothed version
216 : Real aaa, bbb, ccc;
217 302896 : abbo(sin3Lode, sin_angle, aaa, bbb, ccc);
218 : Real daaa, dbbb, dccc;
219 302896 : dabbo(sin3Lode, sin_angle, daaa, dbbb, dccc);
220 302896 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
221 302896 : Real dkk = (daaa + dbbb * sin3Lode + dccc * Utility::pow<2>(sin3Lode)) * dsin_angle;
222 302896 : Real sibar2 = stress.secondInvariant();
223 302896 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
224 302896 : return mean_stress * dsin_angle + sibar2 * kk * dkk / denom - dcohesion(intnl) * cos_angle -
225 302896 : cohesion(intnl) * dcos_angle;
226 : }
227 : }
228 :
229 : RankTwoTensor
230 1259312 : SolidMechanicsPlasticMohrCoulomb::flowPotential(const RankTwoTensor & stress, Real intnl) const
231 : {
232 1259312 : Real sinpsi = std::sin(psi(intnl));
233 1259312 : return df_dsig(stress, sinpsi);
234 : }
235 :
236 : RankFourTensor
237 393216 : SolidMechanicsPlasticMohrCoulomb::dflowPotential_dstress(const RankTwoTensor & stress,
238 : Real intnl) const
239 : {
240 393216 : RankFourTensor dr_dstress;
241 393216 : Real sin_angle = std::sin(psi(intnl));
242 393216 : Real mean_stress = stress.trace() / 3.0;
243 393216 : RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
244 393216 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
245 393216 : 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 90320 : stress.dsymmetricEigenvalues(eigvals, deigvals);
252 90320 : stress.d2symmetricEigenvalues(d2eigvals);
253 :
254 90320 : Real tmp = eigvals[2] - eigvals[0] + (eigvals[2] + eigvals[0] - 2.0 * mean_stress) * sin_angle;
255 : RankTwoTensor dtmp =
256 90320 : deigvals[2] - deigvals[0] + (deigvals[2] + deigvals[0] - 2.0 * dmean_stress) * sin_angle;
257 90320 : Real denom = std::sqrt(smooth(stress) + 0.25 * Utility::pow<2>(tmp));
258 : Real denom3 = Utility::pow<3>(denom);
259 90320 : Real d2smooth_over_denom = d2smooth(stress) / denom;
260 90320 : RankTwoTensor numer = dsmooth(stress) * dmean_stress + 0.5 * tmp * dtmp;
261 :
262 90320 : dr_dstress = 0.25 * tmp *
263 180640 : (d2eigvals[2] - d2eigvals[0] + (d2eigvals[2] + d2eigvals[0]) * sin_angle) / denom;
264 :
265 361280 : for (unsigned i = 0; i < 3; ++i)
266 1083840 : for (unsigned j = 0; j < 3; ++j)
267 3251520 : for (unsigned k = 0; k < 3; ++k)
268 9754560 : for (unsigned l = 0; l < 3; ++l)
269 : {
270 7315920 : dr_dstress(i, j, k, l) +=
271 7315920 : 0.5 * d2smooth_over_denom * dmean_stress(i, j) * dmean_stress(k, l);
272 7315920 : dr_dstress(i, j, k, l) += 0.25 * dtmp(i, j) * dtmp(k, l) / denom;
273 7315920 : dr_dstress(i, j, k, l) -= 0.25 * numer(i, j) * numer(k, l) / denom3;
274 : }
275 : }
276 : else
277 : {
278 : // the edge-smoothed version
279 : Real aaa, bbb, ccc;
280 302896 : abbo(sin3Lode, sin_angle, aaa, bbb, ccc);
281 302896 : RankTwoTensor dsin3Lode = stress.dsin3Lode(_lode_cutoff);
282 302896 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
283 302896 : RankTwoTensor dkk = (bbb + 2 * ccc * sin3Lode) * dsin3Lode;
284 302896 : RankFourTensor d2kk = (bbb + 2 * ccc * sin3Lode) * stress.d2sin3Lode(_lode_cutoff);
285 1211584 : for (unsigned i = 0; i < 3; ++i)
286 3634752 : for (unsigned j = 0; j < 3; ++j)
287 10904256 : for (unsigned k = 0; k < 3; ++k)
288 32712768 : for (unsigned l = 0; l < 3; ++l)
289 24534576 : d2kk(i, j, k, l) += 2 * ccc * dsin3Lode(i, j) * dsin3Lode(k, l);
290 :
291 302896 : Real sibar2 = stress.secondInvariant();
292 302896 : RankTwoTensor dsibar2 = stress.dsecondInvariant();
293 302896 : RankFourTensor d2sibar2 = stress.d2secondInvariant();
294 :
295 302896 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
296 : Real denom3 = Utility::pow<3>(denom);
297 302896 : Real d2smooth_over_denom = d2smooth(stress) / denom;
298 : RankTwoTensor numer_full =
299 302896 : 0.5 * dsmooth(stress) * dmean_stress + 0.5 * dsibar2 * kk * kk + sibar2 * kk * dkk;
300 :
301 302896 : dr_dstress = (0.5 * d2sibar2 * Utility::pow<2>(kk) + sibar2 * kk * d2kk) / denom;
302 1211584 : for (unsigned i = 0; i < 3; ++i)
303 3634752 : for (unsigned j = 0; j < 3; ++j)
304 10904256 : for (unsigned k = 0; k < 3; ++k)
305 32712768 : for (unsigned l = 0; l < 3; ++l)
306 : {
307 24534576 : dr_dstress(i, j, k, l) +=
308 24534576 : 0.5 * d2smooth_over_denom * dmean_stress(i, j) * dmean_stress(k, l);
309 24534576 : dr_dstress(i, j, k, l) +=
310 24534576 : (dsibar2(i, j) * dkk(k, l) * kk + dkk(i, j) * dsibar2(k, l) * kk +
311 24534576 : sibar2 * dkk(i, j) * dkk(k, l)) /
312 : denom;
313 24534576 : dr_dstress(i, j, k, l) -= numer_full(i, j) * numer_full(k, l) / denom3;
314 : }
315 : }
316 393216 : return dr_dstress;
317 : }
318 :
319 : RankTwoTensor
320 393216 : SolidMechanicsPlasticMohrCoulomb::dflowPotential_dintnl(const RankTwoTensor & stress,
321 : Real intnl) const
322 : {
323 393216 : Real sin_angle = std::sin(psi(intnl));
324 393216 : Real dsin_angle = std::cos(psi(intnl)) * dpsi(intnl);
325 :
326 393216 : Real mean_stress = stress.trace() / 3.0;
327 393216 : RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
328 393216 : Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
329 :
330 393216 : if (std::abs(sin3Lode) <= _sin3tt)
331 : {
332 : // the non-edge-smoothed version
333 : std::vector<Real> eigvals;
334 : std::vector<RankTwoTensor> deigvals;
335 90320 : stress.dsymmetricEigenvalues(eigvals, deigvals);
336 90320 : Real tmp = eigvals[2] - eigvals[0] + (eigvals[2] + eigvals[0] - 2.0 * mean_stress) * sin_angle;
337 90320 : Real dtmp_dintnl = (eigvals[2] + eigvals[0] - 2 * mean_stress) * dsin_angle;
338 : RankTwoTensor dtmp_dstress =
339 90320 : deigvals[2] - deigvals[0] + (deigvals[2] + deigvals[0] - 2.0 * dmean_stress) * sin_angle;
340 : RankTwoTensor d2tmp_dstress_dintnl =
341 90320 : (deigvals[2] + deigvals[0] - 2.0 * dmean_stress) * dsin_angle;
342 90320 : Real denom = std::sqrt(smooth(stress) + 0.25 * Utility::pow<2>(tmp));
343 90320 : return dmean_stress * dsin_angle + 0.25 * dtmp_dintnl * dtmp_dstress / denom +
344 90320 : 0.25 * tmp * d2tmp_dstress_dintnl / denom -
345 90320 : 0.5 * (dsmooth(stress) * dmean_stress + 0.5 * tmp * dtmp_dstress) * 0.25 * tmp *
346 90320 : dtmp_dintnl / Utility::pow<3>(denom);
347 : }
348 : else
349 : {
350 : // the edge-smoothed version
351 : Real aaa, bbb, ccc;
352 302896 : abbo(sin3Lode, sin_angle, aaa, bbb, ccc);
353 302896 : Real kk = aaa + bbb * sin3Lode + ccc * Utility::pow<2>(sin3Lode);
354 :
355 : Real daaa, dbbb, dccc;
356 302896 : dabbo(sin3Lode, sin_angle, daaa, dbbb, dccc);
357 302896 : Real dkk_dintnl = (daaa + dbbb * sin3Lode + dccc * Utility::pow<2>(sin3Lode)) * dsin_angle;
358 302896 : RankTwoTensor dkk_dstress = (bbb + 2 * ccc * sin3Lode) * stress.dsin3Lode(_lode_cutoff);
359 : RankTwoTensor d2kk_dstress_dintnl =
360 302896 : (dbbb + 2 * dccc * sin3Lode) * stress.dsin3Lode(_lode_cutoff) * dsin_angle;
361 :
362 302896 : Real sibar2 = stress.secondInvariant();
363 302896 : RankTwoTensor dsibar2 = stress.dsecondInvariant();
364 302896 : Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
365 :
366 302896 : return dmean_stress * dsin_angle +
367 302896 : (dsibar2 * kk * dkk_dintnl + sibar2 * dkk_dintnl * dkk_dstress +
368 302896 : sibar2 * kk * d2kk_dstress_dintnl) /
369 : denom -
370 302896 : (0.5 * dsmooth(stress) * dmean_stress + 0.5 * dsibar2 * Utility::pow<2>(kk) +
371 605792 : sibar2 * kk * dkk_dstress) *
372 302896 : sibar2 * kk * dkk_dintnl / Utility::pow<3>(denom);
373 : }
374 : }
375 :
376 : Real
377 1385304 : SolidMechanicsPlasticMohrCoulomb::cohesion(const Real internal_param) const
378 : {
379 1385304 : return _cohesion.value(internal_param);
380 : }
381 :
382 : Real
383 393216 : SolidMechanicsPlasticMohrCoulomb::dcohesion(const Real internal_param) const
384 : {
385 393216 : return _cohesion.derivative(internal_param);
386 : }
387 :
388 : Real
389 3164190 : SolidMechanicsPlasticMohrCoulomb::phi(const Real internal_param) const
390 : {
391 3164190 : return _phi.value(internal_param);
392 : }
393 :
394 : Real
395 786432 : SolidMechanicsPlasticMohrCoulomb::dphi(const Real internal_param) const
396 : {
397 786432 : return _phi.derivative(internal_param);
398 : }
399 :
400 : Real
401 2439566 : SolidMechanicsPlasticMohrCoulomb::psi(const Real internal_param) const
402 : {
403 2439566 : return _psi.value(internal_param);
404 : }
405 :
406 : Real
407 393216 : SolidMechanicsPlasticMohrCoulomb::dpsi(const Real internal_param) const
408 : {
409 393216 : return _psi.derivative(internal_param);
410 : }
411 :
412 : void
413 2909120 : SolidMechanicsPlasticMohrCoulomb::abbo(
414 : const Real sin3lode, const Real sin_angle, Real & aaa, Real & bbb, Real & ccc) const
415 : {
416 2909120 : Real tmp1 = (sin3lode >= 0 ? _costt - sin_angle * _sintt / std::sqrt(3.0)
417 1152352 : : _costt + sin_angle * _sintt / std::sqrt(3.0));
418 2909120 : Real tmp2 = (sin3lode >= 0 ? _sintt + sin_angle * _costt / std::sqrt(3.0)
419 1152352 : : -_sintt + sin_angle * _costt / std::sqrt(3.0));
420 :
421 2909120 : ccc = -_cos3tt * tmp1;
422 2909120 : ccc += (sin3lode >= 0 ? -3 * _sin3tt * tmp2 : 3 * _sin3tt * tmp2);
423 2909120 : ccc /= 18 * Utility::pow<3>(_cos3tt);
424 :
425 2909120 : bbb = (sin3lode >= 0 ? _sin6tt * tmp1 : -_sin6tt * tmp1);
426 2909120 : bbb -= 6 * _cos6tt * tmp2;
427 2909120 : bbb /= 18 * Utility::pow<3>(_cos3tt);
428 :
429 2909120 : aaa = (sin3lode >= 0 ? -sin_angle * _sintt / std::sqrt(3.0) - bbb * _sin3tt
430 1152352 : : sin_angle * _sintt / std::sqrt(3.0) + bbb * _sin3tt);
431 2909120 : aaa += -ccc * Utility::pow<2>(_sin3tt) + _costt;
432 2909120 : }
433 :
434 : void
435 605792 : SolidMechanicsPlasticMohrCoulomb::dabbo(
436 : const Real sin3lode, const Real /*sin_angle*/, Real & daaa, Real & dbbb, Real & dccc) const
437 : {
438 605792 : Real dtmp1 = (sin3lode >= 0 ? -_sintt / std::sqrt(3.0) : _sintt / std::sqrt(3.0));
439 605792 : Real dtmp2 = _costt / std::sqrt(3.0);
440 :
441 605792 : dccc = -_cos3tt * dtmp1;
442 605792 : dccc += (sin3lode >= 0 ? -3 * _sin3tt * dtmp2 : 3 * _sin3tt * dtmp2);
443 605792 : dccc /= 18 * Utility::pow<3>(_cos3tt);
444 :
445 605792 : dbbb = (sin3lode >= 0 ? _sin6tt * dtmp1 : -_sin6tt * dtmp1);
446 605792 : dbbb -= 6 * _cos6tt * dtmp2;
447 605792 : dbbb /= 18 * Utility::pow<3>(_cos3tt);
448 :
449 605792 : daaa = (sin3lode >= 0 ? -_sintt / std::sqrt(3.0) - dbbb * _sin3tt
450 236160 : : _sintt / std::sqrt(3.0) + dbbb * _sin3tt);
451 605792 : daaa += -dccc * Utility::pow<2>(_sin3tt);
452 605792 : }
453 :
454 : Real
455 3824144 : SolidMechanicsPlasticMohrCoulomb::smooth(const RankTwoTensor & stress) const
456 : {
457 3824144 : Real smoother2 = _small_smoother2;
458 3824144 : if (_tip_scheme == "cap")
459 : {
460 486912 : Real x = stress.trace() / 3.0 - _cap_start;
461 : Real p = 0;
462 486912 : if (x > 0)
463 436704 : p = x * (1 - std::exp(-_cap_rate * x));
464 486912 : smoother2 += Utility::pow<2>(p);
465 : }
466 3824144 : return smoother2;
467 : }
468 :
469 : Real
470 2438960 : SolidMechanicsPlasticMohrCoulomb::dsmooth(const RankTwoTensor & stress) const
471 : {
472 : Real dsmoother2 = 0;
473 2438960 : if (_tip_scheme == "cap")
474 : {
475 312192 : Real x = stress.trace() / 3.0 - _cap_start;
476 : Real p = 0;
477 : Real dp_dx = 0;
478 312192 : if (x > 0)
479 : {
480 279456 : p = x * (1 - std::exp(-_cap_rate * x));
481 279456 : dp_dx = (1 - std::exp(-_cap_rate * x)) + x * _cap_rate * std::exp(-_cap_rate * x);
482 : }
483 312192 : dsmoother2 += 2 * p * dp_dx;
484 : }
485 2438960 : return dsmoother2;
486 : }
487 :
488 : Real
489 393216 : SolidMechanicsPlasticMohrCoulomb::d2smooth(const RankTwoTensor & stress) const
490 : {
491 : Real d2smoother2 = 0;
492 393216 : if (_tip_scheme == "cap")
493 : {
494 49920 : Real x = stress.trace() / 3.0 - _cap_start;
495 : Real p = 0;
496 : Real dp_dx = 0;
497 : Real d2p_dx2 = 0;
498 49920 : if (x > 0)
499 : {
500 44064 : p = x * (1 - std::exp(-_cap_rate * x));
501 44064 : dp_dx = (1 - std::exp(-_cap_rate * x)) + x * _cap_rate * std::exp(-_cap_rate * x);
502 44064 : d2p_dx2 = 2 * _cap_rate * std::exp(-_cap_rate * x) -
503 44064 : x * Utility::pow<2>(_cap_rate) * std::exp(-_cap_rate * x);
504 : }
505 49920 : d2smoother2 += 2 * Utility::pow<2>(dp_dx) + 2 * p * d2p_dx2;
506 : }
507 393216 : return d2smoother2;
508 : }
509 :
510 : std::string
511 0 : SolidMechanicsPlasticMohrCoulomb::modelName() const
512 : {
513 0 : return "MohrCoulomb";
514 : }
|
