Finite-strain shell element formulation

Geometric description

At a given stage in the deformation history of the shell, the position of a material point in the shell is defined by

x (S_{i}) = \bar{x} (S_{α}) + {\bar{f}}_{33} (S_{α}) t_{3} (S_{α}) S_{3},

where the subscript i and other Roman subscripts range from 1 to 3. Subscripts $α$ and other lowercase Greek subscripts which describe the quantities in the reference surface of the shell range from 1 to 2. In the above equation $t_{3}$ is the normal to the reference surface of the shell. The gradient of the position is

\frac{\partial x}{\partial S_{β}} = \frac{\partial \bar{x}}{\partial S_{β}} + {\bar{f}}_{33} \frac{\partial t_{3}}{\partial S_{β}} S_{3}, \frac{\partial x}{\partial S_{3}} = {\bar{f}}_{33} t_{3},

where we have neglected derivatives of ${\bar{f}}_{33}$ with respect to $S_{β}$ . Note that in the above $S_{α}$ are local surface coordinates that are assumed to be orthogonal and distance measuring in the reference state. $S_{3}$ is the coordinate in the thickness direction, distance measuring and orthogonal to $S_{α}$ in the reference state. The thickness increase factor ${\bar{f}}_{33}$ is assumed to be independent of $S_{3}$ .

In the deformed state we define local, orthonormal shell directions $t_{i}$ such that

t_{i} \cdot t_{j} = δ_{i j}, t_{i} t_{i} = I,

where $δ_{i j}$ is the Kronecker delta and $I$ is the identity tensor of rank 2. Summation convention is used for repeated subscripts. The in-plane components of the gradient of the position are obtained as

f_{α β} = t_{α} \cdot \frac{\partial x}{\partial S_{β}} = {\bar{f}}_{α β} + B_{α β} {\bar{f}}_{33} S_{3},

where we have introduced the reference surface deformation gradient

{{\bar{f}}_{α β} \overset{def}{=} t_{α} \cdot \frac{\partial x}{\partial S_{β}} |}_{S_{3} = 0} = t_{α} \cdot \frac{\partial \bar{x}}{\partial S_{β}}

and the reference surface normal gradient

B_{α β} \overset{def}{=} t_{α} \cdot \frac{\partial t_{3}}{\partial S_{β}} .

In the original (reference) configuration we denote the position by $X$ ( $\bar{X}$ for the reference surface) and the direction vectors by $T_{i}$ , which yields

X (S_{i}) = \bar{X} (S_{α}) + T_{3} (S_{α}) S_{3} .

The gradient of the position is

\frac{\partial X}{\partial S_{β}} = \frac{\partial \bar{X}}{\partial S_{β}} + \frac{\partial T_{3}}{\partial S_{β}} S_{3}, \frac{\partial X}{\partial S_{3}} = T_{3},

and the in-plane components of the gradient are obtained as

f_{α β}^{\circ} = T_{α} \cdot \frac{\partial X}{\partial S_{β}} = δ_{α β} + B_{α β}^{\circ} S_{3},

where we have assumed that the in-plane direction vectors follow from the surface coordinates with

T_{β} = {\frac{\partial X}{\partial S_{β}} |}_{S_{3} = 0} = \frac{\partial \bar{X}}{\partial S_{β}}

and defined the original reference surface normal gradient,

B_{α β}^{\circ} \overset{def}{=} T_{α} \cdot \frac{\partial T_{3}}{\partial S_{β}} .

The original reference surface normal gradient is obtained in the finite element formulation from the interpolation of the nodal normals with the shape functions. In the deformed configuration it is not derived from the nodal normals but is updated independently based on the gradient of the incremental rotations.

Parametric interpolation

The position of the points in the shell reference surface is described in terms of discrete nodal positions with parametric interpolation functions $N^{I} (ξ_{α})$ . The functions are $C_{\circ}$ continuous, and $ξ_{α}$ are nonorthogonal, nondistance measuring parametric coordinates. For the reference surface positions one, thus, obtains

\bar{x} (ξ_{α}) = N^{I} (ξ_{α}) {\bar{x}}^{I}, \bar{X} (ξ_{α}) = N^{I} (ξ_{α}) {\bar{X}}^{I} .

The gradients of the position with respect to $ξ_{β}$ are

\frac{\partial \bar{x}}{\partial ξ_{β}} = \frac{\partial N^{I}}{\partial ξ_{β}} {\bar{x}}^{I}, \frac{\partial \bar{X}}{\partial ξ_{β}} = \frac{\partial N^{I}}{\partial ξ_{β}} {\bar{X}}^{I} .

Note that uppercase Roman superscripts such as I denote nodes of an element and that repeated superscripts imply summation over all nodes of an element.

Now consider the original configuration. The unit normal to the shell reference surface is readily obtained as

T_{3} = (\frac{\partial \bar{X}}{\partial ξ_{1}} \times \frac{\partial \bar{X}}{\partial ξ_{2}}) / | \frac{\partial \bar{X}}{\partial ξ_{1}} \times \frac{\partial \bar{X}}{\partial ξ_{2}} | .

Subsequently, we define two orthonormal tangent vectors $T_{α}$ and distance measuring coordinates $S_{α}$ along these vectors. The derivatives of these coordinates with respect to $ξ_{β}$ follow from

\frac{\partial S_{α}}{\partial ξ_{β}} = T_{α} \cdot \frac{\partial \bar{X}}{\partial ξ_{β}} = T_{α} \cdot {\bar{X}}^{I} \frac{\partial N^{I}}{\partial ξ_{β}} .

The gradient of $ξ_{α}$ with respect to $S_{β}$ is readily obtained by inversion:

\frac{\partial ξ_{α}}{\partial S_{β}} = {[\frac{\partial S_{α}}{\partial ξ_{β}}]}^{- 1},

which makes it possible to obtain the gradient operator

\frac{\partial N^{I}}{\partial S_{β}} = \frac{\partial N^{I}}{\partial ξ_{α}} \frac{\partial ξ_{α}}{\partial S_{β}} .

The original reference surface normal gradient is obtained from the nodal normals $T_{3}^{I}$ with

B_{α β}^{\circ} = T_{α} \cdot T_{3}^{I} \frac{\partial N^{I}}{\partial S_{β}} .

Since the original reference surface normal gradient is obtained by taking derivatives with respect to orthogonal distance measuring coordinates, we will call $B_{α β}^{\circ} = b_{α β}^{\circ}$ the original curvature of the reference surface.

Membrane deformation and curvature

It is convenient to define the inverse of the reference surface deformation gradient

{\bar{h}}_{α β} = {[{\bar{f}}_{α β}]}^{- 1} .

With this expression we can define the gradient operator in the current state:

\frac{\partial}{\partial s_{β}} \overset{def}{=} {\bar{h}}_{α β} \frac{\partial}{\partial S_{α}}, or inverted \frac{\partial}{\partial S_{β}} = {\bar{f}}_{α β} \frac{\partial}{\partial s_{α}} .

The gradient operator in the current state can also be defined as the derivative with respect to distance measuring coordinates $s_{α}$ along the base vectors $t_{α}$ , since

t_{α} \cdot \frac{\partial \bar{x}}{\partial s_{β}} = t_{α} \cdot \frac{\partial \bar{x}}{\partial S_{γ}} {\bar{h}}_{γ β} = {\bar{f}}_{α γ} {\bar{h}}_{γ β} = δ_{α β}

and, hence,

t_{α} = \frac{\partial \bar{x}}{\partial s_{α}} .

Hence, it is possible to write for the ${\bar{h}}_{α β}$ :

{\bar{h}}_{α β} = T_{α} \cdot \frac{\partial \bar{X}}{\partial s_{β}}

since

{\bar{f}}_{α γ} {\bar{h}}_{γ β} = t_{α} \cdot \frac{\partial \bar{x}}{\partial \bar{X}} \cdot T_{γ} T_{γ} \frac{\partial \bar{X}}{\partial \bar{x}} \cdot t_{β} = t_{α} \cdot \frac{\partial \bar{x}}{\partial \bar{X}} \cdot \frac{\partial \bar{X}}{\partial \bar{x}} \cdot t_{β} = t_{α} \cdot t_{β} = δ_{α β} .

In an incremental analysis we can also define the incremental deformation tensor

Δ {\bar{f}}_{α β} = t_{α}^{t + Δ t} \cdot \frac{\partial {\bar{x}}^{t + Δ t}}{\partial s_{β}^{t}}

and its inverse

Δ {\bar{h}}_{α β} = t_{α}^{t} \cdot \frac{\partial {\bar{x}}^{t}}{\partial s_{β}^{t + Δ t}} .

With a local coordinate system defined in the current state, the current gradient of the normal can be transformed into the curvature of the surface:

b_{α β} \overset{def}{=} t_{α} \cdot \frac{\partial t_{3}}{\partial s_{β}} = B_{α γ} {\bar{h}}_{γ β} .

Orientation update

The equations given in the earlier sections are valid for any local coordinate system defined in the current state. The $t_{α}$ vectors at the beginning of the analysis are determined following the standard Abaqus conventions. In this section, we outline the way in which the in-plane coordinates are made corotational.

To obtain the updated version of $t_{α}$ , we follow a two-step approach. First, we construct orthogonal vectors ${\hat{t}}_{α}$ tangential to the surface (following Abaqus conventions). Subsequently, we calculate

{\hat{f}}_{α β} = {\hat{t}}_{α} \cdot \frac{\partial {\bar{x}}^{t + Δ t}}{\partial S_{β}} .

We then apply an in-plane rotation $Δ R_{α β}$ to the vectors: ${\hat{t}}_{α}$ :

Δ R_{11} = Δ R_{22} = \cos Δ ψ,

Δ R_{21} = - Δ R_{12} = \sin Δ ψ,

where $Δ ψ$ is to be determined such that the resulting deformation tensor is symmetric, as

{\bar{f}}_{α β} = Δ R_{α γ} {\hat{f}}_{γ β} = {\hat{f}}_{γ α} Δ R_{β γ} = {\bar{f}}_{β α} .

From this follows

\tan Δ ψ = \frac{({\hat{f}}_{12} - {\hat{f}}_{21})}{{\hat{f}}_{11} + {\hat{f}}_{22}} .

Thus, we can calculate the updated local material directions as

{\bar{t}}_{α}^{t + Δ t} = Δ R_{α γ} {\hat{t}}_{γ} .

Curvature change

We assume that the nodal spin will be interpolated with the interpolation functions $N^{I} (ξ_{α})$ . During an increment the nodal spin is assumed to be constant; consequently, the value of the spin at each material point will be constant. Hence, we can use the same interpolation functions for the incremental finite rotation vector $Δ ϕ$ :

Δ ϕ = ω Δ t = N^{I} (ξ_{α}) ω^{I} Δ t = N^{I} (ξ_{α}) Δ ϕ^{I} .

The finite rotation vector can be split in a rotation amplitude $Δ ϕ$ and a rotation axis $p$ :

Δ ϕ = Δ ϕ p, with Δ ϕ \overset{def}{=} | Δ ϕ | and p = Δ ϕ / Δ ϕ .

To rotate the shell normal, we use quaternion algebra. The incremental nodal rotation is represented by the rotation quaternion $Δ q$ , which is defined by

Δ q \overset{def}{=} (\cos \frac{Δ ϕ}{2}, \sin \frac{Δ ϕ}{2} p) .

An updated shell normal is then obtained according to

