www.mooseframework.org
TensorMechanicsPlasticWeakPlaneShear.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  "cohesion",
17  "A TensorMechanicsHardening UserObject that defines hardening of the cohesion. "
18  "Physically the cohesion should not be negative.");
19  params.addRequiredParam<UserObjectName>("tan_friction_angle",
20  "A TensorMechanicsHardening UserObject that defines "
21  "hardening of tan(friction angle). Physically the "
22  "friction angle should be between 0 and 90deg.");
23  params.addRequiredParam<UserObjectName>(
24  "tan_dilation_angle",
25  "A TensorMechanicsHardening UserObject that defines hardening of the "
26  "tan(dilation angle). Usually the dilation angle is not greater than "
27  "the friction angle, and it is between 0 and 90deg.");
28  MooseEnum tip_scheme("hyperbolic cap", "hyperbolic");
29  params.addParam<MooseEnum>(
30  "tip_scheme", tip_scheme, "Scheme by which the cone's tip will be smoothed.");
31  params.addRequiredRangeCheckedParam<Real>(
32  "smoother",
33  "smoother>=0",
34  "For the 'hyperbolic' tip_scheme, the cone vertex at shear-stress "
35  "= 0 will be smoothed by the given amount. For the 'cap' "
36  "tip_scheme, additional smoothing will occur. Typical value is "
37  "0.1*cohesion");
38  params.addParam<Real>(
39  "cap_start",
40  0.0,
41  "For the 'cap' tip_scheme, smoothing is performed in the stress_zz > cap_start region");
42  params.addRangeCheckedParam<Real>("cap_rate",
43  0.0,
44  "cap_rate>=0",
45  "For the 'cap' tip_scheme, this controls how quickly the cap "
46  "degenerates to a hemisphere: small values mean a slow "
47  "degeneration to a hemisphere (and zero means the 'cap' will "
48  "be totally inactive). Typical value is 1/cohesion");
49  params.addClassDescription("Non-associative finite-strain weak-plane shear perfect plasticity. "
50  "Here cohesion = 1, tan(phi) = 1 = tan(psi)");
51 
52  return params;
53 }
54 
56  const InputParameters & parameters)
57  : TensorMechanicsPlasticModel(parameters),
58  _cohesion(getUserObject<TensorMechanicsHardeningModel>("cohesion")),
59  _tan_phi(getUserObject<TensorMechanicsHardeningModel>("tan_friction_angle")),
60  _tan_psi(getUserObject<TensorMechanicsHardeningModel>("tan_dilation_angle")),
61  _tip_scheme(getParam<MooseEnum>("tip_scheme")),
62  _small_smoother2(Utility::pow<2>(getParam<Real>("smoother"))),
63  _cap_start(getParam<Real>("cap_start")),
64  _cap_rate(getParam<Real>("cap_rate"))
65 {
66  // With arbitary UserObjects, it is impossible to check everything, and
67  // I think this is the best I can do
68  if (tan_phi(0) < 0 || tan_psi(0) < 0)
69  mooseError("Weak-Plane-Shear friction and dilation angles must lie in [0, Pi/2]");
70  if (tan_phi(0) < tan_psi(0))
71  mooseError(
72  "Weak-Plane-Shear friction angle must not be less than Weak-Plane-Shear dilation angle");
73  if (cohesion(0) < 0)
74  mooseError("Weak-Plane-Shear cohesion must not be negative");
75 }
76 
77 Real
78 TensorMechanicsPlasticWeakPlaneShear::yieldFunction(const RankTwoTensor & stress, Real intnl) const
79 {
80  // note that i explicitly symmeterise in preparation for Cosserat
81  return std::sqrt(Utility::pow<2>((stress(0, 2) + stress(2, 0)) / 2.0) +
82  Utility::pow<2>((stress(1, 2) + stress(2, 1)) / 2.0) + smooth(stress)) +
83  stress(2, 2) * tan_phi(intnl) - cohesion(intnl);
84 }
85 
86 RankTwoTensor
87 TensorMechanicsPlasticWeakPlaneShear::df_dsig(const RankTwoTensor & stress,
88  Real _tan_phi_or_psi) const
89 {
90  RankTwoTensor deriv; // the constructor zeroes this
91 
92  Real tau = std::sqrt(Utility::pow<2>((stress(0, 2) + stress(2, 0)) / 2.0) +
93  Utility::pow<2>((stress(1, 2) + stress(2, 1)) / 2.0) + smooth(stress));
94  // note that i explicitly symmeterise in preparation for Cosserat
95  if (tau == 0.0)
96  {
97  // the derivative is not defined here, but i want to set it nonzero
98  // because otherwise the return direction might be too crazy
99  deriv(0, 2) = deriv(2, 0) = 0.5;
100  deriv(1, 2) = deriv(2, 1) = 0.5;
101  }
102  else
103  {
104  deriv(0, 2) = deriv(2, 0) = 0.25 * (stress(0, 2) + stress(2, 0)) / tau;
105  deriv(1, 2) = deriv(2, 1) = 0.25 * (stress(1, 2) + stress(2, 1)) / tau;
106  deriv(2, 2) = 0.5 * dsmooth(stress) / tau;
107  }
108  deriv(2, 2) += _tan_phi_or_psi;
109  return deriv;
110 }
111 
112 RankTwoTensor
114  Real intnl) const
115 {
116  return df_dsig(stress, tan_phi(intnl));
117 }
118 
119 Real
121  Real intnl) const
122 {
123  return stress(2, 2) * dtan_phi(intnl) - dcohesion(intnl);
124 }
125 
126 RankTwoTensor
127 TensorMechanicsPlasticWeakPlaneShear::flowPotential(const RankTwoTensor & stress, Real intnl) const
128 {
129  return df_dsig(stress, tan_psi(intnl));
130 }
131 
132 RankFourTensor
134  Real /*intnl*/) const
135 {
136  RankFourTensor dr_dstress;
137  Real tau = std::sqrt(Utility::pow<2>((stress(0, 2) + stress(2, 0)) / 2.0) +
138  Utility::pow<2>((stress(1, 2) + stress(2, 1)) / 2.0) + smooth(stress));
139  if (tau == 0.0)
140  return dr_dstress;
141 
142  // note that i explicitly symmeterise
143  RankTwoTensor dtau;
144  dtau(0, 2) = dtau(2, 0) = 0.25 * (stress(0, 2) + stress(2, 0)) / tau;
145  dtau(1, 2) = dtau(2, 1) = 0.25 * (stress(1, 2) + stress(2, 1)) / tau;
146  dtau(2, 2) = 0.5 * dsmooth(stress) / tau;
147 
148  for (unsigned i = 0; i < 3; ++i)
149  for (unsigned j = 0; j < 3; ++j)
150  for (unsigned k = 0; k < 3; ++k)
151  for (unsigned l = 0; l < 3; ++l)
152  dr_dstress(i, j, k, l) = -dtau(i, j) * dtau(k, l) / tau;
153 
154  // note that i explicitly symmeterise
155  dr_dstress(0, 2, 0, 2) += 0.25 / tau;
156  dr_dstress(0, 2, 2, 0) += 0.25 / tau;
157  dr_dstress(2, 0, 0, 2) += 0.25 / tau;
158  dr_dstress(2, 0, 2, 0) += 0.25 / tau;
159  dr_dstress(1, 2, 1, 2) += 0.25 / tau;
160  dr_dstress(1, 2, 2, 1) += 0.25 / tau;
161  dr_dstress(2, 1, 1, 2) += 0.25 / tau;
162  dr_dstress(2, 1, 2, 1) += 0.25 / tau;
163  dr_dstress(2, 2, 2, 2) += 0.5 * d2smooth(stress) / tau;
164 
165  return dr_dstress;
166 }
167 
168 RankTwoTensor
170  Real intnl) const
171 {
172  RankTwoTensor dr_dintnl;
173  dr_dintnl(2, 2) = dtan_psi(intnl);
174  return dr_dintnl;
175 }
176 
177 Real
178 TensorMechanicsPlasticWeakPlaneShear::cohesion(const Real internal_param) const
179 {
180  return _cohesion.value(internal_param);
181 }
182 
183 Real
184 TensorMechanicsPlasticWeakPlaneShear::dcohesion(const Real internal_param) const
185 {
186  return _cohesion.derivative(internal_param);
187 }
188 
189 Real
190 TensorMechanicsPlasticWeakPlaneShear::tan_phi(const Real internal_param) const
191 {
192  return _tan_phi.value(internal_param);
193 }
194 
195 Real
196 TensorMechanicsPlasticWeakPlaneShear::dtan_phi(const Real internal_param) const
197 {
198  return _tan_phi.derivative(internal_param);
199 }
200 
201 Real
202 TensorMechanicsPlasticWeakPlaneShear::tan_psi(const Real internal_param) const
203 {
204  return _tan_psi.value(internal_param);
205 }
206 
207 Real
208 TensorMechanicsPlasticWeakPlaneShear::dtan_psi(const Real internal_param) const
209 {
210  return _tan_psi.derivative(internal_param);
211 }
212 
213 Real
214 TensorMechanicsPlasticWeakPlaneShear::smooth(const RankTwoTensor & stress) const
215 {
216  Real smoother2 = _small_smoother2;
217  if (_tip_scheme == "cap")
218  {
219  Real x = stress(2, 2) - _cap_start;
220  Real p = 0;
221  if (x > 0)
222  p = x * (1 - std::exp(-_cap_rate * x));
223  smoother2 += Utility::pow<2>(p);
224  }
225  return smoother2;
226 }
227 
228 Real
229 TensorMechanicsPlasticWeakPlaneShear::dsmooth(const RankTwoTensor & stress) const
230 {
231  Real dsmoother2 = 0;
232  if (_tip_scheme == "cap")
233  {
234  Real x = stress(2, 2) - _cap_start;
235  Real p = 0;
236  Real dp_dx = 0;
237  if (x > 0)
238  {
239  const Real exp_cap_rate_x = std::exp(-_cap_rate * x);
240  p = x * (1 - exp_cap_rate_x);
241  dp_dx = (1 - exp_cap_rate_x) + x * _cap_rate * exp_cap_rate_x;
242  }
243  dsmoother2 += 2 * p * dp_dx;
244  }
245  return dsmoother2;
246 }
247 
248 Real
249 TensorMechanicsPlasticWeakPlaneShear::d2smooth(const RankTwoTensor & stress) const
250 {
251  Real d2smoother2 = 0;
252  if (_tip_scheme == "cap")
253  {
254  Real x = stress(2, 2) - _cap_start;
255  Real p = 0;
256  Real dp_dx = 0;
257  Real d2p_dx2 = 0;
258  if (x > 0)
259  {
260  const Real exp_cap_rate_x = std::exp(-_cap_rate * x);
261  p = x * (1.0 - exp_cap_rate_x);
262  dp_dx = (1.0 - exp_cap_rate_x) + x * _cap_rate * exp_cap_rate_x;
263  d2p_dx2 = 2.0 * _cap_rate * exp_cap_rate_x - x * Utility::pow<2>(_cap_rate) * exp_cap_rate_x;
264  }
265  d2smoother2 += 2.0 * Utility::pow<2>(dp_dx) + 2.0 * p * d2p_dx2;
266  }
267  return d2smoother2;
268 }
269 
270 void
272  const RankTwoTensor & stress,
273  Real intnl,
274  const RankFourTensor & Eijkl,
275  std::vector<bool> & act,
276  RankTwoTensor & returned_stress) const
277 {
278  act.assign(1, false);
279  returned_stress = stress;
280 
281  if (f[0] <= _f_tol)
282  return;
283 
284  // in the following i will derive returned_stress for the case smoother=0
285 
286  Real tanpsi = tan_psi(intnl);
287  Real tanphi = tan_phi(intnl);
288 
289  // norm is the normal to the yield surface
290  // with f having psi (dilation angle) instead of phi:
291  // norm(0) = df/dsig(2,0) = df/dsig(0,2)
292  // norm(1) = df/dsig(2,1) = df/dsig(1,2)
293  // norm(2) = df/dsig(2,2)
294  std::vector<Real> norm(3, 0.0);
295  const Real tau = std::sqrt(Utility::pow<2>((stress(0, 2) + stress(2, 0)) / 2.0) +
296  Utility::pow<2>((stress(1, 2) + stress(2, 1)) / 2.0));
297  if (tau > 0.0)
298  {
299  norm[0] = 0.25 * (stress(0, 2) + stress(2, 0)) / tau;
300  norm[1] = 0.25 * (stress(1, 2) + stress(2, 1)) / tau;
301  }
302  else
303  {
304  returned_stress(2, 2) = cohesion(intnl) / tanphi;
305  act[0] = true;
306  return;
307  }
308  norm[2] = tanpsi;
309 
310  // to get the flow directions, we have to multiply norm by Eijkl.
311  // I assume that E(0,2,0,2) = E(1,2,1,2), and E(2,2,0,2) = 0 = E(0,2,1,2), etc
312  // with the usual symmetry. This makes finding the returned_stress
313  // much easier.
314  // returned_stress = stress - alpha*n
315  // where alpha is chosen so that f = 0
316  Real alpha = f[0] / (Eijkl(0, 2, 0, 2) + Eijkl(2, 2, 2, 2) * tanpsi * tanphi);
317 
318  if (1 - alpha * Eijkl(0, 2, 0, 2) / tau >= 0)
319  {
320  // returning to the "surface" of the cone
321  returned_stress(2, 2) = stress(2, 2) - alpha * Eijkl(2, 2, 2, 2) * norm[2];
322  returned_stress(0, 2) = returned_stress(2, 0) =
323  stress(0, 2) - alpha * 2.0 * Eijkl(0, 2, 0, 2) * norm[0];
324  returned_stress(1, 2) = returned_stress(2, 1) =
325  stress(1, 2) - alpha * 2.0 * Eijkl(1, 2, 1, 2) * norm[1];
326  }
327  else
328  {
329  // returning to the "tip" of the cone
330  returned_stress(2, 2) = cohesion(intnl) / tanphi;
331  returned_stress(0, 2) = returned_stress(2, 0) = returned_stress(1, 2) = returned_stress(2, 1) =
332  0.0;
333  }
334  returned_stress(0, 0) =
335  stress(0, 0) - Eijkl(0, 0, 2, 2) * (stress(2, 2) - returned_stress(2, 2)) / Eijkl(2, 2, 2, 2);
336  returned_stress(1, 1) =
337  stress(1, 1) - Eijkl(1, 1, 2, 2) * (stress(2, 2) - returned_stress(2, 2)) / Eijkl(2, 2, 2, 2);
338 
339  act[0] = true;
340 }
341 
342 std::string
344 {
345  return "WeakPlaneShear";
346 }
MooseEnum _tip_scheme
The yield function is modified to f = sqrt(s_xz^2 + s_yz^2 + a) + s_zz*_tan_phi - _cohesion where "a"...
RankTwoTensor flowPotential(const RankTwoTensor &stress, Real intnl) const override
The flow potential.
virtual Real dsmooth(const RankTwoTensor &stress) const
returns the da/dstress(2,2) - see doco for _tip_scheme
TensorMechanicsPlasticWeakPlaneShear(const InputParameters &parameters)
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.
RankTwoTensor dflowPotential_dintnl(const RankTwoTensor &stress, Real intnl) const override
The derivative of the flow potential with respect to the internal parameter.
virtual Real d2smooth(const RankTwoTensor &stress) const
returns the d^2a/dstress(2,2)^2 - see doco for _tip_scheme
virtual Real dtan_psi(const Real internal_param) const
d(tan_psi)/d(internal_param);
virtual Real tan_phi(const Real internal_param) const
tan_phi as a function of internal parameter
const TensorMechanicsHardeningModel & _tan_phi
Hardening model for tan(phi)
virtual Real derivative(Real intnl) const
RankTwoTensor df_dsig(const RankTwoTensor &stress, Real _tan_phi_or_psi) const
Function that&#39;s used in dyieldFunction_dstress and flowPotential.
Real dyieldFunction_dintnl(const RankTwoTensor &stress, Real intnl) const override
The derivative of yield function with respect to the internal parameter.
RankFourTensor dflowPotential_dstress(const RankTwoTensor &stress, Real intnl) const override
The derivative of the flow potential with respect to stress.
Real _cap_start
smoothing parameter dictating when the &#39;cap&#39; will start - see doco for _tip_scheme ...
virtual Real cohesion(const Real internal_param) const
cohesion as a function of internal parameter
const TensorMechanicsHardeningModel & _cohesion
Hardening model for cohesion.
Real yieldFunction(const RankTwoTensor &stress, Real intnl) const override
The following functions are what you should override when building single-plasticity models...
Real _small_smoother2
smoothing parameter for the cone&#39;s tip - see doco for _tip_scheme
InputParameters validParams< TensorMechanicsPlasticWeakPlaneShear >()
InputParameters validParams< TensorMechanicsPlasticModel >()
Real _cap_rate
dictates how quickly the &#39;cap&#39; degenerates to a hemisphere - see doco for _tip_scheme ...
virtual Real dtan_phi(const Real internal_param) const
d(tan_phi)/d(internal_param);
virtual Real value(Real intnl) const
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
const Real _f_tol
Tolerance on yield function.
RankTwoTensor dyieldFunction_dstress(const RankTwoTensor &stress, Real intnl) const override
The derivative of yield function with respect to stress.
Plastic Model base class The virtual functions written below must be over-ridden in derived classes t...
virtual Real tan_psi(const Real internal_param) const
tan_psi as a function of internal parameter
virtual Real smooth(const RankTwoTensor &stress) const
returns the &#39;a&#39; parameter - see doco for _tip_scheme
virtual Real dcohesion(const Real internal_param) const
d(cohesion)/d(internal_param)
const TensorMechanicsHardeningModel & _tan_psi
Hardening model for tan(psi)