Continuum elements with incompatible modes

The lower-order quadrilateral continuum elements in Abaqus of type CPS4I, CPE4I, CAX4I, CPEG4I, and C3D8I, as well as the related hybrid elements, are enhanced by incompatible modes to improve the bending behavior.

The following topics are discussed:

Geometrically linear formulation
Geometrically nonlinear formulation

Geometrically linear formulation

As discussed in the paper by Simo, the incompatible mode formulation can be derived in a rigorous way from the general Hu-Washizu variational principle. In this discussion we will not present this derivation but use only the key results of Simo's work.

In the incompatible mode formulation, the displacement gradient $G = \frac{\partial u}{\partial x}$ is augmented with an additional, incompatible displacement gradient field $\tilde{G}$ :

\bar{G} = G + \tilde{G} = \frac{\partial u}{\partial x} + \tilde{G} .

The incompatible displacement gradient $\tilde{G}$ is chosen internal to an element. The field cannot be selected arbitrarily. It must be independent of the regular displacement gradient.

\tilde{G} \cap G = 0,

which can also be expressed in the form

\tilde{G} - G \neq 0 for all \bar{G} \neq 0.

In addition, it must be orthogonal to any constant gradient field, which yields the condition

\begin{matrix} (1) & \int_{V_{e l}} \tilde{G} d V_{e l} = 0. \end{matrix}

If these conditions are violated, the element does not pass the patch test.

The last condition is used to obtain a suitable general form of the incompatible modes. We describe the incompatible field as a transformation of a parametric gradient field $\tilde{g} (ξ)$ :

\begin{matrix} (2) & \tilde{G} (ξ) = \frac{j (0)}{j (ξ)} \tilde{g} (ξ) \cdot t^{- 1}, \end{matrix}

where $t$ is the parametric transformation at the center of the element

t = {\frac{\partial x}{\partial ξ} |}_{ξ = 0},

$j (ξ)$ is the Jacobian of the parametric transformation at the location $ξ$ , and $j (0)$ is the Jacobian at the center of the element. For planar elements the Jacobian can be written as

j (ξ) = h (ξ) \det (\frac{\partial x}{\partial ξ}),

where h is the thickness; for axisymmetric elements it is

j (ξ) = 2 π r (ξ) \det (\frac{\partial x}{\partial ξ}),

where r is the radius; and for three-dimensional elements it is

j (ξ) = \det (\frac{\partial x}{\partial ξ}) .

Substitution of Equation 2 in Equation 1 allows us to create a simple condition for $\tilde{g} (ξ)$ :

\int_{V_{e l}} \tilde{G} d V_{e l} = \int_{V_{e l}} \frac{j (0)}{j (ξ)} \tilde{g} (ξ) \cdot t^{- 1} d V_{e l} = j (0) (\int_{V_{p a r}} \tilde{g} (ξ) d V_{p a r}) \cdot t^{- 1} = 0.

For two-dimensional elements this yields

\int_{- 1}^{+ 1} \int_{- 1}^{+ 1} \tilde{g} (ξ) d ξ_{1} d ξ_{2} = 0;

and for three-dimensional elements,

\int_{- 1}^{+ 1} \int_{- 1}^{+ 1} \int_{- 1}^{+ 1} \tilde{g} (ξ) d ξ_{1} d ξ_{2} d ξ_{3} = 0.

This makes it possible to write $\tilde{g}$ as a simple polynomial in $ξ$ . The principal contribution to $\tilde{g}$ can be written in the form

{\tilde{g}}_{p} (ξ) = α^{i} ξ_{i},

where $α^{i}$ are vectorial degrees of freedom and $ξ_{i}$ are vectors and the summation i extends over the parametric coordinates. In two-dimensional elements $ξ_{1}$ and $ξ_{2}$ are vectors of the form

ξ_{1} = {\begin{matrix} ξ_{1} \\ 0 \end{matrix}} and ξ_{2} = {\begin{matrix} 0 \\ ξ_{2} \end{matrix}},

and in three-dimensional elements