{\tilde{t}}_{3}^{t + Δ t} = Δ q t_{3}^{t} Δ q^{†} .

This updated shell normal does not actually have to be calculated: it is used only for the derivation of the expression for the curvature change. It is not equal to the shell normal used at the start of the next increment $t_{3}^{t + Δ t}$ , which will again be chosen perpendicular to the reference surface. The updated normal used here will be approximately orthogonal to the reference surface, depending upon the amount of transverse shear deformation. The gradient of the updated shell normal can be obtained by differentiation:

\frac{\partial {\tilde{t}}_{3}^{t + Δ t}}{\partial S_{β}} = \frac{\partial Δ q}{\partial S_{β}} t_{3}^{t} Δ q^{†} + Δ q t_{3}^{t} \frac{\partial Δ q^{†}}{\partial S_{β}} + Δ q \frac{\partial t_{3}^{t}}{\partial S_{β}} Δ q^{†} .

The second term on the right-hand side can be written in the form

Δ q t_{3}^{t} \frac{\partial Δ q^{†}}{\partial S_{β}} = {[\frac{\partial Δ q}{\partial S_{β}} {(t_{3}^{t})}^{†} Δ q^{†}]}^{†} = - {[\frac{\partial Δ q}{\partial S_{β}} t_{3}^{t} Δ q^{†}]}^{†} .

Hence, the scalar parts of the first two terms cancel each other and the vector parts reinforce each other, leading to

\frac{\partial {\tilde{t}}_{3}^{t + Δ t}}{\partial S_{β}} = 2 V (\frac{\partial Δ q}{\partial S_{β}} t_{3}^{t} Δ q^{†}) + Δ q \frac{\partial t_{3}^{t}}{\partial S_{β}} Δ q^{†} .

The inverse of a rotation quaternion such as $Δ q$ is equal to its conjugate ( $Δ q^{- 1} = Δ q^{†}$ ). Hence, we can write

\begin{matrix} \frac{\partial {\tilde{t}}_{3}^{t + Δ t}}{\partial S_{β}} & ​ = 2 V (\frac{\partial Δ q}{\partial S_{β}} Δ q^{†} Δ q t_{3}^{t} Δ q^{†}) + Δ q \frac{\partial t_{3}^{t}}{\partial S_{β}} Δ q^{†} \\ ​ = V (2 \frac{\partial Δ q}{\partial S_{β}} Δ q^{†} {\tilde{t}}_{3}^{t + Δ t}) + Δ q \frac{\partial t_{3}^{t}}{\partial S_{β}} Δ q^{†} \\ ​ = Δ R_{β} \times {\tilde{t}}_{3}^{t + Δ t} + Δ q \frac{\partial t_{3}^{t}}{\partial S_{β}} Δ q^{†}, \end{matrix}

where we have formally defined the incremental gradient update vectors

Δ R_{β} \overset{def}{=} V (2 \frac{\partial Δ q}{\partial S_{β}} Δ q^{†}),

which must be expressed in terms of the gradient of the incremental rotation. From the definition of the incremental quaternion $Δ q$ follows

\frac{\partial Δ q}{\partial S_{β}} = (\frac{1}{2} \sin \frac{Δ ϕ}{2} \frac{\partial Δ ϕ}{\partial S_{β}}, \frac{1}{2} \cos \frac{Δ ϕ}{2} p \frac{\partial Δ ϕ}{\partial S_{β}} + \sin \frac{Δ ϕ}{2} \frac{\partial p}{\partial S_{β}});

thus, for $Δ R_{β}$ , again with use of the incremental quaternion definition

Δ R_{β} = p \frac{\partial Δ ϕ}{\partial S_{β}} + \sin Δ ϕ \frac{\partial p}{\partial S_{β}} + (1 - \cos Δ ϕ) p \times \frac{\partial p}{\partial S_{β}} .

From the definition of $Δ ϕ$ and $p$ follows

\begin{matrix} \frac{\partial Δ ϕ}{\partial S_{β}} & ​ = \frac{1}{Δ ϕ} Δ ϕ \cdot \frac{\partial Δ ϕ}{\partial S_{β}}, \\ \frac{\partial p}{\partial S_{β}} & = \frac{1}{Δ ϕ} \frac{\partial Δ ϕ}{\partial S_{β}} - \frac{1}{Δ ϕ^{3}} Δ ϕ Δ ϕ \cdot \frac{\partial Δ ϕ}{\partial S_{β}} = \frac{1}{Δ ϕ} [\frac{\partial Δ ϕ}{\partial S_{β}} - p p \cdot \frac{\partial Δ ϕ}{\partial S_{β}}] . \end{matrix}

After substitution in the expression for $Δ R_{β}$ and some algebra one obtains

Δ R_{β} = \frac{\sin Δ ϕ}{Δ ϕ} \frac{\partial Δ ϕ}{\partial S_{β}} + \frac{1 - \cos Δ ϕ}{Δ ϕ} p \times \frac{\partial Δ ϕ}{\partial S_{β}} + (1 - \frac{\sin Δ ϕ}{Δ ϕ}) p p \cdot \frac{\partial Δ ϕ}{\partial S_{β}} .

Note that $Δ R_{β} \to \frac{\partial Δ ϕ}{\partial S_{β}}$ when $Δ ϕ \to 0$ .

For the gradient $B_{α β}$ of the updated shell normal we obtain

\begin{matrix} B_{α β}^{t + Δ t} & = t_{α}^{t + Δ t} \cdot \frac{\partial {\tilde{t}}_{3}^{t + Δ t}}{\partial S_{β}} = t_{α}^{t + Δ t} [Δ R_{β} \times {\tilde{t}}_{3}^{t + Δ t} + Δ q t_{γ}^{t} B_{α β}^{t} Δ q^{†}] \\ = Δ R_{β} \cdot ({\tilde{t}}_{3}^{t + Δ t} \times t_{α}^{t + Δ t}) + t_{α}^{t + Δ t} \cdot t_{γ}^{t + Δ t} B_{γ β}^{t} = ϵ_{α}^{γ} t_{γ}^{t + Δ t} \cdot Δ R_{β} + B_{α β}^{t}, \end{matrix}

where we have introduced the two-dimensional alternator $ϵ_{α}^{β}$ :

ϵ_{1}^{1} = ϵ_{2}^{2} = 0, ϵ_{1}^{2} = - ϵ_{2}^{1} = 1.

Note that the change in $B_{α β}$ is independent of $B_{α β}^{t}$ .

Calculation of $B_{α β}$ involves taking the gradient with respect to the reference configuration. It is more convenient to use the reference surface curvature tensor

b_{α β} \overset{def}{=} B_{α γ} {\bar{h}}_{γ β} = t_{α} \cdot \frac{\partial t_{3}}{\partial s_{β}} .

We then introduce the incremental curvature update vectors

Δ r_{β} \overset{def}{=} Δ R_{γ} {\bar{h}}_{γ β} = \frac{\sin Δ ϕ}{Δ ϕ} \frac{\partial Δ ϕ}{\partial s_{β}} + \frac{1 - \cos Δ ϕ}{Δ ϕ} p \times \frac{\partial Δ ϕ}{\partial s_{β}} + (1 - \frac{\sin Δ ϕ}{Δ ϕ}) p p \cdot \frac{\partial Δ ϕ}{\partial s_{β}},

which makes it possible to write the update equation as

b_{α β}^{t + Δ t} = ϵ_{α}^{γ} t_{γ}^{t + Δ t} \cdot Δ r_{β} + b_{α γ}^{t} Δ {\bar{h}}_{γ β} .

This expression makes it feasible to calculate the update in the reference surface curvature by taking gradients in the latest updated state only.

Deformation gradient

We already have obtained an expression for the deformation gradient in the reference surface, and we have assumed that the thickness change is constant:

{\bar{F}}_{α β} = {\bar{f}}_{α β}, {\bar{F}}_{33} = {\bar{f}}_{33} .

At other points in the shell we obtain for the in-plane component

F_{α β} = t_{α} \cdot \frac{\partial x}{\partial S_{γ}} {(T_{β} \cdot \frac{\partial X}{\partial S_{γ}})}^{- 1} = ({\bar{f}}_{α γ} + {\bar{f}}_{33} S_{3} B_{α γ}) {(δ_{γ β} + S_{3} B_{γ β}^{\circ})}^{- 1} .

We neglect terms of order ${(S_{3})}^{2}$ , which yields the simplified relation

F_{α β} = {\bar{F}}_{α β} + S_{3} ({\bar{f}}_{33} B_{α β} - {\bar{f}}_{α γ} B_{γ β}^{\circ}) .

We can write this as the product of a finite-membrane deformation and a bending perturbation:

\begin{matrix} F_{α β} & ​ = [δ_{α γ} + S_{3} ({\bar{f}}_{33} B_{α δ} {\bar{h}}_{δ γ} - {\bar{f}}_{α δ} B_{δ ε}^{\circ} {\bar{h}}_{ε γ})] {\bar{f}}_{γ β} \\ ​ = [δ_{α γ} + S_{3} ({\bar{f}}_{33} b_{α γ} - {\bar{f}}_{α δ} b_{δ ε}^{\circ} {\bar{h}}_{ε γ})] {\bar{f}}_{γ β} . \end{matrix}

It will be assumed that the deformation (strain and rotation) due to bending is small and, therefore,

S_{3} ({\bar{f}}_{33} b_{α γ} - {\bar{f}}_{α δ} b_{δ ε}^{\circ} {\bar{h}}_{ε γ}) ≪ 1.

Membrane strain increment

The membrane strain increment follows from the incremental stretch tensor $Δ V$ , whose components follow from the incremental deformation gradient $Δ {\bar{f}}_{α β}$ by the polar decomposition $Δ {\bar{f}}_{α β} = Δ {\bar{V}}_{α γ} Δ {\bar{R}}_{γ β}$ .

Let ${\bar{f}}_{α β}^{t}$ and ${\bar{f}}_{α β}^{t + Δ t}$ be the deformation gradient at the beginning and the end of the increment, respectively. By definition ${\bar{f}}_{α β}^{t + Δ t} = Δ {\bar{f}}_{α δ} {\bar{f}}_{δ β}^{t}$ . The incremental deformation gradient follows as

Δ {\bar{f}}_{α β} = {\bar{f}}_{α δ}^{t + Δ t} {({\bar{f}}^{t}^{- 1})}_{δ β} .

Since $Δ {\bar{R}}_{γ β}$ are the components of an orthogonal matrix, the square of the incremental stretch tensor can be obtained by

Δ {\bar{f}}_{α γ} Δ {\bar{f}}_{β γ} = Δ {\bar{V}}_{α γ} Δ {\bar{V}}_{β γ} = \overset{2}{\sum_{I = 1}} {(Δ λ_{I})}^{2} a_{α}^{I} a_{β}^{I}

(see Deformation). The logarithmic strain increment is then

Δ ϵ_{α β} = \overset{2}{\sum_{I = 1}} \ln (Δ λ_{I}) a_{α}^{I} a_{β}^{I}

and the average material rotation increment is defined from the polar decomposition:

Δ {\bar{R}}_{α β} = \overset{2}{\sum_{I = 1}} \frac{1}{Δ λ_{I}} a_{α}^{I} a_{γ}^{I} Δ {\bar{f}}_{γ β} .

Due to the choice of the element basis directions, it follows that

Δ {\bar{R}}_{α β} \approx δ_{α β} .

Curvature increment

