15. Element Formulations

This chapter covers the various elements that the Solid Mechanics module has available. We first discuss all of the solid elements, then the shell elements, and finally the beam elements.

15.1. Uniform Gradient Hex8 Solid Element

This element has proven to be the workhorse solid element for Solid Mechanics. The hex8 is an eight-node hexahedron element with a topology and a node numbering convention shown in Fig. 15.1. It also is referred to as the mean-quadrature Hex8 because of the particular manner in which it generates a mean-quadrature representation of the gradient (and divergence) operator. The approach adapted for developing a mean strain rate quadrature for the eight-node hexahedron is that given by Reference [[1]]. While the approach and notation is cumbersome, it provides the structure needed to achieve a closed-form solution for the integration of an arbitrary hexahedron and an explicit and unambiguous identification of the orthogonal hourglass modes.

Fig. 15.1 Isoparametric coordinate representation of the eight-noded hex element.

15.1.1. Kinematics

The eight-node solid hexahedron element relates the spatial coordinates \(x^i\) to the nodal coordinates \(x^i_I\) through the isoparametric shape functions \(N_I\) as follows:

(15.1)\[ x^i = x^i_I N^I (\xi^i).\]

In accordance with index notation convention, repeated subscripts imply summation over the range of that subscript. The lower case subscripts have a range of three, representing the spatial coordinate directions. Upper case subscripts have a range of eight, corresponding to element nodes.

The same shape functions are used to define the element displacement field in terms of the nodal displacements \(u_{iI}\):

(15.2)\[ u_i = u_{iI} N^I (\xi^i).\]

Since these shape functions apply to both spatial coordinates and displacement, their material derivative (represented by a superposed dot) must vanish. Hence, the velocity field is given by:

(15.3)\[ v_i = v_{iI} N^I (\xi^i).\]

The velocity gradient \(v_{i,j}\) is defined as follows:

(15.4)\[ v_{i,j} = v_{iI} N^I_{,j} .\]

By convention, a comma preceding a lower case subscript denotes differentiation with respect to the spatial coordinates, hence \(v_{i,j}\) denotes \(\partial v_i / \partial x^j\).

The shape functions \(N^I\) map a unit cube in the isoparametric coordinates \(\xi^i\) to a general hexahedron in the spatial coordinates \(x^i\). The unit cube is centered at the origin in \(\xi^i\)-space so that the shape functions may be conveniently expanded in terms of an orthogonal set of base vectors, given in Table 15.1, as follows:

(15.5)\[ N^I (\xi^i) = \frac{1}{8} \Sigma^I + \frac{1}{4} \xi^i \Lambda^I_i + \frac{1}{2} \xi^2 \xi^3 \Gamma^I_1 + \frac{1}{2}\xi^1 \xi^3 \Gamma^I_2 + \frac{1}{2} \xi^1 \xi^2 \Gamma^I_3 + \xi^1 \xi^2 \xi^3 \Gamma^I_4.\]

Table 15.1 Deformation modes of the eight-noded hex element

	\(\xi^1\)	\(\xi^2\)	\(\xi^3\)	\(\Sigma^I\)	\(\Lambda^I_1\)	\(\Lambda^I_2\)	\(\Lambda^I_3\)	\(\Gamma^I_1\)	\(\Gamma^I_2\)	\(\Gamma^I_3\)	\(\Gamma^I_4\)
1	\(-\frac{1}{2}\)	\(-\frac{1}{2}\)	\(-\frac{1}{2}\)	1	-1	-1	-1	1	1	1	-1
2	\(\frac{1}{2}\)	\(-\frac{1}{2}\)	\(-\frac{1}{2}\)	1	1	-1	-1	1	-1	-1	1
3	\(\frac{1}{2}\)	\(\frac{1}{2}\)	\(-\frac{1}{2}\)	1	1	1	-1	-1	-1	1	-1
4	\(-\frac{1}{2}\)	\(\frac{1}{2}\)	\(-\frac{1}{2}\)	1	-1	1	-1	-1	1	-1	1
5	\(-\frac{1}{2}\)	\(-\frac{1}{2}\)	\(\frac{1}{2}\)	1	-1	-1	1	-1	-1	1	1
6	\(\frac{1}{2}\)	\(-\frac{1}{2}\)	\(\frac{1}{2}\)	1	1	-1	1	-1	1	-1	-1
7	\(\frac{1}{2}\)	\(\frac{1}{2}\)	\(\frac{1}{2}\)	1	1	1	1	1	1	1	1
8	\(-\frac{1}{2}\)	\(\frac{1}{2}\)	\(\frac{1}{2}\)	1	-1	1	1	1	-1	-1	-1

The above vectors represent the deformation modes of a unit cube, as shown in Fig. 15.2. The first vector, \(\Sigma^I\) accounts for rigid body translation. The linear base vectors \(\Lambda^I_i\) may be readily combined to define three uniform normal strain rate modes, three uniform shear strain rate modes, and three rigid body rotation rates for the unit cube. The last four vectors \(\Gamma^I_{\alpha}\) (Greek subscripts have a range of four) give rise to modes with linear strain variations which are neglected by mean strain quadrature. These vectors define the hourglass patterns for a unit cube. Hence the modes \(\Gamma^I_{\alpha}\) are referred to as the hourglass base vectors.

Fig. 15.2 Deformation modes of the eight-noded hex element.

15.1.2. Mean Quadrature

The variational statement gives the following relationship for the element nodal forces \(f^{iI}\) due to the divergence of the stress field,

(15.6)\[ v_{iI} f^{iI} = \int_V t^{ij} d_{ij} \mathrm{d}v.\]

The integral in (15.6) is evaluated using a constant stress, thereby considering only a mean strain rate within the element:

(15.7)\[ v_{iI} f^{iI} = \int_V t^{ij} d_{ij} \mathrm{d}v = V \bar{t}^{\,ij} \bar{v}_{i,j}.\]

The assumed constant stress field is represented by \(\bar{t}^{\,ij}\), which will be referred to as the mean stress tensor. It is assumed that the mean stress depends only on the mean strain. Mean kinematic quantities are defined by integrating over the element as follows:

