11 #include "libmesh/utility.h" 22 "A SolidMechanicsHardening UserObject that defines hardening of the " 23 "tensile strength. In physical situations this is positive (and always " 24 "must be greater than negative compressive-strength.");
26 "compressive_strength",
27 "A SolidMechanicsHardening UserObject that defines hardening of the " 28 "compressive strength. In physical situations this is positive.");
30 "cohesion",
"A SolidMechanicsHardening UserObject that defines hardening of the cohesion");
32 "A SolidMechanicsHardening UserObject " 33 "that defines hardening of the " 34 "friction angle (in radians)");
37 "A SolidMechanicsHardening UserObject that defines hardening of the " 38 "dilation angle (in radians). Unless you are quite confident, this should " 39 "be set positive and not greater than the friction angle.");
40 params.
addParam<
bool>(
"perfect_guess",
42 "Provide a guess to the Newton-Raphson procedure " 43 "that is the result from perfect plasticity. With " 44 "severe hardening/softening this may be " 47 "tensile (Rankine) and compressive caps, with hardening/softening");
58 _perfect_guess(getParam<bool>(
"perfect_guess")),
64 mooseWarning(
"Usually the Mohr-Coulomb dilation angle is positive and not greater than the " 70 std::vector<Real> & stress_params)
const 76 std::vector<RankTwoTensor>
80 std::vector<RankTwoTensor> dsp;
85 std::vector<RankFourTensor>
88 std::vector<RankFourTensor> d2;
96 const std::vector<Real> & ,
97 const std::vector<Real> & ,
100 std::vector<Real> eigvals;
107 const std::vector<Real> & stress_params,
109 const std::vector<Real> & ,
115 stress =
RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0.0, 0.0, 0.0);
122 const std::vector<Real> & intnl,
123 std::vector<Real> & yf)
const 130 const Real smax = stress_params[2];
131 const Real smid = stress_params[1];
132 const Real smin = stress_params[0];
139 yf[6] = 0.5 * (smax - smin) + 0.5 * (smax + smin) * sinphi - cohcos;
140 yf[7] = 0.5 * (smid - smin) + 0.5 * (smid + smin) * sinphi - cohcos;
141 yf[8] = 0.5 * (smax - smid) + 0.5 * (smax + smid) * sinphi - cohcos;
142 yf[9] = 0.5 * (smid - smax) + 0.5 * (smid + smax) * sinphi - cohcos;
143 yf[10] = 0.5 * (smin - smid) + 0.5 * (smin + smid) * sinphi - cohcos;
144 yf[11] = 0.5 * (smin - smax) + 0.5 * (smin + smax) * sinphi - cohcos;
149 const std::vector<Real> & intnl,
150 std::vector<yieldAndFlow> & all_q)
const 164 const Real smax = stress_params[2];
165 const Real smid = stress_params[1];
166 const Real smin = stress_params[0];
169 all_q[0].f = smax - ts;
170 all_q[1].f = smid - ts;
171 all_q[2].f = smin - ts;
172 all_q[3].f = -smin - cs;
173 all_q[4].f = -smid - cs;
174 all_q[5].f = -smax - cs;
175 all_q[6].f = 0.5 * (smax - smin) + 0.5 * (smax + smin) * sinphi - cohcos;
176 all_q[7].f = 0.5 * (smid - smin) + 0.5 * (smid + smin) * sinphi - cohcos;
177 all_q[8].f = 0.5 * (smax - smid) + 0.5 * (smax + smid) * sinphi - cohcos;
178 all_q[9].f = 0.5 * (smid - smax) + 0.5 * (smid + smax) * sinphi - cohcos;
179 all_q[10].f = 0.5 * (smin - smid) + 0.5 * (smin + smid) * sinphi - cohcos;
180 all_q[11].f = 0.5 * (smin - smax) + 0.5 * (smin + smax) * sinphi - cohcos;
183 for (
unsigned yf = 0; yf <
_num_yf; ++yf)
185 all_q[yf].df[
a] = 0.0;
186 all_q[0].df[2] = 1.0;
187 all_q[1].df[1] = 1.0;
188 all_q[2].df[0] = 1.0;
189 all_q[3].df[0] = -1.0;
190 all_q[4].df[1] = -1.0;
191 all_q[5].df[2] = -1.0;
192 all_q[6].df[2] = 0.5 * (1.0 + sinphi);
193 all_q[6].df[0] = 0.5 * (-1.0 + sinphi);
194 all_q[7].df[1] = 0.5 * (1.0 + sinphi);
195 all_q[7].df[0] = 0.5 * (-1.0 + sinphi);
196 all_q[8].df[2] = 0.5 * (1.0 + sinphi);
197 all_q[8].df[1] = 0.5 * (-1.0 + sinphi);
198 all_q[9].df[1] = 0.5 * (1.0 + sinphi);
199 all_q[9].df[2] = 0.5 * (-1.0 + sinphi);
200 all_q[10].df[0] = 0.5 * (1.0 + sinphi);
201 all_q[10].df[1] = 0.5 * (-1.0 + sinphi);
202 all_q[11].df[0] = 0.5 * (1.0 + sinphi);
203 all_q[11].df[2] = 0.5 * (-1.0 + sinphi);
206 for (
unsigned i = 0; i < 6; ++i)
207 all_q[i].df_di[0] = 0.0;
208 all_q[0].df_di[1] = all_q[1].df_di[1] = all_q[2].df_di[1] = -dts;
209 all_q[3].df_di[1] = all_q[4].df_di[1] = all_q[5].df_di[1] = -dcs;
210 for (
unsigned i = 6; i < 12; ++i)
211 all_q[i].df_di[1] = 0.0;
212 all_q[6].df_di[0] = 0.5 * (smax + smin) * dsinphi - dcohcos;
213 all_q[7].df_di[0] = 0.5 * (smid + smin) * dsinphi - dcohcos;
214 all_q[8].df_di[0] = 0.5 * (smax + smid) * dsinphi - dcohcos;
215 all_q[9].df_di[0] = 0.5 * (smid + smax) * dsinphi - dcohcos;
216 all_q[10].df_di[0] = 0.5 * (smin + smid) * dsinphi - dcohcos;
217 all_q[11].df_di[0] = 0.5 * (smin + smax) * dsinphi - dcohcos;
221 for (
unsigned yf = 0; yf < 6; ++yf)
223 all_q[yf].dg[
a] = all_q[yf].df[
a];
224 all_q[6].dg[2] = all_q[7].dg[1] = all_q[8].dg[2] = all_q[9].dg[1] = all_q[10].dg[0] =
225 all_q[11].dg[0] = 0.5 * (1.0 + sinpsi);
226 all_q[6].dg[0] = all_q[7].dg[0] = all_q[8].dg[1] = all_q[9].dg[2] = all_q[10].dg[1] =
227 all_q[11].dg[2] = 0.5 * (-1.0 + sinpsi);
230 for (
unsigned yf = 0; yf <
_num_yf; ++yf)
233 all_q[yf].d2g_di[
a][i] = 0.0;
234 all_q[6].d2g_di[2][0] = all_q[7].d2g_di[1][0] = all_q[8].d2g_di[2][0] = all_q[9].d2g_di[1][0] =
235 all_q[10].d2g_di[0][0] = all_q[11].d2g_di[0][0] = 0.5 * dsinpsi;
236 all_q[6].d2g_di[0][0] = all_q[7].d2g_di[0][0] = all_q[8].d2g_di[1][0] = all_q[9].d2g_di[2][0] =
237 all_q[10].d2g_di[1][0] = all_q[11].d2g_di[2][0] = 0.5 * dsinpsi;
240 for (
unsigned yf = 0; yf <
_num_yf; ++yf)
243 all_q[yf].d2g[
a][
b] = 0.0;
259 for (
unsigned b = 0;
b <
a; ++
b)
266 const std::vector<Real> & intnl_old,
267 std::vector<Real> & stress_params,
269 std::vector<Real> & intnl)
const 273 for (
unsigned i = 0; i <
_num_sp; ++i)
274 stress_params[i] = trial_stress_params[i];
284 const Real trial_tensile_yf = trial_stress_params[2] - ts;
285 const Real trial_compressive_yf = -trial_stress_params[0] - cs;
286 const Real trial_mc_yf = 0.5 * (trial_stress_params[2] - trial_stress_params[0]) +
287 0.5 * (trial_stress_params[2] + trial_stress_params[0]) * sinphi -
315 bool found_solution =
false;
317 if (trial_tensile_yf <=
_f_tol && trial_compressive_yf <=
_f_tol && trial_mc_yf <=
_f_tol)
321 for (
unsigned i = 0; i <
_num_sp; ++i)
322 stress_params[i] = trial_stress_params[i];
324 found_solution =
true;
327 const bool tensile_possible = (ts < cohcos / sinphi);
328 const bool mc_tip_possible = (cohcos / sinphi < ts);
329 const bool mc_impossible = (0.5 * (ts + cs) + 0.5 * (ts - cs) * sinphi - cohcos <
333 const Real halfplus = 0.5 + 0.5 * sinpsi;
334 const Real neghalfplus = -0.5 + 0.5 * sinpsi;
336 if (!found_solution && tensile_possible && trial_tensile_yf >
_f_tol &&
337 (trial_compressive_yf <=
_f_tol || (trial_compressive_yf >
_f_tol && mc_impossible)))
343 const Real ga = (trial_stress_params[2] - ts) /
_Eij[2][2];
344 stress_params[2] = ts;
345 stress_params[1] = trial_stress_params[1] - ga *
_Eij[1][2];
346 stress_params[0] = trial_stress_params[0] - ga *
_Eij[0][2];
350 const Real to_subtract1 = stress_params[1] - (ts - 0.5 *
_shifter);
351 if (to_subtract1 > 0.0)
353 dist_mod += Utility::pow<2>(to_subtract1 / (trial_stress_params[2] - ts));
354 stress_params[1] -= to_subtract1;
356 const Real to_subtract0 = stress_params[0] - (ts -
_shifter);
357 if (to_subtract0 > 0.0)
359 dist_mod += Utility::pow<2>(to_subtract0 / (trial_stress_params[2] - ts));
360 stress_params[0] -= to_subtract0;
364 const Real to_add0 = -stress_params[0] - cs;
367 dist_mod += Utility::pow<2>(to_add0 / (trial_stress_params[2] - ts));
368 stress_params[0] += to_add0;
370 const Real to_add1 = -cs + 0.5 *
_shifter - stress_params[1];
373 dist_mod += Utility::pow<2>(to_add1 / (trial_stress_params[2] - ts));
374 stress_params[1] += to_add1;
378 const Real new_compressive_yf = -stress_params[0] - cs;
379 const Real new_mc_yf = 0.5 * (stress_params[2] - stress_params[0]) +
380 0.5 * (stress_params[2] + stress_params[0]) * sinphi - cohcos;
381 if (new_mc_yf <=
_f_tol && new_compressive_yf <=
_f_tol)
383 gaE =
std::sqrt(dist_mod) * (trial_stress_params[2] - stress_params[2]);
384 found_solution =
true;
387 if (!found_solution && trial_compressive_yf >
_f_tol &&
388 (trial_tensile_yf <=
_f_tol || (trial_tensile_yf >
_f_tol && mc_impossible)))
393 const Real ga = (trial_stress_params[0] + cs) /
_Eij[0][0];
394 stress_params[0] = -cs;
395 stress_params[1] = trial_stress_params[1] - ga *
_Eij[1][0];
396 stress_params[2] = trial_stress_params[2] - ga *
_Eij[2][0];
400 const Real to_add1 = -cs + 0.5 *
_shifter - stress_params[1];
403 dist_mod += Utility::pow<2>(to_add1 / (trial_stress_params[0] + cs));
404 stress_params[1] += to_add1;
409 dist_mod += Utility::pow<2>(to_add2 / (trial_stress_params[0] + cs));
410 stress_params[2] += to_add2;
414 const Real to_subtract2 = stress_params[2] - ts;
415 if (to_subtract2 > 0.0)
417 dist_mod += Utility::pow<2>(to_subtract2 / (trial_stress_params[0] + cs));
418 stress_params[2] -= to_subtract2;
420 const Real to_subtract1 = stress_params[1] - (ts - 0.5 *
_shifter);
421 if (to_subtract1 > 0.0)
423 dist_mod += Utility::pow<2>(to_subtract1 / (trial_stress_params[0] + cs));
424 stress_params[1] -= to_subtract1;
428 const Real new_tensile_yf = stress_params[2] - ts;
429 const Real new_mc_yf = 0.5 * (stress_params[2] - stress_params[0]) +
430 0.5 * (stress_params[2] + stress_params[0]) * sinphi - cohcos;
433 gaE =
std::sqrt(dist_mod) * (-trial_stress_params[0] + stress_params[0]);
434 found_solution =
true;
437 if (!found_solution && !mc_impossible && trial_mc_yf >
_f_tol)
441 const Real ga = ((trial_stress_params[2] - trial_stress_params[0]) +
442 (trial_stress_params[2] + trial_stress_params[0]) * sinphi - 2.0 * cohcos) /
445 trial_stress_params[2] - ga * (
_Eij[2][2] * halfplus +
_Eij[2][0] * neghalfplus);
446 stress_params[1] = trial_stress_params[1] - ga *
_Eij[1][0] * sinpsi;
448 trial_stress_params[0] - ga * (
_Eij[0][0] * neghalfplus +
_Eij[0][2] * halfplus);
449 const Real f7 = 0.5 * (stress_params[1] - stress_params[0]) +
450 0.5 * (stress_params[1] + stress_params[0]) * sinphi - cohcos;
451 const Real f8 = 0.5 * (stress_params[2] - stress_params[1]) +
452 0.5 * (stress_params[2] + stress_params[1]) * sinphi - cohcos;
453 const Real new_tensile_yf = stress_params[2] - ts;
454 const Real new_compressive_yf = -stress_params[0] - cs;
458 gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
459 (stress_params[2] - stress_params[0])) /
461 found_solution =
true;
464 if (!found_solution && !mc_impossible && trial_mc_yf >
_f_tol)
479 const Real cmax6 =
_Eij[2][2] * halfplus +
_Eij[2][0] * neghalfplus;
480 const Real cmax7 =
_Eij[2][1] * halfplus +
_Eij[2][0] * neghalfplus;
483 const Real cmin6 =
_Eij[0][2] * halfplus +
_Eij[0][0] * neghalfplus;
484 const Real cmin7 =
_Eij[0][1] * halfplus +
_Eij[0][0] * neghalfplus;
488 const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
489 const Real c7 = 0.5 * (cmax7 - cmin7) + 0.5 * (cmax7 + cmin7) * sinphi;
495 const Real f6_trial = trial_mc_yf;
496 const Real f7_trial = 0.5 * (trial_stress_params[1] - trial_stress_params[0]) +
497 0.5 * (trial_stress_params[1] + trial_stress_params[0]) * sinphi -
499 const Real descr = Utility::pow<2>(c6) - Utility::pow<2>(c7);
500 Real ga6 = (c6 * f6_trial - c7 * f7_trial) / descr;
501 Real ga7 = (-c7 * f6_trial + c6 * f7_trial) / descr;
503 stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga7 * cmax7;
504 stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga7 * cmin7;
505 stress_params[1] = stress_params[2] - 0.5 *
_shifter;
507 Real f8 = 0.5 * (stress_params[2] - stress_params[1]) +
508 0.5 * (stress_params[2] + stress_params[1]) * sinphi - cohcos;
510 if (mc_tip_possible && f8 >
_f_tol)
512 stress_params[2] = cohcos / sinphi;
513 stress_params[1] = stress_params[2] - 0.5 *
_shifter;
514 stress_params[0] = stress_params[2] -
_shifter;
520 const Real new_tensile_yf = stress_params[2] - ts;
521 const Real new_compressive_yf = -stress_params[0] - cs;
524 ga6 >= 0.0 && ga7 >= 0.0)
526 gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
527 (stress_params[2] - stress_params[0])) /
529 found_solution =
true;
533 if (!found_solution && !mc_impossible && trial_mc_yf >
_f_tol)
548 const Real cmax6 =
_Eij[2][2] * halfplus +
_Eij[2][0] * neghalfplus;
549 const Real cmax8 =
_Eij[2][2] * halfplus +
_Eij[2][1] * neghalfplus;
552 const Real cmin6 =
_Eij[0][2] * halfplus +
_Eij[0][0] * neghalfplus;
553 const Real cmin8 =
_Eij[0][2] * halfplus +
_Eij[0][1] * neghalfplus;
557 const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
558 const Real c8 = 0.5 * (cmax8 - cmin8) + 0.5 * (cmax8 + cmin8) * sinphi;
564 const Real f6_trial = trial_mc_yf;
565 const Real f8_trial = 0.5 * (trial_stress_params[2] - trial_stress_params[1]) +
566 0.5 * (trial_stress_params[2] + trial_stress_params[1]) * sinphi -
568 const Real descr = Utility::pow<2>(c6) - Utility::pow<2>(c8);
569 Real ga6 = (c6 * f6_trial - c8 * f8_trial) / descr;
570 Real ga8 = (-c8 * f6_trial + c6 * f8_trial) / descr;
572 stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga8 * cmax8;
573 stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga8 * cmin8;
574 stress_params[1] = stress_params[0] + 0.5 *
_shifter;
576 Real f7 = 0.5 * (stress_params[1] - stress_params[0]) +
577 0.5 * (stress_params[1] + stress_params[0]) * sinphi - cohcos;
579 if (mc_tip_possible && f7 >
_f_tol)
581 stress_params[2] = cohcos / sinphi;
582 stress_params[1] = stress_params[2] - 0.5 *
_shifter;
583 stress_params[0] = stress_params[2] -
_shifter;
589 const Real new_tensile_yf = stress_params[2] - ts;
590 const Real new_compressive_yf = -stress_params[0] - cs;
593 ga6 >= 0.0 && ga8 >= 0.0)
595 gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
596 (stress_params[2] - stress_params[0])) /
598 found_solution =
true;
602 if (!found_solution && !mc_impossible && tensile_possible && trial_tensile_yf >
_f_tol)
617 const Real cmax6 =
_Eij[2][2] * halfplus +
_Eij[2][0] * neghalfplus;
619 const Real cmid6 =
_Eij[1][2] * halfplus +
_Eij[1][0] * neghalfplus;
621 const Real cmin6 =
_Eij[0][2] * halfplus +
_Eij[0][0] * neghalfplus;
626 const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
627 const Real c0 = 0.5 * (cmax0 - cmin0) + 0.5 * (cmax0 + cmin0) * sinphi;
631 const Real descr = c0 * cmax6 - c6 * cmax0;
634 const Real ga0 = (-c6 * trial_tensile_yf + cmax6 * trial_mc_yf) / descr;
635 const Real ga6 = (c0 * trial_tensile_yf - cmax0 * trial_mc_yf) / descr;
636 stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga0 * cmax0;
637 stress_params[1] = trial_stress_params[1] - ga6 * cmid6 - ga0 * cmid0;
638 stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga0 * cmin0;
642 std::min(stress_params[0],
647 const bool goto_corner = (stress_params[1] >= stress_params[2] - 0.5 *
_shifter);
648 stress_params[1] = std::max(std::min(stress_params[1], stress_params[2] - 0.5 *
_shifter),
651 const Real new_compressive_yf = -stress_params[0] - cs;
652 if (new_compressive_yf <=
_f_tol &&
653 (goto_corner || (ga0 >= 0.0 && ga6 >= 0.0)))
655 gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
656 (stress_params[2] - stress_params[0])) /
658 (trial_stress_params[2] - stress_params[2]);
659 found_solution =
true;
663 if (!found_solution && !mc_impossible)
678 const Real cmax6 =
_Eij[2][2] * halfplus +
_Eij[2][0] * neghalfplus;
680 const Real cmid6 =
_Eij[1][2] * halfplus +
_Eij[1][0] * neghalfplus;
682 const Real cmin6 =
_Eij[0][2] * halfplus +
_Eij[0][0] * neghalfplus;
687 const Real c6 = 0.5 * (cmax6 - cmin6) + 0.5 * (cmax6 + cmin6) * sinphi;
688 const Real c3 = 0.5 * (cmax3 - cmin3) + 0.5 * (cmax3 + cmin3) * sinphi;
692 const Real descr = c3 * cmin6 - c6 * cmin3;
695 const Real ga3 = (c6 * trial_compressive_yf + cmin6 * trial_mc_yf) / descr;
696 const Real ga6 = (-c3 * trial_compressive_yf - cmin3 * trial_mc_yf) / descr;
697 stress_params[2] = trial_stress_params[2] - ga6 * cmax6 - ga3 * cmax3;
698 stress_params[1] = trial_stress_params[1] - ga6 * cmid6 - ga3 * cmid3;
699 stress_params[0] = trial_stress_params[0] - ga6 * cmin6 - ga3 * cmin3;
701 const Real new_tensile_yf = stress_params[2] - ts;
703 std::max(stress_params[2],
705 stress_params[1] = std::max(std::min(stress_params[1], stress_params[2] - 0.5 *
_shifter),
708 if (new_tensile_yf <= _f_tol && ga6 >= 0.0)
710 gaE = ((trial_stress_params[2] - trial_stress_params[0]) -
711 (stress_params[2] - stress_params[0])) /
713 (-trial_stress_params[0] - stress_params[0]);
715 found_solution =
true;
722 for (
unsigned i = 0; i <
_num_sp; ++i)
723 stress_params[i] = trial_stress_params[i];
725 mooseWarning(
"CappedMohrCoulombStressUpdate cannot initialize from max = ",
738 const std::vector<Real> & current_stress_params,
739 const std::vector<Real> & intnl_old,
740 std::vector<Real> & intnl)
const 743 const Real smax = current_stress_params[2];
744 const Real smin = current_stress_params[0];
745 const Real trial_smax = trial_stress_params[2];
746 const Real trial_smin = trial_stress_params[0];
747 const Real ga_shear = ((trial_smax - trial_smin) - (smax - smin)) / (
_Eij[2][2] -
_Eij[0][2]);
748 intnl[0] = intnl_old[0] + ga_shear;
750 const Real prefactor = (
_Eij[2][2] +
_Eij[0][2]) * sinpsi;
751 const Real shear_correction = prefactor * ga_shear;
753 ((trial_smax + trial_smin) - (smax + smin) - shear_correction) /
755 intnl[1] = intnl_old[1] + ga_tensile;
760 const std::vector<Real> & current_stress_params,
761 const std::vector<Real> & intnl,
762 std::vector<std::vector<Real>> & dintnl)
const 765 const Real smax = current_stress_params[2];
766 const Real smin = current_stress_params[0];
767 const Real trial_smax = trial_stress_params[2];
768 const Real trial_smin = trial_stress_params[0];
769 const Real ga_shear = ((trial_smax - trial_smin) - (smax - smin)) / (
_Eij[2][2] -
_Eij[0][2]);
770 const std::vector<Real> dga_shear = {
773 for (std::size_t i = 0; i <
_num_sp; ++i)
774 dintnl[0][i] = dga_shear[i];
779 const Real prefactor = (
_Eij[2][2] +
_Eij[0][2]) * sinpsi;
780 const Real dprefactor_di0 = (
_Eij[2][2] +
_Eij[0][2]) * dsinpsi_di0;
782 std::vector<Real> dshear_correction(
_num_sp);
783 for (std::size_t i = 0; i <
_num_sp; ++i)
784 dshear_correction[i] = prefactor * dga_shear[i] + dprefactor_di0 * dintnl[0][i] * ga_shear;
789 for (std::size_t i = 0; i <
_num_sp; ++i)
798 const std::vector<Real> & trial_stress_params,
800 const std::vector<Real> & stress_params,
804 bool compute_full_tangent_operator,
805 const std::vector<std::vector<Real>> & dvar_dtrial,
809 if (!compute_full_tangent_operator)
825 if (trial_stress_params[
a] == trial_stress_params[
j])
827 drot_dstress(i,
j,
k, l) +=
830 (trial_stress_params[
j] - trial_stress_params[
a]);
841 dstress_dtrial(i,
j,
k, l) +=
842 drot_dstress(i,
a,
k, l) * stress_params[
a] * eT(
a,
j) +
843 _eigvecs(i,
a) * stress_params[
a] * drot_dstress(
j,
a,
k, l);
852 dstress_dtrial(i,
j,
k, l) +=
const SolidMechanicsHardeningModel & _phi
Hardening model for friction angle.
void symmetricEigenvalues(std::vector< Real > &eigvals) const
CappedMohrCoulombStressUpdate implements rate-independent nonassociative Mohr-Coulomb plus tensile pl...
void computeStressParams(const RankTwoTensor &stress, std::vector< Real > &stress_params) const override
Computes stress_params, given stress.
void computeAllQV(const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< yieldAndFlow > &all_q) const override
Completely fills all_q with correct values.
void dsymmetricEigenvalues(std::vector< Real > &eigvals, std::vector< RankTwoTensorTempl< Real >> &deigvals) const
const bool _perfect_guess
Whether to provide an estimate of the returned stress, based on perfect plasticity.
const unsigned _num_intnl
Number of internal parameters.
virtual void consistentTangentOperatorV(const RankTwoTensor &stress_trial, const std::vector< Real > &trial_stress_params, const RankTwoTensor &stress, const std::vector< Real > &stress_params, Real gaE, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, bool compute_full_tangent_operator, const std::vector< std::vector< Real >> &dvar_dtrial, RankFourTensor &cto) override
Calculates the consistent tangent operator.
Real _En
normalising factor
Struct designed to hold info about a single yield function and its derivatives, as well as the flow d...
const SolidMechanicsHardeningModel & _tensile_strength
Hardening model for tensile strength.
std::vector< std::vector< Real > > _Cij
_Cij[i, j] * _Eij[j, k] = 1 iff j == k
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
virtual Real value(Real intnl) const
void mooseWarning(Args &&... args) const
Real elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
void setEffectiveElasticity(const RankFourTensor &Eijkl) override
Sets _Eij and _En and _Cij.
std::vector< std::vector< Real > > _Eij
E[i, j] in the system of equations to be solved.
const Real _f_tol
The yield-function tolerance.
const unsigned _num_yf
Number of yield functions.
CappedMohrCoulombStressUpdate(const InputParameters ¶meters)
std::vector< RankTwoTensor > dstress_param_dstress(const RankTwoTensor &stress) const override
d(stress_param[i])/d(stress) at given stress
void setIntnlValuesV(const std::vector< Real > &trial_stress_params, const std::vector< Real > ¤t_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &intnl) const override
Sets the internal parameters based on the trial values of stress_params, their current values...
const Real _shifter
When equal-eigenvalues are predicted from the stress initialization routine, shift them by this amoun...
void d2symmetricEigenvalues(std::vector< RankFourTensorTempl< Real >> &deriv) const
virtual void preReturnMapV(const std::vector< Real > &trial_stress_params, 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...
void yieldFunctionValuesV(const std::vector< Real > &stress_params, const std::vector< Real > &intnl, std::vector< Real > &yf) const override
Computes the values of the yield functions, given stress_params and intnl parameters.
registerMooseObject("SolidMechanicsApp", CappedMohrCoulombStressUpdate)
void setIntnlDerivativesV(const std::vector< Real > &trial_stress_params, const std::vector< Real > ¤t_stress_params, 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 stress_params, their current values, and the current values of the internal parameters.
const SolidMechanicsHardeningModel & _compressive_strength
Hardening model for compressive strength.
virtual Real derivative(Real intnl) const
RankTwoTensorTempl< Real > transpose() const
Hardening Model base class.
const SolidMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static InputParameters validParams()
void initializeVarsV(const std::vector< Real > &trial_stress_params, const std::vector< Real > &intnl_old, std::vector< Real > &stress_params, Real &gaE, std::vector< Real > &intnl) const override
Sets (stress_params, intnl) at "good guesses" of the solution to the Return-Map algorithm.
void symmetricEigenvaluesEigenvectors(std::vector< Real > &eigvals, RankTwoTensorTempl< Real > &eigvecs) const
MultiParameterPlasticityStressUpdate performs the return-map algorithm and associated stress updates ...
const SolidMechanicsHardeningModel & _psi
Hardening model for dilation angle.
RankTwoTensor _eigvecs
Eigenvectors of the trial stress as a RankTwoTensor, in order to rotate the returned stress back to s...
static InputParameters validParams()
virtual void setStressAfterReturnV(const RankTwoTensor &stress_trial, const std::vector< Real > &stress_params, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, RankTwoTensor &stress) const override
Sets stress from the admissible parameters.
std::vector< RankFourTensor > d2stress_param_dstress(const RankTwoTensor &stress) const override
d2(stress_param[i])/d(stress)/d(stress) at given stress
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
static constexpr unsigned _tensor_dimensionality
Internal dimensionality of tensors (currently this is 3 throughout tensor_mechanics) ...
static const std::string k
Real _poissons_ratio
Poisson's ratio.
const unsigned _num_sp
Number of stress parameters.