General Formulation of Structural-Element Logic

The low-level implementation details of the structural-element logic are presented here. This information need only be read by users interested in gaining a thorough understanding of the implementation (perhaps, in order to implement nonstandard interaction between structural elements and the grid). The structural-element logic has been developed in a general and flexible fashion in order to support various (as yet unforeseen) application scenarios. The formulation of nonlinear structural elements that provide plasticity in beams and shells is provided in Potyondy (2022).

Finite Elements

The elastic response of each structural element is described by the stiffness that it provides to the surrounding system by means of its nodes. Each structural element has an associated stiffness matrix that is formed using finite-element techniques. The stiffness matrix \([\bf k]\) relates generalized nodal displacements (translations and rotations) \(\{\bf d \}\) to generalized nodal forces (forces and moments) \(\{\bf F \}\) by means of the relation \(\{ {\bf F} \} = [{\bf k}] \{ {\bf d} \}\), in which \([\bf k]\) is a square symmetric matrix of a size equal to the number of degrees of freedom associated with the finite element.

Seven types of finite elements are available: two 2-noded beam finite elements, and five 3-noded shell finite elements. These elements are described in the following two subsections. Beam and pile structural elements utilize the beam finite elements; cable structural elements utilize a degenerate form of the beam finite element in which only the local axially directed degrees-of-freedom are employed; and, shell, geogrid and liner structural elements utilize the shell finite elements. Note that element-to-grid interaction occurs by means of the links that connect the nodes used by elements to the zones.

Beam Finite Elements

The beam finite elements are two-noded, straight, finite elements with six degrees of freedom per node: three translational components and three rotational components (see Beam-Type Coordinate System). There are two beam elements. The Euler-Bernoulli beam element that is appropriate for modeling long thin beams, and the Timoshenko beam element that is appropriate for modeling short thick beams. Both elements are straight; a curved beam can be modeled by using a collection of these elements. Plastic material behavior is supported only by the Euler-Bernoulli beam element. The plasticity is induced by axial and bending deformations; twisting deformation induces an elastic response. For the Euler-Bernoulli beam element, the rotational degree of freedom defines the nodal value of both beam slope and cross-section rotation. For the Timoshenko beam element, beam slope and cross-section rotation differ, and thus, the rotational degree of freedom is regarded as cross-section rotation.

The beam element stiffness matrices (Cook et al., 2002, Eq. 2.3-8; Potyondy, 2022, Eq. 82) model the behavior of a beam that is straight, prismatic and symmetrical about both principal cross-sectional axes (bisymmetrical), and composed of isotropic elastic material. Stiffness matrices are not formed when modeling plastic material behavior; instead, the internal element forces are updated by integrating the internal stresses over the beam volume using a numerical integration scheme. The stiffness-matrix formulation considers displacements resulting from uniform axial deformation, flexural deformation, and twisting deformation. It neglects displacements resulting from the out-of-plane (longitudinal) warping of a cross section that torsional forces may cause. The effect of transverse shearing deformations is also neglected by the Euler-Bernoulli element, but is accounted for by the Timoshenko element. According to McGuire and Gallagher (1979), the equations embodied by the Euler-Bernoulli stiffness matrix are sufficient for the analysis of the large majority of the frameworks encountered in structural-engineering practice.

Shell Finite Elements

There are five 3-noded, triangular finite elements available for shell-type elements: 2 plane-stress elements, 1 plate-bending element and 2 shell elements. The plane-stress elements model membrane action only, the plate-bending element models bending action only, and the shell elements model general shell behavior as a superposition of membrane and bending actions (see Bending Action and Membrance Action). These elements operate upon the 18 degrees of freedom available to the shell-type elements (see Shell-Type Coordinate System). Plastic material behavior is supported by the DKT, CST and DKT-CST finite elements.

We provide the following general recommendation regarding choice of shell finite element. For problems in which a shell element is being connected to zones, use the DKT-CST element. For all other problems, use the DKT-CST hybrid element. These recommendations are explained at the end of the section below on the DKT-CST shell element.

CST Plane-Stress Element (6 degrees of freedom)

The CST plane-stress element is a three-noded plane-stress finite element with two translational degrees-of-freedom per node (shown in Figure 3 (a) by removing the rotational degree-of-freedom). The CST element is described, and an explicit formulation of \([{\bf k}_{\rm CST}]\) is given, by Cook et al. (1989).

CST Hybrid Plane-Stress Element (9 degrees of freedom)

The CST hybrid plane-stress element is a three-noded plane-stress finite element with three degrees of freedom per node: two translational components and one rotational component (shown in Figure 3 (a)). The CST hybrid plane-stress element is formulated by the assumed-stress hybrid approach. It provides a more competent membrane response than a constant strain triangle (CST) through the inclusion of “drilling freedoms,” \(\theta_z\), among its degrees of freedom. The CST hybrid element is described, and an explicit formulation of \([{\bf k}_{\rm CSTH}]\) is given by Cook (1987).

DKT Plate-Bending Element (9 degrees of freedom)

The DKT plate-bending element is a three-noded finite element with three degrees of freedom per node: one translational component and two rotational components (shown in Figure 3 (b)). The DKT plate-bending element enforces the Kirchhoff plate theory assumptions at six discrete locations (corners and mid-sides) within the element; thus, the designation DKT for Discrete Kirchhoff Triangle. The element is described by Batoz et al. (1980), who conclude that it is one of the most efficient and reliable three-node (9 degrees of freedom) plate-bending elements available at the time of publication. An explicit formulation of \([{\bf k}_{\rm DKT}]\) is given by Batoz (1982). The ability of the DKT element to model bending action is quantified in verification problems 1, 2 and 3.

DKT-CST Hybrid Shell Element (18 degrees of freedom)

The DKT-CST hybrid shell element is a three-noded, flat, thin-shell finite element with six degrees of freedom per node: three translational components and three rotational components (shown in Figure 3). It combines a nine degrees-of-freedom CST hybrid plane-stress element to model membrane action in the shell with a nine degrees-of-freedom DKT plate-bending element to model bending action in the shell. This composite flat shell element is referred to as a DKT-CST hybrid, and is described by Cook et al. (1989, pp. 351-352).