Following Koiter-Sanders shell theory, and compensating for the rotation of the base vectors relative to the material, we define the physical curvature increment $Δ κ_{α β}$ as

Δ κ_{α β} = sym [b_{α β}^{t + Δ t} - b_{α γ}^{t} Δ {\bar{R}}_{β γ} + b_{α γ}^{t} Δ {\bar{R}}_{δ γ} Δ {\bar{ε}}_{δ β}] = sym [b_{α β}^{t + Δ t} - b_{α γ}^{t} Δ {\bar{R}}_{δ γ} (δ_{δ β} - Δ {\bar{ε}}_{δ β})] .

Neglecting terms of the order ${(Δ {\bar{ε}}_{α β})}^{2}$ relative to $Δ {\bar{ε}}_{α β}$ , this expression can be rewritten as

Δ κ_{α β} = sym [b_{α β}^{t + Δ t} - b_{α γ}^{t} Δ {\bar{R}}_{δ γ} Δ {\bar{V}}_{δ β}^{- 1}] = sym [b_{α β}^{t + Δ t} - b_{α γ}^{t} Δ {\bar{h}}_{γ β}] = sym [ϵ_{α}^{γ} t_{γ}^{t + Δ t} \cdot Δ r_{β}],

where use was made of the curvature update formula. Observe that the curvature at the beginning of the increment, $b_{α β}^{t}$ , does not appear in this equation. Hence, there is no need to calculate the initial curvature $b_{α β}^{\circ}$ , and we can assume $b_{α β}^{\circ} = 0$ . The deformation gradient can, hence, also be simplified to

F_{α β} = {\bar{f}}_{α β} + S_{3} {\bar{f}}_{33} b_{α γ} {\bar{f}}_{γ β} .

For the material strain increment at a point through the shell thickness Koiter-Sanders theory thus yields

Δ ε_{α β} = Δ {\bar{ε}}_{α β} + {\bar{f}}_{33}^{t + Δ t} S_{3} Δ κ_{α β} .

Virtual work

The virtual work contribution of the stresses is

δ Π = \int_{V} σ_{α β} δ ε_{α β} d V .

We assume that the variations in the strain can be expressed in terms of variation in membrane strain and curvature with the same relations as apply to the increment in strain:

δ ε_{α β} = δ {\bar{ε}}_{α β} + {\bar{f}}_{33}^{t + Δ t} S_{3} δ κ_{α β},

which transforms the virtual work equation into

δ Π = \int_{V} σ_{α β} (δ {\bar{ε}}_{α β} + {\bar{f}}_{33} S_{3} δ κ_{α β}) d V .

We introduce the membrane forces $N_{α β}$ and the bending moments $M_{α β}$ :

N_{α β} \overset{def}{=} \int_{h} σ_{α β} {\bar{f}}_{33} d S_{3}, M_{α β} \overset{def}{=} \int_{h} σ_{α β} {\bar{f}}_{33}^{2} S_{3} d S_{3},

which allows us to write

δ Π = \int_{A} (N_{α β} δ {\bar{ε}}_{α β} + M_{α β} δ κ_{α β}) d A .

The membrane strain variation follows with the usual expression

δ {\bar{ε}}_{α β} = sym (δ {\bar{f}}_{α γ} {\bar{h}}_{γ β}) = sym (δ t_{α} \cdot t_{β} + t_{α} \cdot \frac{\partial δ \bar{x}}{\partial s_{β}}) = sym (t_{α} \cdot \frac{\partial δ \bar{x}}{\partial s_{β}}),

where we have used the identity $sym (δ t_{α} \cdot t_{β}) = δ (t_{α} \cdot t_{β}) = 0$ .

The variation in the curvature is obtained by taking variations in the incremental curvature, which yields

\begin{matrix} δ κ_{α β} & ​ = sym (ϵ_{α}^{γ} t_{γ} \cdot δ Δ r_{β} + ϵ_{α}^{γ} δ t_{γ} \cdot Δ r_{β}) \\ ​ = sym (ϵ_{α}^{γ} t_{γ} \cdot δ Δ R_{δ} {\bar{h}}_{δ β} + ϵ_{α}^{γ} t_{γ} \cdot Δ R_{δ} δ {\bar{h}}_{δ β} + ϵ_{α}^{γ} δ t_{γ} \cdot Δ R_{δ} {\bar{h}}_{δ β}) . \end{matrix}

We neglect the terms of order $Δ κ_{α β}$ and also terms of order $t_{3} \cdot Δ r_{δ}$ , which yields

δ κ_{α β} = sym (ϵ_{α}^{γ} t_{γ} \cdot δ Δ r_{δ} {\bar{h}}_{δ β}) .

We evaluate $δ Δ r_{δ}$ with respect to the current state (at the end of the increment). Hence for the evaluation we can assume $Δ ϕ = 0$ . Moreover, we neglect terms of the order $t_{α} \cdot \frac{\partial δ Δ ϕ}{\partial S_{β}}$ since they are proportional to $Δ κ_{α β}$ . Hence, we obtain

t_{γ} δ Δ r_{δ} \approx t_{γ} \cdot \frac{\partial δ ϕ}{\partial S_{δ}},

which substituted in the expression for $δ κ_{α β}$ yields

δ κ_{α β} = sym (ϵ_{α}^{γ} t_{γ} \cdot \frac{\partial δ ϕ}{\partial S_{δ}} {\bar{h}}_{δ β}) = sym (ϵ_{α}^{γ} t_{γ} \cdot \frac{\partial δ ϕ}{\partial s_{β}}) .

The rate of virtual work

To obtain an expression for the rate of virtual work, we first write the virtual work equation in terms of the reference volume

δ Π = \int_{V^{\circ}} τ_{α β} δ ε_{α β} d V^{\circ} = \int_{A^{\circ}} \int_{h} τ_{α β} δ ε_{α β} d S_{3} d A^{\circ},

where $τ_{α β}$ is the Kirchhoff stress tensor, related to the Cauchy or true stress tensor via

τ_{α β} = J σ_{α β} .

The rate of change then becomes

d δ Π = \int_{A^{\circ}} \int_{h} (d^{\nabla} τ_{α β} δ ε_{α β} + τ_{α β} d^{\nabla} δ ε_{α β}) d S_{3} d A^{\circ} .

Here $d^{\nabla}$ indicates that the rates are taken in a material, corotational coordinate system. The terms involving stress rates are related to the material behavior. We assume constitutive equations of the form

d^{\nabla} τ_{α β} = J C_{α β γ δ} d ε_{δ γ} .

Substituted in the expression $d δ Π$ and transformed back to the current configuration, this yields

d δ Π = \int_{A} \int_{h} (δ ε_{α β} C_{α β γ δ} d ε_{γ δ} + σ_{α β} d^{\nabla} δ ε_{α β}) {\bar{f}}_{33} d S_{3} d A .

Consistent with the derivation of the virtual work equation itself, we neglect terms of the order $d {\bar{f}}_{33} S_{3} δ κ_{α β}$ . Hence, the rate of virtual work can be written as

\begin{matrix} d δ Π = ​ & \int_{A} [\int_{h} (δ {\bar{ε}}_{α β} + {\bar{f}}_{33} S_{3} δ κ_{α β}) C_{α β γ δ} (d {\bar{ε}}_{γ δ} + {\bar{f}}_{33} S_{3} d κ_{γ δ}) {\bar{f}}_{33} d S_{3} + \\ N_{α β} d^{\nabla} δ {\bar{ε}}_{α β} + M_{α β} d^{\nabla} δ κ_{α β}] d A . \end{matrix}

Second variation of the membrane strain

It remains to determine $d^{\nabla} δ {\bar{ε}}_{α β}$ and $d^{\nabla} δ κ_{α β}$ . From the first variation $δ {\bar{ε}}_{α β}$ follows

d^{\nabla} δ {\bar{ε}}_{α β} = sym (d^{\nabla} δ {\bar{f}}_{α γ} {\bar{h}}_{γ β} + δ {\bar{f}}_{α γ} d^{\nabla} {\bar{h}}_{γ β}) .

Since ${\bar{h}}_{α β}$ is the inverse of ${\bar{f}}_{α β}$ , it follows that

d^{\nabla} {\bar{h}}_{α β} = - {\bar{h}}_{α γ} d^{\nabla} {\bar{f}}_{γ δ} {\bar{h}}_{δ β} .

Substitution in the expression for the second variation yields

\begin{matrix} d^{\nabla} δ {\bar{ε}}_{α β} & ​ = sym (d^{\nabla} t_{α} \cdot \frac{\partial δ \bar{x}}{\partial s_{β}} - t_{α} \cdot \frac{\partial δ \bar{x}}{\partial s_{γ}} d^{\nabla} t_{γ} \cdot \frac{\partial \bar{x}}{\partial s_{β}} - t_{α} \cdot \frac{\partial δ \bar{x}}{\partial s_{γ}} t_{γ} \cdot \frac{\partial d^{\nabla} \bar{x}}{\partial s_{β}}) \\ ​ = sym (d^{\nabla} t_{α} \cdot \frac{\partial δ \bar{x}}{\partial s_{β}} - d^{\nabla} t_{γ} \cdot t_{β} t_{α} \cdot \frac{\partial δ \bar{x}}{\partial s_{γ}} - \frac{\partial d^{\nabla} \bar{x}}{\partial s_{β}} \cdot t_{γ} t_{α} \cdot \frac{\partial δ \bar{x}}{\partial s_{γ}}) . \end{matrix}

The corotational rate of the base vectors follows from

d^{\nabla} t_{β} = \frac{\partial d \bar{x}}{\partial s_{β}} + \frac{\partial \bar{x}}{\partial S_{γ}} d^{\nabla} {\bar{h}}_{γ β} = \frac{\partial d \bar{x}}{\partial s_{β}} - d^{\nabla} t_{γ} \cdot t_{β} t_{γ} - \frac{\partial d \bar{x}}{\partial s_{β}} \cdot t_{γ} t_{γ} .

Substituted in the first term of the previous expression yields

\begin{matrix} d^{\nabla} δ {\bar{ε}}_{α β} = ​ & sym [\frac{\partial d \bar{x}}{\partial s_{α}} \cdot \frac{\partial δ \bar{x}}{\partial s_{β}} - d^{\nabla} t_{γ} \cdot t_{α} t_{γ} \cdot \frac{\partial δ \bar{x}}{\partial s_{β}} - d^{\nabla} t_{γ} \cdot t_{β} t_{α} \cdot \frac{\partial δ \bar{x}}{\partial s_{γ}} - \\ \frac{\partial d \bar{x}}{\partial s_{α}} \cdot t_{γ} t_{γ} \cdot \frac{\partial δ \bar{x}}{\partial s_{β}} - \frac{\partial d \bar{x}}{\partial s_{β}} \cdot t_{γ} t_{α} \cdot \frac{\partial δ \bar{x}}{\partial s_{γ}}] \\ = & sym [\frac{\partial d \bar{x}}{\partial s_{α}} \cdot \frac{\partial δ \bar{x}}{\partial s_{β}} - 2 (d^{\nabla} t_{γ} \cdot t_{α} + \frac{\partial d \bar{x}}{\partial s_{α}} \cdot t_{γ}) δ {\bar{ε}}_{γ β}] . \end{matrix}

The in-plane components of the corotational rate of the base vectors can also be expressed in terms of the in-plane material spin in the reference surface:

