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