www.mooseframework.org
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
FiniteStrainPlasticMaterial Class Reference

FiniteStrainPlasticMaterial implements rate-independent associative J2 plasticity with isotropic hardening in the finite-strain framework. More...

#include <FiniteStrainPlasticMaterial.h>

Inheritance diagram for FiniteStrainPlasticMaterial:
[legend]

Public Member Functions

 FiniteStrainPlasticMaterial (const InputParameters &parameters)
 

Protected Member Functions

virtual void computeQpStress ()
 
virtual void initQpStatefulProperties ()
 
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. More...
 
virtual Real yieldFunction (const RankTwoTensor &stress, const Real yield_stress)
 Calculates the yield function. More...
 
virtual RankTwoTensor dyieldFunction_dstress (const RankTwoTensor &stress)
 Derivative of yieldFunction with respect to the stress. More...
 
virtual Real dyieldFunction_dinternal (const Real equivalent_plastic_strain)
 Derivative of yieldFunction with respect to the equivalent plastic strain. More...
 
virtual RankTwoTensor flowPotential (const RankTwoTensor &stress)
 Flow potential, which in this case is just dyieldFunction_dstress because we are doing associative flow, and hence does not depend on the internal hardening parameter equivalent_plastic_strain. More...
 
virtual Real internalPotential ()
 The internal potential. More...
 
Real getSigEqv (const RankTwoTensor &stress)
 Equivalent stress. More...
 
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}(trial_stress - sig))_ij. More...
 
Real getYieldStress (const Real equivalent_plastic_strain)
 yield stress as a function of equivalent plastic strain. More...
 
Real getdYieldStressdPlasticStrain (const Real equivalent_plastic_strain)
 d(yieldstress)/d(equivalent plastic strain) More...
 
virtual void computeQpProperties () override
 

Protected Attributes

std::vector< Real > _yield_stress_vector
 
MaterialProperty< RankTwoTensor > & _plastic_strain
 
const MaterialProperty< RankTwoTensor > & _plastic_strain_old
 
MaterialProperty< Real > & _eqv_plastic_strain
 
const MaterialProperty< Real > & _eqv_plastic_strain_old
 
const MaterialProperty< RankTwoTensor > & _stress_old
 
const MaterialProperty< RankTwoTensor > & _strain_increment
 
const MaterialProperty< RankTwoTensor > & _rotation_increment
 
const MaterialProperty< RankFourTensor > & _elasticity_tensor
 
Real _rtol
 
Real _ftol
 
Real _eptol
 
RankFourTensor _deltaOuter
 
RankFourTensor _deltaMixed
 
const std::string _base_name
 
const std::string _elasticity_tensor_name
 
const MaterialProperty< RankTwoTensor > & _mechanical_strain
 
MaterialProperty< RankTwoTensor > & _stress
 
MaterialProperty< RankTwoTensor > & _elastic_strain
 
const MaterialProperty< RankTwoTensor > & _extra_stress
 Extra stress tensor. More...
 
std::vector< Function * > _initial_stress_fcn
 initial stress components More...
 
MaterialProperty< RankFourTensor > & _Jacobian_mult
 derivative of stress w.r.t. strain (_dstress_dstrain) More...
 

Detailed Description

FiniteStrainPlasticMaterial implements rate-independent associative J2 plasticity with isotropic hardening in the finite-strain framework.

Yield function = sqrt(3*s_ij*s_ij/2) - K(equivalent plastic strain) where s_ij = stress_ij - delta_ij*trace(stress)/3 is the deviatoric stress and K is the yield stress, specified as a piecewise-linear function by the user. Integration is performed in an incremental manner using Newton Raphson.

Definition at line 30 of file FiniteStrainPlasticMaterial.h.

Constructor & Destructor Documentation

◆ FiniteStrainPlasticMaterial()

FiniteStrainPlasticMaterial::FiniteStrainPlasticMaterial ( const InputParameters &  parameters)

