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