www.mooseframework.org
TensorMechanicsPlasticMeanCapTC.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<TensorMechanicsPlasticModel>();
15  params.addRangeCheckedParam<unsigned>("max_iterations",
16  10,
17  "max_iterations>0",
18  "Maximum iterations for custom MeanCapTC return map");
19  params.addParam<bool>(
20  "use_custom_returnMap", true, "Whether to use the custom MeanCapTC returnMap algorithm.");
21  params.addParam<bool>("use_custom_cto",
22  true,
23  "Whether to use the custom consistent tangent operator computations.");
24  params.addRequiredParam<UserObjectName>("tensile_strength",
25  "A TensorMechanicsHardening UserObject that defines "
26  "hardening of the mean-cap tensile strength (it will "
27  "typically be positive). Yield function = trace(stress) "
28  "- tensile_strength for trace(stress)>tensile_strength.");
29  params.addRequiredParam<UserObjectName>(
30  "compressive_strength",
31  "A TensorMechanicsHardening UserObject that defines hardening of the "
32  "mean-cap compressive strength. This should always be less than "
33  "tensile_strength (it will typically be negative). Yield function = "
34  "- (trace(stress) - compressive_strength) for "
35  "trace(stress)<compressive_strength.");
36  params.addClassDescription(
37  "Associative mean-cap tensile and compressive plasticity with hardening/softening");
38 
39  return params;
40 }
41 
43  : TensorMechanicsPlasticModel(parameters),
44  _max_iters(getParam<unsigned>("max_iterations")),
45  _use_custom_returnMap(getParam<bool>("use_custom_returnMap")),
46  _use_custom_cto(getParam<bool>("use_custom_cto")),
47  _strength(getUserObject<TensorMechanicsHardeningModel>("tensile_strength")),
48  _c_strength(getUserObject<TensorMechanicsHardeningModel>("compressive_strength"))
49 {
50  // cannot check the following for all values of the internal parameter, but this will catch most
51  // errors
52  if (_strength.value(0) <= _c_strength.value(0))
53  mooseError("MeanCapTC: tensile strength (which is usually positive) must not be less than "
54  "compressive strength (which is usually negative)");
55 }
56 
57 Real
58 TensorMechanicsPlasticMeanCapTC::yieldFunction(const RankTwoTensor & stress, Real intnl) const
59 {
60  const Real tr = stress.trace();
61  const Real t_str = tensile_strength(intnl);
62  if (tr >= t_str)
63  return tr - t_str;
64 
65  const Real c_str = compressive_strength(intnl);
66  if (tr <= c_str)
67  return -(tr - c_str);
68  // the following is zero at tr = t_str, and at tr = c_str
69  // it also has derivative = 1 at tr = t_str, and derivative = -1 at tr = c_str
70  // it also has second derivative = 0, at these points.
71  // This makes the complete yield function C2 continuous.
72  return (c_str - t_str) / libMesh::pi * std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str));
73 }
74 
75 RankTwoTensor
77  Real intnl) const
78 {
79  return df_dsig(stress, intnl);
80 }
81 
82 Real
84  Real intnl) const
85 {
86  const Real tr = stress.trace();
87  const Real t_str = tensile_strength(intnl);
88  if (tr >= t_str)
89  return -dtensile_strength(intnl);
90 
91  const Real c_str = compressive_strength(intnl);
92  if (tr <= c_str)
93  return dcompressive_strength(intnl);
94 
95  const Real dt = dtensile_strength(intnl);
96  const Real dc = dcompressive_strength(intnl);
97  return (dc - dt) / libMesh::pi * std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) +
98  1.0 / (t_str - c_str) * std::cos(libMesh::pi * (tr - c_str) / (t_str - c_str)) *
99  ((tr - c_str) * dt - (tr - t_str) * dc);
100 }
101 
102 RankTwoTensor
103 TensorMechanicsPlasticMeanCapTC::df_dsig(const RankTwoTensor & stress, Real intnl) const
104 {
105  const Real tr = stress.trace();
106  const Real t_str = tensile_strength(intnl);
107  if (tr >= t_str)
108  return stress.dtrace();
109 
110  const Real c_str = compressive_strength(intnl);
111  if (tr <= c_str)
112  return -stress.dtrace();
113 
114  return -std::cos(libMesh::pi * (tr - c_str) / (t_str - c_str)) * stress.dtrace();
115 }
116 
117 RankTwoTensor
118 TensorMechanicsPlasticMeanCapTC::flowPotential(const RankTwoTensor & stress, Real intnl) const
119 {
120  return df_dsig(stress, intnl);
121 }
122 
123 RankFourTensor
125  Real intnl) const
126 {
127  const Real tr = stress.trace();
128  const Real t_str = tensile_strength(intnl);
129  if (tr >= t_str)
130  return RankFourTensor();
131 
132  const Real c_str = compressive_strength(intnl);
133  if (tr <= c_str)
134  return RankFourTensor();
135 
136  return libMesh::pi / (t_str - c_str) * std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) *
137  stress.dtrace().outerProduct(stress.dtrace());
138 }
139 
140 RankTwoTensor
142  Real intnl) const
143 {
144  const Real tr = stress.trace();
145  const Real t_str = tensile_strength(intnl);
146  if (tr >= t_str)
147  return RankTwoTensor();
148 
149  const Real c_str = compressive_strength(intnl);
150  if (tr <= c_str)
151  return RankTwoTensor();
152 
153  const Real dt = dtensile_strength(intnl);
154  const Real dc = dcompressive_strength(intnl);
155  return std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) * stress.dtrace() * libMesh::pi /
156  Utility::pow<2>(t_str - c_str) * ((tr - t_str) * dc - (tr - c_str) * dt);
157 }
158 
159 Real
160 TensorMechanicsPlasticMeanCapTC::hardPotential(const RankTwoTensor & stress, Real intnl) const
161 {
162  // This is the key for this whole class!
163  const Real tr = stress.trace();
164  const Real t_str = tensile_strength(intnl);
165 
166  if (tr >= t_str)
167  return -1.0; // this will serve to *increase* the internal parameter (so internal parameter will
168  // be a measure of volumetric strain)
169 
170  const Real c_str = compressive_strength(intnl);
171  if (tr <= c_str)
172  return 1.0; // this will serve to *decrease* the internal parameter (so internal parameter will
173  // be a measure of volumetric strain)
174 
175  return std::cos(libMesh::pi * (tr - c_str) /
176  (t_str - c_str)); // this interpolates C2 smoothly between 1 and -1
177 }
178 
179 RankTwoTensor
181  Real intnl) const
182 {
183  const Real tr = stress.trace();
184  const Real t_str = tensile_strength(intnl);
185  if (tr >= t_str)
186  return RankTwoTensor();
187 
188  const Real c_str = compressive_strength(intnl);
189  if (tr <= c_str)
190  return RankTwoTensor();
191 
192  return -std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) * libMesh::pi / (t_str - c_str) *
193  stress.dtrace();
194 }
195 
196 Real
198  Real intnl) const
199 {
200  const Real tr = stress.trace();
201  const Real t_str = tensile_strength(intnl);
202  if (tr >= t_str)
203  return 0.0;
204 
205  const Real c_str = compressive_strength(intnl);
206  if (tr <= c_str)
207  return 0.0;
208 
209  const Real dt = dtensile_strength(intnl);
210  const Real dc = dcompressive_strength(intnl);
211  return -std::sin(libMesh::pi * (tr - c_str) / (t_str - c_str)) * libMesh::pi /
212  Utility::pow<2>(t_str - c_str) * ((tr - t_str) * dc - (tr - c_str) * dt);
213 }
214 
215 Real
217 {
218  return _strength.value(internal_param);
219 }
220 
221 Real
223 {
224  return _strength.derivative(internal_param);
225 }
226 
227 Real
229 {
230  return _c_strength.value(internal_param);
231 }
232 
233 Real
235 {
236  return _c_strength.derivative(internal_param);
237 }
238 
239 void
241  const RankTwoTensor & stress,
242  Real intnl,
243  const RankFourTensor & Eijkl,
244  std::vector<bool> & act,
245  RankTwoTensor & returned_stress) const
246 {
247  act.assign(1, false);
248 
249  if (f[0] <= _f_tol)
250  {
251  returned_stress = stress;
252  return;
253  }
254 
255  const Real tr = stress.trace();
256  const Real t_str = tensile_strength(intnl);
257  Real str;
258  Real dirn;
259  if (tr >= t_str)
260  {
261  str = t_str;
262  dirn = 1.0;
263  }
264  else
265  {
266  str = compressive_strength(intnl);
267  dirn = -1.0;
268  }
269 
270  RankTwoTensor n; // flow direction
271  for (unsigned i = 0; i < 3; ++i)
272  for (unsigned j = 0; j < 3; ++j)
273  for (unsigned k = 0; k < 3; ++k)
274  n(i, j) += dirn * Eijkl(i, j, k, k);
275 
276  // returned_stress = stress - gamma*n
277  // and taking the trace of this and using
278  // Tr(returned_stress) = str, gives
279  // gamma = (Tr(stress) - str)/Tr(n)
280  Real gamma = (stress.trace() - str) / n.trace();
281 
282  for (unsigned i = 0; i < 3; ++i)
283  for (unsigned j = 0; j < 3; ++j)
284  returned_stress(i, j) = stress(i, j) - gamma * n(i, j);
285 
286  act[0] = true;
287 }
288 
289 bool
290 TensorMechanicsPlasticMeanCapTC::returnMap(const RankTwoTensor & trial_stress,
291  Real intnl_old,
292  const RankFourTensor & E_ijkl,
293  Real ep_plastic_tolerance,
294  RankTwoTensor & returned_stress,
295  Real & returned_intnl,
296  std::vector<Real> & dpm,
297  RankTwoTensor & delta_dp,
298  std::vector<Real> & yf,
299  bool & trial_stress_inadmissible) const
300 {
301  if (!(_use_custom_returnMap))
302  return TensorMechanicsPlasticModel::returnMap(trial_stress,
303  intnl_old,
304  E_ijkl,
305  ep_plastic_tolerance,
306  returned_stress,
307  returned_intnl,
308  dpm,
309  delta_dp,
310  yf,
311  trial_stress_inadmissible);
312 
313  yf.resize(1);
314 
315  Real yf_orig = yieldFunction(trial_stress, intnl_old);
316 
317  yf[0] = yf_orig;
318 
319  if (yf_orig < _f_tol)
320  {
321  // the trial_stress is admissible
322  trial_stress_inadmissible = false;
323  return true;
324  }
325 
326  trial_stress_inadmissible = true;
327 
328  // In the following we want to solve
329  // trial_stress - stress = E_ijkl * dpm * r ...... (1)
330  // and either
331  // stress.trace() = tensile_strength(intnl) ...... (2a)
332  // intnl = intnl_old + dpm ...... (3a)
333  // or
334  // stress.trace() = compressive_strength(intnl) ... (2b)
335  // intnl = intnl_old - dpm ...... (3b)
336  // The former (2a and 3a) are chosen if
337  // trial_stress.trace() > tensile_strength(intnl_old)
338  // while the latter (2b and 3b) are chosen if
339  // trial_stress.trace() < compressive_strength(intnl_old)
340  // The variables we want to solve for are stress, dpm
341  // and intnl. We do this using a Newton approach, starting
342  // with stress=trial_stress and intnl=intnl_old and dpm=0
343  const bool tensile_failure = (trial_stress.trace() >= tensile_strength(intnl_old));
344  const Real dirn = (tensile_failure ? 1.0 : -1.0);
345 
346  RankTwoTensor n; // flow direction, which is E_ijkl * r
347  for (unsigned i = 0; i < 3; ++i)
348  for (unsigned j = 0; j < 3; ++j)
349  for (unsigned k = 0; k < 3; ++k)
350  n(i, j) += dirn * E_ijkl(i, j, k, k);
351  const Real n_trace = n.trace();
352 
353  // Perform a Newton-Raphson to find dpm when
354  // residual = trial_stress.trace() - tensile_strength(intnl) - dpm * n.trace() [for
355  // tensile_failure=true]
356  // or
357  // residual = trial_stress.trace() - compressive_strength(intnl) - dpm * n.trace() [for
358  // tensile_failure=false]
359  Real trial_trace = trial_stress.trace();
360  Real residual;
361  Real jac;
362  dpm[0] = 0;
363  unsigned int iter = 0;
364  do
365  {
366  if (tensile_failure)
367  {
368  residual = trial_trace - tensile_strength(intnl_old + dpm[0]) - dpm[0] * n_trace;
369  jac = -dtensile_strength(intnl_old + dpm[0]) - n_trace;
370  }
371  else
372  {
373  residual = trial_trace - compressive_strength(intnl_old - dpm[0]) - dpm[0] * n_trace;
374  jac = -dcompressive_strength(intnl_old - dpm[0]) - n_trace;
375  }
376  dpm[0] += -residual / jac;
377  if (iter > _max_iters) // not converging
378  return false;
379  iter++;
380  } while (residual * residual > _f_tol * _f_tol);
381 
382  // set the returned values
383  yf[0] = 0;
384  returned_intnl = intnl_old + dirn * dpm[0];
385  returned_stress = trial_stress - dpm[0] * n;
386  delta_dp = dpm[0] * dirn * returned_stress.dtrace();
387 
388  return true;
389 }
390 
391 RankFourTensor
393  const RankTwoTensor & trial_stress,
394  Real intnl_old,
395  const RankTwoTensor & stress,
396  Real intnl,
397  const RankFourTensor & E_ijkl,
398  const std::vector<Real> & cumulative_pm) const
399 {
400  if (!_use_custom_cto)
402  trial_stress, intnl_old, stress, intnl, E_ijkl, cumulative_pm);
403 
404  Real df_dq;
405  Real alpha;
406  if (trial_stress.trace() >= tensile_strength(intnl_old))
407  {
408  df_dq = -dtensile_strength(intnl);
409  alpha = 1.0;
410  }
411  else
412  {
413  df_dq = dcompressive_strength(intnl);
414  alpha = -1.0;
415  }
416 
417  RankTwoTensor elas;
418  for (unsigned int i = 0; i < 3; ++i)
419  for (unsigned int j = 0; j < 3; ++j)
420  for (unsigned int k = 0; k < 3; ++k)
421  elas(i, j) += E_ijkl(i, j, k, k);
422 
423  const Real hw = -df_dq + alpha * elas.trace();
424 
425  return E_ijkl - alpha / hw * elas.outerProduct(elas);
426 }
427 
428 bool
430 {
431  return _use_custom_returnMap;
432 }
433 
434 bool
436 {
437  return _use_custom_cto;
438 }
439 
440 std::string
442 {
443  return "MeanCapTC";
444 }
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.
RankTwoTensor flowPotential(const RankTwoTensor &stress, Real intnl) const override
The flow potential.
Real yieldFunction(const RankTwoTensor &stress, Real intnl) const override
The following functions are what you should override when building single-plasticity models...
virtual Real dcompressive_strength(const Real internal_param) const
d(compressive strength)/d(internal_param) as a function of residual value, rate, and internal_param ...
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.
RankTwoTensor dflowPotential_dintnl(const RankTwoTensor &stress, Real intnl) const override
The derivative of the flow potential with respect to the internal parameter.
virtual RankTwoTensor dhardPotential_dstress(const RankTwoTensor &stress, Real intnl) const override
The derivative of the hardening potential with respect to stress.
const TensorMechanicsHardeningModel & _c_strength
the compressive strength
TensorMechanicsPlasticMeanCapTC(const InputParameters &parameters)
virtual Real tensile_strength(const Real internal_param) const
tensile strength as a function of residual value, rate, and internal_param
virtual Real derivative(Real intnl) const
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.
RankTwoTensor df_dsig(const RankTwoTensor &stress, Real intnl) const
Derivative of the yield function with respect to stress.
const bool _use_custom_returnMap
Whether to use the custom return-map algorithm.
virtual Real dtensile_strength(const Real internal_param) const
d(tensile strength)/d(internal_param) as a function of residual value, rate, and internal_param ...
const unsigned _max_iters
max iters for custom return map loop
Real hardPotential(const RankTwoTensor &stress, Real intnl) const override
The hardening potential.
virtual void activeConstraints(const std::vector< Real > &f, const RankTwoTensor &stress, Real intnl, const RankFourTensor &Eijkl, std::vector< bool > &act, RankTwoTensor &returned_stress) const override
The active yield surfaces, given a vector of yield functions.
const bool _use_custom_cto
Whether to use the custom consistent tangent operator algorithm.
RankTwoTensor dyieldFunction_dstress(const RankTwoTensor &stress, Real intnl) const override
The derivative of yield function with respect to stress.
RankFourTensor dflowPotential_dstress(const RankTwoTensor &stress, Real intnl) const override
The derivative of the flow potential with respect to stress.
InputParameters validParams< TensorMechanicsPlasticModel >()
virtual Real value(Real intnl) const
virtual bool useCustomReturnMap() const override
Returns false. You will want to override this in your derived class if you write a custom returnMap f...
virtual Real dhardPotential_dintnl(const RankTwoTensor &stress, Real intnl) const override
The derivative of the hardening potential with respect to the internal parameter. ...
const Real _f_tol
Tolerance on yield function.
Real dyieldFunction_dintnl(const RankTwoTensor &stress, Real intnl) const override
The derivative of yield function with respect to the internal parameter.
virtual Real compressive_strength(const Real internal_param) const
compressive strength as a function of residual value, rate, and internal_param
Plastic Model base class The virtual functions written below must be over-ridden in derived classes t...
InputParameters validParams< TensorMechanicsPlasticMeanCapTC >()
virtual bool useCustomCTO() const override
Returns false. You will want to override this in your derived class if you write a custom consistent ...
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 std::string modelName() const override
const TensorMechanicsHardeningModel & _strength
the tensile strength