# Reactive Plasticity Benchmark Problem

## Overview

The objective of this case study is to present the FEBio analysis of a classical finite plasticity benchmark problem, using the reactive plasticity framework introduced in FEBio 3.5, based on the paper by Zimmerman et al. (2021) ^{1}Zimmerman BK, Jiang D, Weiss JA, Timmins LH, Ateshian GA. On the use of constrained reactive mixtures of solids to model finite deformation isothermal elastoplasticity and elastoplastic damage mechanics. Journal of the Mechanics and Physics of Solids. 2021 Oct 1;155:104534.. The square-cup deep-drawing benchmark problem was developed by the NUMISHEET’93 organizing committee and experimental results were presented at a conference in Japan in September 1993. Details of the experimental results presented in this case study were obtained from the paper by Joachim Danckert ^{2}Danckert J. Experimental investigation of a square-cup deep-drawing process. Journal of Materials Processing Technology. 1995 Mar 1;50(1-4):375-84..

## Reactive Plasticity in FEBio

Starting with FEBio3.5 we have introduced a finite plasticity framework in FEBio, based on constrained reactive mixture theory. This reactive plasticity framework is consistent with classical plasticity, but the presentation of the governing equations differs sufficiently to require a summary presentation, as done here. We start by presenting the concept of reactive elastic-perfect plasticity, where the stress in uniaxial loading of a bar initially increases linearly with strain, then remains constant with increasing strain after yielding. Then, we generalize it to encompass the classical concept of strain hardening.

### Reactive Elastic-Perfect Plasticity Response

Consider that a material consists of molecules that are bonded together and that these bonds store the elastic energy produced by loading the material. Below a certain threshold of loading we may assume that the bonds remain intact as they store more energy. We may label these intact bonds $\alpha=s$ and consider that $\mathbf{X}^{s}$ represents the reference configuration of the molecules involved with these bonds. However, when the load exceeds a certain threshold we assume that the bonds break and (possibly) reform into a new state. If the bonds remain permanently broken, we consider this reaction to model damage mechanics; if the bonds reform in a stress-free state, the reaction models viscoelasticity; if the bonds reform in a stressed state, the reaction models plasticity. In this presentation we only focus on the latter mechanism, in which case the bond-breaking-and-reforming reaction represents *yielding* behavior.

When bonds associated with $\alpha=s$ break and reform into a new stressed state at time $t=u$, we may label this new generation of bonds $\alpha=u$ and consider that the reference configuration of molecules associated with bonds $\alpha=u$ may be represented by $\mathbf{X}^{u}$. The bond-breaking-and-reforming reaction may be denoted by $\mathcal{E}^{s}\to\mathcal{E}^{u}$. At a subsequent time $t=v$, if the loading once again exceeds the ability of bonds $u$ to sustain the applied load, they may break and reform into a new bond generation $\alpha=v$, via the yielding reaction $\mathcal{E}^{u}\to\mathcal{E}^{v}$. Thus, the lifetime of $u-$generation bonds is $u\le t<v$. For an elastic-perfectly plastic response, we assume that all the bonds $u$ present in an elemental region break and reform simultaneously into a new generation $v$. Therefore, only one bond generation is extant at any given time $t$.

The deformation gradient of intact bonds is given by $\mathbf{F}^{s}=∂\boldsymbol{\chi}/∂\mathbf{X}^{s}$, where $\boldsymbol{\chi}\left(\mathbf{X}^{s},t\right)$ is the motion of the continuum. In general, $\mathbf{X}^{\alpha}$ may represent the reference configuration of any bond generation $\alpha$. Thus, using the chain rule of differentiation, we may write

$\mathbf{F}^{s}\left(\mathbf{X}^{s},t\right)=\frac{∂\boldsymbol{\chi}}{∂\mathbf{X}^{s}}=\frac{∂\boldsymbol{\chi}}{∂\mathbf{X}^{\alpha}}\cdot\frac{∂\mathbf{X}^{\alpha}}{∂\mathbf{X}^{s}}=\mathbf{F}^{\alpha}\left(\mathbf{X}^{s},t\right)\cdot\mathbf{F}^{\alpha s}\left(\mathbf{X}^{s}\right)\,.\label{eq:kin-chain-rule}$

