12 #include "libmesh/utility.h" 23 "A SolidMechanicsPlasticDruckerPrager UserObject that defines the " 24 "Drucker-Prager parameters (cohesion, friction angle and dilation angle)");
27 "A SolidMechanicsHardening UserObject that defines hardening of the " 28 "tensile strength. In physical situations this is positive (and always " 29 "must be greater than negative compressive-strength.");
31 "compressive_strength",
32 "A SolidMechanicsHardening UserObject that defines hardening of the " 33 "compressive strength. In physical situations this is positive.");
37 "The cone vertex at J2 = 0 will be smoothed by the given " 38 "amount. Typical value is 0.1*cohesion");
39 params.
addParam<
bool>(
"perfect_guess",
41 "Provide a guess to the Newton-Raphson procedure " 42 "that is the result from perfect plasticity. With " 43 "severe hardening/softening this may be " 45 params.
addParam<
bool>(
"small_dilation",
47 "If true, and if the trial stress exceeds the " 48 "tensile strength, then the user guarantees that " 49 "the returned stress will be independent of the " 50 "compressive strength.");
59 _small_smoother2(Utility::
pow<2>(getParam<
Real>(
"tip_smoother"))),
60 _perfect_guess(getParam<bool>(
"perfect_guess")),
62 _small_dilation(getParam<bool>(
"small_dilation")),
68 mooseError(
"CappedDruckerPragerStressUpdate: Tensile strength plus compressive strength must " 69 "be greater than smoothing_tol");
88 const std::vector<Real> & ,
89 const std::vector<Real> & yf,
114 Eqq = Eijkl(0, 1, 0, 1);
122 const std::vector<Real> & intnl,
123 std::vector<std::vector<Real>> & dintnl)
const 130 dintnl[0][1] = -1.0 /
_Eqq;
131 dintnl[1][0] = -1.0 /
_Epp;
132 dintnl[1][1] = tanpsi /
_Eqq - (q_trial - q) * dtanpsi * dintnl[0][1] /
_Eqq;
138 const std::vector<Real> & intnl,
139 std::vector<Real> & yf)
const 147 yf[1] = std::numeric_limits<Real>::lowest();
152 yf[2] = std::numeric_limits<Real>::lowest();
160 const std::vector<Real> & intnl,
161 std::vector<yieldAndFlow> & all_q)
const 177 all_q[1].f = std::numeric_limits<Real>::lowest();
181 all_q[2].f = std::numeric_limits<Real>::lowest();
187 all_q[0].df[0] = bbb;
188 all_q[1].df[0] = 1.0;
189 all_q[2].df[0] = -1.0;
193 all_q[0].df[1] = 1.0;
196 all_q[1].df[1] = 0.0;
197 all_q[2].df[1] = 0.0;
201 all_q[0].df_di[0] = p * dbbb - daaa;
202 all_q[1].df_di[0] = 0.0;
203 all_q[2].df_di[0] = 0.0;
205 all_q[0].df_di[1] = 0.0;
211 all_q[0].dg[0] = tanpsi;
212 all_q[1].dg[0] = 1.0;
213 all_q[2].dg[0] = -1.0;
216 all_q[0].dg[1] = 1.0;
219 all_q[1].dg[1] = 0.0;
220 all_q[2].dg[1] = 0.0;
224 all_q[0].d2g_di[0][0] = dtanpsi;
225 all_q[1].d2g_di[0][0] = 0.0;
226 all_q[2].d2g_di[0][0] = 0.0;
228 all_q[0].d2g_di[0][1] = 0.0;
229 all_q[1].d2g_di[0][1] = 0.0;
230 all_q[2].d2g_di[0][1] = 0.0;
232 all_q[0].d2g_di[1][0] = 0.0;
233 all_q[1].d2g_di[1][0] = 0.0;
234 all_q[2].d2g_di[1][0] = 0.0;
236 all_q[0].d2g_di[1][1] = 0.0;
237 all_q[1].d2g_di[1][1] = 0.0;
238 all_q[2].d2g_di[1][1] = 0.0;
242 all_q[0].d2g[0][0] = 0.0;
243 all_q[1].d2g[0][0] = 0.0;
244 all_q[2].d2g[0][0] = 0.0;
246 all_q[0].d2g[0][1] = 0.0;
247 all_q[1].d2g[0][1] = 0.0;
248 all_q[2].d2g[0][1] = 0.0;
250 all_q[0].d2g[1][0] = 0.0;
251 all_q[1].d2g[1][0] = 0.0;
252 all_q[2].d2g[1][0] = 0.0;
255 all_q[0].d2g[1][1] = 0.0;
258 all_q[1].d2g[1][1] = 0.0;
259 all_q[2].d2g[1][1] = 0.0;
265 const std::vector<Real> & intnl_old,
269 std::vector<Real> & intnl)
const 286 const Real q_at_T = coh - tens * tanphi;
287 const Real q_at_C = coh + comp * tanphi;
289 if ((p_trial >= tens) && (q_trial <= q_at_T))
296 else if ((p_trial <= -comp) && (q_trial <= q_at_C))
307 const Real ga = (q_trial + p_trial * tanphi - coh) / (
_Eqq +
_Epp * tanphi * tanpsi);
308 if (ga <= 0 && p_trial <= tens && p_trial >= -comp)
318 q = q_trial -
_Eqq * ga;
319 if (q <= 0.0 && q_at_T <= 0.0)
325 gaE = (p_trial - p) / tanpsi;
327 else if (q <= q_at_T)
332 gaE = (p_trial - p) / tanpsi;
334 else if (q >= q_at_C)
339 if (p - p_trial <
_Epp * tanpsi * (q_trial - q) /
_Eqq)
344 gaE = (p - p_trial) / tanpsi;
349 p = p_trial -
_Epp * ga * tanpsi;
363 const std::vector<Real> & intnl_old,
364 std::vector<Real> & intnl)
const 366 intnl[0] = intnl_old[0] + (q_trial - q) /
_Eqq;
369 intnl[1] = intnl_old[1] + (p_trial - p) /
_Epp - (q_trial - q) * tanpsi /
_Eqq;
377 const std::vector<Real> & ,
384 const Real p_trial = stress_trial.
trace();
434 bool compute_full_tangent_operator,
438 if (!compute_full_tangent_operator)
451 (i ==
j) * (1.0 / 3.0) *
453 cto(i,
j,
k, l) -= s_over_q(i,
j) * (
_Epp * (-
_dq_dpt) / 3.0 * (
k == l) +
457 if (smoothed_q.
dg[1] != 0.0)
469 mooseWarning(
"CappedDruckerPragerStressUpdate: Cannot invert 1+T in consistent tangent " 470 "operator computation at quadpoint ",
RankFourTensorTempl< Real > outerProduct(const RankTwoTensorTempl< Real > &b) const
const SolidMechanicsPlasticDruckerPrager & _dp
Hardening model for cohesion, friction and dilation angles.
RankTwoTensorTempl< Real > dsecondInvariant() const
RankFourTensorTempl< Real > d2secondInvariant() const
const Real _small_smoother2
The cone vertex is smoothed by this amount.
virtual RankFourTensor d2qdstress2(const RankTwoTensor &stress) const override
d2(q)/d(stress)/d(stress) Derived classes must override this
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.
void donlyB(Real intnl, int fd, Real &dbbb) const
Calculate d(bbb)/d(intnl) or d(bbb_flow)/d(intnl)
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...
void bothAB(Real intnl, Real &aaa, Real &bbb) const
Calculates aaa and bbb as a function of the internal parameter intnl.
Struct designed to hold info about a single yield function and its derivatives, as well as the flow d...
static InputParameters validParams()
TwoParameterPlasticityStressUpdate performs the return-map algorithm and associated stress updates fo...
Real secondInvariant() const
RankTwoTensorTempl< Real > dtrace() const
virtual void finalizeReturnProcess(const RankTwoTensor &rotation_increment) override
Derived classes may use this to perform calculations after the return-map process has completed succe...
virtual RankTwoTensor dqdstress(const RankTwoTensor &stress) const override
d(q)/d(stress) Derived classes must override this
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
const bool _small_dilation
If true, and if the trial stress exceeds the tensile strength, then the user gaurantees that the retu...
virtual Real value(Real intnl) const
void mooseWarning(Args &&... args) const
virtual RankFourTensor d2pdstress2(const RankTwoTensor &stress) const override
d2(p)/d(stress)/d(stress) Derived classes must override this
CappedDruckerPragerStressUpdate(const InputParameters ¶meters)
void onlyB(Real intnl, int fd, Real &bbb) const
Calculate bbb or bbb_flow.
Rate-independent non-associative Drucker Prager with hardening/softening.
Real _dq_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
virtual RankTwoTensor dpdstress(const RankTwoTensor &stress) const override
d(p)/d(stress) Derived classes must override this
static InputParameters validParams()
RankFourTensorTempl< T > transposeMajor() const
Real _Epp
elasticity tensor in p direction
const Real _smoothing_tol
Smoothing tolerance: edges of the yield surface get smoothed by this amount.
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...
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
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.
virtual Real derivative(Real intnl) const
const SolidMechanicsHardeningModel & _cstrength
Hardening model for compressive strength.
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.
Hardening Model base class.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Real _Eqq
elasticity tensor in q direction
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 SolidMechanicsHardeningModel & _tstrength
Hardening model for tensile strength.
virtual void initializeReturnProcess() override
Derived classes may use this to perform calculations before any return-map process is performed...
StressReturnType
This allows some simplification in the return-map process.
void mooseError(Args &&... args) const
CappedDruckerPragerStressUpdate performs the return-map algorithm and associated stress updates for p...
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.
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.
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
virtual void computePQ(const RankTwoTensor &stress, Real &p, Real &q) const override
Computes p and q, given stress.
registerMooseObject("SolidMechanicsApp", CappedDruckerPragerStressUpdate)
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 bool _perfect_guess
Initialize the NR proceedure from a guess coming from perfect plasticity.
Real _in_q_trial
trial value of q
Real _dp_dqt
derivative of Variable with respect to trial variable (used in consistent-tangent-operator calculatio...
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...
RankFourTensorTempl< T > invSymm() const
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
const Elem *const & _current_elem
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.