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