4.17. Plane Stress Rate Plasticity Model

4.17.1. Theory

The plane stress rate plasticity model is the plane stress formulation of a \(J_2\) plasticity model given by Simo and Taylor [[1]] (and described again in Simo and Hughes [[2]]) extended to include rate-dependent hardening and a failure model for use with shell elements.

Like other plasticity models, the components of the objective stress rate, \(\stackrel{\circ}{\sigma}_{ij}\), are written as,

\[\stackrel{\circ}{\sigma}_{ij} = \mathbb{C}_{ijkl}D_{kl}^{\text{e}}\]

where \(\mathbb{C}_{ijkl}\) are the components of the fourth-order, isotropic elasticity tensor and \(D_{ij}^{\text{e}}\) are the components of the elastic part of the total rate of deformation tensor. An additive split of the total rate of deformation tensor into elastic and plastic contributions is assumed such that,

\[D_{ij}=D_{ij}^{\text{e}}+D_{ij}^{\text{p}}.\]

The plane stress formulation recasts the three-dimensional problem into a constrained subspace with plane stress conditions acting as the constraints. To do this, the plane stress rate plasticity model follows the approach of Simo and Taylor [[1]] to enforce \(\sigma_{13}=\sigma_{23}=\sigma_{33}=0\) and related conditions.

For the plasticity portion of the model, the formulation of Simo and Taylor [[1]] is used (In the work of Simo and Taylor [[1]] and later Simo and Hughes [[2]], hardening is assumed to be rate and temperature independent. Here, such terms are included but do not materially change the formulation. Similarly, the earlier works also introduce kinematic hardening which is not used in the current model.) in which a traditional three-dimensional \(J_2\) plasticity model is recast in reduced subspace. To do this, it is recalled that in three-dimensions the von Mises effective stress, \(\phi\), is written,

\[\phi^2=\frac{3}{2}s_{ij}s_{ij},\]

with \(s_{ij}=\sigma_{ij}-(1/3)\sigma_{kk}\delta_{ij}\) the deviatoric stress. To write an equivalent expression in the reduced subspace, the vector, \(\underline{\sigma}\), and matrix, \(\underline{\underline{\bar{P}}}\), are introduced as (Note, here the \(\underline{x}\) and \(\underline{\underline{X}}\) notations are introduced for vector and matrix objects, respectively, to clearly distinguish that these variables are not tensors. This results from operating in the constrained stress subspace and means that these terms do not have properties of a tensor and act on each other as traditional matrices and vectors.),

\[\begin{split}\underline{\sigma}=\left[\begin{array}{c}\sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{array}\right],~\qquad~;~\qquad~ \underline{\underline{\bar{P}}}=\frac{1}{3}\left[\begin{array}{ccc}2 & -1 & 0 \\ -1 & 2 & 0 \\ 0 & 0 & 3\end{array}\right],\end{split}\]

such that,

\[\begin{split}\underline{s}=\left[\begin{array}{c} s_{11} \\ s_{22} \\ s_{12} \end{array}\right] = \underline{\underline{\bar{P}}}\;\underline{\sigma}.\end{split}\]

In the reduced plane-stress subspace, an alternative effective stress, \(\bar{\phi}\), is given as,

\[\bar{\phi}^2=\underline{\sigma}^T\underline{\underline{P}}\;\underline{\sigma}=\frac{2}{3}\phi^2,\]

where

\[\begin{split}\underline{\underline{P}}=\frac{1}{3}\left[\begin{array}{ccc}2 & -1 & 0 \\ -1 & 2 & 0 \\ 0 & 0 & 6\end{array}\right],\end{split}\]

in which \(\underline{\underline{P}}\) and \(\underline{\underline{\bar{P}}}\) differ by a two in the shear term to reflect Voigt corrections.

A yield function, \(f\), is introduced as,

(4.63)\[f=\bar{\phi}^2-R^2,\]

with,

\[R=\sqrt{\frac{2}{3}}\bar{\sigma}\left(\bar{\varepsilon}^p,\dot{\bar{\varepsilon}}^p,\theta\right),\]