(15.8)\[ \bar{v}_{i,j} = \frac{\displaystyle 1}{ \displaystyle V} \int_V v_{i,j} \mathrm{d}v.\]

The gradient operator \(B^I_i\) is defined by

(15.9)\[ B^I_i = \int_V N^I_{,j} \mathrm{d}v.\]

The mean velocity gradient, applying (15.9) is then given by

(15.10)\[ \bar{v}_{i,j} = \frac{1}{V} v_{iI} B^I_j.\]

The nodal forces are then given by

(15.11)\[ f^{iI} = \bar{t}^{\,ij} B^I_j.\]

Computing nodal forces by this integration scheme requires evaluation of the gradient operator \(B^I_j\) and volume. These two tasks can be linked together by using the relationship \(x^i_{,j} = \delta^i_j\). Therefore (15.9) yields

(15.12)\[ x^i_I B^I_j = \int_V \left( x^i_I N^I \right)_{,j} \mathrm{d}v = V\delta^i_j.\]

Consequently, the gradient operator \(B^I_i\) may alternatively be expressed by

(15.13)\[ B^I_j = \frac{\partial V}{\partial x^i_I}.\]

To integrate the element volume in closed form, the Jacobian of the isoparametric transformation is used transform the integral over the unit cube,

(15.14)\[ V = \int\limits_V \mathrm{d}V = \int\limits_{-1/2}^{1/2} \int\limits_{-1/2}^{1/2} \int\limits_{-1/2}^{1/2} J \mathrm{d}\xi^1 \mathrm{d}\xi^2 \mathrm{d}\xi^3.\]

The Jacobian \(J\) is the determinant of the transformation operator \(\partial x^i / \partial \xi^j\) and may be expressed as

(15.15)\[ J = e^{ijk} \frac{\partial x^1}{\partial \xi^i} \frac{\partial x^2}{\partial \xi^j} \frac{\partial x^3}{\partial \xi^k}.\]

Using (15.1), (15.14), and (15.15), the element volume may be expressed in the following form:

(15.16)\[ V = x^1_I x^2_J x^3_K D^{IJK},\]

where

(15.17)\[ D^{IJK} = e^{ijk} \int\limits_{-1/2}^{1/2} \int\limits_{-1/2}^{1/2} \int\limits_{-1/2}^{1/2} J \frac{\partial N^I}{\partial \xi^i} \frac{\partial N^J}{\partial \xi^j} \frac{\partial N^K}{\partial \xi^k} \mathrm{d}\xi^1 \mathrm{d}\xi^2 \mathrm{d}\xi^3.\]

Observe that the coefficient array \(D^{IJK}\) is identical for all hexahedra. Furthermore, it possesses the alternator properties given by

(15.18)\[ D^{IJK} = D^{JKI} = D^{KIJ} = -D^{IKJ} = -D^{JIK} = -D^{KJI}.\]

Therefore, applying (15.13) and (15.14) to the expression (15.16) yields the following closed-form expression for evaluation the components of the gradient operator, \(B^I_i\):

(15.19)\[\begin{split}\begin{bmatrix} B^I_1 \\ B^I_2 \\ B^I_3 \end{bmatrix} = \begin{bmatrix} x^2_J x^3_K \\ x^3_J x^1_K \\ x^1_J x^2_K \end{bmatrix} D^{IJK} .\end{split}\]

Looking at the form of (15.19), it is evident that evaluating each component of \(D^{IJK}\) involves integrating a polynomial which is at most bi-quadratic. However, since the integration is over a symmetric region, any term with a linear dependence will vanish. The only terms which survive the integration will be the constant, square, double square, and triple square terms. Furthermore, the alternator properties cause half of these remaining terms to drop out. The resulting expression for \(D^{IJK}\) is

(15.20)\[ D^{IJK} = \frac{\displaystyle 1}{\displaystyle 192} e^{ijk} \left( 3\Lambda^I_i \Lambda^J_j \Lambda^K_k + \Lambda^I_i \Gamma^J_k \Gamma^K_j + \Gamma^I_k \Lambda^J_j \Gamma^K_i + \Gamma^I_j \Gamma^J_i \Lambda^K_k \right).\]

The expression in (15.20) is evaluated using Table 15.1, after which practical formula for computing the gradient operator \(B^I_i\) and volume are developed.

The gradient operator component \(B^1_x\) is given explicitly by

(15.21)\[ B^1_x = \left[ y_2 ( z_{63} - z_{45} ) + y_3 z_{24} + y_4 ( z_{38} - z_{52} ) + y_5 ( z_{86} - z_{24} ) + y_6 z_{52} + y_8 z_{45} \right] / 12,\]

where \(\{ x^i_I \} = \{ x_I , y_I , z_I \}\) and \(z_{IJ} = z_I - z_J\). To obtain the balance of the gradient operator components \(B^I_x\), the nodal index permutations contained in Table 15.2 are used. To obtain the components \(B^I_y\) and \(B^I_z\), the coordinate permutations contained in Table 15.3 are used.

Table 15.2 Permutation of Nodal Indices

\(B^1_i\)	2	3	4	5	6	7	8
\(B^2_i\)	3	4	1	6	7	8	5
\(B^3_i\)	4	1	2	7	8	5	6
\(B^4_i\)	1	2	3	8	5	6	7
\(B^5_i\)	8	7	6	1	4	3	2
\(B^6_i\)	5	8	7	2	1	4	3
\(B^7_i\)	6	5	8	3	2	1	4
\(B^8_i\)	7	6	5	4	3	2	1

Table 15.3 Permutation of Nodal Coordinates

\(B^I_x\)	y	z
\(B^I_y\)	z	x
\(B^I_z\)	z	y