d^{\nabla} t_{γ} = d Ω_{γ β} t_{β} = \frac{1}{2} (t_{β} \cdot \frac{\partial d \bar{x}}{\partial s_{γ}} - t_{γ} \cdot \frac{\partial d \bar{x}}{\partial s_{β}}) t_{β} .

Substitution in the last obtained expression for $d^{\nabla} δ {\bar{ε}}_{α β}$ yields

d^{\nabla} δ {\bar{ε}}_{α β} = sym (\frac{\partial d \bar{x}}{\partial s_{α}} \cdot \frac{\partial δ \bar{x}}{\partial s_{β}} - 2 d {\bar{ε}}_{α γ} δ {\bar{ε}}_{γ β}) .

This expression is identical to the one obtained with “standard” continuum elements.

Second variation of the curvature

We need to calculate the second variation of the curvature to calculate the initial stress contribution from the curvature:

m^{α β} d^{\nabla} δ κ_{α β} .

To simplify the computation, we rely on the intrinsic definition of curvature and express the curvature in derivatives with respect to the isoparametric coordinates. Accordingly,

m^{α β} d^{\nabla} δ κ_{α β} = m^{a b} d^{\nabla} δ κ_{a b} = m^{a b} d^{\nabla} (\frac{\partial δ \bar{x}}{\partial ξ^{a}} \cdot \frac{\partial t_{3}}{\partial ξ^{b}} + \frac{\partial \bar{x}}{\partial ξ^{a}} \cdot \frac{\partial δ t_{3}}{\partial ξ^{b}}),

where the bending resultant components $m^{a b}$ are the components expressed in the orthonormal coordinate system ( $m^{α β}$ ) transformed by $\partial ξ^{a} / \partial s_{α}$ .

Denoting derivatives with respect to the isoparametric coordinates as $\partial / \partial ξ^{a} \equiv _{, a}$ , the second variation of the curvature is

d^{\nabla} δ κ_{a b} = δ {\bar{x}}_{, a} \cdot d^{\nabla} t_{3, b} + d {\bar{x}}_{, a} \cdot δ t_{3, b} + {\bar{x}}_{, a} \cdot {(d^{\nabla} δ t_{3})}_{, b} .

Using the fact that $δ t_{3} = δ ϕ \times t_{3}$ and $d^{\nabla} t_{3} = \underset{\cdot}{ϕ} \times t_{3}$ , we find that

\begin{matrix} d^{\nabla} δ κ_{a b} & ​ = δ {\bar{x}}_{, a} \cdot [- {\hat{t}}_{3}] \cdot {\underset{\cdot}{ϕ}}_{, b} + δ {\bar{x}}_{, a} \cdot [- {\hat{t}}_{3, b}] \cdot \underset{\cdot}{ϕ} + δ ϕ_{, b} \cdot [{\hat{t}}_{3}] \cdot d {\bar{x}}_{, a} + δ ϕ \cdot [{\hat{t}}_{3, b}] \cdot d {\bar{x}}_{, a} \\ + δ ϕ_{, b} \cdot (t_{3} {\bar{x}}_{, a}) \cdot \underset{\cdot}{ϕ} + δ ϕ \cdot (t_{3, b} {\bar{x}}_{, a}) \cdot \underset{\cdot}{ϕ} + δ ϕ \cdot (t_{3} {\bar{x}}_{, a}) \cdot {\underset{\cdot}{ϕ}}_{, b} \\ - δ ϕ_{, b} \cdot (t_{3} \cdot {\bar{x}}_{, a} I) \cdot \underset{\cdot}{ϕ} - δ ϕ \cdot (t_{3} \cdot {\bar{x}}_{, a} I) \cdot {\underset{\cdot}{ϕ}}_{, b} - δ ϕ \cdot (t_{3, b} \cdot {\bar{x}}_{, a} I) \cdot \underset{\cdot}{ϕ} . \end{matrix}

Here $[{\hat{t}}_{3}]$ indicates the skew-symmetric tensor with axial vector $t_{3}$ .

Transverse shear treatment

Several interpolation schemes have been proposed to avoid shear-locking, which typically arises as the thickness of a plate or shell goes to zero. Here we employ an assumed strain method based on the Hu-Washizu principle. This scheme derives from that by MacNeal (1978), subsequently extended and reformulated in Hughes and Tezduyar (1981) and MacNeal (1982) and revisited in Bathe and Dvorkin (1984). Computational aspects of the nonlinear theory are investigated in Simo, Fox, and Rifai (1989) for fully integrated quadrilateral shell elements. For reduced integration quadrilateral and triangular shell elements that can be used for both implicit and explicit integration, this assumed strain method needs to be modified. We summarize below the assumed strain method used with fully integrated elements, followed by the modifications required for the one-point integration plus stabilization used in Abaqus.

Construction of the assumed strain field

Consider a typical isoparametric finite element, as depicted in Figure 1, and denote by $A, B, C, D$ the set of midpoints of the element boundaries.

The following assumed transverse shear strain field is used:

\begin{matrix} {\bar{γ}}_{1} & ​ = \frac{1}{2} [(1 - η) γ_{1}^{B} + (1 + η) γ_{1}^{D}] \\ {\bar{γ}}_{2} & ​ = \frac{1}{2} [(1 - ξ) γ_{2}^{A} + (1 + ξ) γ_{2}^{C}], \end{matrix}

where

\begin{matrix} γ_{2}^{A} = t^{A} \cdot {\bar{x}}_{, 2}^{A} - T^{A} \cdot {\bar{X}}_{, 2}^{A}, & γ_{1}^{B} = t^{B} \cdot {\bar{x}}_{, 1}^{B} - T^{B} \cdot {\bar{X}}_{, 1}^{B}, \\ γ_{2}^{C} = t^{C} \cdot {\bar{x}}_{, 2}^{C} - T^{C} \cdot {\bar{X}}_{, 2}^{C}, & γ_{1}^{D} = t^{D} \cdot {\bar{x}}_{, 1}^{D} - T^{D} \cdot {\bar{X}}_{, 1}^{D} \end{matrix}

are the covariant transverse shear strains evaluated at the midpoints of the element boundaries. In the above transverse shear strain definitions, the use of uppercase letters indicates quantities in the reference configuration and the use of lowercase letters indicates the deformed configuration. For readability we have omitted the subscript 3 from the director field. Making use of the bilinear element interpolation, it follows that

\begin{matrix} {\bar{x}}_{, 2}^{A} = \frac{1}{2} ({\bar{x}}_{4} - {\bar{x}}_{1}), & {\bar{x}}_{, 1}^{B} = \frac{1}{2} ({\bar{x}}_{2} - {\bar{x}}_{1}), \\ {\bar{x}}_{, 2}^{C} = \frac{1}{2} ({\bar{x}}_{3} - {\bar{x}}_{2}), & {\bar{x}}_{, 1}^{D} = \frac{1}{2} ({\bar{x}}_{3} - {\bar{x}}_{4}), \end{matrix}

where ${\bar{x}}_{I}$ , for $I = 1, 2, 3, 4$ , are the reference surface position vectors of the element nodes.

By making use of the assumed strain field along with the update formulae for the director field, the assumed covariant transverse shear field can be written concisely in matrix notation. Recall the director field update equation and the corresponding linearized director field:

t_{t + Δ t} = \exp [\hat{Δ ϕ}] \cdot t_{t} \Rightarrow δ t = δ ϕ \times t .

It follows from the element interpolation that

\begin{matrix} δ t^{A} = \frac{1}{2} (δ ϕ_{4} + δ ϕ_{1}) \times t^{A}, & δ t^{B} = \frac{1}{2} (δ ϕ_{2} + δ ϕ_{1}) \times t^{B}, \\ δ t^{C} = \frac{1}{2} (δ ϕ_{3} + δ ϕ_{2}) \times t^{C}, & δ t^{D} = \frac{1}{2} (δ ϕ_{3} + δ ϕ_{4}) \times t^{D} . \end{matrix}

Define the following vectors:

\bar{γ} = {\begin{matrix} {\bar{γ}}_{1} \\ {\bar{γ}}_{2} \end{matrix}}, δ \bar{x} = {\begin{matrix} δ {\bar{x}}_{1} \\ δ {\bar{x}}_{2} \\ δ {\bar{x}}_{3} \\ δ {\bar{x}}_{4} \end{matrix}}, δ ϕ = {\begin{matrix} δ ϕ_{1} \\ δ ϕ_{2} \\ δ ϕ_{3} \\ δ ϕ_{4} \end{matrix}} .

Then, the linearized transverse shear strain is

δ \bar{γ} \equiv {\begin{matrix} δ {\bar{γ}}_{1} \\ δ {\bar{γ}}_{2} \end{matrix}} = {\bar{B}}_{s m} δ \bar{x} + {\bar{B}}_{s b} δ ϕ,

where

{\bar{B}}_{s m} = \frac{1}{4} [\begin{matrix} - (1 - η) t^{B^{T}} & (1 - η) t^{B^{T}} & (1 + η) t^{D^{T}} & - (1 + η) t^{D^{T}} \\ - (1 - ξ) t^{A^{T}} & - (1 + ξ) t^{C^{T}} & (1 + ξ) t^{C^{T}} & (1 - ξ) t^{A^{T}} \end{matrix}] .

Define the four vectors:

\begin{matrix} η_{2}^{A} = t^{A} \times {\bar{x}}_{, 2}^{A}, & η_{1}^{B} = t^{B} \times {\bar{x}}_{, 1}^{B}, \\ η_{2}^{C} = t^{C} \times {\bar{x}}_{, 2}^{C}, & η_{1}^{D} = t^{D} \times {\bar{x}}_{, 1}^{D} . \end{matrix}

Then the rotation or bending part of the strain/displacement operator is written

{\bar{B}}_{s b} = \frac{1}{4} [\begin{matrix} (1 - η) η_{1}^{B}^{T} & (1 - η) η_{1}^{B}^{T} & (1 + η) η_{1}^{D}^{T} & (1 + η) η_{1}^{D}^{T} \\ (1 - ξ) η_{2}^{A}^{T} & (1 + ξ) η_{2}^{C}^{T} & (1 + ξ) η_{2}^{C}^{T} & (1 - ξ) η_{2}^{A}^{T} \end{matrix}] .

Constitutive relations

A St. Venant-Kirchhoff constitutive model for the Kirchhoff curvilinear components of the resultant transverse shear force is written in terms of the transverse shear strains as

{\begin{matrix} Q^{1} \\ Q^{2} \end{matrix}} = C_{s} {\begin{matrix} {\bar{γ}}_{1} \\ {\bar{γ}}_{2} \end{matrix}},

where $C_{s}$ is the transverse shear stiffness in curvilinear coordinates. For a single isotropic layer,

C_{s} = \frac{5}{6} G_{s} h [\begin{matrix} A^{11} & A^{12} \\ A^{21} & A^{22} \end{matrix}] .

The matrix $[A^{α β}]$ is the inverse of the metric $[A_{α β}]$ , where metric components in the reference configuration $A_{α β}$ are defined by the inner product

A_{α β} = {\bar{X}}_{, α} \cdot {\bar{X}}_{, β} .

The Cauchy or true transverse shear force components in the shell orthonormal coordinate system ${q^{1}, q^{2}}^{T}$ are calculated with the coordinate transformation $f_{α}^{a} = \partial s^{a} / \partial ξ^{α}$ as