In this expression, $\mathbf{F}^{\alpha}$ is the deformation gradient of generation $\alpha$ relative to its reference configuration $\mathbf{X}^{\alpha}$, whereas $\mathbf{F}^{\alpha s}=∂\mathbf{X}^{\alpha}/∂\mathbf{X}^{s}$ represents the mapping between $\mathbf{X}^{\alpha}$ and $\mathbf{X}^{s}$. Note that $\mathbf{F}^{\alpha s}$ is not a function of time; the mapping between $\alpha$ and $s$ is determined at the time $t=\alpha$ when generation $\alpha$ is produced in a bond-breaking-and-reforming reaction. In our constrained reactive mixture framework, the mapping $\mathbf{F}^{\alpha s}\left(\mathbf{X}^{s}\right)$ is a function of state for which a constitutive relation must be provided. It follows that $\mathbf{X}^{\alpha}$, which may be calculated by suitably integrating $\mathbf{F}^{\alpha s}$, is not an observable variable (therefore we don’t need to calculate it!). The constitutive relation adopted for $\mathbf{F}^{\alpha s}\left(\mathbf{X}^{s}\right)$ in FEBio is provided below, in the section titled Flow Rule.

### Yield Measure

We may define a scalar yield measure $\Phi\left(\mathbf{F}^{\alpha}\right)$ (such as the von Mises stress) which determines when bonds of generation $\alpha$ will yield. In particular, we may define the *yield surface*

$\varphi\left(\mathbf{F}^{\alpha}\right)=\Phi\left(\mathbf{F}^{\alpha}\right)-\Phi_{m}\label{eq:yield-function}$

where $\Phi_{m}$ is the threshold that must be exceeded in order for generation $\alpha$ to break and reform into a new generation. Let $d\varphi$ represent the change in $\varphi$ over consecutive time increments. Yielding of generation $u$ occurs at time $v$ when $\varphi\left(\mathbf{F}^{u}\left(u\right)\right)=0$ and $\varphi\left(\mathbf{F}^{u}\left(v\right)\right)>0$, producing $d\varphi=\varphi\left(\mathbf{F}^{u}\left(v\right)\right)-\varphi\left(\mathbf{F}^{u}\left(u\right)\right)>0$, at which point generation $u$ breaks and reforms into the new generation $v$ to restore

$d\varphi=\varphi\left(\mathbf{F}^{v}\left(v\right)\right)-\varphi\left(\mathbf{F}^{u}\left(u\right)\right)=0\,.\label{eq:persistency-condition}$

Here, $d\varphi=0$ is known as the persistency condition which must be satisfied by consecutive yielding generations.

### Stress

The stress $\boldsymbol{\sigma}$ in an elastic-perfectly plastic material is evaluated from

$\boldsymbol{\sigma}=\sum_{\alpha}w^{\alpha}\mathbf{T}\left(\mathbf{F}^{\alpha}\right)\label{eq:epp-stress}$

where $\mathbf{T}\left(\mathbf{F}\right)$ is the constitutive model for the elastic response of the material, and $w^{\alpha}$ is the mass fraction of generation $\alpha$ in the mixture. As stated above, only one generation $\alpha$ is extant at any time $t$ in an elastic-perfectly plastic material. For example, during the lifetime $u\le t<v$ of generation $u$, we have $w^{u}=1$ and $w^{\alpha}=0$ for $\alpha\ne u$. Therefore, during that generation, the above expression for the stress reduces to $\boldsymbol{\sigma}=\mathbf{T}\left(\mathbf{F}^{u}\right)$ (Figure 1).

### Flow Rule

The constitutive model for the mapping $\mathbf{F}^{\alpha s}$ of generation $\alpha$ is postulated in a manner that remains consistent with classical finite plasticity. Thus, for historical reasons, we may describe this constitutive model as a flow rule, even though our reactive plasticity framework assumes that $\mathbf{F}^{\alpha s}$ is strictly not a function of time.

First, using the polar decomposition theorem, we evaluate the relative deformation gradient of $\alpha$ as $\mathbf{F}^{\alpha}=\mathbf{R}\cdot\mathbf{U}^{\alpha}$, where $\mathbf{U}^{\alpha}$ is the right stretch tensor of generation $\alpha$, and $\mathbf{R}$ is the rotation tensor, which is assumed to be the same for all generations, thus $\mathbf{F}^{s}=\mathbf{R}\cdot\mathbf{U}^{s}$. In the FEBio finite element implementation recall that the mesh deformation is calculated using $\mathbf{F}^{s}$, thus $\mathbf{R}$ may be extracted from $\mathbf{F}^{s}$ using standard methods.

