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
 
void addQpInitialStress ()
 InitialStress Deprecation: remove this method. More...
 

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...
 
const bool _store_stress_old
 Parameter which decides whether to store old stress. This is required for HHT time integration and Rayleigh damping. More...
 
const bool _initial_stress_provided
 Whether initial stress was provided. InitialStress Deprecation: remove this. More...
 
MaterialProperty< RankTwoTensor > * _initial_stress
 Initial stress, if provided. InitialStress Deprecation: remove this. More...
 
const MaterialProperty< RankTwoTensor > * _initial_stress_old
 Old value of initial stress, which is needed to correctly implement finite-strain rotations. InitialStress Deprecation: remove this. 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 27 of file FiniteStrainPlasticMaterial.h.

Constructor & Destructor Documentation

FiniteStrainPlasticMaterial::FiniteStrainPlasticMaterial ( const InputParameters &  parameters)

Definition at line 28 of file FiniteStrainPlasticMaterial.C.

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 }
const MaterialProperty< RankTwoTensor > & _strain_increment
MaterialProperty< Real > & _eqv_plastic_strain
const MaterialProperty< RankFourTensor > & _elasticity_tensor
const MaterialProperty< Real > & _eqv_plastic_strain_old
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

void ComputeStressBase::addQpInitialStress ( )
protectedinherited

InitialStress Deprecation: remove this method.

Adds initial stress, if it is provided, to _stress[_qp]. This function should NOT be used if you calculate stress using

stress = stress_old + elasticity * strain_increment

because stress_old will already include initial stress. However this function SHOULD be used if your code uses

stress = elasticity * (elastic_strain_old + strain_increment) or stress = elasticity * elastic_strain

since in these cases the elastic_strain and elastic_strain_old will not include any contribution from initial stress.

Definition at line 112 of file ComputeStressBase.C.

Referenced by ComputeLinearElasticStress::computeQpStress(), ComputeCosseratLinearElasticStress::computeQpStress(), ComputeFiniteStrainElasticStress::computeQpStress(), ComputeSmearedCrackingStress::computeQpStress(), ComputeMultipleInelasticStress::computeQpStressIntermediateConfiguration(), ComputeStressBase::initQpStatefulProperties(), ComputeMultipleInelasticStress::updateQpState(), and ComputeMultipleInelasticStress::updateQpStateSingleModel().

113 {
115  _stress[_qp] += (*_initial_stress)[_qp];
116 }
MaterialProperty< RankTwoTensor > & _stress
const bool _initial_stress_provided
Whether initial stress was provided. InitialStress Deprecation: remove this.
void ComputeStressBase::computeQpProperties ( )
overrideprotectedvirtualinherited

Reimplemented in ComputeBirchMurnaghanEquationOfStress, and ComputeStressEosBase.

Definition at line 99 of file ComputeStressBase.C.

100 {
101  // InitialStress Deprecation: remove the following 2 lines
103  (*_initial_stress)[_qp] = (*_initial_stress_old)[_qp];
104 
105  computeQpStress();
106 
107  // Add in extra stress
108  _stress[_qp] += _extra_stress[_qp];
109 }
virtual void computeQpStress()=0
MaterialProperty< RankTwoTensor > & _stress
const bool _initial_stress_provided
Whether initial stress was provided. InitialStress Deprecation: remove this.
const MaterialProperty< RankTwoTensor > & _extra_stress
Extra stress tensor.
void FiniteStrainPlasticMaterial::computeQpStress ( )
protectedvirtual

Implements ComputeStressBase.

Definition at line 56 of file FiniteStrainPlasticMaterial.C.

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 }
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 > & _plastic_strain_old
const MaterialProperty< RankTwoTensor > & _rotation_increment
MaterialProperty< RankTwoTensor > & _plastic_strain
const MaterialProperty< RankTwoTensor > & _stress_old
Real FiniteStrainPlasticMaterial::dyieldFunction_dinternal ( const Real  equivalent_plastic_strain)
protectedvirtual

Derivative of yieldFunction with respect to the equivalent plastic strain.

Definition at line 316 of file FiniteStrainPlasticMaterial.C.

Referenced by returnMap().

317 {
318  return -getdYieldStressdPlasticStrain(equivalent_plastic_strain);
319 }
Real getdYieldStressdPlasticStrain(const Real equivalent_plastic_strain)
d(yieldstress)/d(equivalent plastic strain)
RankTwoTensor FiniteStrainPlasticMaterial::dyieldFunction_dstress ( const RankTwoTensor &  stress)
protectedvirtual

