12 #include "libmesh/utility.h" 23 "A SolidMechanicsHardening UserObject that defines hardening of the cohesion. " 24 "Physically the cohesion should not be negative.");
26 "A SolidMechanicsHardening UserObject that defines " 27 "hardening of tan(friction angle). Physically the " 28 "friction angle should be between 0 and 90deg.");
31 "A SolidMechanicsHardening UserObject that defines hardening of the " 32 "tan(dilation angle). Usually the dilation angle is not greater than " 33 "the friction angle, and it is between 0 and 90deg.");
36 "A SolidMechanicsHardening UserObject that defines hardening of the " 37 "weak-plane tensile strength. In physical situations this is positive " 38 "(and always must be greater than negative compressive-strength.");
40 "A SolidMechanicsHardening UserObject that defines " 41 "hardening of the weak-plane compressive strength. In " 42 "physical situations this is positive.");
46 "The cone vertex at shear-stress = 0 will be smoothed by " 47 "the given amount. Typical value is 0.1*cohesion");
48 params.
addParam<
bool>(
"perfect_guess",
50 "Provide a guess to the Newton-Raphson procedure " 51 "that is the result from perfect plasticity. With " 52 "severe hardening/softening this may be " 64 _small_smoother2(Utility::
pow<2>(getParam<
Real>(
"tip_smoother"))),
65 _perfect_guess(getParam<bool>(
"perfect_guess")),
74 mooseError(
"CappedWeakPlaneStressUpdate: Weak-plane friction and dilation angles must lie in " 77 mooseError(
"CappedWeakPlaneStressUpdate: Weak-plane friction angle must not be less than " 80 mooseError(
"CappedWeakPlaneStressUpdate: Weak-plane cohesion must not be negative");
82 mooseError(
"CappedWeakPlaneStressUpdate: Weak plane tensile strength plus compressive " 83 "strength must be greater than smoothing_tol");
102 const std::vector<Real> & ,
103 const std::vector<Real> & yf,
123 q =
std::sqrt(Utility::pow<2>(stress(0, 2)) + Utility::pow<2>(stress(1, 2)));
129 Epp = Eijkl(2, 2, 2, 2);
130 Eqq = Eijkl(0, 2, 0, 2);
138 const std::vector<Real> & ,
143 stress = stress_trial;
147 stress(0, 0) -= Eijkl(2, 2, 0, 0) * gaE /
_Epp * smoothed_q.
dg[0];
148 stress(1, 1) -= Eijkl(2, 2, 1, 1) * gaE /
_Epp * smoothed_q.
dg[0];
150 stress(2, 0) = stress(2, 1) = stress(0, 2) = stress(1, 2) = 0.0;
163 const std::vector<Real> & intnl,
164 std::vector<std::vector<Real>> & dintnl)
const 168 dintnl[0][1] = -1.0 /
_Eqq;
169 dintnl[1][0] = -1.0 /
_Epp;
184 bool compute_full_tangent_operator,
188 if (!compute_full_tangent_operator)
191 const Real Ezzzz = Eijkl(2, 2, 2, 2);
192 const Real Ezxzx = Eijkl(2, 0, 2, 0);
194 const Real dintnl0_dq = -1.0 / Ezxzx;
195 const Real dintnl0_dqt = 1.0 / Ezxzx;
196 const Real dintnl1_dp = -1.0 / Ezzzz;
197 const Real dintnl1_dpt = 1.0 / Ezzzz;
198 const Real dintnl1_dq =
200 const Real dintnl1_dqt =
205 const Real dpt_depii = Eijkl(2, 2, i, i);
206 cto(2, 2, i, i) =
_dp_dpt * dpt_depii;
207 const Real poisson_effect =
208 Eijkl(2, 2, 0, 0) / Ezzzz *
212 gaE * smoothed_q.
d2g_di[0][1] *
215 cto(0, 0, i, i) -= poisson_effect;
216 cto(1, 1, i, i) -= poisson_effect;
224 const Real poisson_effect =
225 -Eijkl(2, 2, 0, 0) / Ezzzz *
228 gaE * smoothed_q.
d2g_di[0][0] * (dintnl0_dqt + dintnl0_dq *
_dq_dqt) +
231 const Real dqt_dep20 = (q_trial == 0.0 ? 1.0 :
_in_trial02 / q_trial) * Eijkl(2, 0, 2, 0);
232 cto(2, 2, 2, 0) = cto(2, 2, 0, 2) =
_dp_dqt * dqt_dep20;
233 cto(0, 0, 2, 0) = cto(0, 0, 0, 2) = cto(1, 1, 2, 0) = cto(1, 1, 0, 2) =
234 poisson_effect * dqt_dep20;
238 cto(0, 2, 0, 2) = cto(2, 0, 0, 2) = cto(0, 2, 2, 0) = cto(2, 0, 2, 0) =
239 Eijkl(2, 0, 2, 0) * q / q_trial +
241 cto(1, 2, 0, 2) = cto(2, 1, 0, 2) = cto(1, 2, 2, 0) = cto(2, 1, 2, 0) =
245 const Real dqt_dep21 = (q_trial == 0.0 ? 1.0 :
_in_trial12 / q_trial) * Eijkl(2, 1, 2, 1);
246 cto(2, 2, 2, 1) = cto(2, 2, 1, 2) =
_dp_dqt * dqt_dep21;
247 cto(0, 0, 2, 1) = cto(0, 0, 1, 2) = cto(1, 1, 2, 1) = cto(1, 1, 1, 2) =
248 poisson_effect * dqt_dep21;
252 cto(0, 2, 1, 2) = cto(2, 0, 1, 2) = cto(0, 2, 2, 1) = cto(2, 0, 2, 1) =
254 cto(1, 2, 1, 2) = cto(2, 1, 1, 2) = cto(1, 2, 2, 1) = cto(2, 1, 2, 1) =
255 Eijkl(2, 1, 2, 1) * q / q_trial +
263 const std::vector<Real> & intnl,
264 std::vector<Real> & yf)
const 270 yf[1] = std::numeric_limits<Real>::lowest();
275 yf[2] = std::numeric_limits<Real>::lowest();
283 const std::vector<Real> & intnl,
284 std::vector<yieldAndFlow> & all_q)
const 290 all_q[1].f = std::numeric_limits<Real>::lowest();
294 all_q[2].f = std::numeric_limits<Real>::lowest();
301 all_q[1].df[0] = 1.0;
302 all_q[2].df[0] = -1.0;
306 all_q[0].df[1] = 1.0;
309 all_q[1].df[1] = 0.0;
310 all_q[2].df[1] = 0.0;
315 all_q[1].df_di[0] = 0.0;
316 all_q[2].df_di[0] = 0.0;
318 all_q[0].df_di[1] = 0.0;
325 all_q[1].dg[0] = 1.0;
326 all_q[2].dg[0] = -1.0;
329 all_q[0].dg[1] = 1.0;
332 all_q[1].dg[1] = 0.0;
333 all_q[2].dg[1] = 0.0;
338 all_q[1].d2g_di[0][0] = 0.0;
339 all_q[2].d2g_di[0][0] = 0.0;
341 all_q[0].d2g_di[0][1] = 0.0;
342 all_q[1].d2g_di[0][1] = 0.0;
343 all_q[2].d2g_di[0][1] = 0.0;
345 all_q[0].d2g_di[1][0] = 0.0;
346 all_q[1].d2g_di[1][0] = 0.0;
347 all_q[2].d2g_di[1][0] = 0.0;
349 all_q[0].d2g_di[1][1] = 0.0;
350 all_q[1].d2g_di[1][1] = 0.0;
351 all_q[2].d2g_di[1][1] = 0.0;
355 all_q[0].d2g[0][0] = 0.0;
356 all_q[1].d2g[0][0] = 0.0;
357 all_q[2].d2g[0][0] = 0.0;
359 all_q[0].d2g[0][1] = 0.0;
360 all_q[1].d2g[0][1] = 0.0;
361 all_q[2].d2g[0][1] = 0.0;
363 all_q[0].d2g[1][0] = 0.0;
364 all_q[1].d2g[1][0] = 0.0;
365 all_q[2].d2g[1][0] = 0.0;
368 all_q[0].d2g[1][1] = 0.0;
371 all_q[1].d2g[1][1] = 0.0;
372 all_q[2].d2g[1][1] = 0.0;
378 const std::vector<Real> & intnl_old,
382 std::vector<Real> & intnl)
const 398 const Real q_at_T = coh - tens * tanphi;
399 const Real q_at_C = coh + comp * tanphi;
401 if ((p_trial >= tens) && (q_trial <= q_at_T))
408 else if ((p_trial <= -comp) && (q_trial <= q_at_C))
419 const Real ga = (q_trial + p_trial * tanphi - coh) / (
_Eqq +
_Epp * tanphi * tanpsi);
420 if (ga <= 0 && p_trial <= tens && p_trial >= -comp)
430 q = q_trial -
_Eqq * ga;
431 if (q <= 0.0 && q_at_T <= 0.0)
437 gaE = (p_trial - p) / tanpsi;
439 else if (q <= q_at_T)
444 gaE = (p_trial - p) / tanpsi;
446 else if (q >= q_at_C)
451 if (p - p_trial <
_Epp * tanpsi * (q_trial - q) /
_Eqq)
456 gaE = (p - p_trial) / tanpsi;
461 p = p_trial -
_Epp * ga * tanpsi;
475 const std::vector<Real> & intnl_old,
476 std::vector<Real> & intnl)
const 478 intnl[0] = intnl_old[0] + (q_trial - q) /
_Eqq;
480 intnl[1] = intnl_old[1] + (p_trial - p) /
_Epp - (q_trial - q) * tanpsi /
_Eqq;
501 const Real q =
std::sqrt(Utility::pow<2>(stress(2, 0)) + Utility::pow<2>(stress(2, 1)));
504 deriv(2, 0) =
deriv(0, 2) = 0.5 * stress(2, 0) / q;
505 deriv(2, 1) =
deriv(1, 2) = 0.5 * stress(2, 1) / q;
521 const Real q =
std::sqrt(Utility::pow<2>(stress(2, 0)) + Utility::pow<2>(stress(2, 1)));
526 dq(2, 0) = dq(0, 2) = 0.25 * (stress(2, 0) + stress(0, 2)) / q;
527 dq(2, 1) = dq(1, 2) = 0.25 * (stress(2, 1) + stress(1, 2)) / q;
533 d2(i,
j,
k, l) = -dq(i,
j) * dq(
k, l) / q;
535 d2(0, 2, 0, 2) += 0.25 / q;
536 d2(0, 2, 2, 0) += 0.25 / q;
537 d2(2, 0, 0, 2) += 0.25 / q;
538 d2(2, 0, 2, 0) += 0.25 / q;
539 d2(1, 2, 1, 2) += 0.25 / q;
540 d2(1, 2, 2, 1) += 0.25 / q;
541 d2(2, 1, 1, 2) += 0.25 / q;
542 d2(2, 1, 2, 1) += 0.25 / q;
const SolidMechanicsHardeningModel & _cstrength
Hardening model for compressive strength.
virtual RankTwoTensor dqdstress(const RankTwoTensor &stress) const override
d(q)/d(stress) Derived classes must override this
std::vector< std::vector< Real > > d2g
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.
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...
virtual void yieldFunctionValues(Real p, Real q, const std::vector< Real > &intnl, std::vector< Real > &yf) const override
Computes the values of the yield functions, given p, q and intnl parameters.
Real _dp_dpt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
Real _dq_dpt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
Struct designed to hold info about a single yield function and its derivatives, as well as the flow d...
TwoParameterPlasticityStressUpdate performs the return-map algorithm and associated stress updates fo...
Real _dgaE_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
virtual RankFourTensor d2pdstress2(const RankTwoTensor &stress) const override
d2(p)/d(stress)/d(stress) Derived classes must override this
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
virtual void computePQ(const RankTwoTensor &stress, Real &p, Real &q) const override
Computes p and q, given stress.
Real _dgaE_dpt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
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.
virtual Real value(Real intnl) const
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.
virtual RankTwoTensor dpdstress(const RankTwoTensor &stress) const override
d(p)/d(stress) Derived classes must override this
const SolidMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
MaterialProperty< std::vector< Real > > & _intnl
internal parameters
Real _in_trial02
trial value of stress(0, 2)
virtual void setIntnlValues(Real p_trial, Real q_trial, Real p, Real q, const std::vector< Real > &intnl_old, std::vector< Real > &intnl) const override
Sets the internal parameters based on the trial values of p and q, their current values, and the old values of the internal parameters.
Real _dq_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
Real deriv(unsigned n, unsigned alpha, unsigned beta, Real x)
const SolidMechanicsHardeningModel & _tan_psi
Hardening model for tan(psi)
virtual void finalizeReturnProcess(const RankTwoTensor &rotation_increment) override
Derived classes may use this to perform calculations after the return-map process has completed succe...
StressReturnType
This allows some simplification in the return-map process.
static InputParameters validParams()
registerMooseObject("SolidMechanicsApp", CappedWeakPlaneStressUpdate)
const bool _perfect_guess
Initialize the NR proceedure from a guess coming from perfect plasticity.
Real _Epp
elasticity tensor in p direction
virtual void setIntnlDerivatives(Real p_trial, Real q_trial, Real p, Real q, const std::vector< Real > &intnl, std::vector< std::vector< Real >> &dintnl) const override
Sets the derivatives of internal parameters, based on the trial values of p and q, their current values, and the old values of the internal parameters.
const Real _smoothing_tol
Smoothing tolerance: edges of the yield surface get smoothed by this amount.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
virtual Real derivative(Real intnl) const
Hardening Model base class.
Real _in_q_trial
trial value of q
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Real _Eqq
elasticity tensor in q direction
enum CappedWeakPlaneStressUpdate::StressReturnType _stress_return_type
const SolidMechanicsHardeningModel & _tstrength
Hardening model for tensile strength.
void mooseError(Args &&... args) const
static InputParameters validParams()
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
std::vector< std::vector< Real > > d2g_di
CappedWeakPlaneStressUpdate(const InputParameters ¶meters)
virtual RankFourTensor d2qdstress2(const RankTwoTensor &stress) const override
d2(q)/d(stress)/d(stress) Derived classes must override this
const SolidMechanicsHardeningModel & _tan_phi
Hardening model for tan(phi)
CappedWeakPlaneStressUpdate performs the return-map algorithm and associated stress updates for plast...
virtual void initializeReturnProcess() override
Derived classes may use this to perform calculations before any return-map process is performed...
virtual void setStressAfterReturn(const RankTwoTensor &stress_trial, Real p_ok, Real q_ok, Real gaE, const std::vector< Real > &intnl, const yieldAndFlow &smoothed_q, const RankFourTensor &Eijkl, RankTwoTensor &stress) const override
Sets stress from the admissible parameters.
const Real _small_smoother2
The cone vertex is smoothed by this amount.
Real _dp_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
MooseUnits pow(const MooseUnits &, int)
static constexpr unsigned _tensor_dimensionality
Internal dimensionality of tensors (currently this is 3 throughout tensor_mechanics) ...
static const std::string k
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.
Real _in_trial12
trial value of stress(1, 2)