ξ_{1} = {\begin{matrix} ξ_{1} \\ 0 \\ 0 \end{matrix}}, ξ_{2} = {\begin{matrix} 0 \\ ξ_{2} \\ 0 \end{matrix}}, and ξ_{3} = {\begin{matrix} 0 \\ 0 \\ ξ_{3} \end{matrix}} .

The principal contribution to the incompatible displacement gradient hence becomes

{\tilde{G}}_{p} (ξ) = α^{i} \frac{j (0)}{j (ξ)} ξ_{i} \cdot t^{- 1} .

With the addition of these terms, parasitic shear and Poisson's effect in bending are eliminated. Note that the vectors $α^{i}$ in $\tilde{G}$ appear in a similar form as the nodal displacement vectors $u^{i}$ in $G$ ,

G = u^{i} \frac{\partial N_{i}}{\partial x},

and can be treated similar to displacement degrees of freedom.

For approximately incompressible material behavior, bilinear patterns in the hydrostatic stress can still be observed in all elements except CPS4I. These patterns can be eliminated by the introduction of additional incompatible modes of the form

{\tilde{G}}_{I} (ξ) = β^{j} \frac{j (0)}{j (ξ)} ϑ_{j} (ξ) I,

where $β^{j}$ are additional scalar degrees of freedom. For two-dimensional elements a single term $β^{1}$ is added with

ϑ_{1} = ξ_{1} ξ_{2} .

In the three-dimensional elements four additional terms $β^{j}$ are added with

ϑ_{1} = ξ_{1} ξ_{2}, ϑ_{2} = ξ_{1} ξ_{3}, ϑ_{3} = ξ_{2} ξ_{3}, and ϑ_{4} = ξ_{1} ξ_{2} ξ_{3} .

Thus, the incompatible displacement gradient $\tilde{G}$ takes the final form

\tilde{G} (ξ) = α^{i} \frac{j (0)}{j (ξ)} ξ_{i} \cdot t^{- 1} + β^{j} \frac{j (0)}{j (ξ)} ϑ_{j} (ξ) I .

The symmetric part of the incompatible displacement gradient contributes to incompatible strains:

{\begin{matrix} {\tilde{ε}}_{11} \\ {\tilde{ε}}_{22} \\ {\tilde{ε}}_{33} \\ {\tilde{γ}}_{12} \\ {\tilde{γ}}_{13} \\ {\tilde{γ}}_{23} \end{matrix}} = {\begin{matrix} {\tilde{G}}_{11} \\ {\tilde{G}}_{22} \\ {\tilde{G}}_{33} \\ {\tilde{G}}_{21} + {\tilde{G}}_{12} \\ {\tilde{G}}_{31} + {\tilde{G}}_{13} \\ {\tilde{G}}_{32} + {\tilde{G}}_{23} \end{matrix}} .

The skew-symmetric part plays no role in the geometrically linear formulation.

Geometrically nonlinear formulation

Since we want to use a formulation that can be used for any material model, we want to express the incompatible modes as a modification of the deformation gradient $F$ . The most obvious approach is to add the incompatible modes to the deformation gradient:

\begin{matrix} (3) & \bar{F} = F + \tilde{F} = \frac{\partial x}{\partial X} + \tilde{F} . \end{matrix}

This approach has been used successfully by Simo and Armero. Elements formulated on this basis satisfy the large-strain patch test; i.e., any patch of elements will be able to represent homogeneous deformations exactly. However, once the elements become distorted due to deformation, the patch test will no longer be satisfied in an incremental sense; that is, subsequent homogeneous deformations will not be represented exactly. This turns out to be a fatal flaw in the formulation for problems involving large distortions in compression.

Satisfaction of the instantaneous patch test requires the addition of an incompatible deformation rate tensor to the standard rate of deformation:

\dot{\bar{ε}} = \dot{ε} + \dot{\tilde{ε}} .

To obtain an approximate relation of this type, we write the total deformation gradient as the product of a series of incremental deformation gradients:

