Porous metal plasticity

The porous metal plasticity model is intended for use with mildly voided metals. Even though the material that contains the voids (also known as the matrix material) is assumed to be plastically incompressible, the plastic behavior of the bulk material is pressure-dependent due to the presence of voids. The model is described, followed by a brief description of the material point calculations.

The following topics are discussed:

Yield condition
Flow rule
Evolution of the plastic strains and f
Integration of the elastoplastic equations
Computing the linearization moduli

Yield condition

For a metal containing a dilute concentration of voids, based on a rigid-plastic upper bound solution for spherically symmetric deformations of a single spherical void, Gurson (1977) proposed a yield condition of the form

ϕ = {(\frac{q}{σ_{y}})}^{2} + 2 q_{1} f \cosh (- \frac{3}{2} \frac{q_{2} p}{σ_{y}}) - (1 + q_{3} f^{2}) = 0,

where

S = p I + σ

is the deviatoric part of the macroscopic Cauchy stress tensor $σ$ ;

q = \sqrt{\frac{3}{2} S : S}

is the Mises stress;

p = - \frac{1}{3} σ : I

is the hydrostatic pressure; f is the volume fraction of the voids in the material; and $σ_{y} ({\bar{ε}}_{m}^{p l})$ is the yield stress of the fully dense matrix material as a function of ${\bar{ε}}_{m}^{p l}$ , the equivalent plastic strain in the matrix. Tvergaard (1981) introduced the constants $q_{1}$ , $q_{2}$ , and $q_{3} = q_{1}^{2}$ (as coefficients of the void volume fraction and pressure terms) to make the predictions of the Gurson model agree with numerical studies of “ordered” voided materials in plane strain tensile fields; one can recover the original Gurson model by setting $q_{1} = q_{2} = q_{3} = 1$ .

It should be noted that $f = 0$ implies that the material is fully dense, and the Gurson yield condition reduces to that of von Mises; $f = 1$ implies that the material is fully voided and has no stress carrying capacity. This is illustrated in Figure 1, where the yield surfaces for different levels of porosity are shown in the p–q plane.

Figure 1. Schematic of the yield surface in the p–q plane.

Figure 2 compares the behavior of a porous metal (which has an initial void volume fraction of $f_{0}$ ) in tension and compression against the perfectly plastic matrix material; the initial yield stress of the porous metal is $σ_{y_{0}}$ .

Figure 2. Schematic of uniaxial behavior of a porous metal.

In compression the porous material “hardens” due to closing of the voids, and in tension it “softens” due to growth and nucleation of the voids.

Flow rule

The plastic strains are derived from the yield potential; the presence of the first invariant of the stress tensor in the yield condition results in nondeviatoric plastic strains:

{\dot{ε}}^{p l} = \dot{λ} \frac{\partial ϕ}{\partial σ} = \dot{λ} (- \frac{1}{3} \frac{\partial ϕ}{\partial p} I + \frac{3}{2 q} \frac{\partial ϕ}{\partial q} S),

where $\dot{λ}$ is the nonnegative plastic flow multiplier.

Evolution of the plastic strains and f

The hardening of the (fully dense) matrix material is described through $σ_{y} = σ_{y} ({\bar{ε}}_{m}^{p l})$ . The evolution of ${\bar{ε}}_{m}^{p l}$ is assumed to be governed by the equivalent plastic work expression; i.e.,

(1 - f) σ_{y} {\dot{\bar{ε}}}_{m}^{p l} = σ : {\dot{ε}}^{p l} .

The change in volume fraction of the voids is due partly to the growth of existing voids and partly to the nucleation of voids. Growth of existing voids is based on the law of conservation of mass and is expressed as

{\dot{f}}_{gr} = (1 - f) {\dot{ε}}^{p l} : I .

Nucleation of voids can occur due to micro-cracking and/or decohesion of the particle-matrix interface. Abaqus assumes that the nucleation of new voids is plastic strain controlled (see Chu and Needleman, 1980), so that

{\dot{f}}_{nucl} = A {\dot{\bar{ε}}}_{m}^{p l},

where

A = \frac{f_{N}}{s_{N} \sqrt{2 π}} \exp [- \frac{1}{2} {(\frac{{\bar{ε}}_{m}^{p l} - ε_{N}}{s_{N}})}^{2}] .

The normal distribution of the nucleation strain has a mean value $ε_{N}$ , a standard deviation $s_{N}$ , and nucleates voids with volume fraction $f_{N}$ . The total rate of change of f is given as

\dot{f} = {\dot{f}}_{gr} + {\dot{f}}_{nucl} .

Voids are nucleated only in tension; Abaqus will not consider the nucleation term at a material point if the stress state is compressive.

The nucleation function $A / f_{N}$ , which is assumed to have a normal distribution, is shown in Figure 3 for different values of the parameter $s_{N}$ . Figure 4 shows the extent of softening in a uniaxial tension test of a porous material for different values of $f_{N}$ .

Figure 3. Nucleation function

A / f_{N}

Figure 4. Softening (in uniaxial tension) as a function of

f_{N}

Integration of the elastoplastic equations

The integration of the elastoplastic equations for the porous plasticity model is carried out using the backward Euler scheme proposed by Aravas (1987). This method is briefly discussed in the following paragraphs; the user can refer to the paper for further details.

During the constitutive calculations in an increment, the stress and state variables are known at time t (beginning of the increment). Given a total incremental strain $Δ ε$ , the stress and state variables need to be updated at $t + Δ t$ (end of the increment) so that they satisfy the yield condition, flow rules, and evolution equation of the state variables. To do this, consider the elasticity equations