Definition at line 33 of file FiniteStrainPlasticMaterial.C.

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(getMaterialProperty<RankFourTensor>(_base_name + "elasticity_tensor")),
44  _rtol(getParam<Real>("rtol")),
45  _ftol(getParam<Real>("ftol")),
46  _eptol(getParam<Real>("eptol")),
47  _deltaOuter(RankTwoTensor::Identity().outerProduct(RankTwoTensor::Identity())),
48  _deltaMixed(RankTwoTensor::Identity().mixedProductIkJl(RankTwoTensor::Identity()))
49 {
50 }
const MaterialProperty< RankTwoTensor > & _strain_increment
MaterialProperty< Real > & _eqv_plastic_strain
const MaterialProperty< RankFourTensor > & _elasticity_tensor
const MaterialProperty< Real > & _eqv_plastic_strain_old
const std::string _base_name
ComputeStressBase(const InputParameters &parameters)
const MaterialProperty< RankTwoTensor > & _plastic_strain_old
const MaterialProperty< RankTwoTensor > & _rotation_increment
MaterialProperty< RankTwoTensor > & _plastic_strain
const MaterialProperty< RankTwoTensor > & _stress_old

Member Function Documentation

◆ computeQpProperties()

void ComputeStressBase::computeQpProperties ( )
overrideprotectedvirtualinherited

Definition at line 51 of file ComputeStressBase.C.

52 {
54 
55  // Add in extra stress
56  _stress[_qp] += _extra_stress[_qp];
57 }
virtual void computeQpStress()=0
MaterialProperty< RankTwoTensor > & _stress
const MaterialProperty< RankTwoTensor > & _extra_stress
Extra stress tensor.

◆ computeQpStress()

void FiniteStrainPlasticMaterial::computeQpStress ( )
protectedvirtual

Implements ComputeStressBase.

Definition at line 61 of file FiniteStrainPlasticMaterial.C.

62 {
63 
64  // perform the return-mapping algorithm
68  _strain_increment[_qp],
69  _elasticity_tensor[_qp],
70  _stress[_qp],
72  _plastic_strain[_qp]);
73 
74  // Rotate the stress tensor to the current configuration
75  _stress[_qp] = _rotation_increment[_qp] * _stress[_qp] * _rotation_increment[_qp].transpose();
76 
77  // Rotate plastic strain tensor to the current configuration
78  _plastic_strain[_qp] =
79  _rotation_increment[_qp] * _plastic_strain[_qp] * _rotation_increment[_qp].transpose();
80 
81  // Calculate the elastic strain_increment
83 
85 }
MaterialProperty< RankFourTensor > & _Jacobian_mult
derivative of stress w.r.t. strain (_dstress_dstrain)
const MaterialProperty< RankTwoTensor > & _strain_increment
MaterialProperty< Real > & _eqv_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.
MaterialProperty< RankTwoTensor > & _stress
const MaterialProperty< RankFourTensor > & _elasticity_tensor
const MaterialProperty< Real > & _eqv_plastic_strain_old
const MaterialProperty< RankTwoTensor > & _mechanical_strain
const MaterialProperty< RankTwoTensor > & _plastic_strain_old
MaterialProperty< RankTwoTensor > & _elastic_strain
const MaterialProperty< RankTwoTensor > & _rotation_increment
MaterialProperty< RankTwoTensor > & _plastic_strain
const MaterialProperty< RankTwoTensor > & _stress_old

◆ dyieldFunction_dinternal()

Real FiniteStrainPlasticMaterial::dyieldFunction_dinternal ( const Real  equivalent_plastic_strain)
protectedvirtual

Derivative of yieldFunction with respect to the equivalent plastic strain.

Definition at line 324 of file FiniteStrainPlasticMaterial.C.

Referenced by returnMap().

325 {
326  return -getdYieldStressdPlasticStrain(equivalent_plastic_strain);
327 }
Real getdYieldStressdPlasticStrain(const Real equivalent_plastic_strain)
d(yieldstress)/d(equivalent plastic strain)

◆ dyieldFunction_dstress()

RankTwoTensor FiniteStrainPlasticMaterial::dyieldFunction_dstress ( const RankTwoTensor &  stress)
protectedvirtual

Derivative of yieldFunction with respect to the stress.

Definition at line 316 of file FiniteStrainPlasticMaterial.C.

Referenced by flowPotential(), getJac(), and returnMap().

317 {
318  RankTwoTensor deriv = sig.dsecondInvariant();
319  deriv *= std::sqrt(3.0 / sig.secondInvariant()) / 2.0;
320  return deriv;
321 }

◆ flowPotential()

RankTwoTensor FiniteStrainPlasticMaterial::flowPotential ( const RankTwoTensor &  stress)
protectedvirtual

Flow potential, which in this case is just dyieldFunction_dstress because we are doing associative flow, and hence does not depend on the internal hardening parameter equivalent_plastic_strain.

