www.mooseframework.org
TensorMechanicsPlasticTensileMulti.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 
12 template <>
13 InputParameters
15 {
16  InputParameters params = validParams<TensorMechanicsPlasticModel>();
17  params.addClassDescription("Associative tensile plasticity with hardening/softening");
18  params.addRequiredParam<UserObjectName>(
19  "tensile_strength",
20  "A TensorMechanicsHardening UserObject that defines hardening of the tensile strength");
21  params.addParam<Real>("shift",
22  "Yield surface is shifted by this amount to avoid problems with "
23  "defining derivatives when eigenvalues are equal. If this is "
24  "larger than f_tol, a warning will be issued. Default = f_tol.");
25  params.addParam<unsigned int>("max_iterations",
26  10,
27  "Maximum number of Newton-Raphson iterations "
28  "allowed in the custom return-map algorithm. "
29  " For highly nonlinear hardening this may "
30  "need to be higher than 10.");
31  params.addParam<bool>("use_custom_returnMap",
32  true,
33  "Whether to use the custom returnMap "
34  "algorithm. Set to true if you are using "
35  "isotropic elasticity.");
36  params.addParam<bool>("use_custom_cto",
37  true,
38  "Whether to use the custom consistent tangent "
39  "operator computations. Set to true if you are "
40  "using isotropic elasticity.");
41  return params;
42 }
43 
45  const InputParameters & parameters)
46  : TensorMechanicsPlasticModel(parameters),
47  _strength(getUserObject<TensorMechanicsHardeningModel>("tensile_strength")),
48  _max_iters(getParam<unsigned int>("max_iterations")),
49  _shift(parameters.isParamValid("shift") ? getParam<Real>("shift") : _f_tol),
50  _use_custom_returnMap(getParam<bool>("use_custom_returnMap")),
51  _use_custom_cto(getParam<bool>("use_custom_cto"))
52 {
53  if (_shift < 0)
54  mooseError("Value of 'shift' in TensorMechanicsPlasticTensileMulti must not be negative\n");
55  if (_shift > _f_tol)
56  _console << "WARNING: value of 'shift' in TensorMechanicsPlasticTensileMulti is probably set "
57  "too high\n";
58  if (LIBMESH_DIM != 3)
59  mooseError("TensorMechanicsPlasticTensileMulti is only defined for LIBMESH_DIM=3");
60  MooseRandom::seed(0);
61 }
62 
63 unsigned int
65 {
66  return 3;
67 }
68 
69 void
71  Real intnl,
72  std::vector<Real> & f) const
73 {
74  std::vector<Real> eigvals;
75  stress.symmetricEigenvalues(eigvals);
76  const Real str = tensile_strength(intnl);
77 
78  f.resize(3);
79  f[0] = eigvals[0] + _shift - str;
80  f[1] = eigvals[1] - str;
81  f[2] = eigvals[2] - _shift - str;
82 }
83 
84 void
86  const RankTwoTensor & stress, Real /*intnl*/, std::vector<RankTwoTensor> & df_dstress) const
87 {
88  std::vector<Real> eigvals;
89  stress.dsymmetricEigenvalues(eigvals, df_dstress);
90 
91  if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
92  {
93  Real small_perturbation;
94  RankTwoTensor shifted_stress = stress;
95  while (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
96  {
97  for (unsigned i = 0; i < 3; ++i)
98  for (unsigned j = 0; j <= i; ++j)
99  {
100  small_perturbation = 0.1 * _shift * 2 * (MooseRandom::rand() - 0.5);
101  shifted_stress(i, j) += small_perturbation;
102  shifted_stress(j, i) += small_perturbation;
103  }
104  shifted_stress.dsymmetricEigenvalues(eigvals, df_dstress);
105  }
106  }
107 }
108 
109 void
111  Real intnl,
112  std::vector<Real> & df_dintnl) const
113 {
114  df_dintnl.assign(3, -dtensile_strength(intnl));
115 }
116 
117 void
119  Real intnl,
120  std::vector<RankTwoTensor> & r) const
121 {
122  // This plasticity is associative so
123  dyieldFunction_dstressV(stress, intnl, r);
124 }
125 
126 void
128  const RankTwoTensor & stress, Real /*intnl*/, std::vector<RankFourTensor> & dr_dstress) const
129 {
130  stress.d2symmetricEigenvalues(dr_dstress);
131 }
132 
133 void
135  const RankTwoTensor & /*stress*/, Real /*intnl*/, std::vector<RankTwoTensor> & dr_dintnl) const
136 {
137  dr_dintnl.assign(3, RankTwoTensor());
138 }
139 
140 Real
142 {
143  return _strength.value(internal_param);
144 }
145 
146 Real
148 {
149  return _strength.derivative(internal_param);
150 }
151 
152 void
154  const RankTwoTensor & stress,
155  Real intnl,
156  const RankFourTensor & Eijkl,
157  std::vector<bool> & act,
158  RankTwoTensor & returned_stress) const
159 {
160  act.assign(3, false);
161 
162  if (f[0] <= _f_tol && f[1] <= _f_tol && f[2] <= _f_tol)
163  {
164  returned_stress = stress;
165  return;
166  }
167 
168  Real returned_intnl;
169  std::vector<Real> dpm(3);
170  RankTwoTensor delta_dp;
171  std::vector<Real> yf(3);
172  bool trial_stress_inadmissible;
173  doReturnMap(stress,
174  intnl,
175  Eijkl,
176  0.0,
177  returned_stress,
178  returned_intnl,
179  dpm,
180  delta_dp,
181  yf,
182  trial_stress_inadmissible);
183 
184  for (unsigned i = 0; i < 3; ++i)
185  act[i] = (dpm[i] > 0);
186 }
187 
188 Real
189 TensorMechanicsPlasticTensileMulti::dot(const std::vector<Real> & a,
190  const std::vector<Real> & b) const
191 {
192  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
193 }
194 
195 Real
196 TensorMechanicsPlasticTensileMulti::triple(const std::vector<Real> & a,
197  const std::vector<Real> & b,
198  const std::vector<Real> & c) const
199 {
200  return a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) +
201  a[2] * (b[0] * c[1] - b[1] * c[0]);
202 }
203 
204 std::string
206 {
207  return "TensileMulti";
208 }
209 
210 bool
211 TensorMechanicsPlasticTensileMulti::returnMap(const RankTwoTensor & trial_stress,
212  Real intnl_old,
213  const RankFourTensor & E_ijkl,
214  Real ep_plastic_tolerance,
215  RankTwoTensor & returned_stress,
216  Real & returned_intnl,
217  std::vector<Real> & dpm,
218  RankTwoTensor & delta_dp,
219  std::vector<Real> & yf,
220  bool & trial_stress_inadmissible) const
221 {
223  return TensorMechanicsPlasticModel::returnMap(trial_stress,
224  intnl_old,
225  E_ijkl,
226  ep_plastic_tolerance,
227  returned_stress,
228  returned_intnl,
229  dpm,
230  delta_dp,
231  yf,
232  trial_stress_inadmissible);
233 
234  return doReturnMap(trial_stress,
235  intnl_old,
236  E_ijkl,
237  ep_plastic_tolerance,
238  returned_stress,
239  returned_intnl,
240  dpm,
241  delta_dp,
242  yf,
243  trial_stress_inadmissible);
244 }
245 
246 bool
247 TensorMechanicsPlasticTensileMulti::doReturnMap(const RankTwoTensor & trial_stress,
248  Real intnl_old,
249  const RankFourTensor & E_ijkl,
250  Real ep_plastic_tolerance,
251  RankTwoTensor & returned_stress,
252  Real & returned_intnl,
253  std::vector<Real> & dpm,
254  RankTwoTensor & delta_dp,
255  std::vector<Real> & yf,
256  bool & trial_stress_inadmissible) const
257 {
258  mooseAssert(dpm.size() == 3,
259  "TensorMechanicsPlasticTensileMulti size of dpm should be 3 but it is "
260  << dpm.size());
261 
262  std::vector<Real> eigvals;
263  RankTwoTensor eigvecs;
264  trial_stress.symmetricEigenvaluesEigenvectors(eigvals, eigvecs);
265  eigvals[0] += _shift;
266  eigvals[2] -= _shift;
267 
268  Real str = tensile_strength(intnl_old);
269 
270  yf.resize(3);
271  yf[0] = eigvals[0] - str;
272  yf[1] = eigvals[1] - str;
273  yf[2] = eigvals[2] - str;
274 
275  if (yf[0] <= _f_tol && yf[1] <= _f_tol && yf[2] <= _f_tol)
276  {
277  // purely elastic (trial_stress, intnl_old)
278  trial_stress_inadmissible = false;
279  return true;
280  }
281 
282  trial_stress_inadmissible = true;
283  delta_dp.zero();
284  returned_stress.zero();
285 
286  // In the following i often assume that E_ijkl is
287  // for an isotropic situation. This reduces FLOPS
288  // substantially which is important since the returnMap
289  // is potentially the most compute-intensive function
290  // of a simulation.
291  // In many comments i write the general expression, and
292  // i hope that might guide future coders if they are
293  // generalising to a non-istropic E_ijkl
294 
295  // n[alpha] = E_ijkl*r[alpha]_kl expressed in principal stress space
296  // (alpha = 0, 1, 2, corresponding to the three surfaces)
297  // Note that in principal stress space, the flow
298  // directions are, expressed in 'vector' form,
299  // r[0] = (1,0,0), r[1] = (0,1,0), r[2] = (0,0,1).
300  // Similar for _n:
301  // so _n[0] = E_ij00*r[0], _n[1] = E_ij11*r[1], _n[2] = E_ij22*r[2]
302  // In the following I assume that the E_ijkl is
303  // for an isotropic situation.
304  // In the anisotropic situation, we couldn't express
305  // the flow directions as vectors in the same principal
306  // stress space as the stress: they'd be full rank-2 tensors
307  std::vector<RealVectorValue> n(3);
308  n[0](0) = E_ijkl(0, 0, 0, 0);
309  n[0](1) = E_ijkl(1, 1, 0, 0);
310  n[0](2) = E_ijkl(2, 2, 0, 0);
311  n[1](0) = E_ijkl(0, 0, 1, 1);
312  n[1](1) = E_ijkl(1, 1, 1, 1);
313  n[1](2) = E_ijkl(2, 2, 1, 1);
314  n[2](0) = E_ijkl(0, 0, 2, 2);
315  n[2](1) = E_ijkl(1, 1, 2, 2);
316  n[2](2) = E_ijkl(2, 2, 2, 2);
317 
318  // With non-zero Poisson's ratio and hardening
319  // it is not computationally cheap to know whether
320  // the trial stress will return to the tip, edge,
321  // or plane. The following is correct for zero
322  // Poisson's ratio and no hardening, and at least
323  // gives a not-completely-stupid guess in the
324  // more general case.
325  // trial_order[0] = type of return to try first
326  // trial_order[1] = type of return to try second
327  // trial_order[2] = type of return to try third
328  const unsigned int number_of_return_paths = 3;
329  std::vector<int> trial_order(number_of_return_paths);
330  if (yf[0] > _f_tol) // all the yield functions are positive, since eigvals are ordered eigvals[0]
331  // <= eigvals[1] <= eigvals[2]
332  {
333  trial_order[0] = tip;
334  trial_order[1] = edge;
335  trial_order[2] = plane;
336  }
337  else if (yf[1] > _f_tol) // two yield functions are positive
338  {
339  trial_order[0] = edge;
340  trial_order[1] = tip;
341  trial_order[2] = plane;
342  }
343  else
344  {
345  trial_order[0] = plane;
346  trial_order[1] = edge;
347  trial_order[2] = tip;
348  }
349 
350  unsigned trial;
351  bool nr_converged = false;
352  for (trial = 0; trial < number_of_return_paths; ++trial)
353  {
354  switch (trial_order[trial])
355  {
356  case tip:
357  nr_converged = returnTip(eigvals, n, dpm, returned_stress, intnl_old, 0);
358  break;
359  case edge:
360  nr_converged = returnEdge(eigvals, n, dpm, returned_stress, intnl_old, 0);
361  break;
362  case plane:
363  nr_converged = returnPlane(eigvals, n, dpm, returned_stress, intnl_old, 0);
364  break;
365  }
366 
367  str = tensile_strength(intnl_old + dpm[0] + dpm[1] + dpm[2]);
368 
369  if (nr_converged && KuhnTuckerOK(returned_stress, dpm, str, ep_plastic_tolerance))
370  break;
371  }
372 
373  if (trial == number_of_return_paths)
374  {
375  Moose::err << "Trial stress = \n";
376  trial_stress.print(Moose::err);
377  Moose::err << "Internal parameter = " << intnl_old << "\n";
378  mooseError("TensorMechanicsPlasticTensileMulti: FAILURE! You probably need to implement a "
379  "line search\n");
380  // failure - must place yield function values at trial stress into yf
381  str = tensile_strength(intnl_old);
382  yf[0] = eigvals[0] - str;
383  yf[1] = eigvals[1] - str;
384  yf[2] = eigvals[2] - str;
385  return false;
386  }
387 
388  // success
389 
390  returned_intnl = intnl_old;
391  for (unsigned i = 0; i < 3; ++i)
392  {
393  yf[i] = returned_stress(i, i) - str;
394  delta_dp(i, i) = dpm[i];
395  returned_intnl += dpm[i];
396  }
397  returned_stress = eigvecs * returned_stress * (eigvecs.transpose());
398  delta_dp = eigvecs * delta_dp * (eigvecs.transpose());
399  return true;
400 }
401 
402 bool
403 TensorMechanicsPlasticTensileMulti::returnTip(const std::vector<Real> & eigvals,
404  const std::vector<RealVectorValue> & n,
405  std::vector<Real> & dpm,
406  RankTwoTensor & returned_stress,
407  Real intnl_old,
408  Real initial_guess) const
409 {
410  // The returned point is defined by f0=f1=f2=0.
411  // that is, returned_stress = diag(str, str, str), where
412  // str = tensile_strength(intnl),
413  // where intnl = intnl_old + dpm[0] + dpm[1] + dpm[2]
414  // The 3 plastic multipliers, dpm, are defiend by the normality condition
415  // eigvals - str = dpm[0]*n[0] + dpm[1]*n[1] + dpm[2]*n[2]
416  // (Kuhn-Tucker demands that all dpm are non-negative, but we leave
417  // that checking for later.)
418  // This is a vector equation with solution (A):
419  // dpm[0] = triple(eigvals - str, n[1], n[2])/trip;
420  // dpm[1] = triple(eigvals - str, n[2], n[0])/trip;
421  // dpm[2] = triple(eigvals - str, n[0], n[1])/trip;
422  // where trip = triple(n[0], n[1], n[2]).
423  // By adding the three components of that solution together
424  // we can get an equation for x = dpm[0] + dpm[1] + dpm[2],
425  // and then our Newton-Raphson only involves one variable (x).
426  // In the following, i specialise to the isotropic situation.
427 
428  Real x = initial_guess;
429  const Real denom = (n[0](0) - n[0](1)) * (n[0](0) + 2 * n[0](1));
430  Real str = tensile_strength(intnl_old + x);
431 
432  if (_strength.modelName().compare("Constant") != 0)
433  {
434  // Finding x is expensive. Therefore
435  // although x!=0 for Constant Hardening, solution (A)
436  // demonstrates that we don't
437  // actually need to know x to find the dpm for
438  // Constant Hardening.
439  //
440  // However, for nontrivial Hardening, the following
441  // is necessary
442  const Real eig = eigvals[0] + eigvals[1] + eigvals[2];
443  const Real bul = (n[0](0) + 2 * n[0](1));
444 
445  // and finally, the equation we want to solve is:
446  // bul*x - eig + 3*str = 0
447  // where str=tensile_strength(intnl_old + x)
448  // and x = dpm[0] + dpm[1] + dpm[2]
449  // (Note this has units of stress, so using _f_tol as a convergence check is reasonable.)
450  // Use Netwon-Raphson with initial guess x = 0
451  Real residual = bul * x - eig + 3 * str;
452  Real jacobian;
453  unsigned int iter = 0;
454  do
455  {
456  jacobian = bul + 3 * dtensile_strength(intnl_old + x);
457  x += -residual / jacobian;
458  if (iter > _max_iters) // not converging
459  return false;
460  str = tensile_strength(intnl_old + x);
461  residual = bul * x - eig + 3 * str;
462  iter++;
463  } while (residual * residual > _f_tol * _f_tol);
464  }
465 
466  // The following is the solution (A) written above
467  // (dpm[0] = triple(eigvals - str, n[1], n[2])/trip, etc)
468  // in the isotropic situation
469  dpm[0] = (n[0](0) * (eigvals[0] - str) + n[0](1) * (eigvals[0] - eigvals[1] - eigvals[2] + str)) /
470  denom;
471  dpm[1] = (n[0](0) * (eigvals[1] - str) + n[0](1) * (eigvals[1] - eigvals[2] - eigvals[0] + str)) /
472  denom;
473  dpm[2] = (n[0](0) * (eigvals[2] - str) + n[0](1) * (eigvals[2] - eigvals[0] - eigvals[1] + str)) /
474  denom;
475  returned_stress(0, 0) = returned_stress(1, 1) = returned_stress(2, 2) = str;
476  return true;
477 }
478 
479 bool
480 TensorMechanicsPlasticTensileMulti::returnEdge(const std::vector<Real> & eigvals,
481  const std::vector<RealVectorValue> & n,
482  std::vector<Real> & dpm,
483  RankTwoTensor & returned_stress,
484  Real intnl_old,
485  Real initial_guess) const
486 {
487  // work out the point to which we would return, "a". It is defined by
488  // f1 = 0 = f2, and the normality condition:
489  // (eigvals - a).(n1 x n2) = 0,
490  // where eigvals is the starting position
491  // (it is a vector in principal stress space).
492  // To get f1=0=f2, we need a = (a0, str, str), and a0 is found
493  // by expanding the normality condition to yield:
494  // a0 = (-str*n1xn2[1] - str*n1xn2[2] + edotn1xn2)/n1xn2[0];
495  // where edotn1xn2 = eigvals.(n1 x n2)
496  //
497  // We need to find the plastic multipliers, dpm, defined by
498  // eigvals - a = dpm[1]*n1 + dpm[2]*n2
499  // For the isotropic case, and defining eminusa = eigvals - a,
500  // the solution is easy:
501  // dpm[0] = 0;
502  // dpm[1] = (eminusa[1] - eminusa[0])/(n[1][1] - n[1][0]);
503  // dpm[2] = (eminusa[2] - eminusa[0])/(n[2][2] - n[2][0]);
504  //
505  // Now specialise to the isotropic case. Define
506  // x = dpm[1] + dpm[2] = (eigvals[1] + eigvals[2] - 2*str)/(n[0][0] + n[0][1])
507  // Notice that the RHS is a function of x, so we solve using
508  // Newton-Raphson starting with x=initial_guess
509  Real x = initial_guess;
510  const Real denom = n[0](0) + n[0](1);
511  Real str = tensile_strength(intnl_old + x);
512 
513  if (_strength.modelName().compare("Constant") != 0)
514  {
515  // Finding x is expensive. Therefore
516  // although x!=0 for Constant Hardening, solution
517  // for dpm above demonstrates that we don't
518  // actually need to know x to find the dpm for
519  // Constant Hardening.
520  //
521  // However, for nontrivial Hardening, the following
522  // is necessary
523  const Real eig = eigvals[1] + eigvals[2];
524  Real residual = denom * x - eig + 2 * str;
525  Real jacobian;
526  unsigned int iter = 0;
527  do
528  {
529  jacobian = denom + 2 * dtensile_strength(intnl_old + x);
530  x += -residual / jacobian;
531  if (iter > _max_iters) // not converging
532  return false;
533  str = tensile_strength(intnl_old + x);
534  residual = denom * x - eig + 2 * str;
535  iter++;
536  } while (residual * residual > _f_tol * _f_tol);
537  }
538 
539  dpm[0] = 0;
540  dpm[1] = ((eigvals[1] * n[0](0) - eigvals[2] * n[0](1)) / (n[0](0) - n[0](1)) - str) / denom;
541  dpm[2] = ((eigvals[2] * n[0](0) - eigvals[1] * n[0](1)) / (n[0](0) - n[0](1)) - str) / denom;
542 
543  returned_stress(0, 0) = eigvals[0] - n[0](1) * (dpm[1] + dpm[2]);
544  returned_stress(1, 1) = returned_stress(2, 2) = str;
545  return true;
546 }
547 
548 bool
549 TensorMechanicsPlasticTensileMulti::returnPlane(const std::vector<Real> & eigvals,
550  const std::vector<RealVectorValue> & n,
551  std::vector<Real> & dpm,
552  RankTwoTensor & returned_stress,
553  Real intnl_old,
554  Real initial_guess) const
555 {
556  // the returned point, "a", is defined by f2=0 and
557  // a = p - dpm[2]*n2.
558  // This is a vector equation in
559  // principal stress space, and dpm[2] is the third
560  // plasticity multiplier (dpm[0]=0=dpm[1] for return
561  // to the plane) and "p" is the starting
562  // position (p=eigvals).
563  // (Kuhn-Tucker demands that dpm[2]>=0, but we leave checking
564  // that condition for later.)
565  // Since f2=0, we must have a[2]=tensile_strength,
566  // so we can just look at the [2] component of the
567  // equation, which yields
568  // n[2][2]*dpm[2] - eigvals[2] + str = 0
569  // For hardening, str=tensile_strength(intnl_old+dpm[2]),
570  // and we want to solve for dpm[2].
571  // Use Newton-Raphson with initial guess dpm[2] = initial_guess
572  dpm[2] = initial_guess;
573  Real residual = n[2](2) * dpm[2] - eigvals[2] + tensile_strength(intnl_old + dpm[2]);
574  Real jacobian;
575  unsigned int iter = 0;
576  do
577  {
578  jacobian = n[2](2) + dtensile_strength(intnl_old + dpm[2]);
579  dpm[2] += -residual / jacobian;
580  if (iter > _max_iters) // not converging
581  return false;
582  residual = n[2](2) * dpm[2] - eigvals[2] + tensile_strength(intnl_old + dpm[2]);
583  iter++;
584  } while (residual * residual > _f_tol * _f_tol);
585 
586  dpm[0] = 0;
587  dpm[1] = 0;
588  returned_stress(0, 0) = eigvals[0] - dpm[2] * n[2](0);
589  returned_stress(1, 1) = eigvals[1] - dpm[2] * n[2](1);
590  returned_stress(2, 2) = eigvals[2] - dpm[2] * n[2](2);
591  return true;
592 }
593 
594 bool
595 TensorMechanicsPlasticTensileMulti::KuhnTuckerOK(const RankTwoTensor & returned_diagonal_stress,
596  const std::vector<Real> & dpm,
597  Real str,
598  Real ep_plastic_tolerance) const
599 {
600  for (unsigned i = 0; i < 3; ++i)
602  returned_diagonal_stress(i, i) - str, dpm[i], ep_plastic_tolerance))
603  return false;
604  return true;
605 }
606 
607 RankFourTensor
609  const RankTwoTensor & trial_stress,
610  Real intnl_old,
611  const RankTwoTensor & stress,
612  Real intnl,
613  const RankFourTensor & E_ijkl,
614  const std::vector<Real> & cumulative_pm) const
615 {
616  if (!_use_custom_cto)
618  trial_stress, intnl_old, stress, intnl, E_ijkl, cumulative_pm);
619 
620  mooseAssert(cumulative_pm.size() == 3,
621  "TensorMechanicsPlasticTensileMulti size of cumulative_pm should be 3 but it is "
622  << cumulative_pm.size());
623 
624  if (cumulative_pm[2] <= 0) // All cumulative_pm are non-positive, so this is admissible
625  return E_ijkl;
626 
627  // Need the eigenvalues at the returned configuration
628  std::vector<Real> eigvals;
629  stress.symmetricEigenvalues(eigvals);
630 
631  // need to rotate to and from principal stress space
632  // using the eigenvectors of the trial configuration
633  // (not the returned configuration).
634  std::vector<Real> trial_eigvals;
635  RankTwoTensor trial_eigvecs;
636  trial_stress.symmetricEigenvaluesEigenvectors(trial_eigvals, trial_eigvecs);
637 
638  // The returnMap will have returned to the Tip, Edge or
639  // Plane. The consistentTangentOperator describes the
640  // change in stress for an arbitrary change in applied
641  // strain. I assume that the change in strain will not
642  // change the type of return (Tip remains Tip, Edge remains
643  // Edge, Plane remains Plane).
644  // I assume isotropic elasticity.
645  //
646  // The consistent tangent operator is a little different
647  // than cases where no rotation to principal stress space
648  // is made during the returnMap. Let S_ij be the stress
649  // in original coordinates, and s_ij be the stress in the
650  // principal stress coordinates, so that
651  // s_ij = diag(eigvals[0], eigvals[1], eigvals[2])
652  // We want dS_ij under an arbitrary change in strain (ep->ep+dep)
653  // dS = S(ep+dep) - S(ep)
654  // = R(ep+dep) s(ep+dep) R(ep+dep)^T - R(ep) s(ep) R(ep)^T
655  // Here R = the rotation to principal-stress space, ie
656  // R_ij = eigvecs[i][j] = i^th component of j^th eigenvector
657  // Expanding to first order in dep,
658  // dS = R(ep) (s(ep+dep) - s(ep)) R(ep)^T
659  // + dR/dep s(ep) R^T + R(ep) s(ep) dR^T/dep
660  // The first line is all that is usually calculated in the
661  // consistent tangent operator calculation, and is called
662  // cto below.
663  // The second line involves changes in the eigenvectors, and
664  // is called sec below.
665 
666  RankFourTensor cto;
667  const Real hard = dtensile_strength(intnl);
668  const Real la = E_ijkl(0, 0, 1, 1);
669  const Real mu = 0.5 * (E_ijkl(0, 0, 0, 0) - la);
670 
671  if (cumulative_pm[1] <= 0)
672  {
673  // only cumulative_pm[2] is positive, so this is return to the Plane
674  const Real denom = hard + la + 2 * mu;
675  const Real al = la * la / denom;
676  const Real be = la * (la + 2 * mu) / denom;
677  const Real ga = hard * (la + 2 * mu) / denom;
678  std::vector<Real> comps(9);
679  comps[0] = comps[4] = la + 2 * mu - al;
680  comps[1] = comps[3] = la - al;
681  comps[2] = comps[5] = comps[6] = comps[7] = la - be;
682  comps[8] = ga;
683  cto.fillFromInputVector(comps, RankFourTensor::principal);
684  }
685  else if (cumulative_pm[0] <= 0)
686  {
687  // both cumulative_pm[2] and cumulative_pm[1] are positive, so Edge
688  const Real denom = 2 * hard + 2 * la + 2 * mu;
689  const Real al = hard * 2 * la / denom;
690  const Real be = hard * (2 * la + 2 * mu) / denom;
691  std::vector<Real> comps(9);
692  comps[0] = la + 2 * mu - 2 * la * la / denom;
693  comps[1] = comps[2] = al;
694  comps[3] = comps[6] = al;
695  comps[4] = comps[5] = comps[7] = comps[8] = be;
696  cto.fillFromInputVector(comps, RankFourTensor::principal);
697  }
698  else
699  {
700  // all cumulative_pm are positive, so Tip
701  const Real denom = 3 * hard + 3 * la + 2 * mu;
702  std::vector<Real> comps(2);
703  comps[0] = hard * (3 * la + 2 * mu) / denom;
704  comps[1] = 0;
705  cto.fillFromInputVector(comps, RankFourTensor::symmetric_isotropic);
706  }
707 
708  cto.rotate(trial_eigvecs);
709 
710  // drdsig = change in eigenvectors under a small stress change
711  // drdsig(i,j,m,n) = dR(i,j)/dS_mn
712  // The formula below is fairly easily derived:
713  // S R = R s, so taking the variation
714  // dS R + S dR = dR s + R ds, and multiplying by R^T
715  // R^T dS R + R^T S dR = R^T dR s + ds .... (eqn 1)
716  // I demand that RR^T = 1 = R^T R, and also that
717  // (R+dR)(R+dR)^T = 1 = (R+dT)^T (R+dR), which means
718  // that dR = R*c, for some antisymmetric c, so Eqn1 reads
719  // R^T dS R + s c = c s + ds
720  // Grabbing the components of this gives ds/dS (already
721  // in RankTwoTensor), and c, which is:
722  // dR_ik/dS_mn = drdsig(i, k, m, n) = trial_eigvecs(m, b)*trial_eigvecs(n, k)*trial_eigvecs(i,
723  // b)/(trial_eigvals[k] - trial_eigvals[b]);
724  // (sum over b!=k).
725 
726  RankFourTensor drdsig;
727  for (unsigned k = 0; k < 3; ++k)
728  for (unsigned b = 0; b < 3; ++b)
729  {
730  if (b == k)
731  continue;
732  for (unsigned m = 0; m < 3; ++m)
733  for (unsigned n = 0; n < 3; ++n)
734  for (unsigned i = 0; i < 3; ++i)
735  drdsig(i, k, m, n) += trial_eigvecs(m, b) * trial_eigvecs(n, k) * trial_eigvecs(i, b) /
736  (trial_eigvals[k] - trial_eigvals[b]);
737  }
738 
739  // With diagla = diag(eigvals[0], eigvals[1], digvals[2])
740  // The following implements
741  // ans(i, j, a, b) += (drdsig(i, k, m, n)*trial_eigvecs(j, l)*diagla(k, l) + trial_eigvecs(i,
742  // k)*drdsig(j, l, m, n)*diagla(k, l))*E_ijkl(m, n, a, b);
743  // (sum over k, l, m and n)
744 
745  RankFourTensor ans;
746  for (unsigned i = 0; i < 3; ++i)
747  for (unsigned j = 0; j < 3; ++j)
748  for (unsigned a = 0; a < 3; ++a)
749  for (unsigned k = 0; k < 3; ++k)
750  for (unsigned m = 0; m < 3; ++m)
751  ans(i, j, a, a) += (drdsig(i, k, m, m) * trial_eigvecs(j, k) +
752  trial_eigvecs(i, k) * drdsig(j, k, m, m)) *
753  eigvals[k] * la; // E_ijkl(m, n, a, b) = la*(m==n)*(a==b);
754 
755  for (unsigned i = 0; i < 3; ++i)
756  for (unsigned j = 0; j < 3; ++j)
757  for (unsigned a = 0; a < 3; ++a)
758  for (unsigned b = 0; b < 3; ++b)
759  for (unsigned k = 0; k < 3; ++k)
760  {
761  ans(i, j, a, b) += (drdsig(i, k, a, b) * trial_eigvecs(j, k) +
762  trial_eigvecs(i, k) * drdsig(j, k, a, b)) *
763  eigvals[k] * mu; // E_ijkl(m, n, a, b) = mu*(m==a)*(n==b)
764  ans(i, j, a, b) += (drdsig(i, k, b, a) * trial_eigvecs(j, k) +
765  trial_eigvecs(i, k) * drdsig(j, k, b, a)) *
766  eigvals[k] * mu; // E_ijkl(m, n, a, b) = mu*(m==b)*(n==a)
767  }
768 
769  return cto + ans;
770 }
771 
772 bool
774 {
775  return _use_custom_returnMap;
776 }
777 
778 bool
780 {
781  return _use_custom_cto;
782 }
bool returnPlane(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real initial_guess) const
Tries to return-map to the Tensile plane The return value is true if the internal Newton-Raphson proc...
virtual Real tensile_strength(const Real internal_param) const
tensile strength as a function of residual value, rate, and internal_param
virtual void yieldFunctionV(const RankTwoTensor &stress, Real intnl, std::vector< Real > &f) const override
Calculates the yield functions.
virtual std::string modelName() const =0
virtual RankFourTensor consistentTangentOperator(const RankTwoTensor &trial_stress, Real intnl_old, const RankTwoTensor &stress, Real intnl, const RankFourTensor &E_ijkl, const std::vector< Real > &cumulative_pm) const
Calculates a custom consistent tangent operator.
bool returnTip(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real initial_guess) const
Tries to return-map to the Tensile tip.
virtual bool useCustomCTO() const override
Returns false. You will want to override this in your derived class if you write a custom consistent ...
virtual 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
Just like returnMap, but a protected interface that definitely uses the algorithm, since returnMap itself does not use the algorithm if _use_returnMap=false.
virtual Real dtensile_strength(const Real internal_param) const
d(tensile strength)/d(internal_param) as a function of residual value, rate, and internal_param ...
virtual RankFourTensor consistentTangentOperator(const RankTwoTensor &trial_stress, Real intnl_old, const RankTwoTensor &stress, Real intnl, const RankFourTensor &E_ijkl, const std::vector< Real > &cumulative_pm) const override
Calculates a custom consistent tangent operator.
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 std::string modelName() const override
const bool _use_custom_returnMap
Whether to use the custom return-map algorithm.
bool KuhnTuckerOK(const RankTwoTensor &returned_diagonal_stress, const std::vector< Real > &dpm, Real str, Real ep_plastic_tolerance) const
Returns true if the Kuhn-Tucker conditions are satisfied.
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...
Real dot(const std::vector< Real > &a, const std::vector< Real > &b) const
dot product of two 3-dimensional vectors
Real triple(const std::vector< Real > &a, const std::vector< Real > &b, const std::vector< Real > &c) const
triple product of three 3-dimensional vectors
bool KuhnTuckerSingleSurface(Real yf, Real dpm, Real dpm_tol) const
Returns true if the Kuhn-Tucker conditions for the single surface are satisfied.
const unsigned int _max_iters
maximum iterations allowed in the custom return-map algorithm
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< TensorMechanicsPlasticModel >()
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.
virtual Real value(Real intnl) const
InputParameters validParams< TensorMechanicsPlasticTensileMulti >()
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.
virtual unsigned int numberSurfaces() const override
The number of yield surfaces for this plasticity model.
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.
const Real _f_tol
Tolerance on yield function.
virtual bool useCustomReturnMap() const override
Returns false. You will want to override this in your derived class if you write a custom returnMap f...
virtual void flowPotentialV(const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &r) const override
The flow potentials.
const bool _use_custom_cto
Whether to use the custom consistent tangent operator calculation.
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.
Plastic Model base class The virtual functions written below must be over-ridden in derived classes t...
TensorMechanicsPlasticTensileMulti(const InputParameters &parameters)
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.
const TensorMechanicsHardeningModel & _strength
bool returnEdge(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real initial_guess) const
Tries to return-map to the Tensile edge.