www.mooseframework.org
FiniteStrainPlasticMaterial.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<ComputeStressBase>();
15 
16  params.addRequiredParam<std::vector<Real>>(
17  "yield_stress",
18  "Input data as pairs of equivalent plastic strain and yield stress: Should "
19  "start with equivalent plastic strain 0");
20  params.addParam<Real>("rtol", 1e-8, "Plastic strain NR tolerance");
21  params.addParam<Real>("ftol", 1e-4, "Consistency condition NR tolerance");
22  params.addParam<Real>("eptol", 1e-7, "Equivalent plastic strain NR tolerance");
23  params.addClassDescription("Associative J2 plasticity with isotropic hardening.");
24 
25  return params;
26 }
27 
29  : ComputeStressBase(parameters),
30  _yield_stress_vector(getParam<std::vector<Real>>("yield_stress")), // Read from input file
31  _plastic_strain(declareProperty<RankTwoTensor>("plastic_strain")),
32  _plastic_strain_old(getMaterialPropertyOld<RankTwoTensor>("plastic_strain")),
33  _eqv_plastic_strain(declareProperty<Real>("eqv_plastic_strain")),
34  _eqv_plastic_strain_old(getMaterialPropertyOld<Real>("eqv_plastic_strain")),
35  _stress_old(getMaterialPropertyOld<RankTwoTensor>("stress")),
36  _strain_increment(getMaterialProperty<RankTwoTensor>("strain_increment")),
37  _rotation_increment(getMaterialProperty<RankTwoTensor>("rotation_increment")),
38  _elasticity_tensor(getMaterialProperty<RankFourTensor>("elasticity_tensor")),
39  _rtol(getParam<Real>("rtol")),
40  _ftol(getParam<Real>("ftol")),
41  _eptol(getParam<Real>("eptol")),
42  _deltaOuter(RankTwoTensor::Identity().outerProduct(RankTwoTensor::Identity())),
43  _deltaMixed(RankTwoTensor::Identity().mixedProductIkJl(RankTwoTensor::Identity()))
44 {
45 }
46 
47 void
49 {
51  _plastic_strain[_qp].zero();
52  _eqv_plastic_strain[_qp] = 0.0;
53 }
54 
55 void
57 {
58 
59  // perform the return-mapping algorithm
63  _strain_increment[_qp],
64  _elasticity_tensor[_qp],
65  _stress[_qp],
67  _plastic_strain[_qp]);
68 
69  // Rotate the stress tensor to the current configuration
70  _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
71 
72  // Rotate plastic strain tensor to the current configuration
73  _plastic_strain[_qp] =
74  _rotation_increment[_qp] * _plastic_strain[_qp] * _rotation_increment[_qp].transpose();
75 
77 }
78 
140 void
141 FiniteStrainPlasticMaterial::returnMap(const RankTwoTensor & sig_old,
142  const Real eqvpstrain_old,
143  const RankTwoTensor & plastic_strain_old,
144  const RankTwoTensor & delta_d,
145  const RankFourTensor & E_ijkl,
146  RankTwoTensor & sig,
147  Real & eqvpstrain,
148  RankTwoTensor & plastic_strain)
149 {
150  // the yield function, must be non-positive
151  // Newton-Raphson sets this to zero if trial stress enters inadmissible region
152  Real f;
153 
154  // the consistency parameter, must be non-negative
155  // change in plastic strain in this timestep = flow_incr*flow_potential
156  Real flow_incr = 0.0;
157 
158  // direction of flow defined by the potential
159  RankTwoTensor flow_dirn;
160 
161  // Newton-Raphson sets this zero
162  // resid_ij = flow_incr*flow_dirn_ij - (plastic_strain - plastic_strain_old)
163  RankTwoTensor resid;
164 
165  // Newton-Raphson sets this zero
166  // rep = -flow_incr*internalPotential - (eqvpstrain - eqvpstrain_old)
167  Real rep;
168 
169  // change in the stress (sig) in a Newton-Raphson iteration
170  RankTwoTensor ddsig;
171 
172  // change in the consistency parameter in a Newton-Raphson iteration
173  Real dflow_incr = 0.0;
174 
175  // change in equivalent plastic strain in one Newton-Raphson iteration
176  Real deqvpstrain = 0.0;
177 
178  // convenience variable that holds the change in plastic strain incurred during the return
179  // delta_dp = plastic_strain - plastic_strain_old
180  // delta_dp = E^{-1}*(trial_stress - sig), where trial_stress = E*(strain - plastic_strain_old)
181  RankTwoTensor delta_dp;
182 
183  // d(yieldFunction)/d(stress)
184  RankTwoTensor df_dsig;
185 
186  // d(resid_ij)/d(sigma_kl)
187  RankFourTensor dr_dsig;
188 
189  // dr_dsig_inv_ijkl*dr_dsig_klmn = 0.5*(de_ij de_jn + de_ij + de_jm), where de_ij = 1 if i=j, but
190  // zero otherwise
191  RankFourTensor dr_dsig_inv;
192 
193  // d(yieldFunction)/d(eqvpstrain)
194  Real fq;
195 
196  // yield stress at the start of this "time step" (ie, evaluated with
197  // eqvpstrain_old). It is held fixed during the Newton-Raphson return,
198  // even if eqvpstrain != eqvpstrain_old.
199  Real yield_stress;
200 
201  // measures of whether the Newton-Raphson process has converged
202  Real err1, err2, err3;
203 
204  // number of Newton-Raphson iterations performed
205  unsigned int iter = 0;
206 
207  // maximum number of Newton-Raphson iterations allowed
208  unsigned int maxiter = 100;
209 
210  // plastic loading occurs if yieldFunction > toly
211  Real toly = 1.0e-8;
212 
213  // Assume this strain increment does not induce any plasticity
214  // This is the elastic-predictor
215  sig = sig_old + E_ijkl * delta_d; // the trial stress
216  eqvpstrain = eqvpstrain_old;
217  plastic_strain = plastic_strain_old;
218 
219  yield_stress = getYieldStress(eqvpstrain); // yield stress at this equivalent plastic strain
220  if (yieldFunction(sig, yield_stress) > toly)
221  {
222  // the sig just calculated is inadmissable. We must return to the yield surface.
223  // This is done iteratively, using a Newton-Raphson process.
224 
225  delta_dp.zero();
226 
227  sig = sig_old + E_ijkl * delta_d; // this is the elastic predictor
228 
229  flow_dirn = flowPotential(sig);
230 
231  resid = flow_dirn * flow_incr - delta_dp; // Residual 1 - refer Hughes Simo
232  f = yieldFunction(sig, yield_stress);
233  rep = -eqvpstrain + eqvpstrain_old - flow_incr * internalPotential(); // Residual 3 rep=0
234 
235  err1 = resid.L2norm();
236  err2 = std::abs(f);
237  err3 = std::abs(rep);
238 
239  while ((err1 > _rtol || err2 > _ftol || err3 > _eptol) &&
240  iter < maxiter) // Stress update iteration (hardness fixed)
241  {
242  iter++;
243 
244  df_dsig = dyieldFunction_dstress(sig);
245  getJac(sig, E_ijkl, flow_incr, dr_dsig); // gets dr_dsig = d(resid_ij)/d(sig_kl)
246  fq = dyieldFunction_dinternal(eqvpstrain); // d(f)/d(eqvpstrain)
247 
259  dr_dsig_inv = dr_dsig.invSymm();
260 
268  dflow_incr = (f - df_dsig.doubleContraction(dr_dsig_inv * resid) + fq * rep) /
269  (df_dsig.doubleContraction(dr_dsig_inv * flow_dirn) - fq);
270  ddsig =
271  dr_dsig_inv *
272  (-resid -
273  flow_dirn * dflow_incr); // from solving the top row of linear system, given dflow_incr
274  deqvpstrain =
275  rep + dflow_incr; // from solving the bottom row of linear system, given dflow_incr
276 
277  // update the variables
278  flow_incr += dflow_incr;
279  delta_dp -= E_ijkl.invSymm() * ddsig;
280  sig += ddsig;
281  eqvpstrain += deqvpstrain;
282 
283  // evaluate the RHS equations ready for next Newton-Raphson iteration
284  flow_dirn = flowPotential(sig);
285  resid = flow_dirn * flow_incr - delta_dp;
286  f = yieldFunction(sig, yield_stress);
287  rep = -eqvpstrain + eqvpstrain_old - flow_incr * internalPotential();
288 
289  err1 = resid.L2norm();
290  err2 = std::abs(f);
291  err3 = std::abs(rep);
292  }
293 
294  if (iter >= maxiter)
295  mooseError("Constitutive failure");
296 
297  plastic_strain += delta_dp;
298  }
299 }
300 
301 Real
302 FiniteStrainPlasticMaterial::yieldFunction(const RankTwoTensor & stress, const Real yield_stress)
303 {
304  return getSigEqv(stress) - yield_stress;
305 }
306 
307 RankTwoTensor
309 {
310  RankTwoTensor deriv = sig.dsecondInvariant();
311  deriv *= std::sqrt(3.0 / sig.secondInvariant()) / 2.0;
312  return deriv;
313 }
314 
315 Real
316 FiniteStrainPlasticMaterial::dyieldFunction_dinternal(const Real equivalent_plastic_strain)
317 {
318  return -getdYieldStressdPlasticStrain(equivalent_plastic_strain);
319 }
320 
321 RankTwoTensor
323 {
324  return dyieldFunction_dstress(sig); // this plasticity model assumes associative flow
325 }
326 
327 Real
329 {
330  return -1;
331 }
332 
333 Real
334 FiniteStrainPlasticMaterial::getSigEqv(const RankTwoTensor & stress)
335 {
336  return std::sqrt(3 * stress.secondInvariant());
337 }
338 
339 // Jacobian for stress update algorithm
340 void
341 FiniteStrainPlasticMaterial::getJac(const RankTwoTensor & sig,
342  const RankFourTensor & E_ijkl,
343  Real flow_incr,
344  RankFourTensor & dresid_dsig)
345 {
346  RankTwoTensor sig_dev, df_dsig, flow_dirn;
347  RankTwoTensor dfi_dft, dfi_dsig;
348  RankFourTensor dft_dsig, dfd_dft, dfd_dsig;
349  Real sig_eqv;
350  Real f1, f2, f3;
351  RankFourTensor temp;
352 
353  sig_dev = sig.deviatoric();
354  sig_eqv = getSigEqv(sig);
355  df_dsig = dyieldFunction_dstress(sig);
356  flow_dirn = flowPotential(sig);
357 
358  f1 = 3.0 / (2.0 * sig_eqv);
359  f2 = f1 / 3.0;
360  f3 = 9.0 / (4.0 * Utility::pow<3>(sig_eqv));
361 
362  dft_dsig = f1 * _deltaMixed - f2 * _deltaOuter - f3 * sig_dev.outerProduct(sig_dev);
363 
364  dfd_dsig = dft_dsig;
365  dresid_dsig = E_ijkl.invSymm() + dfd_dsig * flow_incr;
366 }
367 
368 // Obtain yield stress for a given equivalent plastic strain (input)
369 Real
371 {
372  unsigned nsize;
373 
374  nsize = _yield_stress_vector.size();
375 
376  if (_yield_stress_vector[0] > 0.0 || nsize % 2 > 0) // Error check for input inconsitency
377  mooseError("Error in yield stress input: Should be a vector with eqv plastic strain and yield "
378  "stress pair values.\n");
379 
380  unsigned int ind = 0;
381  Real tol = 1e-8;
382 
383  while (ind < nsize)
384  {
385  if (std::abs(eqpe - _yield_stress_vector[ind]) < tol)
386  return _yield_stress_vector[ind + 1];
387 
388  if (ind + 2 < nsize)
389  {
390  if (eqpe > _yield_stress_vector[ind] && eqpe < _yield_stress_vector[ind + 2])
391  return _yield_stress_vector[ind + 1] +
392  (eqpe - _yield_stress_vector[ind]) /
393  (_yield_stress_vector[ind + 2] - _yield_stress_vector[ind]) *
394  (_yield_stress_vector[ind + 3] - _yield_stress_vector[ind + 1]);
395  }
396  else
397  return _yield_stress_vector[nsize - 1];
398 
399  ind += 2;
400  }
401 
402  return 0.0;
403 }
404 
405 Real
407 {
408  unsigned nsize;
409 
410  nsize = _yield_stress_vector.size();
411 
412  if (_yield_stress_vector[0] > 0.0 || nsize % 2 > 0) // Error check for input inconsitency
413  mooseError("Error in yield stress input: Should be a vector with eqv plastic strain and yield "
414  "stress pair values.\n");
415 
416  unsigned int ind = 0;
417 
418  while (ind < nsize)
419  {
420  if (ind + 2 < nsize)
421  {
422  if (eqpe >= _yield_stress_vector[ind] && eqpe < _yield_stress_vector[ind + 2])
423  return (_yield_stress_vector[ind + 3] - _yield_stress_vector[ind + 1]) /
424  (_yield_stress_vector[ind + 2] - _yield_stress_vector[ind]);
425  }
426  else
427  return 0.0;
428 
429  ind += 2;
430  }
431 
432  return 0.0;
433 }
MaterialProperty< RankFourTensor > & _Jacobian_mult
derivative of stress w.r.t. strain (_dstress_dstrain)
ComputeStressBase is the base class for stress tensors.
const MaterialProperty< RankTwoTensor > & _strain_increment
MaterialProperty< Real > & _eqv_plastic_strain
virtual Real dyieldFunction_dinternal(const Real equivalent_plastic_strain)
Derivative of yieldFunction with respect to the equivalent plastic strain.
virtual void returnMap(const RankTwoTensor &sig_old, const Real eqvpstrain_old, const RankTwoTensor &plastic_strain_old, const RankTwoTensor &delta_d, const RankFourTensor &E_ijkl, RankTwoTensor &sig, Real &eqvpstrain, RankTwoTensor &plastic_strain)
Implements the return map.
virtual RankTwoTensor flowPotential(const RankTwoTensor &stress)
Flow potential, which in this case is just dyieldFunction_dstress because we are doing associative fl...
MaterialProperty< RankTwoTensor > & _stress
const MaterialProperty< RankFourTensor > & _elasticity_tensor
const MaterialProperty< Real > & _eqv_plastic_strain_old
virtual Real yieldFunction(const RankTwoTensor &stress, const Real yield_stress)
Calculates the yield function.
Real getYieldStress(const Real equivalent_plastic_strain)
yield stress as a function of equivalent plastic strain.
static const double tol
Definition: XFEMFuncs.h:26
FiniteStrainPlasticMaterial(const InputParameters &parameters)
InputParameters validParams< ComputeStressBase >()
virtual void initQpStatefulProperties() override
const MaterialProperty< RankTwoTensor > & _plastic_strain_old
Real getSigEqv(const RankTwoTensor &stress)
Equivalent stress.
virtual RankTwoTensor dyieldFunction_dstress(const RankTwoTensor &stress)
Derivative of yieldFunction with respect to the stress.
virtual Real internalPotential()
The internal potential.
InputParameters validParams< FiniteStrainPlasticMaterial >()
const MaterialProperty< RankTwoTensor > & _rotation_increment
Real getdYieldStressdPlasticStrain(const Real equivalent_plastic_strain)
d(yieldstress)/d(equivalent plastic strain)
MaterialProperty< RankTwoTensor > & _plastic_strain
virtual void getJac(const RankTwoTensor &sig, const RankFourTensor &E_ijkl, Real flow_incr, RankFourTensor &dresid_dsig)
Evaluates the derivative d(resid_ij)/d(sig_kl), where resid_ij = flow_incr*flowPotential_ij - (E^{-1}...
const MaterialProperty< RankTwoTensor > & _stress_old