Now, the flow rule adopted in FEBio for the reaction $\mathcal{E}^{u}\to\mathcal{E}^{v}$ is the recursive relation

$\left(\mathbf{F}^{vs}\right)^{-1}=\xi\left(\mathbf{F}^{us}\right)^{-1}\cdot\left(\mathbf{I}-\lambda\hat{\mathbf{N}}^{v}\right)\,,\label{eq:flow-rule}$

where

$\hat{\mathbf{N}}^{\alpha}=\frac{\mathbf{N}^{\alpha}}{\sqrt{\mathbf{N}^{\alpha}:\mathbf{N}^{\alpha}}}\,,\quad\mathbf{N}^{\alpha}=\frac{∂\varphi}{∂\mathbf{U}^{\alpha}}\,.$

Here, $\mathbf{N}^{\alpha}$ is the tensorial normal to the yield surface and $\hat{\mathbf{N}}^{\alpha}$ is its normalized value. The scalar $\lambda$ is obtained by satisfying the persistency condition, namely $\varphi\left(\mathbf{U}^{v}\right)=0$ at time $v$. The scalar $\xi$ is used to optionally enforce isochoric plastic flow, or equivalently $\det\mathbf{F}^{vs}=\det\mathbf{F}^{us}=1$. Thus,

$\xi=\begin{cases}\left[\det\left(\mathbf{I}-\lambda\hat{\mathbf{N}}^{v}\right)\right]^{-1/3} & \text{isochoric plastic flow} 1 & \text{otherwise}\end{cases}\,.\label{eq:flow-condition}$

For the earliest yielded generation $u$, the preceding $s-$generation is in the elastic regime; therefore, $\mathbf{F}^{us}=\mathbf{I}$ at the start of the above recursive relation.

## Kinematic “Hardening”

Real materials exhibit a “hardening” response following the initial yielding, with the stress subsequently increasing nonlinearly with increasing strain. In the reactive plasticity framework of FEBio, it is assumed that a multitude of bond families $\beta$ maintain the cohesiveness of molecular species within an elemental volume, with each bond family having a different yield threshold $\Phi_{m\beta}$. As the first family yields at $\Phi_{m0}$, other ones remain intact until the next threshold $\Phi_{m1}$ is reached, causing the second bond family to yield, and so on until all bond families have yielded, at which point a perfectly plastic response may be achieved. Optionally, it is possible to include a bond family that never yields (i.e., remains intact), causing a perpetual hardening behavior instead of the constant stress of a perfectly plastic response. In other words, the “hardening” response is simulated by causing bond families to yield in sequence, instead of all bonds yielding at once.

In FEBio it is assumed that there are $n_{f}$ yielding bond families, with $\beta=0,\ldots,n_{f}-1$, where $n_{f}$ is a user-defined parameter. The mass fraction of each of these bond families is $w_{\beta}$. If we choose to include an additional bond family $\beta=n_{f}$ that never yields, we denote its mass fraction by $w_{e}$. These mass fractions necessarily satisfy

$\sum_{\beta=0}^{n_{f}}w_{\beta}=1\,,\label{eq:bond-family-saturation}$

where the summation includes $w_{e}$. The bond mass fractions $w_\beta$ and thresholds $\Phi_{m\beta}$ are calculated internally based on user-defined parameters, including the lowest yield measure $\Upsilon_{0}=\Phi_{m0}$ (e.g., the yield stress), the ultimate yield measure $\Upsilon_{\text{max}}$ (e.g., the ultimate stress), and the number of bond families $n_f$. The stress in a hardening material is evaluated from

$\boldsymbol{\sigma}=\sum_{\beta=0}^{n_{f}}w_{\beta}\sum_{\alpha}w_{\beta}^{\alpha}\mathbf{T}\left(\mathbf{F}_{\beta}^{\alpha}\right)\,,\label{eq:hardening-stress}$