It is worth noting at this point the difference between the mean quadrature (alt. mean strain rate, mean stress) approach and one-point Gauss quadrature. The latter method would effectively neglect the last three terms of (15.21). In a parallelepiped, the nodal coordinates contain no component of the hourglass base vectors, consequently, only the first term of (15.21) is necessary to compute the gradient operator and volume. In such a case, one-point quadrature is equivalent to the mean quadrature formula. However, for a general hexahedron shape, one-point quadrature does not correctly assess a state of uniform stress and strain, thus, may not be convergent [Zienkiewicz, 1977]. In view of the requirements of the Iron’s patch test, it is likely that (15.20) is unique.

15.1.3. Orthogonal Hourglass Control

The mean stress - mean strain rate formulation considers only the linear part of the velocity field. The remaining portion of the velocity field is the so-called hourglass field. Excitation of these modes may lead to sever, unresisted mesh distortion. A method for isolating the hourglass modes so that they may be treated independently of the rigid body and uniform strain modes is required. This is accomplished by developing an hourglass gradient operator connected with hourglass restoring forces. The linear velocity field on which the mean strain rates are based is given by

(15.22)\[ v^{\text{LIN}}_i = v_{iI} \left( \frac{1}{8} \Sigma^I + \frac{1}{V} ( x^j - \frac{1}{8} x^j_J \Sigma^J ) B^I_j \right) .\]

The hourglass velocity field \(v^{\text{HG}}_i\) may be defined by removing the linear portion of the velocity field. Thus,

(15.23)\[ v^{\text{HG}}_i = v_i - v^{\text{LIN}}_i ,\]

or in terms of the nodal velocities,

(15.24)\[ v^{\text{HG}}_i = v_{iI} - v_{i0} \Sigma_I - \frac{1}{V} \left( x^j_I - x^j_0 \Sigma_I \right) v_{iJ} B^J_j ,\]

where \(v_{i0} = \frac{1}{8} v_{iI} \Sigma_I\) and \(x^i_0 = \frac{1}{8} x^i_I \Sigma_I\).

The hourglass velocity field, (15.24), is in the improper null space of the gradient operator \(B^I_i\). The linear velocity field, (15.22), spans 12 degrees of freedom: 3 rates of rigid body translation, 3 rates of rigid body rotation, and 6 uniform strain rates, which means that the hourglass subspace is remaining 12 degrees of freedom.

An hourglass gradient operator is constructed from the hourglass basis vectors \(\Lambda^I_{\alpha}\) as follows:

(15.25)\[G^I_{\alpha} =\frac{V}{\delta} \bigg[ \Lambda^I_1 , ~ \Lambda^I_2 , ~ \Lambda^I_3 , ~ \Lambda^I_4 \bigg] ,\]

where \(\delta\) is a generalized element dimension developed below. This scaling provides the hourglass gradient operator with the same dimensional characteristics as the uniform gradient operator. While \(G^I_{\alpha}\) is orthogonal to \(B^I_i\), the following property:

(15.26)\[B^I_i G^I_{\alpha} \neq 0 ,\]

means that \(G^I_{\alpha}\) used with the full velocity field \(v_{iI}\) will couple the hourglass behavior to the uniform strain rate behavior. Thus, hourglass strain rates \(\dot{q}_{i \alpha}\) are developed with \(G^I_{\alpha}\) operating on only the hourglass velocities \(v^{\text{HG}}_{iI}\),

(15.27)\[\dot{q}_{i \alpha} = \frac{1}{V} v^{\text{HG}}_{iI} G^I_{\alpha} .\]

Alternatively, an unrestricted operator may be developed by requiring is to satisfy the following condition:

(15.28)\[v_{iI} \gamma^I_{\alpha} = v^{\text{HG}}_{iI} G^I_{\alpha} .\]

Using the hourglass velocity, (15.24), provides

(15.29)\[v_{iI} \gamma^I_{\alpha} = \Big[ v_{iI} - v_{i0} \Sigma^I - \frac{1}{V} \left( x^j_I - x^j_0 \Sigma_I \right) v_{iJ} B^J_j \Big] G^I_{\alpha} ,\]

which, when rearranged and using the orthogonality of the mode shapes \(\Sigma^I\) and \(\Gamma^I_{\alpha}\) (i.e., \(\Sigma^I \Gamma^I_{\alpha}=0\)) gives

(15.30)\[v_{iI} \gamma^I_{\alpha} = v_{iI} \left( G^I_{\alpha} - \frac{1}{V} x^j_J G^J_{\alpha} B^I_j \right).\]

The condition for the unrestricted operator is satisfied if the hourglass operator \(\gamma^I_{\alpha}\) is defined as

(15.31)\[\gamma^I_{\alpha} = \frac{V}{\delta} \left( \Gamma^I_{\alpha} - \frac{1}{V} x^j_J \Gamma^J_{\alpha} B^I_j \right) ,\]

and the hourglass strain rates are defined as

(15.32)\[\dot{q}_{i \alpha} = \frac{1}{V} v_{iI} \gamma^I_{\alpha} .\]

To control the hourglass modes, generalized forces \(Q^{i \alpha}\) are defined which are conjugate to \(\dot{q}_{i \alpha}\), so that the work rate is given by

(15.33)\[v_{iI} f^{iI}_{\text{HG}} = V Q^{i \alpha} \dot{q}_{i \alpha} .\]

Utilizing (15.31), the contribution to the nodal forces due to hourglass resistance is given be

(15.34)\[f^{iI}_{\text{HG}} = Q^{i \alpha} \gamma^I_{\alpha}.\]

The hourglass restoring forces are calculated from

(15.35)\[\check{Q}^{i \alpha} = \epsilon ~2 \mu_{tan} \delta^{ij} \delta^{\alpha \beta} \dot{q}_{j \beta} ,\]

where \(2 \mu_{tan}\) is the tangent shear stiffness obtained from the deviatoric constitutive behavior of the mean stress and mean strain rate in the element, and \(\epsilon\) is a scaling parameter. The scaling \(\epsilon\) assures the level of the hourglass restoring forces remains below that of the mean deviatoric stress state.

The deviatoric behavior is used since the hourglass modes are constant volume, higher order straining modes of the element. The tangent modulus assures that the evolution of the hourglass restoring forces parallels that of the mean deviatoric stress state.