q = {\begin{matrix} q^{1} \\ q^{2} \end{matrix}} = \frac{A}{a} [\begin{matrix} f_{1}^{1} & f_{2}^{1} \\ f_{1}^{2} & f_{2}^{2} \end{matrix}] {\begin{matrix} Q^{1} \\ Q^{2} \end{matrix}},

where A is the element's reference area and a is the current area.

Initial stress stiffness

The calculation of the initial stress stiffness matrix requires the second variation of the assumed transverse strain field. This calculation can be summarized in matrix notation as follows. Define vectors of variations of the nodal displacement quantities:

δ u = {\begin{matrix} δ {\bar{x}}_{1} \\ δ ϕ_{1} \\ δ {\bar{x}}_{2} \\ δ ϕ_{2} \\ δ {\bar{x}}_{3} \\ δ ϕ_{3} \\ δ {\bar{x}}_{4} \\ δ ϕ_{4} \end{matrix}}, d u = {\begin{matrix} d {\bar{x}}_{1} \\ d ϕ_{1} \\ d {\bar{x}}_{2} \\ d ϕ_{2} \\ d {\bar{x}}_{3} \\ d ϕ_{3} \\ d {\bar{x}}_{4} \\ d ϕ_{4} \end{matrix}} .

Then the initial stress contribution is written

q \cdot d δ \bar{γ} = \int δ u \cdot K_{s} \cdot Δ u d a,

where $d a$ is the area measure in the current configuration and $K_{s}$ is the (symmetric) transverse shear contribution to the initial stress, defined as follows. Let $I$ be the $3 \times 3$ identity matrix; then define the symmetric matrices

\begin{matrix} A = Q^{2} (1 - ξ) [sym {t^{A} {({\bar{x}}_{, 2}^{A})}^{T}} - γ_{2}^{A} I], & B = Q^{1} (1 - η) [sym {t^{B} {({\bar{x}}_{, 1}^{B})}^{T}} - γ_{1}^{B} I], \\ C = Q^{2} (1 + ξ) [sym {t^{C} {({\bar{x}}_{, 2}^{C})}^{T}} - γ_{2}^{C} I], & D = Q^{1} (1 + η) [sym {t^{D} {({\bar{x}}_{, 1}^{D})}^{T}} - γ_{1}^{D} I] . \end{matrix}

Also define the skew-symmetric matrices

\begin{matrix} \hat{a} = Q^{2} (1 - ξ) {\hat{t}}^{A}, & \hat{b} = Q^{1} (1 - η) {\hat{t}}^{B}, \\ \hat{c} = Q^{2} (1 + ξ) {\hat{t}}^{C}, & \hat{d} = Q^{1} (1 + η) {\hat{t}}^{D} . \end{matrix}

Also, let $0$ be the $3 \times 3$ zero matrix. Then $K_{s}$ is written

K_{s} = \frac{1}{8} [\begin{matrix} 0 & \hat{a} + \hat{b} & 0 & \hat{b} & 0 & 0 & \hat{a} \\ - \hat{a} - \hat{b} & A + B & \hat{b} & B & 0 & 0 & \hat{a} & A \\ 0 & - \hat{b} & 0 & - \hat{b} + \hat{c} & 0 & \hat{c} & 0 & 0 \\ - \hat{b} & B & \hat{b} - \hat{c} & B + C & \hat{c} & C & 0 & 0 \\ 0 & 0 & 0 & - \hat{c} & 0 & - \hat{c} - \hat{d} & 0 & - \hat{d} \\ 0 & 0 & - \hat{c} & C & \hat{c} + \hat{d} & C + D & - \hat{d} & D \\ 0 & - \hat{a} & 0 & 0 & 0 & \hat{d} & 0 & \hat{d} - \hat{a} \\ - \hat{a} & A & 0 & 0 & \hat{d} & D & - \hat{d} + \hat{a} & D + A \end{matrix}] .

One point integration plus stabilization

For reduced-integration elements the transverse shear force components need to be evaluated at the center of the elements. Consider $π_{s}$ the transverse shear contribution to the internal energy:

π_{s} = \frac{1}{2} \int \bar{γ} \cdot C_{s} \cdot \bar{γ} d A .

The reference area measure $d A$ is written in terms of the isoparametric coordinates as $d A = J d ξ d η$ , where $J = \sqrt{A_{11} A_{22} - {(A_{12})}^{2}}$ and $A_{α β}$ are the components of the reference surface metric in the undeformed configuration.

This transverse shear energy can be approximated in many ways to produce a one point integration at the center of the element plus hourglass stabilization. It is important that this treatment yield accurate representation of transverse shear deformation in thick shell problems and provide robust performance for skewed elements. The treatment should collapse smoothly to a triangle, which should be insensitive to the node numbering during collapse; that is, the triangle's response should not depend on the nodal connectivity. For an entire mesh of triangular elements, the treatment should give convergent results (that is, the element should not lock). Furthermore, the high frequency response of the transverse shear treatment should be controlled so that transverse shear response does not dominate the stable time increment for explicit dynamic analysis (including for skewed triangular or quadrilateral geometries). All of these requirements are embodied in the following transverse shear treatment.

Define the transverse shear strain at the center of the element (the homogeneous part) and the “hourglass” transverse shear strain vectors as

{\bar{γ}}_{0} = \frac{1}{2} {\begin{matrix} γ_{1}^{B} + γ_{1}^{D} \\ γ_{2}^{A} + γ_{2}^{C} \end{matrix}} + γ_{c c} {\begin{matrix} c_{ξ} \\ c_{η} \end{matrix}} and γ_{h g} = {\begin{matrix} γ_{b f} \\ γ_{c c} \end{matrix}} .

The element distortion coefficients $c_{ξ}$ and $c_{η}$ are constants determined by the element reference geometry. For geometries with constant Jacobian transformation, $c_{ξ} = c_{η} = 0$ . The components of the hourglass strain vector $γ_{h g}$ are defined in terms of the edge strains as

γ_{b f} = - Γ_{A} γ_{2}^{A} + Γ_{B} γ_{1}^{B} + Γ_{C} γ_{2}^{C} - Γ_{D} γ_{1}^{D} and γ_{c c} = - γ_{2}^{A} + γ_{1}^{B} + γ_{2}^{C} - γ_{1}^{D} .

The coefficients $Γ_{A}$ , $Γ_{B}$ , $Γ_{C}$ , and $Γ_{D}$ are constants determined from the reference geometry of the element. For rectangular elements $Γ_{A} = \frac{1}{2}$ , $Γ_{B} = - \frac{1}{2}$ , $Γ_{C} = \frac{1}{2}$ , $Γ_{D} = - \frac{1}{2},$ and $γ_{b f}$ can be identified as the strain associated with the rotational “butterfly” deformation pattern. We call $γ_{c c}$ the “crop circle” mode strain since it corresponds to a deformation pattern that resembles the sweeping over the element normals in a circular pattern.

The inclusion of the crop circle strain $γ_{c c}$ in the homogeneous part of the transverse shear strain ${\bar{γ}}_{0}$ has two important consequences. First, it makes the transverse shear response insensitive to the nodal connectivity for a triangular element. That is, when a side of a quadrilateral element is collapsed to form a triangle, the element's response is independent of the choice of node numbering on the element. Second, for explicit dynamic analyses the coefficients $c_{ξ}$ and $c_{η}$ are chosen to minimize the highest frequencies associated with the homogeneous part of the transverse shear response.

To illustrate the crop circle and butterfly transverse shear patterns, consider a square, initially flat element. Furthermore, consider plate theory kinematics; that is, two rotations and a vertical deflection at the nodes. The crop circle pattern has zero vertical deflection at the nodes and a nodal rotation vector pattern as illustrated in Figure 2.

The butterfly pattern has vertical deflections that correspond to cross-diagonal bending; that is, two equal deflections at two nodes across a diagonal, with equal and opposite deflections at the remaining two nodes. The nodal rotations develop in a way that opposes the bending motion of the reference surface; that is, the rotations are opposite the rotations that would develop for this displacement pattern to produce pure bending. The butterfly mode's nodal vertical deflection and rotation vector pattern are illustrated in Figure 3.

Let the reference element area be $A_{0} = 4 J_{0}$ . The transverse shear energy can be approximated as a center point value plus a stabilization term:

π_{s} = \frac{1}{2} A_{0} {\bar{γ}}_{0} \cdot C_{s}_{0} \cdot {\bar{γ}}_{0} + \frac{1}{2} A_{0} γ_{h g} \cdot H \cdot γ_{h g},

where $C_{s}_{0}$ is the transverse shear stiffness evaluated at the center of the element and the hourglass stiffness $H$ is the diagonal matrix

H = \frac{C_{s 0}^{e f f}}{12} [\begin{matrix} 1 & 0 \\ 0 & 0.001 \end{matrix}] .

The effective stiffness $C_{s 0}^{e f f}$ is the average direct component of the transverse shear stiffness, $C_{s 0}^{e f f} = (C_{s 0}^{11} + C_{s 0}^{22}) / 2$ .

The formulation of the homogeneous part of the transverse shear has two contributions: the average edge strain across the element, plus the element distortion term. The average strain treatment is essentially the same as that for the assumed strain formulation of MacNeal and others presented earlier, with expressions evaluated at the center of the element ( $ξ = 0$ and $η = 0$ ). The details of this part are omitted; only the element distortion term is presented in detail. The variation of the homogeneous transverse shear strain can be written

δ {\bar{γ}}_{0} = {\bar{B}}_{s m 0} δ \bar{x} + {\bar{B}}_{s b 0} δ ϕ + {\bar{B}}_{c c d} δ \bar{x} + {\bar{B}}_{c c r} δ ϕ,

where ${\bar{B}}_{s m 0}$ and ${\bar{B}}_{s b 0}$ are ${\bar{B}}_{s m}$ and ${\bar{B}}_{s b}$ evaluated at the center of the element,

{\bar{B}}_{c c d} = \frac{1}{2} [\begin{matrix} c_{ξ} {(t^{A} - t^{B})}^{T} & c_{ξ} {(t^{B} - t^{C})}^{T} & c_{ξ} {(t^{C} - t^{D})}^{T} & c_{ξ} {(t^{D} - t^{A})}^{T} \\ c_{η} {(t^{A} - t^{B})}^{T} & c_{η} {(t^{B} - t^{C})}^{T} & c_{η} {(t^{C} - t^{D})}^{T} & c_{η} {(t^{D} - t^{A})}^{T} \end{matrix}]

and

{\bar{B}}_{c c r} = \frac{1}{2} [\begin{matrix} c_{ξ} {(η_{1}^{B} - η_{2}^{A})}^{T} & c_{ξ} {(η_{1}^{B} + η_{2}^{C})}^{T} & c_{ξ} {(η_{2}^{C} - η_{1}^{D})}^{T} & c_{ξ} {(- η_{2}^{A} - η_{1}^{D})}^{T} \\ c_{η} {(η_{1}^{B} - η_{2}^{A})}^{T} & c_{η} {(η_{1}^{B} + η_{2}^{C})}^{T} & c_{η} {(η_{2}^{C} - η_{1}^{D})}^{T} & c_{η} {(- η_{2}^{A} - η_{1}^{D})}^{T} \end{matrix}] .

The stabilization term has a similar formulation. The variation of the hourglass strain is

δ γ_{h g} = {\bar{B}}_{d} δ \bar{x} + {\bar{B}}_{r} δ ϕ,

where

