www.mooseframework.org
TensorMechanicsPlasticDruckerPragerHyperbolic.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
11 #include "RankFourTensor.h"
12 #include "libmesh/utility.h"
13 
15 
17 
18 InputParameters
20 {
21  InputParameters params = TensorMechanicsPlasticDruckerPrager::validParams();
22  params.addParam<bool>("use_custom_returnMap",
23  true,
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",
28  true,
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",
34  0.0,
35  "smoother>=0",
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 "
40  "parameter");
41  params.addRangeCheckedParam<unsigned>(
42  "max_iterations",
43  40,
44  "max_iterations>0",
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.");
48  return params;
49 }
50 
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"))
58 {
59 }
60 
61 Real
63  Real intnl) const
64 {
65  Real aaa;
66  Real bbb;
67  bothAB(intnl, aaa, bbb);
68  return std::sqrt(stress.secondInvariant() + _smoother2) + stress.trace() * bbb - aaa;
69 }
70 
73 {
74  return 0.5 * stress.dsecondInvariant() / std::sqrt(stress.secondInvariant() + _smoother2) +
75  stress.dtrace() * bbb;
76 }
77 
80  Real /*intnl*/) const
81 {
82  RankFourTensor dr_dstress;
83  dr_dstress = 0.5 * stress.d2secondInvariant() / std::sqrt(stress.secondInvariant() + _smoother2);
84  dr_dstress += -0.5 * 0.5 * stress.dsecondInvariant().outerProduct(stress.dsecondInvariant()) /
85  std::pow(stress.secondInvariant() + _smoother2, 1.5);
86  return dr_dstress;
87 }
88 
89 std::string
91 {
92  return "DruckerPragerHyperbolic";
93 }
94 
95 bool
97  Real intnl_old,
98  const RankFourTensor & E_ijkl,
99  Real ep_plastic_tolerance,
100  RankTwoTensor & returned_stress,
101  Real & returned_intnl,
102  std::vector<Real> & dpm,
103  RankTwoTensor & delta_dp,
104  std::vector<Real> & yf,
105  bool & trial_stress_inadmissible) const
106 {
107  if (!(_use_custom_returnMap))
108  return TensorMechanicsPlasticModel::returnMap(trial_stress,
109  intnl_old,
110  E_ijkl,
111  ep_plastic_tolerance,
112  returned_stress,
113  returned_intnl,
114  dpm,
115  delta_dp,
116  yf,
117  trial_stress_inadmissible);
118 
119  yf.resize(1);
120 
121  yf[0] = yieldFunction(trial_stress, intnl_old);
122 
123  if (yf[0] < _f_tol)
124  {
125  // the trial_stress is admissible
126  trial_stress_inadmissible = false;
127  return true;
128  }
129 
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();
136 
137  // Perform a Newton-Raphson to find dpm when
138  // residual = (1 + dpm*mu/ll)sqrt(ll^2 - s^2) - sqrt(J2trial) = 0, with s=smoother
139  // with ll = sqrt(J2 + s^2) = aaa - bbb*Tr(stress) = aaa - bbb(Tr_trial - p*3*bulky*bbb_flow)
140  Real aaa;
141  Real daaa;
142  Real bbb;
143  Real dbbb;
144  Real bbb_flow;
145  Real dbbb_flow;
146  Real ll;
147  Real dll;
148  Real residual;
149  Real jac;
150  dpm[0] = 0;
151  unsigned int iter = 0;
152  do
153  {
154  bothAB(intnl_old + dpm[0], aaa, bbb);
155  dbothAB(intnl_old + dpm[0], daaa, dbbb);
156  onlyB(intnl_old + dpm[0], dilation, bbb_flow);
157  donlyB(intnl_old + dpm[0], dilation, dbbb_flow);
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;
169  if (iter > _max_iters) // not converging
170  return false;
171  iter++;
172  } while (residual * residual > _f_tol * _f_tol);
173 
174  // set the returned values
175  yf[0] = 0;
176  returned_intnl = intnl_old + dpm[0];
177 
178  bothAB(returned_intnl, aaa, bbb);
179  onlyB(returned_intnl, dilation, bbb_flow);
180  ll = aaa - bbb * (Tr_trial - dpm[0] * bulky * 3.0 * bbb_flow);
181  returned_stress =
182  trial_stress.deviatoric() / (1.0 + dpm[0] * mu / ll); // this is the deviatoric part only
183 
184  RankTwoTensor rij = 0.5 * returned_stress.deviatoric() /
185  ll; // this is the derivatoric part the flow potential only
186 
187  // form the returned stress and the full flow potential
188  const Real returned_trace_over_3 = (aaa - ll) / bbb / 3.0;
189  for (unsigned i = 0; i < 3; ++i)
190  {
191  returned_stress(i, i) += returned_trace_over_3;
192  rij(i, i) += bbb_flow;
193  }
194 
195  delta_dp = rij * dpm[0];
196 
197  return true;
198 }
199 
202  const RankTwoTensor & trial_stress,
203  Real intnl_old,
204  const RankTwoTensor & stress,
205  Real intnl,
206  const RankFourTensor & E_ijkl,
207  const std::vector<Real> & cumulative_pm) const
208 {
209  if (!_use_custom_cto)
211  trial_stress, intnl_old, stress, intnl, E_ijkl, cumulative_pm);
212 
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;
216  Real bbb;
217  onlyB(intnl, friction, bbb);
218  Real bbb_flow;
219  onlyB(intnl, dilation, bbb_flow);
220  Real dbbb_flow;
221  donlyB(intnl, dilation, dbbb_flow);
222  const Real bbb_flow_mod = bbb_flow + cumulative_pm[0] * dbbb_flow;
223  const Real J2 = stress.secondInvariant();
224  const RankTwoTensor sij = stress.deviatoric();
225  const Real sq = std::sqrt(J2 + _smoother2);
226 
227  const Real one_over_h =
228  1.0 / (-dyieldFunction_dintnl(stress, intnl) + mu * J2 / Utility::pow<2>(sq) +
229  3.0 * bbb * bbb_flow_mod * bulky); // correct up to hard
230 
231  const RankTwoTensor df_dsig_timesE =
232  mu * sij / sq + bulky * bbb * RankTwoTensor(RankTwoTensor::initIdentity); // correct
233  const RankTwoTensor rmod_timesE =
234  mu * sij / sq +
235  bulky * bbb_flow_mod * RankTwoTensor(RankTwoTensor::initIdentity); // correct up to hard
236 
237  const RankFourTensor coeff_ep =
238  E_ijkl - one_over_h * rmod_timesE.outerProduct(df_dsig_timesE); // correct
239 
240  const RankFourTensor dr_dsig_timesE = -0.5 * mu * sij.outerProduct(sij) / Utility::pow<3>(sq) +
241  mu * stress.d2secondInvariant() / sq; // correct
242  const RankTwoTensor df_dsig_E_dr_dsig =
243  0.5 * mu * _smoother2 * sij / Utility::pow<4>(sq); // correct
244 
245  const RankFourTensor coeff_si =
246  RankFourTensor(RankFourTensor::initIdentitySymmetricFour) +
247  cumulative_pm[0] *
248  (dr_dsig_timesE - one_over_h * rmod_timesE.outerProduct(df_dsig_E_dr_dsig));
249 
250  RankFourTensor s_inv;
251  try
252  {
253  s_inv = coeff_si.invSymm();
254  }
255  catch (MooseException & e)
256  {
257  return coeff_ep; // when coeff_si is singular return the "linear" tangent operator
258  }
259 
260  return s_inv * coeff_ep;
261 }
262 
263 bool
265 {
266  return _use_custom_returnMap;
267 }
268 
269 bool
271 {
272  return _use_custom_cto;
273 }
TensorMechanicsPlasticDruckerPragerHyperbolic::useCustomCTO
virtual bool useCustomCTO() const override
Returns false. You will want to override this in your derived class if you write a custom consistent ...
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.C:270
TensorMechanicsPlasticDruckerPragerHyperbolic::returnMap
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.
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.C:96
TensorMechanicsPlasticDruckerPrager::onlyB
void onlyB(Real intnl, int fd, Real &bbb) const
Calculate bbb or bbb_flow.
Definition: TensorMechanicsPlasticDruckerPrager.C:161
pow
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
Definition: ExpressionBuilder.h:673
TensorMechanicsPlasticDruckerPrager::donlyB
void donlyB(Real intnl, int fd, Real &dbbb) const
Calculate d(bbb)/d(intnl) or d(bbb_flow)/d(intnl)
Definition: TensorMechanicsPlasticDruckerPrager.C:177
TensorMechanicsPlasticDruckerPrager
Rate-independent non-associative Drucker Prager with hardening/softening.
Definition: TensorMechanicsPlasticDruckerPrager.h:27
defineLegacyParams
defineLegacyParams(TensorMechanicsPlasticDruckerPragerHyperbolic)
TensorMechanicsPlasticDruckerPragerHyperbolic::_use_custom_returnMap
const bool _use_custom_returnMap
whether to use the custom returnMap function
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.h:71
TensorMechanicsPlasticDruckerPragerHyperbolic::_use_custom_cto
const bool _use_custom_cto
Whether to use the custom consistent tangent operator calculation.
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.h:74
TensorMechanicsPlasticDruckerPragerHyperbolic.h
TensorMechanicsPlasticDruckerPrager::dilation
Definition: TensorMechanicsPlasticDruckerPrager.h:49
TensorMechanicsPlasticDruckerPragerHyperbolic::_max_iters
const unsigned _max_iters
max iters for custom return map loop
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.h:77
TensorMechanicsPlasticDruckerPragerHyperbolic::yieldFunction
Real yieldFunction(const RankTwoTensor &stress, Real intnl) const override
The following functions are what you should override when building single-plasticity models.
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.C:62
TensorMechanicsPlasticDruckerPragerHyperbolic
Rate-independent non-associative Drucker Prager with hardening/softening.
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.h:26
TensorMechanicsPlasticDruckerPrager::friction
Definition: TensorMechanicsPlasticDruckerPrager.h:48
TensorMechanicsPlasticModel::_f_tol
const Real _f_tol
Tolerance on yield function.
Definition: TensorMechanicsPlasticModel.h:175
TensorMechanicsPlasticDruckerPragerHyperbolic::useCustomReturnMap
virtual bool useCustomReturnMap() const override
Returns false. You will want to override this in your derived class if you write a custom returnMap f...
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.C:264
TensorMechanicsPlasticDruckerPrager::dyieldFunction_dintnl
virtual Real dyieldFunction_dintnl(const RankTwoTensor &stress, Real intnl) const override
The derivative of yield function with respect to the internal parameter.
Definition: TensorMechanicsPlasticDruckerPrager.C:105
registerMooseObject
registerMooseObject("TensorMechanicsApp", TensorMechanicsPlasticDruckerPragerHyperbolic)
TensorMechanicsPlasticModel::returnMap
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.
Definition: TensorMechanicsPlasticModel.C:221
TensorMechanicsPlasticDruckerPragerHyperbolic::validParams
static InputParameters validParams()
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.C:19
RankTwoTensor
RankTwoTensorTempl< Real > RankTwoTensor
Definition: ACGrGrElasticDrivingForce.h:17
TensorMechanicsPlasticModel::consistentTangentOperator
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.
Definition: TensorMechanicsPlasticModel.C:254
TensorMechanicsPlasticDruckerPragerHyperbolic::_smoother2
const Real _smoother2
smoothing parameter for the cone's tip
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.h:68
TensorMechanicsPlasticDruckerPragerHyperbolic::modelName
virtual std::string modelName() const override
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.C:90
RankFourTensorTempl
Definition: ACGrGrElasticDrivingForce.h:20
TensorMechanicsPlasticDruckerPragerHyperbolic::consistentTangentOperator
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.
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.C:201
TensorMechanicsPlasticDruckerPragerHyperbolic::TensorMechanicsPlasticDruckerPragerHyperbolic
TensorMechanicsPlasticDruckerPragerHyperbolic(const InputParameters &parameters)
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.C:51
TensorMechanicsPlasticDruckerPragerHyperbolic::df_dsig
virtual RankTwoTensor df_dsig(const RankTwoTensor &stress, Real bbb) const override
Function that's used in dyieldFunction_dstress and flowPotential.
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.C:72
RankTwoTensorTempl< Real >
TensorMechanicsPlasticDruckerPrager::validParams
static InputParameters validParams()
Definition: TensorMechanicsPlasticDruckerPrager.C:19
TensorMechanicsPlasticDruckerPragerHyperbolic::dflowPotential_dstress
RankFourTensor dflowPotential_dstress(const RankTwoTensor &stress, Real intnl) const override
The derivative of the flow potential with respect to stress.
Definition: TensorMechanicsPlasticDruckerPragerHyperbolic.C:79
TensorMechanicsPlasticDruckerPrager::dbothAB
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.
Definition: TensorMechanicsPlasticDruckerPrager.C:219
TensorMechanicsPlasticDruckerPrager::bothAB
void bothAB(Real intnl, Real &aaa, Real &bbb) const
Calculates aaa and bbb as a function of the internal parameter intnl.
Definition: TensorMechanicsPlasticDruckerPrager.C:149
RankFourTensor
RankFourTensorTempl< Real > RankFourTensor
Definition: ACGrGrElasticDrivingForce.h:20