The invariant time derivative of the generalized forces \(Q^{i \alpha}\) accounts for the finite rotations expected in use of the element in applications. The derivative is given by

(15.36)\[\check{Q}^{i \alpha} = \dot{Q}^{i \alpha} - \omega_{ij} Q^{j \alpha} ,\]

where \(\omega_{ij}\) is the spin tensor.

The hourglass restoring forces are added to those obtained from the divergence of the mean stress state so that the complete result is

(15.37)\[f^{iI} = \left( \bar{t}^{\,ij} B^I_j + Q^{i \alpha} \gamma^I_{\alpha} \right) .\]

15.1.4. Linear Hyperelastic Hourglass Control

The traditional hourglass formulation is in rate form, and as such is not guaranteed to be energetically reversible and may not result in symmetric contributions to the finite element stiffness matrix. Here we describe a hyperelastic hourglass formulation that overcomes these limitations and also allows us to define Lagrangian hourglass strains which are valuable for extending the formulation to nonlinear hourglass response, as discussed in the next section.

Consider an element’s total hourglass energy defined as follows:

(15.38)\[\psi_e(\mathbf{x}) = V \epsilon \mu_{tan} \left( \sum_{\alpha=1}^4 \varepsilon_\alpha^2 \right),\]

where \(V\) is the reference volume of the element, \(\mathbf{x}\) is the element’s current nodal coordinates, and \(\varepsilon_\alpha\) is a measure of the hourglass strain for hourglass mode \(\alpha\). The hourglass strain here is based on an hourglass operator defined in the model/reference coordinates (this is in contrast to the total and incremental hourglass formulations which use hourglass operators in the current configuration). The hourglass strains are given by

\[\varepsilon_\alpha = \sqrt{ \sum_{i=1}^3 {\epsilon}_{i \alpha} {\epsilon}_{i \alpha} } \quad (\textrm{no sum over }\alpha),\]

with hourglass strain vector \(\epsilon\):

\[\epsilon_{i \alpha} = \sum_{I=1}^8 {H}_{\alpha}^I {x}_{i}^I,\]

where \(H_\alpha^I\) is the unrestricted hourglass operator defined in the reference configuration:

\[{H}_{\alpha}^I = \frac{1}{\delta} \left[ {\Gamma}_{\alpha}^I - \frac{1}{V} {X}_{j}^J {\Gamma}_{\alpha}^J {B}_{j}^I \right],\]

where \({B}\) is the gradient operator with respect to the reference configuration, and \(\mathbf{X}\) is the element’s reference coordinates. Note that this definition of the unrestricted hourglass is analogous to \(\gamma^I_{\alpha}\) from Section 15.1.3, but includes a factor of \(\frac{1}{V}\) to simplify notation. The value of \(\delta\), a characteristic element length scale, is chosen so that the hyperelastic hourglass forces match the classical formulation at small stains. The proposed hourglass energy is objective due to the definition of the hourglass strain, which is invariant to rigid body rotations of the current coordinates, \(\mathbf{x}\). The resulting hourglass forces follow from work conjugacy:

\[{f}_{\text{HG}}^{iI} = -\frac{\partial\psi_e}{\partial{x}_i^I}.\]

Being derived from an objective potential energy, these forces are objective, path independent. The resulting stiffness matrix is guaranteed to be symmetric and hourglass deformations are elastically reversible. Because we are assuming an energy which is quadratic in the hourglass strain, the resulting forces are linear with hourglass displacement. In addition, the unrestricted hourglass operator \(H\) has both rigid body modes and affine deformation modes in its null-space, meaning the hourglass forces will be orthogonal to affine motions of the element.

An extension to a non-linear hyperelastic hourglass formulation is described below.

15.1.5. Nonlinear Hyperelastic Hourglass Control

An additional limitation of the traditional hourglass control is that the hourglass resistance is typically formulated to be (incrementally) linearly proportional to hourglass deformation increments. A hyperelastic hourglass control formulation can overcome this limitation by using a nonlinear hyperelastic energy which is a function of the four hourglass strains. We use a generalized definition of the energy function from (15.38):

\[\hat{\psi}_e(\mathbf{x}) = V \epsilon \mu_{tan} \varepsilon_0^2\left( \sum_{\alpha=1}^4 \left[ e_m\left(\frac{\varepsilon_\alpha}{\varepsilon_0}\right)\right]^2 \right),\]

where \(\varepsilon_0\) is called the transition strain, and \(e_m(\cdot)\) is a function which takes a strain and returns an alternative Seth-Hill strain measure. In particular, it is defined by

\[e_m(\varepsilon) = \frac{1}{m} \left( \left(\varepsilon+1 \right)^m -1\right),\]

if \(\varepsilon\) is a strain, \(e_m(\varepsilon)\) satisfies all the requirements of a strain measure. In particular, it satisfies \(e_m(0)\) = 0, \(e_m^\prime(0) = 1\). As a result, at small strains, \(e_m \approx \varepsilon\) for any \(m\). A value of \(m=2\) corresponds to a Green-Lagrange strain measure, while \(m=1\) is an identity map. The variable \(\varepsilon_0\) is called the transition strain because it sets the strain level at which the nonlinearity of the strain measure begins to become dominant. For a very large transition strain, the hourglass force response will remain linear up to large hourglass strains as the term \(\frac{\varepsilon_\alpha}{\varepsilon_0}\) remains small. For small transition strains, the nonlinearity become noticeable earlier. The hourglass energy for this model scales as \(\varepsilon^{2m}\) as strains get large, meaning the hourglass force scales as the hourglass displacement to the \(2m-1\) power at large deformations.

Fig. 15.3 shows the hourglass resistance force vs displacement for varying transition strains and Seth-Hill exponent \(m\).

Fig. 15.3 Nonlinear hourglass force versus displacement.

15.2. Tet4 Solid Element

This element is the standard 4-noded tetrahedral element. It is notoriously stiff and prone to locking, but included for completeness. More information on this element can be found in [[2]].