{\bar{B}}_{d} = \frac{1}{2} [\begin{matrix} {(Γ_{A} t^{A} - Γ_{B} t^{B})}^{T} & {(Γ_{B} t^{B} - Γ_{C} t^{C})}^{T} & {(Γ_{C} t^{C} - Γ_{D} t^{D})}^{T} & {(Γ_{D} t^{D} - Γ_{A} t^{A})}^{T} \\ {(t^{A} - t^{B})}^{T} & {(t^{B} - t^{C})}^{T} & {(t^{C} - t^{D})}^{T} & {(t^{D} - t^{A})}^{T} \end{matrix}]

and

{\bar{B}}_{r} = \frac{1}{2} [\begin{matrix} {(Γ_{B} η_{1}^{B} - Γ_{A} η_{2}^{A})}^{T} & {(Γ_{B} η_{1}^{B} + Γ_{C} η_{2}^{C})}^{T} & {(Γ_{C} η_{2}^{C} - Γ_{D} η_{1}^{D})}^{T} & {(- Γ_{A} η_{2}^{A} - Γ_{D} η_{1}^{D})}^{T} \\ {(η_{1}^{B} - η_{2}^{A})}^{T} & {(η_{1}^{B} + η_{2}^{C})}^{T} & {(η_{2}^{C} - η_{1}^{D})}^{T} & {(- η_{2}^{A} - η_{1}^{D})}^{T} \end{matrix}] .

The hourglass force components $h^{1}$ and $h^{2}$ are given by the constitutive relations

h^{1} = \frac{C_{s 0}^{e f f}}{12} γ_{b f} and h^{2} = 0.001 \frac{C_{s 0}^{e f f}}{12} γ_{c c} .

Comments on stabilization

(1) The butterfly mode $γ_{b f}$ is applied with a “large” or physical hourglass stiffness. For a reference geometry with constant Jacobian, the butterfly stabilization term can be derived from an exact integration of the assumed strain formulation of the transverse shear energy. It is important to apply this constraint with a high stiffness to prevent overly flexible response for quadrilateral elements. The crop circle mode is applied with a “small” or weak stiffness. Although this mode can propagate, it is rarely problematic and is often prevented with boundary conditions.

(2) As the quadrilateral element is degenerated to a triangle, the two hourglass constraints converge into a single constraint: the crop circle constraint. However, as is well-known, for a constant strain triangle the element will lock for certain meshes with three transverse shear constraints per element. Therefore, in the case of a triangular element, the (strong) butterfly mode stabilization is not applied. Only the (weak) crop circle mode stabilization is applied. Thus, in addition to the two homogeneous transverse shear strains, the triangle has a weak constraint to prevent spurious zero energy modes, yet avoids locking in most situations.

The initial stress contribution from the stabilization terms takes the following form:

A_{0} (h^{1} d δ γ_{b f} + h^{2} d δ γ_{c c}) = A_{0} δ u \cdot K_{h} \cdot d u,

where $K_{h}$ is the (symmetric) transverse shear stabilization contribution to the initial stress. Define the symmetric matrices

\begin{matrix} H^{A} = (h^{1} Γ_{A} + h^{2}) [sym {t^{A} {({\bar{x}}_{, 2}^{A})}^{T}} - γ_{2}^{A} I], & H^{B} = (h^{1} Γ_{B} + h^{2}) [sym {t^{B} {({\bar{x}}_{, 1}^{B})}^{T}} - γ_{1}^{B} I], \\ H^{C} = (h^{1} Γ_{C} + h^{2}) [sym {t^{C} {({\bar{x}}_{, 2}^{C})}^{T}} - γ_{2}^{C} I], & H^{D} = (h^{1} Γ_{D} + h^{2}) [sym {t^{D} {({\bar{x}}_{, 1}^{D})}^{T}} - γ_{1}^{D} I] . \end{matrix}

Also define the skew-symmetric matrices

\begin{matrix} {\hat{a}}_{h} = (h^{1} Γ_{A} + h^{2}) {\hat{t}}^{A}, & {\hat{b}}_{h} = (h^{1} Γ_{B} + h^{2}) {\hat{t}}^{B}, \\ {\hat{c}}_{h} = (h^{1} Γ_{C} + h^{2}) {\hat{t}}^{C}, & {\hat{d}}_{h} = (h^{1} Γ_{D} + h^{2}) {\hat{t}}^{D} . \end{matrix}

Then $K_{h}$ is written

\begin{matrix} K_{h} = ​ \\ \frac{1}{4} [\begin{matrix} 0 & - {\hat{a}}_{h} + {\hat{b}}_{h} & 0 & {\hat{b}}_{h} & 0 & 0 & 0 & - {\hat{a}}_{h} \\ {\hat{a}}_{h} - {\hat{b}}_{h} & - H^{A} + H^{B} & {\hat{b}}_{h} & H^{B} & 0 & 0 & - {\hat{a}}_{h} & - H^{A} \\ 0 & - {\hat{b}}_{h} & 0 & - {\hat{b}}_{h} + {\hat{c}}_{h} & 0 & {\hat{c}}_{h} & 0 & 0 \\ - {\hat{b}}_{h} & H^{B} & {\hat{b}}_{h} - {\hat{c}}_{h} & H^{B} + H^{C} & {\hat{c}}_{h} & H^{C} & 0 & 0 \\ 0 & 0 & 0 & - {\hat{c}}_{h} & 0 & - {\hat{c}}_{h} + {\hat{d}}_{h} & 0 & {\hat{d}}_{h} \\ 0 & 0 & - {\hat{c}}_{h} & H^{C} & {\hat{c}}_{h} - {\hat{d}}_{h} & H^{C} - H^{D} & {\hat{d}}_{h} & - H^{D} \\ 0 & {\hat{a}}_{h} & 0 & 0 & 0 & - {\hat{d}}_{h} & 0 & - {\hat{d}}_{h} + {\hat{a}}_{h} \\ {\hat{a}}_{h} & - H^{A} & 0 & 0 & - {\hat{d}}_{h} & - H^{D} & {\hat{d}}_{h} - {\hat{a}}_{h} & - H^{D} - H^{A} \end{matrix}] . \end{matrix}

Note that once the matrix entries in $K_{h}$ are defined, $K_{h}$ is filled just as $K_{s}$ .

The initial stress contribution from the homogeneous part consists of two terms, one from the assumed strain formulation (evaluated at $ξ = η = 0$ ) as detailed earlier, and the other from the crop circle mode addition. These two terms can be written

{A_{0} q |}_{0} \cdot d δ {\bar{γ}}_{0} = {A_{0} δ u \cdot K_{s} |}_{0} \cdot d u + A_{0} (q^{1} c_{ξ} + q^{2} c_{η}) d δ γ_{c c},

where ${q |}_{0}$ and ${K_{s} |}_{0}$ are the shear force and matrix $K_{s}$ evaluated at the element center. The matrix expression for $(q^{1} c_{ξ} + q^{2} c_{η}) d δ γ_{c c}$ is analogous to $K_{h}$ from the stabilization terms.

In-plane displacement hourglass control

The in-plane displacement hourglass control is applied in the same way as in the Abaqus membrane elements. The hourglass strains are defined by

z_{α} = \frac{\partial \bar{x}}{\partial S_{α}} \cdot {\bar{x}}^{I} γ^{I} - \frac{\partial \bar{X}}{\partial S_{α}} \cdot {\bar{X}}^{I} γ^{I} = \frac{\partial \bar{x}}{\partial S_{α}} \cdot {\bar{x}}^{I} γ^{I} - T_{α} \cdot {\bar{X}}^{I} γ^{I},

where $γ^{I}$ is the hourglass mode. This mode is obtained by making the “regular” hourglass mode $Γ^{I} = {1, - 1, 1, - 1}$ orthogonal to the homogeneous deformation mode in the undeformed shape of the element. This last condition can be written as

γ^{I} = Γ^{I} - Γ^{J} X^{J} \cdot T_{β} \frac{\partial N^{I}}{\partial S_{β}} .

Observe that

T_{α} \cdot {\bar{X}}^{I} γ^{I} = T_{α} \cdot {\bar{X}}^{I} Γ^{I} - Γ^{J} {\bar{X}}^{J} \cdot T_{β} \frac{\partial \bar{X}}{\partial S_{β}} \cdot T_{α} = T_{α} \cdot {\bar{X}}^{I} Γ^{I} - Γ^{J} {\bar{X}}^{J} \cdot T_{α} = 0,

and consequently

z_{α} = \frac{\partial x}{\partial S_{α}} \cdot x^{I} γ^{I} .

This expression can be worked out further. We define the projected nodal coordinates

S_{α}^{I} \overset{def}{=} T_{α} \cdot {\bar{X}}^{I}

and the projected element area

A = \frac{1}{2} [(S_{1}^{3} - S_{1}^{1}) (S_{2}^{4} - S_{2}^{2}) + (S_{1}^{2} - S_{1}^{4}) (S_{2}^{3} - S_{2}^{1})] .

The hourglass mode can then be written in the form

\begin{matrix} γ^{1} = ​ & \frac{1}{A} [S_{1}^{2} (S_{2}^{3} - S_{2}^{4}) + S_{1}^{3} (S_{2}^{4} - S_{2}^{2}) + S_{1}^{4} (S_{2}^{2} - S_{2}^{3})], \\ γ^{2} = ​ & \frac{1}{A} [S_{1}^{3} (S_{2}^{1} - S_{2}^{4}) + S_{1}^{4} (S_{2}^{3} - S_{2}^{1}) + S_{1}^{1} (S_{2}^{4} - S_{2}^{3})], \\ γ^{3} = ​ & \frac{1}{A} [S_{1}^{4} (S_{2}^{1} - S_{2}^{2}) + S_{1}^{1} (S_{2}^{2} - S_{2}^{4}) + S_{1}^{2} (S_{2}^{4} - S_{2}^{1})], \\ γ^{4} = ​ & \frac{1}{A} [S_{1}^{1} (S_{2}^{3} - S_{2}^{2}) + S_{1}^{2} (S_{2}^{1} - S_{2}^{3}) + S_{1}^{3} (S_{2}^{2} - S_{2}^{1})] . \end{matrix}

The hourglass stiffness is chosen equal to

K_{h} = (r_{F} G) \frac{\partial N^{I}}{\partial S_{β}} \frac{\partial N^{I}}{\partial S_{β}} h A,

where G is the shear modulus and $r_{F}$ is a small number chosen to be 0.005 in Abaqus/Standard and 0.05 in Abaqus/Explicit. When the hourglass control is based on assumed enhanced strain, the artificial stiffness factor is replaced by coefficients derived from a three-field variational principle. The hourglass force Z conjugate to z is then equal to

Z_{α} = K_{h} z_{α} .

For virtual work we need the first variation of the hourglass strain. From the expression for the strain follows immediately

δ z_{α} = \frac{\partial \bar{x}}{\partial S_{α}} \cdot δ {\bar{x}}^{I} γ^{I} + \frac{\partial δ \bar{x}}{\partial S_{α}} \cdot {\bar{x}}^{I} γ^{I} .

Note that the second term vanishes in the initial configuration since ${\bar{X}}^{I} γ^{I} = 0$ . The second variation is needed for the Jacobian. From the first variation follows right away

d δ z_{α} = \frac{\partial d \bar{x}}{\partial S_{α}} \cdot δ {\bar{x}}^{I} γ^{I} + \frac{\partial δ \bar{x}}{\partial S_{α}} \cdot d {\bar{x}}^{I} γ^{I} .

