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