{\bar{F}}^{(n + 1)} = Δ {\bar{F}}^{(n + 1)} \cdot Δ {\bar{F}}^{(n)} \cdot \dots \cdot Δ {\bar{F}}^{(1)} .

Principal incompatible modes are then added to the incremental deformation gradient:

Δ {\bar{F}}_{p} = Δ F + Δ \tilde{F} = \frac{\partial x_{t + Δ t}}{\partial x_{t}} + Δ \tilde{F} .

Similar to Equation 2 the principal incompatible modes are described as a transformation of the parametric gradient field $Δ \tilde{g} (ξ)$ :

\begin{matrix} (4) & Δ \tilde{F} (ξ) = {\frac{j (0)}{j (ξ)} |}_{t} Δ \tilde{g} (ξ) \cdot t_{t}^{- 1}, \end{matrix}

where $t_{t}$ is the parametric transformation at the center of the element in the state at the start of the increment,

t_{t} = {\frac{\partial x_{t}}{\partial ξ} |}_{ξ = 0};

$j (ξ)$ is the Jacobian of the parametric transformation at the location $ξ$ ; and $j (0)$ is the Jacobian at the centroid at the start of the increment. Note that $j (ξ)$ is evaluated based on the deformation caused by the displacement degrees of freedom only and does not include the volume change due to the incompatible modes.

The incremental parametric gradient field $Δ \tilde{g} (ξ)$ has exactly the same form as in the linear formulation,

Δ \tilde{g} (ξ) = Δ α^{i} ξ_{i},

which yields the principal incremental incompatible deformation gradient

Δ \tilde{F} (ξ) = {Δ α^{i} \frac{j (0)}{j (ξ)} |}_{t} ξ_{i} \cdot t_{t}^{- 1} .

Bilinear volumetric terms are added to the principal terms in a multiplicative way:

\begin{matrix} (5) & Δ \bar{F} = Δ {\bar{F}}_{p} Δ F_{I} = (Δ F + Δ \tilde{F}) \exp ({Δ β^{j} \frac{j (0)}{j (ξ)} |}_{t} ϑ_{j} (ξ)) . \end{matrix}

The variation in the gradient of the position with respect to the current state is obtained from the fundamental relation

\begin{matrix} δ \bar{L} = δ \bar{F} \cdot {\bar{F}}^{- 1} = δ Δ \bar{F} \cdot Δ {\bar{F}}^{- 1} & = [(δ Δ F + δ Δ \tilde{F}) Δ F_{I} + Δ {\bar{F}}_{p} δ Δ F_{I}] \cdot Δ {\bar{F}}^{- 1} \\ = (δ Δ F + δ Δ \tilde{F}) \cdot Δ {\bar{F}}_{p}^{- 1} + δ Δ F_{I} Δ F_{I}^{- 1}, \end{matrix}

with

\begin{matrix} δ Δ F & = \frac{\partial δ x}{\partial x_{t}}, \\ δ Δ \tilde{F} & = \frac{j (0)}{j (ξ)} |_{t} δ \tilde{g} (ξ) \cdot t_{t}^{- 1} = δ α^{i} \frac{j (0)}{j (ξ)} |_{t} ξ_{i} \cdot t_{t}^{- 1}, \\ δ Δ F_{I} & = {δ β^{j} \frac{j (0)}{j (ξ)} |}_{t} ϑ_{j} (ξ) Δ F_{I} . \end{matrix}

This allows us to write for $δ \bar{L}$

δ \bar{L} = \frac{\partial δ x}{\partial x_{t}} \cdot Δ {\bar{F}}_{p}^{- 1} + {δ α^{i} \frac{j (0)}{j (ξ)} |}_{t} ξ_{i} \cdot t_{t}^{- 1} \cdot Δ {\bar{F}}_{p}^{- 1} + {δ β^{j} \frac{j (0)}{j (ξ)} |}_{t} ϑ_{j} (ξ) .

The integral of the principal incompatible modes over the element volume at the start of the increment is, hence, equal to

