www.mooseframework.org
TensorMechanicsPlasticDruckerPragerHyperbolic.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
8 #include "libmesh/utility.h"
9 
10 template <>
11 InputParameters
13 {
14  InputParameters params = validParams<TensorMechanicsPlasticDruckerPrager>();
15  params.addParam<bool>("use_custom_returnMap",
16  true,
17  "Whether to use the custom returnMap "
18  "algorithm. Set to true if you are using "
19  "isotropic elasticity.");
20  params.addParam<bool>("use_custom_cto",
21  true,
22  "Whether to use the custom consistent tangent "
23  "operator computations. Set to true if you are "
24  "using isotropic elasticity.");
25  params.addClassDescription("J2 plasticity, associative, with hardening");
26  params.addRangeCheckedParam<Real>("smoother",
27  0.0,
28  "smoother>=0",
29  "The cone vertex at J2=0 is smoothed. The maximum mean "
30  "stress possible, which is Cohesion*Cot(friction_angle) for "
31  "smoother=0, becomes (Cohesion - "
32  "smoother/3)*Cot(friction_angle). This is a non-hardening "
33  "parameter");
34  params.addRangeCheckedParam<unsigned>(
35  "max_iterations",
36  40,
37  "max_iterations>0",
38  "Maximum iterations to use in the custom return map function");
39  params.addClassDescription(
40  "Non-associative Drucker Prager plasticity with hyperbolic smoothing of the cone tip.");
41  return params;
42 }
43 
45  const InputParameters & parameters)
47  _smoother2(Utility::pow<2>(getParam<Real>("smoother"))),
48  _use_custom_returnMap(getParam<bool>("use_custom_returnMap")),
49  _use_custom_cto(getParam<bool>("use_custom_cto")),
50  _max_iters(getParam<unsigned>("max_iterations"))
51 {
52 }
53 
54 Real
56  Real intnl) const
57 {
58  Real aaa;
59  Real bbb;
60  bothAB(intnl, aaa, bbb);
61  return std::sqrt(stress.secondInvariant() + _smoother2) + stress.trace() * bbb - aaa;
62 }
63 
64 RankTwoTensor
65 TensorMechanicsPlasticDruckerPragerHyperbolic::df_dsig(const RankTwoTensor & stress, Real bbb) const
66 {
67  return 0.5 * stress.dsecondInvariant() / std::sqrt(stress.secondInvariant() + _smoother2) +
68  stress.dtrace() * bbb;
69 }
70 
71 RankFourTensor
73  Real /*intnl*/) const
74 {
75  RankFourTensor dr_dstress;
76  dr_dstress = 0.5 * stress.d2secondInvariant() / std::sqrt(stress.secondInvariant() + _smoother2);
77  dr_dstress += -0.5 * 0.5 * stress.dsecondInvariant().outerProduct(stress.dsecondInvariant()) /
78  std::pow(stress.secondInvariant() + _smoother2, 1.5);
79  return dr_dstress;
80 }
81 
82 std::string
84 {
85  return "DruckerPragerHyperbolic";
86 }
87 
88 bool
90  Real intnl_old,
91  const RankFourTensor & E_ijkl,
92  Real ep_plastic_tolerance,
93  RankTwoTensor & returned_stress,
94  Real & returned_intnl,
95  std::vector<Real> & dpm,
96  RankTwoTensor & delta_dp,
97  std::vector<Real> & yf,
98  bool & trial_stress_inadmissible) const
99 {
100  if (!(_use_custom_returnMap))
101  return TensorMechanicsPlasticModel::returnMap(trial_stress,
102  intnl_old,
103  E_ijkl,
104  ep_plastic_tolerance,
105  returned_stress,
106  returned_intnl,
107  dpm,
108  delta_dp,
109  yf,
110  trial_stress_inadmissible);
111 
112  yf.resize(1);
113 
114  yf[0] = yieldFunction(trial_stress, intnl_old);
115 
116  if (yf[0] < _f_tol)
117  {
118  // the trial_stress is admissible
119  trial_stress_inadmissible = false;
120  return true;
121  }
122 
123  trial_stress_inadmissible = true;
124  const Real mu = E_ijkl(0, 1, 0, 1);
125  const Real lambda = E_ijkl(0, 0, 0, 0) - 2.0 * mu;
126  const Real bulky = 3.0 * lambda + 2.0 * mu;
127  const Real Tr_trial = trial_stress.trace();
128  const Real J2trial = trial_stress.secondInvariant();
129 
130  // Perform a Newton-Raphson to find dpm when
131  // residual = (1 + dpm*mu/ll)sqrt(ll^2 - s^2) - sqrt(J2trial) = 0, with s=smoother
132  // with ll = sqrt(J2 + s^2) = aaa - bbb*Tr(stress) = aaa - bbb(Tr_trial - p*3*bulky*bbb_flow)
133  Real aaa;
134  Real daaa;
135  Real bbb;
136  Real dbbb;
137  Real bbb_flow;
138  Real dbbb_flow;
139  Real ll;
140  Real dll;
141  Real residual;
142  Real jac;
143  dpm[0] = 0;
144  unsigned int iter = 0;
145  do
146  {
147  bothAB(intnl_old + dpm[0], aaa, bbb);
148  dbothAB(intnl_old + dpm[0], daaa, dbbb);
149  onlyB(intnl_old + dpm[0], dilation, bbb_flow);
150  donlyB(intnl_old + dpm[0], dilation, dbbb_flow);
151  ll = aaa - bbb * (Tr_trial - dpm[0] * bulky * 3 * bbb_flow);
152  dll = daaa - dbbb * (Tr_trial - dpm[0] * bulky * 3 * bbb_flow) +
153  bbb * bulky * 3 * (bbb_flow + dpm[0] * dbbb_flow);
154  residual = bbb * (Tr_trial - dpm[0] * bulky * 3 * bbb_flow) - aaa +
155  std::sqrt(J2trial / Utility::pow<2>(1 + dpm[0] * mu / ll) + _smoother2);
156  jac = dbbb * (Tr_trial - dpm[0] * bulky * 3 * bbb_flow) -
157  bbb * bulky * 3 * (bbb_flow + dpm[0] * dbbb_flow) - daaa +
158  0.5 * J2trial * (-2.0) * (mu / ll - dpm[0] * mu * dll / ll / ll) /
159  Utility::pow<3>(1 + dpm[0] * mu / ll) /
160  std::sqrt(J2trial / Utility::pow<2>(1.0 + dpm[0] * mu / ll) + _smoother2);
161  dpm[0] += -residual / jac;
162  if (iter > _max_iters) // not converging
163  return false;
164  iter++;
165  } while (residual * residual > _f_tol * _f_tol);
166 
167  // set the returned values
168  yf[0] = 0;
169  returned_intnl = intnl_old + dpm[0];
170 
171  bothAB(returned_intnl, aaa, bbb);
172  onlyB(returned_intnl, dilation, bbb_flow);
173  ll = aaa - bbb * (Tr_trial - dpm[0] * bulky * 3.0 * bbb_flow);
174  returned_stress =
175  trial_stress.deviatoric() / (1.0 + dpm[0] * mu / ll); // this is the deviatoric part only
176 
177  RankTwoTensor rij = 0.5 * returned_stress.deviatoric() /
178  ll; // this is the derivatoric part the flow potential only
179 
180  // form the returned stress and the full flow potential
181  const Real returned_trace_over_3 = (aaa - ll) / bbb / 3.0;
182  for (unsigned i = 0; i < 3; ++i)
183  {
184  returned_stress(i, i) += returned_trace_over_3;
185  rij(i, i) += bbb_flow;
186  }
187 
188  delta_dp = rij * dpm[0];
189 
190  return true;
191 }
192 
193 RankFourTensor
195  const RankTwoTensor & trial_stress,
196  Real intnl_old,
197  const RankTwoTensor & stress,
198  Real intnl,
199  const RankFourTensor & E_ijkl,
200  const std::vector<Real> & cumulative_pm) const
201 {
202  if (!_use_custom_cto)
204  trial_stress, intnl_old, stress, intnl, E_ijkl, cumulative_pm);
205 
206  const Real mu = E_ijkl(0, 1, 0, 1);
207  const Real la = E_ijkl(0, 0, 0, 0) - 2.0 * mu;
208  const Real bulky = 3.0 * la + 2.0 * mu;
209  Real bbb;
210  onlyB(intnl, friction, bbb);
211  Real bbb_flow;
212  onlyB(intnl, dilation, bbb_flow);
213  Real dbbb_flow;
214  donlyB(intnl, dilation, dbbb_flow);
215  const Real bbb_flow_mod = bbb_flow + cumulative_pm[0] * dbbb_flow;
216  const Real J2 = stress.secondInvariant();
217  const RankTwoTensor sij = stress.deviatoric();
218  const Real sq = std::sqrt(J2 + _smoother2);
219 
220  const Real one_over_h =
221  1.0 / (-dyieldFunction_dintnl(stress, intnl) + mu * J2 / Utility::pow<2>(sq) +
222  3.0 * bbb * bbb_flow_mod * bulky); // correct up to hard
223 
224  const RankTwoTensor df_dsig_timesE =
225  mu * sij / sq + bulky * bbb * RankTwoTensor(RankTwoTensor::initIdentity); // correct
226  const RankTwoTensor rmod_timesE =
227  mu * sij / sq +
228  bulky * bbb_flow_mod * RankTwoTensor(RankTwoTensor::initIdentity); // correct up to hard
229 
230  const RankFourTensor coeff_ep =
231  E_ijkl - one_over_h * rmod_timesE.outerProduct(df_dsig_timesE); // correct
232 
233  const RankFourTensor dr_dsig_timesE = -0.5 * mu * sij.outerProduct(sij) / Utility::pow<3>(sq) +
234  mu * stress.d2secondInvariant() / sq; // correct
235  const RankTwoTensor df_dsig_E_dr_dsig =
236  0.5 * mu * _smoother2 * sij / Utility::pow<4>(sq); // correct
237 
238  const RankFourTensor coeff_si =
239  RankFourTensor(RankFourTensor::initIdentitySymmetricFour) +
240  cumulative_pm[0] *
241  (dr_dsig_timesE - one_over_h * rmod_timesE.outerProduct(df_dsig_E_dr_dsig));
242 
243  RankFourTensor s_inv;
244  try
245  {
246  s_inv = coeff_si.invSymm();
247  }
248  catch (MooseException & e)
249  {
250  return coeff_ep; // when coeff_si is singular return the "linear" tangent operator
251  }
252 
253  return s_inv * coeff_ep;
254 }
255 
256 bool
258 {
259  return _use_custom_returnMap;
260 }
261 
262 bool
264 {
265  return _use_custom_cto;
266 }
virtual Real dyieldFunction_dintnl(const RankTwoTensor &stress, Real intnl) const override
The derivative of yield function with respect to the internal parameter.
const Real _smoother2
smoothing parameter for the cone&#39;s tip
virtual RankFourTensor consistentTangentOperator(const RankTwoTensor &trial_stress, Real intnl_old, const RankTwoTensor &stress, Real intnl, const RankFourTensor &E_ijkl, const std::vector< Real > &cumulative_pm) const override
Calculates a custom consistent tangent operator.
virtual RankFourTensor consistentTangentOperator(const RankTwoTensor &trial_stress, Real intnl_old, const RankTwoTensor &stress, Real intnl, const RankFourTensor &E_ijkl, const std::vector< Real > &cumulative_pm) const
Calculates a custom consistent tangent operator.
virtual bool useCustomCTO() const override
Returns false. You will want to override this in your derived class if you write a custom consistent ...
Real yieldFunction(const RankTwoTensor &stress, Real intnl) const override
The following functions are what you should override when building single-plasticity models...
virtual bool returnMap(const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const
Performs a custom return-map.
InputParameters validParams< TensorMechanicsPlasticDruckerPrager >()
void bothAB(Real intnl, Real &aaa, Real &bbb) const
Calculates aaa and bbb as a function of the internal parameter intnl.
RankFourTensor dflowPotential_dstress(const RankTwoTensor &stress, Real intnl) const override
The derivative of the flow potential with respect to stress.
virtual bool returnMap(const RankTwoTensor &trial_stress, Real intnl_old, const RankFourTensor &E_ijkl, Real ep_plastic_tolerance, RankTwoTensor &returned_stress, Real &returned_intnl, std::vector< Real > &dpm, RankTwoTensor &delta_dp, std::vector< Real > &yf, bool &trial_stress_inadmissible) const override
Performs a custom return-map.
virtual RankTwoTensor df_dsig(const RankTwoTensor &stress, Real bbb) const override
Function that&#39;s used in dyieldFunction_dstress and flowPotential.
Rate-independent non-associative Drucker Prager with hardening/softening.
const unsigned _max_iters
max iters for custom return map loop
const bool _use_custom_cto
Whether to use the custom consistent tangent operator calculation.
virtual bool useCustomReturnMap() const override
Returns false. You will want to override this in your derived class if you write a custom returnMap f...
void donlyB(Real intnl, int fd, Real &dbbb) const
Calculate d(bbb)/d(intnl) or d(bbb_flow)/d(intnl)
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...
InputParameters validParams< TensorMechanicsPlasticDruckerPragerHyperbolic >()
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
const Real _f_tol
Tolerance on yield function.
void onlyB(Real intnl, int fd, Real &bbb) const
Calculate bbb or bbb_flow.
const bool _use_custom_returnMap
whether to use the custom returnMap function