11 #include "RankFourTensor.h"
12 #include "libmesh/utility.h"
22 params.addParam<
bool>(
"use_custom_returnMap",
24 "Whether to use the custom returnMap "
25 "algorithm. Set to true if you are using "
26 "isotropic elasticity.");
27 params.addParam<
bool>(
"use_custom_cto",
29 "Whether to use the custom consistent tangent "
30 "operator computations. Set to true if you are "
31 "using isotropic elasticity.");
32 params.addClassDescription(
"J2 plasticity, associative, with hardening");
33 params.addRangeCheckedParam<Real>(
"smoother",
36 "The cone vertex at J2=0 is smoothed. The maximum mean "
37 "stress possible, which is Cohesion*Cot(friction_angle) for "
38 "smoother=0, becomes (Cohesion - "
39 "smoother/3)*Cot(friction_angle). This is a non-hardening "
41 params.addRangeCheckedParam<
unsigned>(
45 "Maximum iterations to use in the custom return map function");
46 params.addClassDescription(
47 "Non-associative Drucker Prager plasticity with hyperbolic smoothing of the cone tip.");
52 const InputParameters & parameters)
54 _smoother2(Utility::
pow<2>(getParam<Real>(
"smoother"))),
55 _use_custom_returnMap(getParam<bool>(
"use_custom_returnMap")),
56 _use_custom_cto(getParam<bool>(
"use_custom_cto")),
57 _max_iters(getParam<unsigned>(
"max_iterations"))
68 return std::sqrt(stress.secondInvariant() +
_smoother2) + stress.trace() * bbb - aaa;
74 return 0.5 * stress.dsecondInvariant() / std::sqrt(stress.secondInvariant() +
_smoother2) +
75 stress.dtrace() * bbb;
83 dr_dstress = 0.5 * stress.d2secondInvariant() / std::sqrt(stress.secondInvariant() +
_smoother2);
84 dr_dstress += -0.5 * 0.5 * stress.dsecondInvariant().outerProduct(stress.dsecondInvariant()) /
92 return "DruckerPragerHyperbolic";
99 Real ep_plastic_tolerance,
101 Real & returned_intnl,
102 std::vector<Real> & dpm,
104 std::vector<Real> & yf,
105 bool & trial_stress_inadmissible)
const
111 ep_plastic_tolerance,
117 trial_stress_inadmissible);
126 trial_stress_inadmissible =
false;
130 trial_stress_inadmissible =
true;
131 const Real mu = E_ijkl(0, 1, 0, 1);
132 const Real lambda = E_ijkl(0, 0, 0, 0) - 2.0 * mu;
133 const Real bulky = 3.0 * lambda + 2.0 * mu;
134 const Real Tr_trial = trial_stress.trace();
135 const Real J2trial = trial_stress.secondInvariant();
151 unsigned int iter = 0;
154 bothAB(intnl_old + dpm[0], aaa, bbb);
155 dbothAB(intnl_old + dpm[0], daaa, dbbb);
158 ll = aaa - bbb * (Tr_trial - dpm[0] * bulky * 3 * bbb_flow);
159 dll = daaa - dbbb * (Tr_trial - dpm[0] * bulky * 3 * bbb_flow) +
160 bbb * bulky * 3 * (bbb_flow + dpm[0] * dbbb_flow);
161 residual = bbb * (Tr_trial - dpm[0] * bulky * 3 * bbb_flow) - aaa +
162 std::sqrt(J2trial / Utility::pow<2>(1 + dpm[0] * mu / ll) +
_smoother2);
163 jac = dbbb * (Tr_trial - dpm[0] * bulky * 3 * bbb_flow) -
164 bbb * bulky * 3 * (bbb_flow + dpm[0] * dbbb_flow) - daaa +
165 0.5 * J2trial * (-2.0) * (mu / ll - dpm[0] * mu * dll / ll / ll) /
166 Utility::pow<3>(1 + dpm[0] * mu / ll) /
167 std::sqrt(J2trial / Utility::pow<2>(1.0 + dpm[0] * mu / ll) +
_smoother2);
168 dpm[0] += -residual / jac;
176 returned_intnl = intnl_old + dpm[0];
178 bothAB(returned_intnl, aaa, bbb);
180 ll = aaa - bbb * (Tr_trial - dpm[0] * bulky * 3.0 * bbb_flow);
182 trial_stress.deviatoric() / (1.0 + dpm[0] * mu / ll);
188 const Real returned_trace_over_3 = (aaa - ll) / bbb / 3.0;
189 for (
unsigned i = 0; i < 3; ++i)
191 returned_stress(i, i) += returned_trace_over_3;
192 rij(i, i) += bbb_flow;
195 delta_dp = rij * dpm[0];
207 const std::vector<Real> & cumulative_pm)
const
211 trial_stress, intnl_old, stress, intnl, E_ijkl, cumulative_pm);
213 const Real mu = E_ijkl(0, 1, 0, 1);
214 const Real la = E_ijkl(0, 0, 0, 0) - 2.0 * mu;
215 const Real bulky = 3.0 * la + 2.0 * mu;
222 const Real bbb_flow_mod = bbb_flow + cumulative_pm[0] * dbbb_flow;
223 const Real J2 = stress.secondInvariant();
227 const Real one_over_h =
229 3.0 * bbb * bbb_flow_mod * bulky);
232 mu * sij / sq + bulky * bbb *
RankTwoTensor(RankTwoTensor::initIdentity);
235 bulky * bbb_flow_mod *
RankTwoTensor(RankTwoTensor::initIdentity);
238 E_ijkl - one_over_h * rmod_timesE.outerProduct(df_dsig_timesE);
240 const RankFourTensor dr_dsig_timesE = -0.5 * mu * sij.outerProduct(sij) / Utility::pow<3>(sq) +
241 mu * stress.d2secondInvariant() / sq;
243 0.5 * mu *
_smoother2 * sij / Utility::pow<4>(sq);
248 (dr_dsig_timesE - one_over_h * rmod_timesE.outerProduct(df_dsig_E_dr_dsig));
253 s_inv = coeff_si.invSymm();
255 catch (MooseException & e)
260 return s_inv * coeff_ep;