www.mooseframework.org
CappedDruckerPragerStressUpdate.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 Drucker-Prager plasticity stress calculator");
17  params.addRequiredParam<UserObjectName>(
18  "DP_model",
19  "A TensorMechanicsPlasticDruckerPrager UserObject that defines the "
20  "Drucker-Prager parameters (cohesion, friction angle and dilation angle)");
21  params.addRequiredParam<UserObjectName>(
22  "tensile_strength",
23  "A TensorMechanicsHardening UserObject that defines hardening of the "
24  "tensile strength. In physical situations this is positive (and always "
25  "must be greater than negative compressive-strength.");
26  params.addRequiredParam<UserObjectName>(
27  "compressive_strength",
28  "A TensorMechanicsHardening UserObject that defines hardening of the "
29  "compressive strength. In physical situations this is positive.");
30  params.addRequiredRangeCheckedParam<Real>(
31  "tip_smoother",
32  "tip_smoother>=0",
33  "The cone vertex at J2 = 0 will be smoothed by the given "
34  "amount. Typical value is 0.1*cohesion");
35  params.addParam<bool>("perfect_guess",
36  true,
37  "Provide a guess to the Newton-Raphson proceedure "
38  "that is the result from perfect plasticity. With "
39  "severe hardening/softening this may be "
40  "suboptimal.");
41  params.addParam<bool>("small_dilation",
42  true,
43  "If true, and if the trial stress exceeds the "
44  "tensile strength, then the user gaurantees that "
45  "the returned stress will be independent of the "
46  "compressive strength.");
47  return params;
48 }
49 
51  : TwoParameterPlasticityStressUpdate(parameters, 3, 2),
52  _dp(getUserObject<TensorMechanicsPlasticDruckerPrager>("DP_model")),
53  _tstrength(getUserObject<TensorMechanicsHardeningModel>("tensile_strength")),
54  _cstrength(getUserObject<TensorMechanicsHardeningModel>("compressive_strength")),
55  _small_smoother2(Utility::pow<2>(getParam<Real>("tip_smoother"))),
56  _perfect_guess(getParam<bool>("perfect_guess")),
57  _stress_return_type(StressReturnType::nothing_special),
58  _small_dilation(getParam<bool>("small_dilation")),
59  _in_q_trial(0.0)
60 {
61  // With arbitary UserObjects, it is impossible to check everything,
62  // but this will catch the common errors
64  mooseError("CappedDruckerPragerStressUpdate: Tensile strength plus compressive strength must "
65  "be greater than smoothing_tol");
66 }
67 
68 void
70 {
72 }
73 
74 void
75 CappedDruckerPragerStressUpdate::finalizeReturnProcess(const RankTwoTensor & /*rotation_increment*/)
76 {
78 }
79 
80 void
82  Real q_trial,
83  const RankTwoTensor & /*stress_trial*/,
84  const std::vector<Real> & /*intnl_old*/,
85  const std::vector<Real> & yf,
86  const RankFourTensor & /*Eijkl*/)
87 {
88  // If it's obvious, then simplify the return-type
89  if (yf[2] >= 0)
91  else if (_small_dilation && yf[1] >= 0)
93 
94  _in_q_trial = q_trial;
95 }
96 
97 void
98 CappedDruckerPragerStressUpdate::computePQ(const RankTwoTensor & stress, Real & p, Real & q) const
99 {
100  p = stress.trace();
101  q = std::sqrt(stress.secondInvariant());
102 }
103 
104 void
105 CappedDruckerPragerStressUpdate::setEppEqq(const RankFourTensor & Eijkl,
106  Real & Epp,
107  Real & Eqq) const
108 {
109  Epp = Eijkl.sum3x3();
110  Eqq = Eijkl(0, 1, 0, 1);
111 }
112 
113 void
115  Real q_trial,
116  Real /*p*/,
117  Real q,
118  const std::vector<Real> & intnl,
119  std::vector<std::vector<Real>> & dintnl) const
120 {
121  Real tanpsi;
122  _dp.onlyB(intnl[0], _dp.dilation, tanpsi);
123  Real dtanpsi;
124  _dp.donlyB(intnl[0], _dp.dilation, dtanpsi);
125  dintnl[0][0] = 0.0;
126  dintnl[0][1] = -1.0 / _Eqq;
127  dintnl[1][0] = -1.0 / _Epp;
128  dintnl[1][1] = tanpsi / _Eqq - (q_trial - q) * dtanpsi * dintnl[0][1] / _Eqq;
129 }
130 
131 void
133  Real q,
134  const std::vector<Real> & intnl,
135  std::vector<Real> & yf) const
136 {
137  Real aaa;
138  Real bbb;
139  _dp.bothAB(intnl[0], aaa, bbb);
140  yf[0] = std::sqrt(Utility::pow<2>(q) + _small_smoother2) + p * bbb - aaa;
141 
143  yf[1] = std::numeric_limits<Real>::lowest();
144  else
145  yf[1] = p - _tstrength.value(intnl[1]);
146 
148  yf[2] = std::numeric_limits<Real>::lowest();
149  else
150  yf[2] = -p - _cstrength.value(intnl[1]);
151 }
152 
153 void
155  Real q,
156  const std::vector<Real> & intnl,
157  std::vector<yieldAndFlow> & all_q) const
158 {
159  Real aaa;
160  Real bbb;
161  _dp.bothAB(intnl[0], aaa, bbb);
162  Real daaa;
163  Real dbbb;
164  _dp.dbothAB(intnl[0], daaa, dbbb);
165  Real tanpsi;
166  _dp.onlyB(intnl[0], _dp.dilation, tanpsi);
167  Real dtanpsi;
168  _dp.donlyB(intnl[0], _dp.dilation, dtanpsi);
169 
170  // yield function values
171  all_q[0].f = std::sqrt(Utility::pow<2>(q) + _small_smoother2) + p * bbb - aaa;
173  all_q[1].f = std::numeric_limits<Real>::lowest();
174  else
175  all_q[1].f = p - _tstrength.value(intnl[1]);
177  all_q[2].f = std::numeric_limits<Real>::lowest();
178  else
179  all_q[2].f = -p - _cstrength.value(intnl[1]);
180 
181  // d(yield Function)/d(p, q)
182  // derivatives wrt p
183  all_q[0].df[0] = bbb;
184  all_q[1].df[0] = 1.0;
185  all_q[2].df[0] = -1.0;
186 
187  // derivatives wrt q
188  if (_small_smoother2 == 0.0)
189  all_q[0].df[1] = 1.0;
190  else
191  all_q[0].df[1] = q / std::sqrt(Utility::pow<2>(q) + _small_smoother2);
192  all_q[1].df[1] = 0.0;
193  all_q[2].df[1] = 0.0;
194 
195  // d(yield Function)/d(intnl)
196  // derivatives wrt intnl[0] (shear plastic strain)
197  all_q[0].df_di[0] = p * dbbb - daaa;
198  all_q[1].df_di[0] = 0.0;
199  all_q[2].df_di[0] = 0.0;
200  // derivatives wrt intnl[q] (tensile plastic strain)
201  all_q[0].df_di[1] = 0.0;
202  all_q[1].df_di[1] = -_tstrength.derivative(intnl[1]);
203  all_q[2].df_di[1] = -_cstrength.derivative(intnl[1]);
204 
205  // d(flowPotential)/d(p, q)
206  // derivatives wrt p
207  all_q[0].dg[0] = tanpsi;
208  all_q[1].dg[0] = 1.0;
209  all_q[2].dg[0] = -1.0;
210  // derivatives wrt q
211  if (_small_smoother2 == 0.0)
212  all_q[0].dg[1] = 1.0;
213  else
214  all_q[0].dg[1] = q / std::sqrt(Utility::pow<2>(q) + _small_smoother2);
215  all_q[1].dg[1] = 0.0;
216  all_q[2].dg[1] = 0.0;
217 
218  // d2(flowPotential)/d(p, q)/d(intnl)
219  // d(dg/dp)/dintnl[0]
220  all_q[0].d2g_di[0][0] = dtanpsi;
221  all_q[1].d2g_di[0][0] = 0.0;
222  all_q[2].d2g_di[0][0] = 0.0;
223  // d(dg/dp)/dintnl[1]
224  all_q[0].d2g_di[0][1] = 0.0;
225  all_q[1].d2g_di[0][1] = 0.0;
226  all_q[2].d2g_di[0][1] = 0.0;
227  // d(dg/dq)/dintnl[0]
228  all_q[0].d2g_di[1][0] = 0.0;
229  all_q[1].d2g_di[1][0] = 0.0;
230  all_q[2].d2g_di[1][0] = 0.0;
231  // d(dg/dq)/dintnl[1]
232  all_q[0].d2g_di[1][1] = 0.0;
233  all_q[1].d2g_di[1][1] = 0.0;
234  all_q[2].d2g_di[1][1] = 0.0;
235 
236  // d2(flowPotential)/d(p, q)/d(p, q)
237  // d(dg/dp)/dp
238  all_q[0].d2g[0][0] = 0.0;
239  all_q[1].d2g[0][0] = 0.0;
240  all_q[2].d2g[0][0] = 0.0;
241  // d(dg/dp)/dq
242  all_q[0].d2g[0][1] = 0.0;
243  all_q[1].d2g[0][1] = 0.0;
244  all_q[2].d2g[0][1] = 0.0;
245  // d(dg/dq)/dp
246  all_q[0].d2g[1][0] = 0.0;
247  all_q[1].d2g[1][0] = 0.0;
248  all_q[2].d2g[1][0] = 0.0;
249  // d(dg/dq)/dq
250  if (_small_smoother2 == 0.0)
251  all_q[0].d2g[1][1] = 0.0;
252  else
253  all_q[0].d2g[1][1] = _small_smoother2 / std::pow(Utility::pow<2>(q) + _small_smoother2, 1.5);
254  all_q[1].d2g[1][1] = 0.0;
255  all_q[2].d2g[1][1] = 0.0;
256 }
257 
258 void
260  Real q_trial,
261  const std::vector<Real> & intnl_old,
262  Real & p,
263  Real & q,
264  Real & gaE,
265  std::vector<Real> & intnl) const
266 {
267  if (!_perfect_guess)
268  {
269  p = p_trial;
270  q = q_trial;
271  gaE = 0.0;
272  }
273  else
274  {
275  Real coh;
276  Real tanphi;
277  _dp.bothAB(intnl[0], coh, tanphi);
278  Real tanpsi;
279  _dp.onlyB(intnl_old[0], _dp.dilation, tanpsi);
280  const Real tens = _tstrength.value(intnl_old[1]);
281  const Real comp = _cstrength.value(intnl_old[1]);
282  const Real q_at_T = coh - tens * tanphi;
283  const Real q_at_C = coh + comp * tanphi;
284 
285  if ((p_trial >= tens) && (q_trial <= q_at_T))
286  {
287  // pure tensile failure
288  p = tens;
289  q = q_trial;
290  gaE = p_trial - p;
291  }
292  else if ((p_trial <= -comp) && (q_trial <= q_at_C))
293  {
294  // pure compressive failure
295  p = -comp;
296  q = q_trial;
297  gaE = p - p_trial;
298  }
299  else
300  {
301  // shear failure or a mixture
302  // Calculate ga assuming a pure shear failure
303  const Real ga = (q_trial + p_trial * tanphi - coh) / (_Eqq + _Epp * tanphi * tanpsi);
304  if (ga <= 0 && p_trial <= tens && p_trial >= -comp)
305  {
306  // very close to one of the rounded corners: there is no advantage to guessing the
307  // solution, so:
308  p = p_trial;
309  q = q_trial;
310  gaE = 0.0;
311  }
312  else
313  {
314  q = q_trial - _Eqq * ga;
315  if (q <= 0.0 && q_at_T <= 0.0)
316  {
317  // user has set tensile strength so large that it is irrelevant: return to the tip of the
318  // shear cone
319  q = 0.0;
320  p = coh / tanphi;
321  gaE = (p_trial - p) / tanpsi; // just a guess, based on the angle to the corner
322  }
323  else if (q <= q_at_T)
324  {
325  // pure shear is incorrect: mixture of tension and shear is correct
326  q = q_at_T;
327  p = tens;
328  gaE = (p_trial - p) / tanpsi; // just a guess, based on the angle to the corner
329  }
330  else if (q >= q_at_C)
331  {
332  // pure shear is incorrect: mixture of compression and shear is correct
333  q = q_at_C;
334  p = -comp;
335  if (p - p_trial < _Epp * tanpsi * (q_trial - q) / _Eqq)
336  // trial point is sitting almost directly above corner
337  gaE = (q_trial - q) * _Epp / _Eqq;
338  else
339  // trial point is sitting to the left of the corner
340  gaE = (p - p_trial) / tanpsi;
341  }
342  else
343  {
344  // pure shear was correct
345  p = p_trial - _Epp * ga * tanpsi;
346  gaE = ga * _Epp;
347  }
348  }
349  }
350  }
351  setIntnlValues(p_trial, q_trial, p, q, intnl_old, intnl);
352 }
353 
354 void
356  Real q_trial,
357  Real p,
358  Real q,
359  const std::vector<Real> & intnl_old,
360  std::vector<Real> & intnl) const
361 {
362  intnl[0] = intnl_old[0] + (q_trial - q) / _Eqq;
363  Real tanpsi;
364  _dp.onlyB(intnl[0], _dp.dilation, tanpsi);
365  intnl[1] = intnl_old[1] + (p_trial - p) / _Epp - (q_trial - q) * tanpsi / _Eqq;
366 }
367 
368 void
369 CappedDruckerPragerStressUpdate::setStressAfterReturn(const RankTwoTensor & stress_trial,
370  Real p_ok,
371  Real q_ok,
372  Real /*gaE*/,
373  const std::vector<Real> & /*intnl*/,
374  const yieldAndFlow & /*smoothed_q*/,
375  const RankFourTensor & /*Eijkl*/,
376  RankTwoTensor & stress) const
377 {
378  // stress = s_ij + de_ij tr(stress) / 3 = q / q_trial * s_ij^trial + de_ij p / 3 = q / q_trial *
379  // (stress_ij^trial - de_ij tr(stress^trial) / 3) + de_ij p / 3
380  const Real p_trial = stress_trial.trace();
381  stress = RankTwoTensor(RankTwoTensor::initIdentity) / 3.0 *
382  (p_ok - (_in_q_trial == 0.0 ? 0.0 : p_trial * q_ok / _in_q_trial));
383  if (_in_q_trial > 0)
384  stress += q_ok / _in_q_trial * stress_trial;
385 }
386 
387 RankTwoTensor
388 CappedDruckerPragerStressUpdate::dpdstress(const RankTwoTensor & stress) const
389 {
390  return stress.dtrace();
391 }
392 
393 RankFourTensor
394 CappedDruckerPragerStressUpdate::d2pdstress2(const RankTwoTensor & /*stress*/) const
395 {
396  return RankFourTensor();
397 }
398 
399 RankTwoTensor
400 CappedDruckerPragerStressUpdate::dqdstress(const RankTwoTensor & stress) const
401 {
402  const Real j2 = stress.secondInvariant();
403  if (j2 == 0.0)
404  return RankTwoTensor();
405  return 0.5 * stress.dsecondInvariant() / std::sqrt(j2);
406 }
407 
408 RankFourTensor
409 CappedDruckerPragerStressUpdate::d2qdstress2(const RankTwoTensor & stress) const
410 {
411  const Real j2 = stress.secondInvariant();
412  if (j2 == 0.0)
413  return RankFourTensor();
414 
415  const RankTwoTensor dj2 = stress.dsecondInvariant();
416  return -0.25 * dj2.outerProduct(dj2) / std::pow(j2, 1.5) +
417  0.5 * stress.d2secondInvariant() / std::sqrt(j2);
418 }
419 
420 void
421 CappedDruckerPragerStressUpdate::consistentTangentOperator(const RankTwoTensor & /*stress_trial*/,
422  Real /*p_trial*/,
423  Real /*q_trial*/,
424  const RankTwoTensor & stress,
425  Real /*p*/,
426  Real q,
427  Real gaE,
428  const yieldAndFlow & smoothed_q,
429  const RankFourTensor & Eijkl,
430  bool compute_full_tangent_operator,
431  RankFourTensor & cto) const
432 {
433  cto = Eijkl;
434  if (!compute_full_tangent_operator)
435  return;
436 
437  const RankTwoTensor s_over_q =
438  (q == 0.0 ? RankTwoTensor()
439  : (stress - stress.trace() * RankTwoTensor(RankTwoTensor::initIdentity) / 3.0) / q);
440 
441  for (unsigned i = 0; i < _tensor_dimensionality; ++i)
442  for (unsigned j = 0; j < _tensor_dimensionality; ++j)
443  for (unsigned k = 0; k < _tensor_dimensionality; ++k)
444  for (unsigned l = 0; l < _tensor_dimensionality; ++l)
445  {
446  cto(i, j, k, l) -= (i == j) * (1.0 / 3.0) * (_Epp * (1.0 - _dp_dpt) / 3.0 * (k == l) +
447  s_over_q(k, l) * _Eqq * (-_dp_dqt));
448  cto(i, j, k, l) -= s_over_q(i, j) * (_Epp * (-_dq_dpt) / 3.0 * (k == l) +
449  s_over_q(k, l) * _Eqq * (1.0 - _dq_dqt));
450  }
451 
452  if (smoothed_q.dg[1] != 0.0)
453  {
454  const RankFourTensor Tijab = Eijkl * (gaE / _Epp) * smoothed_q.dg[1] * d2qdstress2(stress);
455  RankFourTensor inv = RankFourTensor(RankFourTensor::initIdentityFour) + Tijab;
456  try
457  {
458  inv = inv.transposeMajor().invSymm();
459  }
460  catch (const MooseException & e)
461  {
462  // Cannot form the inverse, so probably at some degenerate place in stress space.
463  // Just return with the "best estimate" of the cto.
464  mooseWarning("CappedDruckerPragerStressUpdate: Cannot invert 1+T in consistent tangent "
465  "operator computation at quadpoint ",
466  _qp,
467  " of element ",
468  _current_elem->id());
469  return;
470  }
471  cto = (cto.transposeMajor() * inv).transposeMajor();
472  }
473 }
const Real _small_smoother2
The cone vertex is smoothed by this amount.
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.
virtual RankTwoTensor dpdstress(const RankTwoTensor &stress) const override
d(p)/d(stress) Derived classes must override this
virtual void computePQ(const RankTwoTensor &stress, Real &p, Real &q) const override
Computes p and q, given stress.
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...
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.
Real _dq_dpt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
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...
TwoParameterPlasticityStressUpdate performs the return-map algorithm and associated stress updates fo...
virtual void finalizeReturnProcess(const RankTwoTensor &rotation_increment) override
Derived classes may use this to perform calculations after the return-map process has completed succe...
const bool _small_dilation
If true, and if the trial stress exceeds the tensile strength, then the user gaurantees that the retu...
void bothAB(Real intnl, Real &aaa, Real &bbb) const
Calculates aaa and bbb as a function of the internal parameter intnl.
const TensorMechanicsHardeningModel & _tstrength
Hardening model for tensile strength.
CappedDruckerPragerStressUpdate(const InputParameters &parameters)
const TensorMechanicsHardeningModel & _cstrength
Hardening model for compressive strength.
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.
InputParameters validParams< CappedDruckerPragerStressUpdate >()
Real _dq_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
Rate-independent non-associative Drucker Prager with hardening/softening.
virtual RankFourTensor d2pdstress2(const RankTwoTensor &stress) const override
d2(p)/d(stress)/d(stress) Derived classes must override this
InputParameters validParams< TwoParameterPlasticityStressUpdate >()
void donlyB(Real intnl, int fd, Real &dbbb) const
Calculate d(bbb)/d(intnl) or d(bbb_flow)/d(intnl)
virtual RankTwoTensor dqdstress(const RankTwoTensor &stress) const override
d(q)/d(stress) Derived classes must override this
Real _Epp
elasticity tensor in p direction
virtual RankFourTensor d2qdstress2(const RankTwoTensor &stress) const override
d2(q)/d(stress)/d(stress) Derived classes must override this
void dbothAB(Real intnl, Real &daaa, Real &dbbb) const
Calculates d(aaa)/d(intnl) and d(bbb)/d(intnl) as a function of the internal parameter intnl...
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)
const TensorMechanicsPlasticDruckerPrager & _dp
Hardening model for cohesion, friction and dilation angles.
void onlyB(Real intnl, int fd, Real &bbb) const
Calculate bbb or bbb_flow.
Real _Eqq
elasticity tensor in q direction
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.
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.
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.
StressReturnType
This allows some simplification in the return-map process.
enum CappedDruckerPragerStressUpdate::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.
const bool _perfect_guess
Initialize the NR proceedure from a guess coming from perfect plasticity.
Real _dp_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
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...
static constexpr unsigned _tensor_dimensionality
Internal dimensionality of tensors (currently this is 3 throughout tensor_mechanics) ...