www.mooseframework.org
CappedWeakPlaneStressUpdate.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
8 
9 #include "libmesh/utility.h"
10 
11 template <>
12 InputParameters
14 {
15  InputParameters params = validParams<TwoParameterPlasticityStressUpdate>();
16  params.addClassDescription("Capped weak-plane plasticity stress calculator");
17  params.addRequiredParam<UserObjectName>(
18  "cohesion",
19  "A TensorMechanicsHardening UserObject that defines hardening of the cohesion. "
20  "Physically the cohesion should not be negative.");
21  params.addRequiredParam<UserObjectName>("tan_friction_angle",
22  "A TensorMechanicsHardening UserObject that defines "
23  "hardening of tan(friction angle). Physically the "
24  "friction angle should be between 0 and 90deg.");
25  params.addRequiredParam<UserObjectName>(
26  "tan_dilation_angle",
27  "A TensorMechanicsHardening UserObject that defines hardening of the "
28  "tan(dilation angle). Usually the dilation angle is not greater than "
29  "the friction angle, and it is between 0 and 90deg.");
30  params.addRequiredParam<UserObjectName>(
31  "tensile_strength",
32  "A TensorMechanicsHardening UserObject that defines hardening of the "
33  "weak-plane tensile strength. In physical situations this is positive "
34  "(and always must be greater than negative compressive-strength.");
35  params.addRequiredParam<UserObjectName>("compressive_strength",
36  "A TensorMechanicsHardening UserObject that defines "
37  "hardening of the weak-plane compressive strength. In "
38  "physical situations this is positive.");
39  params.addRequiredRangeCheckedParam<Real>(
40  "tip_smoother",
41  "tip_smoother>=0",
42  "The cone vertex at shear-stress = 0 will be smoothed by "
43  "the given amount. Typical value is 0.1*cohesion");
44  params.addParam<bool>("perfect_guess",
45  true,
46  "Provide a guess to the Newton-Raphson proceedure "
47  "that is the result from perfect plasticity. With "
48  "severe hardening/softening this may be "
49  "suboptimal.");
50  return params;
51 }
52 
54  : TwoParameterPlasticityStressUpdate(parameters, 3, 2),
55  _cohesion(getUserObject<TensorMechanicsHardeningModel>("cohesion")),
56  _tan_phi(getUserObject<TensorMechanicsHardeningModel>("tan_friction_angle")),
57  _tan_psi(getUserObject<TensorMechanicsHardeningModel>("tan_dilation_angle")),
58  _tstrength(getUserObject<TensorMechanicsHardeningModel>("tensile_strength")),
59  _cstrength(getUserObject<TensorMechanicsHardeningModel>("compressive_strength")),
60  _small_smoother2(Utility::pow<2>(getParam<Real>("tip_smoother"))),
61  _perfect_guess(getParam<bool>("perfect_guess")),
62  _stress_return_type(StressReturnType::nothing_special),
63  _in_trial02(0.0),
64  _in_trial12(0.0),
65  _in_q_trial(0.0)
66 {
67  // With arbitary UserObjects, it is impossible to check everything,
68  // but this will catch the common errors
69  if (_tan_phi.value(0) < 0 || _tan_psi.value(0) < 0)
70  mooseError("CappedWeakPlaneStressUpdate: Weak-plane friction and dilation angles must lie in "
71  "[0, Pi/2]");
72  if (_tan_phi.value(0) < _tan_psi.value(0))
73  mooseError("CappedWeakPlaneStressUpdate: Weak-plane friction angle must not be less than "
74  "dilation angle");
75  if (_cohesion.value(0) < 0)
76  mooseError("CappedWeakPlaneStressUpdate: Weak-plane cohesion must not be negative");
78  mooseError("CappedWeakPlaneStressUpdate: Weak plane tensile strength plus compressive "
79  "strength must be greater than smoothing_tol");
80 }
81 
82 void
84 {
86 }
87 
88 void
89 CappedWeakPlaneStressUpdate::finalizeReturnProcess(const RankTwoTensor & /*rotation_increment*/)
90 {
92 }
93 
94 void
96  Real q_trial,
97  const RankTwoTensor & stress_trial,
98  const std::vector<Real> & /*intnl_old*/,
99  const std::vector<Real> & yf,
100  const RankFourTensor & /*Eijkl*/)
101 {
102  // If it's obvious, then simplify the return-type
103  if (yf[1] >= 0)
105  else if (yf[2] >= 0)
107 
108  // The following are useful for the Cosserat case too
109  _in_trial02 = stress_trial(0, 2);
110  _in_trial12 = stress_trial(1, 2);
111  _in_q_trial = q_trial;
112 }
113 
114 void
115 CappedWeakPlaneStressUpdate::computePQ(const RankTwoTensor & stress, Real & p, Real & q) const
116 {
117  p = stress(2, 2);
118  // Because the following is not explicitly symmeterised, it is useful for the Cosserat case too
119  q = std::sqrt(Utility::pow<2>(stress(0, 2)) + Utility::pow<2>(stress(1, 2)));
120 }
121 
122 void
123 CappedWeakPlaneStressUpdate::setEppEqq(const RankFourTensor & Eijkl, Real & Epp, Real & Eqq) const
124 {
125  Epp = Eijkl(2, 2, 2, 2);
126  Eqq = Eijkl(0, 2, 0, 2);
127 }
128 
129 void
130 CappedWeakPlaneStressUpdate::setStressAfterReturn(const RankTwoTensor & stress_trial,
131  Real p_ok,
132  Real q_ok,
133  Real gaE,
134  const std::vector<Real> & /*intnl*/,
135  const yieldAndFlow & smoothed_q,
136  const RankFourTensor & Eijkl,
137  RankTwoTensor & stress) const
138 {
139  stress = stress_trial;
140  stress(2, 2) = p_ok;
141  // stress_xx and stress_yy are sitting at their trial-stress values
142  // so need to bring them back via Poisson's ratio
143  stress(0, 0) -= Eijkl(2, 2, 0, 0) * gaE / _Epp * smoothed_q.dg[0];
144  stress(1, 1) -= Eijkl(2, 2, 1, 1) * gaE / _Epp * smoothed_q.dg[0];
145  if (_in_q_trial == 0.0)
146  stress(2, 0) = stress(2, 1) = stress(0, 2) = stress(1, 2) = 0.0;
147  else
148  {
149  stress(2, 0) = stress(0, 2) = _in_trial02 * q_ok / _in_q_trial;
150  stress(2, 1) = stress(1, 2) = _in_trial12 * q_ok / _in_q_trial;
151  }
152 }
153 
154 void
156  Real q_trial,
157  Real /*p*/,
158  Real q,
159  const std::vector<Real> & intnl,
160  std::vector<std::vector<Real>> & dintnl) const
161 {
162  const Real tanpsi = _tan_psi.value(intnl[0]);
163  dintnl[0][0] = 0.0;
164  dintnl[0][1] = -1.0 / _Eqq;
165  dintnl[1][0] = -1.0 / _Epp;
166  dintnl[1][1] =
167  tanpsi / _Eqq - (q_trial - q) * _tan_psi.derivative(intnl[0]) * dintnl[0][1] / _Eqq;
168 }
169 
170 void
171 CappedWeakPlaneStressUpdate::consistentTangentOperator(const RankTwoTensor & /*stress_trial*/,
172  Real /*p_trial*/,
173  Real q_trial,
174  const RankTwoTensor & /*stress*/,
175  Real /*p*/,
176  Real q,
177  Real gaE,
178  const yieldAndFlow & smoothed_q,
179  const RankFourTensor & Eijkl,
180  bool compute_full_tangent_operator,
181  RankFourTensor & cto) const
182 {
183  cto = Eijkl;
184  if (!compute_full_tangent_operator)
185  return;
186 
187  const Real Ezzzz = Eijkl(2, 2, 2, 2);
188  const Real Ezxzx = Eijkl(2, 0, 2, 0);
189  const Real tanpsi = _tan_psi.value(_intnl[_qp][0]);
190  const Real dintnl0_dq = -1.0 / Ezxzx;
191  const Real dintnl0_dqt = 1.0 / Ezxzx;
192  const Real dintnl1_dp = -1.0 / Ezzzz;
193  const Real dintnl1_dpt = 1.0 / Ezzzz;
194  const Real dintnl1_dq =
195  tanpsi / Ezxzx - (q_trial - q) * _tan_psi.derivative(_intnl[_qp][0]) * dintnl0_dq / Ezxzx;
196  const Real dintnl1_dqt =
197  -tanpsi / Ezxzx - (q_trial - q) * _tan_psi.derivative(_intnl[_qp][0]) * dintnl0_dqt / Ezxzx;
198 
199  for (unsigned i = 0; i < _tensor_dimensionality; ++i)
200  {
201  const Real dpt_depii = Eijkl(2, 2, i, i);
202  cto(2, 2, i, i) = _dp_dpt * dpt_depii;
203  const Real poisson_effect =
204  Eijkl(2, 2, 0, 0) / Ezzzz *
205  (_dgaE_dpt * smoothed_q.dg[0] + gaE * smoothed_q.d2g[0][0] * _dp_dpt +
206  gaE * smoothed_q.d2g[0][1] * _dq_dpt +
207  gaE * smoothed_q.d2g_di[0][0] * (dintnl0_dq * _dq_dpt) +
208  gaE * smoothed_q.d2g_di[0][1] *
209  (dintnl1_dpt + dintnl1_dp * _dp_dpt + dintnl1_dq * _dq_dpt)) *
210  dpt_depii;
211  cto(0, 0, i, i) -= poisson_effect;
212  cto(1, 1, i, i) -= poisson_effect;
213  if (q_trial > 0.0)
214  {
215  cto(2, 0, i, i) = cto(0, 2, i, i) = _in_trial02 / q_trial * _dq_dpt * dpt_depii;
216  cto(2, 1, i, i) = cto(1, 2, i, i) = _in_trial12 / q_trial * _dq_dpt * dpt_depii;
217  }
218  }
219 
220  const Real poisson_effect =
221  -Eijkl(2, 2, 0, 0) / Ezzzz *
222  (_dgaE_dqt * smoothed_q.dg[0] + gaE * smoothed_q.d2g[0][0] * _dp_dqt +
223  gaE * smoothed_q.d2g[0][1] * _dq_dqt +
224  gaE * smoothed_q.d2g_di[0][0] * (dintnl0_dqt + dintnl0_dq * _dq_dqt) +
225  gaE * smoothed_q.d2g_di[0][1] * (dintnl1_dqt + dintnl1_dp * _dp_dqt + dintnl1_dq * _dq_dqt));
226 
227  const Real dqt_dep20 = (q_trial == 0.0 ? 1.0 : _in_trial02 / q_trial) * Eijkl(2, 0, 2, 0);
228  cto(2, 2, 2, 0) = cto(2, 2, 0, 2) = _dp_dqt * dqt_dep20;
229  cto(0, 0, 2, 0) = cto(0, 0, 0, 2) = cto(1, 1, 2, 0) = cto(1, 1, 0, 2) =
230  poisson_effect * dqt_dep20;
231  if (q_trial > 0.0)
232  {
233  // for q_trial=0, Jacobian_mult is just given by the elastic case
234  cto(0, 2, 0, 2) = cto(2, 0, 0, 2) = cto(0, 2, 2, 0) = cto(2, 0, 2, 0) =
235  Eijkl(2, 0, 2, 0) * q / q_trial +
236  _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep20;
237  cto(1, 2, 0, 2) = cto(2, 1, 0, 2) = cto(1, 2, 2, 0) = cto(2, 1, 2, 0) =
238  _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep20;
239  }
240 
241  const Real dqt_dep21 = (q_trial == 0.0 ? 1.0 : _in_trial12 / q_trial) * Eijkl(2, 1, 2, 1);
242  cto(2, 2, 2, 1) = cto(2, 2, 1, 2) = _dp_dqt * dqt_dep21;
243  cto(0, 0, 2, 1) = cto(0, 0, 1, 2) = cto(1, 1, 2, 1) = cto(1, 1, 1, 2) =
244  poisson_effect * dqt_dep21;
245  if (q_trial > 0.0)
246  {
247  // for q_trial=0, Jacobian_mult is just given by the elastic case
248  cto(0, 2, 1, 2) = cto(2, 0, 1, 2) = cto(0, 2, 2, 1) = cto(2, 0, 2, 1) =
249  _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep21;
250  cto(1, 2, 1, 2) = cto(2, 1, 1, 2) = cto(1, 2, 2, 1) = cto(2, 1, 2, 1) =
251  Eijkl(2, 1, 2, 1) * q / q_trial +
252  _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep21;
253  }
254 }
255 
256 void
258  Real q,
259  const std::vector<Real> & intnl,
260  std::vector<Real> & yf) const
261 {
262  yf[0] = std::sqrt(Utility::pow<2>(q) + _small_smoother2) + p * _tan_phi.value(intnl[0]) -
263  _cohesion.value(intnl[0]);
264 
266  yf[1] = std::numeric_limits<Real>::lowest();
267  else
268  yf[1] = p - _tstrength.value(intnl[1]);
269 
271  yf[2] = std::numeric_limits<Real>::lowest();
272  else
273  yf[2] = -p - _cstrength.value(intnl[1]);
274 }
275 
276 void
278  Real q,
279  const std::vector<Real> & intnl,
280  std::vector<yieldAndFlow> & all_q) const
281 {
282  // yield function values
283  all_q[0].f = std::sqrt(Utility::pow<2>(q) + _small_smoother2) + p * _tan_phi.value(intnl[0]) -
284  _cohesion.value(intnl[0]);
286  all_q[1].f = std::numeric_limits<Real>::lowest();
287  else
288  all_q[1].f = p - _tstrength.value(intnl[1]);
290  all_q[2].f = std::numeric_limits<Real>::lowest();
291  else
292  all_q[2].f = -p - _cstrength.value(intnl[1]);
293 
294  // d(yield Function)/d(p, q)
295  // derivatives wrt p
296  all_q[0].df[0] = _tan_phi.value(intnl[0]);
297  all_q[1].df[0] = 1.0;
298  all_q[2].df[0] = -1.0;
299 
300  // derivatives wrt q
301  if (_small_smoother2 == 0.0)
302  all_q[0].df[1] = 1.0;
303  else
304  all_q[0].df[1] = q / std::sqrt(Utility::pow<2>(q) + _small_smoother2);
305  all_q[1].df[1] = 0.0;
306  all_q[2].df[1] = 0.0;
307 
308  // d(yield Function)/d(intnl)
309  // derivatives wrt intnl[0] (shear plastic strain)
310  all_q[0].df_di[0] = p * _tan_phi.derivative(intnl[0]) - _cohesion.derivative(intnl[0]);
311  all_q[1].df_di[0] = 0.0;
312  all_q[2].df_di[0] = 0.0;
313  // derivatives wrt intnl[q] (tensile plastic strain)
314  all_q[0].df_di[1] = 0.0;
315  all_q[1].df_di[1] = -_tstrength.derivative(intnl[1]);
316  all_q[2].df_di[1] = -_cstrength.derivative(intnl[1]);
317 
318  // d(flowPotential)/d(p, q)
319  // derivatives wrt p
320  all_q[0].dg[0] = _tan_psi.value(intnl[0]);
321  all_q[1].dg[0] = 1.0;
322  all_q[2].dg[0] = -1.0;
323  // derivatives wrt q
324  if (_small_smoother2 == 0.0)
325  all_q[0].dg[1] = 1.0;
326  else
327  all_q[0].dg[1] = q / std::sqrt(Utility::pow<2>(q) + _small_smoother2);
328  all_q[1].dg[1] = 0.0;
329  all_q[2].dg[1] = 0.0;
330 
331  // d2(flowPotential)/d(p, q)/d(intnl)
332  // d(dg/dp)/dintnl[0]
333  all_q[0].d2g_di[0][0] = _tan_psi.derivative(intnl[0]);
334  all_q[1].d2g_di[0][0] = 0.0;
335  all_q[2].d2g_di[0][0] = 0.0;
336  // d(dg/dp)/dintnl[1]
337  all_q[0].d2g_di[0][1] = 0.0;
338  all_q[1].d2g_di[0][1] = 0.0;
339  all_q[2].d2g_di[0][1] = 0.0;
340  // d(dg/dq)/dintnl[0]
341  all_q[0].d2g_di[1][0] = 0.0;
342  all_q[1].d2g_di[1][0] = 0.0;
343  all_q[2].d2g_di[1][0] = 0.0;
344  // d(dg/dq)/dintnl[1]
345  all_q[0].d2g_di[1][1] = 0.0;
346  all_q[1].d2g_di[1][1] = 0.0;
347  all_q[2].d2g_di[1][1] = 0.0;
348 
349  // d2(flowPotential)/d(p, q)/d(p, q)
350  // d(dg/dp)/dp
351  all_q[0].d2g[0][0] = 0.0;
352  all_q[1].d2g[0][0] = 0.0;
353  all_q[2].d2g[0][0] = 0.0;
354  // d(dg/dp)/dq
355  all_q[0].d2g[0][1] = 0.0;
356  all_q[1].d2g[0][1] = 0.0;
357  all_q[2].d2g[0][1] = 0.0;
358  // d(dg/dq)/dp
359  all_q[0].d2g[1][0] = 0.0;
360  all_q[1].d2g[1][0] = 0.0;
361  all_q[2].d2g[1][0] = 0.0;
362  // d(dg/dq)/dq
363  if (_small_smoother2 == 0.0)
364  all_q[0].d2g[1][1] = 0.0;
365  else
366  all_q[0].d2g[1][1] = _small_smoother2 / std::pow(Utility::pow<2>(q) + _small_smoother2, 1.5);
367  all_q[1].d2g[1][1] = 0.0;
368  all_q[2].d2g[1][1] = 0.0;
369 }
370 
371 void
373  Real q_trial,
374  const std::vector<Real> & intnl_old,
375  Real & p,
376  Real & q,
377  Real & gaE,
378  std::vector<Real> & intnl) const
379 {
380  const Real tanpsi = _tan_psi.value(intnl_old[0]);
381 
382  if (!_perfect_guess)
383  {
384  p = p_trial;
385  q = q_trial;
386  gaE = 0.0;
387  }
388  else
389  {
390  const Real coh = _cohesion.value(intnl_old[0]);
391  const Real tanphi = _tan_phi.value(intnl_old[0]);
392  const Real tens = _tstrength.value(intnl_old[1]);
393  const Real comp = _cstrength.value(intnl_old[1]);
394  const Real q_at_T = coh - tens * tanphi;
395  const Real q_at_C = coh + comp * tanphi;
396 
397  if ((p_trial >= tens) && (q_trial <= q_at_T))
398  {
399  // pure tensile failure
400  p = tens;
401  q = q_trial;
402  gaE = p_trial - p;
403  }
404  else if ((p_trial <= -comp) && (q_trial <= q_at_C))
405  {
406  // pure compressive failure
407  p = -comp;
408  q = q_trial;
409  gaE = p - p_trial;
410  }
411  else
412  {
413  // shear failure or a mixture
414  // Calculate ga assuming a pure shear failure
415  const Real ga = (q_trial + p_trial * tanphi - coh) / (_Eqq + _Epp * tanphi * tanpsi);
416  if (ga <= 0 && p_trial <= tens && p_trial >= -comp)
417  {
418  // very close to one of the rounded corners: there is no advantage to guessing the
419  // solution, so:
420  p = p_trial;
421  q = q_trial;
422  gaE = 0.0;
423  }
424  else
425  {
426  q = q_trial - _Eqq * ga;
427  if (q <= 0.0 && q_at_T <= 0.0)
428  {
429  // user has set tensile strength so large that it is irrelevant: return to the tip of the
430  // shear cone
431  q = 0.0;
432  p = coh / tanphi;
433  gaE = (p_trial - p) / tanpsi; // just a guess, based on the angle to the corner
434  }
435  else if (q <= q_at_T)
436  {
437  // pure shear is incorrect: mixture of tension and shear is correct
438  q = q_at_T;
439  p = tens;
440  gaE = (p_trial - p) / tanpsi; // just a guess, based on the angle to the corner
441  }
442  else if (q >= q_at_C)
443  {
444  // pure shear is incorrect: mixture of compression and shear is correct
445  q = q_at_C;
446  p = -comp;
447  if (p - p_trial < _Epp * tanpsi * (q_trial - q) / _Eqq)
448  // trial point is sitting almost directly above corner
449  gaE = (q_trial - q) * _Epp / _Eqq;
450  else
451  // trial point is sitting to the left of the corner
452  gaE = (p - p_trial) / tanpsi;
453  }
454  else
455  {
456  // pure shear was correct
457  p = p_trial - _Epp * ga * tanpsi;
458  gaE = ga * _Epp;
459  }
460  }
461  }
462  }
463  setIntnlValues(p_trial, q_trial, p, q, intnl_old, intnl);
464 }
465 
466 void
468  Real q_trial,
469  Real p,
470  Real q,
471  const std::vector<Real> & intnl_old,
472  std::vector<Real> & intnl) const
473 {
474  intnl[0] = intnl_old[0] + (q_trial - q) / _Eqq;
475  const Real tanpsi = _tan_psi.value(intnl[0]);
476  intnl[1] = intnl_old[1] + (p_trial - p) / _Epp - (q_trial - q) * tanpsi / _Eqq;
477 }
478 
479 RankTwoTensor
480 CappedWeakPlaneStressUpdate::dpdstress(const RankTwoTensor & /*stress*/) const
481 {
482  RankTwoTensor deriv = RankTwoTensor();
483  deriv(2, 2) = 1.0;
484  return deriv;
485 }
486 
487 RankFourTensor
488 CappedWeakPlaneStressUpdate::d2pdstress2(const RankTwoTensor & /*stress*/) const
489 {
490  return RankFourTensor();
491 }
492 
493 RankTwoTensor
494 CappedWeakPlaneStressUpdate::dqdstress(const RankTwoTensor & stress) const
495 {
496  RankTwoTensor deriv = RankTwoTensor();
497  const Real q = std::sqrt(Utility::pow<2>(stress(2, 0)) + Utility::pow<2>(stress(2, 1)));
498  if (q > 0.0)
499  {
500  deriv(2, 0) = deriv(0, 2) = 0.5 * stress(2, 0) / q;
501  deriv(2, 1) = deriv(1, 2) = 0.5 * stress(2, 1) / q;
502  }
503  else
504  {
505  // derivative is not defined here. For now i'll set:
506  deriv(2, 0) = deriv(0, 2) = 0.5;
507  deriv(2, 1) = deriv(1, 2) = 0.5;
508  }
509  return deriv;
510 }
511 
512 RankFourTensor
513 CappedWeakPlaneStressUpdate::d2qdstress2(const RankTwoTensor & stress) const
514 {
515  RankFourTensor d2 = RankFourTensor();
516 
517  const Real q = std::sqrt(Utility::pow<2>(stress(2, 0)) + Utility::pow<2>(stress(2, 1)));
518  if (q == 0.0)
519  return d2;
520 
521  RankTwoTensor dq = RankTwoTensor();
522  dq(2, 0) = dq(0, 2) = 0.25 * (stress(2, 0) + stress(0, 2)) / q;
523  dq(2, 1) = dq(1, 2) = 0.25 * (stress(2, 1) + stress(1, 2)) / q;
524 
525  for (unsigned i = 0; i < _tensor_dimensionality; ++i)
526  for (unsigned j = 0; j < _tensor_dimensionality; ++j)
527  for (unsigned k = 0; k < _tensor_dimensionality; ++k)
528  for (unsigned l = 0; l < _tensor_dimensionality; ++l)
529  d2(i, j, k, l) = -dq(i, j) * dq(k, l) / q;
530 
531  d2(0, 2, 0, 2) += 0.25 / q;
532  d2(0, 2, 2, 0) += 0.25 / q;
533  d2(2, 0, 0, 2) += 0.25 / q;
534  d2(2, 0, 2, 0) += 0.25 / q;
535  d2(1, 2, 1, 2) += 0.25 / q;
536  d2(1, 2, 2, 1) += 0.25 / q;
537  d2(2, 1, 1, 2) += 0.25 / q;
538  d2(2, 1, 2, 1) += 0.25 / q;
539 
540  return d2;
541 }
const TensorMechanicsHardeningModel & _cstrength
Hardening model for compressive strength.
virtual RankTwoTensor dpdstress(const RankTwoTensor &stress) const override
d(p)/d(stress) Derived classes must override this
virtual void preReturnMap(Real p_trial, Real q_trial, 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...
const TensorMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
Real _dp_dpt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
Real _dq_dpt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
virtual RankFourTensor d2pdstress2(const RankTwoTensor &stress) const override
d2(p)/d(stress)/d(stress) Derived classes must override this
virtual Real derivative(Real intnl) const
Struct designed to hold info about a single yield function and its derivatives, as well as the flow d...
virtual RankTwoTensor dqdstress(const RankTwoTensor &stress) const override
d(q)/d(stress) Derived classes must override this
TwoParameterPlasticityStressUpdate performs the return-map algorithm and associated stress updates fo...
virtual void consistentTangentOperator(const RankTwoTensor &stress_trial, Real p_trial, Real q_trial, const RankTwoTensor &stress, Real p, Real q, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, bool compute_full_tangent_operator, RankFourTensor &cto) const override
Calculates the consistent tangent operator.
Real _dgaE_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
Real _dgaE_dpt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
MaterialProperty< std::vector< Real > > & _intnl
internal parameters
const TensorMechanicsHardeningModel & _tan_phi
Hardening model for tan(phi)
InputParameters validParams< CappedWeakPlaneStressUpdate >()
Real _in_trial02
trial value of stress(0, 2)
Real _dq_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
virtual void yieldFunctionValues(Real p, Real q, const std::vector< Real > &intnl, std::vector< Real > &yf) const override
Computes the values of the yield functions, given p, q and intnl parameters.
InputParameters validParams< TwoParameterPlasticityStressUpdate >()
virtual void setStressAfterReturn(const RankTwoTensor &stress_trial, Real p_ok, Real q_ok, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, RankTwoTensor &stress) const override
Sets stress from the admissible parameters.
virtual void finalizeReturnProcess(const RankTwoTensor &rotation_increment) override
Derived classes may use this to perform calculations after the return-map process has completed succe...
StressReturnType
This allows some simplification in the return-map process.
const TensorMechanicsHardeningModel & _tan_psi
Hardening model for tan(psi)
const bool _perfect_guess
Initialize the NR proceedure from a guess coming from perfect plasticity.
virtual void computeAllQ(Real p, Real q, const std::vector< Real > &intnl, std::vector< yieldAndFlow > &all_q) const override
Completely fills all_q with correct values.
Real _Epp
elasticity tensor in p direction
virtual Real value(Real intnl) const
const Real _smoothing_tol
Smoothing tolerance: edges of the yield surface get smoothed by this amount.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
Real _Eqq
elasticity tensor in q direction
enum CappedWeakPlaneStressUpdate::StressReturnType _stress_return_type
virtual void setEppEqq(const RankFourTensor &Eijkl, Real &Epp, Real &Eqq) const override
Set Epp and Eqq based on the elasticity tensor Derived classes must override this.
virtual void setIntnlValues(Real p_trial, Real q_trial, Real p, Real q, const std::vector< Real > &intnl_old, std::vector< Real > &intnl) const override
Sets the internal parameters based on the trial values of p and q, their current values, and the old values of the internal parameters.
virtual RankFourTensor d2qdstress2(const RankTwoTensor &stress) const override
d2(q)/d(stress)/d(stress) Derived classes must override this
virtual void initializeVars(Real p_trial, Real q_trial, const std::vector< Real > &intnl_old, Real &p, Real &q, Real &gaE, std::vector< Real > &intnl) const override
Sets (p, q, gaE, intnl) at "good guesses" of the solution to the Return-Map algorithm.
virtual void computePQ(const RankTwoTensor &stress, Real &p, Real &q) const override
Computes p and q, given stress.
CappedWeakPlaneStressUpdate(const InputParameters &parameters)
virtual void initializeReturnProcess() override
Derived classes may use this to perform calculations before any return-map process is performed...
const Real _small_smoother2
The cone vertex is smoothed by this amount.
Real _dp_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
const TensorMechanicsHardeningModel & _tstrength
Hardening model for tensile strength.
static constexpr unsigned _tensor_dimensionality
Internal dimensionality of tensors (currently this is 3 throughout tensor_mechanics) ...
virtual void setIntnlDerivatives(Real p_trial, Real q_trial, Real p, Real q, 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 p and q, their current values, and the old values of the internal parameters.
Real _in_trial12
trial value of stress(1, 2)