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