Definition at line 330 of file FiniteStrainPlasticMaterial.C.

Referenced by getJac(), and returnMap().

331 {
332  return dyieldFunction_dstress(sig); // this plasticity model assumes associative flow
333 }
virtual RankTwoTensor dyieldFunction_dstress(const RankTwoTensor &stress)
Derivative of yieldFunction with respect to the stress.

◆ getdYieldStressdPlasticStrain()

Real FiniteStrainPlasticMaterial::getdYieldStressdPlasticStrain ( const Real  equivalent_plastic_strain)
protected

d(yieldstress)/d(equivalent plastic strain)

Definition at line 414 of file FiniteStrainPlasticMaterial.C.

Referenced by dyieldFunction_dinternal().

415 {
416  unsigned nsize;
417 
418  nsize = _yield_stress_vector.size();
419 
420  if (_yield_stress_vector[0] > 0.0 || nsize % 2 > 0) // Error check for input inconsitency
421  mooseError("Error in yield stress input: Should be a vector with eqv plastic strain and yield "
422  "stress pair values.\n");
423 
424  unsigned int ind = 0;
425 
426  while (ind < nsize)
427  {
428  if (ind + 2 < nsize)
429  {
430  if (eqpe >= _yield_stress_vector[ind] && eqpe < _yield_stress_vector[ind + 2])
431  return (_yield_stress_vector[ind + 3] - _yield_stress_vector[ind + 1]) /
432  (_yield_stress_vector[ind + 2] - _yield_stress_vector[ind]);
433  }
434  else
435  return 0.0;
436 
437  ind += 2;
438  }
439 
440  return 0.0;
441 }

◆ getJac()

void FiniteStrainPlasticMaterial::getJac ( const RankTwoTensor &  sig,
const RankFourTensor &  E_ijkl,
Real  flow_incr,
RankFourTensor &  dresid_dsig 
)
protectedvirtual

Evaluates the derivative d(resid_ij)/d(sig_kl), where resid_ij = flow_incr*flowPotential_ij - (E^{-1}(trial_stress - sig))_ij.

Parameters
sigstress
E_ijklelasticity tensor (sig = E*(strain - plastic_strain))
flow_incrconsistency parameter
dresid_dsigthe required derivative (this is an output variable)

Definition at line 349 of file FiniteStrainPlasticMaterial.C.

Referenced by returnMap().

353 {
354  RankTwoTensor sig_dev, df_dsig, flow_dirn;
355  RankTwoTensor dfi_dft, dfi_dsig;
356  RankFourTensor dft_dsig, dfd_dft, dfd_dsig;
357  Real sig_eqv;
358  Real f1, f2, f3;
359  RankFourTensor temp;
360 
361  sig_dev = sig.deviatoric();
362  sig_eqv = getSigEqv(sig);
363  df_dsig = dyieldFunction_dstress(sig);
364  flow_dirn = flowPotential(sig);
365 
366  f1 = 3.0 / (2.0 * sig_eqv);
367  f2 = f1 / 3.0;
368  f3 = 9.0 / (4.0 * Utility::pow<3>(sig_eqv));
369 
370  dft_dsig = f1 * _deltaMixed - f2 * _deltaOuter - f3 * sig_dev.outerProduct(sig_dev);
371 
372  dfd_dsig = dft_dsig;
373  dresid_dsig = E_ijkl.invSymm() + dfd_dsig * flow_incr;
374 }
virtual RankTwoTensor flowPotential(const RankTwoTensor &stress)
Flow potential, which in this case is just dyieldFunction_dstress because we are doing associative fl...
Real getSigEqv(const RankTwoTensor &stress)
Equivalent stress.
virtual RankTwoTensor dyieldFunction_dstress(const RankTwoTensor &stress)
Derivative of yieldFunction with respect to the stress.

◆ getSigEqv()

Real FiniteStrainPlasticMaterial::getSigEqv ( const RankTwoTensor &  stress)
protected

Equivalent stress.

Definition at line 342 of file FiniteStrainPlasticMaterial.C.

Referenced by getJac(), and yieldFunction().

343 {
344  return std::sqrt(3 * stress.secondInvariant());
345 }

◆ getYieldStress()

Real FiniteStrainPlasticMaterial::getYieldStress ( const Real  equivalent_plastic_strain)
protected

yield stress as a function of equivalent plastic strain.

This is a piecewise linear function entered by the user in the yield_stress vector