Derivative of yieldFunction with respect to the stress.

Definition at line 308 of file FiniteStrainPlasticMaterial.C.

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

309 {
310  RankTwoTensor deriv = sig.dsecondInvariant();
311  deriv *= std::sqrt(3.0 / sig.secondInvariant()) / 2.0;
312  return deriv;
313 }
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 322 of file FiniteStrainPlasticMaterial.C.

Referenced by getJac(), and returnMap().

323 {
324  return dyieldFunction_dstress(sig); // this plasticity model assumes associative flow
325 }
virtual RankTwoTensor dyieldFunction_dstress(const RankTwoTensor &stress)
Derivative of yieldFunction with respect to the stress.
Real FiniteStrainPlasticMaterial::getdYieldStressdPlasticStrain ( const Real  equivalent_plastic_strain)
protected

d(yieldstress)/d(equivalent plastic strain)

Definition at line 406 of file FiniteStrainPlasticMaterial.C.

Referenced by dyieldFunction_dinternal().

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 }
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 341 of file FiniteStrainPlasticMaterial.C.

Referenced by returnMap().

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 }
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.
Real FiniteStrainPlasticMaterial::getSigEqv ( const RankTwoTensor &  stress)
protected

Equivalent stress.

Definition at line 334 of file FiniteStrainPlasticMaterial.C.

Referenced by getJac(), and yieldFunction().

335 {
336  return std::sqrt(3 * stress.secondInvariant());
337 }
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 370 of file FiniteStrainPlasticMaterial.C.

Referenced by returnMap().

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 }
static const double tol
Definition: XFEMFuncs.h:26
void FiniteStrainPlasticMaterial::initQpStatefulProperties ( )
protectedvirtual

Reimplemented from ComputeStressBase.

Definition at line 48 of file FiniteStrainPlasticMaterial.C.

49 {
51  _plastic_strain[_qp].zero();
52  _eqv_plastic_strain[_qp] = 0.0;
53 }
MaterialProperty< Real > & _eqv_plastic_strain
virtual void initQpStatefulProperties() override
MaterialProperty< RankTwoTensor > & _plastic_strain
Real FiniteStrainPlasticMaterial::internalPotential ( )
protectedvirtual

The internal potential.

For associative J2 plasticity this is just -1

Definition at line 328 of file FiniteStrainPlasticMaterial.C.

Referenced by returnMap().

329 {
330  return -1;
331 }
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 141 of file FiniteStrainPlasticMaterial.C.

Referenced by computeQpStress().

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 }
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}...
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 302 of file FiniteStrainPlasticMaterial.C.

Referenced by returnMap().

303 {
304  return getSigEqv(stress) - yield_stress;
305 }
Real getSigEqv(const RankTwoTensor &stress)
Equivalent stress.

Member Data Documentation

const std::string ComputeStressBase::_base_name
protectedinherited
RankFourTensor FiniteStrainPlasticMaterial::_deltaMixed
protected

Definition at line 50 of file FiniteStrainPlasticMaterial.h.

Referenced by getJac().

RankFourTensor FiniteStrainPlasticMaterial::_deltaOuter
protected

Definition at line 50 of file FiniteStrainPlasticMaterial.h.

Referenced by getJac().

MaterialProperty<RankTwoTensor>& ComputeStressBase::_elastic_strain
protectedinherited
const MaterialProperty<RankFourTensor>& FiniteStrainPlasticMaterial::_elasticity_tensor
protected

Definition at line 44 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress().

const std::string ComputeStressBase::_elasticity_tensor_name
protectedinherited
Real FiniteStrainPlasticMaterial::_eptol
protected

Definition at line 47 of file FiniteStrainPlasticMaterial.h.

Referenced by returnMap().

MaterialProperty<Real>& FiniteStrainPlasticMaterial::_eqv_plastic_strain
protected

Definition at line 39 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

const MaterialProperty<Real>& FiniteStrainPlasticMaterial::_eqv_plastic_strain_old
protected

Definition at line 40 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress().

const MaterialProperty<RankTwoTensor>& ComputeStressBase::_extra_stress
protectedinherited
Real FiniteStrainPlasticMaterial::_ftol
protected

Definition at line 46 of file FiniteStrainPlasticMaterial.h.

Referenced by returnMap().

MaterialProperty<RankTwoTensor>* ComputeStressBase::_initial_stress
protectedinherited