15.3. Tet10 Solid Element

The default 10-noded tetrahedral solid element is the composite tetrahedron, as given in [[3]], which is an extension of the composite tetrahedron and triangle formulations in [[4]] and [[5]]. This 10-noded tetrahedron consists of 12 linear, 4-noded sub-tetrahedra. The nodal fields, including displacement, are linear within each sub-tetrahedron and, therefore, piecewise linear within the parent 10-noded tetrahedron. The deformation gradient and stress fields are formulated to be linear over this 10-noded tetrahedron; the gradient operator projects the piecewise discontinuous gradients among the 12 sub-tetrahedra into a linear basis on the parent tetrahedron.

As stated in [[3]], [[4]], and [[5]], there are several advantages of this composite tetrahedron formulation over commonly used alternatives. Tetrahedral elements provide generally more robust and efficient finite element meshing than hexahedral elements. As opposed to the traditional formulation of the quadratic 10-noded tetrahedron, this formulation has a well-balanced mass lumping to all 10 nodes lending to improved performance in explicit transient dynamics, contact, and other solid mechanics capabilities that rely upon the nodal mass distribution. Volumetric locking and unrealistic pressure oscillations are still possible for this element when modeling isochoric deformation (plasticity) and nearly incompressible materials, but this behavior can be alleviated further by volume-averaging the dilation over the element [[3]]. In Sierra/SM, this option is called VOLUME AVERAGE J = ON, which is a default setting for the composite tetrahedron.

A 10-noded, quadratic, fully-integrated tetrahedral element is also available in Sierra/SM. For more information on this element, refer to [[2]].

15.4. Belytschko-Tsay Shell Element

The 4-noded Belytschko-Tsay shell (or BT-shell4) is the simplest of the shell elements offered. The original reference can be found in [[6]]. It should be considered as the minimal 5-parameter Mindlin-type formulation that includes a constant transverse shear contribution.

15.5. Key-Hoff Shell Element

The 4-noded Key-Hoff shell (or KH-shell4) is slightly more involved than the BT-shell4, in that it includes a term for a linear-varying transverse shear in its formulation. The inclusion of this term is an improvement on the BT-shell4 because it properly models warped shell geometry - albeit in a low-order way. The original reference for this element can be found in [[7]].

15.6. Belytschko-Leviathan Shell Element

The 4-noded Belytschko-Leviathan shell (or BL-shell4) is slightly more involved than the KH-shell4, in that it includes additional shear terms as well as additional hourglass controls. The inclusion of the hourglass terms (also known as the physical stabilization parameter) is an improvement on the KH-shell in that it eliminates some of the over-stiffness sometimes observed in the KH-shell4. The BL-shell4 also includes a projection of the angular velocities and the internal forces. The original references for this element can be found in [[8]] and [[9]].

15.7. Shear Correction for Layered Shell Elements

For sandwich composite plates with a low modulus core, the effects of transverse shear deformation can be significant. Thus, the results of first-order shear deformation theory, as applied for layered shell elements in Sierra/SM, are affected by the choice of shear correction factor (\(K\)). In reference [[10]], an expression is derived for the variation of transverse shear through the thickness of a laminated plate. The expression given for the shear correction factor is:

(15.39)\[K = \left[ \left( A_{44}-\frac{A^{2}_{45}}{A_{55}} \right) \int_{-h/2}^{h/2} \frac{\left[ \int_{-h/2}^{z} (\bar{Q}_{1i}\beta_{1i}+z\bar{Q}_{1i}\delta_{1i})dz \right]^2}{\left[\bar{q}_{44}-\frac{\bar{Q}_{45}^2}{\bar{Q}_{55}}\right]} \right]^{-1}, \mbox{ for } i = 1,2,6.\]

Here, the \(A_{ij}\) are the corresponding terms in the laminate extensional stiffness matrix, the \(\bar{Q}_{ij}\) are the terms of the reduced stiffness matrix, \(\beta_{ij}\) and \(\delta_{ij}\) are terms of the compliance sub-matrices, and \(h\) is the thickness of the section. When expressed in algebraic form, for a laminate of \(N\) layers and a coordinate system centered at the centroidal axis of the section, the shear correction factor can be expressed as

(15.40)\[K = \frac{\left[ A_{44}-\frac{A^{2}_{45}}{A_{55}} \right]^{-1}}{\sum^{N}_{k=1} \frac{1}{\left( \bar{Q}^{k}_{44}-\frac{\bar{Q}^{k^2}_{45}}{\bar{Q}^k_{55}} \right)} \left[ P_k(z_{k+1}-z_k) + \frac{R_k}{2}(z^{2}_{k+1}-z^{2}_{k}) + \frac{V_k}{3}(z^{3}_{k+1}-z^{3}_{k}) + \frac{W_k}{4}(z^{4}_{k+1}-z^{4}_{k}) + \frac{X_k}{5}(z^{5}_{k+1}-z^{5}_{k}) \right]},\]

where

(15.41)\[\begin{split}P_{k} = T^{2}_{k} + H^{2}_{k}z^{2}_{k} - 2T_kH_kz_k + U^{2}_{k} + \frac{J^{2}_{k}z^{4}_{k}}{4} - U_kJ_kz^{2}_{k} + 2T_kU_k - T_kJ_kz^{2}_{k} - 2H_kU_kz_k + H_kJ_kz^{3}_{k}, \\ R_{k} = 2T_kH_k - 2H^{2}_{k}z_k + 2H_kU_k - H_kJ_kz^{2}_{k}, \\ V_k = H^{2}_{k} - \frac{J^{2}_{k}z^{2}_{k}}{2} + U_kJ_k + T_kJ_k - z_kH_kJ_k, \\ W_k = H_kJ_k, \\ X_k = \frac{J^{2}_{k}}{4}, \\ T_k = \sum^{k-1}_{m=1} H_m(z_{m+1}-z_m), \\ U_k = \sum^{k-1}_{m=1} \frac{J_m}{2}(z^{2}_{m+1}-z^{2}_m), \\ H_k = \bar{Q}^{k}_{1i}\beta_{1i} \mbox{ , for } i=1,2,6,\end{split}\]

