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