\begin{matrix} \int_{v_{e l}} δ Δ \tilde{F} \cdot Δ {\bar{F}}_{p}^{- 1} d v_{e l} & = {\int_{v_{e l}} \frac{j (0)}{j (ξ)} |}_{t} δ \tilde{g} (ξ) \cdot t_{t}^{- 1} \cdot Δ {\bar{F}}_{p}^{- 1} d v_{e l} \\ = {\int_{V_{p a r}} j (0) |}_{t} δ \tilde{g} (ξ) \cdot t_{t}^{- 1} \cdot Δ {\bar{F}}_{p}^{- 1} d V_{p a r} . \end{matrix}

Note that the integral will vanish if the incremental deformation is homogeneous; that is, $Δ {\bar{F}}_{p} (ξ) = Δ {\bar{F}}_{0}$ , since in that case the integral can be written as

\int_{v_{e l}} δ Δ \tilde{F} \cdot Δ {\bar{F}}_{p}^{- 1} d v_{e l} = {j (0) |}_{t} (\int_{V_{p a r}} δ {\tilde{g}}_{p} (ξ) d V_{p a r}) \cdot t_{t}^{- 1} \cdot Δ {\bar{F}}_{0}^{- 1} = 0.

Hence, the incremental patch test will be satisfied. The rate of the gradient of the position with respect to the current state is obtained similar to the variation

\begin{matrix} d \bar{L} = d \bar{F} \cdot {\bar{F}}^{- 1} = d Δ \bar{F} \cdot Δ {\bar{F}}^{- 1} & = [(d Δ F + d Δ \tilde{F}) Δ F_{I} + Δ {\bar{F}}_{p} d Δ F_{I}] \cdot Δ {\bar{F}}^{- 1} \\ = (d Δ F + d Δ \tilde{F}) \cdot Δ {\bar{F}}_{p}^{- 1} + d Δ F_{I} Δ F_{I}^{- 1}, \end{matrix}

with

\begin{matrix} d Δ F & = \frac{\partial d x}{\partial x_{t}}, \\ d Δ \tilde{F} & = \frac{j (0)}{j (ξ)} |_{t} d \tilde{g} (ξ) \cdot t_{t}^{- 1} = d α^{i} \frac{j (0)}{j (ξ)} |_{t} ξ_{i} \cdot t_{t}^{- 1}, \\ d Δ F_{I} & = {d β^{j} \frac{j (0)}{j (ξ)} |}_{t} ϑ_{j} (ξ) Δ F_{I} . \end{matrix}

For a finite-strain increment we use the midincrement approach proposed by Hughes and Winget. This yields

Δ \bar{L} \approx Δ \bar{F} \cdot Δ {\bar{F}}_{t + Δ t / 2}^{- 1} .

The second variation is obtained in the usual way. Since the regular deformation gradient and the principal incompatible modes are purely displacement based, the initial stress stiffness terms are readily obtained as

\int_{v_{e l}} σ : (δ {\bar{L}}_{p}^{T} \cdot d {\bar{L}}_{p} - 2 δ {\bar{D}}_{p} \cdot d {\bar{D}}_{p}) d v_{e l},

where $d {\bar{L}}_{p} = d {\bar{F}}_{p} \cdot Δ {\bar{F}}_{p}^{- 1}$ is the gradient of the velocity and $δ {\bar{L}}_{p} = δ {\bar{F}}_{p} \cdot Δ {\bar{F}}_{p}^{- 1}$ is the gradient of the displacement variation including the primary incompatible modes. Further, $d {\bar{D}}_{p} = \frac{1}{2} (d {\bar{L}}_{p} + d {\bar{L}}_{p}^{T})$ is the rate of deformation, and $δ {\bar{D}}_{p} = \frac{1}{2} (δ {\bar{L}}_{p} + δ {\bar{L}}_{p}^{T})$ is the variation in the deformation.

The additional bilinear modes appear in the first variation only as variations (the values $β^{i}$ themselves do not appear); hence, the contributions to the second variation can be neglected.