and

\[J_k = \bar{Q}^{k}_{1i}\delta_{1i} \mbox{ , for } i=1,2,6.\]

For the laminate cross-section geometry, see Figure \(2\) in reference [[10]]. This correction factor has been coded into a subroutine and is used for the layered shell formulation in Sierra/SM.

15.8. 3D Beam Element

The two-noded beam in Sierra/SM is based on conventional Timoshenko beam theory in which the functional form of the deformation is made explicit on a cross section normal to the reference axis. Thus, the deformation is described in terms of kinematic variables that depend on the coordinate along the reference axis. As shown in Fig. 15.4, the axis connecting node 1 and node 2 labeled \(\xi_1\) is this reference axis. The beam is defined by a cross-section of fixed-shape existing uniformly along the reference axis and is formulated using isoparametric coordinates.

As will be apparent, the assumptions about the deformation of the beam are those of a Timoshenko beam theory. In particular, the transverse shear deformation is modeled. Planar cross-sections originally perpendicular to the reference axis remain flat and undeformed though not necessarily remain perpendicular to the reference axis under deformation.

When initially curved beams are modeled with straight beam segments, the global curvature properties are represented by the change in orientation from one beam element to the next. In effect, the smooth variation in curvature of the original reference axis is approximated by discrete changes in orientation occurring at the element ends; the elements are chord approximations to the original curved beam, much like linear shell elements when modeling a curved structure. This approximation is the same order as the constant membrane and bending stress approximations introduced in the element integration.

Fig. 15.4 Isoparametric coordinate representation of the two-noded beam element.

15.8.1. Kinematics

The motion of the beam is defined in terms of the velocity of the reference axis and the additional rotation of the region within the cross-section defined by \(A(\xi_2,\xi_3)\),

(15.42)\[v_i (x^j) = v_i (\xi_1) - \varepsilon_{imn} \rho^m \omega^n (\xi_1) .\]

Here, \(\rho^m\) is the position vector from the reference axis to a point in the cross-section \(A(\xi_2,\xi_3)\). The position vector is perpendicular to the reference axis and has the units of length.

Based on (15.42), the spatial gradient of the velocity is given by

(15.43)\[v_{i,j} (x^k) = v_{i,j} (\xi_1) - \varepsilon_{imn} \rho^m \omega^n_{,j} (\xi_1).\]

[Note: In the special case when the isoparametric coordinates \(\xi_i\) coincide with the spatial coordinates \(x_i\), the velocity of the beam is given by:

(15.44)\[\{v_i \} = \{ (v_x + z \omega_y - y \omega_z), \ (v_y - z \omega_x), \ (v_z + y \omega_x)\}.\]

The stretching (symmetric part of the velocity gradient) is then given by

(15.45)\[\begin{split}\begin{aligned} d_{xx} &= v_{x,x} + z\omega_{y,z} - y \omega{z,x}, \\ d_{yy} &= 0, \\ d_{zz} &= 0, \\ 2d_{xy} &= -\omega_z + v_{y,x} - z\omega{x,x}, \\ 2d_{xz} &= \omega_y + v_{z,x} - y\omega{x,x}, \\ 2d_{yz} &= 0, \end{aligned}\end{split}\]

and the spin (skew-symmetric part of the velocity gradient) is given by

(15.46)\[\begin{split}\begin{aligned} 2\omega_{xy} &= -\omega_z - v_{y,x} + z\omega{x,x}, \\ 2\omega_{xz} &= \omega_y - v_{z,x} - y\omega{x,x}, \\ 2\omega_{yz} &= -2\omega_x, \end{aligned}\end{split}\]

where it is now apparent that Timoshenko beam theory allowing transverse shear deformation is considered, see Reference [[11]]. Using Timoshenko beam theory allows the rotation rates \(\omega_y\) and \(\omega_z\) to be described separately, rather than defined by \(-v_{z,x}\) and \(v_{y,x}\), respectively. Consequently, the (separate) finite element assumptions on the velocity and rotation rates are required to be no more than continuous represented. In the event that the slender beam limit of vanishing transverse shear strains holds, classical beam theory is recovered, though special considerations in the element formulation (introduced below) are needed to prevent shear locking.]

Returning to our description of the more general case, the two-noded beam relates the spatial coordinates \(x_{iI}\) through the isoparametric shape functions \(N_I, \ I=1,2\) as follows:

(15.47)\[x_i = x_{iI} N_I (\xi_1).\]

The shape functions map a unit interval in the isoparametric coordinate \(\xi_i\) to a general beam segment in the spatial coordinates, \(x_i\). The unit interval is centered at the origin in the \(\xi_1\)-space so that the shape functions may be conveniently expanded in terms of an orthogonal set of base vectors:

(15.48)\[N_I (\xi_1) = \frac{1}{2} \Sigma_I + \xi_1 \Lambda_I\]

where at node 1: \(\xi_1 = -\frac{1}{2}\), \(\Sigma_1 = 1\), \(\Lambda_1=-1\), and at node 2: \(\xi_1 = +\frac{1}{2}\), \(\Sigma_2 = 1\), \(\Lambda_2=1\). As shown in Fig. 15.5, these two modes represent the deformation modes of a unit interval \(-\frac{1}{2} \leq \xi_1 \leq \frac{1}{2}\).

Fig. 15.5 Deformation modes of a unit interval.

Although the velocity gradient of the two-noded beam is quite complex in description when using a Timoshenko beam theory (15.43), the modes \(\Sigma_I\) and \(\Lambda_I\) combine to represent rates of rigid body translation and rotation, and the uniform strain rates, with no hourglass mode of deformation.

The same shape functions are used to define the reference axis displacement in terms of the nodal displacements, \(u_{iI}\):

(15.49)\[u_i = u_{iI} N_I (\xi_1) .\]

Since these shape functions apply to spatial coordinates and displacements, their material derivatives must vanish. Hence, the velocity field and rotational rate are given by

