www.mooseframework.org
CappedMohrCoulombStressUpdate.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
11 #include "libmesh/utility.h"
12 #include "ElasticityTensorTools.h"
13 
15 
18 {
20  params.addRequiredParam<UserObjectName>(
21  "tensile_strength",
22  "A SolidMechanicsHardening UserObject that defines hardening of the "
23  "tensile strength. In physical situations this is positive (and always "
24  "must be greater than negative compressive-strength.");
25  params.addRequiredParam<UserObjectName>(
26  "compressive_strength",
27  "A SolidMechanicsHardening UserObject that defines hardening of the "
28  "compressive strength. In physical situations this is positive.");
29  params.addRequiredParam<UserObjectName>(
30  "cohesion", "A SolidMechanicsHardening UserObject that defines hardening of the cohesion");
31  params.addRequiredParam<UserObjectName>("friction_angle",
32  "A SolidMechanicsHardening UserObject "
33  "that defines hardening of the "
34  "friction angle (in radians)");
35  params.addRequiredParam<UserObjectName>(
36  "dilation_angle",
37  "A SolidMechanicsHardening UserObject that defines hardening of the "
38  "dilation angle (in radians). Unless you are quite confident, this should "
39  "be set positive and not greater than the friction angle.");
40  params.addParam<bool>("perfect_guess",
41  true,
42  "Provide a guess to the Newton-Raphson procedure "
43  "that is the result from perfect plasticity. With "
44  "severe hardening/softening this may be "
45  "suboptimal.");
46  params.addClassDescription("Nonassociative, smoothed, Mohr-Coulomb plasticity capped with "
47  "tensile (Rankine) and compressive caps, with hardening/softening");
48  return params;
49 }
50 
52  : MultiParameterPlasticityStressUpdate(parameters, 3, 12, 2),
53  _tensile_strength(getUserObject<SolidMechanicsHardeningModel>("tensile_strength")),
54  _compressive_strength(getUserObject<SolidMechanicsHardeningModel>("compressive_strength")),
55  _cohesion(getUserObject<SolidMechanicsHardeningModel>("cohesion")),
56  _phi(getUserObject<SolidMechanicsHardeningModel>("friction_angle")),
57  _psi(getUserObject<SolidMechanicsHardeningModel>("dilation_angle")),
58  _perfect_guess(getParam<bool>("perfect_guess")),
59  _poissons_ratio(0.0),
60  _shifter(_f_tol),
61  _eigvecs(RankTwoTensor())
62 {
63  if (_psi.value(0.0) <= 0.0 || _psi.value(0.0) > _phi.value(0.0))
64  mooseWarning("Usually the Mohr-Coulomb dilation angle is positive and not greater than the "
65  "friction angle");
66 }
67 
68 void
70  std::vector<Real> & stress_params) const
71 {
72  // stress_params[0] = smallest eigenvalue, stress_params[2] = largest eigenvalue
73  stress.symmetricEigenvalues(stress_params);
74 }
75 
76 std::vector<RankTwoTensor>
78 {
79  std::vector<Real> sp;
80  std::vector<RankTwoTensor> dsp;
81  stress.dsymmetricEigenvalues(sp, dsp);
82  return dsp;
83 }
84 
85 std::vector<RankFourTensor>
87 {
88  std::vector<RankFourTensor> d2;
89  stress.d2symmetricEigenvalues(d2);
90  return d2;
91 }
92 
93 void
94 CappedMohrCoulombStressUpdate::preReturnMapV(const std::vector<Real> & /*trial_stress_params*/,
95  const RankTwoTensor & stress_trial,
96  const std::vector<Real> & /*intnl_old*/,
97  const std::vector<Real> & /*yf*/,
98  const RankFourTensor & Eijkl)
99 {
100  std::vector<Real> eigvals;
101  stress_trial.symmetricEigenvaluesEigenvectors(eigvals, _eigvecs);
103 }
104 
105 void
107  const std::vector<Real> & stress_params,
108  Real /*gaE*/,
109  const std::vector<Real> & /*intnl*/,
110  const yieldAndFlow & /*smoothed_q*/,
111  const RankFourTensor & /*Eijkl*/,
112  RankTwoTensor & stress) const
113 {
114  // form the diagonal stress
115  stress = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0.0, 0.0, 0.0);
116  // rotate to the original frame
117  stress = _eigvecs * stress * (_eigvecs.transpose());
118 }
119 
120 void
121 CappedMohrCoulombStressUpdate::yieldFunctionValuesV(const std::vector<Real> & stress_params,
122  const std::vector<Real> & intnl,
123  std::vector<Real> & yf) const
124 {
125  // intnl[0] = shear, intnl[1] = tensile
126  const Real ts = _tensile_strength.value(intnl[1]);
127  const Real cs = _compressive_strength.value(intnl[1]);
128  const Real sinphi = std::sin(_phi.value(intnl[0]));
129  const Real cohcos = _cohesion.value(intnl[0]) * std::cos(_phi.value(intnl[0]));
130  const Real smax = stress_params[2]; // largest eigenvalue
131  const Real smid = stress_params[1];
132  const Real smin = stress_params[0]; // smallest eigenvalue
133  yf[0] = smax - ts;
134  yf[1] = smid - ts;
135  yf[2] = smin - ts;
136  yf[3] = -smin - cs;
137  yf[4] = -smid - cs;
138  yf[5] = -smax - cs;
139  yf[6] = 0.5 * (smax - smin) + 0.5 * (smax + smin) * sinphi - cohcos;
140  yf[7] = 0.5 * (smid - smin) + 0.5 * (smid + smin) * sinphi - cohcos;
141  yf[8] = 0.5 * (smax - smid) + 0.5 * (smax + smid) * sinphi - cohcos;
142  yf[9] = 0.5 * (smid - smax) + 0.5 * (smid + smax) * sinphi - cohcos;
143  yf[10] = 0.5 * (smin - smid) + 0.5 * (smin + smid) * sinphi - cohcos;
144  yf[11] = 0.5 * (smin - smax) + 0.5 * (smin + smax) * sinphi - cohcos;
145 }
146 
147 void
148 CappedMohrCoulombStressUpdate::computeAllQV(const std::vector<Real> & stress_params,
149  const std::vector<Real> & intnl,
150  std::vector<yieldAndFlow> & all_q) const
151 {
152  const Real ts = _tensile_strength.value(intnl[1]);
153  const Real dts = _tensile_strength.derivative(intnl[1]);
154  const Real cs = _compressive_strength.value(intnl[1]);
155  const Real dcs = _compressive_strength.derivative(intnl[1]);
156  const Real sinphi = std::sin(_phi.value(intnl[0]));
157  const Real dsinphi = std::cos(_phi.value(intnl[0])) * _phi.derivative(intnl[0]);
158  const Real sinpsi = std::sin(_psi.value(intnl[0]));
159  const Real dsinpsi = std::cos(_psi.value(intnl[0])) * _psi.derivative(intnl[0]);
160  const Real cohcos = _cohesion.value(intnl[0]) * std::cos(_phi.value(intnl[0]));
161  const Real dcohcos =
162  _cohesion.derivative(intnl[0]) * std::cos(_phi.value(intnl[0])) -
163  _cohesion.value(intnl[0]) * std::sin(_phi.value(intnl[0])) * _phi.derivative(intnl[0]);
164  const Real smax = stress_params[2]; // largest eigenvalue
165  const Real smid = stress_params[1];
166  const Real smin = stress_params[0]; // smallest eigenvalue
167 
168  // yield functions. See comment in yieldFunctionValuesV
169  all_q[0].f = smax - ts;
170  all_q[1].f = smid - ts;
171  all_q[2].f = smin - ts;
172  all_q[3].f = -smin - cs;
173  all_q[4].f = -smid - cs;
174  all_q[5].f = -smax - cs;
175  all_q[6].f = 0.5 * (smax - smin) + 0.5 * (smax + smin) * sinphi - cohcos;
176  all_q[7].f = 0.5 * (smid - smin) + 0.5 * (smid + smin) * sinphi - cohcos;
177  all_q[8].f = 0.5 * (smax - smid) + 0.5 * (smax + smid) * sinphi - cohcos;
178  all_q[9].f = 0.5 * (smid - smax) + 0.5 * (smid + smax) * sinphi - cohcos;
179  all_q[10].f = 0.5 * (smin - smid) + 0.5 * (smin + smid) * sinphi - cohcos;
180  all_q[11].f = 0.5 * (smin - smax) + 0.5 * (smin + smax) * sinphi - cohcos;
181 
182  // d(yield function)/d(stress_params)
183  for (unsigned yf = 0; yf < _num_yf; ++yf)
184  for (unsigned a = 0; a < _num_sp; ++a)
185  all_q[yf].df[a] = 0.0;
186  all_q[0].df[2] = 1.0;
187  all_q[1].df[1] = 1.0;
188  all_q[2].df[0] = 1.0;
189  all_q[3].df[0] = -1.0;
190  all_q[4].df[1] = -1.0;
191  all_q[5].df[2] = -1.0;
192  all_q[6].df[2] = 0.5 * (1.0 + sinphi);
193  all_q[6].df[0] = 0.5 * (-1.0 + sinphi);
194  all_q[7].df[1] = 0.5 * (1.0 + sinphi);
195  all_q[7].df[0] = 0.5 * (-1.0 + sinphi);
196  all_q[8].df[2] = 0.5 * (1.0 + sinphi);
197  all_q[8].df[1] = 0.5 * (-1.0 + sinphi);
198  all_q[9].df[1] = 0.5 * (1.0 + sinphi);
199  all_q[9].df[2] = 0.5 * (-1.0 + sinphi);
200  all_q[10].df[0] = 0.5 * (1.0 + sinphi);
201  all_q[10].df[1] = 0.5 * (-1.0 + sinphi);
202  all_q[11].df[0] = 0.5 * (1.0 + sinphi);
203  all_q[11].df[2] = 0.5 * (-1.0 + sinphi);
204 
205  // d(yield function)/d(intnl)
206  for (unsigned i = 0; i < 6; ++i)
207  all_q[i].df_di[0] = 0.0;
208  all_q[0].df_di[1] = all_q[1].df_di[1] = all_q[2].df_di[1] = -dts;
209  all_q[3].df_di[1] = all_q[4].df_di[1] = all_q[5].df_di[1] = -dcs;
210  for (unsigned i = 6; i < 12; ++i)
211  all_q[i].df_di[1] = 0.0;
212  all_q[6].df_di[0] = 0.5 * (smax + smin) * dsinphi - dcohcos;
213  all_q[7].df_di[0] = 0.5 * (smid + smin) * dsinphi - dcohcos;
214  all_q[8].df_di[0] = 0.5 * (smax + smid) * dsinphi - dcohcos;
215  all_q[9].df_di[0] = 0.5 * (smid + smax) * dsinphi - dcohcos;
216  all_q[10].df_di[0] = 0.5 * (smin + smid) * dsinphi - dcohcos;
217  all_q[11].df_di[0] = 0.5 * (smin + smax) * dsinphi - dcohcos;
218 
219  // the flow potential is just the yield function with phi->psi
220  // d(flow potential)/d(stress_params)
221  for (unsigned yf = 0; yf < 6; ++yf)
222  for (unsigned a = 0; a < _num_sp; ++a)
223  all_q[yf].dg[a] = all_q[yf].df[a];
224  all_q[6].dg[2] = all_q[7].dg[1] = all_q[8].dg[2] = all_q[9].dg[1] = all_q[10].dg[0] =
225  all_q[11].dg[0] = 0.5 * (1.0 + sinpsi);
226  all_q[6].dg[0] = all_q[7].dg[0] = all_q[8].dg[1] = all_q[9].dg[2] = all_q[10].dg[1] =
227  all_q[11].dg[2] = 0.5 * (-1.0 + sinpsi);
228 
229  // d(flow potential)/d(stress_params)/d(intnl)
230  for (unsigned yf = 0; yf < _num_yf; ++yf)
231  for (unsigned a = 0; a < _num_sp; ++a)
232  for (unsigned i = 0; i < _num_intnl; ++i)
233  all_q[yf].d2g_di[a][i] = 0.0;
234  all_q[6].d2g_di[2][0] = all_q[7].d2g_di[1][0] = all_q[8].d2g_di[2][0] = all_q[9].d2g_di[1][0] =
235  all_q[10].d2g_di[0][0] = all_q[11].d2g_di[0][0] = 0.5 * dsinpsi;
236  all_q[6].d2g_di[0][0] = all_q[7].d2g_di[0][0] = all_q[8].d2g_di[1][0] = all_q[9].d2g_di[2][0] =
237  all_q[10].d2g_di[1][0] = all_q[11].d2g_di[2][0] = 0.5 * dsinpsi;
238 
239  // d(flow potential)/d(stress_params)/d(stress_params)
240  for (unsigned yf = 0; yf < _num_yf; ++yf)
241  for (unsigned a = 0; a < _num_sp; ++a)
242  for (unsigned b = 0; b < _num_sp; ++b)
243  all_q[yf].d2g[a][b] = 0.0;
244 }
245 
246 void
248 {
249  // Eijkl is required to be isotropic, so we can use the
250  // frame where stress is diagonal
251  for (unsigned a = 0; a < _num_sp; ++a)
252  for (unsigned b = 0; b < _num_sp; ++b)
253  _Eij[a][b] = Eijkl(a, a, b, b);
254  _En = _Eij[2][2];
255  const Real denom = _Eij[0][0] * (_Eij[0][0] + _Eij[0][1]) - 2 * Utility::pow<2>(_Eij[0][1]);
256  for (unsigned a = 0; a < _num_sp; ++a)
257  {
258  _Cij[a][a] = (_Eij[0][0] + _Eij[0][1]) / denom;
259  for (unsigned b = 0; b < a; ++b)
260  _Cij[a][b] = _Cij[b][a] = -_Eij[0][1] / denom;
261  }
262 }
263 
264 void
265 CappedMohrCoulombStressUpdate::initializeVarsV(const std::vector<Real> & trial_stress_params,
266  const std::vector<Real> & intnl_old,
267  std::vector<Real> & stress_params,
268  Real & gaE,
269  std::vector<Real> & intnl) const
270 {
271  if (!_perfect_guess)
272  {
273  for (unsigned i = 0; i < _num_sp; ++i)
274  stress_params[i] = trial_stress_params[i];
275  gaE = 0.0;
276  }
277  else
278  {
279  const Real ts = _tensile_strength.value(intnl_old[1]);
280  const Real cs = _compressive_strength.value(intnl_old[1]);
281  const Real sinphi = std::sin(_phi.value(intnl_old[0]));
282  const Real cohcos = _cohesion.value(intnl_old[0]) * std::cos(_phi.value(intnl_old[0]));
283 
284  const Real trial_tensile_yf = trial_stress_params[2] - ts;
285  const Real trial_compressive_yf = -trial_stress_params[0] - cs;
286  const Real trial_mc_yf = 0.5 * (trial_stress_params[2] - trial_stress_params[0]) +
287  0.5 * (trial_stress_params[2] + trial_stress_params[0]) * sinphi -
288  cohcos;
289 
315  bool found_solution = false;
316 
317  if (trial_tensile_yf <= _f_tol && trial_compressive_yf <= _f_tol && trial_mc_yf <= _f_tol)
318  {
319  // this is needed because although we know the smoothed yield function is
320  // positive, the individual yield functions may not be
321  for (unsigned i = 0; i < _num_sp; ++i)
322  stress_params[i] = trial_stress_params[i];
323  gaE = 0.0;
324  found_solution = true;
325  }
326 
327  const bool tensile_possible = (ts < cohcos / sinphi); // tensile chops MC tip
328  const bool mc_tip_possible = (cohcos / sinphi < ts); // MC tip pokes through tensile
329  const bool mc_impossible = (0.5 * (ts + cs) + 0.5 * (ts - cs) * sinphi - cohcos <
330  _f_tol); // MC outside tensile and compressive
331 
332  const Real sinpsi = std::sin(_psi.value(intnl_old[0]));
333  const Real halfplus = 0.5 + 0.5 * sinpsi;
334  const Real neghalfplus = -0.5 + 0.5 * sinpsi;
335 
336  if (!found_solution && tensile_possible && trial_tensile_yf > _f_tol &&
337  (trial_compressive_yf <= _f_tol || (trial_compressive_yf > _f_tol && mc_impossible)))
338  {
339  // try pure tensile failure, return to the plane
340  // This involves solving yf[0] = 0 and the three flow-direction equations
341  // Don't try this if there is compressive failure, since returning to
342  // the tensile yield surface will only make compressive failure worse
343  const Real ga = (trial_stress_params[2] - ts) / _Eij[2][2];
344  stress_params[2] = ts; // largest eigenvalue
345  stress_params[1] = trial_stress_params[1] - ga * _Eij[1][2];
346  stress_params[0] = trial_stress_params[0] - ga * _Eij[0][2];
347 
348  // if we have to return to the edge, or tip, do that
349  Real dist_mod = 1.0;
350  const Real to_subtract1 = stress_params[1] - (ts - 0.5 * _shifter);
351  if (to_subtract1 > 0.0)
352  {
353  dist_mod += Utility::pow<2>(to_subtract1 / (trial_stress_params[2] - ts));
354  stress_params[1] -= to_subtract1;
355  }
356  const Real to_subtract0 = stress_params[0] - (ts - _shifter);
357  if (to_subtract0 > 0.0)
358  {
359  dist_mod += Utility::pow<2>(to_subtract0 / (trial_stress_params[2] - ts));
360  stress_params[0] -= to_subtract0;
361  }
362  if (mc_impossible) // might have to shift up to the compressive yield surface
363  {
364  const Real to_add0 = -stress_params[0] - cs;
365  if (to_add0 > 0.0)
366  {
367  dist_mod += Utility::pow<2>(to_add0 / (trial_stress_params[2] - ts));
368  stress_params[0] += to_add0;
369  }
370  const Real to_add1 = -cs + 0.5 * _shifter - stress_params[1];
371  if (to_add1 > 0.0)
372  {
373  dist_mod += Utility::pow<2>(to_add1 / (trial_stress_params[2] - ts));
374  stress_params[1] += to_add1;
375  }
376  }
377 
378  const Real new_compressive_yf = -stress_params[0] - cs;
379  const Real new_mc_yf = 0.5 * (stress_params[2] - stress_params[0]) +
380  0.5 * (stress_params[2] + stress_params[0]) * sinphi - cohcos;
381  if (new_mc_yf <= _f_tol && new_compressive_yf <= _f_tol)
382  {
383  gaE = std::sqrt(dist_mod) * (trial_stress_params[2] - stress_params[2]);
384  found_solution = true;
385  }
386  }
387  if (!found_solution && trial_compressive_yf > _f_tol &&
388  (trial_tensile_yf <= _f_tol || (trial_tensile_yf > _f_tol && mc_impossible)))
389  {
390  // try pure compressive failure
391  // Don't try this if there is tensile failure, since returning to
392  // the compressive yield surface will only make tensile failure worse
393  const Real ga = (trial_stress_params[0] + cs) / _Eij[0][0]; // this is negative
394  stress_params[0] = -cs;
395  stress_params[1] = trial_stress_params[1] - ga * _Eij[1][0];
396  stress_params[2] = trial_stress_params[2] - ga * _Eij[2][0];
397 
398  // if we have to return to the edge, or tip, do that
399  Real dist_mod = 1.0;
400  const Real to_add1 = -cs + 0.5 * _shifter - stress_params[1];
401  if (to_add1 > 0.0)
402  {
403  dist_mod += Utility::pow<2>(to_add1 / (trial_stress_params[0] + cs));
404  stress_params[1] += to_add1;
405  }
406  const Real to_add2 = -cs + _shifter - stress_params[2];
407  if (to_add2 > 0.0)
408  {
409  dist_mod += Utility::pow<2>(to_add2 / (trial_stress_params[0] + cs));
410  stress_params[2] += to_add2;
411  }
412  if (mc_impossible) // might have to shift down to the tensile yield surface
413  {
414  const Real to_subtract2 = stress_params[2] - ts;
415  if (to_subtract2 > 0.0)
416  {
417  dist_mod += Utility::pow<2>(to_subtract2 / (trial_stress_params[0] + cs));
418  stress_params[2] -= to_subtract2;
419  }
420  const Real to_subtract1 = stress_params[1] - (ts - 0.5 * _shifter);
421  if (to_subtract1 > 0.0)
422  {
423  dist_mod += Utility::pow<2>(to_subtract1 / (trial_stress_params[0] + cs));
424  stress_params[1] -= to_subtract1;
425  }
426  }
427 
428  const Real new_tensile_yf = stress_params[2] - ts;
429  const Real new_mc_yf = 0.5 * (stress_params[2] - stress_params[0]) +
430  0.5 * (stress_params[2] + stress_params[0]) * sinphi - cohcos;
431  if (new_mc_yf <= _f_tol && new_tensile_yf <= _f_tol)
432  {
433  gaE = std::sqrt(dist_mod) * (-trial_stress_params[0] + stress_params[0]);
434  found_solution = true;
435  }
436  }
437  if (!found_solution && !mc_impossible && trial_mc_yf > _f_tol)
438  {
439  // try pure shear failure, return to the plane
440  // This involves solving yf[6]=0 and the three flow-direction equations
441  const Real ga = ((trial_stress_params[2] - trial_stress_params[0]) +
442  (trial_stress_params[2] + trial_stress_params[0]) * sinphi - 2.0 * cohcos) /
443  (_Eij[2][2] - _Eij[0][2] + (_Eij[2][2] + _Eij[0][2]) * sinpsi * sinphi);
444  stress_params[2] =
445  trial_stress_params[2] - ga * (_Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus);
446  stress_params[1] = trial_stress_params[1] - ga * _Eij[1][0] * sinpsi;
447  stress_params[0] =
448  trial_stress_params[0] - ga * (_Eij[0][0] * neghalfplus + _Eij[0][2] * halfplus);
449  const Real f7 = 0.5 * (stress_params[1] - stress_params[0]) +
450  0.5 * (stress_params[1] + stress_params[0]) * sinphi - cohcos;
451  const Real f8 = 0.5 * (stress_params[2] - stress_params[1]) +
452  0.5 * (stress_params[2] + stress_params[1]) * sinphi - cohcos;
453  const Real new_tensile_yf = stress_params[2] - ts;
454  const Real new_compressive_yf = -stress_params[0] - cs;
455 
456  if (f7 <= _f_tol && f8 <= _f_tol && new_tensile_yf <= _f_tol && new_compressive_yf <= _f_tol)
457  {
458  gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
459  (stress_params[2] - stress_params[0])) /
460  (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2];
461  found_solution = true;
462  }
463  }
464  if (!found_solution && !mc_impossible && trial_mc_yf > _f_tol)
465  {
466  // Try return to the max=mid MC line.
467  // To return to the max=mid line, we need to solve f6 = 0 = f7 and
468  // the three flow-direction equations. In the flow-direction equations
469  // there are two plasticity multipliers, which i call ga6 and ga7,
470  // corresponding to the amounts of strain normal to the f6 and f7
471  // directions, respectively.
472  // So:
473  // Smax = Smax^trial - ga6 Emax,a dg6/dSa - ga7 Emax,a dg7/dSa
474  // = Smax^trial - ga6 cmax6 - ga7 cmax7 (with cmax6 and cmax7 evaluated below)
475  // Smid = Smid^trial - ga6 Emid,a dg6/dSa - ga7 Emid,a dg7/dSa
476  // = Smid^trial - ga6 cmid6 - ga7 cmid7
477  // Smin = Smin^trial - ga6 Emin,a dg6/dSa - ga7 Emin,a dg7/dSa
478  // = Smin^trial - ga6 cmin6 - ga7 cmin7
479  const Real cmax6 = _Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus;
480  const Real cmax7 = _Eij[2][1] * halfplus + _Eij[2][0] * neghalfplus;
481  // const Real cmid6 = _Eij[1][2] * halfplus + _Eij[1][0] * neghalfplus;
482  // const Real cmid7 = _Eij[1][1] * halfplus + _Eij[1][0] * neghalfplus;
483  const Real cmin6 = _Eij[0][2] * halfplus + _Eij[0][0] * neghalfplus;
484  const Real cmin7 = _Eij[0][1] * halfplus + _Eij[0][0] * neghalfplus;
485  // Substituting these into f6 = 0 yields
486  // 0 = f6_trial - ga6 (0.5(cmax6 - cmin6) + 0.5(cmax6 + cmin6)sinphi) - ga7 (0.5(cmax6 -
487  // cmin6) + 0.5(cmax6 + cmin6)sinphi) = f6_trial - ga6 c6 - ga7 c7, where
488  const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
489  const Real c7 = 0.5 * (cmax7 - cmin7) + 0.5 * (cmax7 + cmin7) * sinphi;
490  // It isn't too hard to check that the other equation is
491  // 0 = f7_trial - ga6 c7 - ga7 c6
492  // These equations may be inverted to yield
493  if (c6 != c7)
494  {
495  const Real f6_trial = trial_mc_yf;
496  const Real f7_trial = 0.5 * (trial_stress_params[1] - trial_stress_params[0]) +
497  0.5 * (trial_stress_params[1] + trial_stress_params[0]) * sinphi -
498  cohcos;
499  const Real descr = Utility::pow<2>(c6) - Utility::pow<2>(c7);
500  Real ga6 = (c6 * f6_trial - c7 * f7_trial) / descr;
501  Real ga7 = (-c7 * f6_trial + c6 * f7_trial) / descr;
502  // and finally
503  stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga7 * cmax7;
504  stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga7 * cmin7;
505  stress_params[1] = stress_params[2] - 0.5 * _shifter;
506 
507  Real f8 = 0.5 * (stress_params[2] - stress_params[1]) +
508  0.5 * (stress_params[2] + stress_params[1]) * sinphi - cohcos;
509 
510  if (mc_tip_possible && f8 > _f_tol)
511  {
512  stress_params[2] = cohcos / sinphi;
513  stress_params[1] = stress_params[2] - 0.5 * _shifter;
514  stress_params[0] = stress_params[2] - _shifter;
515  f8 = 0.0;
516  ga6 = 1.0;
517  ga7 = 1.0;
518  }
519 
520  const Real new_tensile_yf = stress_params[2] - ts;
521  const Real new_compressive_yf = -stress_params[0] - cs;
522 
523  if (f8 <= _f_tol && new_tensile_yf <= _f_tol && new_compressive_yf <= _f_tol &&
524  ga6 >= 0.0 && ga7 >= 0.0)
525  {
526  gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
527  (stress_params[2] - stress_params[0])) /
528  (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2];
529  found_solution = true;
530  }
531  }
532  }
533  if (!found_solution && !mc_impossible && trial_mc_yf > _f_tol)
534  {
535  // Try return to the mid=min line.
536  // To return to the mid=min line, we need to solve f6 = 0 = f8 and
537  // the three flow-direction equations. In the flow-direction equations
538  // there are two plasticity multipliers, which i call ga6 and ga8,
539  // corresponding to the amounts of strain normal to the f6 and f8
540  // directions, respectively.
541  // So:
542  // Smax = Smax^trial - ga6 Emax,a dg6/dSa - ga8 Emax,a dg8/dSa
543  // = Smax^trial - ga6 cmax6 - ga8 cmax8 (with cmax6 and cmax8 evaluated below)
544  // Smid = Smid^trial - ga6 Emid,a dg6/dSa - ga8 Emid,a dg8/dSa
545  // = Smid^trial - ga6 cmid6 - ga8 cmid8
546  // Smin = Smin^trial - ga6 Emin,a dg6/dSa - ga8 Emin,a dg8/dSa
547  // = Smin^trial - ga6 cmin6 - ga8 cmin8
548  const Real cmax6 = _Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus;
549  const Real cmax8 = _Eij[2][2] * halfplus + _Eij[2][1] * neghalfplus;
550  // const Real cmid6 = _Eij[1][2] * halfplus + _Eij[1][0] * neghalfplus;
551  // const Real cmid8 = _Eij[1][2] * halfplus + _Eij[1][1] * neghalfplus;
552  const Real cmin6 = _Eij[0][2] * halfplus + _Eij[0][0] * neghalfplus;
553  const Real cmin8 = _Eij[0][2] * halfplus + _Eij[0][1] * neghalfplus;
554  // Substituting these into f6 = 0 yields
555  // 0 = f6_trial - ga6 (0.5(cmax6 - cmin6) + 0.5(cmax6 + cmin6)sinphi) - ga8 (0.5(cmax6 -
556  // cmin6) + 0.5(cmax6 + cmin6)sinphi) = f6_trial - ga6 c6 - ga8 c8, where
557  const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
558  const Real c8 = 0.5 * (cmax8 - cmin8) + 0.5 * (cmax8 + cmin8) * sinphi;
559  // It isn't too hard to check that the other equation is
560  // 0 = f8_trial - ga6 c8 - ga8 c6
561  // These equations may be inverted to yield
562  if (c6 != c8)
563  {
564  const Real f6_trial = trial_mc_yf;
565  const Real f8_trial = 0.5 * (trial_stress_params[2] - trial_stress_params[1]) +
566  0.5 * (trial_stress_params[2] + trial_stress_params[1]) * sinphi -
567  cohcos;
568  const Real descr = Utility::pow<2>(c6) - Utility::pow<2>(c8);
569  Real ga6 = (c6 * f6_trial - c8 * f8_trial) / descr;
570  Real ga8 = (-c8 * f6_trial + c6 * f8_trial) / descr;
571  // and finally
572  stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga8 * cmax8;
573  stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga8 * cmin8;
574  stress_params[1] = stress_params[0] + 0.5 * _shifter;
575 
576  Real f7 = 0.5 * (stress_params[1] - stress_params[0]) +
577  0.5 * (stress_params[1] + stress_params[0]) * sinphi - cohcos;
578 
579  if (mc_tip_possible && f7 > _f_tol)
580  {
581  stress_params[2] = cohcos / sinphi;
582  stress_params[1] = stress_params[2] - 0.5 * _shifter;
583  stress_params[0] = stress_params[2] - _shifter;
584  f7 = 0.0;
585  ga6 = 1.0;
586  ga8 = 1.0;
587  }
588 
589  const Real new_tensile_yf = stress_params[2] - ts;
590  const Real new_compressive_yf = -stress_params[0] - cs;
591 
592  if (f7 <= _f_tol && new_tensile_yf <= _f_tol && new_compressive_yf <= _f_tol &&
593  ga6 >= 0.0 && ga8 >= 0.0)
594  {
595  gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
596  (stress_params[2] - stress_params[0])) /
597  (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2];
598  found_solution = true;
599  }
600  }
601  }
602  if (!found_solution && !mc_impossible && tensile_possible && trial_tensile_yf > _f_tol)
603  {
604  // Return to the line where yf[0] = 0 = yf[6].
605  // To return to this line, we need to solve f0 = 0 = f6 and
606  // the three flow-direction equations. In the flow-direction equations
607  // there are two plasticity multipliers, which i call ga0 and ga6
608  // corresponding to the amounts of strain normal to the f0 and f6
609  // directions, respectively.
610  // So:
611  // Smax = Smax^trial - ga6 Emax,a dg6/dSa - ga0 Emax,a dg0/dSa
612  // = Smax^trial - ga6 cmax6 - ga0 cmax0 (with cmax6 and cmax0 evaluated below)
613  // Smid = Smid^trial - ga6 Emid,a dg6/dSa - ga0 Emid,a dg0/dSa
614  // = Smid^trial - ga6 cmid6 - ga0 cmid0
615  // Smin = Smin^trial - ga6 Emin,a dg6/dSa - ga0 Emin,a dg0/dSa
616  // = Smin^trial - ga6 cmin6 - ga0 cmin0
617  const Real cmax6 = _Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus;
618  const Real cmax0 = _Eij[2][2];
619  const Real cmid6 = _Eij[1][2] * halfplus + _Eij[1][0] * neghalfplus;
620  const Real cmid0 = _Eij[1][2];
621  const Real cmin6 = _Eij[0][2] * halfplus + _Eij[0][0] * neghalfplus;
622  const Real cmin0 = _Eij[0][2];
623  // Substituting these into f6 = 0 yields
624  // 0 = f6_trial - ga6 (0.5(cmax6 - cmin6) + 0.5(cmax6 + cmin6)sinphi) - ga0 (0.5(cmax0 -
625  // cmin0) + 0.5(cmax0 + cmin0)sinphi) = f6_trial - ga6 c6 - ga0 c0, where
626  const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
627  const Real c0 = 0.5 * (cmax0 - cmin0) + 0.5 * (cmax0 + cmin0) * sinphi;
628  // Substituting these into f0 = 0 yields
629  // 0 = f0_trial - ga6 cmax6 - ga0 cmax0
630  // These equations may be inverted to yield the following
631  const Real descr = c0 * cmax6 - c6 * cmax0;
632  if (descr != 0.0)
633  {
634  const Real ga0 = (-c6 * trial_tensile_yf + cmax6 * trial_mc_yf) / descr;
635  const Real ga6 = (c0 * trial_tensile_yf - cmax0 * trial_mc_yf) / descr;
636  stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga0 * cmax0;
637  stress_params[1] = trial_stress_params[1] - ga6 * cmid6 - ga0 * cmid0;
638  stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga0 * cmin0;
639 
640  // enforce physicality (go to corners if necessary)
641  stress_params[0] =
642  std::min(stress_params[0],
643  stress_params[2] - _shifter); // account for poor choice from user
644  // if goto_corner then the max(min()) in the subsequent line will force the solution to lie
645  // at the corner where max = mid = tensile. This means the signs of ga0 and ga6 become
646  // irrelevant in the check below
647  const bool goto_corner = (stress_params[1] >= stress_params[2] - 0.5 * _shifter);
648  stress_params[1] = std::max(std::min(stress_params[1], stress_params[2] - 0.5 * _shifter),
649  stress_params[0] + 0.5 * _shifter);
650 
651  const Real new_compressive_yf = -stress_params[0] - cs;
652  if (new_compressive_yf <= _f_tol &&
653  (goto_corner || (ga0 >= 0.0 && ga6 >= 0.0))) // enforce ga>=0 unless going to a corner
654  {
655  gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
656  (stress_params[2] - stress_params[0])) /
657  (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2] +
658  (trial_stress_params[2] - stress_params[2]);
659  found_solution = true;
660  }
661  }
662  }
663  if (!found_solution && !mc_impossible)
664  {
665  // Return to the line where yf[3] = 0 = yf[6].
666  // To return to this line, we need to solve f3 = 0 = f6 and
667  // the three flow-direction equations. In the flow-direction equations
668  // there are two plasticity multipliers, which i call ga3 and ga6
669  // corresponding to the amounts of strain normal to the f3 and f6
670  // directions, respectively.
671  // So:
672  // Smax = Smax^trial - ga6 Emax,a dg6/dSa - ga3 Emax,a dg3/dSa
673  // = Smax^trial - ga6 cmax6 - ga3 cmax3 (with cmax6 and cmax3 evaluated below)
674  // Smid = Smid^trial - ga6 Emid,a dg6/dSa - ga3 Emid,a dg3/dSa
675  // = Smid^trial - ga6 cmid6 - ga3 cmid3
676  // Smin = Smin^trial - ga6 Emin,a dg6/dSa - ga3 Emin,a dg3/dSa
677  // = Smin^trial - ga6 cmin6 - ga3 cmin3
678  const Real cmax6 = _Eij[2][2] * halfplus + _Eij[2][0] * neghalfplus;
679  const Real cmax3 = -_Eij[2][0];
680  const Real cmid6 = _Eij[1][2] * halfplus + _Eij[1][0] * neghalfplus;
681  const Real cmid3 = -_Eij[1][0];
682  const Real cmin6 = _Eij[0][2] * halfplus + _Eij[0][0] * neghalfplus;
683  const Real cmin3 = -_Eij[0][0];
684  // Substituting these into f6 = 0 yields
685  // 0 = f6_trial - ga6 (0.5(cmax6 - cmin6) + 0.5(cmax6 + cmin6)sinphi) - ga3 (0.5(cmax3 -
686  // cmin3) + 0.5(cmax3 + cmin3)sinphi) = f6_trial - ga6 c6 - ga3 c3, where
687  const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
688  const Real c3 = 0.5 * (cmax3 - cmin3) + 0.5 * (cmax3 + cmin3) * sinphi;
689  // Substituting these into f3 = 0 yields
690  // 0 = - f3_trial - ga6 cmin6 - ga3 cmin3
691  // These equations may be inverted to yield the following
692  const Real descr = c3 * cmin6 - c6 * cmin3;
693  if (descr != 0.0)
694  {
695  const Real ga3 = (c6 * trial_compressive_yf + cmin6 * trial_mc_yf) / descr;
696  const Real ga6 = (-c3 * trial_compressive_yf - cmin3 * trial_mc_yf) / descr;
697  stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga3 * cmax3;
698  stress_params[1] = trial_stress_params[1] - ga6 * cmid6 - ga3 * cmid3;
699  stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga3 * cmin3;
700 
701  const Real new_tensile_yf = stress_params[2] - ts;
702  stress_params[2] =
703  std::max(stress_params[2],
704  stress_params[0] + _shifter); // account for poor choice from user
705  stress_params[1] = std::max(std::min(stress_params[1], stress_params[2] - 0.5 * _shifter),
706  stress_params[0] + 0.5 * _shifter);
707 
708  if (new_tensile_yf <= _f_tol && ga6 >= 0.0)
709  {
710  gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
711  (stress_params[2] - stress_params[0])) /
712  (_Eij[2][2] - _Eij[0][2]) * _Eij[2][2] +
713  (-trial_stress_params[0] - stress_params[0]);
714 
715  found_solution = true;
716  }
717  }
718  }
719  if (!found_solution)
720  {
721  // Cannot find an acceptable initialisation
722  for (unsigned i = 0; i < _num_sp; ++i)
723  stress_params[i] = trial_stress_params[i];
724  gaE = 0.0;
725  mooseWarning("CappedMohrCoulombStressUpdate cannot initialize from max = ",
726  stress_params[2],
727  " mid = ",
728  stress_params[1],
729  " min = ",
730  stress_params[0]);
731  }
732  }
733  setIntnlValuesV(trial_stress_params, stress_params, intnl_old, intnl);
734 }
735 
736 void
737 CappedMohrCoulombStressUpdate::setIntnlValuesV(const std::vector<Real> & trial_stress_params,
738  const std::vector<Real> & current_stress_params,
739  const std::vector<Real> & intnl_old,
740  std::vector<Real> & intnl) const
741 {
742  // intnl[0] = shear, intnl[1] = tensile
743  const Real smax = current_stress_params[2]; // largest eigenvalue
744  const Real smin = current_stress_params[0]; // smallest eigenvalue
745  const Real trial_smax = trial_stress_params[2]; // largest eigenvalue
746  const Real trial_smin = trial_stress_params[0]; // smallest eigenvalue
747  const Real ga_shear = ((trial_smax - trial_smin) - (smax - smin)) / (_Eij[2][2] - _Eij[0][2]);
748  intnl[0] = intnl_old[0] + ga_shear;
749  const Real sinpsi = std::sin(_psi.value(intnl[0]));
750  const Real prefactor = (_Eij[2][2] + _Eij[0][2]) * sinpsi;
751  const Real shear_correction = prefactor * ga_shear;
752  const Real ga_tensile = (1.0 - _poissons_ratio) *
753  ((trial_smax + trial_smin) - (smax + smin) - shear_correction) /
754  _Eij[2][2];
755  intnl[1] = intnl_old[1] + ga_tensile;
756 }
757 
758 void
759 CappedMohrCoulombStressUpdate::setIntnlDerivativesV(const std::vector<Real> & trial_stress_params,
760  const std::vector<Real> & current_stress_params,
761  const std::vector<Real> & intnl,
762  std::vector<std::vector<Real>> & dintnl) const
763 {
764  // intnl[0] = shear, intnl[1] = tensile
765  const Real smax = current_stress_params[2]; // largest eigenvalue
766  const Real smin = current_stress_params[0]; // smallest eigenvalue
767  const Real trial_smax = trial_stress_params[2]; // largest eigenvalue
768  const Real trial_smin = trial_stress_params[0]; // smallest eigenvalue
769  const Real ga_shear = ((trial_smax - trial_smin) - (smax - smin)) / (_Eij[2][2] - _Eij[0][2]);
770  const std::vector<Real> dga_shear = {
771  1.0 / (_Eij[2][2] - _Eij[0][2]), 0.0, -1.0 / (_Eij[2][2] - _Eij[0][2])};
772  // intnl[0] = intnl_old[0] + ga_shear;
773  for (std::size_t i = 0; i < _num_sp; ++i)
774  dintnl[0][i] = dga_shear[i];
775 
776  const Real sinpsi = std::sin(_psi.value(intnl[0]));
777  const Real dsinpsi_di0 = _psi.derivative(intnl[0]) * std::cos(_psi.value(intnl[0]));
778 
779  const Real prefactor = (_Eij[2][2] + _Eij[0][2]) * sinpsi;
780  const Real dprefactor_di0 = (_Eij[2][2] + _Eij[0][2]) * dsinpsi_di0;
781  // const Real shear_correction = prefactor * ga_shear;
782  std::vector<Real> dshear_correction(_num_sp);
783  for (std::size_t i = 0; i < _num_sp; ++i)
784  dshear_correction[i] = prefactor * dga_shear[i] + dprefactor_di0 * dintnl[0][i] * ga_shear;
785  // const Real ga_tensile = (1 - _poissons_ratio) * ((trial_smax + trial_smin) - (smax + smin) -
786  // shear_correction) /
787  // _Eij[2][2];
788  // intnl[1] = intnl_old[1] + ga_tensile;
789  for (std::size_t i = 0; i < _num_sp; ++i)
790  dintnl[1][i] = -(1.0 - _poissons_ratio) * dshear_correction[i] / _Eij[2][2];
791  dintnl[1][2] += -(1.0 - _poissons_ratio) / _Eij[2][2];
792  dintnl[1][0] += -(1.0 - _poissons_ratio) / _Eij[2][2];
793 }
794 
795 void
797  const RankTwoTensor & stress_trial,
798  const std::vector<Real> & trial_stress_params,
799  const RankTwoTensor & /*stress*/,
800  const std::vector<Real> & stress_params,
801  Real /*gaE*/,
802  const yieldAndFlow & /*smoothed_q*/,
803  const RankFourTensor & elasticity_tensor,
804  bool compute_full_tangent_operator,
805  const std::vector<std::vector<Real>> & dvar_dtrial,
806  RankFourTensor & cto)
807 {
808  cto = elasticity_tensor;
809  if (!compute_full_tangent_operator)
810  return;
811 
812  // dvar_dtrial has been computed already, so
813  // d(stress)/d(trial_stress) = d(eigvecs * stress_params * eigvecs.transpose())/d(trial_stress)
814  // eigvecs is a rotation matrix, rot(i, j) = e_j(i) = i^th component of j^th eigenvector
815  // d(rot_ij)/d(stress_kl) = d(e_j(i))/d(stress_kl)
816  // = sum_a 0.5 * e_a(i) * (e_a(k)e_j(l) + e_a(l)e_j(k)) / (la_j - la_a)
817  // = sum_a 0.5 * rot(i,a) * (rot(k,a)rot(l,j) + rot(l,a)*rot(k,j)) / (la_j - la_a)
818  RankFourTensor drot_dstress;
819  for (unsigned i = 0; i < _tensor_dimensionality; ++i)
820  for (unsigned j = 0; j < _tensor_dimensionality; ++j)
821  for (unsigned k = 0; k < _tensor_dimensionality; ++k)
822  for (unsigned l = 0; l < _tensor_dimensionality; ++l)
823  for (unsigned a = 0; a < _num_sp; ++a)
824  {
825  if (trial_stress_params[a] == trial_stress_params[j])
826  continue;
827  drot_dstress(i, j, k, l) +=
828  0.5 * _eigvecs(i, a) *
829  (_eigvecs(k, a) * _eigvecs(l, j) + _eigvecs(l, a) * _eigvecs(k, j)) /
830  (trial_stress_params[j] - trial_stress_params[a]);
831  }
832 
833  const RankTwoTensor eT = _eigvecs.transpose();
834 
835  RankFourTensor dstress_dtrial;
836  for (unsigned i = 0; i < _tensor_dimensionality; ++i)
837  for (unsigned j = 0; j < _tensor_dimensionality; ++j)
838  for (unsigned k = 0; k < _tensor_dimensionality; ++k)
839  for (unsigned l = 0; l < _tensor_dimensionality; ++l)
840  for (unsigned a = 0; a < _num_sp; ++a)
841  dstress_dtrial(i, j, k, l) +=
842  drot_dstress(i, a, k, l) * stress_params[a] * eT(a, j) +
843  _eigvecs(i, a) * stress_params[a] * drot_dstress(j, a, k, l);
844 
845  const std::vector<RankTwoTensor> dsp_trial = dstress_param_dstress(stress_trial);
846  for (unsigned i = 0; i < _tensor_dimensionality; ++i)
847  for (unsigned j = 0; j < _tensor_dimensionality; ++j)
848  for (unsigned k = 0; k < _tensor_dimensionality; ++k)
849  for (unsigned l = 0; l < _tensor_dimensionality; ++l)
850  for (unsigned a = 0; a < _num_sp; ++a)
851  for (unsigned b = 0; b < _num_sp; ++b)
852  dstress_dtrial(i, j, k, l) +=
853  _eigvecs(i, a) * dvar_dtrial[a][b] * dsp_trial[b](k, l) * eT(a, j);
854 
855  cto = dstress_dtrial * elasticity_tensor;
856 }
const SolidMechanicsHardeningModel & _phi
Hardening model for friction angle.
void symmetricEigenvalues(std::vector< Real > &eigvals) const
CappedMohrCoulombStressUpdate implements rate-independent nonassociative Mohr-Coulomb plus tensile pl...
void computeStressParams(const RankTwoTensor &stress, std::vector< Real > &stress_params) const override
Computes stress_params, given stress.
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
void computeAllQV(const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< yieldAndFlow > &all_q) const override
Completely fills all_q with correct values.
void dsymmetricEigenvalues(std::vector< Real > &eigvals, std::vector< RankTwoTensorTempl< Real >> &deigvals) const
const bool _perfect_guess
Whether to provide an estimate of the returned stress, based on perfect plasticity.
const unsigned _num_intnl
Number of internal parameters.
virtual void consistentTangentOperatorV(const RankTwoTensor &stress_trial, const std::vector< Real > &trial_stress_params, const RankTwoTensor &stress, const std::vector< Real > &stress_params, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, bool compute_full_tangent_operator, const std::vector< std::vector< Real >> &dvar_dtrial, RankFourTensor &cto) override
Calculates the consistent tangent operator.
Struct designed to hold info about a single yield function and its derivatives, as well as the flow d...
const SolidMechanicsHardeningModel & _tensile_strength
Hardening model for tensile strength.
std::vector< std::vector< Real > > _Cij
_Cij[i, j] * _Eij[j, k] = 1 iff j == k
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
virtual Real value(Real intnl) const
void mooseWarning(Args &&... args) const
void addRequiredParam(const std::string &name, const std::string &doc_string)
Real elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
void setEffectiveElasticity(const RankFourTensor &Eijkl) override
Sets _Eij and _En and _Cij.
std::vector< std::vector< Real > > _Eij
E[i, j] in the system of equations to be solved.
const Real _f_tol
The yield-function tolerance.
const unsigned _num_yf
Number of yield functions.
CappedMohrCoulombStressUpdate(const InputParameters &parameters)
std::vector< RankTwoTensor > dstress_param_dstress(const RankTwoTensor &stress) const override
d(stress_param[i])/d(stress) at given stress
void setIntnlValuesV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &intnl) const override
Sets the internal parameters based on the trial values of stress_params, their current values...
const Real _shifter
When equal-eigenvalues are predicted from the stress initialization routine, shift them by this amoun...
void d2symmetricEigenvalues(std::vector< RankFourTensorTempl< Real >> &deriv) const
virtual void preReturnMapV(const std::vector< Real > &trial_stress_params, const RankTwoTensor &stress_trial, const std::vector< Real > &intnl_old, const std::vector< Real > &yf, const RankFourTensor &Eijkl) override
Derived classes may employ this function to record stuff or do other computations prior to the return...
void yieldFunctionValuesV(const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< Real > &yf) const override
Computes the values of the yield functions, given stress_params and intnl parameters.
registerMooseObject("SolidMechanicsApp", CappedMohrCoulombStressUpdate)
void setIntnlDerivativesV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &current_stress_params, const std::vector< Real > &intnl, std::vector< std::vector< Real >> &dintnl) const override
Sets the derivatives of internal parameters, based on the trial values of stress_params, their current values, and the current values of the internal parameters.
const SolidMechanicsHardeningModel & _compressive_strength
Hardening model for compressive strength.
virtual Real derivative(Real intnl) const
RankTwoTensorTempl< Real > transpose() const
const SolidMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
T getIsotropicPoissonsRatio(const RankFourTensorTempl< T > &elasticity_tensor)
Get the Poisson&#39;s modulus for an isotropic elasticity tensor param elasticity_tensor the tensor (must...
void initializeVarsV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &stress_params, Real &gaE, std::vector< Real > &intnl) const override
Sets (stress_params, intnl) at "good guesses" of the solution to the Return-Map algorithm.
void symmetricEigenvaluesEigenvectors(std::vector< Real > &eigvals, RankTwoTensorTempl< Real > &eigvecs) const
MultiParameterPlasticityStressUpdate performs the return-map algorithm and associated stress updates ...
const SolidMechanicsHardeningModel & _psi
Hardening model for dilation angle.
RankTwoTensor _eigvecs
Eigenvectors of the trial stress as a RankTwoTensor, in order to rotate the returned stress back to s...
void addClassDescription(const std::string &doc_string)
virtual void setStressAfterReturnV(const RankTwoTensor &stress_trial, const std::vector< Real > &stress_params, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, RankTwoTensor &stress) const override
Sets stress from the admissible parameters.
std::vector< RankFourTensor > d2stress_param_dstress(const RankTwoTensor &stress) const override
d2(stress_param[i])/d(stress)/d(stress) at given stress
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
static constexpr unsigned _tensor_dimensionality
Internal dimensionality of tensors (currently this is 3 throughout tensor_mechanics) ...
static const std::string k
Definition: NS.h:124
const unsigned _num_sp
Number of stress parameters.