/FAIL/GURSON
Block Format Keyword A GursonNahshonHutchinson failure model describing the damage in terms of void nucleation and growth in metal plasticity.
The modified Gurson formulation adds additional damage accumulation terms for shear dominated loads, specific treatment under compressive loading, and elastic stiffness loss with damage.
Format
(1)  (2)  (3)  (4)  (5)  (6)  (7)  (8)  (9)  (10) 

/FAIL/GURSON/mat_ID/unit_ID  
${q}_{1}$  ${q}_{2}$  I_{loc}  
${\epsilon}_{n}$  A_{s}  K_{w}  
${f}_{c}$  ${f}_{F}$  ${f}_{I}$  
R_{len}  H_{chi}  ${L}_{e}^{MAX}$ 
(1)  (2)  (3)  (4)  (5)  (6)  (7)  (8)  (9)  (10) 

fail_ID 
Definition
Field  Contents  SI Unit Example 

mat_ID  Material identifier. (Integer, maximum 10 digits) 

unit_ID  (Optional) Unit Identifier. (Integer, maximum 10 digits) 

${q}_{1}$  First Gurson damage
coefficient. Default = 1.5 (Real) 

${q}_{2}$  Second Gurson damage coefficient,
maximum value = 1.02. Default = 1.0 (Real) 

I_{loc}  Damage variable accumulation method flag.
(Integer) 

${\epsilon}_{n}$  Equivalent plastic strain at void
nucleation. (Real) 

A_{s}  Linear void nucleation
slope. (Real) 

K_{w}  Shear damage growth
coefficient. (Real) 

${f}_{c}$  Critical void volume fraction at
void coalescence. (Real) 

${f}_{F}$  Void volume fraction at ductile
failure. (Real) 

${f}_{\mathrm{I}}$  Initial void volume
fraction. (Real) 

R_{len}  Radius of nonlocal variable
influence (I_{loc}
> 1). (Real) 
$\left[\text{m}\right]$ 
H_{chi}  Nonlocal penalty parameter
(Micromorphic method only,
I_{loc} =
2). (Real) 
$\left[\text{Pa}\right]$ 
${L}_{e}^{MAX}$  Mesh convergence element length
target. 5
(Real) 
$\left[\text{m}\right]$ 
fail_ID  (Optional) Failure criteria
identifier. (Integer, maximum 10 digits) 
Comments
 The Gurson damage model
may only be used with the elastoplastic material
/MAT/LAW104. The yield surface definition of the
material law is modified by adding the damage evolution terms:
$$\varphi =\frac{{\sigma}_{eq}^{2}}{{\sigma}_{yld}^{2}}1+2{q}_{1}{f}^{*}cosh\left(\frac{{\eta}_{t}{q}_{2}Tr\left(\sigma \right)}{2{\sigma}_{yld}}\right){\left({q}_{1}{f}^{*}\right)}^{2}=0$$Where,
 ${q}_{1}$ , ${q}_{2}$
 Two GursonTveergardNeedleman parameters.
 ${f}^{*}$
 Effective damage
 ${\eta}_{t}$
 Factor defined as:
 ${f}_{t}$
 Total void volume fraction that is computed incrementally.