Definition at line 378 of file FiniteStrainPlasticMaterial.C.

Referenced by returnMap().

379 {
380  unsigned nsize;
381 
382  nsize = _yield_stress_vector.size();
383 
384  if (_yield_stress_vector[0] > 0.0 || nsize % 2 > 0) // Error check for input inconsitency
385  mooseError("Error in yield stress input: Should be a vector with eqv plastic strain and yield "
386  "stress pair values.\n");
387 
388  unsigned int ind = 0;
389  Real tol = 1e-8;
390 
391  while (ind < nsize)
392  {
393  if (std::abs(eqpe - _yield_stress_vector[ind]) < tol)
394  return _yield_stress_vector[ind + 1];
395 
396  if (ind + 2 < nsize)
397  {
398  if (eqpe > _yield_stress_vector[ind] && eqpe < _yield_stress_vector[ind + 2])
399  return _yield_stress_vector[ind + 1] +
400  (eqpe - _yield_stress_vector[ind]) /
401  (_yield_stress_vector[ind + 2] - _yield_stress_vector[ind]) *
402  (_yield_stress_vector[ind + 3] - _yield_stress_vector[ind + 1]);
403  }
404  else
405  return _yield_stress_vector[nsize - 1];
406 
407  ind += 2;
408  }
409 
410  return 0.0;
411 }
const double tol
Definition: Setup.h:19

◆ initQpStatefulProperties()

void FiniteStrainPlasticMaterial::initQpStatefulProperties ( )
protectedvirtual

Reimplemented from ComputeStressBase.

Definition at line 53 of file FiniteStrainPlasticMaterial.C.

54 {
56  _plastic_strain[_qp].zero();
57  _eqv_plastic_strain[_qp] = 0.0;
58 }
MaterialProperty< Real > & _eqv_plastic_strain
virtual void initQpStatefulProperties() override
MaterialProperty< RankTwoTensor > & _plastic_strain

◆ internalPotential()

Real FiniteStrainPlasticMaterial::internalPotential ( )
protectedvirtual

The internal potential.

For associative J2 plasticity this is just -1

Definition at line 336 of file FiniteStrainPlasticMaterial.C.

Referenced by returnMap().

337 {
338  return -1;
339 }

◆ returnMap()

void FiniteStrainPlasticMaterial::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 
)
protectedvirtual

Implements the return map.

Implements the return-map algorithm via a Newton-Raphson process.

Parameters
sig_oldThe stress at the previous "time" step
eqvpstrain_oldThe equivalent plastic strain at the previous "time" step
plastic_strain_oldThe value of plastic strain at the previous "time" step
delta_dThe total strain increment for this "time" step
E_ijklThe elasticity tensor. If no plasiticity then sig_new = sig_old + E_ijkl*delta_d
sigThe stress after returning to the yield surface (this is an output variable)
eqvpstrainThe equivalent plastic strain after returning to the yield surface (this is an output variable)
plastic_strainThe value of plastic strain after returning to the yield surface (this is an output variable) Note that this algorithm doesn't do any rotations. In order to find the final stress and plastic_strain, sig and plastic_strain must be rotated using _rotation_increment.

This idea is fully explained in Simo and Hughes "Computational Inelasticity" Springer 1997, for instance, as well as many other books on plasticity. The basic idea is as follows. Given: sig_old - the stress at the start of the "time step" plastic_strain_old - the plastic strain at the start of the "time step" eqvpstrain_old - equivalent plastic strain at the start of the "time step" (In general we would be given some number of internal parameters that the yield function depends upon.) delta_d - the prescribed strain increment for this "time step" we want to determine the following parameters at the end of this "time step": sig - the stress plastic_strain - the plastic strain eqvpstrain - the equivalent plastic strain (again, in general, we would have an arbitrary number of internal parameters).

To determine these parameters, introduce the "yield function", f the "consistency parameter", flow_incr the "flow potential", flow_dirn_ij the "internal potential", internalPotential (in general there are as many internalPotential functions as there are internal parameters). All three of f, flow_dirn_ij, and internalPotential, are functions of sig and eqvpstrain. To find sig, plastic_strain and eqvpstrain, we need to solve the following resid_ij = 0 f = 0 rep = 0 This is done by using Newton-Raphson. There are 8 equations here: six from resid_ij=0 (more generally there are nine but in this case resid is symmetric); one from f=0; one from rep=0 (more generally, for N internal parameters there are N of these equations).