The second variation does not contribute in the initial configuration since initially $Z_{α} = 0$ .

Rotational hourglass control

The expressions for the curvature change, the transverse shear constraints, and the drilling mode constraints still leave three nonhomogeneous rotational modes unconstrained. These modes correspond to zero rotation at the midedges and zero gradient at the centroid. Hence, they correspond to the familiar $Γ^{I} = {1, - 1, 1, - 1}$ hourglass pattern. To pass curvature patch tests exactly, it is necessary to use orthogonalized hourglass patterns as derived for in-plane hourglass control.

This last aspect implies that the rotational hourglass mode corresponds to the mixed derivative of the rotation at the centroid:

γ^{I} Δ ϕ^{I} = 4 \times \frac{\partial^{2} Δ ϕ}{\partial ξ_{1} \partial ξ_{2}} .

We cannot use the above formulation directly in a formulation suitable for multiple finite rotation increments. Hence, we use the same approach as for the calculation of the curvature change. For the purpose of the calculation we define the updated shell direction vectors

{\tilde{t}}_{i}^{t + Δ t} = Δ q t_{i}^{t} Δ q^{†} .

The updated shell direction vectors do not actually have to be stored: they are used only for the derivation of the expression for the hourglass strain. We now formally define the hourglass strain tensor as

ϑ_{i j} = t_{i} \cdot \frac{\partial^{2} t_{j}}{\partial ξ_{1} \partial ξ_{2}} .

Observe that

\frac{\partial^{2} (t_{i} \cdot t_{j})}{\partial ξ_{1} \partial ξ_{2}} = t_{i} \cdot \frac{\partial^{2} t_{j}}{\partial ξ_{1} \partial ξ_{2}} + \frac{\partial^{2} t_{i}}{\partial ξ_{1} \partial ξ_{2}} \cdot t_{j} + \frac{\partial t_{i}}{\partial ξ_{1}} \cdot \frac{\partial t_{j}}{\partial ξ_{2}} + \frac{\partial t_{i}}{\partial ξ_{2}} \cdot \frac{\partial t_{j}}{\partial ξ_{1}} = 0.

For the purpose of hourglass strain calculation we assume that all products of first-order derivatives with respect to $ξ_{1}$ and $ξ_{2}$ can be neglected. Consequently,

\frac{\partial^{2} (t_{i} \cdot t_{j})}{\partial ξ_{1} \partial ξ_{2}} = t_{i} \cdot \frac{\partial^{2} t_{j}}{\partial ξ_{1} \partial ξ_{2}} + \frac{\partial^{2} t_{i}}{\partial ξ_{1} \partial ξ_{2}} \cdot t_{j} = 0;

and, hence,

ϑ_{i j} = - ϑ_{j i}

is skew-symmetric. Observe that the mixed derivative of $t_{i}$ can be expressed in terms of the hourglass strain tensor with

t_{i} ϑ_{i j} = t_{i} t_{i} \cdot \frac{\partial^{2} t_{j}}{\partial ξ_{1} \partial ξ_{2}} = \frac{\partial^{2} t_{j}}{\partial ξ_{1} \partial ξ_{2}} .

In the undeformed configuration, we assume that $ϑ_{i j} = 0$ . Subsequent values of $ϑ_{i j}$ are obtained incrementally. From the expression for ${\tilde{t}}_{i}^{t + Δ t}$ we obtain

\begin{matrix} \frac{\partial^{2} {\tilde{t}}_{i}^{t + Δ t}}{\partial ξ_{1} \partial ξ_{2}} = ​ & \frac{\partial^{2} Δ q}{\partial ξ_{1} \partial ξ_{2}} t_{i}^{t} Δ q^{†} + Δ q t_{i}^{t} \frac{\partial Δ q^{†}}{\partial ξ_{1} \partial ξ_{2}} + Δ q \frac{\partial t_{i}^{t}}{\partial ξ_{1} \partial ξ_{2}} Δ q^{†} + \\ \frac{\partial Δ q}{\partial ξ_{1}} t_{i}^{t} \frac{\partial Δ q^{†}}{\partial ξ_{2}} + \frac{\partial Δ q}{\partial ξ_{2}} t_{i}^{t} \frac{\partial Δ q^{†}}{\partial ξ_{1}} + \frac{\partial Δ q}{\partial ξ_{1}} \frac{\partial t_{i}^{t}}{\partial ξ_{2}} Δ q^{†} + \\ \frac{\partial Δ q}{\partial ξ_{2}} \frac{\partial t_{i}^{t}}{\partial ξ_{1}} Δ q^{†} + Δ q \frac{\partial t_{i}^{t}}{\partial ξ_{1}} \frac{\partial Δ q^{†}}{\partial ξ_{2}} + Δ q \frac{\partial t_{i}^{t}}{\partial ξ_{2}} \frac{\partial Δ q^{†}}{\partial ξ_{1}} . \end{matrix}

In this expression we also ignore all terms with products of derivatives with respect to $ξ_{1}$ and $ξ_{2}$ . Hence, the above expression simplifies to

\frac{\partial {\tilde{t}}_{i}^{t + Δ t}}{\partial ξ_{1} \partial ξ_{2}} = \frac{\partial^{2} Δ q}{\partial ξ_{1} \partial ξ_{2}} t_{i}^{t} Δ q^{†} + Δ q t_{i}^{t} \frac{\partial Δ q^{†}}{\partial ξ_{1} \partial ξ_{2}} + Δ q \frac{\partial t_{i}^{t}}{\partial ξ_{1} \partial ξ_{2}} Δ q^{†} .

The second term on the right-hand side can be written in the form

Δ q t_{i}^{t} \frac{\partial Δ q^{†}}{\partial ξ_{1} \partial ξ_{2}} = {[\frac{\partial Δ q}{\partial ξ_{1} \partial ξ_{2}} {(t_{i}^{t})}^{†} Δ q^{†}]}^{†} = - {[\frac{\partial Δ q}{\partial ξ_{1} \partial ξ_{2}} t_{i}^{t} Δ q^{†}]}^{†} .

Hence, the scalar parts of the first two terms cancel each other and the vector parts reinforce each other, leading to

\frac{\partial {\tilde{t}}_{i}^{t + Δ t}}{\partial ξ_{1} \partial ξ_{2}} = 2 V (\frac{\partial Δ q}{\partial ξ_{1} \partial ξ_{2}} t_{i}^{t} Δ q^{†}) + Δ q \frac{\partial t_{i}^{t}}{\partial ξ_{1} \partial ξ_{2}} Δ q^{†} .

The inverse of a rotation quaternion such as $Δ q$ is equal to its conjugate ( $Δ q^{- 1} = Δ q^{†}$ ); hence, we can write

\begin{matrix} \frac{\partial {\tilde{t}}_{i}^{t + Δ t}}{\partial ξ_{1} \partial ξ_{2}} & ​ = 2 V (\frac{\partial Δ q}{\partial ξ_{1} \partial ξ_{2}} Δ q^{†} Δ q t_{i}^{t} Δ q^{†}) + Δ q \frac{\partial t_{i}^{t}}{\partial ξ_{1} \partial ξ_{2}} Δ q^{†} \\ ​ = V (2 \frac{\partial Δ q}{\partial ξ_{1} \partial ξ_{2}} Δ q^{†} {\tilde{t}}_{i}^{t + Δ t}) + Δ q \frac{\partial t_{i}^{t}}{\partial ξ_{1} \partial ξ_{2}} Δ q^{†} \\ ​ = Δ R_{Γ} \times {\tilde{t}}_{i}^{t + Δ t} + Δ q \frac{\partial t_{i}^{t}}{\partial ξ_{1} \partial ξ_{2}} Δ q^{†}, \end{matrix}

where we have formally defined the incremental hourglass update vector

Δ R_{Γ} \overset{def}{=} V (2 \frac{\partial Δ q}{\partial ξ_{1} \partial ξ_{2}} Δ q^{†}),

which must be expressed in terms of the incremental rotation hourglass mode. From the definition of the incremental quaternion $Δ q$ follows, while neglecting the products of $ξ_{1}$ and $ξ_{2}$ derivatives:

\frac{\partial Δ q}{\partial ξ_{1} \partial ξ_{2}} = (\frac{1}{2} \sin \frac{Δ ϕ}{2} \frac{\partial Δ ϕ}{\partial ξ_{1} \partial ξ_{2}}, \frac{1}{2} \cos \frac{Δ ϕ}{2} p \frac{\partial Δ ϕ}{\partial ξ_{1} \partial ξ_{2}} + \sin \frac{Δ ϕ}{2} \frac{\partial p}{\partial ξ_{1} \partial ξ_{2}});

thus, for $Δ R_{Γ}$ again with use of the incremental quaternion definition

Δ R_{Γ} = p \frac{\partial Δ ϕ}{\partial ξ_{1} \partial ξ_{2}} + \sin Δ ϕ \frac{\partial p}{\partial ξ_{1} \partial ξ_{2}} + (1 - \cos Δ ϕ) p \times \frac{\partial p}{\partial ξ_{1} \partial ξ_{2}} .

From the definition of $Δ ϕ$ and $p$ follows, again neglecting the products of first derivatives

\begin{matrix} \frac{\partial Δ ϕ}{\partial ξ_{1} \partial ξ_{2}} & ​ = \frac{1}{Δ ϕ} Δ ϕ \cdot \frac{\partial Δ ϕ}{\partial ξ_{1} \partial ξ_{2}}, \\ \frac{\partial p}{\partial ξ_{1} \partial ξ_{2}} & = \frac{1}{Δ ϕ} \frac{\partial Δ ϕ}{\partial ξ_{1} \partial ξ_{2}} - \frac{1}{Δ ϕ^{3}} Δ ϕ Δ ϕ \cdot \frac{\partial Δ ϕ}{\partial ξ_{1} \partial ξ_{2}} = \frac{1}{Δ ϕ} [\frac{\partial Δ ϕ}{\partial ξ_{1} \partial ξ_{2}} - p p \cdot \frac{\partial Δ ϕ}{\partial ξ_{1} \partial ξ_{2}}] . \end{matrix}

After substitution in the expression for $Δ R_{Γ}$ and some algebra one obtains

Δ R_{Γ} = \frac{\sin Δ ϕ}{Δ ϕ} \frac{\partial Δ ϕ}{\partial ξ_{1} \partial ξ_{2}} + \frac{1 - \cos Δ ϕ}{Δ ϕ} p \times \frac{\partial Δ ϕ}{\partial ξ_{1} \partial ξ_{2}} + (1 - \frac{\sin Δ ϕ}{Δ ϕ}) p p \cdot \frac{\partial Δ ϕ}{\partial ξ_{1} \partial ξ_{2}} .

Note that $Δ R_{Γ} \to \frac{\partial Δ ϕ}{\partial ξ_{1} \partial ξ_{2}}$ when $Δ ϕ \to 0$ . For the updated hourglass tensor one readily obtains