where $\mathbf{T}\left(\mathbf{F}\right)$ is the constitutive model for the elastic response of the material and $w_{\beta}^{\alpha}$ is the bond mass fraction of generation $\alpha$ in bond family $\beta$. As before, only one generation $\alpha$ is extant at any time during the response of bond family $\beta$. Thus, during the time interval $u\le t<v$ of the $u-$generation of family $\beta$, we have $w_{\beta}^{u}=1$ and $w_{\beta}^{\alpha}=0$ for $\alpha\ne u$ (Figure 2). Similarly, $\mathbf{F}_{\beta}^{\alpha}$ is the deformation gradient of generation $\alpha$ of bond family $\beta$, relative to its reference configuration $\mathbf{X}_{\beta}^{\alpha}$, which satisfies $\mathbf{F}^{s}\left(\mathbf{X}^{s},t\right)=\mathbf{F}_{\beta}^{\alpha}\left(\mathbf{X}^{s},t\right)\cdot\mathbf{F}_{\beta}^{\alpha s}\left(\mathbf{X}^{s}\right)$. Here again, $\mathbf{F}_{\beta}^{\alpha s}\left(\mathbf{X}^{s}\right)$ is a function of state describing the “flow rule” for yielding bond family $\beta$. In FEBio we assume that all yielding bond families obey the same flow rule, as given above.

## Deep Drawing of a Square Cup

This benchmark problem concerns the deep drawing of a cup from a sheet metal blank, an important application in many metal forming processes. In this example, a mild steel blank is sandwiched between a die and a holder and kept in place by application of a 19.6 kN holding force while a punch draws the material. Friction between all components is modeled as well.

The specific geometry used for this validation comes from the experimental benchmark devised by the organizers of the NUMISHEET 1993 international conference^{3}Makinouchi A, Nakamachi E, Onate E, Wagoner RH. Numerical simulation of 3-D sheet metal forming processes. Verification of simulation with experiments, NUMISHEET. 1993;93. and is shown in Figure 3; a discussion of the experimental procedure may be found in Danckert (1995)^{2}Danckert J. Experimental investigation of a square-cup deep-drawing process. Journal of Materials Processing Technology. 1995 Mar 1;50(1-4):375-84.. This example represents a challenging finite strain, finite rotation elastoplastic frictional contact problem with damage involved and thus serves to thoroughly validate our novel theory.

### Material Properties of Mild Steel

Material properties of the mild steel used in this validation were $E=206\,\text{GPa}$ and $\nu=0.3$ according to Table 1 of Danckert (1995)^{2}Danckert J. Experimental investigation of a square-cup deep-drawing process. Journal of Materials Processing Technology. 1995 Mar 1;50(1-4):375-84.; the mean value of the yield stress was given as $\Upsilon_0=173.1\,\text{MPa}$. This study also reported an equation for the mean true stress-true plastic strain response of the steel, $\sigma=565.32\left(0.007117+\eps_p\right)^{0.2589}\,\text{MPa}$, where $\sigma$ is the true stress and $\eps_p$ is the true plastic strain. According to this formula, $\sigma=157.1\,\text{MPa}$ when $\eps_p=0$, which is inconsistent with the given value of $\Upsilon_0$. Therefore we first adjust the formula to produce the correct value of $\Upsilon_0$: $\sigma=565.32\left(0.010344+\eps_p\right)^{0.2589}\,\text{MPa}$. Then, since FEBio uses the total strain $\eps = \eps_0 +\eps_p$ instead of the plastic strain when specifying the flow curve, we evaluate the yield strain $\eps_0 = \Upsilon_0/E = 173.1/206\times 10^3$, or $\eps_0=8.4029\times 10^{-4}$. Now, the stress-strain response reduces to $\sigma=565.32\left(0.0095037+\eps\right)^{0.2589}\,\text{MPa}$.

In FEBio we use the `"reactive plasticity"`

material to represent this mild steel, with isochoric plastic flow, `"natural neo-Hookean"`

isotropic hyperelastic response, and a yield criterion based on the von Mises stress. The flow curve is represented using a mathematical expression:

```
<material id="1" name="mild steel" type="reactive plasticity">
<density>1</density>
<isochoric>1</isochoric>
<elastic type="natural neo-Hookean">
<E>206000</E>
<v>0.3</v>
</elastic>
<yield_criterion type="DC von Mises stress"/>
<flow_curve type="PFC math">
<nf>15</nf>
<e0>8.4029e-4</e0>
<emax>1.5</emax>
<plastic_response>565.32*(0.0095037+eps)^0.2589</plastic_response>
</flow_curve?
</material>
```

The flow curve employed here is of type `"PFC math"`

