www.mooseframework.org
TensorMechanicsPlasticMohrCoulombMulti.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 
9 // Following is for perturbing eigvenvalues. This looks really bodgy, but works quite well!
10 #include "MooseRandom.h"
11 #include "libmesh/utility.h"
12 
13 template <>
14 InputParameters
16 {
17  InputParameters params = validParams<TensorMechanicsPlasticModel>();
18  params.addClassDescription("Non-associative Mohr-Coulomb plasticity with hardening/softening");
19  params.addRequiredParam<UserObjectName>(
20  "cohesion", "A TensorMechanicsHardening UserObject that defines hardening of the cohesion");
21  params.addRequiredParam<UserObjectName>("friction_angle",
22  "A TensorMechanicsHardening UserObject "
23  "that defines hardening of the "
24  "friction angle (in radians)");
25  params.addRequiredParam<UserObjectName>("dilation_angle",
26  "A TensorMechanicsHardening UserObject "
27  "that defines hardening of the "
28  "dilation angle (in radians)");
29  params.addParam<unsigned int>("max_iterations",
30  10,
31  "Maximum number of Newton-Raphson iterations "
32  "allowed in the custom return-map algorithm. "
33  " For highly nonlinear hardening this may "
34  "need to be higher than 10.");
35  params.addParam<Real>("shift",
36  "Yield surface is shifted by this amount to avoid problems with "
37  "defining derivatives when eigenvalues are equal. If this is "
38  "larger than f_tol, a warning will be issued. This may be set "
39  "very small when using the custom returnMap. Default = f_tol.");
40  params.addParam<bool>("use_custom_returnMap",
41  true,
42  "Use a custom return-map algorithm for this "
43  "plasticity model, which may speed up "
44  "computations considerably. Set to true "
45  "only for isotropic elasticity with no "
46  "hardening of the dilation angle. In this "
47  "case you may set 'shift' very small.");
48 
49  return params;
50 }
51 
53  const InputParameters & parameters)
54  : TensorMechanicsPlasticModel(parameters),
55  _cohesion(getUserObject<TensorMechanicsHardeningModel>("cohesion")),
56  _phi(getUserObject<TensorMechanicsHardeningModel>("friction_angle")),
57  _psi(getUserObject<TensorMechanicsHardeningModel>("dilation_angle")),
58  _max_iters(getParam<unsigned int>("max_iterations")),
59  _shift(parameters.isParamValid("shift") ? getParam<Real>("shift") : _f_tol),
60  _use_custom_returnMap(getParam<bool>("use_custom_returnMap"))
61 {
62  if (_shift < 0)
63  mooseError("Value of 'shift' in TensorMechanicsPlasticMohrCoulombMulti must not be negative\n");
64  if (_shift > _f_tol)
65  _console << "WARNING: value of 'shift' in TensorMechanicsPlasticMohrCoulombMulti is probably "
66  "set too high\n";
67  if (LIBMESH_DIM != 3)
68  mooseError("TensorMechanicsPlasticMohrCoulombMulti is only defined for LIBMESH_DIM=3");
69  MooseRandom::seed(0);
70 }
71 
72 unsigned int
74 {
75  return 6;
76 }
77 
78 void
80  Real intnl,
81  std::vector<Real> & f) const
82 {
83  std::vector<Real> eigvals;
84  stress.symmetricEigenvalues(eigvals);
85  eigvals[0] += _shift;
86  eigvals[2] -= _shift;
87 
88  const Real sinphi = std::sin(phi(intnl));
89  const Real cosphi = std::cos(phi(intnl));
90  const Real cohcos = cohesion(intnl) * cosphi;
91 
92  yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, f);
93 }
94 
95 void
97  Real e0, Real e1, Real e2, Real sinphi, Real cohcos, std::vector<Real> & f) const
98 {
99  // Naively it seems a shame to have 6 yield functions active instead of just
100  // 3. But 3 won't do. Eg, think of a loading with eigvals[0]=eigvals[1]=eigvals[2]
101  // Then to return to the yield surface would require 2 positive plastic multipliers
102  // and one negative one. Boo hoo.
103 
104  f.resize(6);
105  f[0] = 0.5 * (e0 - e1) + 0.5 * (e0 + e1) * sinphi - cohcos;
106  f[1] = 0.5 * (e1 - e0) + 0.5 * (e0 + e1) * sinphi - cohcos;
107  f[2] = 0.5 * (e0 - e2) + 0.5 * (e0 + e2) * sinphi - cohcos;
108  f[3] = 0.5 * (e2 - e0) + 0.5 * (e0 + e2) * sinphi - cohcos;
109  f[4] = 0.5 * (e1 - e2) + 0.5 * (e1 + e2) * sinphi - cohcos;
110  f[5] = 0.5 * (e2 - e1) + 0.5 * (e1 + e2) * sinphi - cohcos;
111 }
112 
113 void
115  std::vector<Real> & eigvals,
116  std::vector<RankTwoTensor> & deigvals) const
117 {
118  Real small_perturbation;
119  RankTwoTensor shifted_stress = stress;
120  while (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
121  {
122  for (unsigned i = 0; i < 3; ++i)
123  for (unsigned j = 0; j <= i; ++j)
124  {
125  small_perturbation = 0.1 * _shift * 2 * (MooseRandom::rand() - 0.5);
126  shifted_stress(i, j) += small_perturbation;
127  shifted_stress(j, i) += small_perturbation;
128  }
129  shifted_stress.dsymmetricEigenvalues(eigvals, deigvals);
130  }
131 }
132 
133 void
135  Real sin_angle,
136  std::vector<RankTwoTensor> & df) const
137 {
138  std::vector<Real> eigvals;
139  std::vector<RankTwoTensor> deigvals;
140  stress.dsymmetricEigenvalues(eigvals, deigvals);
141 
142  if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
143  perturbStress(stress, eigvals, deigvals);
144 
145  df.resize(6);
146  df[0] = 0.5 * (deigvals[0] - deigvals[1]) + 0.5 * (deigvals[0] + deigvals[1]) * sin_angle;
147  df[1] = 0.5 * (deigvals[1] - deigvals[0]) + 0.5 * (deigvals[0] + deigvals[1]) * sin_angle;
148  df[2] = 0.5 * (deigvals[0] - deigvals[2]) + 0.5 * (deigvals[0] + deigvals[2]) * sin_angle;
149  df[3] = 0.5 * (deigvals[2] - deigvals[0]) + 0.5 * (deigvals[0] + deigvals[2]) * sin_angle;
150  df[4] = 0.5 * (deigvals[1] - deigvals[2]) + 0.5 * (deigvals[1] + deigvals[2]) * sin_angle;
151  df[5] = 0.5 * (deigvals[2] - deigvals[1]) + 0.5 * (deigvals[1] + deigvals[2]) * sin_angle;
152 }
153 
154 void
156  const RankTwoTensor & stress, Real intnl, std::vector<RankTwoTensor> & df_dstress) const
157 {
158  const Real sinphi = std::sin(phi(intnl));
159  df_dsig(stress, sinphi, df_dstress);
160 }
161 
162 void
164  Real intnl,
165  std::vector<Real> & df_dintnl) const
166 {
167  std::vector<Real> eigvals;
168  stress.symmetricEigenvalues(eigvals);
169  eigvals[0] += _shift;
170  eigvals[2] -= _shift;
171 
172  const Real sin_angle = std::sin(phi(intnl));
173  const Real cos_angle = std::cos(phi(intnl));
174  const Real dsin_angle = cos_angle * dphi(intnl);
175  const Real dcos_angle = -sin_angle * dphi(intnl);
176  const Real dcohcos = dcohesion(intnl) * cos_angle + cohesion(intnl) * dcos_angle;
177 
178  df_dintnl.resize(6);
179  df_dintnl[0] = df_dintnl[1] = 0.5 * (eigvals[0] + eigvals[1]) * dsin_angle - dcohcos;
180  df_dintnl[2] = df_dintnl[3] = 0.5 * (eigvals[0] + eigvals[2]) * dsin_angle - dcohcos;
181  df_dintnl[4] = df_dintnl[5] = 0.5 * (eigvals[1] + eigvals[2]) * dsin_angle - dcohcos;
182 }
183 
184 void
186  Real intnl,
187  std::vector<RankTwoTensor> & r) const
188 {
189  const Real sinpsi = std::sin(psi(intnl));
190  df_dsig(stress, sinpsi, r);
191 }
192 
193 void
195  const RankTwoTensor & stress, Real intnl, std::vector<RankFourTensor> & dr_dstress) const
196 {
197  std::vector<RankFourTensor> d2eigvals;
198  stress.d2symmetricEigenvalues(d2eigvals);
199 
200  const Real sinpsi = std::sin(psi(intnl));
201 
202  dr_dstress.resize(6);
203  dr_dstress[0] =
204  0.5 * (d2eigvals[0] - d2eigvals[1]) + 0.5 * (d2eigvals[0] + d2eigvals[1]) * sinpsi;
205  dr_dstress[1] =
206  0.5 * (d2eigvals[1] - d2eigvals[0]) + 0.5 * (d2eigvals[0] + d2eigvals[1]) * sinpsi;
207  dr_dstress[2] =
208  0.5 * (d2eigvals[0] - d2eigvals[2]) + 0.5 * (d2eigvals[0] + d2eigvals[2]) * sinpsi;
209  dr_dstress[3] =
210  0.5 * (d2eigvals[2] - d2eigvals[0]) + 0.5 * (d2eigvals[0] + d2eigvals[2]) * sinpsi;
211  dr_dstress[4] =
212  0.5 * (d2eigvals[1] - d2eigvals[2]) + 0.5 * (d2eigvals[1] + d2eigvals[2]) * sinpsi;
213  dr_dstress[5] =
214  0.5 * (d2eigvals[2] - d2eigvals[1]) + 0.5 * (d2eigvals[1] + d2eigvals[2]) * sinpsi;
215 }
216 
217 void
219  const RankTwoTensor & stress, Real intnl, std::vector<RankTwoTensor> & dr_dintnl) const
220 {
221  const Real cos_angle = std::cos(psi(intnl));
222  const Real dsin_angle = cos_angle * dpsi(intnl);
223 
224  std::vector<Real> eigvals;
225  std::vector<RankTwoTensor> deigvals;
226  stress.dsymmetricEigenvalues(eigvals, deigvals);
227 
228  if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
229  perturbStress(stress, eigvals, deigvals);
230 
231  dr_dintnl.resize(6);
232  dr_dintnl[0] = dr_dintnl[1] = 0.5 * (deigvals[0] + deigvals[1]) * dsin_angle;
233  dr_dintnl[2] = dr_dintnl[3] = 0.5 * (deigvals[0] + deigvals[2]) * dsin_angle;
234  dr_dintnl[4] = dr_dintnl[5] = 0.5 * (deigvals[1] + deigvals[2]) * dsin_angle;
235 }
236 
237 void
239  const RankTwoTensor & stress,
240  Real intnl,
241  const RankFourTensor & Eijkl,
242  std::vector<bool> & act,
243  RankTwoTensor & returned_stress) const
244 {
245  act.assign(6, false);
246 
247  if (f[0] <= _f_tol && f[1] <= _f_tol && f[2] <= _f_tol && f[3] <= _f_tol && f[4] <= _f_tol &&
248  f[5] <= _f_tol)
249  {
250  returned_stress = stress;
251  return;
252  }
253 
254  Real returned_intnl;
255  std::vector<Real> dpm(6);
256  RankTwoTensor delta_dp;
257  std::vector<Real> yf(6);
258  bool trial_stress_inadmissible;
259  doReturnMap(stress,
260  intnl,
261  Eijkl,
262  0.0,
263  returned_stress,
264  returned_intnl,
265  dpm,
266  delta_dp,
267  yf,
268  trial_stress_inadmissible);
269 
270  for (unsigned i = 0; i < 6; ++i)
271  act[i] = (dpm[i] > 0);
272 }
273 
274 Real
275 TensorMechanicsPlasticMohrCoulombMulti::cohesion(const Real internal_param) const
276 {
277  return _cohesion.value(internal_param);
278 }
279 
280 Real
282 {
283  return _cohesion.derivative(internal_param);
284 }
285 
286 Real
287 TensorMechanicsPlasticMohrCoulombMulti::phi(const Real internal_param) const
288 {
289  return _phi.value(internal_param);
290 }
291 
292 Real
293 TensorMechanicsPlasticMohrCoulombMulti::dphi(const Real internal_param) const
294 {
295  return _phi.derivative(internal_param);
296 }
297 
298 Real
299 TensorMechanicsPlasticMohrCoulombMulti::psi(const Real internal_param) const
300 {
301  return _psi.value(internal_param);
302 }
303 
304 Real
305 TensorMechanicsPlasticMohrCoulombMulti::dpsi(const Real internal_param) const
306 {
307  return _psi.derivative(internal_param);
308 }
309 
310 std::string
312 {
313  return "MohrCoulombMulti";
314 }
315 
316 bool
317 TensorMechanicsPlasticMohrCoulombMulti::returnMap(const RankTwoTensor & trial_stress,
318  Real intnl_old,
319  const RankFourTensor & E_ijkl,
320  Real ep_plastic_tolerance,
321  RankTwoTensor & returned_stress,
322  Real & returned_intnl,
323  std::vector<Real> & dpm,
324  RankTwoTensor & delta_dp,
325  std::vector<Real> & yf,
326  bool & trial_stress_inadmissible) const
327 {
329  return TensorMechanicsPlasticModel::returnMap(trial_stress,
330  intnl_old,
331  E_ijkl,
332  ep_plastic_tolerance,
333  returned_stress,
334  returned_intnl,
335  dpm,
336  delta_dp,
337  yf,
338  trial_stress_inadmissible);
339 
340  return doReturnMap(trial_stress,
341  intnl_old,
342  E_ijkl,
343  ep_plastic_tolerance,
344  returned_stress,
345  returned_intnl,
346  dpm,
347  delta_dp,
348  yf,
349  trial_stress_inadmissible);
350 }
351 
352 bool
353 TensorMechanicsPlasticMohrCoulombMulti::doReturnMap(const RankTwoTensor & trial_stress,
354  Real intnl_old,
355  const RankFourTensor & E_ijkl,
356  Real ep_plastic_tolerance,
357  RankTwoTensor & returned_stress,
358  Real & returned_intnl,
359  std::vector<Real> & dpm,
360  RankTwoTensor & delta_dp,
361  std::vector<Real> & yf,
362  bool & trial_stress_inadmissible) const
363 {
364  mooseAssert(dpm.size() == 6,
365  "TensorMechanicsPlasticMohrCoulombMulti size of dpm should be 6 but it is "
366  << dpm.size());
367 
368  std::vector<Real> eigvals;
369  RankTwoTensor eigvecs;
370  trial_stress.symmetricEigenvaluesEigenvectors(eigvals, eigvecs);
371  eigvals[0] += _shift;
372  eigvals[2] -= _shift;
373 
374  Real sinphi = std::sin(phi(intnl_old));
375  Real cosphi = std::cos(phi(intnl_old));
376  Real coh = cohesion(intnl_old);
377  Real cohcos = coh * cosphi;
378 
379  yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, yf);
380 
381  if (yf[0] <= _f_tol && yf[1] <= _f_tol && yf[2] <= _f_tol && yf[3] <= _f_tol && yf[4] <= _f_tol &&
382  yf[5] <= _f_tol)
383  {
384  // purely elastic (trial_stress, intnl_old)
385  trial_stress_inadmissible = false;
386  return true;
387  }
388 
389  trial_stress_inadmissible = true;
390  delta_dp.zero();
391  returned_stress = RankTwoTensor();
392 
393  // these are the normals to the 6 yield surfaces, which are const because of the assumption of no
394  // psi hardening
395  std::vector<RealVectorValue> norm(6);
396  const Real sinpsi = std::sin(psi(intnl_old));
397  const Real oneminus = 0.5 * (1 - sinpsi);
398  const Real oneplus = 0.5 * (1 + sinpsi);
399  norm[0](0) = oneplus;
400  norm[0](1) = -oneminus;
401  norm[0](2) = 0;
402  norm[1](0) = -oneminus;
403  norm[1](1) = oneplus;
404  norm[1](2) = 0;
405  norm[2](0) = oneplus;
406  norm[2](1) = 0;
407  norm[2](2) = -oneminus;
408  norm[3](0) = -oneminus;
409  norm[3](1) = 0;
410  norm[3](2) = oneplus;
411  norm[4](0) = 0;
412  norm[4](1) = oneplus;
413  norm[4](2) = -oneminus;
414  norm[5](0) = 0;
415  norm[5](1) = -oneminus;
416  norm[5](2) = oneplus;
417 
418  // the flow directions are these norm multiplied by Eijkl.
419  // I call the flow directions "n".
420  // In the following I assume that the Eijkl is
421  // for an isotropic situation. Then I don't have to
422  // rotate to the principal-stress frame, and i don't
423  // have to worry about strange off-diagonal things
424  std::vector<RealVectorValue> n(6);
425  for (unsigned ys = 0; ys < 6; ++ys)
426  for (unsigned i = 0; i < 3; ++i)
427  for (unsigned j = 0; j < 3; ++j)
428  n[ys](i) += E_ijkl(i, i, j, j) * norm[ys](j);
429  const Real mag_E = E_ijkl(0, 0, 0, 0);
430 
431  // With non-zero Poisson's ratio and hardening
432  // it is not computationally cheap to know whether
433  // the trial stress will return to the tip, edge,
434  // or plane. The following at least
435  // gives a not-completely-stupid guess
436  // trial_order[0] = type of return to try first
437  // trial_order[1] = type of return to try second
438  // trial_order[2] = type of return to try third
439  // trial_order[3] = type of return to try fourth
440  // trial_order[4] = type of return to try fifth
441  // In the following the "binary" stuff indicates the
442  // deactive (0) and active (1) surfaces, eg
443  // 110100 means that surfaces 0, 1 and 3 are active
444  // and 2, 4 and 5 are deactive
445  const unsigned int number_of_return_paths = 5;
446  std::vector<int> trial_order(number_of_return_paths);
447  if (yf[1] > _f_tol && yf[3] > _f_tol && yf[5] > _f_tol)
448  {
449  trial_order[0] = tip110100;
450  trial_order[1] = edge010100;
451  trial_order[2] = plane000100;
452  trial_order[3] = edge000101;
453  trial_order[4] = tip010101;
454  }
455  else if (yf[1] <= _f_tol && yf[3] > _f_tol && yf[5] > _f_tol)
456  {
457  trial_order[0] = edge000101;
458  trial_order[1] = plane000100;
459  trial_order[2] = tip110100;
460  trial_order[3] = tip010101;
461  trial_order[4] = edge010100;
462  }
463  else if (yf[1] <= _f_tol && yf[3] > _f_tol && yf[5] <= _f_tol)
464  {
465  trial_order[0] = plane000100;
466  trial_order[1] = edge000101;
467  trial_order[2] = edge010100;
468  trial_order[3] = tip110100;
469  trial_order[4] = tip010101;
470  }
471  else
472  {
473  trial_order[0] = edge010100;
474  trial_order[1] = plane000100;
475  trial_order[2] = edge000101;
476  trial_order[3] = tip110100;
477  trial_order[4] = tip010101;
478  }
479 
480  unsigned trial;
481  bool nr_converged = false;
482  bool kt_success = false;
483  std::vector<RealVectorValue> ntip(3);
484  std::vector<Real> dpmtip(3);
485 
486  for (trial = 0; trial < number_of_return_paths; ++trial)
487  {
488  switch (trial_order[trial])
489  {
490  case tip110100:
491  for (unsigned int i = 0; i < 3; ++i)
492  {
493  ntip[0](i) = n[0](i);
494  ntip[1](i) = n[1](i);
495  ntip[2](i) = n[3](i);
496  }
497  kt_success = returnTip(eigvals,
498  ntip,
499  dpmtip,
500  returned_stress,
501  intnl_old,
502  sinphi,
503  cohcos,
504  0,
505  nr_converged,
506  ep_plastic_tolerance,
507  yf);
508  if (nr_converged && kt_success)
509  {
510  dpm[0] = dpmtip[0];
511  dpm[1] = dpmtip[1];
512  dpm[3] = dpmtip[2];
513  dpm[2] = dpm[4] = dpm[5] = 0;
514  }
515  break;
516 
517  case tip010101:
518  for (unsigned int i = 0; i < 3; ++i)
519  {
520  ntip[0](i) = n[1](i);
521  ntip[1](i) = n[3](i);
522  ntip[2](i) = n[5](i);
523  }
524  kt_success = returnTip(eigvals,
525  ntip,
526  dpmtip,
527  returned_stress,
528  intnl_old,
529  sinphi,
530  cohcos,
531  0,
532  nr_converged,
533  ep_plastic_tolerance,
534  yf);
535  if (nr_converged && kt_success)
536  {
537  dpm[1] = dpmtip[0];
538  dpm[3] = dpmtip[1];
539  dpm[5] = dpmtip[2];
540  dpm[0] = dpm[2] = dpm[4] = 0;
541  }
542  break;
543 
544  case edge000101:
545  kt_success = returnEdge000101(eigvals,
546  n,
547  dpm,
548  returned_stress,
549  intnl_old,
550  sinphi,
551  cohcos,
552  0,
553  mag_E,
554  nr_converged,
555  ep_plastic_tolerance,
556  yf);
557  break;
558 
559  case edge010100:
560  kt_success = returnEdge010100(eigvals,
561  n,
562  dpm,
563  returned_stress,
564  intnl_old,
565  sinphi,
566  cohcos,
567  0,
568  mag_E,
569  nr_converged,
570  ep_plastic_tolerance,
571  yf);
572  break;
573 
574  case plane000100:
575  kt_success = returnPlane(eigvals,
576  n,
577  dpm,
578  returned_stress,
579  intnl_old,
580  sinphi,
581  cohcos,
582  0,
583  nr_converged,
584  ep_plastic_tolerance,
585  yf);
586  break;
587  }
588 
589  if (nr_converged && kt_success)
590  break;
591  }
592 
593  if (trial == number_of_return_paths)
594  {
595  sinphi = std::sin(phi(intnl_old));
596  cosphi = std::cos(phi(intnl_old));
597  coh = cohesion(intnl_old);
598  cohcos = coh * cosphi;
599  yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, yf);
600  Moose::err << "Trial stress = \n";
601  trial_stress.print(Moose::err);
602  Moose::err << "which has eigenvalues = " << eigvals[0] << " " << eigvals[1] << " " << eigvals[2]
603  << "\n";
604  Moose::err << "and yield functions = " << yf[0] << " " << yf[1] << " " << yf[2] << " " << yf[3]
605  << " " << yf[4] << " " << yf[5] << "\n";
606  Moose::err << "Internal parameter = " << intnl_old << "\n";
607  mooseError("TensorMechanicsPlasticMohrCoulombMulti: FAILURE! You probably need to implement a "
608  "line search if your hardening is too severe, or you need to tune your tolerances "
609  "(eg, yield_function_tolerance should be a little smaller than (young "
610  "modulus)*ep_plastic_tolerance).\n");
611  return false;
612  }
613 
614  // success
615 
616  returned_intnl = intnl_old;
617  for (unsigned i = 0; i < 6; ++i)
618  returned_intnl += dpm[i];
619  for (unsigned i = 0; i < 6; ++i)
620  for (unsigned j = 0; j < 3; ++j)
621  delta_dp(j, j) += dpm[i] * norm[i](j);
622  returned_stress = eigvecs * returned_stress * (eigvecs.transpose());
623  delta_dp = eigvecs * delta_dp * (eigvecs.transpose());
624  return true;
625 }
626 
627 bool
628 TensorMechanicsPlasticMohrCoulombMulti::returnTip(const std::vector<Real> & eigvals,
629  const std::vector<RealVectorValue> & n,
630  std::vector<Real> & dpm,
631  RankTwoTensor & returned_stress,
632  Real intnl_old,
633  Real & sinphi,
634  Real & cohcos,
635  Real initial_guess,
636  bool & nr_converged,
637  Real ep_plastic_tolerance,
638  std::vector<Real> & yf) const
639 {
640  // This returns to the Mohr-Coulomb tip using the THREE directions
641  // given in n, and yields the THREE dpm values. Note that you
642  // must supply THREE suitable n vectors out of the total of SIX
643  // flow directions, and then interpret the THREE dpm values appropriately.
644  //
645  // Eg1. You supply the flow directions n[0], n[1] and n[3] as
646  // the "n" vectors. This is return-to-the-tip via 110100.
647  // Then the three returned dpm values will be dpm[0], dpm[1] and dpm[3].
648 
649  // Eg2. You supply the flow directions n[1], n[3] and n[5] as
650  // the "n" vectors. This is return-to-the-tip via 010101.
651  // Then the three returned dpm values will be dpm[1], dpm[3] and dpm[5].
652 
653  // The returned point is defined by the three yield functions (corresonding
654  // to the three supplied flow directions) all being zero.
655  // that is, returned_stress = diag(cohcot, cohcot, cohcot), where
656  // cohcot = cohesion*cosphi/sinphi
657  // where intnl = intnl_old + dpm[0] + dpm[1] + dpm[2]
658  // The 3 plastic multipliers, dpm, are defiend by the normality condition
659  // eigvals - cohcot = dpm[0]*n[0] + dpm[1]*n[1] + dpm[2]*n[2]
660  // (Kuhn-Tucker demands that all dpm are non-negative, but we leave
661  // that checking for the end.)
662  // (Remember that these "n" vectors and "dpm" values must be interpreted
663  // differently depending on what you pass into this function.)
664  // This is a vector equation with solution (A):
665  // dpm[0] = triple(eigvals - cohcot, n[1], n[2])/trip;
666  // dpm[1] = triple(eigvals - cohcot, n[2], n[0])/trip;
667  // dpm[2] = triple(eigvals - cohcot, n[0], n[1])/trip;
668  // where trip = triple(n[0], n[1], n[2]).
669  // By adding the three components of that solution together
670  // we can get an equation for x = dpm[0] + dpm[1] + dpm[2],
671  // and then our Newton-Raphson only involves one variable (x).
672  // In the following, i specialise to the isotropic situation.
673 
674  mooseAssert(n.size() == 3,
675  "TensorMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
676  "supplied with n of size 3, whereas yours is "
677  << n.size());
678  mooseAssert(dpm.size() == 3,
679  "TensorMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
680  "supplied with dpm of size 3, whereas yours is "
681  << dpm.size());
682  mooseAssert(yf.size() == 6,
683  "TensorMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
684  "supplied with yf of size 6, whereas yours is "
685  << yf.size());
686 
687  Real x = initial_guess;
688  const Real trip = triple_product(n[0], n[1], n[2]);
689  sinphi = std::sin(phi(intnl_old + x));
690  Real cosphi = std::cos(phi(intnl_old + x));
691  Real coh = cohesion(intnl_old + x);
692  cohcos = coh * cosphi;
693  Real cohcot = cohcos / sinphi;
694 
695  if (_cohesion.modelName().compare("Constant") != 0 || _phi.modelName().compare("Constant") != 0)
696  {
697  // Finding x is expensive. Therefore
698  // although x!=0 for Constant Hardening, solution (A)
699  // demonstrates that we don't
700  // actually need to know x to find the dpm for
701  // Constant Hardening.
702  //
703  // However, for nontrivial Hardening, the following
704  // is necessary
705  // cohcot_coeff = [1,1,1].(Cross[n[1], n[2]] + Cross[n[2], n[0]] + Cross[n[0], n[1]])/trip
706  Real cohcot_coeff =
707  (n[0](0) * (n[1](1) - n[1](2) - n[2](1)) + (n[1](2) - n[1](1)) * n[2](0) +
708  (n[1](0) - n[1](2)) * n[2](1) + n[0](2) * (n[1](0) - n[1](1) - n[2](0) + n[2](1)) +
709  n[0](1) * (n[1](2) - n[1](0) + n[2](0) - n[2](2)) +
710  (n[0](0) - n[1](0) + n[1](1)) * n[2](2)) /
711  trip;
712  // eig_term = eigvals.(Cross[n[1], n[2]] + Cross[n[2], n[0]] + Cross[n[0], n[1]])/trip
713  Real eig_term = eigvals[0] * (-n[0](2) * n[1](1) + n[0](1) * n[1](2) + n[0](2) * n[2](1) -
714  n[1](2) * n[2](1) - n[0](1) * n[2](2) + n[1](1) * n[2](2)) /
715  trip;
716  eig_term += eigvals[1] * (n[0](2) * n[1](0) - n[0](0) * n[1](2) - n[0](2) * n[2](0) +
717  n[1](2) * n[2](0) + n[0](0) * n[2](2) - n[1](0) * n[2](2)) /
718  trip;
719  eig_term += eigvals[2] * (n[0](0) * n[1](1) - n[1](1) * n[2](0) + n[0](1) * n[2](0) -
720  n[0](1) * n[1](0) - n[0](0) * n[2](1) + n[1](0) * n[2](1)) /
721  trip;
722  // and finally, the equation we want to solve is:
723  // x - eig_term + cohcot*cohcot_coeff = 0
724  // but i divide by cohcot_coeff so the result has the units of
725  // stress, so using _f_tol as a convergence check is reasonable
726  eig_term /= cohcot_coeff;
727  Real residual = x / cohcot_coeff - eig_term + cohcot;
728  Real jacobian;
729  Real deriv_phi;
730  Real deriv_coh;
731  unsigned int iter = 0;
732  do
733  {
734  deriv_phi = dphi(intnl_old + x);
735  deriv_coh = dcohesion(intnl_old + x);
736  jacobian = 1.0 / cohcot_coeff + deriv_coh * cosphi / sinphi -
737  coh * deriv_phi / Utility::pow<2>(sinphi);
738  x += -residual / jacobian;
739 
740  if (iter > _max_iters) // not converging
741  {
742  nr_converged = false;
743  return false;
744  }
745 
746  sinphi = std::sin(phi(intnl_old + x));
747  cosphi = std::cos(phi(intnl_old + x));
748  coh = cohesion(intnl_old + x);
749  cohcos = coh * cosphi;
750  cohcot = cohcos / sinphi;
751  residual = x / cohcot_coeff - eig_term + cohcot;
752  iter++;
753  } while (residual * residual > _f_tol * _f_tol / 100);
754  }
755 
756  // so the NR process converged, but we must
757  // calculate the individual dpm values and
758  // check Kuhn-Tucker
759  nr_converged = true;
760  if (x < -3 * ep_plastic_tolerance)
761  // obviously at least one of the dpm are < -ep_plastic_tolerance. No point in proceeding. This
762  // is a potential weak-point: if the user has set _f_tol quite large, and ep_plastic_tolerance
763  // quite small, the above NR process will quickly converge, but the solution may be wrong and
764  // violate Kuhn-Tucker.
765  return false;
766 
767  // The following is the solution (A) written above
768  // (dpm[0] = triple(eigvals - cohcot, n[1], n[2])/trip, etc)
769  // in the isotropic situation
770  RealVectorValue v;
771  v(0) = eigvals[0] - cohcot;
772  v(1) = eigvals[1] - cohcot;
773  v(2) = eigvals[2] - cohcot;
774  dpm[0] = triple_product(v, n[1], n[2]) / trip;
775  dpm[1] = triple_product(v, n[2], n[0]) / trip;
776  dpm[2] = triple_product(v, n[0], n[1]) / trip;
777 
778  if (dpm[0] < -ep_plastic_tolerance || dpm[1] < -ep_plastic_tolerance ||
779  dpm[2] < -ep_plastic_tolerance)
780  // Kuhn-Tucker failure. No point in proceeding
781  return false;
782 
783  // Kuhn-Tucker has succeeded: just need returned_stress and yf values
784  // I do not use the dpm to calculate returned_stress, because that
785  // might add error (and computational effort), simply:
786  returned_stress(0, 0) = returned_stress(1, 1) = returned_stress(2, 2) = cohcot;
787  // So by construction the yield functions are all zero
788  yf[0] = yf[1] = yf[2] = yf[3] = yf[4] = yf[5] = 0;
789  return true;
790 }
791 
792 bool
793 TensorMechanicsPlasticMohrCoulombMulti::returnPlane(const std::vector<Real> & eigvals,
794  const std::vector<RealVectorValue> & n,
795  std::vector<Real> & dpm,
796  RankTwoTensor & returned_stress,
797  Real intnl_old,
798  Real & sinphi,
799  Real & cohcos,
800  Real initial_guess,
801  bool & nr_converged,
802  Real ep_plastic_tolerance,
803  std::vector<Real> & yf) const
804 {
805  // This returns to the Mohr-Coulomb plane using n[3] (ie 000100)
806  //
807  // The returned point is defined by the f[3]=0 and
808  // a = eigvals - dpm[3]*n[3]
809  // where "a" is the returned point and dpm[3] is the plastic multiplier.
810  // This equation is a vector equation in principal stress space.
811  // (Kuhn-Tucker also demands that dpm[3]>=0, but we leave checking
812  // that condition for the end.)
813  // Since f[3]=0, we must have
814  // a[2]*(1+sinphi) + a[0]*(-1+sinphi) - 2*coh*cosphi = 0
815  // which gives dpm[3] as the solution of
816  // alpha*dpm[3] + eigvals[2] - eigvals[0] + beta*sinphi - 2*coh*cosphi = 0
817  // with alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0))*sinphi
818  // beta = eigvals[2] + eigvals[0]
819 
820  mooseAssert(n.size() == 6,
821  "TensorMechanicsPlasticMohrCoulombMulti: Custom plane-return algorithm must be "
822  "supplied with n of size 6, whereas yours is "
823  << n.size());
824  mooseAssert(dpm.size() == 6,
825  "TensorMechanicsPlasticMohrCoulombMulti: Custom plane-return algorithm must be "
826  "supplied with dpm of size 6, whereas yours is "
827  << dpm.size());
828  mooseAssert(yf.size() == 6,
829  "TensorMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
830  "supplied with yf of size 6, whereas yours is "
831  << yf.size());
832 
833  dpm[3] = initial_guess;
834  sinphi = std::sin(phi(intnl_old + dpm[3]));
835  Real cosphi = std::cos(phi(intnl_old + dpm[3]));
836  Real coh = cohesion(intnl_old + dpm[3]);
837  cohcos = coh * cosphi;
838 
839  Real alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0)) * sinphi;
840  Real deriv_phi;
841  Real dalpha;
842  const Real beta = eigvals[2] + eigvals[0];
843  Real deriv_coh;
844 
845  Real residual =
846  alpha * dpm[3] + eigvals[2] - eigvals[0] + beta * sinphi - 2.0 * cohcos; // this is 2*yf[3]
847  Real jacobian;
848 
849  const Real f_tol2 = Utility::pow<2>(_f_tol);
850  unsigned int iter = 0;
851  do
852  {
853  deriv_phi = dphi(intnl_old + dpm[3]);
854  dalpha = -(n[3](2) + n[3](0)) * cosphi * deriv_phi;
855  deriv_coh = dcohesion(intnl_old + dpm[3]);
856  jacobian = alpha + dalpha * dpm[3] + beta * cosphi * deriv_phi - 2.0 * deriv_coh * cosphi +
857  2.0 * coh * sinphi * deriv_phi;
858 
859  dpm[3] -= residual / jacobian;
860  if (iter > _max_iters) // not converging
861  {
862  nr_converged = false;
863  return false;
864  }
865 
866  sinphi = std::sin(phi(intnl_old + dpm[3]));
867  cosphi = std::cos(phi(intnl_old + dpm[3]));
868  coh = cohesion(intnl_old + dpm[3]);
869  cohcos = coh * cosphi;
870  alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0)) * sinphi;
871  residual = alpha * dpm[3] + eigvals[2] - eigvals[0] + beta * sinphi - 2.0 * cohcos;
872  iter++;
873  } while (residual * residual > f_tol2);
874 
875  // so the NR process converged, but we must
876  // check Kuhn-Tucker
877  nr_converged = true;
878  if (dpm[3] < -ep_plastic_tolerance)
879  // Kuhn-Tucker failure
880  return false;
881 
882  for (unsigned i = 0; i < 3; ++i)
883  returned_stress(i, i) = eigvals[i] - dpm[3] * n[3](i);
885  returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
886 
887  // by construction abs(yf[3]) = abs(residual/2) < _f_tol/2
888  if (yf[0] > _f_tol || yf[1] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol || yf[5] > _f_tol)
889  // Kuhn-Tucker failure
890  return false;
891 
892  // success!
893  dpm[0] = dpm[1] = dpm[2] = dpm[4] = dpm[5] = 0;
894  return true;
895 }
896 
897 bool
899  const std::vector<RealVectorValue> & n,
900  std::vector<Real> & dpm,
901  RankTwoTensor & returned_stress,
902  Real intnl_old,
903  Real & sinphi,
904  Real & cohcos,
905  Real initial_guess,
906  Real mag_E,
907  bool & nr_converged,
908  Real ep_plastic_tolerance,
909  std::vector<Real> & yf) const
910 {
911  // This returns to the Mohr-Coulomb edge
912  // with 000101 being active. This means that
913  // f3=f5=0. Denoting the returned stress by "a"
914  // (in principal stress space), this means that
915  // a0=a1 = (2Ccosphi + a2(1+sinphi))/(sinphi-1)
916  //
917  // Also, a = eigvals - dpm3*n3 - dpm5*n5,
918  // where the dpm are the plastic multipliers
919  // and the n are the flow directions.
920  //
921  // Hence we have 5 equations and 5 unknowns (a,
922  // dpm3 and dpm5). Eliminating all unknowns
923  // except for x = dpm3+dpm5 gives the residual below
924  // (I used mathematica)
925 
926  Real x = initial_guess;
927  sinphi = std::sin(phi(intnl_old + x));
928  Real cosphi = std::cos(phi(intnl_old + x));
929  Real coh = cohesion(intnl_old + x);
930  cohcos = coh * cosphi;
931 
932  const Real numer_const =
933  -n[3](2) * eigvals[0] - n[5](1) * eigvals[0] + n[5](2) * eigvals[0] + n[3](2) * eigvals[1] +
934  n[5](0) * eigvals[1] - n[5](2) * eigvals[1] - n[5](0) * eigvals[2] + n[5](1) * eigvals[2] +
935  n[3](0) * (-eigvals[1] + eigvals[2]) - n[3](1) * (-eigvals[0] + eigvals[2]);
936  const Real numer_coeff1 = 2 * (-n[3](0) + n[3](1) + n[5](0) - n[5](1));
937  const Real numer_coeff2 =
938  n[5](2) * (eigvals[0] - eigvals[1]) + n[3](2) * (-eigvals[0] + eigvals[1]) +
939  n[5](1) * (eigvals[0] + eigvals[2]) + (n[3](0) - n[5](0)) * (eigvals[1] + eigvals[2]) -
940  n[3](1) * (eigvals[0] + eigvals[2]);
941  Real numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
942  const Real denom_const = -n[3](2) * (n[5](0) - n[5](1)) - n[3](1) * (-n[5](0) + n[5](2)) +
943  n[3](0) * (-n[5](1) + n[5](2));
944  const Real denom_coeff = -n[3](2) * (n[5](0) - n[5](1)) - n[3](1) * (n[5](0) + n[5](2)) +
945  n[3](0) * (n[5](1) + n[5](2));
946  Real denom = denom_const + denom_coeff * sinphi;
947  Real residual = -x + numer / denom;
948 
949  Real deriv_phi;
950  Real deriv_coh;
951  Real jacobian;
952  const Real tol = Utility::pow<2>(_f_tol / (mag_E * 10.0));
953  unsigned int iter = 0;
954  do
955  {
956  do
957  {
958  deriv_phi = dphi(intnl_old + dpm[3]);
959  deriv_coh = dcohesion(intnl_old + dpm[3]);
960  jacobian = -1 +
961  (numer_coeff1 * deriv_coh * cosphi - numer_coeff1 * coh * sinphi * deriv_phi +
962  numer_coeff2 * cosphi * deriv_phi) /
963  denom -
964  numer * denom_coeff * cosphi * deriv_phi / denom / denom;
965  x -= residual / jacobian;
966  if (iter > _max_iters) // not converging
967  {
968  nr_converged = false;
969  return false;
970  }
971 
972  sinphi = std::sin(phi(intnl_old + x));
973  cosphi = std::cos(phi(intnl_old + x));
974  coh = cohesion(intnl_old + x);
975  cohcos = coh * cosphi;
976  numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
977  denom = denom_const + denom_coeff * sinphi;
978  residual = -x + numer / denom;
979  iter++;
980  } while (residual * residual > tol);
981 
982  // now must ensure that yf[3] and yf[5] are both "zero"
983  const Real dpm3minusdpm5 =
984  (2.0 * (eigvals[0] - eigvals[1]) + x * (n[3](1) - n[3](0) + n[5](1) - n[5](0))) /
985  (n[3](0) - n[3](1) + n[5](1) - n[5](0));
986  dpm[3] = (x + dpm3minusdpm5) / 2.0;
987  dpm[5] = (x - dpm3minusdpm5) / 2.0;
988 
989  for (unsigned i = 0; i < 3; ++i)
990  returned_stress(i, i) = eigvals[i] - dpm[3] * n[3](i) - dpm[5] * n[5](i);
992  returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
993  } while (yf[3] * yf[3] > _f_tol * _f_tol && yf[5] * yf[5] > _f_tol * _f_tol);
994 
995  // so the NR process converged, but we must
996  // check Kuhn-Tucker
997  nr_converged = true;
998 
999  if (dpm[3] < -ep_plastic_tolerance || dpm[5] < -ep_plastic_tolerance)
1000  // Kuhn-Tucker failure. This is a potential weak-point: if the user has set _f_tol quite
1001  // large, and ep_plastic_tolerance quite small, the above NR process will quickly converge, but
1002  // the solution may be wrong and violate Kuhn-Tucker.
1003  return false;
1004 
1005  if (yf[0] > _f_tol || yf[1] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol)
1006  // Kuhn-Tucker failure
1007  return false;
1008 
1009  // success
1010  dpm[0] = dpm[1] = dpm[2] = dpm[4] = 0;
1011  return true;
1012 }
1013 
1014 bool
1016  const std::vector<RealVectorValue> & n,
1017  std::vector<Real> & dpm,
1018  RankTwoTensor & returned_stress,
1019  Real intnl_old,
1020  Real & sinphi,
1021  Real & cohcos,
1022  Real initial_guess,
1023  Real mag_E,
1024  bool & nr_converged,
1025  Real ep_plastic_tolerance,
1026  std::vector<Real> & yf) const
1027 {
1028  // This returns to the Mohr-Coulomb edge
1029  // with 010100 being active. This means that
1030  // f1=f3=0. Denoting the returned stress by "a"
1031  // (in principal stress space), this means that
1032  // a1=a2 = (2Ccosphi + a0(1-sinphi))/(sinphi+1)
1033  //
1034  // Also, a = eigvals - dpm1*n1 - dpm3*n3,
1035  // where the dpm are the plastic multipliers
1036  // and the n are the flow directions.
1037  //
1038  // Hence we have 5 equations and 5 unknowns (a,
1039  // dpm1 and dpm3). Eliminating all unknowns
1040  // except for x = dpm1+dpm3 gives the residual below
1041  // (I used mathematica)
1042 
1043  Real x = initial_guess;
1044  sinphi = std::sin(phi(intnl_old + x));
1045  Real cosphi = std::cos(phi(intnl_old + x));
1046  Real coh = cohesion(intnl_old + x);
1047  cohcos = coh * cosphi;
1048 
1049  const Real numer_const = -n[1](2) * eigvals[0] - n[3](1) * eigvals[0] + n[3](2) * eigvals[0] -
1050  n[1](0) * eigvals[1] + n[1](2) * eigvals[1] + n[3](0) * eigvals[1] -
1051  n[3](2) * eigvals[1] + n[1](0) * eigvals[2] - n[3](0) * eigvals[2] +
1052  n[3](1) * eigvals[2] - n[1](1) * (-eigvals[0] + eigvals[2]);
1053  const Real numer_coeff1 = 2 * (n[1](1) - n[1](2) - n[3](1) + n[3](2));
1054  const Real numer_coeff2 =
1055  n[3](2) * (-eigvals[0] - eigvals[1]) + n[1](2) * (eigvals[0] + eigvals[1]) +
1056  n[3](1) * (eigvals[0] + eigvals[2]) + (n[1](0) - n[3](0)) * (eigvals[1] - eigvals[2]) -
1057  n[1](1) * (eigvals[0] + eigvals[2]);
1058  Real numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
1059  const Real denom_const = -n[1](0) * (n[3](1) - n[3](2)) + n[1](2) * (-n[3](0) + n[3](1)) +
1060  n[1](1) * (-n[3](2) + n[3](0));
1061  const Real denom_coeff =
1062  n[1](0) * (n[3](1) - n[3](2)) + n[1](2) * (n[3](0) + n[3](1)) - n[1](1) * (n[3](0) + n[3](2));
1063  Real denom = denom_const + denom_coeff * sinphi;
1064  Real residual = -x + numer / denom;
1065 
1066  Real deriv_phi;
1067  Real deriv_coh;
1068  Real jacobian;
1069  const Real tol = Utility::pow<2>(_f_tol / (mag_E * 10.0));
1070  unsigned int iter = 0;
1071  do
1072  {
1073  do
1074  {
1075  deriv_phi = dphi(intnl_old + dpm[3]);
1076  deriv_coh = dcohesion(intnl_old + dpm[3]);
1077  jacobian = -1 +
1078  (numer_coeff1 * deriv_coh * cosphi - numer_coeff1 * coh * sinphi * deriv_phi +
1079  numer_coeff2 * cosphi * deriv_phi) /
1080  denom -
1081  numer * denom_coeff * cosphi * deriv_phi / denom / denom;
1082  x -= residual / jacobian;
1083  if (iter > _max_iters) // not converging
1084  {
1085  nr_converged = false;
1086  return false;
1087  }
1088 
1089  sinphi = std::sin(phi(intnl_old + x));
1090  cosphi = std::cos(phi(intnl_old + x));
1091  coh = cohesion(intnl_old + x);
1092  cohcos = coh * cosphi;
1093  numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
1094  denom = denom_const + denom_coeff * sinphi;
1095  residual = -x + numer / denom;
1096  iter++;
1097  } while (residual * residual > tol);
1098 
1099  // now must ensure that yf[1] and yf[3] are both "zero"
1100  Real dpm1minusdpm3 =
1101  (2 * (eigvals[1] - eigvals[2]) + x * (n[1](2) - n[1](1) + n[3](2) - n[3](1))) /
1102  (n[1](1) - n[1](2) + n[3](2) - n[3](1));
1103  dpm[1] = (x + dpm1minusdpm3) / 2.0;
1104  dpm[3] = (x - dpm1minusdpm3) / 2.0;
1105 
1106  for (unsigned i = 0; i < 3; ++i)
1107  returned_stress(i, i) = eigvals[i] - dpm[1] * n[1](i) - dpm[3] * n[3](i);
1109  returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
1110  } while (yf[1] * yf[1] > _f_tol * _f_tol && yf[3] * yf[3] > _f_tol * _f_tol);
1111 
1112  // so the NR process converged, but we must
1113  // check Kuhn-Tucker
1114  nr_converged = true;
1115 
1116  if (dpm[1] < -ep_plastic_tolerance || dpm[3] < -ep_plastic_tolerance)
1117  // Kuhn-Tucker failure. This is a potential weak-point: if the user has set _f_tol quite
1118  // large, and ep_plastic_tolerance quite small, the above NR process will quickly converge, but
1119  // the solution may be wrong and violate Kuhn-Tucker
1120  return false;
1121 
1122  if (yf[0] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol || yf[5] > _f_tol)
1123  // Kuhn-Tucker failure
1124  return false;
1125 
1126  // success
1127  dpm[0] = dpm[2] = dpm[4] = dpm[5] = 0;
1128  return true;
1129 }
1130 
1131 bool
1133  const std::vector<Real> & dpm,
1134  Real ep_plastic_tolerance) const
1135 {
1136  mooseAssert(
1137  yf.size() == 6,
1138  "TensorMechanicsPlasticMohrCoulomb::KuhnTuckerOK requires yf to be size 6, but your is "
1139  << yf.size());
1140  mooseAssert(
1141  dpm.size() == 6,
1142  "TensorMechanicsPlasticMohrCoulomb::KuhnTuckerOK requires dpm to be size 6, but your is "
1143  << dpm.size());
1144  for (unsigned i = 0; i < 6; ++i)
1145  if (!TensorMechanicsPlasticModel::KuhnTuckerSingleSurface(yf[i], dpm[i], ep_plastic_tolerance))
1146  return false;
1147  return true;
1148 }
1149 
1150 bool
1152 {
1153  return _use_custom_returnMap;
1154 }
virtual void dyieldFunction_dstressV(const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &df_dstress) const override
The derivative of yield functions with respect to stress.
const unsigned int _max_iters
Maximum Newton-Raphison iterations in the custom returnMap algorithm.
virtual bool returnMap(const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const override
Performs a custom return-map.
bool returnPlane(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
Tries to return-map to the MC plane using the n[3] direction The return value is true if the internal...
virtual bool useCustomReturnMap() const override
Returns false. You will want to override this in your derived class if you write a custom returnMap f...
bool returnEdge010100(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, Real mag_E, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
Tries to return-map to the MC edge using the n[1] and n[3] directions The return value is true if the...
virtual Real phi(const Real internal_param) const
phi as a function of residual value, rate, and internal_param
virtual std::string modelName() const =0
virtual void flowPotentialV(const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &r) const override
The flow potentials.
virtual void dflowPotential_dstressV(const RankTwoTensor &stress, Real intnl, std::vector< RankFourTensor > &dr_dstress) const override
The derivative of the flow potential with respect to stress.
InputParameters validParams< TensorMechanicsPlasticMohrCoulombMulti >()
virtual Real derivative(Real intnl) const
virtual bool returnMap(const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const
Performs a custom return-map.
virtual void activeConstraints(const std::vector< Real > &f, const RankTwoTensor &stress, Real intnl, const RankFourTensor &Eijkl, std::vector< bool > &act, RankTwoTensor &returned_stress) const override
The active yield surfaces, given a vector of yield functions.
const TensorMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
const TensorMechanicsHardeningModel & _phi
Hardening model for phi.
virtual Real dphi(const Real internal_param) const
d(phi)/d(internal_param) as a function of residual value, rate, and internal_param ...
const TensorMechanicsHardeningModel & _psi
Hardening model for psi.
virtual void yieldFunctionV(const RankTwoTensor &stress, Real intnl, std::vector< Real > &f) const override
Calculates the yield functions.
virtual Real psi(const Real internal_param) const
psi as a function of residual value, rate, and internal_param
bool KuhnTuckerSingleSurface(Real yf, Real dpm, Real dpm_tol) const
Returns true if the Kuhn-Tucker conditions for the single surface are satisfied.
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...
bool KuhnTuckerOK(const std::vector< Real > &yf, const std::vector< Real > &dpm, Real ep_plastic_tolerance) const
Returns true if the Kuhn-Tucker conditions are satisfied.
virtual void dyieldFunction_dintnlV(const RankTwoTensor &stress, Real intnl, std::vector< Real > &df_dintnl) const override
The derivative of yield functions with respect to the internal parameter.
InputParameters validParams< TensorMechanicsPlasticModel >()
virtual Real dcohesion(const Real internal_param) const
d(cohesion)/d(internal_param) as a function of residual value, rate, and internal_param ...
void yieldFunctionEigvals(Real e0, Real e1, Real e2, Real sinphi, Real cohcos, std::vector< Real > &f) const
Calculates the yield functions given the eigenvalues of stress.
bool doReturnMap(const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const
See doco for returnMap function.
static const double tol
Definition: XFEMFuncs.h:26
void df_dsig(const RankTwoTensor &stress, Real sin_angle, std::vector< RankTwoTensor > &df) const
this is exactly dyieldFunction_dstress, or flowPotential, depending on whether sin_angle = sin(phi)...
virtual Real dpsi(const Real internal_param) const
d(psi)/d(internal_param) as a function of residual value, rate, and internal_param ...
virtual Real value(Real intnl) const
const Real _f_tol
Tolerance on yield function.
TensorMechanicsPlasticMohrCoulombMulti(const InputParameters &parameters)
virtual unsigned int numberSurfaces() const override
The number of yield surfaces for this plasticity model.
virtual Real cohesion(const Real internal_param) const
cohesion as a function of residual value, rate, and internal_param
Plastic Model base class The virtual functions written below must be over-ridden in derived classes t...
const bool _use_custom_returnMap
Whether to use the custom return-map algorithm.
virtual void dflowPotential_dintnlV(const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &dr_dintnl) const override
The derivative of the flow potential with respect to the internal parameter.
bool returnTip(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
Tries to return-map to the MC tip using the THREE directions given in n, and THREE dpm values are ret...
void perturbStress(const RankTwoTensor &stress, std::vector< Real > &eigvals, std::vector< RankTwoTensor > &deigvals) const
perturbs the stress tensor in the case of almost-equal eigenvalues.
bool returnEdge000101(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, Real mag_E, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
Tries to return-map to the MC edge using the n[4] and n[6] directions The return value is true if the...