\begin{matrix} (1) & σ = D^{e l} : ε^{e l} = D^{e l} : ({ε^{e l} |}_{t} + Δ ε^{e l}) = σ^{e l} - D^{e l} : Δ ε^{p l}, \end{matrix}

where

σ^{e l} = D^{e l} : ({ε^{e l} |}_{t} + Δ ε)

is the elastic predictor and

D^{e l} = 2 G ℑ + (K - \frac{2}{3} G) I I

is the linear isotropic elasticity tensor with G and K being the shear and bulk modulus, respectively, and $ℑ$ and $I$ being the fourth- and second-order identity tensors, respectively. Also, in the above, the additive decomposition of strain is used to write the total incremental strain as the sum of the elastic and plastic parts. All of the stress and state variables are evaluated at $t + Δ t$ unless indicated otherwise.

The yield condition, the flow rule, and the evolution of the state variables are rewritten as

ϕ (σ, H^{α}) = 0,

\begin{matrix} (2) & Δ ε^{p l} = \frac{1}{3} Δ ε_{p} I + Δ ε_{q} n, \end{matrix}

and

\begin{matrix} (3) & Δ H^{α} = {\bar{H}}^{α} (Δ ε^{p l}, σ, H^{β}), \end{matrix}

where

\begin{matrix} n & = \frac{3}{2 q} S, \\ Δ ε_{p} & = - Δ λ \frac{\partial ϕ}{\partial p}, \\ Δ ε_{q} & = Δ λ \frac{\partial ϕ}{\partial q} . \end{matrix}

In the above $H^{α}$ , $α = 1, 2$ are the state variables ${\bar{ε}}_{m}^{p l}$ and f, respectively. The plastic multiplier $Δ λ$ is eliminated from the last two equations to give

Δ ε_{p} \frac{\partial ϕ}{\partial q} + Δ ε_{q} \frac{\partial ϕ}{\partial p} = 0.

Equation 2 is used in the elasticity Equation 1 to yield

\begin{matrix} (4) & σ = σ^{e l} - K Δ ε_{p} I - 2 G Δ ε_{q} n . \end{matrix}

$S^{e l}$ and $n$ are coaxial; therefore, $n$ is determined as

\begin{matrix} (5) & n = \frac{3}{2 q^{e l}} S^{e l} . \end{matrix}

Once $n$ is known, it is easily seen that consistent determination of the scalars $Δ ε_{p}$ and $Δ ε_{q}$ completes the solution. Therefore, the problem of integrating the pressure-dependent elastoplastic constitutive equations is reduced to solving the following two nonlinear equations for the scalars $Δ ε_{p}$ and $Δ ε_{q}$ :

\begin{matrix} (6) & ϕ (p, q, H^{α}) = 0, \end{matrix}

\begin{matrix} (7) & Δ ε_{p} \frac{\partial ϕ}{\partial q} + Δ ε_{q} \frac{\partial ϕ}{\partial p} = 0. \end{matrix}

In the above equations p, q, and $H^{α}$ are defined by

\begin{matrix} (8) & p = p^{e l} + K Δ ε_{p}, \end{matrix}

\begin{matrix} (9) & q = q^{e l} - 3 G Δ ε_{q}, \end{matrix}

\begin{matrix} (10) & Δ H^{α} = h^{α} (Δ ε_{p}, Δ ε_{q}, p, q, H^{β}), \end{matrix}

where Equation 8 and Equation 9 are obtained by projecting Equation 4 onto $I$ and $n$ , respectively, and Equation 10 is an alternate form of Equation 3. Solving the above system of equations for the unknowns $p, q, Δ ε_{p}, Δ ε_{q},$ and $H^{α}$ completes the integration algorithm for the porous plasticity model.

Equation 6 and Equation 7 are solved for $Δ ε_{p}$ and $Δ ε_{q}$ using Newton's method; p and q are updated using Equation 8 and Equation 9; the state variables are updated using Equation 10.

Computing the linearization moduli

In the implicit finite element method of solving large-deformation problems, the discretized equilibrium equations result in a set of nonlinear equations for the nodal unknowns at the end of the increment. Abaqus/Standard uses Newton's method to solve these equations, which requires the computation of linearization moduli

J = \frac{\partial Δ σ}{\partial Δ ε} = {(\frac{\partial σ}{\partial ε})}_{t + Δ t} .

To compute $J$ (also known as the Jacobian), we start with the elasticity equations (Equation 4), which can be rewritten as

σ = 2 G ({e^{e l} |}_{t} + Δ e - Δ ε_{q} n) + K ({ε_{k k}^{e l} |}_{t} + Δ ε_{k k} - Δ ε_{p}) I,

where $e = ε - (1 / 3) ε_{k k} I$ is the deviatoric part of $ε$ . From the above equation you find that

\begin{matrix} (11) & \partial σ = 2 G (\partial e - \partial Δ ε_{q} n - Δ ε_{q} \partial n) + K (\partial ε_{k k} - \partial Δ ε_{p}) I . \end{matrix}

Once again Equation 6 and Equation 7 are used in computing the variations $\partial Δ ε_{p}$ and $\partial Δ ε_{q}$ .

Aravas (1987) After some lengthy algebraic calculations, a set of linear equations is obtained that can be solved for $\partial Δ ε_{p}$ and $\partial Δ ε_{q}$ . These derivatives are substituted into Equation 11 to obtain the linearization moduli. In general this linearization moduli is not symmetric. Further details of the derivation of the Jacobian can be found in Aravas (1987).