resid_ij = flow_incr*flow_dirn_ij - (plastic_strain - plastic_strain_old)_ij = flow_incr*flow_dirn_ij - (E^{-1}(trial_stress - sig))_ij Here trial_stress = E*(strain - plastic_strain_old) sig = E*(strain - plastic_strain) Note: flow_dirn_ij is evaluated at sig and eqvpstrain (not the old values).

f is the yield function, evaluated at sig and eqvpstrain

rep = -flow_incr*internalPotential - (eqvpstrain - eqvpstrain_old) Here internalPotential are evaluated at sig and eqvpstrain (not the old values).

The Newton-Raphson procedure has sig, flow_incr, and eqvpstrain as its variables. Therefore we need the derivatives of resid_ij, f, and rep with respect to these parameters

In this associative J2 with isotropic hardening, things are a little more specialised. (1) f = sqrt(3*s_ij*s_ij/2) - K(eqvpstrain) (this is called "isotropic hardening") (2) associativity means that flow_dirn_ij = df/d(sig_ij) = s_ij*sqrt(3/2/(s_ij*s_ij)), and this means that flow_dirn_ij*flow_dirn_ij = 3/2, so when resid_ij=0, we get (plastic_strain_dot)_ij*(plastic_strain_dot)_ij = (3/2)*flow_incr^2, where plastic_strain_dot = plastic_strain - plastic_strain_old (3) The definition of equivalent plastic strain is through eqvpstrain_dot = sqrt(2*plastic_strain_dot_ij*plastic_strain_dot_ij/3), so using (2), we obtain eqvpstrain_dot = flow_incr, and this yields internalPotential = -1 in the "rep" equation.

The linear system is ( dr_dsig flow_dirn 0 )( ddsig ) ( - resid ) ( df_dsig 0 fq )( dflow_incr ) = ( - f ) ( 0 1 -1 )( deqvpstrain ) ( - rep ) The zeroes are: d(resid_ij)/d(eqvpstrain) = flow_dirn*d(df/d(sig_ij))/d(eqvpstrain) = 0 and df/d(flow_dirn) = 0 (this is always true, even for general hardening and non-associative) and d(rep)/d(sig_ij) = -flow_incr*d(internalPotential)/d(sig_ij) = 0

Because of the zeroes and ones, the linear system is not impossible to solve by hand. NOTE: andy believes there was originally a sign-error in the next line. The next line is unchanged, however andy's definition of fq is negative of the original definition of fq. andy can't see any difference in any tests!

Definition at line 149 of file FiniteStrainPlasticMaterial.C.

Referenced by computeQpStress().

