www.mooseframework.org
TensorMechanicsPlasticTensile.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.addRequiredParam<UserObjectName>(
16  "tensile_strength",
17  "A TensorMechanicsHardening UserObject that defines hardening of the tensile strength");
18  params.addRangeCheckedParam<Real>(
19  "tensile_edge_smoother",
20  25.0,
21  "tensile_edge_smoother>=0 & tensile_edge_smoother<=30",
22  "Smoothing parameter: the edges of the cone are smoothed by the given amount.");
23  MooseEnum tip_scheme("hyperbolic cap", "hyperbolic");
24  params.addParam<MooseEnum>(
25  "tip_scheme", tip_scheme, "Scheme by which the pyramid's tip will be smoothed.");
26  params.addRequiredRangeCheckedParam<Real>("tensile_tip_smoother",
27  "tensile_tip_smoother>=0",
28  "For the 'hyperbolic' tip_scheme, the pyramid vertex "
29  "will be smoothed by the given amount. For the 'cap' "
30  "tip_scheme, additional smoothing will occur. Typical "
31  "value is 0.1*tensile_strength");
32  params.addParam<Real>(
33  "cap_start",
34  0.0,
35  "For the 'cap' tip_scheme, smoothing is performed in the stress_mean > cap_start region");
36  params.addRangeCheckedParam<Real>("cap_rate",
37  0.0,
38  "cap_rate>=0",
39  "For the 'cap' tip_scheme, this controls how quickly the cap "
40  "degenerates to a hemisphere: small values mean a slow "
41  "degeneration to a hemisphere (and zero means the 'cap' will "
42  "be totally inactive). Typical value is 1/tensile_strength");
43  params.addParam<Real>("tensile_lode_cutoff",
44  "If the second invariant of stress is less than "
45  "this amount, the Lode angle is assumed to be zero. "
46  " This is to gaurd against precision-loss problems, "
47  "and this parameter should be set small. Default = "
48  "0.00001*((yield_Function_tolerance)^2)");
49  params.addClassDescription(
50  "Associative tensile plasticity with hardening/softening, and tensile_strength = 1");
51 
52  return params;
53 }
54 
56  : TensorMechanicsPlasticModel(parameters),
57  _strength(getUserObject<TensorMechanicsHardeningModel>("tensile_strength")),
58  _tip_scheme(getParam<MooseEnum>("tip_scheme")),
59  _small_smoother2(Utility::pow<2>(getParam<Real>("tensile_tip_smoother"))),
60  _cap_start(getParam<Real>("cap_start")),
61  _cap_rate(getParam<Real>("cap_rate")),
62  _tt(getParam<Real>("tensile_edge_smoother") * libMesh::pi / 180.0),
63  _sin3tt(std::sin(3.0 * _tt)),
64  _lode_cutoff(parameters.isParamValid("tensile_lode_cutoff")
65  ? getParam<Real>("tensile_lode_cutoff")
66  : 1.0e-5 * Utility::pow<2>(_f_tol))
67 
68 {
69  if (_lode_cutoff < 0)
70  mooseError("tensile_lode_cutoff must not be negative");
71  _ccc = (-std::cos(3.0 * _tt) * (std::cos(_tt) - std::sin(_tt) / std::sqrt(3.0)) -
72  3.0 * _sin3tt * (std::sin(_tt) + std::cos(_tt) / std::sqrt(3.0))) /
73  (18.0 * Utility::pow<3>(std::cos(3.0 * _tt)));
74  _bbb = (std::sin(6.0 * _tt) * (std::cos(_tt) - std::sin(_tt) / std::sqrt(3.0)) -
75  6.0 * std::cos(6.0 * _tt) * (std::sin(_tt) + std::cos(_tt) / std::sqrt(3.0))) /
76  (18.0 * Utility::pow<3>(std::cos(3.0 * _tt)));
77  _aaa = -std::sin(_tt) / std::sqrt(3.0) - _bbb * _sin3tt - _ccc * Utility::pow<2>(_sin3tt) +
78  std::cos(_tt);
79 }
80 
81 Real
82 TensorMechanicsPlasticTensile::yieldFunction(const RankTwoTensor & stress, Real intnl) const
83 {
84  Real mean_stress = stress.trace() / 3.0;
85  Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
86  if (sin3Lode <= _sin3tt)
87  {
88  // the non-edge-smoothed version
89  std::vector<Real> eigvals;
90  stress.symmetricEigenvalues(eigvals);
91  return mean_stress + std::sqrt(smooth(stress) + Utility::pow<2>(eigvals[2] - mean_stress)) -
92  tensile_strength(intnl);
93  }
94  else
95  {
96  // the edge-smoothed version
97  Real kk = _aaa + _bbb * sin3Lode + _ccc * Utility::pow<2>(sin3Lode);
98  Real sibar2 = stress.secondInvariant();
99  return mean_stress + std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk)) -
100  tensile_strength(intnl);
101  }
102 }
103 
104 RankTwoTensor
106  Real /*intnl*/) const
107 {
108  Real mean_stress = stress.trace() / 3.0;
109  RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
110  Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
111  if (sin3Lode <= _sin3tt)
112  {
113  // the non-edge-smoothed version
114  std::vector<Real> eigvals;
115  std::vector<RankTwoTensor> deigvals;
116  stress.dsymmetricEigenvalues(eigvals, deigvals);
117  Real denom = std::sqrt(smooth(stress) + Utility::pow<2>(eigvals[2] - mean_stress));
118  return dmean_stress +
119  (0.5 * dsmooth(stress) * dmean_stress +
120  (eigvals[2] - mean_stress) * (deigvals[2] - dmean_stress)) /
121  denom;
122  }
123  else
124  {
125  // the edge-smoothed version
126  Real kk = _aaa + _bbb * sin3Lode + _ccc * Utility::pow<2>(sin3Lode);
127  RankTwoTensor dkk = (_bbb + 2.0 * _ccc * sin3Lode) * stress.dsin3Lode(_lode_cutoff);
128  Real sibar2 = stress.secondInvariant();
129  RankTwoTensor dsibar2 = stress.dsecondInvariant();
130  Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
131  return dmean_stress +
132  (0.5 * dsmooth(stress) * dmean_stress + 0.5 * dsibar2 * Utility::pow<2>(kk) +
133  sibar2 * kk * dkk) /
134  denom;
135  }
136 }
137 
138 Real
140  Real intnl) const
141 {
142  return -dtensile_strength(intnl);
143 }
144 
145 RankTwoTensor
146 TensorMechanicsPlasticTensile::flowPotential(const RankTwoTensor & stress, Real intnl) const
147 {
148  // This plasticity is associative so
149  return dyieldFunction_dstress(stress, intnl);
150 }
151 
152 RankFourTensor
154  Real /*intnl*/) const
155 {
156  Real mean_stress = stress.trace() / 3.0;
157  RankTwoTensor dmean_stress = stress.dtrace() / 3.0;
158  Real sin3Lode = stress.sin3Lode(_lode_cutoff, 0);
159  if (sin3Lode <= _sin3tt)
160  {
161  // the non-edge-smoothed version
162  std::vector<Real> eigvals;
163  std::vector<RankTwoTensor> deigvals;
164  std::vector<RankFourTensor> d2eigvals;
165  stress.dsymmetricEigenvalues(eigvals, deigvals);
166  stress.d2symmetricEigenvalues(d2eigvals);
167 
168  Real denom = std::sqrt(smooth(stress) + Utility::pow<2>(eigvals[2] - mean_stress));
169  Real denom3 = Utility::pow<3>(denom);
170  RankTwoTensor numer_part = deigvals[2] - dmean_stress;
171  RankTwoTensor numer_full =
172  0.5 * dsmooth(stress) * dmean_stress + (eigvals[2] - mean_stress) * numer_part;
173  Real d2smooth_over_denom = d2smooth(stress) / denom;
174 
175  RankFourTensor dr_dstress = (eigvals[2] - mean_stress) * d2eigvals[2] / denom;
176  for (unsigned i = 0; i < 3; ++i)
177  for (unsigned j = 0; j < 3; ++j)
178  for (unsigned k = 0; k < 3; ++k)
179  for (unsigned l = 0; l < 3; ++l)
180  {
181  dr_dstress(i, j, k, l) +=
182  0.5 * d2smooth_over_denom * dmean_stress(i, j) * dmean_stress(k, l);
183  dr_dstress(i, j, k, l) += numer_part(i, j) * numer_part(k, l) / denom;
184  dr_dstress(i, j, k, l) -= numer_full(i, j) * numer_full(k, l) / denom3;
185  }
186  return dr_dstress;
187  }
188  else
189  {
190  // the edge-smoothed version
191  RankTwoTensor dsin3Lode = stress.dsin3Lode(_lode_cutoff);
192  Real kk = _aaa + _bbb * sin3Lode + _ccc * Utility::pow<2>(sin3Lode);
193  RankTwoTensor dkk = (_bbb + 2.0 * _ccc * sin3Lode) * dsin3Lode;
194  RankFourTensor d2kk = (_bbb + 2.0 * _ccc * sin3Lode) * stress.d2sin3Lode(_lode_cutoff);
195  for (unsigned i = 0; i < 3; ++i)
196  for (unsigned j = 0; j < 3; ++j)
197  for (unsigned k = 0; k < 3; ++k)
198  for (unsigned l = 0; l < 3; ++l)
199  d2kk(i, j, k, l) += 2.0 * _ccc * dsin3Lode(i, j) * dsin3Lode(k, l);
200 
201  Real sibar2 = stress.secondInvariant();
202  RankTwoTensor dsibar2 = stress.dsecondInvariant();
203  RankFourTensor d2sibar2 = stress.d2secondInvariant();
204 
205  Real denom = std::sqrt(smooth(stress) + sibar2 * Utility::pow<2>(kk));
206  Real denom3 = Utility::pow<3>(denom);
207  Real d2smooth_over_denom = d2smooth(stress) / denom;
208  RankTwoTensor numer_full =
209  0.5 * dsmooth(stress) * dmean_stress + 0.5 * dsibar2 * kk * kk + sibar2 * kk * dkk;
210 
211  RankFourTensor dr_dstress = (0.5 * d2sibar2 * Utility::pow<2>(kk) + sibar2 * kk * d2kk) / denom;
212  for (unsigned i = 0; i < 3; ++i)
213  for (unsigned j = 0; j < 3; ++j)
214  for (unsigned k = 0; k < 3; ++k)
215  for (unsigned l = 0; l < 3; ++l)
216  {
217  dr_dstress(i, j, k, l) +=
218  0.5 * d2smooth_over_denom * dmean_stress(i, j) * dmean_stress(k, l);
219  dr_dstress(i, j, k, l) +=
220  (dsibar2(i, j) * dkk(k, l) * kk + dkk(i, j) * dsibar2(k, l) * kk +
221  sibar2 * dkk(i, j) * dkk(k, l)) /
222  denom;
223  dr_dstress(i, j, k, l) -= numer_full(i, j) * numer_full(k, l) / denom3;
224  }
225  return dr_dstress;
226  }
227 }
228 
229 RankTwoTensor
231  Real /*intnl*/) const
232 {
233  return RankTwoTensor();
234 }
235 
236 Real
237 TensorMechanicsPlasticTensile::tensile_strength(const Real internal_param) const
238 {
239  return _strength.value(internal_param);
240 }
241 
242 Real
243 TensorMechanicsPlasticTensile::dtensile_strength(const Real internal_param) const
244 {
245  return _strength.derivative(internal_param);
246 }
247 
248 Real
249 TensorMechanicsPlasticTensile::smooth(const RankTwoTensor & stress) const
250 {
251  Real smoother2 = _small_smoother2;
252  if (_tip_scheme == "cap")
253  {
254  Real x = stress.trace() / 3.0 - _cap_start;
255  Real p = 0;
256  if (x > 0)
257  p = x * (1 - std::exp(-_cap_rate * x));
258  smoother2 += Utility::pow<2>(p);
259  }
260  return smoother2;
261 }
262 
263 Real
264 TensorMechanicsPlasticTensile::dsmooth(const RankTwoTensor & stress) const
265 {
266  Real dsmoother2 = 0;
267  if (_tip_scheme == "cap")
268  {
269  Real x = stress.trace() / 3.0 - _cap_start;
270  Real p = 0;
271  Real dp_dx = 0;
272  if (x > 0)
273  {
274  p = x * (1 - std::exp(-_cap_rate * x));
275  dp_dx = (1 - std::exp(-_cap_rate * x)) + x * _cap_rate * std::exp(-_cap_rate * x);
276  }
277  dsmoother2 += 2.0 * p * dp_dx;
278  }
279  return dsmoother2;
280 }
281 
282 Real
283 TensorMechanicsPlasticTensile::d2smooth(const RankTwoTensor & stress) const
284 {
285  Real d2smoother2 = 0;
286  if (_tip_scheme == "cap")
287  {
288  Real x = stress.trace() / 3.0 - _cap_start;
289  Real p = 0;
290  Real dp_dx = 0;
291  Real d2p_dx2 = 0;
292  if (x > 0)
293  {
294  p = x * (1 - std::exp(-_cap_rate * x));
295  dp_dx = (1 - std::exp(-_cap_rate * x)) + x * _cap_rate * std::exp(-_cap_rate * x);
296  d2p_dx2 = 2.0 * _cap_rate * std::exp(-_cap_rate * x) -
297  x * Utility::pow<2>(_cap_rate) * std::exp(-_cap_rate * x);
298  }
299  d2smoother2 += 2.0 * Utility::pow<2>(dp_dx) + 2.0 * p * d2p_dx2;
300  }
301  return d2smoother2;
302 }
303 
304 std::string
306 {
307  return "Tensile";
308 }
virtual Real d2smooth(const RankTwoTensor &stress) const
returns the d^2a/dstress_mean^2 - see doco for _tip_scheme
Real dyieldFunction_dintnl(const RankTwoTensor &stress, Real intnl) const override
The derivative of yield function with respect to the internal parameter.
Real _bbb
Abbo et al&#39;s B parameter.
virtual Real dsmooth(const RankTwoTensor &stress) const
returns the da/dstress_mean - see doco for _tip_scheme
TensorMechanicsPlasticTensile(const InputParameters &parameters)
Real _sin3tt
sin(3*_tt) - useful for making comparisons with Lode angle
RankTwoTensor dyieldFunction_dstress(const RankTwoTensor &stress, Real intnl) const override
The derivative of yield function with respect to stress.
virtual Real derivative(Real intnl) const
RankTwoTensor dflowPotential_dintnl(const RankTwoTensor &stress, Real intnl) const override
The derivative of the flow potential with respect to the internal parameter.
Real _ccc
Abbo et al&#39;s C parameter.
virtual Real tensile_strength(const Real internal_param) const
tensile strength as a function of residual value, rate, and internal_param
virtual std::string modelName() const override
RankTwoTensor flowPotential(const RankTwoTensor &stress, Real intnl) const override
The flow potential.
InputParameters validParams< TensorMechanicsPlasticModel >()
Real _tt
edge smoothing parameter, in radians
Real _small_smoother2
Square of tip smoothing parameter to smooth the cone at mean_stress = T.
Real yieldFunction(const RankTwoTensor &stress, Real intnl) const override
The following functions are what you should override when building single-plasticity models...
virtual Real value(Real intnl) const
Real _cap_start
smoothing parameter dictating when the &#39;cap&#39; will start - see doco for _tip_scheme ...
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
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 ...
Real _aaa
Abbo et al&#39;s A parameter.
Plastic Model base class The virtual functions written below must be over-ridden in derived classes t...
InputParameters validParams< TensorMechanicsPlasticTensile >()
Real _cap_rate
dictates how quickly the &#39;cap&#39; degenerates to a hemisphere - see doco for _tip_scheme ...
MooseEnum _tip_scheme
The yield function is modified to f = s_m + sqrt(a + s_bar^2 K^2) - tensile_strength where "a" depend...
RankFourTensor dflowPotential_dstress(const RankTwoTensor &stress, Real intnl) const override
The derivative of the flow potential with respect to stress.
virtual Real smooth(const RankTwoTensor &stress) const
returns the &#39;a&#39; parameter - see doco for _tip_scheme
const TensorMechanicsHardeningModel & _strength
Real _lode_cutoff
if secondInvariant < _lode_cutoff then set Lode angle to zero. This is to guard against precision-los...