where \(\bar{\varepsilon}^p\) and \(\dot{\bar{\varepsilon}}^p\) are the equivalent plastic strain (isotropic hardening variable) and its rate, respectively. Various hardening options may be used with this model. In general, the current flow stress is written as,

(4.64)\[\bar{\sigma}\left(\bar{\varepsilon}^p,\dot{\bar{\varepsilon}}^p,\theta\right)=\left(\sigma_y + K\left(\bar{\varepsilon}^p\right)\right) \hat{\sigma}\left(\dot{\bar{\varepsilon}}^p\right)\left(1-\left(\frac{\theta-\theta_{\text{ref}}}{\theta_{\text{melt}}-\theta_{\text{ref}}}\right)^M\right),\]

in which \(\sigma_y\) is the original yield stress, \(K\) is the isotropic hardening function that may take linear, power-law, or multilinear form, \(\hat{\sigma}\) the rate multiplier whose specification will be defined later, and the right-most term is the Johnson-Cook temperature dependence term that may be optionally used to give temperature dependence of the flow-stress.

To complete the theoretical formulation, the flow rules are specified as,

\[\begin{split}d\underline{\varepsilon}^p & = \lambda\underline{\underline{P}} \; \underline{\sigma}, \\ d\bar{\varepsilon}^p & = \dot{\bar{\varepsilon}}^{p}\Delta t = \lambda\sqrt{\frac{2}{3}}\phi,\end{split}\]

where \(\lambda\) is the consistency multiplier enforcing \(f=0\) during plastic deformation and \(d\underline{\varepsilon}^p\) is the plastic strain increment in the constrained subspace. It is emphasized here that the yield surface described in (4.63) is not homogeneous of degree one like in other three-dimensional formulations presented in this manual. As such, the consistency multiplier and equivalent plastic strain increment are not equivalent. As an example of this, by consideration of the preceding relations, it is apparent that \(\lambda\) has units of one over stress.

The specification of the rate dependence, \(\hat{\sigma}\), is important as it enables the consideration of two different model responses. These behaviors are controlled via the USER RATE DEPENDENCE command. If this input parameter is zero, then either an analytical or user-defined rate-dependence may be given. Importantly, in this case failure is not modeled. For the analytical case, the Johnson-Cook [[3], [4]] rate-multiplier is used such that,