157 {
158  // the yield function, must be non-positive
159  // Newton-Raphson sets this to zero if trial stress enters inadmissible region
160  Real f;
161 
162  // the consistency parameter, must be non-negative
163  // change in plastic strain in this timestep = flow_incr*flow_potential
164  Real flow_incr = 0.0;
165 
166  // direction of flow defined by the potential
167  RankTwoTensor flow_dirn;
168 
169  // Newton-Raphson sets this zero
170  // resid_ij = flow_incr*flow_dirn_ij - (plastic_strain - plastic_strain_old)
171  RankTwoTensor resid;
172 
173  // Newton-Raphson sets this zero
174  // rep = -flow_incr*internalPotential - (eqvpstrain - eqvpstrain_old)
175  Real rep;
176 
177  // change in the stress (sig) in a Newton-Raphson iteration
178  RankTwoTensor ddsig;
179 
180  // change in the consistency parameter in a Newton-Raphson iteration
181  Real dflow_incr = 0.0;
182 
183  // change in equivalent plastic strain in one Newton-Raphson iteration
184  Real deqvpstrain = 0.0;
185 
186  // convenience variable that holds the change in plastic strain incurred during the return
187  // delta_dp = plastic_strain - plastic_strain_old
188  // delta_dp = E^{-1}*(trial_stress - sig), where trial_stress = E*(strain - plastic_strain_old)
189  RankTwoTensor delta_dp;
190 
191  // d(yieldFunction)/d(stress)
192  RankTwoTensor df_dsig;
193 
194  // d(resid_ij)/d(sigma_kl)
195  RankFourTensor dr_dsig;
196 
197  // 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
198  // zero otherwise
199  RankFourTensor dr_dsig_inv;
200 
201  // d(yieldFunction)/d(eqvpstrain)
202  Real fq;
203 
204  // yield stress at the start of this "time step" (ie, evaluated with
205  // eqvpstrain_old). It is held fixed during the Newton-Raphson return,
206  // even if eqvpstrain != eqvpstrain_old.
207  Real yield_stress;
208 
209  // measures of whether the Newton-Raphson process has converged
210  Real err1, err2, err3;
211 
212  // number of Newton-Raphson iterations performed
213  unsigned int iter = 0;
214 
215  // maximum number of Newton-Raphson iterations allowed
216  unsigned int maxiter = 100;
217 
218  // plastic loading occurs if yieldFunction > toly
219  Real toly = 1.0e-8;
220 
221  // Assume this strain increment does not induce any plasticity
222  // This is the elastic-predictor
223  sig = sig_old + E_ijkl * delta_d; // the trial stress
224  eqvpstrain = eqvpstrain_old;
225  plastic_strain = plastic_strain_old;
226 
227  yield_stress = getYieldStress(eqvpstrain); // yield stress at this equivalent plastic strain
228  if (yieldFunction(sig, yield_stress) > toly)
229  {
230  // the sig just calculated is inadmissable. We must return to the yield surface.
231  // This is done iteratively, using a Newton-Raphson process.
232 
233  delta_dp.zero();
234 
235  sig = sig_old + E_ijkl * delta_d; // this is the elastic predictor
236 
237  flow_dirn = flowPotential(sig);
238 
239  resid = flow_dirn * flow_incr - delta_dp; // Residual 1 - refer Hughes Simo
240  f = yieldFunction(sig, yield_stress);
241  rep = -eqvpstrain + eqvpstrain_old - flow_incr * internalPotential(); // Residual 3 rep=0
242 
243  err1 = resid.L2norm();
244  err2 = std::abs(f);
245  err3 = std::abs(rep);
246 
247  while ((err1 > _rtol || err2 > _ftol || err3 > _eptol) &&
248  iter < maxiter) // Stress update iteration (hardness fixed)
249  {
250  iter++;
251 
252  df_dsig = dyieldFunction_dstress(sig);
253  getJac(sig, E_ijkl, flow_incr, dr_dsig); // gets dr_dsig = d(resid_ij)/d(sig_kl)
254  fq = dyieldFunction_dinternal(eqvpstrain); // d(f)/d(eqvpstrain)
255 
267  dr_dsig_inv = dr_dsig.invSymm();
268 
276  dflow_incr = (f - df_dsig.doubleContraction(dr_dsig_inv * resid) + fq * rep) /
277  (df_dsig.doubleContraction(dr_dsig_inv * flow_dirn) - fq);
278  ddsig =
279  dr_dsig_inv *
280  (-resid -
281  flow_dirn * dflow_incr); // from solving the top row of linear system, given dflow_incr
282  deqvpstrain =
283  rep + dflow_incr; // from solving the bottom row of linear system, given dflow_incr
284 
285  // update the variables
286  flow_incr += dflow_incr;
287  delta_dp -= E_ijkl.invSymm() * ddsig;
288  sig += ddsig;
289  eqvpstrain += deqvpstrain;
290 
291  // evaluate the RHS equations ready for next Newton-Raphson iteration
292  flow_dirn = flowPotential(sig);
293  resid = flow_dirn * flow_incr - delta_dp;
294  f = yieldFunction(sig, yield_stress);
295  rep = -eqvpstrain + eqvpstrain_old - flow_incr * internalPotential();
296 
297  err1 = resid.L2norm();
298  err2 = std::abs(f);
299  err3 = std::abs(rep);
300  }
301 
302  if (iter >= maxiter)
303  mooseError("Constitutive failure");
304 
305  plastic_strain += delta_dp;
306  }
307 }
virtual Real dyieldFunction_dinternal(const Real equivalent_plastic_strain)
Derivative of yieldFunction with respect to the equivalent plastic strain.
virtual RankTwoTensor flowPotential(const RankTwoTensor &stress)
Flow potential, which in this case is just dyieldFunction_dstress because we are doing associative fl...
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.
virtual RankTwoTensor dyieldFunction_dstress(const RankTwoTensor &stress)
Derivative of yieldFunction with respect to the stress.
virtual Real internalPotential()
The internal potential.
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}...

◆ yieldFunction()

Real FiniteStrainPlasticMaterial::yieldFunction ( const RankTwoTensor &  stress,
const Real  yield_stress 
)
protectedvirtual

