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