(15.50)\[\begin{split}\begin{aligned} v_i &= v_{iI} N_I (\xi_1), \\ \omega_i &= \omega_{iI} N_I. (\xi_1) \end{aligned}\end{split}\]

The velocity gradient and the gradient of the rotational rate are defined as follows:

(15.51)\[\begin{split}\begin{aligned} v_{i,j} &= v_{iI} N_{I,j}, \\ \omega_{i,j} &= \omega_{iI} N_{I,j}. \end{aligned}\end{split}\]

15.8.2. Mean Quadrature

In order to introduce the concept of a mean (constant) strain and stress in the beam, we need to deal with the explicit dependence of the velocity on the coordinates \(\xi_2\) and \(\xi_3\) normal to the reference axis. The divergence of the stress field in the variational statement is expanded for the beam as follows:

(15.52)\[\begin{split}\begin{aligned} \int_V t^{ij} d_{ij} \mathrm{d}v & \int\limits_{-1/2}^{+1/2} \int\limits_A t^{ij} v_{i,j}(\xi_1) l \mathrm{d}a \mathrm{d}\xi_1 \\ &-\int\limits_{-1/2}^{+1/2} \int\limits_A t^{ij} \Bigg[ \varepsilon_{imn} \rho^m_{,j} \omega^n (\xi_1) + \varepsilon_{imn} \rho^m \omega^n_{,j} (\xi_1) \Bigg] \,l\,\mathrm{d}a\mathrm{d}\xi_1 . \end{aligned}\end{split}\]

The dependence on \(\xi_2\) and \(\xi_3\) is explicit since \(J\) specializes to

(15.53)\[J = A \varepsilon_{rst} \frac{\partial x_r}{\partial \xi_1} m_s n_t = Al ,\]

where \(l\) is the length of the beam, \(A\) its (constant) cross-sectional area, and \(m_s\) and \(n_t\) are the unit vectors along the \(\xi_2\) and \(\xi_3\) axes, respectively.

At this point, we can write the classical force and bending stress resultants \(\mathcal{N}^{ij}\) and \(\mathcal{M}^i_j\) as:

(15.54)\[\begin{split}\begin{aligned} \mathcal{N}^{ij} &= \int\limits_A t^{ij} \mathrm{d}a ,\\ \mathcal{M}^i_j &= \int\limits_A t^{in} \varepsilon_{nmj} \rho^m \mathrm{d}a . \end{aligned}\end{split}\]

Now, we introduce the average stresses \(\tau^i{j}\) and average bending stresses \(\mu^i_j\) as:

(15.55)\[\begin{split}\begin{aligned} \tau^{ij} &= \frac{1}{A} \mathcal{N}^{ij} , \\ \mu^I_j &= \frac{1}{A} \mathcal{M}^i_j . \end{aligned}\end{split}\]

Combining (15.52) through (15.55) yields a reduced divergence of the stress field:

(15.56)\[\begin{split}\begin{aligned} \int_V t^{ij} d_{ij} \mathrm{d}v & \int\limits_{-1/2}^{+1/2} \tau^{ij} v_{i,j}(\xi_1) l A \mathrm{d}\xi_1 \\ &-\int\limits_{-1/2}^{+1/2} \left[ \tau^{ij} \varepsilon_{imn} \rho^m_{,j} \omega^n (\xi_1) + \mu^j_n \omega^n_{,j} (\xi_1) \right] l A \mathrm{d}\xi_1 . \end{aligned}\end{split}\]

The integrals in the reduced divergence of the stress field in the element are evaluated using a mean stress, thereby considering only a state of constant axial, bending, and torsional stress within the element. Expressing (15.56) explicitly with mean kinematic quantities \(\bar{v}_{i,j}\), \(\bar{\omega}^n\) and \(\bar{\omega}^n_{,j}\), and mean stresses \(\bar{\tau}^{ij}\) and \(\bar{\mu}^j_n\), yields:

(15.57)\[\int_V t^{ij} d_{ij} \mathrm{d}v = V \left( \bar{\tau}^{ij} \bar{v}_{i,j} + \bar{\tau}^{ij} \varepsilon_{imn} \rho^m_{,j} \bar{\omega}^n + \bar{\mu}^j_n \bar{\omega}^n_{,j} \right) ,\]

where the mean kinematic quantities are defined by integrating over the element as follows:

(15.58)\[\begin{split}\begin{aligned} \bar{v}_{i,j} &= \frac{1}{V} \int_V v_{i,j} \mathrm{d}v ,\\ \bar{\omega}^n &= \frac{1}{V} \int_V \omega^n \mathrm{d}v , \\ \bar{\omega}^n_{,j} &= \frac{1}{V} \int_V \omega^n_{,j} \mathrm{d}v . \end{aligned}\end{split}\]

The gradient operator is defined by

\[ \begin{align}\begin{aligned}:label: element_formulations:eq:b18\\B_{iI} =\int\limits_{-1/2}^{+1/2} N_{I,j} J \mathrm{d}\xi_1 ,\end{aligned}\end{align} \]

and an averaging operator is defined by

\[ \begin{align}\begin{aligned}:label: element_formulations:eq:b19\\A_I =\int\limits_{-1/2}^{+1/2} N_I J \mathrm{d}\xi_1 .\end{aligned}\end{align} \]

With these definitions, the mean velocity gradient, mean rotational rate, and mean rotational rate gradient can be expressed in the more convenient form:

(15.59)\[\begin{split}\begin{aligned} \bar{v}_{i,j} &= \frac{1}{V} v_{iI} B_{jI} , \\ \bar{\omega}^n &= \frac{1}{V} \omega^n_I A_I , \\ \bar{\omega}^n_{,j} &= \frac{1}{V} \omega^n_I B_{jI} . \end{aligned}\end{split}\]

Thus, the divergence of the stress field becomes:

(15.60)\[\int_V t^{ij} d_{ij} \mathrm{d}v = v_{iI} \left( \bar{\tau}^{ij} B_{jI} \right) - \varepsilon_{imn} \rho^m_{,j} \omega^n_I \left( \bar{\tau}^{ij} A_I \right) - \omega^n_I \left( \bar{\mu}^j_n B_{jI} \right) ,\]

