Go to the documentation of this file.
12 #include "libmesh/utility.h"
22 params.addClassDescription(
"Capped Drucker-Prager plasticity stress calculator");
23 params.addRequiredParam<UserObjectName>(
25 "A TensorMechanicsPlasticDruckerPrager UserObject that defines the "
26 "Drucker-Prager parameters (cohesion, friction angle and dilation angle)");
27 params.addRequiredParam<UserObjectName>(
29 "A TensorMechanicsHardening UserObject that defines hardening of the "
30 "tensile strength. In physical situations this is positive (and always "
31 "must be greater than negative compressive-strength.");
32 params.addRequiredParam<UserObjectName>(
33 "compressive_strength",
34 "A TensorMechanicsHardening UserObject that defines hardening of the "
35 "compressive strength. In physical situations this is positive.");
36 params.addRequiredRangeCheckedParam<Real>(
39 "The cone vertex at J2 = 0 will be smoothed by the given "
40 "amount. Typical value is 0.1*cohesion");
41 params.addParam<
bool>(
"perfect_guess",
43 "Provide a guess to the Newton-Raphson procedure "
44 "that is the result from perfect plasticity. With "
45 "severe hardening/softening this may be "
47 params.addParam<
bool>(
"small_dilation",
49 "If true, and if the trial stress exceeds the "
50 "tensile strength, then the user guarantees that "
51 "the returned stress will be independent of the "
52 "compressive strength.");
61 _small_smoother2(Utility::
pow<2>(getParam<Real>(
"tip_smoother"))),
62 _perfect_guess(getParam<bool>(
"perfect_guess")),
64 _small_dilation(getParam<bool>(
"small_dilation")),
70 mooseError(
"CappedDruckerPragerStressUpdate: Tensile strength plus compressive strength must "
71 "be greater than smoothing_tol");
90 const std::vector<Real> & ,
91 const std::vector<Real> & yf,
107 q = std::sqrt(stress.secondInvariant());
115 Epp = Eijkl.sum3x3();
116 Eqq = Eijkl(0, 1, 0, 1);
124 const std::vector<Real> & intnl,
125 std::vector<std::vector<Real>> & dintnl)
const
132 dintnl[0][1] = -1.0 /
_Eqq;
133 dintnl[1][0] = -1.0 /
_Epp;
134 dintnl[1][1] = tanpsi /
_Eqq - (q_trial - q) * dtanpsi * dintnl[0][1] /
_Eqq;
140 const std::vector<Real> & intnl,
141 std::vector<Real> & yf)
const
149 yf[1] = std::numeric_limits<Real>::lowest();
154 yf[2] = std::numeric_limits<Real>::lowest();
162 const std::vector<Real> & intnl,
163 std::vector<yieldAndFlow> & all_q)
const
177 all_q[0].f = std::sqrt(Utility::pow<2>(q) +
_small_smoother2) + p * bbb - aaa;
179 all_q[1].f = std::numeric_limits<Real>::lowest();
183 all_q[2].f = std::numeric_limits<Real>::lowest();
189 all_q[0].df[0] = bbb;
190 all_q[1].df[0] = 1.0;
191 all_q[2].df[0] = -1.0;
195 all_q[0].df[1] = 1.0;
198 all_q[1].df[1] = 0.0;
199 all_q[2].df[1] = 0.0;
203 all_q[0].df_di[0] = p * dbbb - daaa;
204 all_q[1].df_di[0] = 0.0;
205 all_q[2].df_di[0] = 0.0;
207 all_q[0].df_di[1] = 0.0;
213 all_q[0].dg[0] = tanpsi;
214 all_q[1].dg[0] = 1.0;
215 all_q[2].dg[0] = -1.0;
218 all_q[0].dg[1] = 1.0;
221 all_q[1].dg[1] = 0.0;
222 all_q[2].dg[1] = 0.0;
226 all_q[0].d2g_di[0][0] = dtanpsi;
227 all_q[1].d2g_di[0][0] = 0.0;
228 all_q[2].d2g_di[0][0] = 0.0;
230 all_q[0].d2g_di[0][1] = 0.0;
231 all_q[1].d2g_di[0][1] = 0.0;
232 all_q[2].d2g_di[0][1] = 0.0;
234 all_q[0].d2g_di[1][0] = 0.0;
235 all_q[1].d2g_di[1][0] = 0.0;
236 all_q[2].d2g_di[1][0] = 0.0;
238 all_q[0].d2g_di[1][1] = 0.0;
239 all_q[1].d2g_di[1][1] = 0.0;
240 all_q[2].d2g_di[1][1] = 0.0;
244 all_q[0].d2g[0][0] = 0.0;
245 all_q[1].d2g[0][0] = 0.0;
246 all_q[2].d2g[0][0] = 0.0;
248 all_q[0].d2g[0][1] = 0.0;
249 all_q[1].d2g[0][1] = 0.0;
250 all_q[2].d2g[0][1] = 0.0;
252 all_q[0].d2g[1][0] = 0.0;
253 all_q[1].d2g[1][0] = 0.0;
254 all_q[2].d2g[1][0] = 0.0;
257 all_q[0].d2g[1][1] = 0.0;
260 all_q[1].d2g[1][1] = 0.0;
261 all_q[2].d2g[1][1] = 0.0;
267 const std::vector<Real> & intnl_old,
271 std::vector<Real> & intnl)
const
288 const Real q_at_T = coh - tens * tanphi;
289 const Real q_at_C = coh + comp * tanphi;
291 if ((p_trial >= tens) && (q_trial <= q_at_T))
298 else if ((p_trial <= -comp) && (q_trial <= q_at_C))
309 const Real ga = (q_trial + p_trial * tanphi - coh) / (
_Eqq +
_Epp * tanphi * tanpsi);
310 if (ga <= 0 && p_trial <= tens && p_trial >= -comp)
320 q = q_trial -
_Eqq * ga;
321 if (q <= 0.0 && q_at_T <= 0.0)
327 gaE = (p_trial - p) / tanpsi;
329 else if (q <= q_at_T)
334 gaE = (p_trial - p) / tanpsi;
336 else if (q >= q_at_C)
341 if (p - p_trial <
_Epp * tanpsi * (q_trial - q) /
_Eqq)
346 gaE = (p - p_trial) / tanpsi;
351 p = p_trial -
_Epp * ga * tanpsi;
365 const std::vector<Real> & intnl_old,
366 std::vector<Real> & intnl)
const
368 intnl[0] = intnl_old[0] + (q_trial - q) /
_Eqq;
371 intnl[1] = intnl_old[1] + (p_trial - p) /
_Epp - (q_trial - q) * tanpsi /
_Eqq;
379 const std::vector<Real> & ,
386 const Real p_trial = stress_trial.trace();
396 return stress.dtrace();
408 const Real j2 = stress.secondInvariant();
411 return 0.5 * stress.dsecondInvariant() / std::sqrt(j2);
417 const Real j2 = stress.secondInvariant();
422 return -0.25 * dj2.outerProduct(dj2) /
std::pow(j2, 1.5) +
423 0.5 * stress.d2secondInvariant() / std::sqrt(j2);
436 bool compute_full_tangent_operator,
440 if (!compute_full_tangent_operator)
445 : (stress - stress.trace() *
RankTwoTensor(RankTwoTensor::initIdentity) / 3.0) / q);
453 (i == j) * (1.0 / 3.0) *
455 cto(i, j, k, l) -= s_over_q(i, j) * (
_Epp * (-
_dq_dpt) / 3.0 * (k == l) +
459 if (smoothed_q.
dg[1] != 0.0)
465 inv = inv.transposeMajor().invSymm();
467 catch (
const MooseException & e)
471 mooseWarning(
"CappedDruckerPragerStressUpdate: Cannot invert 1+T in consistent tangent "
472 "operator computation at quadpoint ",
475 _current_elem->id());
478 cto = (cto.transposeMajor() * inv).transposeMajor();
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.
const TensorMechanicsHardeningModel & _cstrength
Hardening model for compressive strength.
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 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.
virtual void finalizeReturnProcess(const RankTwoTensor &rotation_increment) override
Derived classes may use this to perform calculations after the return-map process has completed succe...
static InputParameters validParams()
StressReturnType
This allows some simplification in the return-map process.
Struct designed to hold info about a single yield function and its derivatives, as well as the flow d...
virtual RankFourTensor d2qdstress2(const RankTwoTensor &stress) const override
d2(q)/d(stress)/d(stress) Derived classes must override this
const Real _smoothing_tol
Smoothing tolerance: edges of the yield surface get smoothed by this amount.
void onlyB(Real intnl, int fd, Real &bbb) const
Calculate bbb or bbb_flow.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
void donlyB(Real intnl, int fd, Real &dbbb) const
Calculate d(bbb)/d(intnl) or d(bbb_flow)/d(intnl)
Rate-independent non-associative Drucker Prager with hardening/softening.
virtual RankFourTensor d2pdstress2(const RankTwoTensor &stress) const override
d2(p)/d(stress)/d(stress) Derived classes must override this
static InputParameters validParams()
Real _dp_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
Real _Epp
elasticity tensor in p direction
CappedDruckerPragerStressUpdate performs the return-map algorithm and associated stress updates for p...
virtual void initializeReturnProcess() override
Derived classes may use this to perform calculations before any return-map process is performed,...
Real _dq_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
const bool _small_dilation
If true, and if the trial stress exceeds the tensile strength, then the user gaurantees that the retu...
TwoParameterPlasticityStressUpdate performs the return-map algorithm and associated stress updates fo...
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.
enum CappedDruckerPragerStressUpdate::StressReturnType _stress_return_type
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.
const TensorMechanicsPlasticDruckerPrager & _dp
Hardening model for cohesion, friction and dilation angles.
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.
defineLegacyParams(CappedDruckerPragerStressUpdate)
virtual RankTwoTensor dqdstress(const RankTwoTensor &stress) const override
d(q)/d(stress) Derived classes must override this
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,...
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 _in_q_trial
trial value of q
virtual Real derivative(Real intnl) const
virtual void computePQ(const RankTwoTensor &stress, Real &p, Real &q) const override
Computes p and q, given stress.
virtual RankTwoTensor dpdstress(const RankTwoTensor &stress) const override
d(p)/d(stress) Derived classes must override this
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,...
const Real _small_smoother2
The cone vertex is smoothed by this amount.
virtual Real value(Real intnl) const
RankTwoTensorTempl< Real > RankTwoTensor
const bool _perfect_guess
Initialize the NR proceedure from a guess coming from perfect plasticity.
const TensorMechanicsHardeningModel & _tstrength
Hardening model for tensile strength.
Real _dq_dpt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
Real _Eqq
elasticity tensor in q direction
registerMooseObject("TensorMechanicsApp", CappedDruckerPragerStressUpdate)
void dbothAB(Real intnl, Real &daaa, Real &dbbb) const
Calculates d(aaa)/d(intnl) and d(bbb)/d(intnl) as a function of the internal parameter intnl.
CappedDruckerPragerStressUpdate(const InputParameters ¶meters)
void bothAB(Real intnl, Real &aaa, Real &bbb) const
Calculates aaa and bbb as a function of the internal parameter intnl.
RankFourTensorTempl< Real > RankFourTensor
constexpr static unsigned _tensor_dimensionality
Internal dimensionality of tensors (currently this is 3 throughout tensor_mechanics)
Hardening Model base class.