The kinetic equations of the damage factor increments are: Void nucleation (creation of microcavities), decreasing at low
triaxiality.$$d{f}_{n}=\{\begin{array}{c}{A}_{s}d{\epsilon}_{p},\text{}{\epsilon}_{p}\ge {\epsilon}_{n}\text{and}\tau \ge 0\\ {A}_{s}\left(1+3\tau \right)d{\epsilon}_{p},\text{}{\epsilon}_{p}\ge {\epsilon}_{n}\text{and}\frac{1}{3}\le \tau 0\\ 0,\text{}{\epsilon}_{p}{\epsilon}_{n}\text{and}\tau \frac{1}{3}\end{array}$$Where, $\tau $ is the stress triaxiality defined as:$$\tau =\frac{Tr\left(\sigma \right)}{3{\sigma}_{eq}}$$
 Void growth at high triaxiality:$$d{f}_{g}=\left(1{f}_{t}\right)Tr\left({\epsilon}_{p}\right)$$
 Additional shear void growth at low triaxiality which is shear
dominated: $$d{f}_{sh}={K}_{w}{f}_{t}w\left(\theta \right)\frac{S:{\epsilon}_{p}}{{\sigma}_{eq}}$$Where, $w\left(\theta \right)$ is a weight function depending on the Lode angle:$$w\left(\theta \right)=1{\left(\text{cos}\left(3\theta \right)\right)}^{2}$$
To represent the cavities coalescence when a critical void volume fraction ${f}_{c}$ is reached by ${f}_{t}$ , the effective damage (which has an influence on the stress computation) ${f}^{*}$ is introduced in the model and its expression depends on ${f}_{t}$ :
$${f}^{*}=f\left({f}_{t}\right)=\{\begin{array}{c}{f}_{t},{f}_{t}{f}_{c}\\ {f}_{c}+\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{1ex}{${q}_{1}$}\right.{f}_{c}\right)\frac{\left({f}_{t}{f}_{c}\right)}{\left({f}_{R}{f}_{c}\right)},{f}_{t}\ge {f}_{c}\end{array}$$Where, ${f}_{R}$ is the total void volume fraction at rupture for which ${f}^{*}=1/{q}_{1}$ .To take into account the effect of the stiffness loss, a damage variable $D$ is computed as:
$$D={q}_{1}{f}^{*}$$The effective damage ${f}^{*}$ is normalized by its rupture value $1/{q}_{1}$ which gives $0\le D\le 1$ . The stress tensor is then computed as:
$$\sigma =\left(1D\right)\tilde{C}{\epsilon}_{e}$$Where, $\tilde{C}$ is the elastic stiffness matrix.
The material fails when the cumulated total damage factor reaches the limit value ${f}_{R}$ . The element is then deleted.
 By default
${I}_{loc}=1$
, the damage variable is calculated step by
step using the local plastic strain values at each integration point.
However, one may want to use nonlocal regularization which offers mesh size
and the mesh orientation independent results (mesh convergence) for all
meshes using the mesh size
${L}_{e}$
less than equal to the maximum value set by
you
${L}_{e}\le {L}_{e}^{max}$
. This maximum mesh size
${L}_{e}^{max}$
is then the highest mesh size used for which
results are mesh convergent.
If one of the nonlocal formulations is used, $({I}_{loc}>1)$ , the damage increments depend on a regularized nodal “nonlocal” plastic strain calculated on the entire mesh. The nonlocal plastic strain at nodes denoted ${\epsilon}_{p}^{nl}$ is computed accounting for its own gradient and its local counterpart ${\epsilon}_{p}$ computed at the Gauss points following the set of equations:
$$\begin{array}{c}{R}_{len}^{2}\text{\Delta}{\epsilon}_{p}^{nl}\gamma \dot{{\epsilon}_{p}^{nl}}+\left({\epsilon}_{p}{\epsilon}_{p}^{nl}\right)=\zeta \ddot{{\epsilon}_{p}^{nl}}\\ \overrightarrow{\nabla}{\epsilon}_{p}^{nl}.\overrightarrow{n}=0\end{array}\begin{array}{c}on\\ on\end{array}\begin{array}{c}\Omega \\ \Gamma \end{array}$$The parameters $\gamma $ and $\zeta $ are automatically set. You have to set the parameter R_{len} (or ${L}_{e}^{MAX}$  Comment 5) which defines a nonlocal “internal length” which corresponds to a radius of influence in the nonlocal variable computation. This defines the size of the nonlocal regularization band ${L}_{r}=f\left({R}_{len}\right)$ (Figure 5).To help choose a value for the parameter R_{len}, one may follow the following expression:
$${R}_{len}\approx \frac{3{L}_{e}^{max}}{\sqrt{\pi}}$$  If
${I}_{loc}=2$
, the nonlocal Micromorphic method will be
used. For this method, the parameter is required,
H_{chi}. This
parameter and the nonlocal plastic strain
${\epsilon}_{p}^{nl}$
are introduced in the constitutive equation
as: $${R}_{chi}\left({\epsilon}_{p},{\epsilon}_{p}^{nl}\right)=R\left({\epsilon}_{p}\right){H}_{chi}\left({\epsilon}_{p}^{nl}{\epsilon}_{p}\right)$$
Where, $R\left({\epsilon}_{p}\right)$ is the classical workhardening function. This newly defined micromorphic workhardening function R_{chi} is introduced in the flow stress computation ${\sigma}_{yld}$ . The parameter H_{chi} becomes a penalty parameter and if ${H}_{chi}\to \infty $ , then ${\epsilon}_{p}\to {\epsilon}_{p}^{nl}$ and ${\epsilon}_{p}^{nl}\to {\epsilon}_{p}$ and so ${\epsilon}_{p}\approx {\epsilon}_{p}^{nl}$ . This method is thermodynamically well defined. However, it is hard to identify the input values and it changes the plastic behavior of the model. This is why it is recommended to use the Peerlings method ${I}_{loc}=3$ .
 If
${I}_{loc}=3$
, the nonlocal Peerlings method will be
used. For this method, the parameter
H_{chi} is used. Only
the nonlocal length R_{len}
is used. This method is simpler than the Micromorphic one. It introduces the
nonlocal plastic strain in the softening variable kinetic equation (damage
and temperature if thermal effects are considered): $${\dot{f}}_{t}=\underset{\text{Voidnucleation}}{\underbrace{A\dot{{\epsilon}_{p}^{nl}}}}+\underset{\begin{array}{c}\text{Voidgrowth}\\ \left(\text{hightriaxiality}\right)\end{array}}{\underbrace{\left(1{f}_{t}\right)Tr\left(\dot{{\epsilon}_{p}^{nl}}\right)}}\text{+}\underset{\begin{array}{l}\text{Shearnucleation}\\ \text{(lowtriaxility)}\end{array}}{\underbrace{{K}_{w}{f}_{t}w\left(\theta \right)\frac{S:\dot{{\epsilon}_{p}^{nl}}}{{\sigma}_{eq}}}}$$$$\dot{T}=\omega \left(\dot{{\epsilon}_{p}^{nl}}\right)\frac{\eta}{\rho {C}_{p}}\overline{\overline{\sigma}}:\dot{{\epsilon}_{p}^{nl}}$$
This method is recommended since it is simple to identify the input parameters and does not modify the plastic behavior of the material.
 To set the nonlocal
length parameter R_{len}, you
can choose:
 Directly input the value of R_{len} in the input card, if a direct control on this parameter is wanted. In this case, the parameter ${L}_{e}^{MAX}$ must be ignored and set to none.
 Input the maximum mesh size ${L}_{e}^{MAX}$ for which results have reached mesh convergence. The nonlocal regularization will then be effective for all mesh sizes ${L}_{e}^{}$ such as ${L}_{e}\le {L}_{e}^{max}$ . In this case R_{len} is calculated automatically according to the value of ${L}_{e}^{MAX}$ , and the input value of R_{len} is ignored. For instance, if you want to get converged and meshindependent results for a mesh size of 5 mm, ${L}_{e}^{max}=5$ mm. In this case, the results will be converged, meshsize and mesh orientation independent for ${L}_{e}\le 5$ mm.
 When a nonlocal regularization is used for shell elements, an
additional regularization is made on the thickness variation computation
avoiding an additional localization issue. In the common local case (Figure 6), the compatibility of thickness
between shell elements is not ensured due to the lack of kinematic equations
in the zdirection, and the thickness variation is locally computed at Gauss
points. By introducing the nonlocal plastic strain in the “inthickness”
strain increment, the compatibility is restored, (Figure 7). $$\text{\Delta}{\epsilon}_{zz}=\frac{\nu}{1\nu}\left(\text{\Delta}{\epsilon}_{xx}\text{\Delta}{\lambda}_{nl}{n}_{xx}+\text{\Delta}{\epsilon}_{yy}\text{\Delta}{\lambda}_{nl}{n}_{yy}\right)+\text{\Delta}{\lambda}_{nl}{n}_{zz}$$Where, $\text{\Delta}{\lambda}_{nl}=f\left({\epsilon}_{p}^{nl}\right)$ is the nonlocal plastic multiplier.Note: This last point implies that the identified parameters can be used on solid and shells, as results will be identical within the same range of stress triaxiality $\raisebox{1ex}{$2$}\!\left/ \!\raisebox{1ex}{$3$}\right.\le \eta \le \raisebox{1ex}{$2$}\!\left/ \!\raisebox{1ex}{$3$}\right.$ .
 To create a specific damage output DAMA in ANIM and H3D files, the
total damage is normalized by its rupture value: $$D=\frac{{f}_{t}}{{f}_{R}}$$