where, evident by inspection of (15.60), the nodal forces due to the divergence of the stress field are then given by:

(15.61)\[f_{iI}= \bar{\tau}^{ij} B_{jI} ,\]

and the nodal torques by:

(15.62)\[m_{nI} = - \varepsilon_{imn} \rho^m_{,j} \bar{\tau}^{ij} A_I - \bar{\mu}^j_n B_{jI} .\]

15.8.3. Evaluation of Stress Resultants

The constant axial, bending, and torsional stress resultant assumptions result in a mean gradient operator and an averaging operator that select mean strain rates linear over the cross-section of the beam. Material non-linearities, though, will result in the stress distribution over the cross to be anything but linear, e.g., in the case of plasticity. As a consequence, the integrals for the force and bending stress resultants are implemented using numerical quadrature over the cross-section. The location of these integration points for the different cross-sections supported are shown in the Theusersguide.

At each integration point, the strain rate is computed from the nodal velocities and rotation rates. The material constitutive behavior is also incrementally evaluated. With a weighting factor and distance from the reference axis for each integration point, the stress resultant integrals are computed simply as a weighted-sum over all integration points. Finally, the stress resultant integrals include the optional offset of the reference axis from the geometric centroid of the cross-section. Details of how the cross-section is specified and how the reference axis is offset are discussed in the Sierra/SM User Manual.

15.8.4. Bending Performance

A correction of the strain energy in the bending of thick beams is necessary due to the overestimation of the transverse shear contributions. This correction of \(\frac{4}{5}\) (Reference [[11]]) is related to the discrepancy between the constant distributions of transverse shear strains implied by the displacement assumptions of the beam and the parabolic through distribution.

In the limit of reducing cross-sectional area, a beam theory with transverse shear becomes arbitrarily stiff in transverse shear response and the transverse shear strains should vanish. Without any correction, the result is a \(\frac{1}{h^2}\) (and \(\frac{1}{w^2}\) ) growth in the transverse shear strain energy, known as shear locking. If \(l\) is the length of the beam element, transverse shear strain energy scale factors of \(\frac{6h^2}{l^2}\) (and \(\frac{6w^2}{l^2}\)) provide, in the limiting case of slender beam behavior, quadratic displacement convergence to the Kirchhoff bending result without the shear locking in the element. Implementation of the shear locking correction factors is done by considering the minimum of \(\frac{4}{5}\) and \(\frac{6h^2}{l^2}\) (and \(\frac{4}{5}\) and \(\frac{6w^2}{l^2}\)), thus allowing a transition from the transverse shear corrected thick beam to the vanishing transverse shear strains \(d_{xz}\), \(d_{yz}\) (implying \(v_{z,x}=\omega_y\) and \(-v_{y,x} = \omega_x\)) required for the thin slender behavior.

15.9. 3D Spring Element

The 2-noded 3D spring element is a simple beam formulation that includes concepts embodied in Timoshenko beam theory.

15.10. Superelement

The superelement formulation in Sierra/SM conforms to the Craig-Bampton reduction capability in Sierra/SD. An option in Sierra/SM is a corotational formulation, which uses the Kabsch algorithm [[12]] to minimize the root mean square deviations between model coordinates and current coordinates of the superelement connection nodes. Additionally, this superelement formulation supports uniform gravity load and uniform initial velocity, with the latter satisfying a zero modal velocity condition, \(\nu_{i}=0\).

---

[1]

D.P. Flanagan and T. Belytschko. Simultaneous relaxation in structural dynamics. Journal of the Engineering Mechanics Division, ASCE, 107:1039–1055, 1981.

[2] (1,2)

T.J.R. Hughes. The Finite Element Method: Linear static and dynamic finite element analysis. Dover, 2000. Reprint of “The Finite Element Method”, Prentice-Hall, 1987.

[3] (1,2,3)

J.T. Ostien, J.W. Foulk, A. Mota, and M.G. Veilleux. A 10-node composite tetrahedral finite element for solid mechanics. International Journal for Numerical Methods in Engineering, 2016.

[4] (1,2)

P. Thoutireddy, J.F. Molinari, E.A. Repetto, and M. Ortiz. Tetrahedral composite finite elements. International Journal for Numerical Methods in Engineering, 53(6):1337–1351, 2002.

[5] (1,2)

Y.Guo, M. Ortiz, T. Belytschko, and E.A. Repetto. Triangular composite finite elements. International Journal for Numerical Methods in Engineering, 47(1-3):287–316, 2000.

[6]

T. Belytschko and C.-S Tsay. Explicit algorithms for nonlinear dynamics of shells. Nonlinear finite element analysis of plates and shells, pages 209–231, 1981.

[7]

S.W. Key and C.C. Hoff. An improved constant membrane and bending stress shell element for explicit transient dynamics. Computer Methods in Applied Mechanics and Engineering, 124(1–2):33–47, 1995. doi:10.1016/0045-7825(95)00785-Y.

[8]

T. Belytschko and I. Leviathan. Physical stabilization of the 4-node shell element with one point quadrature. Computer Methods in Applied Mechanics and Engineering, 113:321–350, 1994.

[9]

T. Belytschko and I. Leviathan. Projection schemes for one-point quadrature shell elements. Computer Methods in Applied Mechanics and Engineering, 115:277–286, 1994.

[10] (1,2)

P. Madabhusi-Raman and J.F. Davalos. Static shear correction factor for laminated rectangular beams. Composites Part B: Engineering, 27(3):285–293, 1996. doi:10.1016/1359-8368(95)00014-3.

[11] (1,2)

G.R. Cowper. The shear coefficient in Timoshenko’s beam theory. Journal of Applied Mechanics, 33:335–340, 1966.

[12]

W. Kabsch. A solution for the best rotation to relate two sets of vectors. Acta Crystallographica Section A, 32(5):922–923, Sep 1976. doi:10.1107/S0567739476001873.