Let a typical element lie in the \(xy\)-plane of a local coordinate system \(xyz\). Nodal degrees-of-freedom are shown in Figure 3, where \(\lfloor \theta_x \quad \theta_y \rfloor = \lfloor {\partial w \over \partial y} \quad -{\partial w \over \partial x} \rfloor\). Let these degrees of freedom be called \(\{{\bf d'} \}\), and be arranged in the order

\[\begin{Bmatrix} {\bf d^{\prime}} \end{Bmatrix} = {\lfloor {\bf u}_i \; {\bf v}_i \; {\boldsymbol\theta}_{zi} \; {\bf w}_i \; {\boldsymbol\theta}_{xi} \; {\boldsymbol\theta}_{yi} \rfloor}^T\]

where \(\lfloor {\bf u}_i \rfloor = \lfloor u_1\quad u_2 \quad u_3 \rfloor\), \(\lfloor {\boldsymbol\theta}_{xi} \rfloor = \lfloor \theta_{x1}\quad \theta_{x2} \quad \theta_{x3} \rfloor\), and so on. For the DKT-CST hybrid shell element, whose stiffness matrix in local \(xyz\)-coordinates is called \([{\bf k'}]\), we have

(1)\[\begin{split}\begin{bmatrix} {\bf k^{\prime}} \end{bmatrix} \begin{Bmatrix} {\bf d^{\prime}} \end{Bmatrix} = \begin{bmatrix} \begin{matrix} \begin{bmatrix} {\bf k}_{CSTH} \end{bmatrix} \\ 9\times 9 \end{matrix} & \begin{matrix} \begin{bmatrix} {\bf 0} \end{bmatrix} \\ 9\times 9 \end{matrix} \\ \begin{matrix} \begin{bmatrix} {\bf 0} \end{bmatrix} \\ 9\times 9 \end{matrix} & \begin{matrix} \begin{bmatrix} {\bf k}_{DKT} \end{bmatrix} \\ 9\times 9 \end{matrix} \end{bmatrix} \begin{Bmatrix} {\bf u}_i \\ {\bf v}_i \\ {\boldsymbol\theta}_{zi} \\ {\bf w}_i \\ {\boldsymbol\theta}_{xi} \\ {\boldsymbol\theta}_{yi} \end{Bmatrix}\end{split}\]

where \([{\bf k}_{\rm CSTH}]\) is the stiffness matrix of a CST hybrid plane-stress element, and \([{\bf k}_{\rm DKT}]\) is the stiffness matrix of a DKT plate-bending element.

The following observations about the DKT-CST hybrid element are appropriate (also see the observations about the DKT-CST element in the next section).

  1. \([{\bf k}_{\rm CSTH}]\) only contains degrees of freedom associated with membrane action, and \([{\bf k}_{\rm DKT}]\) only contains degrees of freedom associated with bending action. Although membrane-bending coupling is present throughout an actual curved shell, it is absent in an individual element (as evidenced by the 9 by 9 null matrices of (1)). Nevertheless, the element works well enough to be competitive with curved elements, as demonstrated by Carpenter et al. (1986).

  2. The element provides a solution which converges to the classical Kirchhoff thin-plate solution. Thus, the element is appropriate for modeling thin shells in which the displacements caused by transverse-shearing deformations can be neglected.

../../../../../../_images/dkt-cst-element.png

Figure 3: DKT-CST hybrid flat shell finite element in a local xy-plane: (a) degrees of freedom associated with membrane action (CST hybrid sub-element); (b) degrees of freedom associated with bending action (DKT sub-element)

DKT-CST Shell Element (15 degrees of freedom)

The DKT-CST shell element is a three-noded, flat, thin-shell finite element with five degrees of freedom per node: three translational components and two rotational components. It combines a six degrees-of-freedom CST plane-stress element to model membrane action in the shell with a nine degrees-of-freedom DKT plate-bending element to model bending action in the shell. This composite flat shell element is referred to as a DKT-CST, and is described by Cook et al. (1989, pp. 351-352).

Let a typical element lie in the \(xy\)-plane of a local coordinate system \(xyz\). Nodal degrees-of-freedom are shown in Figure 3 by removing \(\theta_z\), where \(\lfloor \theta_x \quad \theta_y \rfloor = \lfloor {\partial w \over \partial y} \quad -{\partial w \over \partial x} \rfloor\). Let these degrees of freedom be called \(\{{\bf d^\star} \}\), and be arranged in the order

\[\begin{Bmatrix} {\bf d^{\star}} \end{Bmatrix} = {\lfloor {\bf u}_i \, {\bf v}_i \, {\bf w}_i \, {\boldsymbol\theta}_{xi} \, {\boldsymbol\theta}_{yi} \rfloor}^T\]

where \(\lfloor {\bf u}_i \rfloor = \lfloor u_1\quad u_2 \quad u_3 \rfloor\), \(\lfloor {\boldsymbol\theta}_{xi} \rfloor = \lfloor \theta_{x1}\quad \theta_{x2} \quad \theta_{x3} \rfloor\), and so on. For the DKT-CST shell element, whose stiffness matrix in local \(xyz\)-coordinates is called \([{\bf k^\star}]\), we have

(2)\[\begin{split}\begin{bmatrix} {\bf k^{\star}} \end{bmatrix} \begin{Bmatrix} {\bf d^{\star}} \end{Bmatrix} = \begin{bmatrix} \begin{matrix} \begin{bmatrix} {\bf k}_{CST} \end{bmatrix} \\ 6\times 6 \end{matrix} & \begin{matrix} \begin{bmatrix} {\bf 0} \end{bmatrix} \\ 6\times 9 \end{matrix} \\ \begin{matrix} \begin{bmatrix} {\bf 0} \end{bmatrix} \\ 9\times 6 \end{matrix} & \begin{matrix} \begin{bmatrix} {\bf k}_{DKT} \end{bmatrix} \\ 9\times 9 \end{matrix} \end{bmatrix} \begin{Bmatrix} {\bf u}_i \\ {\bf v}_i \\ {\bf w}_i \\ {\boldsymbol\theta}_{xi} \\ {\boldsymbol\theta}_{yi} \end{Bmatrix}\end{split}\]

where \([{\bf k}_{\rm CST}]\) is the stiffness matrix of a CST plane-stress element, and \([{\bf k}_{\rm DKT}]\) is the stiffness matrix of a DKT plate-bending element.

The following observations about the DKT-CST element are appropriate.

  1. \([{\bf k}_{\rm CST}]\) only contains degrees of freedom associated with membrane action, and \([{\bf k}_{\rm DKT}]\) only contains degrees of freedom associated with bending action. Although membrane-bending coupling is present throughout an actual curved shell, it is absent in an individual element (as evidenced by the null matrices of (2).

  2. The element provides a solution which converges to the classical Kirchhoff thin-plate solution. Thus, the element is appropriate for modeling thin shells in which the displacements caused by transverse-shearing deformations can be neglected.

  3. The element differs from the DKT-CST hybrid element by not having drilling degrees-of-freedom. This makes the element overly stiff for problems with large membrane-stress gradients. For such problems, the DKT-CST hybrid element will perform better. For problems without large membrane-stress gradients, the performance of the DKT-CST element will be much improved, making it competitive with the DKT-CST hybrid element. The performance of the DKT-CST and DKT-CST hybrid elements is compared in this verification problem.

  4. The lack of drilling degrees-of-freedom makes this element more suitable for coupling with the model grid than the DKT-CST hybrid element. The zones have a linear displacement variation, and the 3 drilling degrees-of-freedom give the DKT-CST hybrid a quadratic in-plane displacement variation along its edges. When a shell-type element is rigidly connected to the grid, only the translational degrees-of-freedom are slaved to move with the grid, while the rotational degrees-of-freedom are left free. Because we are not constraining the drilling degrees-of-freedom, it is possible for gaps to form between the element and the grid. These incompatibilities seem to be activated when shell-type elements of different sizes meet along a common line. The DKT-CST element has a linear in-plane displacement variation along its edges, which is fully compatible with that of the zones and, thus, gaps cannot form when these elements are rigidly connected to zones.

  1. The two preceding observations lead to the following general recommendation regarding choice of shell finite element. For problems in which a shell element is being connected to zones, use the DKT-CST element. For all other problems, use the DKT-CST hybrid element.

Stress Recovery Procedure (in detail)

The stress-recovery procedure is outlined in Stress Recovery Procedure. A more detailed description of the process is provided here and in Potyondy (2022). Stress recovery is applied to a group of shell-type elements. Nodal averaging of bending and membrane stress resultants only occurs for elements within this group. Stress recovery consists of two steps:

  1. Compute stress-resultant field.

    1. Define the surface coordinate system at each node.

      This system has its \(z\)-direction equal to the average \(z\)-direction of all shell-type elements using the node, and has its \(x\)-direction equal to the projection of a user-specified surface vector onto the \(z\)-plane.

    2. Compute the bending and membrane stress resultants

      (\(M_x\), \(M_y\), \(M_{xy}\), \(N_x\), \(N_y\), and \(N_{xy}\)) at each node. These values are found by first computing them within each shell-type element at its nodes, and then averaging the contributions from all shell-type elements that use each node.

      Note that these six stress resultants vary linearly over each element.

      The bending and membrane stress resultants within each shell-type element with a plastic material model are found by integrating the stress through the element thickness using Gauss quadrature. The bending and membrane stress resultants within each shell-type element with an elastic material model are found using the flexibility-stiffness transformation from McGuire and Gallagher (1979, pp. 76-77) to determine which displacements are consistent with a given set of nodal forces acting on the element. This is done by forming a flexibility matrix that is derived from the element stiffness matrix. The flexibilities are derived from the stiffnesses by defining a stable, statically determinate support system, removing from the stiffness matrix the rows and columns corresponding to the support components, and inverting the remainder. These displacements are then fed into stress-recovery procedures appropriate to the particular finite element.

  1. Compute the transverse-shear stress resultants (\(Q_x\) and \(Q_y\)) within each shell-type element by enforcing the equilibrium equations

    \[ \begin{align}\begin{aligned}Q_x = {\partial M_x \over \partial x} + {\partial M_{xy} \over \partial y}\\Q_y = {\partial M_y \over \partial y} + {\partial M_{xy} \over \partial x}\end{aligned}\end{align} \]

    by assuming that the bending resultant field varies linearly over each element. Note that \(Q_x\) and \(Q_y\) are constant over each shell-type element, and are not averaged at nodes.

  1. Compute stress field at specified shell depth using Stresses From Resultants, and then rotating the stresses into the global system. This step is only done for elements with an elastic material model.

Governing Equations

The discretization and time integration of the equations that govern the dynamic response of the structural elements and their associated entities (i.e., nodes and links) utilize standard finite-element practice. The description given here is similar to that given in Chapter 13 of Cook et al. (1989).

The equations that govern the dynamic response of a structure can be derived by requiring the work of external forces to be absorbed by the work of internal, inertial and damping forces for any small, kinematically admissible motion. For a single element, this work balance becomes

(3)\[\begin{split}\int_{V_e} \begin{Bmatrix} \delta {\bf u} \end{Bmatrix}^T \begin{Bmatrix} {\bf F} \end{Bmatrix} dV + \int_{S_e} \begin{Bmatrix} \delta {\bf u} \end{Bmatrix}^T \begin{Bmatrix} {\boldsymbol\Phi} \end{Bmatrix} dS +\ \sum_{i=1}^n \begin{Bmatrix} \delta {\bf u} \end{Bmatrix}_i^T \begin{Bmatrix} {\bf p} \end{Bmatrix}_i = \\ \qquad \qquad \int_{V_e} \bigg( \begin{Bmatrix} \delta {\boldsymbol\varepsilon} \end{Bmatrix}^T \begin{Bmatrix} {\boldsymbol\sigma} \end{Bmatrix} + \begin{Bmatrix} \delta {\bf u} \end{Bmatrix}^T \rho \begin{Bmatrix} {\ddot{\bf u}} \end{Bmatrix} + \begin{Bmatrix} \delta {\bf u} \end{Bmatrix}^T \kappa_d \begin{Bmatrix} {\dot{\bf u}} \end{Bmatrix} \bigg) dV\end{split}\]

where \(\{ \delta {\bf u} \}\) and \(\{ \delta {\boldsymbol\varepsilon} \}\) are, respectively, small arbitrary displacements and their corresponding strains, \(\{ {\bf F} \}\) are body forces, \(\{ {\boldsymbol\Phi} \}\) are prescribed surface tractions that act over a portion of the element surface \(S_e\), \(\{ {\bf p}_i \}\) are concentrated loads that act at a total of \(n\) points on the element, \(\{ {\bf u} \}_i^T\) is the displacement of the point at which load \(\{ {\bf p}_i \}\) is applied, \(\{ {\boldsymbol\sigma} \}\) are stresses, \(\rho\) is the mass density of the material, \(\kappa_d\) is a material-damping parameter analogous to viscosity, and volume integration is carried out over the element volume, \(V_e\).

If we express the displacement field \(\{ {\bf u} \}\) (which is a function of both space and time) and its first two time derivatives as

(4)\[\{ {\bf u} \} = [{\bf N}] \{ {\bf d} \} \quad\quad \{ \dot{\bf u} \} = [{\bf N}] \{ \dot{\bf d} \} \quad\quad \{ \ddot{\bf u} \} = [{\bf N}] \{ \ddot{\bf d} \}\]

where shape functions \([{\bf N}]\) are functions of space only, and nodal degrees-of-freedom \(\{{\bf d}\}\) are functions of time only, then (4) represent a local separation of variables. By substituting (4) into (3), and noting that \(\{ \delta{\bf u} \}\), and thus \(\{ \delta{\bf d} \}\), are arbitrary, one obtains

(5)\[[{\bf m}] \{ \ddot{\bf d} \} + [{\bf c}] \{ \dot{\bf d} \} + \{ {\bf r}^{\rm int} \} = \{ {\bf r}^{\rm ext} \}\]

where the element mass and damping matrices are defined as

\[[{\bf m}] = \int_{V\!{\scriptscriptstyle e}} \rho [{\bf N}]^T [{\bf N}]\,dV\]
(6)\[[{\bf c}] = \int_{V\!{\scriptscriptstyle e}} \kappa_{\scriptscriptstyle d} [{\bf N}]^T [{\bf N}]\,dV\]

and the element internal force and external load vectors are defined as

(7)\[\begin{Bmatrix} {\bf r}^{int} \end{Bmatrix} = \int_{V_e} \begin{bmatrix} {\bf B} \end{bmatrix}^T \begin{Bmatrix} {\boldsymbol\sigma} \end{Bmatrix} dV\]
(8)\[\begin{Bmatrix} {\bf r}^{ext} \end{Bmatrix} = \int_{V_e} \begin{bmatrix} {\bf N} \end{bmatrix}^T \begin{Bmatrix} {\bf F} \end{Bmatrix} dV + \int_{S_e} \begin{bmatrix} {\bf N} \end{bmatrix}^T \begin{Bmatrix} {\boldsymbol\Phi} \end{Bmatrix} dS + \sum_{i=1}^{n} \begin{bmatrix} {\bf p} \end{bmatrix}_i\]

where \([{\bf B}]\) is the strain-displacement matrix that satisfies \(\{\delta {\boldsymbol\varepsilon}\} = [{\bf B}]\{{\bf d}\}\).

(5) is a system of coupled, second-order, ordinary differential equations in time, and is called a finite-element semi-discretization because, although displacements \([{\bf d}]\) are discrete functions of space, they are still continuous functions of time. The structural-element logic employs an explicit, direct-integration method which discretizes (5) in time, using the central-difference method to obtain a sequence of uncoupled algebraic equations.

Recall that (5) was written for a single element. Equations for the assembled structural system can be written as

(9)\[\begin{bmatrix} {\bf M} \end{bmatrix} \begin{Bmatrix} \ddot{\bf D} \end{Bmatrix} + \begin{bmatrix} {\bf C} \end{bmatrix} \begin{Bmatrix} \dot{\bf D} \end{Bmatrix} = \begin{Bmatrix} {\bf R}^{ext} \end{Bmatrix} - \begin{Bmatrix} {\bf R}^{int} \end{Bmatrix} = \begin{Bmatrix} {\bf A} \end{Bmatrix}\]

where \(\{{\bf A}\}\) designates the accumulated out-of-balance force, and uppercase English alphabet characters designate structure matrices which are constructed by the conceptual expansion of element matrices to “structure size,” followed by addition of overlapping coefficients. See any standard finite-element text (e.g., Cook et al. 1989) for a detailed description of the assembly procedure.

The equations of motion, (9), are written for a specific instant of time:

(10)\[\begin{bmatrix} {\bf M} \end{bmatrix} \begin{Bmatrix} \ddot{\bf D} \end{Bmatrix}_n + \begin{bmatrix} {\bf C} \end{bmatrix} \begin{Bmatrix} \dot{\bf D} \end{Bmatrix}_n = \begin{Bmatrix} {\bf A} \end{Bmatrix}_n\]

where subscript \(n\) denotes time \(n\Delta t\), and \(\Delta t\) is the size of the time increment or timestep.

Central-Difference Approximations

The central-difference approximations for acceleration at time \(n\Delta t\), velocity at time \((n - {\scriptstyle 1/2})\Delta t\), and velocity at time \(n\Delta t\) are given by

(11)\[\{ \ddot{\bf D} \}_n = {1\over\Delta t} \biggl( \{ \dot{\bf D} \}_{n \scriptscriptstyle + 1/2} - \{ \dot{\bf D} \}_{n \scriptscriptstyle - 1/2} \biggr) + O( \Delta t^2 )\]
(12)\[\{ \dot{\bf D} \}_{n \scriptscriptstyle - 1/2} = {1\over\Delta t} \biggl( \{ {\bf D} \}_{n} - \{ {\bf D} \}_{n \scriptscriptstyle - 1} \biggr) + O( \Delta t^2 )\]
(13)\[\{ \dot{\bf D} \}_n = {1\over2} \biggl( \{ \dot{\bf D} \}_{n \scriptscriptstyle - 1/2} + \{ \dot{\bf D} \}_{n \scriptscriptstyle + 1/2} \biggr) + O( \Delta t^2 )\]

These expressions are second-order accurate, which implies that halving the timestep should approximately quarter the error in the approximation.

Accumulated Out-of-Balance Force Vector

The structure internal force vector that makes up one portion of the accumulated out-of-balance force vector \(\{{\bf A}\}_n\) in (10) is obtained by assembling the contributions from each structural element and from each deformable direction of each link

\[\{ {\bf R}^{\rm int} \}_n = \sum_{j} \left( \{ {\bf r}^{\rm int} \}_n \right)_j\]

with

\[ \begin{align}\begin{aligned}\{ {\bf r}^{\rm int} \}_n &= \{ {\bf r}^{\rm int} \}_{n \scriptscriptstyle - 1} + \{ \Delta {\bf r}^{\rm int} \}\\= \{ {\bf r}^{\rm int} \}_{n \scriptscriptstyle - 1} + [{\bf k}]_{n \scriptscriptstyle - 1/2} \bigl( \{ {\bf d} \}_n - \{ {\bf d} \}_{n \scriptscriptstyle - 1} \bigr)\\= \{ {\bf r}^{\rm int} \}_{n \scriptscriptstyle - 1} + [{\bf k}]_{n \scriptscriptstyle - 1/2} \{ \dot{\bf d} \}_{n \scriptscriptstyle - 1/2} \Delta t\end{aligned}\end{align} \]

in which the central-difference approximation of (12) has been employed. The internal-force vector \(\{ {\bf r}^{\rm int} \}\) that is maintained by each structural element and each deformable direction of each link is also updated during the computation of \(\{ {\bf R}^{\rm int} \}\).

After distributing all of the internal forces to the appropriate nodes, each rigidly attached direction, \({\bf n}^{[{\rm S}]}\), of each link contributes the following portion of the accumulated out-of-balance force of its source node \({\bf A}^{[{\rm S}]}\) to the accumulated out-of-balance force \({\bf A}^{[{\rm T}]}\) of its target entity as

\[{\bf A}^{[{\rm T}]} \leftarrow {\bf A}^{[{\rm T}]} + \left( {\bf A}^{[{\rm S}]} \cdot {\bf n}^{[{\rm S}]} \right) {\bf n}^{[{\rm S}]}\]

Mass Matrix

The structure mass matrix in (10) is obtained by utilizing a nodal-lumping scheme, and assembling the contributions from each structural element

\[[ {\bf M} ] = \lceil {\bf M} \rfloor = \sum_{j=1}^{\rm numel} \lceil {\bf m} \rfloor_j\]

with (in \(3D\))

\[\begin{split}\lceil {\bf m} \rfloor_j = \begin{cases} \frac{\rho_j V_j}{n_j}, & \text{all 3 translational d.o.f.} \\ \frac{2}{5}\Bigg(\frac{\rho_j V_j}{n_j}\Bigg)\Bigg(\frac{3V_j}{4\pi n_j}\Bigg)^{2/3}, & \text{all 3 rotational d.o.f.} \end{cases}\end{split}\]

and (in \(2D\))

\[\begin{split}\lceil {\bf m} \rfloor_j = \begin{cases} \frac{\rho_j V_j}{n_j}, & \text{2 translational d.o.f.} \\ \frac{1}{24}(\rho_j V_j {L_j}^2), & \text{1 rotational d.o.f.} \end{cases}\end{split}\]

where \(\rho_j\), \(V_j\) and \(n_j\) are the mass density, volume and number of nodes of element \(j\). The rotational inertias are found by assuming that the portion of element mass seen by each node is distributed as a sphere.

After distributing the mass associated with all of the structural elements to the appropriate nodes, each rigidly attached direction of each link contributes the following portion of its source mass in direction- \(k\) \(m_k^{[{\rm S}]}\) to the accumulated mass in target direction-\(i\) \(m_i^{[{\rm T}]}\), using the weighting function

(14)\[m_i^{[{\rm T}]} \leftarrow m_i^{[{\rm T}]} + m_k^{[{\rm S}]} \left( { \left| {\bf n}_k^{[{\rm S}]} \cdot {\bf n}_i^{[{\rm T}]} \right| \over \sum^3_{l=1} \left| {\bf n}_k^{[{\rm S}]} \cdot {\bf n}_l^{[{\rm T}]} \right| } \right)\]

where \({\bf n}_i^{[{\rm S}]}\) and \({\bf n}_i^{[{\rm T}]}\) are, respectively, the unit-vector triads of the source and target local coordinate systems.

Damping Formulation

The derivation of equation (6) utilized a viscous-damping parameter that is suitable for modeling a Newtonian fluid; however, when modeling structural systems, we are less interested in viscous damping than in dry friction and hysteresis loss. These energy-loss mechanisms are not well-understood and, from a practical standpoint, (6) does not correctly represent structural damping. Thus, three ad hoc damping schemes (local damping, combined local damping and Rayleigh or proportional damping) are provided. The implementation details of the three damping schemes are presented here; refer to Mechanical Damping and Rayleigh Damping for a description of the rationale behind each scheme.

The starting point for all three damping schemes is given by rewriting (10) as

(15)\[[{\bf M}] \{ \ddot{\bf D} \}_n + \{ {\bf R}^{\rm d} \}_n = \{ {\bf A} \}_n\]

where \(\{ {\bf R}^{\rm d} \}_n\) represents the force produced by the damping mechanism at time \(n\Delta t\).

Local Damping

Substitution of (11) into (15) yields

(16)\[\{ \dot{\bf D} \}_{n \scriptscriptstyle + 1/2} = \{ \dot{\bf D} \}_{n \scriptscriptstyle - 1/2} + \Delta t [{\bf M}]^{-1} \left( \{ {\bf A} \}_n - \{ {\bf R}^{\rm d} \}_n \right)\]

The use of a diagonal mass matrix (as described in Mass Matrix) uncouples this set of equations such that we can rewrite them for a single degree of freedom \(i\) as

\[( \dot D_i )_{n \scriptscriptstyle + 1/2} = ( \dot D_i )_{n \scriptscriptstyle - 1/2} + {\Delta t \over M_i} \left( (A_i)_n - ( R_i^{\rm d} )_n \right)\]

Local damping is provided by setting

(17)\[( R_i^{\rm d} )_n = \lambda_i \> \bigl| (A_i)_n \bigr| \> {\rm sgn} ( \dot D_i )_n\]

where \(\lambda_i\) is the local-damping factor for degree-of-freedom \(i\), and

\[\begin{split}{\rm sgn}(x) = \begin{cases} +1, & \text{if} x \geq 0 \\ -1, & \text{if} x < 0 \end{cases}\end{split}\]

The local-damping formulation is implemented by computing an undamped \((\dot D_i )^{\rm ud}\) and a damped \(( \dot D_i )^{\rm d}\) new velocity:

\[ \begin{align}\begin{aligned}( \dot D_i )^{\rm ud} = ( \dot D_i )_{n \scriptscriptstyle - 1/2} + {\Delta t \over M_i} (A_i)_n\\( \dot D_i )^{\rm d} = ( \dot D_i )_{n \scriptscriptstyle - 1/2} + {\Delta t \over M_i} \biggl( (A_i)_n - \lambda_i \> | (A_i)_n | \> {\rm sgn} ( \dot D_i )_{n \scriptscriptstyle - 1/2} \biggr)\end{aligned}\end{align} \]

and setting the actual new velocity equal to the damped velocity (only if the velocity has not changed sign) as

(18)\[\begin{split}\begin{pmatrix} \dot{D}_i \end{pmatrix}_{n+1/2} = \begin{cases} \begin{pmatrix} \dot{D}_i \end{pmatrix}^{ud}, & \text{if} \begin{pmatrix} \dot{D}_i \end{pmatrix}_{n-1/2}\begin{pmatrix} \dot{D}_i \end{pmatrix}^d \leq 0 \\ \begin{pmatrix} \dot{D}_i \end{pmatrix}^{d}, & \text{otherwise} \end{cases}\end{split}\]

It can be shown that this implementation maintains second-order accuracy as follows. Assume that

(19)\[( \dot D_i )_{n \scriptscriptstyle + 1/2} ( \dot D_i )_{n \scriptscriptstyle - 1/2} > 0\]

This implies that

\[{\rm sgn}( \dot D_i )_n = {\rm sgn}( \dot D_i )_{n \scriptscriptstyle - 1/2}\]

which makes (18) an exact expression of (17). Note that if the assumption of (19) is not satisfied (i.e., the velocity either changes sign or is zero), then no damping is applied to this degree of freedom.

Combined Local Damping

Combined local damping is provided by utilizing the same procedure as described for local damping, but with (17) replaced by

\[( R_i^{\rm d} )_n = \lambda_i \> \bigl| (A_i)_n \bigr| \> {1\over2} \bigl( {\rm sgn} ( \dot D_i )_n - {\rm sgn} ( \dot A_i )_n \bigr)\]

where

\[{\rm sgn} ( \dot A_i )_n = {\rm sgn} \bigl( ( A_i )_n - ( A_i )_{n-1} \bigr)\]

Rayleigh Damping

Rayleigh damping results by forming the damping matrix \([ {\bf C} ]\) as a linear combination of the stiffness and mass matrices:

\[[ {\bf C} ] = \beta [ {\bf K} ] + \alpha [ {\bf M} ]\]

where \(\beta\) and \(\alpha\) are called, respectively, the stiffness- and mass-proportional damping constants. The relation between \(\beta\), \(\alpha\) and the fraction of critical damping \(\xi\) at frequency \(\omega\) is given by

\[\xi = {1\over2} \left( \beta\omega + {\alpha\over\omega} \right)\]

Damping attributable to stiffness increases with increasing frequency, whereas damping attributable to mass increases with decreasing frequency.

Damping constants \(\beta\) and \(\alpha\) are related to the fractions of critical damping (\(\xi_1\) and \(\xi_2\)) at frequencies (\(\omega_1\) and \(\omega_2\)) by

\[ \begin{align}\begin{aligned}\beta &= 2(\xi_2\omega_2 - \xi_1\omega_1) / (\omega_2^2 - \omega_1^2)\\\alpha &= 2\omega_1\omega_2(\xi_1\omega_2 - \xi_2\omega_1) / (\omega_2^2 - \omega_1^2)\end{aligned}\end{align} \]

The governing equation of motion used when Rayleigh damping is active is derived as follows. Express the force produced by the damping mechanism (\(\{ {\bf R}^{\rm d} \}_n\) in (15)) as the sum of the stiffness and mass contributions:

\[\{ {\bf R}^{\rm d} \}_n = \{ {\bf R}^{\beta} \}_n + \{ {\bf R}^{\alpha} \}_n = \{ {\bf R}^{\beta} \}_n + \alpha [{\bf M}] \{ \dot{\bf D} \}_n\]

and substitute this expression into (15) to obtain

(20)\[[{\bf M}] \{ \ddot{\bf D} \}_n + \alpha [{\bf M}] \{ \dot{\bf D} \}_n = \{ {\bf A} \}_n - \{ {\bf R}^{\alpha} \}_n\]

Substitution of the central-difference approximations given by (11) and (13) into (20), followed by rearrangement of terms, yields

(21)\[\{ \dot{\bf D} \}_{n \scriptscriptstyle + 1/2} = \left( {2 - \alpha\Delta t \over 2 + \alpha\Delta t} \right) \{ \dot{\bf D} \}_{n \scriptscriptstyle - 1/2} + \left( {2\Delta t \over 2 + \alpha\Delta t} \right) [{\bf M}]^{-1} \biggl( \{ {\bf A} \}_n - \{ {\bf R}^{\beta} \}_n \biggr)\]

in which the stiffness contribution to the damping force is obtained by assembling the contributions from each structural element and from each deformable direction of each link:

(22)\[\{ {\bf R}^{\beta} \}_n = \sum_j \left( \{ {\bf r}^{\beta} \}_n \right)_j\]

with

(23)\[ \begin{align}\begin{aligned}\{ {\bf r}^{\beta} \}_n & = \beta [{\bf k}]_{n \scriptscriptstyle - 1/2} \bigl( \{ {\bf d} \}_n - \{ {\bf d} \}_{n \scriptscriptstyle - 1} \bigr) / \Delta t\\& = \beta [{\bf k}]_{n \scriptscriptstyle - 1/2} \{ \dot{\bf d} \}_{n \scriptscriptstyle - 1/2}\end{aligned}\end{align} \]

Computation of Stable Timestep

Time integration carried out using the central-difference method is a conditionally stable process. The stable timestep for our central-difference equations ((16) if local or combined local damping is active, or (20) if Rayleigh damping is active) is given by (Belytschko 1983)

(24)\[\Delta t_{\rm stab} = {2\over\omega_{\rm max}} F(\xi_{\rm max})\quad\hbox{where } F(\xi) = \sqrt{1+\xi^2} - \xi\]

where \(\omega_{\rm max}\) is the highest natural frequency of the discretized system, and \(\xi_{\rm max}\) is the fraction of critical damping at \(\omega_{\rm max}\). The value of \(\xi_{\rm max}\) is dependent upon the damping scheme,

(25)\[\begin{split}\xi_{max} = \begin{cases} \frac{1}{2} ( \beta\omega_{max} + \alpha / \omega_{max} ), & \text{(Rayleigh Damping)} \\ \text{max} \lambda_i / \pi, & \text{(local or combined local damping)} \end{cases}\end{split}\]

The natural frequencies, \(\omega_n\), of the discretized system are given by solution of the eigenproblem,

(26)\[\left( [{\bf K}] - \omega^2 [{\bf M}] \right) \{ \tilde{\bf D} \} = \{ {\bf 0} \}\]

The maximum frequency can be bounded by applying Gerschgorin’s Theorem (Isaacson and Keller 1963) to (26) which, for a diagonal mass matrix, states that

(27)\[\omega_{\rm max}^2 \leq \max_i \left( G_i / M_i \right) \qquad\hbox{with}\quad i = 1,2,\ldots,n_e\]

where \(n_e\) is the number of degrees of freedom, and \(G_i\) is the stiffness sum for the \(i\)th row of the assembled stiffness matrix

\[G_i = K_{ii} + \sum_{j=1 \atop j \neq i}^{n_e} \left| K_{ij} \right|\]

The six components of \(G_i\) are stored at each structural-element node and expressed in terms of the node-local system. Stiffness contributions are provided by all structural elements and by all deformable directions of all links. After distributing all of the stiffnesses to the appropriate nodes, each rigidly attached direction of each link contributes the following portion of its source stiffness in direction-\(k\) \(G_k^{[{\rm S}]}\) to the accumulated stiffness in target direction-\(i\) \(G_i^{[{\rm T}]}\), using the same weighting function as given by (14), namely

\[G_i^{[{\rm T}]} \leftarrow G_i^{[{\rm T}]} + G_k^{[{\rm S}]} \left( { \left| {\bf n}_k^{[{\rm S}]} \cdot {\bf n}_i^{[{\rm T}]} \right| \over \sum^3_{l=1} \left| {\bf n}_k^{[{\rm S}]} \cdot {\bf n}_l^{[{\rm T}]} \right| } \right)\]

where \({\bf n}_i^{[{\rm S}]}\) and \({\bf n}_i^{[{\rm T}]}\) are, respectively, the unit-vector triads of the source and target local coordinate systems.

The final timestep computed by the structural-element portion of the program is

\[\Delta t = \beta_u \Delta t_{\rm stab}\]

where \(\beta_u\) is a user-specified safety factor, and \(\Delta t_{\rm stab}\) is computed using (24), (25) and (27), in which (27) is only applied to active degrees-of-freedom.

Mass Scaling

The structural-element portion of the program may be run in one of three modes (fully dynamic, fully static or partially dynamic):

  • When run in the fully dynamic mode (structure scale-rotational-mass off), mass scaling is not applied; instead, the mass matrix reflects the actual mass of the structural system and is assembled as described in Mass Matrix. The timestep necessary to ensure a stable integration of the equations of motion is determined as described in Computation of Stable Timestep.

  • When run in the fully static mode, the values in the mass matrix are adjusted to ensure that all degrees of freedom have the same frequency response — one that is appropriate to ensure a stable integration of the equations of motion for a fixed value of timestep. The procedure by which this is accomplished is described below.

  • When run in the partially dynamic mode (structure scale-rotational-mass on), the mass matrix reflects the actual mass of the translational degrees-of-freedom, and only the translational degrees-of-freedom are considered when determining the timestep necessary to ensure a stable integration of the equations of motion. The masses of the rotational degrees-of-freedom are adjusted to ensure that all rotational degrees-of-freedom have the same frequency response — one that is appropriate to ensure a stable integration of the equations of motion for the fixed value of timestep determined for the translational degrees-of-freedom. The procedure by which this is accomplished is described below.

Fully Static Mode

For a given value of timestep \(\Delta\tilde t\), stability requires that

\[\Delta t_{\rm stab} = {2\Delta\tilde t \over \beta_u} = {2\over\omega_{\rm max}} F(\xi_{\rm max})\]

where the estimated stable timestep has been multiplied by a safety factor of one-half.

By setting the maximum frequency equal to its bounding value provided by (27), one can write

\[\omega_{\rm max}^2 = \left( {F(\xi_{\rm max}) \beta_u \over \Delta\tilde t} \right)^2 \equiv \psi^2 = \max_i \left( G_i / M_i \right)\]

in which \(\psi\) has been introduced to represent the second term. One can ensure that all degrees of freedom have the same frequency response by setting each \(M_i\) such that

\[M_i = G_i / \psi^2\]

\(\psi\) can be computed directly because \(\xi_{\rm max}\) is given by (25).

Partially Dynamic Mode

  1. Set all translational masses using the procedure described in Mass Matrix.

  2. Determine the timestep necessary to ensure a stable integration of the equations of motion by only considering the translational degrees-of-freedom when using the procedure described in Computation of Stable Timestep.

  3. Determine the highest natural frequency of the discretized system by applying (27) to only the translational degrees-of-freedom, and call this value \(\omega^\star\). Set all masses associated with rotational degrees-of-freedom such that they satisfy

    \[M_i = G_i / {\omega^\star}^2\]

    in order to ensure that the rotational degrees-of-freedom do not produce a larger natural frequency than \(\omega^\star\).

Computational Algorithm

  1. Assume that boundary conditions consisting of externally applied loads, \(\{{\bf R}^{\rm ext}\}\), as well as structural element, link and nodal properties have been specified.

  2. Set initial conditions:

\[\begin{split}\begin{split} n &= 0 \\ \{ {\bf D} \}_0 &= \{ {\bf D}(t = 0) \} \\ \{ \dot{\bf D} \}_{\scriptscriptstyle - 1/2} &= \{ \dot{\bf D}(t = 0) \} \\ [{\bf k}]_{\scriptscriptstyle - 1/2} &= [{\bf k}(t=0)] \quad \hbox{for all elements} \\ \{ {\bf r}^{\rm int} \}_{\scriptscriptstyle - 1} &= \{ {\bf 0} \} \quad\qquad\hbox{for all elements} \end{split}\end{split}\]

  1. Assemble the structure mass matrix \(\lceil {\bf M} \rfloor\) by considering the mode in which the structural-element portion of the program is being run. For the fully dynamic case, employ the procedure described in Mass Matrix. For the fully static or partially dynamic cases, employ the procedures described in Mass Scaling.

  2. Compute stable timestep \(\Delta t\), as described in Computation of Stable Timestep.

  3. Compute the accumulated out-of-balance force vector \(\{ {\bf A} \}_n\):

    1. Set \(\{ {\bf A} \}_n = \{ {\bf R}^{\rm ext} \}_n\).

    2. Assemble the contributions from each structural element, each deformable direction of each link, and each rigidly attached direction of each link, as described in Accumulated Out-of-Balance Force Vector.

    3. If Rayleigh damping is active, assemble the stiffness contribution to the damping force \(\{ {\bf R}^{\beta} \}_n\) using (22) and (23).

  4. Update all nodal velocities \(\{ \dot{\bf D} \}_{n \scriptscriptstyle + 1/2}\):

    1. For all active degrees-of-freedom (i.e., those that are not fixed and are not slaved), update velocities using (16) if local or combined local damping is active, or (20) if Rayleigh damping is active.

    2. For each rigidly attached direction \({\bf n}^{[{\rm S}]}\) of each link, set velocity of source node \(\dot{\bf D}^{[{\rm S}]}\) equal to component of velocity of target entity in the attached direction: \(\quad\dot{\bf D}^{[{\rm S}]} = \left( \dot{\bf D}^{[{\rm T}]} \cdot {\bf n}^{[{\rm S}]} \right) {\bf n}^{[{\rm S}]}\)

  5. Update all nodal positions. Each node stores:

    1. its reference position \({\bf x}^{\rm ref}\), which corresponds to its position at the time that all related stiffness matrices were formed,

    2. its current position \({\bf x}^{\rm cur}\), and

    3. its accumulated rotation \({\boldsymbol\theta}^{\text{cur}}\) since the last update of \({\bf x}^{\text{ref}}\).

    The displacement that has occurred since the last update of \({\bf x}^{\text{ref}}\) is given by

    \[\begin{split}\begin{Bmatrix} {\bf d} \end{Bmatrix} = \begin{Bmatrix} {\bf d}^{\text{tran}} \\ {\bf d}^{\text{rot}} \end{Bmatrix} = \begin{Bmatrix} {\bf x}^{\text{cur}} - {\bf x}^{\text{ref}} \\ {\boldsymbol\theta}^{\text{cur}} \end{Bmatrix}\end{split}\]

    and the current positions and rotations are updated using the central-difference approximation of (12):

    \[ \begin{align}\begin{aligned}\begin{Bmatrix} {\bf x}^{\text{cur}} \end{Bmatrix}_{n+1} = \begin{Bmatrix} {\bf x}^{\text{cur}} \end{Bmatrix}_n + \begin{Bmatrix} {\bf \dot{d}}^{\text{tran}} \end{Bmatrix}_{n+1/2} \Delta t\\\begin{Bmatrix} {\boldsymbol\theta}^{\text{cur}} \end{Bmatrix}_{n+1} = \begin{Bmatrix} {\boldsymbol\theta}^{\text{cur}} \end{Bmatrix}_n + \begin{Bmatrix} {\bf \dot{d}}^{\text{rot}} \end{Bmatrix}_{n+1/2} \Delta t\end{aligned}\end{align} \]

  6. Update all structural element stiffness matrices if a large-strain update is due; otherwise, set

    \[[{\bf k}]_{n \scriptscriptstyle + 1/2} = [{\bf k}]_{n \scriptscriptstyle - 1/2} \; \text{for all elements}\]

  7. Set \(n \leftarrow n + 1\), and go to step 3.

References

AISC. Manual of Steel Construction, Eighth Edition. Chicago: American Institute of Steel Construction Inc. (1980).

Batoz, J. L. “An Explicit Formulation for an Efficient Triangular Plate-Bending Element,” Int. J. Num. Meth. Engng., 18(7), 1077–1089 (1982).

Batoz, J. L., K. J. Bathe, and L. W. Lo. “A Study of Three-Node Triangular Plate Bending Elements,” Int. J. Num. Meth. Engng., 15(12), 1771–1812 (1980).

Belytschko, T. “An Overview of Semidiscretization and Time Integration Procedures,” in Computational Methods for Transient Analysis, Chapter 1, pp. 1–65. T. Belytschko and T. J. R. Hughes, eds. New York: Elsevier Science Publishers, B. V. (1983).

Beneito, C., and Ph. Gotteland. “Three-Dimensional Numerical Modeling of Geosynthetics Mechanical Behavior,” in FLAC and Numerical Modeling in Geomechanics (Proceedings of the Second International FLAC Symposium on Numerical Modeling in Geomechanics, Lyon, France, 29-31 October 2001). D. Billaux et al., eds. Lisse: Balkema (2001).

Carpenter, N., H. Stolarski, and T. Belytschko. “Improvements in 3-Node Triangular Shell Elements,” Int. J. Num. Meth. Engng., 23(9), 1643–1667 (1986).

Cernica, J. N. Geotechnical Engineering: Foundation Design, New York: John Wiley & Sons Inc. (1995).

Clayton, C. R. I., J. Milititsky, and R. T. Woods. Earth Pressure and Earth-Retaining Structures, pp. 145–148. London: Blackie Academic & Professional (1993).

Cook, R. D. “A Plane Hybrid Element with Rotational D. O. F. and Adjustable Stiffness,” Int. J. Num. Meth. Engng., 24(8), 1499–1508 (1987).

Cook, R. D., D. S. Malkus, and M. E. Plesha. Concepts and Applications of Finite Element Analysis, Third Edition. New York: John Wiley & Sons Inc. (1989).

Cook, R. D., D. S. Malkus, M. E. Plesha, and R. J. Witt. Concepts and Applications of Finite Element Analysis, Fourth Edition. New York: John Wiley & Sons Inc. (2002).

Isaacson, E., and H. B. Keller. Analysis of Numerical Methods, pp. 135–136. New York: John Wiley & Sons Inc. (1963).

Massonnet, C. E. Résistance Des Matériaux. Sciences Et Lettres, Liége (1960).

McGuire, W., and R. H. Gallagher. Matrix Structural Analysis. New York: John Wiley & Sons Inc. (1979).

Potyondy, D. “Nonlinear Structural Elements,” Itasca Consulting Group, Inc., Technical Memorandum 5-8121:22TM36 (September 9, 2022), Minneapolis, Minnesota (2022).

St. John, C. M., and D. E. Van Dillen. “Rockbolts: A New Numerical Representation and Its Application in Tunnel Design,” in Rock Mechanics — Theory - Experiment - Practice (Proceedings of the 24th U.S. Symposium on Rock Mechanics, Texas A&M University, June 1983), pp. 13–26. New York: Association of Engineering Geologists (1983).

Ugural, A. C. Stresses in Plates and Shells, pp. 1–16. New York: McGraw-Hill Publishing Company Inc. (1981).