https://mooseframework.inl.gov
Loading...
Searching...
No Matches
SolidMechanicsPlasticMohrCoulombMulti.C
Go to the documentation of this file.
1//* This file is part of the MOOSE framework
2//* https://mooseframework.inl.gov
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 "SolidMechanicsPlasticMohrCoulombMulti.h"
11#include "RankFourTensor.h"
12
13// Following is for perturbing eigvenvalues. This looks really bodgy, but works quite well!
14#include "MooseRandom.h"
15#include "libmesh/utility.h"
16
17registerMooseObject("SolidMechanicsApp", SolidMechanicsPlasticMohrCoulombMulti);
18registerMooseObjectRenamed("SolidMechanicsApp",
19 TensorMechanicsPlasticMohrCoulombMulti,
20 "01/01/2025 00:00",
21 SolidMechanicsPlasticMohrCoulombMulti);
22
23InputParameters
24SolidMechanicsPlasticMohrCoulombMulti::validParams()
25{
26 InputParameters params = SolidMechanicsPlasticModel::validParams();
27 params.addClassDescription("Non-associative Mohr-Coulomb plasticity with hardening/softening");
28 params.addRequiredParam<UserObjectName>(
29 "cohesion", "A SolidMechanicsHardening UserObject that defines hardening of the cohesion");
30 params.addRequiredParam<UserObjectName>("friction_angle",
31 "A SolidMechanicsHardening UserObject "
32 "that defines hardening of the "
33 "friction angle (in radians)");
34 params.addRequiredParam<UserObjectName>("dilation_angle",
35 "A SolidMechanicsHardening UserObject "
36 "that defines hardening of the "
37 "dilation angle (in radians)");
38 params.addParam<unsigned int>("max_iterations",
39 10,
40 "Maximum number of Newton-Raphson iterations "
41 "allowed in the custom return-map algorithm. "
42 " For highly nonlinear hardening this may "
43 "need to be higher than 10.");
44 params.addParam<Real>("shift",
45 "Yield surface is shifted by this amount to avoid problems with "
46 "defining derivatives when eigenvalues are equal. If this is "
47 "larger than f_tol, a warning will be issued. This may be set "
48 "very small when using the custom returnMap. Default = f_tol.");
49 params.addParam<bool>("use_custom_returnMap",
50 true,
51 "Use a custom return-map algorithm for this "
52 "plasticity model, which may speed up "
53 "computations considerably. Set to true "
54 "only for isotropic elasticity with no "
55 "hardening of the dilation angle. In this "
56 "case you may set 'shift' very small.");
57
58 return params;
59}
60
61SolidMechanicsPlasticMohrCoulombMulti::SolidMechanicsPlasticMohrCoulombMulti(
62 const InputParameters & parameters)
63 : SolidMechanicsPlasticModel(parameters),
64 _cohesion(getUserObject<SolidMechanicsHardeningModel>("cohesion")),
65 _phi(getUserObject<SolidMechanicsHardeningModel>("friction_angle")),
66 _psi(getUserObject<SolidMechanicsHardeningModel>("dilation_angle")),
67 _max_iters(getParam<unsigned int>("max_iterations")),
68 _shift(parameters.isParamValid("shift") ? getParam<Real>("shift") : _f_tol),
69 _use_custom_returnMap(getParam<bool>("use_custom_returnMap"))
70{
71 if (_shift < 0)
72 mooseError("Value of 'shift' in SolidMechanicsPlasticMohrCoulombMulti must not be negative\n");
73 if (_shift > _f_tol)
74 _console << "WARNING: value of 'shift' in SolidMechanicsPlasticMohrCoulombMulti is probably "
75 "set too high"
76 << std::endl;
77 if (LIBMESH_DIM != 3)
78 mooseError("SolidMechanicsPlasticMohrCoulombMulti is only defined for LIBMESH_DIM=3");
79 MooseRandom::seed(0);
80}
81
82unsigned int
87
88void
89SolidMechanicsPlasticMohrCoulombMulti::yieldFunctionV(const RankTwoTensor & stress,
90 Real intnl,
91 std::vector<Real> & f) const
92{
93 std::vector<Real> eigvals;
94 stress.symmetricEigenvalues(eigvals);
95 eigvals[0] += _shift;
96 eigvals[2] -= _shift;
97
98 const Real sinphi = std::sin(phi(intnl));
99 const Real cosphi = std::cos(phi(intnl));
100 const Real cohcos = cohesion(intnl) * cosphi;
101
102 yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, f);
103}
104
105void
106SolidMechanicsPlasticMohrCoulombMulti::yieldFunctionEigvals(
107 Real e0, Real e1, Real e2, Real sinphi, Real cohcos, std::vector<Real> & f) const
108{
109 // Naively it seems a shame to have 6 yield functions active instead of just
110 // 3. But 3 won't do. Eg, think of a loading with eigvals[0]=eigvals[1]=eigvals[2]
111 // Then to return to the yield surface would require 2 positive plastic multipliers
112 // and one negative one. Boo hoo.
113
114 f.resize(6);
115 f[0] = 0.5 * (e0 - e1) + 0.5 * (e0 + e1) * sinphi - cohcos;
116 f[1] = 0.5 * (e1 - e0) + 0.5 * (e0 + e1) * sinphi - cohcos;
117 f[2] = 0.5 * (e0 - e2) + 0.5 * (e0 + e2) * sinphi - cohcos;
118 f[3] = 0.5 * (e2 - e0) + 0.5 * (e0 + e2) * sinphi - cohcos;
119 f[4] = 0.5 * (e1 - e2) + 0.5 * (e1 + e2) * sinphi - cohcos;
120 f[5] = 0.5 * (e2 - e1) + 0.5 * (e1 + e2) * sinphi - cohcos;
121}
122
123void
124SolidMechanicsPlasticMohrCoulombMulti::perturbStress(const RankTwoTensor & stress,
125 std::vector<Real> & eigvals,
126 std::vector<RankTwoTensor> & deigvals) const
127{
128 Real small_perturbation;
129 RankTwoTensor shifted_stress = stress;
130 while (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
131 {
132 for (unsigned i = 0; i < 3; ++i)
133 for (unsigned j = 0; j <= i; ++j)
134 {
135 small_perturbation = 0.1 * _shift * 2 * (MooseRandom::rand() - 0.5);
136 shifted_stress(i, j) += small_perturbation;
137 shifted_stress(j, i) += small_perturbation;
138 }
139 shifted_stress.dsymmetricEigenvalues(eigvals, deigvals);
140 }
141}
142
143void
144SolidMechanicsPlasticMohrCoulombMulti::df_dsig(const RankTwoTensor & stress,
145 Real sin_angle,
146 std::vector<RankTwoTensor> & df) const
147{
148 std::vector<Real> eigvals;
149 std::vector<RankTwoTensor> deigvals;
150 stress.dsymmetricEigenvalues(eigvals, deigvals);
151
152 if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
153 perturbStress(stress, eigvals, deigvals);
154
155 df.resize(6);
156 df[0] = 0.5 * (deigvals[0] - deigvals[1]) + 0.5 * (deigvals[0] + deigvals[1]) * sin_angle;
157 df[1] = 0.5 * (deigvals[1] - deigvals[0]) + 0.5 * (deigvals[0] + deigvals[1]) * sin_angle;
158 df[2] = 0.5 * (deigvals[0] - deigvals[2]) + 0.5 * (deigvals[0] + deigvals[2]) * sin_angle;
159 df[3] = 0.5 * (deigvals[2] - deigvals[0]) + 0.5 * (deigvals[0] + deigvals[2]) * sin_angle;
160 df[4] = 0.5 * (deigvals[1] - deigvals[2]) + 0.5 * (deigvals[1] + deigvals[2]) * sin_angle;
161 df[5] = 0.5 * (deigvals[2] - deigvals[1]) + 0.5 * (deigvals[1] + deigvals[2]) * sin_angle;
162}
163
164void
165SolidMechanicsPlasticMohrCoulombMulti::dyieldFunction_dstressV(
166 const RankTwoTensor & stress, Real intnl, std::vector<RankTwoTensor> & df_dstress) const
167{
168 const Real sinphi = std::sin(phi(intnl));
169 df_dsig(stress, sinphi, df_dstress);
170}
171
172void
173SolidMechanicsPlasticMohrCoulombMulti::dyieldFunction_dintnlV(const RankTwoTensor & stress,
174 Real intnl,
175 std::vector<Real> & df_dintnl) const
176{
177 std::vector<Real> eigvals;
178 stress.symmetricEigenvalues(eigvals);
179 eigvals[0] += _shift;
180 eigvals[2] -= _shift;
181
182 const Real sin_angle = std::sin(phi(intnl));
183 const Real cos_angle = std::cos(phi(intnl));
184 const Real dsin_angle = cos_angle * dphi(intnl);
185 const Real dcos_angle = -sin_angle * dphi(intnl);
186 const Real dcohcos = dcohesion(intnl) * cos_angle + cohesion(intnl) * dcos_angle;
187
188 df_dintnl.resize(6);
189 df_dintnl[0] = df_dintnl[1] = 0.5 * (eigvals[0] + eigvals[1]) * dsin_angle - dcohcos;
190 df_dintnl[2] = df_dintnl[3] = 0.5 * (eigvals[0] + eigvals[2]) * dsin_angle - dcohcos;
191 df_dintnl[4] = df_dintnl[5] = 0.5 * (eigvals[1] + eigvals[2]) * dsin_angle - dcohcos;
192}
193
194void
195SolidMechanicsPlasticMohrCoulombMulti::flowPotentialV(const RankTwoTensor & stress,
196 Real intnl,
197 std::vector<RankTwoTensor> & r) const
198{
199 const Real sinpsi = std::sin(psi(intnl));
200 df_dsig(stress, sinpsi, r);
201}
202
203void
204SolidMechanicsPlasticMohrCoulombMulti::dflowPotential_dstressV(
205 const RankTwoTensor & stress, Real intnl, std::vector<RankFourTensor> & dr_dstress) const
206{
207 std::vector<RankFourTensor> d2eigvals;
208 stress.d2symmetricEigenvalues(d2eigvals);
209
210 const Real sinpsi = std::sin(psi(intnl));
211
212 dr_dstress.resize(6);
213 dr_dstress[0] =
214 0.5 * (d2eigvals[0] - d2eigvals[1]) + 0.5 * (d2eigvals[0] + d2eigvals[1]) * sinpsi;
215 dr_dstress[1] =
216 0.5 * (d2eigvals[1] - d2eigvals[0]) + 0.5 * (d2eigvals[0] + d2eigvals[1]) * sinpsi;
217 dr_dstress[2] =
218 0.5 * (d2eigvals[0] - d2eigvals[2]) + 0.5 * (d2eigvals[0] + d2eigvals[2]) * sinpsi;
219 dr_dstress[3] =
220 0.5 * (d2eigvals[2] - d2eigvals[0]) + 0.5 * (d2eigvals[0] + d2eigvals[2]) * sinpsi;
221 dr_dstress[4] =
222 0.5 * (d2eigvals[1] - d2eigvals[2]) + 0.5 * (d2eigvals[1] + d2eigvals[2]) * sinpsi;
223 dr_dstress[5] =
224 0.5 * (d2eigvals[2] - d2eigvals[1]) + 0.5 * (d2eigvals[1] + d2eigvals[2]) * sinpsi;
225}
226
227void
228SolidMechanicsPlasticMohrCoulombMulti::dflowPotential_dintnlV(
229 const RankTwoTensor & stress, Real intnl, std::vector<RankTwoTensor> & dr_dintnl) const
230{
231 const Real cos_angle = std::cos(psi(intnl));
232 const Real dsin_angle = cos_angle * dpsi(intnl);
233
234 std::vector<Real> eigvals;
235 std::vector<RankTwoTensor> deigvals;
236 stress.dsymmetricEigenvalues(eigvals, deigvals);
237
238 if (eigvals[0] > eigvals[1] - 0.1 * _shift || eigvals[1] > eigvals[2] - 0.1 * _shift)
239 perturbStress(stress, eigvals, deigvals);
240
241 dr_dintnl.resize(6);
242 dr_dintnl[0] = dr_dintnl[1] = 0.5 * (deigvals[0] + deigvals[1]) * dsin_angle;
243 dr_dintnl[2] = dr_dintnl[3] = 0.5 * (deigvals[0] + deigvals[2]) * dsin_angle;
244 dr_dintnl[4] = dr_dintnl[5] = 0.5 * (deigvals[1] + deigvals[2]) * dsin_angle;
245}
246
247void
248SolidMechanicsPlasticMohrCoulombMulti::activeConstraints(const std::vector<Real> & f,
249 const RankTwoTensor & stress,
250 Real intnl,
251 const RankFourTensor & Eijkl,
252 std::vector<bool> & act,
253 RankTwoTensor & returned_stress) const
254{
255 act.assign(6, false);
256
257 if (f[0] <= _f_tol && f[1] <= _f_tol && f[2] <= _f_tol && f[3] <= _f_tol && f[4] <= _f_tol &&
258 f[5] <= _f_tol)
259 {
260 returned_stress = stress;
261 return;
262 }
263
264 Real returned_intnl;
265 std::vector<Real> dpm(6);
266 RankTwoTensor delta_dp;
267 std::vector<Real> yf(6);
268 bool trial_stress_inadmissible;
269 doReturnMap(stress,
270 intnl,
271 Eijkl,
272 0.0,
273 returned_stress,
274 returned_intnl,
275 dpm,
276 delta_dp,
277 yf,
278 trial_stress_inadmissible);
279
280 for (unsigned i = 0; i < 6; ++i)
281 act[i] = (dpm[i] > 0);
282}
283
284Real
285SolidMechanicsPlasticMohrCoulombMulti::cohesion(const Real internal_param) const
286{
287 return _cohesion.value(internal_param);
288}
289
290Real
291SolidMechanicsPlasticMohrCoulombMulti::dcohesion(const Real internal_param) const
292{
293 return _cohesion.derivative(internal_param);
294}
295
296Real
297SolidMechanicsPlasticMohrCoulombMulti::phi(const Real internal_param) const
298{
299 return _phi.value(internal_param);
300}
301
302Real
303SolidMechanicsPlasticMohrCoulombMulti::dphi(const Real internal_param) const
304{
305 return _phi.derivative(internal_param);
306}
307
308Real
309SolidMechanicsPlasticMohrCoulombMulti::psi(const Real internal_param) const
310{
311 return _psi.value(internal_param);
312}
313
314Real
315SolidMechanicsPlasticMohrCoulombMulti::dpsi(const Real internal_param) const
316{
317 return _psi.derivative(internal_param);
318}
319
320std::string
321SolidMechanicsPlasticMohrCoulombMulti::modelName() const
322{
323 return "MohrCoulombMulti";
324}
325
326bool
327SolidMechanicsPlasticMohrCoulombMulti::returnMap(const RankTwoTensor & trial_stress,
328 Real intnl_old,
329 const RankFourTensor & E_ijkl,
330 Real ep_plastic_tolerance,
331 RankTwoTensor & returned_stress,
332 Real & returned_intnl,
333 std::vector<Real> & dpm,
334 RankTwoTensor & delta_dp,
335 std::vector<Real> & yf,
336 bool & trial_stress_inadmissible) const
337{
338 if (!_use_custom_returnMap)
339 return SolidMechanicsPlasticModel::returnMap(trial_stress,
340 intnl_old,
341 E_ijkl,
342 ep_plastic_tolerance,
343 returned_stress,
344 returned_intnl,
345 dpm,
346 delta_dp,
347 yf,
348 trial_stress_inadmissible);
349
350 return doReturnMap(trial_stress,
351 intnl_old,
352 E_ijkl,
353 ep_plastic_tolerance,
354 returned_stress,
355 returned_intnl,
356 dpm,
357 delta_dp,
358 yf,
359 trial_stress_inadmissible);
360}
361
362bool
363SolidMechanicsPlasticMohrCoulombMulti::doReturnMap(const RankTwoTensor & trial_stress,
364 Real intnl_old,
365 const RankFourTensor & E_ijkl,
366 Real ep_plastic_tolerance,
367 RankTwoTensor & returned_stress,
368 Real & returned_intnl,
369 std::vector<Real> & dpm,
370 RankTwoTensor & delta_dp,
371 std::vector<Real> & yf,
372 bool & trial_stress_inadmissible) const
373{
374 mooseAssert(dpm.size() == 6,
375 "SolidMechanicsPlasticMohrCoulombMulti size of dpm should be 6 but it is "
376 << dpm.size());
377
378 std::vector<Real> eigvals;
379 RankTwoTensor eigvecs;
380 trial_stress.symmetricEigenvaluesEigenvectors(eigvals, eigvecs);
381 eigvals[0] += _shift;
382 eigvals[2] -= _shift;
383
384 Real sinphi = std::sin(phi(intnl_old));
385 Real cosphi = std::cos(phi(intnl_old));
386 Real coh = cohesion(intnl_old);
387 Real cohcos = coh * cosphi;
388
389 yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, yf);
390
391 if (yf[0] <= _f_tol && yf[1] <= _f_tol && yf[2] <= _f_tol && yf[3] <= _f_tol && yf[4] <= _f_tol &&
392 yf[5] <= _f_tol)
393 {
394 // purely elastic (trial_stress, intnl_old)
395 trial_stress_inadmissible = false;
396 return true;
397 }
398
399 trial_stress_inadmissible = true;
400 delta_dp.zero();
401 returned_stress = RankTwoTensor();
402
403 // these are the normals to the 6 yield surfaces, which are const because of the assumption of no
404 // psi hardening
405 std::vector<RealVectorValue> norm(6);
406 const Real sinpsi = std::sin(psi(intnl_old));
407 const Real oneminus = 0.5 * (1 - sinpsi);
408 const Real oneplus = 0.5 * (1 + sinpsi);
409 norm[0](0) = oneplus;
410 norm[0](1) = -oneminus;
411 norm[0](2) = 0;
412 norm[1](0) = -oneminus;
413 norm[1](1) = oneplus;
414 norm[1](2) = 0;
415 norm[2](0) = oneplus;
416 norm[2](1) = 0;
417 norm[2](2) = -oneminus;
418 norm[3](0) = -oneminus;
419 norm[3](1) = 0;
420 norm[3](2) = oneplus;
421 norm[4](0) = 0;
422 norm[4](1) = oneplus;
423 norm[4](2) = -oneminus;
424 norm[5](0) = 0;
425 norm[5](1) = -oneminus;
426 norm[5](2) = oneplus;
427
428 // the flow directions are these norm multiplied by Eijkl.
429 // I call the flow directions "n".
430 // In the following I assume that the Eijkl is
431 // for an isotropic situation. Then I don't have to
432 // rotate to the principal-stress frame, and i don't
433 // have to worry about strange off-diagonal things
434 std::vector<RealVectorValue> n(6);
435 for (unsigned ys = 0; ys < 6; ++ys)
436 for (unsigned i = 0; i < 3; ++i)
437 for (unsigned j = 0; j < 3; ++j)
438 n[ys](i) += E_ijkl(i, i, j, j) * norm[ys](j);
439 const Real mag_E = E_ijkl(0, 0, 0, 0);
440
441 // With non-zero Poisson's ratio and hardening
442 // it is not computationally cheap to know whether
443 // the trial stress will return to the tip, edge,
444 // or plane. The following at least
445 // gives a not-completely-stupid guess
446 // trial_order[0] = type of return to try first
447 // trial_order[1] = type of return to try second
448 // trial_order[2] = type of return to try third
449 // trial_order[3] = type of return to try fourth
450 // trial_order[4] = type of return to try fifth
451 // In the following the "binary" stuff indicates the
452 // deactive (0) and active (1) surfaces, eg
453 // 110100 means that surfaces 0, 1 and 3 are active
454 // and 2, 4 and 5 are deactive
455 const unsigned int number_of_return_paths = 5;
456 std::vector<int> trial_order(number_of_return_paths);
457 if (yf[1] > _f_tol && yf[3] > _f_tol && yf[5] > _f_tol)
458 {
459 trial_order[0] = tip110100;
460 trial_order[1] = edge010100;
461 trial_order[2] = plane000100;
462 trial_order[3] = edge000101;
463 trial_order[4] = tip010101;
464 }
465 else if (yf[1] <= _f_tol && yf[3] > _f_tol && yf[5] > _f_tol)
466 {
467 trial_order[0] = edge000101;
468 trial_order[1] = plane000100;
469 trial_order[2] = tip110100;
470 trial_order[3] = tip010101;
471 trial_order[4] = edge010100;
472 }
473 else if (yf[1] <= _f_tol && yf[3] > _f_tol && yf[5] <= _f_tol)
474 {
475 trial_order[0] = plane000100;
476 trial_order[1] = edge000101;
477 trial_order[2] = edge010100;
478 trial_order[3] = tip110100;
479 trial_order[4] = tip010101;
480 }
481 else
482 {
483 trial_order[0] = edge010100;
484 trial_order[1] = plane000100;
485 trial_order[2] = edge000101;
486 trial_order[3] = tip110100;
487 trial_order[4] = tip010101;
488 }
489
490 unsigned trial;
491 bool nr_converged = false;
492 bool kt_success = false;
493 std::vector<RealVectorValue> ntip(3);
494 std::vector<Real> dpmtip(3);
495
496 for (trial = 0; trial < number_of_return_paths; ++trial)
497 {
498 switch (trial_order[trial])
499 {
500 case tip110100:
501 for (unsigned int i = 0; i < 3; ++i)
502 {
503 ntip[0](i) = n[0](i);
504 ntip[1](i) = n[1](i);
505 ntip[2](i) = n[3](i);
506 }
507 kt_success = returnTip(eigvals,
508 ntip,
509 dpmtip,
510 returned_stress,
511 intnl_old,
512 sinphi,
513 cohcos,
514 0,
515 nr_converged,
516 ep_plastic_tolerance,
517 yf);
518 if (nr_converged && kt_success)
519 {
520 dpm[0] = dpmtip[0];
521 dpm[1] = dpmtip[1];
522 dpm[3] = dpmtip[2];
523 dpm[2] = dpm[4] = dpm[5] = 0;
524 }
525 break;
526
527 case tip010101:
528 for (unsigned int i = 0; i < 3; ++i)
529 {
530 ntip[0](i) = n[1](i);
531 ntip[1](i) = n[3](i);
532 ntip[2](i) = n[5](i);
533 }
534 kt_success = returnTip(eigvals,
535 ntip,
536 dpmtip,
537 returned_stress,
538 intnl_old,
539 sinphi,
540 cohcos,
541 0,
542 nr_converged,
543 ep_plastic_tolerance,
544 yf);
545 if (nr_converged && kt_success)
546 {
547 dpm[1] = dpmtip[0];
548 dpm[3] = dpmtip[1];
549 dpm[5] = dpmtip[2];
550 dpm[0] = dpm[2] = dpm[4] = 0;
551 }
552 break;
553
554 case edge000101:
555 kt_success = returnEdge000101(eigvals,
556 n,
557 dpm,
558 returned_stress,
559 intnl_old,
560 sinphi,
561 cohcos,
562 0,
563 mag_E,
564 nr_converged,
565 ep_plastic_tolerance,
566 yf);
567 break;
568
569 case edge010100:
570 kt_success = returnEdge010100(eigvals,
571 n,
572 dpm,
573 returned_stress,
574 intnl_old,
575 sinphi,
576 cohcos,
577 0,
578 mag_E,
579 nr_converged,
580 ep_plastic_tolerance,
581 yf);
582 break;
583
584 case plane000100:
585 kt_success = returnPlane(eigvals,
586 n,
587 dpm,
588 returned_stress,
589 intnl_old,
590 sinphi,
591 cohcos,
592 0,
593 nr_converged,
594 ep_plastic_tolerance,
595 yf);
596 break;
597 }
598
599 if (nr_converged && kt_success)
600 break;
601 }
602
603 if (trial == number_of_return_paths)
604 {
605 sinphi = std::sin(phi(intnl_old));
606 cosphi = std::cos(phi(intnl_old));
607 coh = cohesion(intnl_old);
608 cohcos = coh * cosphi;
609 yieldFunctionEigvals(eigvals[0], eigvals[1], eigvals[2], sinphi, cohcos, yf);
610 Moose::err << "Trial stress = \n";
611 trial_stress.print(Moose::err);
612 Moose::err << "which has eigenvalues = " << eigvals[0] << " " << eigvals[1] << " " << eigvals[2]
613 << "\n";
614 Moose::err << "and yield functions = " << yf[0] << " " << yf[1] << " " << yf[2] << " " << yf[3]
615 << " " << yf[4] << " " << yf[5] << "\n";
616 Moose::err << "Internal parameter = " << intnl_old << std::endl;
617 mooseError("SolidMechanicsPlasticMohrCoulombMulti: FAILURE! You probably need to implement a "
618 "line search if your hardening is too severe, or you need to tune your tolerances "
619 "(eg, yield_function_tolerance should be a little smaller than (young "
620 "modulus)*ep_plastic_tolerance).\n");
621 return false;
622 }
623
624 // success
625
626 returned_intnl = intnl_old;
627 for (unsigned i = 0; i < 6; ++i)
628 returned_intnl += dpm[i];
629 for (unsigned i = 0; i < 6; ++i)
630 for (unsigned j = 0; j < 3; ++j)
631 delta_dp(j, j) += dpm[i] * norm[i](j);
632 returned_stress = eigvecs * returned_stress * (eigvecs.transpose());
633 delta_dp = eigvecs * delta_dp * (eigvecs.transpose());
634 return true;
635}
636
637bool
638SolidMechanicsPlasticMohrCoulombMulti::returnTip(const std::vector<Real> & eigvals,
639 const std::vector<RealVectorValue> & n,
640 std::vector<Real> & dpm,
641 RankTwoTensor & returned_stress,
642 Real intnl_old,
643 Real & sinphi,
644 Real & cohcos,
645 Real initial_guess,
646 bool & nr_converged,
647 Real ep_plastic_tolerance,
648 std::vector<Real> & yf) const
649{
650 // This returns to the Mohr-Coulomb tip using the THREE directions
651 // given in n, and yields the THREE dpm values. Note that you
652 // must supply THREE suitable n vectors out of the total of SIX
653 // flow directions, and then interpret the THREE dpm values appropriately.
654 //
655 // Eg1. You supply the flow directions n[0], n[1] and n[3] as
656 // the "n" vectors. This is return-to-the-tip via 110100.
657 // Then the three returned dpm values will be dpm[0], dpm[1] and dpm[3].
658
659 // Eg2. You supply the flow directions n[1], n[3] and n[5] as
660 // the "n" vectors. This is return-to-the-tip via 010101.
661 // Then the three returned dpm values will be dpm[1], dpm[3] and dpm[5].
662
663 // The returned point is defined by the three yield functions (corresonding
664 // to the three supplied flow directions) all being zero.
665 // that is, returned_stress = diag(cohcot, cohcot, cohcot), where
666 // cohcot = cohesion*cosphi/sinphi
667 // where intnl = intnl_old + dpm[0] + dpm[1] + dpm[2]
668 // The 3 plastic multipliers, dpm, are defiend by the normality condition
669 // eigvals - cohcot = dpm[0]*n[0] + dpm[1]*n[1] + dpm[2]*n[2]
670 // (Kuhn-Tucker demands that all dpm are non-negative, but we leave
671 // that checking for the end.)
672 // (Remember that these "n" vectors and "dpm" values must be interpreted
673 // differently depending on what you pass into this function.)
674 // This is a vector equation with solution (A):
675 // dpm[0] = triple(eigvals - cohcot, n[1], n[2])/trip;
676 // dpm[1] = triple(eigvals - cohcot, n[2], n[0])/trip;
677 // dpm[2] = triple(eigvals - cohcot, n[0], n[1])/trip;
678 // where trip = triple(n[0], n[1], n[2]).
679 // By adding the three components of that solution together
680 // we can get an equation for x = dpm[0] + dpm[1] + dpm[2],
681 // and then our Newton-Raphson only involves one variable (x).
682 // In the following, i specialise to the isotropic situation.
683
684 mooseAssert(n.size() == 3,
685 "SolidMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
686 "supplied with n of size 3, whereas yours is "
687 << n.size());
688 mooseAssert(dpm.size() == 3,
689 "SolidMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
690 "supplied with dpm of size 3, whereas yours is "
691 << dpm.size());
692 mooseAssert(yf.size() == 6,
693 "SolidMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
694 "supplied with yf of size 6, whereas yours is "
695 << yf.size());
696
697 Real x = initial_guess;
698 const Real trip = triple_product(n[0], n[1], n[2]);
699 sinphi = std::sin(phi(intnl_old + x));
700 Real cosphi = std::cos(phi(intnl_old + x));
701 Real coh = cohesion(intnl_old + x);
702 cohcos = coh * cosphi;
703 Real cohcot = cohcos / sinphi;
704
705 if (_cohesion.modelName().compare("Constant") != 0 || _phi.modelName().compare("Constant") != 0)
706 {
707 // Finding x is expensive. Therefore
708 // although x!=0 for Constant Hardening, solution (A)
709 // demonstrates that we don't
710 // actually need to know x to find the dpm for
711 // Constant Hardening.
712 //
713 // However, for nontrivial Hardening, the following
714 // is necessary
715 // cohcot_coeff = [1,1,1].(Cross[n[1], n[2]] + Cross[n[2], n[0]] + Cross[n[0], n[1]])/trip
716 Real cohcot_coeff =
717 (n[0](0) * (n[1](1) - n[1](2) - n[2](1)) + (n[1](2) - n[1](1)) * n[2](0) +
718 (n[1](0) - n[1](2)) * n[2](1) + n[0](2) * (n[1](0) - n[1](1) - n[2](0) + n[2](1)) +
719 n[0](1) * (n[1](2) - n[1](0) + n[2](0) - n[2](2)) +
720 (n[0](0) - n[1](0) + n[1](1)) * n[2](2)) /
721 trip;
722 // eig_term = eigvals.(Cross[n[1], n[2]] + Cross[n[2], n[0]] + Cross[n[0], n[1]])/trip
723 Real eig_term = eigvals[0] *
724 (-n[0](2) * n[1](1) + n[0](1) * n[1](2) + n[0](2) * n[2](1) -
725 n[1](2) * n[2](1) - n[0](1) * n[2](2) + n[1](1) * n[2](2)) /
726 trip;
727 eig_term += eigvals[1] *
728 (n[0](2) * n[1](0) - n[0](0) * n[1](2) - n[0](2) * n[2](0) + n[1](2) * n[2](0) +
729 n[0](0) * n[2](2) - n[1](0) * n[2](2)) /
730 trip;
731 eig_term += eigvals[2] *
732 (n[0](0) * n[1](1) - n[1](1) * n[2](0) + n[0](1) * n[2](0) - n[0](1) * n[1](0) -
733 n[0](0) * n[2](1) + n[1](0) * n[2](1)) /
734 trip;
735 // and finally, the equation we want to solve is:
736 // x - eig_term + cohcot*cohcot_coeff = 0
737 // but i divide by cohcot_coeff so the result has the units of
738 // stress, so using _f_tol as a convergence check is reasonable
739 eig_term /= cohcot_coeff;
740 Real residual = x / cohcot_coeff - eig_term + cohcot;
741 Real jacobian;
742 Real deriv_phi;
743 Real deriv_coh;
744 unsigned int iter = 0;
745 do
746 {
747 deriv_phi = dphi(intnl_old + x);
748 deriv_coh = dcohesion(intnl_old + x);
749 jacobian = 1.0 / cohcot_coeff + deriv_coh * cosphi / sinphi -
750 coh * deriv_phi / Utility::pow<2>(sinphi);
751 x += -residual / jacobian;
752
753 if (iter > _max_iters) // not converging
754 {
755 nr_converged = false;
756 return false;
757 }
758
759 sinphi = std::sin(phi(intnl_old + x));
760 cosphi = std::cos(phi(intnl_old + x));
761 coh = cohesion(intnl_old + x);
762 cohcos = coh * cosphi;
763 cohcot = cohcos / sinphi;
764 residual = x / cohcot_coeff - eig_term + cohcot;
765 iter++;
766 } while (residual * residual > _f_tol * _f_tol / 100);
767 }
768
769 // so the NR process converged, but we must
770 // calculate the individual dpm values and
771 // check Kuhn-Tucker
772 nr_converged = true;
773 if (x < -3 * ep_plastic_tolerance)
774 // obviously at least one of the dpm are < -ep_plastic_tolerance. No point in proceeding. This
775 // is a potential weak-point: if the user has set _f_tol quite large, and ep_plastic_tolerance
776 // quite small, the above NR process will quickly converge, but the solution may be wrong and
777 // violate Kuhn-Tucker.
778 return false;
779
780 // The following is the solution (A) written above
781 // (dpm[0] = triple(eigvals - cohcot, n[1], n[2])/trip, etc)
782 // in the isotropic situation
783 RealVectorValue v;
784 v(0) = eigvals[0] - cohcot;
785 v(1) = eigvals[1] - cohcot;
786 v(2) = eigvals[2] - cohcot;
787 dpm[0] = triple_product(v, n[1], n[2]) / trip;
788 dpm[1] = triple_product(v, n[2], n[0]) / trip;
789 dpm[2] = triple_product(v, n[0], n[1]) / trip;
790
791 if (dpm[0] < -ep_plastic_tolerance || dpm[1] < -ep_plastic_tolerance ||
792 dpm[2] < -ep_plastic_tolerance)
793 // Kuhn-Tucker failure. No point in proceeding
794 return false;
795
796 // Kuhn-Tucker has succeeded: just need returned_stress and yf values
797 // I do not use the dpm to calculate returned_stress, because that
798 // might add error (and computational effort), simply:
799 returned_stress(0, 0) = returned_stress(1, 1) = returned_stress(2, 2) = cohcot;
800 // So by construction the yield functions are all zero
801 yf[0] = yf[1] = yf[2] = yf[3] = yf[4] = yf[5] = 0;
802 return true;
803}
804
805bool
806SolidMechanicsPlasticMohrCoulombMulti::returnPlane(const std::vector<Real> & eigvals,
807 const std::vector<RealVectorValue> & n,
808 std::vector<Real> & dpm,
809 RankTwoTensor & returned_stress,
810 Real intnl_old,
811 Real & sinphi,
812 Real & cohcos,
813 Real initial_guess,
814 bool & nr_converged,
815 Real ep_plastic_tolerance,
816 std::vector<Real> & yf) const
817{
818 // This returns to the Mohr-Coulomb plane using n[3] (ie 000100)
819 //
820 // The returned point is defined by the f[3]=0 and
821 // a = eigvals - dpm[3]*n[3]
822 // where "a" is the returned point and dpm[3] is the plastic multiplier.
823 // This equation is a vector equation in principal stress space.
824 // (Kuhn-Tucker also demands that dpm[3]>=0, but we leave checking
825 // that condition for the end.)
826 // Since f[3]=0, we must have
827 // a[2]*(1+sinphi) + a[0]*(-1+sinphi) - 2*coh*cosphi = 0
828 // which gives dpm[3] as the solution of
829 // alpha*dpm[3] + eigvals[2] - eigvals[0] + beta*sinphi - 2*coh*cosphi = 0
830 // with alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0))*sinphi
831 // beta = eigvals[2] + eigvals[0]
832
833 mooseAssert(n.size() == 6,
834 "SolidMechanicsPlasticMohrCoulombMulti: Custom plane-return algorithm must be "
835 "supplied with n of size 6, whereas yours is "
836 << n.size());
837 mooseAssert(dpm.size() == 6,
838 "SolidMechanicsPlasticMohrCoulombMulti: Custom plane-return algorithm must be "
839 "supplied with dpm of size 6, whereas yours is "
840 << dpm.size());
841 mooseAssert(yf.size() == 6,
842 "SolidMechanicsPlasticMohrCoulombMulti: Custom tip-return algorithm must be "
843 "supplied with yf of size 6, whereas yours is "
844 << yf.size());
845
846 dpm[3] = initial_guess;
847 sinphi = std::sin(phi(intnl_old + dpm[3]));
848 Real cosphi = std::cos(phi(intnl_old + dpm[3]));
849 Real coh = cohesion(intnl_old + dpm[3]);
850 cohcos = coh * cosphi;
851
852 Real alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0)) * sinphi;
853 Real deriv_phi;
854 Real dalpha;
855 const Real beta = eigvals[2] + eigvals[0];
856 Real deriv_coh;
857
858 Real residual =
859 alpha * dpm[3] + eigvals[2] - eigvals[0] + beta * sinphi - 2.0 * cohcos; // this is 2*yf[3]
860 Real jacobian;
861
862 const Real f_tol2 = Utility::pow<2>(_f_tol);
863 unsigned int iter = 0;
864 do
865 {
866 deriv_phi = dphi(intnl_old + dpm[3]);
867 dalpha = -(n[3](2) + n[3](0)) * cosphi * deriv_phi;
868 deriv_coh = dcohesion(intnl_old + dpm[3]);
869 jacobian = alpha + dalpha * dpm[3] + beta * cosphi * deriv_phi - 2.0 * deriv_coh * cosphi +
870 2.0 * coh * sinphi * deriv_phi;
871
872 dpm[3] -= residual / jacobian;
873 if (iter > _max_iters) // not converging
874 {
875 nr_converged = false;
876 return false;
877 }
878
879 sinphi = std::sin(phi(intnl_old + dpm[3]));
880 cosphi = std::cos(phi(intnl_old + dpm[3]));
881 coh = cohesion(intnl_old + dpm[3]);
882 cohcos = coh * cosphi;
883 alpha = n[3](0) - n[3](2) - (n[3](2) + n[3](0)) * sinphi;
884 residual = alpha * dpm[3] + eigvals[2] - eigvals[0] + beta * sinphi - 2.0 * cohcos;
885 iter++;
886 } while (residual * residual > f_tol2);
887
888 // so the NR process converged, but we must
889 // check Kuhn-Tucker
890 nr_converged = true;
891 if (dpm[3] < -ep_plastic_tolerance)
892 // Kuhn-Tucker failure
893 return false;
894
895 for (unsigned i = 0; i < 3; ++i)
896 returned_stress(i, i) = eigvals[i] - dpm[3] * n[3](i);
897 yieldFunctionEigvals(
898 returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
899
900 // by construction abs(yf[3]) = abs(residual/2) < _f_tol/2
901 if (yf[0] > _f_tol || yf[1] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol || yf[5] > _f_tol)
902 // Kuhn-Tucker failure
903 return false;
904
905 // success!
906 dpm[0] = dpm[1] = dpm[2] = dpm[4] = dpm[5] = 0;
907 return true;
908}
909
910bool
911SolidMechanicsPlasticMohrCoulombMulti::returnEdge000101(const std::vector<Real> & eigvals,
912 const std::vector<RealVectorValue> & n,
913 std::vector<Real> & dpm,
914 RankTwoTensor & returned_stress,
915 Real intnl_old,
916 Real & sinphi,
917 Real & cohcos,
918 Real initial_guess,
919 Real mag_E,
920 bool & nr_converged,
921 Real ep_plastic_tolerance,
922 std::vector<Real> & yf) const
923{
924 // This returns to the Mohr-Coulomb edge
925 // with 000101 being active. This means that
926 // f3=f5=0. Denoting the returned stress by "a"
927 // (in principal stress space), this means that
928 // a0=a1 = (2Ccosphi + a2(1+sinphi))/(sinphi-1)
929 //
930 // Also, a = eigvals - dpm3*n3 - dpm5*n5,
931 // where the dpm are the plastic multipliers
932 // and the n are the flow directions.
933 //
934 // Hence we have 5 equations and 5 unknowns (a,
935 // dpm3 and dpm5). Eliminating all unknowns
936 // except for x = dpm3+dpm5 gives the residual below
937 // (I used mathematica)
938
939 Real x = initial_guess;
940 sinphi = std::sin(phi(intnl_old + x));
941 Real cosphi = std::cos(phi(intnl_old + x));
942 Real coh = cohesion(intnl_old + x);
943 cohcos = coh * cosphi;
944
945 const Real numer_const =
946 -n[3](2) * eigvals[0] - n[5](1) * eigvals[0] + n[5](2) * eigvals[0] + n[3](2) * eigvals[1] +
947 n[5](0) * eigvals[1] - n[5](2) * eigvals[1] - n[5](0) * eigvals[2] + n[5](1) * eigvals[2] +
948 n[3](0) * (-eigvals[1] + eigvals[2]) - n[3](1) * (-eigvals[0] + eigvals[2]);
949 const Real numer_coeff1 = 2 * (-n[3](0) + n[3](1) + n[5](0) - n[5](1));
950 const Real numer_coeff2 =
951 n[5](2) * (eigvals[0] - eigvals[1]) + n[3](2) * (-eigvals[0] + eigvals[1]) +
952 n[5](1) * (eigvals[0] + eigvals[2]) + (n[3](0) - n[5](0)) * (eigvals[1] + eigvals[2]) -
953 n[3](1) * (eigvals[0] + eigvals[2]);
954 Real numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
955 const Real denom_const = -n[3](2) * (n[5](0) - n[5](1)) - n[3](1) * (-n[5](0) + n[5](2)) +
956 n[3](0) * (-n[5](1) + n[5](2));
957 const Real denom_coeff = -n[3](2) * (n[5](0) - n[5](1)) - n[3](1) * (n[5](0) + n[5](2)) +
958 n[3](0) * (n[5](1) + n[5](2));
959 Real denom = denom_const + denom_coeff * sinphi;
960 Real residual = -x + numer / denom;
961
962 Real deriv_phi;
963 Real deriv_coh;
964 Real jacobian;
965 const Real tol = Utility::pow<2>(_f_tol / (mag_E * 10.0));
966 unsigned int iter = 0;
967 do
968 {
969 do
970 {
971 deriv_phi = dphi(intnl_old + dpm[3]);
972 deriv_coh = dcohesion(intnl_old + dpm[3]);
973 jacobian = -1 +
974 (numer_coeff1 * deriv_coh * cosphi - numer_coeff1 * coh * sinphi * deriv_phi +
975 numer_coeff2 * cosphi * deriv_phi) /
976 denom -
977 numer * denom_coeff * cosphi * deriv_phi / denom / denom;
978 x -= residual / jacobian;
979 if (iter > _max_iters) // not converging
980 {
981 nr_converged = false;
982 return false;
983 }
984
985 sinphi = std::sin(phi(intnl_old + x));
986 cosphi = std::cos(phi(intnl_old + x));
987 coh = cohesion(intnl_old + x);
988 cohcos = coh * cosphi;
989 numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
990 denom = denom_const + denom_coeff * sinphi;
991 residual = -x + numer / denom;
992 iter++;
993 } while (residual * residual > tol);
994
995 // now must ensure that yf[3] and yf[5] are both "zero"
996 const Real dpm3minusdpm5 =
997 (2.0 * (eigvals[0] - eigvals[1]) + x * (n[3](1) - n[3](0) + n[5](1) - n[5](0))) /
998 (n[3](0) - n[3](1) + n[5](1) - n[5](0));
999 dpm[3] = (x + dpm3minusdpm5) / 2.0;
1000 dpm[5] = (x - dpm3minusdpm5) / 2.0;
1001
1002 for (unsigned i = 0; i < 3; ++i)
1003 returned_stress(i, i) = eigvals[i] - dpm[3] * n[3](i) - dpm[5] * n[5](i);
1004 yieldFunctionEigvals(
1005 returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
1006 } while (yf[3] * yf[3] > _f_tol * _f_tol && yf[5] * yf[5] > _f_tol * _f_tol);
1007
1008 // so the NR process converged, but we must
1009 // check Kuhn-Tucker
1010 nr_converged = true;
1011
1012 if (dpm[3] < -ep_plastic_tolerance || dpm[5] < -ep_plastic_tolerance)
1013 // Kuhn-Tucker failure. This is a potential weak-point: if the user has set _f_tol quite
1014 // large, and ep_plastic_tolerance quite small, the above NR process will quickly converge, but
1015 // the solution may be wrong and violate Kuhn-Tucker.
1016 return false;
1017
1018 if (yf[0] > _f_tol || yf[1] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol)
1019 // Kuhn-Tucker failure
1020 return false;
1021
1022 // success
1023 dpm[0] = dpm[1] = dpm[2] = dpm[4] = 0;
1024 return true;
1025}
1026
1027bool
1028SolidMechanicsPlasticMohrCoulombMulti::returnEdge010100(const std::vector<Real> & eigvals,
1029 const std::vector<RealVectorValue> & n,
1030 std::vector<Real> & dpm,
1031 RankTwoTensor & returned_stress,
1032 Real intnl_old,
1033 Real & sinphi,
1034 Real & cohcos,
1035 Real initial_guess,
1036 Real mag_E,
1037 bool & nr_converged,
1038 Real ep_plastic_tolerance,
1039 std::vector<Real> & yf) const
1040{
1041 // This returns to the Mohr-Coulomb edge
1042 // with 010100 being active. This means that
1043 // f1=f3=0. Denoting the returned stress by "a"
1044 // (in principal stress space), this means that
1045 // a1=a2 = (2Ccosphi + a0(1-sinphi))/(sinphi+1)
1046 //
1047 // Also, a = eigvals - dpm1*n1 - dpm3*n3,
1048 // where the dpm are the plastic multipliers
1049 // and the n are the flow directions.
1050 //
1051 // Hence we have 5 equations and 5 unknowns (a,
1052 // dpm1 and dpm3). Eliminating all unknowns
1053 // except for x = dpm1+dpm3 gives the residual below
1054 // (I used mathematica)
1055
1056 Real x = initial_guess;
1057 sinphi = std::sin(phi(intnl_old + x));
1058 Real cosphi = std::cos(phi(intnl_old + x));
1059 Real coh = cohesion(intnl_old + x);
1060 cohcos = coh * cosphi;
1061
1062 const Real numer_const = -n[1](2) * eigvals[0] - n[3](1) * eigvals[0] + n[3](2) * eigvals[0] -
1063 n[1](0) * eigvals[1] + n[1](2) * eigvals[1] + n[3](0) * eigvals[1] -
1064 n[3](2) * eigvals[1] + n[1](0) * eigvals[2] - n[3](0) * eigvals[2] +
1065 n[3](1) * eigvals[2] - n[1](1) * (-eigvals[0] + eigvals[2]);
1066 const Real numer_coeff1 = 2 * (n[1](1) - n[1](2) - n[3](1) + n[3](2));
1067 const Real numer_coeff2 =
1068 n[3](2) * (-eigvals[0] - eigvals[1]) + n[1](2) * (eigvals[0] + eigvals[1]) +
1069 n[3](1) * (eigvals[0] + eigvals[2]) + (n[1](0) - n[3](0)) * (eigvals[1] - eigvals[2]) -
1070 n[1](1) * (eigvals[0] + eigvals[2]);
1071 Real numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
1072 const Real denom_const = -n[1](0) * (n[3](1) - n[3](2)) + n[1](2) * (-n[3](0) + n[3](1)) +
1073 n[1](1) * (-n[3](2) + n[3](0));
1074 const Real denom_coeff =
1075 n[1](0) * (n[3](1) - n[3](2)) + n[1](2) * (n[3](0) + n[3](1)) - n[1](1) * (n[3](0) + n[3](2));
1076 Real denom = denom_const + denom_coeff * sinphi;
1077 Real residual = -x + numer / denom;
1078
1079 Real deriv_phi;
1080 Real deriv_coh;
1081 Real jacobian;
1082 const Real tol = Utility::pow<2>(_f_tol / (mag_E * 10.0));
1083 unsigned int iter = 0;
1084 do
1085 {
1086 do
1087 {
1088 deriv_phi = dphi(intnl_old + dpm[3]);
1089 deriv_coh = dcohesion(intnl_old + dpm[3]);
1090 jacobian = -1 +
1091 (numer_coeff1 * deriv_coh * cosphi - numer_coeff1 * coh * sinphi * deriv_phi +
1092 numer_coeff2 * cosphi * deriv_phi) /
1093 denom -
1094 numer * denom_coeff * cosphi * deriv_phi / denom / denom;
1095 x -= residual / jacobian;
1096 if (iter > _max_iters) // not converging
1097 {
1098 nr_converged = false;
1099 return false;
1100 }
1101
1102 sinphi = std::sin(phi(intnl_old + x));
1103 cosphi = std::cos(phi(intnl_old + x));
1104 coh = cohesion(intnl_old + x);
1105 cohcos = coh * cosphi;
1106 numer = numer_const + numer_coeff1 * cohcos + numer_coeff2 * sinphi;
1107 denom = denom_const + denom_coeff * sinphi;
1108 residual = -x + numer / denom;
1109 iter++;
1110 } while (residual * residual > tol);
1111
1112 // now must ensure that yf[1] and yf[3] are both "zero"
1113 Real dpm1minusdpm3 =
1114 (2 * (eigvals[1] - eigvals[2]) + x * (n[1](2) - n[1](1) + n[3](2) - n[3](1))) /
1115 (n[1](1) - n[1](2) + n[3](2) - n[3](1));
1116 dpm[1] = (x + dpm1minusdpm3) / 2.0;
1117 dpm[3] = (x - dpm1minusdpm3) / 2.0;
1118
1119 for (unsigned i = 0; i < 3; ++i)
1120 returned_stress(i, i) = eigvals[i] - dpm[1] * n[1](i) - dpm[3] * n[3](i);
1121 yieldFunctionEigvals(
1122 returned_stress(0, 0), returned_stress(1, 1), returned_stress(2, 2), sinphi, cohcos, yf);
1123 } while (yf[1] * yf[1] > _f_tol * _f_tol && yf[3] * yf[3] > _f_tol * _f_tol);
1124
1125 // so the NR process converged, but we must
1126 // check Kuhn-Tucker
1127 nr_converged = true;
1128
1129 if (dpm[1] < -ep_plastic_tolerance || dpm[3] < -ep_plastic_tolerance)
1130 // Kuhn-Tucker failure. This is a potential weak-point: if the user has set _f_tol quite
1131 // large, and ep_plastic_tolerance quite small, the above NR process will quickly converge, but
1132 // the solution may be wrong and violate Kuhn-Tucker
1133 return false;
1134
1135 if (yf[0] > _f_tol || yf[2] > _f_tol || yf[4] > _f_tol || yf[5] > _f_tol)
1136 // Kuhn-Tucker failure
1137 return false;
1138
1139 // success
1140 dpm[0] = dpm[2] = dpm[4] = dpm[5] = 0;
1141 return true;
1142}
1143
1144bool
1145SolidMechanicsPlasticMohrCoulombMulti::KuhnTuckerOK(const std::vector<Real> & yf,
1146 const std::vector<Real> & dpm,
1147 Real ep_plastic_tolerance) const
1148{
1149 mooseAssert(
1150 yf.size() == 6,
1151 "SolidMechanicsPlasticMohrCoulomb::KuhnTuckerOK requires yf to be size 6, but your is "
1152 << yf.size());
1153 mooseAssert(
1154 dpm.size() == 6,
1155 "SolidMechanicsPlasticMohrCoulomb::KuhnTuckerOK requires dpm to be size 6, but your is "
1156 << dpm.size());
1157 for (unsigned i = 0; i < 6; ++i)
1158 if (!SolidMechanicsPlasticModel::KuhnTuckerSingleSurface(yf[i], dpm[i], ep_plastic_tolerance))
1159 return false;
1160 return true;
1161}
1162
1163bool
f
Real f(Real x)
Test function for Brents method.
Definition BrentsMethodTest.C:21
tol
const double tol
Definition EquilibriumConstantFitTest.C:16
x
const std::vector< double > x
Definition EquilibriumConstantFitTest.C:17
v
const double v
Definition GasLiquidMassTransferTest.C:23
registerMooseObject
registerMooseObject("SolidMechanicsApp", SolidMechanicsPlasticMohrCoulombMulti)
registerMooseObjectRenamed
registerMooseObjectRenamed("SolidMechanicsApp", TensorMechanicsPlasticMohrCoulombMulti, "01/01/2025 00:00", SolidMechanicsPlasticMohrCoulombMulti)
int
void ErrorVector unsigned int
GeneralUserObject::_console
const ConsoleStream _console
InputParameters::addRequiredParam
void addRequiredParam(const std::string &name, const std::string &doc_string)
InputParameters::addParam
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
InputParameters::addClassDescription
void addClassDescription(const std::string &doc_string)
GeneralUserObject::mooseError
void mooseError(Args &&... args) const
MooseRandom::rand
static Real rand()
MooseRandom::seed
void seed(std::size_t i, unsigned int seed)
RankTwoTensorTempl::dsymmetricEigenvalues
void dsymmetricEigenvalues(std::vector< T > &eigvals, std::vector< RankTwoTensorTempl< T > > &deigvals) const
RankTwoTensorTempl::transpose
RankTwoTensorTempl< T > transpose() const
RankTwoTensorTempl::print
void print(std::ostream &stm=Moose::out) const
RankTwoTensorTempl::d2symmetricEigenvalues
void d2symmetricEigenvalues(std::vector< RankFourTensorTempl< T > > &deriv) const
RankTwoTensorTempl::symmetricEigenvalues
void symmetricEigenvalues(std::vector< T > &eigvals) const
RankTwoTensorTempl::symmetricEigenvaluesEigenvectors
void symmetricEigenvaluesEigenvectors(std::vector< T > &eigvals, RankTwoTensorTempl< T > &eigvecs) const
SolidMechanicsHardeningModel::derivative
virtual Real derivative(Real intnl) const
Definition SolidMechanicsHardeningModel.C:50
SolidMechanicsHardeningModel::modelName
virtual std::string modelName() const =0
SolidMechanicsHardeningModel::value
virtual Real value(Real intnl) const
Definition SolidMechanicsHardeningModel.C:44
SolidMechanicsPlasticModel
Plastic Model base class The virtual functions written below must be over-ridden in derived classes t...
Definition SolidMechanicsPlasticModel.h:38
SolidMechanicsPlasticModel::returnMap
virtual bool returnMap(const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const
Performs a custom return-map.
Definition SolidMechanicsPlasticModel.C:219
SolidMechanicsPlasticModel::KuhnTuckerSingleSurface
bool KuhnTuckerSingleSurface(Real yf, Real dpm, Real dpm_tol) const
Returns true if the Kuhn-Tucker conditions for the single surface are satisfied.
Definition SolidMechanicsPlasticModel.C:246
SolidMechanicsPlasticModel::_f_tol
const Real _f_tol
Tolerance on yield function.
Definition SolidMechanicsPlasticModel.h:170
SolidMechanicsPlasticModel::validParams
static InputParameters validParams()
Definition SolidMechanicsPlasticModel.C:16
SolidMechanicsPlasticMohrCoulombMulti
FiniteStrainMohrCoulombMulti implements rate-independent non-associative mohr-coulomb with hardening/...
Definition SolidMechanicsPlasticMohrCoulombMulti.h:20
SolidMechanicsPlasticMohrCoulombMulti::dyieldFunction_dstressV
virtual void dyieldFunction_dstressV(const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &df_dstress) const override
The derivative of yield functions with respect to stress.
Definition SolidMechanicsPlasticMohrCoulombMulti.C:165
SolidMechanicsPlasticMohrCoulombMulti::_shift
const Real _shift
yield function is shifted by this amount to avoid problems with stress-derivatives at equal eigenvalu...
Definition SolidMechanicsPlasticMohrCoulombMulti.h:107
SolidMechanicsPlasticMohrCoulombMulti::perturbStress
void perturbStress(const RankTwoTensor &stress, std::vector< Real > &eigvals, std::vector< RankTwoTensor > &deigvals) const
perturbs the stress tensor in the case of almost-equal eigenvalues.
Definition SolidMechanicsPlasticMohrCoulombMulti.C:124
SolidMechanicsPlasticMohrCoulombMulti::SolidMechanicsPlasticMohrCoulombMulti
SolidMechanicsPlasticMohrCoulombMulti(const InputParameters &parameters)
Definition SolidMechanicsPlasticMohrCoulombMulti.C:61
SolidMechanicsPlasticMohrCoulombMulti::returnEdge010100
bool returnEdge010100(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, Real mag_E, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
Tries to return-map to the MC edge using the n[1] and n[3] directions The return value is true if the...
Definition SolidMechanicsPlasticMohrCoulombMulti.C:1028
SolidMechanicsPlasticMohrCoulombMulti::df_dsig
void df_dsig(const RankTwoTensor &stress, Real sin_angle, std::vector< RankTwoTensor > &df) const
this is exactly dyieldFunction_dstress, or flowPotential, depending on whether sin_angle = sin(phi),...
Definition SolidMechanicsPlasticMohrCoulombMulti.C:144
SolidMechanicsPlasticMohrCoulombMulti::dpsi
virtual Real dpsi(const Real internal_param) const
d(psi)/d(internal_param) as a function of residual value, rate, and internal_param
Definition SolidMechanicsPlasticMohrCoulombMulti.C:315
SolidMechanicsPlasticMohrCoulombMulti::_max_iters
const unsigned int _max_iters
Maximum Newton-Raphison iterations in the custom returnMap algorithm.
Definition SolidMechanicsPlasticMohrCoulombMulti.h:104
SolidMechanicsPlasticMohrCoulombMulti::dcohesion
virtual Real dcohesion(const Real internal_param) const
d(cohesion)/d(internal_param) as a function of residual value, rate, and internal_param
Definition SolidMechanicsPlasticMohrCoulombMulti.C:291
SolidMechanicsPlasticMohrCoulombMulti::dflowPotential_dstressV
virtual void dflowPotential_dstressV(const RankTwoTensor &stress, Real intnl, std::vector< RankFourTensor > &dr_dstress) const override
The derivative of the flow potential with respect to stress.
Definition SolidMechanicsPlasticMohrCoulombMulti.C:204
SolidMechanicsPlasticMohrCoulombMulti::numberSurfaces
virtual unsigned int numberSurfaces() const override
The number of yield surfaces for this plasticity model.
Definition SolidMechanicsPlasticMohrCoulombMulti.C:83
SolidMechanicsPlasticMohrCoulombMulti::cohesion
virtual Real cohesion(const Real internal_param) const
cohesion as a function of residual value, rate, and internal_param
Definition SolidMechanicsPlasticMohrCoulombMulti.C:285
SolidMechanicsPlasticMohrCoulombMulti::_cohesion
const SolidMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
Definition SolidMechanicsPlasticMohrCoulombMulti.h:95
SolidMechanicsPlasticMohrCoulombMulti::_phi
const SolidMechanicsHardeningModel & _phi
Hardening model for phi.
Definition SolidMechanicsPlasticMohrCoulombMulti.h:98
SolidMechanicsPlasticMohrCoulombMulti::doReturnMap
bool doReturnMap(const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const
See doco for returnMap function.
Definition SolidMechanicsPlasticMohrCoulombMulti.C:363
SolidMechanicsPlasticMohrCoulombMulti::KuhnTuckerOK
bool KuhnTuckerOK(const std::vector< Real > &yf, const std::vector< Real > &dpm, Real ep_plastic_tolerance) const
Returns true if the Kuhn-Tucker conditions are satisfied.
Definition SolidMechanicsPlasticMohrCoulombMulti.C:1145
SolidMechanicsPlasticMohrCoulombMulti::yieldFunctionV
virtual void yieldFunctionV(const RankTwoTensor &stress, Real intnl, std::vector< Real > &f) const override
Calculates the yield functions.
Definition SolidMechanicsPlasticMohrCoulombMulti.C:89
SolidMechanicsPlasticMohrCoulombMulti::useCustomReturnMap
virtual bool useCustomReturnMap() const override
Returns false. You will want to override this in your derived class if you write a custom returnMap f...
Definition SolidMechanicsPlasticMohrCoulombMulti.C:1164
SolidMechanicsPlasticMohrCoulombMulti::psi
virtual Real psi(const Real internal_param) const
psi as a function of residual value, rate, and internal_param
Definition SolidMechanicsPlasticMohrCoulombMulti.C:309
SolidMechanicsPlasticMohrCoulombMulti::returnTip
bool returnTip(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
Tries to return-map to the MC tip using the THREE directions given in n, and THREE dpm values are ret...
Definition SolidMechanicsPlasticMohrCoulombMulti.C:638
SolidMechanicsPlasticMohrCoulombMulti::returnPlane
bool returnPlane(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
Tries to return-map to the MC plane using the n[3] direction The return value is true if the internal...
Definition SolidMechanicsPlasticMohrCoulombMulti.C:806
SolidMechanicsPlasticMohrCoulombMulti::_use_custom_returnMap
const bool _use_custom_returnMap
Whether to use the custom return-map algorithm.
Definition SolidMechanicsPlasticMohrCoulombMulti.h:110
SolidMechanicsPlasticMohrCoulombMulti::returnEdge000101
bool returnEdge000101(const std::vector< Real > &eigvals, const std::vector< RealVectorValue > &n, std::vector< Real > &dpm, RankTwoTensor &returned_stress, Real intnl_old, Real &sinphi, Real &cohcos, Real initial_guess, Real mag_E, bool &nr_converged, Real ep_plastic_tolerance, std::vector< Real > &yf) const
Tries to return-map to the MC edge using the n[4] and n[6] directions The return value is true if the...
Definition SolidMechanicsPlasticMohrCoulombMulti.C:911
SolidMechanicsPlasticMohrCoulombMulti::returnMap
virtual bool returnMap(const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const override
Performs a custom return-map.
Definition SolidMechanicsPlasticMohrCoulombMulti.C:327
SolidMechanicsPlasticMohrCoulombMulti::_psi
const SolidMechanicsHardeningModel & _psi
Hardening model for psi.
Definition SolidMechanicsPlasticMohrCoulombMulti.h:101
SolidMechanicsPlasticMohrCoulombMulti::phi
virtual Real phi(const Real internal_param) const
phi as a function of residual value, rate, and internal_param
Definition SolidMechanicsPlasticMohrCoulombMulti.C:297
SolidMechanicsPlasticMohrCoulombMulti::dflowPotential_dintnlV
virtual void dflowPotential_dintnlV(const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &dr_dintnl) const override
The derivative of the flow potential with respect to the internal parameter.
Definition SolidMechanicsPlasticMohrCoulombMulti.C:228
SolidMechanicsPlasticMohrCoulombMulti::yieldFunctionEigvals
void yieldFunctionEigvals(Real e0, Real e1, Real e2, Real sinphi, Real cohcos, std::vector< Real > &f) const
Calculates the yield functions given the eigenvalues of stress.
Definition SolidMechanicsPlasticMohrCoulombMulti.C:106
SolidMechanicsPlasticMohrCoulombMulti::dyieldFunction_dintnlV
virtual void dyieldFunction_dintnlV(const RankTwoTensor &stress, Real intnl, std::vector< Real > &df_dintnl) const override
The derivative of yield functions with respect to the internal parameter.
Definition SolidMechanicsPlasticMohrCoulombMulti.C:173
SolidMechanicsPlasticMohrCoulombMulti::flowPotentialV
virtual void flowPotentialV(const RankTwoTensor &stress, Real intnl, std::vector< RankTwoTensor > &r) const override
The flow potentials.
Definition SolidMechanicsPlasticMohrCoulombMulti.C:195
SolidMechanicsPlasticMohrCoulombMulti::dphi
virtual Real dphi(const Real internal_param) const
d(phi)/d(internal_param) as a function of residual value, rate, and internal_param
Definition SolidMechanicsPlasticMohrCoulombMulti.C:303
SolidMechanicsPlasticMohrCoulombMulti::activeConstraints
virtual void activeConstraints(const std::vector< Real > &f, const RankTwoTensor &stress, Real intnl, const RankFourTensor &Eijkl, std::vector< bool > &act, RankTwoTensor &returned_stress) const override
The active yield surfaces, given a vector of yield functions.
Definition SolidMechanicsPlasticMohrCoulombMulti.C:248
Generated on Mon Aug 24 2026 19:10:41 for https://mooseframework.inl.gov by doxygen 1.9.8