\begin{matrix} ϑ_{i j}^{t + Δ t} & = {\tilde{t}}_{i}^{t + Δ t} \cdot \frac{\partial^{2} {\tilde{t}}_{j}^{t + Δ t}}{\partial ξ_{1} \partial ξ_{2}} = {\tilde{t}}_{i}^{t + Δ t} \cdot [Δ R_{Γ} \times {\tilde{t}}_{j}^{t + Δ t} + Δ q \frac{\partial^{2} t_{j}^{t}}{\partial ξ_{1} \partial ξ_{2}} Δ q^{†}] \\ ​ = - Δ R_{Γ} \cdot ({\tilde{t}}_{i}^{t + Δ t} \times {\tilde{t}}_{j}^{t + Δ t}) + {\tilde{t}}_{i}^{t + Δ t} \cdot (Δ q t_{k}^{t} ϑ_{k j}^{t} Δ q^{†}) \\ = - Δ R_{Γ} \cdot ({\tilde{t}}_{i}^{t + Δ t} \times {\tilde{t}}_{j}^{t + Δ t}) + {\tilde{t}}_{i}^{t + Δ t} \cdot {\tilde{t}}_{k}^{t + Δ t} ϑ_{k j}^{t} = - Δ R_{Γ} \cdot ({\tilde{t}}_{i}^{t + Δ t} \times {\tilde{t}}_{j}^{t + Δ t}) + ϑ_{i j}^{t} . \end{matrix}

This expression simplifies further with the introduction of the hourglass vector

ϑ_{i} = - \frac{1}{2} e_{i j k} ϑ_{j k},

which yields the update formula

ϑ_{i}^{t + Δ t} = Δ R_{Γ} \cdot {\tilde{t}}_{i}^{t + Δ t} + ϑ_{i}^{t} .

The first and second variation are obtained in entirely the same way as the first and second variation of the curvature change. For the first variation we neglect terms of order $Δ R_{Γ}$ and obtain

δ ϑ_{i} = δ Δ R_{Γ} \cdot t_{i} + Δ R_{Γ} \cdot δ t_{i} = δ Δ R_{Γ} \cdot t_{i} = t_{i} \cdot \frac{\partial^{2} δ ϕ}{\partial ξ_{1} \partial ξ_{2}} .

For the second variation we ignore in addition the terms of order $d R_{Γ}$ and $δ R_{Γ}$ with as final result

d δ ϑ_{i} = 0.

Degenerate elements

In general meshes it will be desirable to collapse at least some of the quadrilateral elements to triangles or to use the triangular element S3 or S3R, which is in fact an internally collapsed S4R element. For this case the calculation of the membrane strains and the curvature changes proceeds along the same lines as before. The transverse shears will now be zero at the degenerate edges. Finally, calculation of all hourglass constraints will be omitted.

Rotary inertia scaling for explicit dynamics

For numerical efficiency in explicit dynamic analysis, it is desirable to have the stable time increment determined by the membrane response of the structure. For this reason scaling of the rotary inertia based on the element's reference geometry is included in Abaqus/Explicit.

In explicit dynamic analyses the stable time increment is proportional to the inverse of the highest frequency of the element. Therefore, we must ensure that the highest frequency associated with the transverse shear response does not exceed the highest frequency associated with the membrane response. For thick elements (that is, for elements whose thickness is order unity relative to a characteristic length in the element), the membrane frequencies are dominant. The primary consideration in choosing appropriate scalings is that in the limit as the element's thickness goes to zero, the transverse shear frequencies remain below the membrane ones. Recall that for a one-dimensional spring-mass oscillator, the natural frequency $ω$ can be written in terms of the stiffness K and the mass M as

ω = \sqrt{\frac{K}{M}} .

For the transverse shear response the rotary inertia, which is proportional to the cube of the thickness, plays the role of the mass of the system. All other quantities—the membrane stiffness, the mass associated with membrane deformation, and the transverse shear stiffness—are proportional to the thickness. Hence as the thickness of the element goes to zero, the frequencies associated with transverse shear go to infinity proportional to the inverse of the thickness, while the membrane frequencies remain constant. Without scaling, the stable time increment would go to zero as the thickness becomes small.

Rotary inertia scaling

For thin elements the rotary inertia is small (negligible) relative to the rotational inertia of the mass at the nodes rotating about an axis through the element center. Therefore, we choose a scaling on the rotary inertia such that it never becomes smaller than a fixed (small) percentage of the rotational inertia of the mass at the nodes rotating about an axis through the center of the element.

Let R be the nondimensional rotary inertia scaling, where $R \geq 1$ . When $R = 1$ , the true rotary inertia is used. Consider a lumped mass matrix for a 4-node element, and let the element be flat and square. For rotations about an axis in the plane of the element, parallel to an element edge, passing through the center, the contribution to the rotational inertia of the mass at the nodes is

\frac{1}{4} ρ A h L^{2},

where A is the area of the element, L is the characteristic element length, and $ρ$ is the mass density. The sum of the rotary inertia at the four nodes is

R ρ A \frac{h^{3}}{12} .

The ratio of the in-plane contribution to the rotary contribution is

\frac{3 L^{2}}{R h^{2}} .

For the rotary inertia to remain a fixed fraction—say $ϵ$ —of the mass contribution as the thickness goes to zero, asymptotically R must be proportional to $L^{2} / h^{2}$ ; that is,

R \to ϵ \frac{3 L^{2}}{h^{2}} .

For planar geometries with element directors along the normal direction, closed-form expressions are possible for the highest membrane and transverse shear frequencies. In such cases the length parameter L is interpreted as a characteristic element length that depends on the element distortion. To handle arbitrarily shaped curved elements, exact calculation of the element frequencies becomes difficult. However, we can safely bound the frequencies by an appropriate choice of L in the following scaling:

R = {[16 + ϵ^{4} {(\frac{3 L^{2}}{h^{2}})}^{4}]}^{\frac{1}{4}},

where $ϵ = 1 / 4$ . For quadrilateral elements the characteristic length, $L,$ is,

L = \frac{2 ∥ x_{, ξ} \times x_{, η} ∥}{max (∥ x_{, ξ} ∥, ∥ x_{, η} ∥)} .

The factor 16 in the definition of R is used to protect against bending frequencies determining the stable time increment in very fine meshes subjected to loads that cause an increase in thickness of the shell.

Rotational bulk viscosity for explicit dynamics

For the displacement degrees of freedom, bulk viscosity introduces damping associated with volumetric straining. Linear bulk viscosity or truncation frequency damping is used to damp the high frequency ringing that leads to unwanted noise in the solution or spurious overshoot in the response amplitude. For the same reason, in shells we need to damp the high frequency ringing in the rotational degrees of freedom with linear bulk viscosity acting on the mean curvature strain rate. This damping generates a bulk viscosity “pressure moment,” m, which is linear in the mean curvature strain rate:

m = b \frac{h_{0}^{2}}{12} ρ c_{d} L \frac{Δ κ}{Δ t},

where b is a damping coefficient (default = 0.06), $h_{0}$ is the original thickness, $ρ$ is the mass density, $c_{d}$ is the current dilatational wave speed, L is the characteristic length used for rotary inertia scaling, and $Δ κ = Δ κ_{11} + Δ κ_{22}$ is twice the increment in mean curvature. The dilatational wave speed is given in terms of the effective Lamé's constants as

c_{d} = \sqrt{\frac{(λ + 2 μ)}{ρ}} .

The resultant pressure moment $m h$ , where h is the current thickness, is added to the direct components of the moment resultant.

Fully integrated finite-membrane-strain shell formulation

Element S4 is a fully integrated finite-membrane-strain shell element. Since the element's stiffness is fully integrated, no spurious membrane or bending zero energy modes exist and no membrane or bending mode hourglass stabilization is used. Drill rotation control, however, is required. Element S4 uses the same drill stiffness formulation as used for element S4R. Similarly, element S4 assumes that the transverse shear strain (and force, since the transverse shear treatment is elastic based on the initial elastic modulus of the material) is constant over the element. Therefore, all four stiffness integration locations will have the same transverse shear strain, transverse shear section force, and transverse shear stress distribution. The transverse shear treatment for S4 is identical to that for S4R.

It is well known that a standard displacement formulation will exhibit shear locking for applications dominated by in-plane bending deformation. However, a standard displacement formulation for the out-of-plane bending stiffness is not subject to similar locking response. Hence, S4 uses a standard displacement formulation for the element's bending stiffness, and the theory presented above for the rotation kinematics and bending strain measures applies to S4. The primary difference between the element formulations for S4 and S4R is the treatment of the membrane strain field. This formulation is the topic of the following discussion.

The membrane formulation used for S4 does not rely on the fact that S4 is a shell element. Hence, the discussion below details the formulation from the point of view that the membrane response is governed by the equilibrium for a three-dimensional body in a state of plane stress.

Consider an enhancement $\tilde{D}$ to the rate of deformation tensor $D$ . We introduce the enhanced rate of deformation tensor, $\bar{D}$ , as

\bar{D} \overset{def}{=} D + \tilde{D},

where $\tilde{D}$ is defined subsequently.

Admissible variations in the rate of deformation are also introduced as

\bar{δ} D \overset{def}{=} δ D + δ \tilde{D},

where

δ D \overset{def}{=} sym (\frac{\partial δ x}{\partial x}) .

We now introduce constraints on the enhancements $D$ and $\tilde{D}$ :

\int_{V} (\bar{D} - D) : \overset{*}{σ} d V = \int_{V} (δ \bar{D} - δ D) : \overset{*}{σ} d V = 0

so that the modified virtual work statement can be written in the form

\int_{V} σ (\bar{D}) : δ \bar{D} d V + \int_{V} b \cdot δ x d V - \int_{S_{T}} t \cdot δ x d S_{T} = 0,

where $t$ is the specified traction on $S_{T}$ and $δ x = 0$ on $S_{u}$ . $\overset{*}{σ}$ is an arbitrary stress field, and the constitutive equation $σ = σ$ ( $\bar{D}$ ) is enforced pointwise.

In the modified virtual work statement all kinematic quantities and corresponding variations ( $\bar{D}$ , $D$ , $δ \bar{D}$ , and $δ D$ ) are known functions of $x$ , $δ x$ , and the reference configuration. A fundamental requirement for the validity of the formulation is that the modified virtual work statement leads to the proper equilibrium equations. If $\overset{*}{σ}$ is arbitrary, the constraint equations can be rewritten as

\bar{D} = D and δ \bar{D} = δ D .

Substituting these two relations in the modified virtual work equation yields

\int_{V} σ : δ D d V + \int_{V} b \cdot δ x d V - \int_{S_{T}} t \cdot δ x d S_{T} = 0

where we have used the constitutive equation $σ = σ$ $(D)$ . We recognize this variational statement as the usual virtual work equation, and a straightforward application of the divergence theorem leads to the standard equilibrium equations.

In the actual implementation we choose to satisfy the constraints only for piecewise constant stress fields $\overset{*}{σ}$ . Hence, over the element domain $V_{e}$ we require

\int_{V^{e}} \tilde{D} d V = \int_{V^{e}} δ \tilde{D} d V = 0.

The enhancements $\tilde{D}$ and $δ \tilde{D}$ are chosen such that they eliminate the shear locking for in-plane bending. In addition, the direct strain field is enhanced to approximate the strains caused by Poisson's effect in bending.

Patch test

To pass the patch test, the choice of enhancements $\tilde{D}$ cannot be arbitrary. A sufficient condition for the satisfaction of the patch test is that for homogeneous deformations we have $\bar{D} = D$ or $\tilde{D} = 0$ pointwise. In that case $σ (\bar{D}) = σ (D)$ and the stress is homogeneous. Since the stress is homogeneous, it can be moved outside the volume integral in the modified virtual work statement. The volume integral condition on $δ \tilde{D}$ implies that the expression is independent of the enhancement and leads to the standard displacement formulation, which is known to satisfy the patch test.