\[\begin{split}\hat{\sigma}\left(\dot{\bar{\varepsilon}}^p\right)=\left\{\begin{array}{cc} 1+C\ln\left(\frac{\dot{\bar{\varepsilon}}^p}{\dot{\bar{\varepsilon}}_0}\right) & \dot{\bar{\varepsilon}}^p > \dot{\bar{\varepsilon}}_0 \\ 1 & \dot{\bar{\varepsilon}}^p \leq \dot{\bar{\varepsilon}}_0 \end{array}\right. ,\end{split}\]

with \(C\) being the rate dependence multiplier and \(\dot{\bar{\varepsilon}}_0\) is a reference rate. Note, while other models allow user specification of the reference rate, the plane stress rate plasticity model uses the value set in the original work of Johnson-Cook [[3]] such that \(\dot{\bar{\varepsilon}}_0=1 \text{s}^{-1}\). Alternatively, a user function may be specified for the rate multiplier, \(\hat{\sigma}\).

If USER RATE DEPENDENCE is set to one, both rate dependence and failure may be modeled. With respect to the rate dependence, (4.64) is rewritten,

\[\bar{\sigma}\left(\bar{\varepsilon}^p,\dot{\bar{\varepsilon}}^p,\theta\right)=\tilde{\sigma}\left(\bar{\varepsilon}^p,\dot{\bar{\varepsilon}}^p\right)\left(1-\left(\frac{\theta-\theta_{\text{ref}}}{\theta_{\text{melt}}-\theta_{\text{ref}}}\right)^M\right),\]

in which both isotropic hardening and rate dependence are described via definition of \(\tilde{\sigma}\). In this case, \(\tilde{\sigma}\) cannot be specified through analytical expressions and must instead be given as a series of isotropic hardening curves; each at a different strain rate. For rates not explicitly given, interpolation is performed between relevant curves. Note, no extrapolation is performed with respect to the rates. If a rate is determined outside any specified curves, the hardening is calculated with respect to the bounding curve.

For failure, a failure parameter, \(\alpha\), is calculated as

(4.65)\[\alpha=\int_0^t \frac{1}{\varepsilon_f\left(\eta,\dot{\bar{\varepsilon}}^p\right)}\dot{\bar{\varepsilon}}^pdt = \sum_{t_0}^{t}\frac{d\bar{\varepsilon}^p}{\varepsilon_f\left(\eta,\dot{\bar{\varepsilon}}^p\right)},\]

in which the summation is used to imply the discrete calculation of the damage variable over a series of loadsteps and \(\varepsilon_f\) is the rate and triaxiality, \(\eta\), dependent failure strain. The failure strain, \(\varepsilon_f\), is specified in a fashion similar to \(\tilde{\sigma}\). Specifically, a series of triaxiality dependent functions are defined each at a given strain rate. Interpolation is used at rates between those specified. Extrapolation outside the defined bounds is not done and the extremum curves are instead used.

The onset of damage is assumed to occur when \(\alpha = 1\) and the current failure strain is taken to be the critical one such that \(\varepsilon_f^{cr}=\varepsilon_f\left(t=t_{cr}\right)\) with \(t_{cr}\) being the time at which \(\alpha=1\). Subsequent damage calculation is performed via,

(4.66)\[\alpha\left(t>t_{cr}\right)=\int_{0}^{t_{cr}}\frac{1}{\varepsilon_f\left(\eta,\dot{\bar{\varepsilon}}^p\right)}\dot{\bar{\varepsilon}}^pdt + \int_{t_{cr}}^{t}\frac{1}{\varepsilon_f^{cr}}\dot{\bar{\varepsilon}}^pdt.\]

After the critical failure parameter has been reached, an exponential decay relation is used to decrease the strength of the material. In this fashion, a decay relation of the form,

(4.67)\[\bar{\sigma}=\bar{\sigma}e^{C_1\left(1-\alpha\right)},\]

is used in which \(C_1\) is the decay coefficient.

For more information about the plane stress rate plasticity model, consult [[1], [2]].

4.17.2. Implementation

The plane stress rate plasticity model encapsulates both a plasticity and failure model. These features are implemented in a decoupled, sequential sense. As such, the implementation of these features will also be presented and discussed in a sequential fashion.

For the plasticity portion, the approach of Simo and Taylor [[1]] (and Simo and Hughes [[2]]) in developing a single scalar equation to solve is adopted. As will be discussed, a slightly different approach will be used to solve this equation versus that used previously. To get to this single scalar equation, an elastic-predictor inelastic corrector scheme is adopted. In this scheme, an elastic predictor is calculated by assuming all deformation is elastic such that,

(4.68)\[\underline{\sigma}^{tr}=\underline{\sigma}^{n}+\Delta t\underline{\underline{\bar{C}}}\;d\underline{\varepsilon}^{n+1},\]

in which n and n+1 denote the material states at \(t=t_n\) and \(t=t_{n+1}\), respectively, with \(\Delta t = t_{n+1}-t_{n}\). The plane stress stiffness matrix, \(\underline{\underline{\bar{C}}}\), is given as,

\[\begin{split}\underline{\underline{\bar{C}}}=\frac{E}{1-\nu^2}\left[\begin{array}{ccc} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{array}\right],\end{split}\]

with \(E\) and \(\nu\) being the Youngs Modulus and Poisson’s ratio, respectively, and \(d\underline{\varepsilon}^{n+1}\) is the plane stress total strain increment that is written,

\[\begin{split}d\underline{\varepsilon}=\Delta t\left[\begin{array}{c} d_{11} \\ d_{22} \\ 2d_{12} \end{array}\right],\end{split}\]

where \(d_{ij}\) are the components of the rate of deformation tensor.

The stress at time \(t=t_{n+1}\) may be given as,

\[\underline{\sigma}_{n+1}=\underline{\underline{\bar{C}}}\left[\underline{\varepsilon}^{n+1}-\underline{\varepsilon}^{\text{p}}\right],\]

which noting the definition of the trial stress in (4.68) may be implicitly rewritten,

\[\underline{\sigma}^{n+1} = \underline{\sigma}^{tr}-\lambda\,\underline{\underline{\bar{C}}}\;\underline{\underline{P}}\;\underline{\sigma}^{n+1}.\]

Rearranging yields,

(4.69)\[\left[\underline{\underline{I}}+\lambda\,\underline{\underline{\bar{C}}}\;\underline{\underline{P}}\right]\underline{\sigma}^{n+1}=\underline{\sigma}^{tr},\]

with \(\underline{\underline{I}}\) the identity matrix. Importantly, by noting that \(\underline{\sigma}^{tr}\) is known it is clear that (4.69) is an equation for the updated stress vector in terms of only the unknown scalar consistency parameter, \(\lambda\). To further simplify the problem, it can be shown that \(\underline{\underline{\bar{C}}}\) and \(\underline{\underline{P}}\) share the same principal subspaces such that (see Simo and Hughes [[2]] for details),

\[\underline{\underline{P}}=\underline{\underline{Q}}\;\underline{\underline{\Lambda}}^P\underline{\underline{Q}}^T, ~\qquad~;~\qquad~ \underline{\underline{\bar{C}}}=\underline{\underline{Q}}\;\underline{\underline{\Lambda}}^C\underline{\underline{Q}}^T,\]

where \(\underline{\underline{Q}}^T\) is an orthogonal matrix such that \(\underline{\underline{Q}}^T=\underline{\underline{Q}}^{-1}\) and the matrices \(\underline{\underline{Q}},~\underline{\underline{\Lambda}}^P\) and \(\underline{\underline{\Lambda}}^C\) are given as,

\[\begin{split}\underline{\underline{\Lambda}}^P & = \left[\begin{array}{ccc} 1/3 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 2\end{array}\right], ~~;~~ \underline{\underline{\Lambda}}^C=\left[\begin{array}{ccc} E/\left(1-\nu\right) & 0 & 0 \\ 0 & 2\mu & 0 \\ 0 & 0 & 2\mu\end{array}\right] \\ \qquad\qquad & \underline{\underline{Q}} = \frac{1}{\sqrt{2}}\left[\begin{array}{ccc} 1 & -1 & 0 \\ 1 & 1 & 0 \\ 0 & 0 & \sqrt{2} \end{array}\right].\end{split}\]

By introducing a transformed stress vector in the principal space of \(\underline{\underline{\bar{C}}}\) and \(\underline{\underline{P}}\), \(\eta\), such that,

\[\begin{split}\underline{\eta} = \left[\begin{array}{c}\eta_{11} \\ \eta_{12} \\ \eta_{12} \end{array}\right]=\underline{\underline{Q}}^T\underline{\sigma},\end{split}\]

the effective stress may be rewritten,

\[\bar{\phi}^2=\underline{\eta}^T\underline{\underline{\Lambda}}^P\underline{\eta}.\]

By rewriting (4.69) in the transformed space the premultiplying matrix on the left-hand side can be analytically inverted such that,

\[\begin{split}\underline{\eta}^{n+1}=\left[\begin{array}{c}\frac{\eta^{tr}_{11}}{1+\lambda\frac{E}{3\left(1-\nu\right)}} \\ \frac{\eta_{22}^{tr}}{1+\lambda2\mu} \\ \frac{\eta_{12}^{tr}}{1+\lambda2\mu}\end{array}\right],\end{split}\]

and the effective stress may be written as a scalar function of \(\lambda\),

\[\bar{\phi}^2\left(\lambda\right)=\frac{\frac{1}{3}\left(\eta_{11}^{tr}\right)^2}{\left[1+\lambda\frac{E}{2\left(1-\nu\right)}\right]^2} \frac{\left(\eta_{22}^{tr}\right)^2+2\left(\eta_{12}^{tr}\right)^2}{\left[1+\lambda2\mu\right]^2}.\]

Noting that the equivalent plastic strain and rate may be written,

\[\bar{\varepsilon}^{p\left(n+1\right)}=\bar{\varepsilon}^{p\left(n\right)}+\lambda\sqrt{\frac{2}{3}}\phi,~\qquad~;~\qquad~ \dot{\bar{\varepsilon}}^p=\frac{\lambda}{\Delta t}\sqrt{\frac{2}{3}}\phi,\]

means determining the updated states reduces to solving the scalar consistency equation,

\[f\left(\lambda\right)=\bar{\phi}^2\left(\lambda\right)-\bar{\sigma}\left(\lambda\right),\]

for \(\lambda\). This is done iteratively by using a line-search augmented Newton-Raphson method like that described in [[5]].

Failure is handled separately from plasticity and in a straight-forward fashion. Specifically, if \(\alpha^n>1\) (above the critical value) then a decay coefficient, \(\beta^{n+1}\), is calculated via

\[\beta^{n+1}=e^{C_1\left(1-\alpha^n\right)},\]

and the yield stress is scaled accordingly such that,

\[\bar{\sigma}=\beta\bar{\sigma}.\]

Such corrections are done prior to performing the plasticity calculation. Updating the damage variable, \(\alpha^{n+1}\), is done via relations (4.65) (or (4.66)) after convergence is achieved for the inelastic correction.

4.17.3. User Guide

BEGIN PARAMETERS FOR MODEL PLANE_STRESS_RATE_PLASTICITY
   #
   # Elastic constants
   #
   YOUNGS MODULUS = <real>
   POISSONS RATIO = <real>
   SHEAR MODULUS  = <real>
   BULK MODULUS   = <real>
   LAMBDA         = <real>
   TWO MU         = <real>
   #
   # Optional parameters related to inelastic correction criteria
   #
   TOLERANCE          = <real> tolerance (1.0e-10)
   MAX_INEL_CORR_ITER = <int> maximum_correction_iterations (100)
   MAX_LS_CORR_ITER   = <int> maximum_line_search_cutbacks  (20)
   #
   USER RATE DEPENDENCE = 0|1(0)
   YIELD STRESS         = <real>
   #
   FORMULATION          = <int> formulation (1)
   #
   # Input Options for USER RATE DEPENDENCE = 0
   #
   # linear hardening
   HARDENING MODULUS = <real> hardening_modulus
   # power law hardening
   HARDENING CONSTANT = <real> hardening_constant
   HARDENING EXPONENT = <real> hardening_exponent (0.5)
   # multilinear hardening
   HARDENING FUNCTION = <string>
   #
   # Rate dependence
   #
   # Johnson-Cook rate dependence
   RATE CONSTANT = <real>
   # multilinear rate dependence
   RATE FUNCTION = <string> rate_function_name
   #
   # Input Options for USER RATE DEPENDENCE = 1
   #
   # rate-dependent yield
   YIELD STRAIN RATES = <real_list> yield_strain_rates
   YIELD CURVES       = <string_list> yield_function_names
   # rate-dependent damage
   FRACTURE STRAIN RATES = <real_list> fracture_strain_rates
   FRACTURE CURVES       = <string_list> fracture_function_names
   DECAY COEFFICIENT     = <real> (1.0)
   #
   # Thermal softening commands (Johnson-Cook)
   INITIAL TEMPERATURE   = <real>
   MELT TEMPERATURE      = <real>
   REFERENCE TEMPERATURE = <real>
   THERMAL EXPONENT      = <real>
   #
END [PARAMETERS FOR MODEL PLANE_STRESS_RATE_PLASTICITY]

In the command blocks that define the Plane Stress Rate Plasticity model:

• TOLERANCE specifies the numerical value used for assessing convergence of the plastic correction routine.

• Optionally, the user can specify the maximum number of inelastic correction iterations for the plasticity inelastic correction routines. By default this value is 100.

• Optionally, the user can specify the maximum number of line search cutbacks in the plasticity correction routine. By default this value is 20 and only impacts the formulation not equal to 0 (plastic strain rate) case.

• The reference nominal yield stress, \(\bar{\sigma}\), is defined with the YIELD STRESS command line.

• The formulation parameter defines whether the total strain rate (formulation = 0) or equivalent plastic strain rate (anything else) is used for calculating rate dependence.

• INITIAL TEMPERATURE defines the initial temperature at \(t=0\).

• MELT TEMPERATURE defines \(T_{melt}\) in (4.64).

• REFERENCE TEMPERATURE defines \(T_{ref}\) in (4.64).

• THERMAL EXPONENT defines \(M\) in (4.64).

• The USER RATE DEPENDENCE is used to control the way hardening may be specified and whether or not failure is calculated.

  • For USER RATE DEPENDENCE = 0, plastic hardening may be specified as linear, power-law, OR multilinear. Failure cannot be used with USER RATE DEPENDENCE = 0:

    • For linear hardening, a non-zero HARDENING MODULUS should be specified. Do not give if using power-law or multilinear hardening.

    • For power-law hardening, the HARDENING CONSTANT should be specified. Optionally, the HARDENING EXPONENT parameter should be specified if the default value (0.5) is not to be used. Do not give for linear or multilinear hardening.

    • For multilinear hardening, a function name should be given for HARDENING FUNCTION command. Do not specify for linear or power-law hardening.

    • For Johnson-Cook rate-dependence, a rate constant must be specified via the RATE CONSTANT command. Do not use if using functionally specified rate-dependence.

    • For functionally defined rate-dependence, a function name should be given via the RATE FUNCTION command. Do not specify if using Johnson-Cook type rate dependence.

  • For USER RATE DEPENDENCE = 1, plastic hardening is specified through a series of user functions. Failure can be modeled

    • Rate-dependent plastic hardening is specified jointly via the YIELD STRAIN RATES and YIELD CURVES commands. YIELD STRAIN RATES is a list of strain rates corresponding one-to-one to functions specified in the YIELD CURVES list of strings giving user-defined function names. Each YIELD CURVES function should give the plastic isotropic plastic hardening curve at the corresponding rate given in YIELD STRAIN RATES.

    • Failure strains used in calculating damage are specified via the FRACTURE STRAIN RATES and FRACTURE CURVES commands. Similarly to the plastic hardening, a list of strain rates should be given with the FRACTURE STRAIN RATES input. Each rate should correspond one-to-one with a user function listed via FRACTURE CURVES. These functions are specified as a function of triaxiality and should give the failure strain at the specified rate.

    • The decay coefficient, \(C_1\), that controls the exponential decay of the yield stress related to failure should be specified via the DECAY COEFFICIENT command.

Output variables available for this model are listed in Table 4.25.

Table 4.25 State Variables for PLANE STRESS RATE PLASTICITY Model

Name	Description
EQPS	equivalent plastic strain, \(\bar{\varepsilon}^{p}\)
EQDOT	equivalent plastic strain rate, \(\dot{\bar{\varepsilon}}^{p}\)
SEFF	effective stress, \(\phi\)

---

[1] (1,2,3,4,5,6)

J. C. Simo and R. L. Taylor. A return mapping algorithm for plane stress elastoplasticity. International Journal for Numerical Methods in Engineering, 22:649–670, 1986.

[2] (1,2,3,4,5)

J.C. Simo and T.J.R. Hughes. Computational Inelasticity. Springer-Verlag, New York, NY, 1998.

[3] (1,2)

G.R. Johnson and W.H. Cook. A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures. In Proc. 7th. Int. Symp. on Ballistics, 541–547. The Hague, The Netherlands, 1983.

[4]

G.R. Johnson and W.H. Cook. Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures. Engineering Fracture Mechanics, 21(1):31–48, 1985.

[5]

W. M. Scherzinger. A return mapping algorithm for isotropic and anisotropic plasticity models using a line search method. Computer Methods in Applied Mechanics and Engineering, 317:526–553, 2017.