Initial stress, if provided. InitialStress Deprecation: remove this.

Definition at line 74 of file ComputeStressBase.h.

Referenced by ComputeStressBase::initQpStatefulProperties().

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

initial stress components

Definition at line 62 of file ComputeStressBase.h.

Referenced by ComputeStressBase::ComputeStressBase(), and ComputeStressBase::initQpStatefulProperties().

const MaterialProperty<RankTwoTensor>* ComputeStressBase::_initial_stress_old
protectedinherited

Old value of initial stress, which is needed to correctly implement finite-strain rotations. InitialStress Deprecation: remove this.

Definition at line 77 of file ComputeStressBase.h.

const bool ComputeStressBase::_initial_stress_provided
protectedinherited
MaterialProperty<RankFourTensor>& ComputeStressBase::_Jacobian_mult
protectedinherited
const MaterialProperty<RankTwoTensor>& ComputeStressBase::_mechanical_strain
protectedinherited
MaterialProperty<RankTwoTensor>& FiniteStrainPlasticMaterial::_plastic_strain
protected

Definition at line 37 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress(), and initQpStatefulProperties().

const MaterialProperty<RankTwoTensor>& FiniteStrainPlasticMaterial::_plastic_strain_old
protected

Definition at line 38 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress().

const MaterialProperty<RankTwoTensor>& FiniteStrainPlasticMaterial::_rotation_increment
protected

Definition at line 43 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress().

Real FiniteStrainPlasticMaterial::_rtol
protected

Definition at line 45 of file FiniteStrainPlasticMaterial.h.

Referenced by returnMap().

const bool ComputeStressBase::_store_stress_old
protectedinherited

Parameter which decides whether to store old stress. This is required for HHT time integration and Rayleigh damping.

Definition at line 68 of file ComputeStressBase.h.

const MaterialProperty<RankTwoTensor>& FiniteStrainPlasticMaterial::_strain_increment
protected

Definition at line 42 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress().

MaterialProperty<RankTwoTensor>& ComputeStressBase::_stress
protectedinherited

Definition at line 53 of file ComputeStressBase.h.

Referenced by ComputeStressBase::addQpInitialStress(), ComputeMultipleInelasticCosseratStress::computeAdmissibleState(), ComputeMultipleInelasticStress::computeAdmissibleState(), ComputeStressBase::computeQpProperties(), ComputeStressEosBase::computeQpProperties(), ComputeBirchMurnaghanEquationOfStress::computeQpProperties(), ComputeLinearElasticStress::computeQpStress(), ComputeStrainIncrementBasedStress::computeQpStress(), ComputeLinearElasticPFFractureStress::computeQpStress(), ComputeCosseratLinearElasticStress::computeQpStress(), ComputeFiniteStrainElasticStress::computeQpStress(), ComputeIsotropicLinearElasticPFFractureStress::computeQpStress(), ComputeFiniteStrainElasticStressBirchMurnaghan::computeQpStress(), computeQpStress(), ComputeMultiPlasticityStress::computeQpStress(), ComputeSmearedCrackingStress::computeQpStress(), ComputeLinearViscoelasticStress::computeQpStress(), ComputeMultipleInelasticStress::computeQpStressIntermediateConfiguration(), ComputeSmearedCrackingStress::crackingStressRotation(), ComputeMultipleInelasticStress::finiteStrainRotation(), ComputeStressBase::initQpStatefulProperties(), FiniteStrainCrystalPlasticity::initQpStatefulProperties(), FiniteStrainUObasedCP::initQpStatefulProperties(), FiniteStrainHyperElasticViscoPlastic::initQpStatefulProperties(), ComputeMultiPlasticityStress::postReturnMap(), FiniteStrainUObasedCP::postSolveQp(), FiniteStrainHyperElasticViscoPlastic::postSolveQp(), FiniteStrainCrystalPlasticity::postSolveQp(), ComputeMultipleInelasticStress::updateQpState(), ComputeMultipleInelasticStress::updateQpStateSingleModel(), and LinearIsoElasticPFDamage::updateVar().

const MaterialProperty<RankTwoTensor>& FiniteStrainPlasticMaterial::_stress_old
protected

Definition at line 41 of file FiniteStrainPlasticMaterial.h.

Referenced by computeQpStress().

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

Definition at line 36 of file FiniteStrainPlasticMaterial.h.

Referenced by getdYieldStressdPlasticStrain(), and getYieldStress().


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