Hybrid incompressible solid element formulation

ProductsAbaqus/StandardAbaqus/Explicit

Many problems involve the prediction of the response of almost incompressible materials. This is especially true at large strains, since most solid materials show relatively incompressible behavior under large deformations. In this section we describe the augmented virtual work basis provided in Abaqus/Standard for such cases. The method is described in the context of incompressible elasticity theory, since that is where it is most likely to be used.

When the material response is incompressible, the solution to a problem cannot be obtained in terms of the displacement history only, since a purely hydrostatic pressure can be added without changing the displacements. The nearly incompressible case (that is, when the bulk modulus is much larger than the shear modulus or Poisson's ratio, $ν$ , is greater than 0.4999999) exhibits behavior approaching this limit, in that a very small change in displacement produces extremely large changes in pressure, so that a purely displacement-based solution is too sensitive to be useful numerically (for example, round-off on the computer may cause the method to fail). We remove this singular behavior in the system by treating the pressure stress as an independently interpolated basic solution variable, coupled to the displacement solution through the constitutive theory and the compatibility condition, with this coupling implemented by a Lagrange multiplier. This independent interpolation of pressure stress is the basis of these “hybrid” elements. More precisely, they are “mixed formulation” elements, using a mixture of displacement and stress variables with an augmented variational principle to approximate the equilibrium equations and compatibility conditions. The hybrid elements also remedy the problem of volume strain “locking,” which can occur at much lower values of $ν$ (i.e., $ν = $ 0.49). Volume strain locking occurs if the finite element mesh cannot properly represent incompressible deformations. Volume strain locking can be avoided in regular displacement elements by fully or selectively reduced integration, as described in Solid isoparametric quadrilaterals and hexahedra.

We begin by writing the internal virtual work:

\begin{matrix} (1) & δ W = \int_{V} σ : δ ε d V, \end{matrix}

where $δ ε$ is the virtual strain:

δ ε = sym (\frac{\partial δ u}{\partial x}),

where $δ u$ is the virtual displacement field, $σ$ is the true (Cauchy) stress, V is the current volume, and $δ W$ is the virtual work as defined by this equation. See Equilibrium and virtual work for a detailed discussion of the virtual work concept.

In a displacement-based formulation the Cauchy stress, $σ$ , is obtained with the constitutive equations from the deformation, usually in rate form:

\begin{matrix} (2) & d σ = C : d ε + d W \cdot σ - σ \cdot d W, \end{matrix}

where $C$ is the “material stiffness matrix” and $d W$ is the rate of rotation (spin) of the material.

We modify the Cauchy stress by introducing an independent hydrostatic pressure field $\hat{p}$ as follows:

\begin{matrix} (3) & \bar{σ} = σ + (1 - ρ) I (p - \hat{p}), \end{matrix}

where

p = - \frac{1}{3} trace (σ)

is the hydrostatic pressure stress and $ρ$ is a small number. If $ρ$ was set equal to zero, the hydrostatic component in $\bar{σ}$ would be identical to the independent pressure field $\hat{p}$ , corresponding to a pure “mixed” formulation. The small nonzero value ( $10^{- 9}$ ) is chosen to avoid equation solver difficulties. This relation is used in incremental form:

\begin{matrix} (4) & \bar{σ} = {\bar{σ}}^{0} + Δ σ + (1 - ρ) I (Δ p - Δ \hat{p}), \end{matrix}

where ${\bar{σ}}^{0}$ is the modified Cauchy stress at the start of the increment. We use the modified Cauchy stress in the virtual work expression and augment the expression with the Lagrange multiplier enforced constraint $Δ p - Δ \hat{p} = 0$ :

\begin{matrix} (5) & δ \bar{W} = \int_{V} [\bar{σ} : δ ε + J^{- 1} δ λ (Δ p - Δ \hat{p})] d V, \end{matrix}

with J the volume change ratio (Jacobian) and $δ λ$ a Lagrange multiplier whose interpolation must still be determined. $Δ \hat{p}$ will be interpolated over each element so that the constraint is satisfied in an integrated (average) sense. Since $Δ p$ is the value of the equivalent pressure stress increment computed from the kinematic solution, Equation 4 does not make sense if the material is fully incompressible because then $Δ p$ cannot be computed. For the purpose of development we regard the bulk modulus as finite, and we will be able to show that the final formulation approaches a usable limit as we allow the bulk modulus to approach infinity.

For the formulation of the tangent stiffness (the Jacobian), we need to define the rate of change of $δ \bar{W}$ . Therefore, we rewrite the virtual work equation in terms of the reference volume $V^{0}$ :

\begin{matrix} (6) & δ \bar{W} = \int_{V^{0}} [J \bar{σ} : δ ε + δ λ (Δ p - Δ \hat{p})] d V^{0} . \end{matrix}

The rate of change $d δ \bar{W}$ is then readily obtained as

\begin{matrix} (7) & d δ \bar{W} = \int_{V^{0}} [d J \bar{σ} : δ ε + J d \bar{σ} : δ ε + J \bar{σ} : d δ ε + δ λ (d p - d \hat{p}) + d δ λ (Δ p - Δ \hat{p})] d V^{0} . \end{matrix}

We rewrite this expression in terms of the current volume:

\begin{matrix} (8) & d δ \bar{W} = \int_{V} [d \bar{σ} : δ ε + d ε : I \bar{σ} : δ ε + \bar{σ} : d δ ε + J^{- 1} δ λ (d p - d \hat{p}) + J^{- 1} d δ λ (Δ p - Δ \hat{p})] d V, \end{matrix}

where we used the identity $J^{- 1} d J = I : d ε$ .

The rate of the modified stress follows from Equation 4 and the constitutive equations:

\begin{matrix} (9) & d \bar{σ} = C : d ε + (1 - ρ) I (d p - d \hat{p}) + d W \cdot \bar{σ} - \bar{σ} \cdot d W, \end{matrix}

where

d p = - \frac{1}{3} trace (d σ) = - \frac{1}{3} I : C : d ε;

and we used the fact that $d W \cdot \bar{σ} - \bar{σ} \cdot d W = d W \cdot σ - σ \cdot d W$ , since $\bar{σ}$ and $σ$ differ only in the hydrostatic part. Substituting these expressions into the expression for the rate of virtual work yields

\begin{matrix} (10) & \begin{matrix} d δ \bar{W} = \int_{V} & {δ ε : C : d ε + δ ε : \bar{σ} I : d ε - \frac{1}{3} (1 - ρ) δ ε : I I : C : d ε \\ - (1 - ρ) δ ε : I d \hat{p} - \frac{1}{3} J^{- 1} δ λ I : C : d ε - J^{- 1} δ λ d \hat{p} \\ + J^{- 1} d δ λ (Δ p - Δ \hat{p}) + \bar{σ} : d δ ε + δ ε : (d W \cdot \bar{σ} - \bar{σ} \cdot d W)} d V . \end{matrix} \end{matrix}

It remains to choose $δ λ$ . To get a symmetric expression for the rate of virtual work, we choose

\begin{matrix} (11) & δ λ = (1 - ρ) J (\frac{1}{3 K} I : C : δ ε - I : δ ε + \frac{1}{K} δ \hat{p}), \end{matrix}

where

\begin{matrix} (12) & K = \frac{1}{9} I : C : I \end{matrix}

is the (instantaneous) bulk modulus. This is a suitable choice for $δ λ$ , because the (independent) term proportional to $δ \hat{p}$ ensures that the modified incremental pressure field, $Δ \hat{p}$ , is properly constrained to the incremental pressure, $Δ p$ . If we assume that the volumetric moduli $I : C$ and K change slowly with strain and ignore changes in volume, we can write for the second variation $d δ λ$ :

\begin{matrix} (13) & d δ λ = (1 - ρ) J (\frac{1}{3 K} I : C - I) : d δ ε . \end{matrix}

Hence, we find for the virtual work expression:

\begin{matrix} (14) & δ \bar{W} = \int_{V} [\tilde{σ} : δ ε + (1 - ρ) (Δ p - Δ \hat{p}) \frac{1}{K} δ \hat{p}] d V, \end{matrix}

where

\begin{matrix} (15) & \tilde{σ} = \bar{σ} + (1 - ρ) (Δ p - Δ \hat{p}) (\frac{1}{3 K} I : C - I) . \end{matrix}

For the rate of change of virtual work we find

\begin{matrix} (16) & \begin{matrix} d δ \bar{W} = \int_{V} & {δ ε : C : d ε - \frac{1}{9 K} (1 - ρ) δ ε : C^{T} : I I : C : d ε + δ ε : \bar{σ} I : d ε \\ - \frac{1}{3 K} (1 - ρ) (δ \hat{p} I : C : d ε + δ ε : C^{T} : I d \hat{p}) - \frac{1}{K} (1 - ρ) δ \hat{p} d \hat{p} \\ + \tilde{σ} : d δ ε + \bar{σ} : (δ ε \cdot d W - d W \cdot δ ε)} d V . \end{matrix} \end{matrix}

The initial stress term can be approximated by

\tilde{σ} : (d δ ε + δ ε \cdot d W - d W \cdot δ ε),

which can be written as

\tilde{σ} : [{(\frac{\partial δ u}{\partial x})}^{T} \cdot \frac{\partial d u}{\partial x} - 2 δ ε \cdot d ε],

so that the final expression for the rate of virtual work becomes

\begin{matrix} (17) & \begin{matrix} d δ \bar{W} = \int_{V} & {δ ε : C : d ε - \frac{1}{9 K} (1 - ρ) δ ε : C^{T} : I I : C : d ε + δ ε : \bar{σ} I : d ε \\ - \frac{1}{3 K} (1 - ρ) (δ \hat{p} I : C : d ε + δ ε : C^{T} : I d \hat{p}) - \frac{1}{K} (1 - ρ) δ \hat{p} d \hat{p} \\ + \tilde{σ} : [{(\frac{\partial δ u}{\partial x})}^{T} \cdot \frac{\partial d u}{\partial x} - 2 δ ε \cdot d ε]} d V . \end{matrix} \end{matrix}

The asymmetric term $δ ε : \bar{σ} I : d ε$ is significant only if large volume changes occur. Hence, the term is ignored except for material models with volumetric plasticity, such as the (capped) Drucker-Prager model and the Cam-clay model. For these models the constitutive matrix $C$ is usually asymmetric anyway so that the addition of this nonsymmetric term does not affect the cost of the analysis. It was assumed in the expression for $d δ λ$ that the (volumetric) moduli change only slowly with strain. This is not the case for material models with volumetric plasticity, in which these moduli can change abruptly. This may lead to slow convergence or even convergence failures. Failures usually occur only in higher-order elements, since in lower-order elements $Δ p - Δ \hat{p}$ approaches zero at every point and the error in $d δ λ$ has no impact.