. It requires the user to specify the number of bond families, `nf`

($n_f$), the yield strain `e0`

($\eps_0$), the maximum strain `emax`

expected in an analysis, and the formula for the plastic flow response, which was evaluated above. Internally, this formula is used to calcuate the bond mass fractions $w_\beta$, using a linear approximation for the hyperelastic response of the intact (non-yielded) material to evaluate $\Phi_{m\beta}=\Phi\left(\eps_\beta\right)=E\,\eps_\beta$, where $\eps_\beta$ is the strain at which bond family $\beta$ yields. However, when the material undergoes finite deformations, the actual flow curve may deviate from the prescribed formula, due to the nonlinearity of the hyperelastic response in the large strain regime. The `"natural neo-Hookean"`

material comes closest to approximating the linear hyperelastic response $\Phi\left(\eps\right)=E\,\eps$, but deviations of the plastic flow curve from the requested formula may persist even with this choice of hyperelastic material. As a result, it may be beneficial to run a separate reactive plasticity analysis of this material on a single finite element representing a unit cube, loaded uniaxially in tension. This analysis is performed in the file named `mild_steel_stress_strain_math.feb`

(available from the FEBioStudio Repository), for the strain range $0 \le \eps \le 1.25$ (equivalent to prescribing a displacement up to 2.5 on one of the cube faces, or equivalently a stretch ratio of 3.5, where $\ln{3.5} = 1.25$). A plot of $\sigma$ versus $\eps$ from the results of this analysis is presented in Figure 4 (blue), along with the desired flow curve (gray). A small discrepancy can be seen, where the results of the reactive plasticity analysis slightly overestimate the desired flow curve. In this example, a simple scaling of the flow curve formula may be used to adjust this deviation: Based on the final values of both curves, a scale factor of $0.96485$ is used on the coefficient $565.32$, whereas the offset strain value is adjusted to properly produce the desired $\Upsilon_0$ when $\eps=\eps_0$. The adjusted formula becomes `<plastic_response>545.46*(0.011024+eps)^0.2589</plastic_response>`

, and the resulting flow curve predicted by the finite element analysis (filename `mild_steel_stress_strain_final.feb`

) becomes indistinguishable from the desired response in Figure 4.

### FEA of Deep Drawing of Steel Plate

The finite element model file for this problem (`CupDrawingNUMISHEET93.fsm`

) can be found on the FEBioStudio repository. The finite element meshes of all components are shown in Figure 5, where only one-quarter of the problem was considered due to symmetry. The punch, die, and blank holder are all modeled as rigid bodies, and both the punch and die are represented by shell elements for simplicity. Due to the aspect ratio and significant rotation undergone by the blank, the blank is modeled with 625 isoparametric quadratic (20-node) hexahedral brick elements, as a convergence study revealed mesh locking with eight-node hexahedral elements. All surfaces were initially in perfect (grazing) contact, except the blank holder was slightly displaced into the mesh of the blank. This initial overlap was required for the load-control constraint on the blank holder to be computationally enforced. Our recently-developed robust surface-to-surface frictional contact algorithm ^{4}Zimmerman, B. K., and Ateshian, G. A. (July 2, 2018). “A Surface-to-Surface Finite Element Algorithm for Large Deformation Frictional Contact in febio.” ASME. *J Biomech Eng*. August 2018; 140(8): 081013. was employed between all contacting surfaces, with the friction coefficient $\mu=0.144$ given in the original problem specification^{2}Danckert J. Experimental investigation of a square-cup deep-drawing process. Journal of Materials Processing Technology. 1995 Mar 1;50(1-4):375-84.^{3}Makinouchi A, Nakamachi E, Onate E, Wagoner RH. Numerical simulation of 3-D sheet metal forming processes. Verification of simulation with experiments, NUMISHEET. 1993;93.. The contact penalty scale factor was chosen to be $\alpha_{c}=0.01$, with the auto-penalty calculation turned on.

**Figure 5.** Finite element meshes of all components for the deep drawing of a square cup. Dimensions are in mm.

A punch displacement of $u_{z}=-40$ mm was prescribed in 300 time steps, and the punch was then returned to the initial position to evaluate the elastic springback.

**Figure 6.** Development of effective plastic strain during deep drawing of a square cup. One-quarter of the cup was modeled; the images above use reflection planes in the post-processor.