Calculates the yield function.

Parameters
stressthe stress at which to calculate the yield function
yield_stressthe current value of the yield stress
Returns
equivalentstress - yield_stress

Definition at line 310 of file FiniteStrainPlasticMaterial.C.

Referenced by returnMap().

311 {
312  return getSigEqv(stress) - yield_stress;
313 }
Real getSigEqv(const RankTwoTensor &stress)
Equivalent stress.

Member Data Documentation

◆ _base_name

const std::string ComputeStressBase::_base_name
protectedinherited

◆ _deltaMixed

RankFourTensor FiniteStrainPlasticMaterial::_deltaMixed
protected

Definition at line 52 of file FiniteStrainPlasticMaterial.h.

Referenced by getJac().

◆ _deltaOuter

RankFourTensor FiniteStrainPlasticMaterial::_deltaOuter
protected

Definition at line 52 of file FiniteStrainPlasticMaterial.h.

Referenced by getJac().

◆ _elastic_strain

MaterialProperty<RankTwoTensor>& ComputeStressBase::_elastic_strain
protectedinherited

◆ _elasticity_tensor

const MaterialProperty<RankFourTensor>& FiniteStrainPlasticMaterial::_elasticity_tensor
protected

Definition at line 46 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress().

◆ _elasticity_tensor_name

const std::string ComputeStressBase::_elasticity_tensor_name
protectedinherited

◆ _eptol

Real FiniteStrainPlasticMaterial::_eptol
protected

Definition at line 49 of file FiniteStrainPlasticMaterial.h.

Referenced by returnMap().

◆ _eqv_plastic_strain

MaterialProperty<Real>& FiniteStrainPlasticMaterial::_eqv_plastic_strain
protected

Definition at line 41 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

◆ _eqv_plastic_strain_old

const MaterialProperty<Real>& FiniteStrainPlasticMaterial::_eqv_plastic_strain_old
protected

Definition at line 42 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress().

◆ _extra_stress

const MaterialProperty<RankTwoTensor>& ComputeStressBase::_extra_stress
protectedinherited

Extra stress tensor.

Definition at line 47 of file ComputeStressBase.h.

Referenced by ComputeStressBase::computeQpProperties().

◆ _ftol

Real FiniteStrainPlasticMaterial::_ftol
protected

Definition at line 48 of file FiniteStrainPlasticMaterial.h.

Referenced by returnMap().

◆ _initial_stress_fcn

std::vector<Function *> ComputeStressBase::_initial_stress_fcn
protectedinherited

initial stress components

Definition at line 50 of file ComputeStressBase.h.

◆ _Jacobian_mult

MaterialProperty<RankFourTensor>& ComputeStressBase::_Jacobian_mult
protectedinherited

◆ _mechanical_strain

const MaterialProperty<RankTwoTensor>& ComputeStressBase::_mechanical_strain
protectedinherited

◆ _plastic_strain

MaterialProperty<RankTwoTensor>& FiniteStrainPlasticMaterial::_plastic_strain
protected

Definition at line 39 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

◆ _plastic_strain_old

const MaterialProperty<RankTwoTensor>& FiniteStrainPlasticMaterial::_plastic_strain_old
protected

Definition at line 40 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress().

◆ _rotation_increment

const MaterialProperty<RankTwoTensor>& FiniteStrainPlasticMaterial::_rotation_increment
protected

Definition at line 45 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress().

◆ _rtol

Real FiniteStrainPlasticMaterial::_rtol
protected

Definition at line 47 of file FiniteStrainPlasticMaterial.h.

Referenced by returnMap().

◆ _strain_increment

const MaterialProperty<RankTwoTensor>& FiniteStrainPlasticMaterial::_strain_increment
protected

Definition at line 44 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress().

◆ _stress

MaterialProperty<RankTwoTensor>& ComputeStressBase::_stress
protectedinherited

◆ _stress_old

const MaterialProperty<RankTwoTensor>& FiniteStrainPlasticMaterial::_stress_old
protected

Definition at line 43 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress().

◆ _yield_stress_vector

std::vector<Real> FiniteStrainPlasticMaterial::_yield_stress_vector
protected

Definition at line 38 of file FiniteStrainPlasticMaterial.h.

Referenced by getdYieldStressdPlasticStrain(), and getYieldStress().


The documentation for this class was generated from the following files: