Element Stiffness

To generate a stiffness matrix for a curvilinear quadrilateral or triangular element a new approach must be used. Most finite element codes used an approach based on isoparametric or similar elements.

In an isoparametric element the element displacements are interpolated in the same way as the geometry, e.g. a plane stress element. In a superparametric or degenerate isoparametric element the interpolation on the geometry is of a higher order than the interpolation of the displacements, e.g. a plate element. In a subparametric element the interpolation of the geometry is of a lower order than the interpolation of the displacements, e.g. an eight noded straight sided quadrilateral element, where a different geometric interpolation function is used for the geometry from that for the displacements. The term isoparametric is often used as a general term to cover all these element types.

For a plane stress problem, we can establish a material matrix $\mathbf{C}$ that relates stress and strain. The displacements in a local coordinate system are

$$\mathbf{u} = \left\{ u,v \right\}$$

the strains are

$$\mathbf{\varepsilon} = \left\{ \varepsilon_{xx},\varepsilon_{yy},\gamma_{xy} \right\}$$

and the stresses are

$$\mathbf{\sigma} = \left\{ \sigma_{xx},\sigma_{yy},\tau_{xy} \right\}$$

For an elastic-isotropic material the material matrix is

$$\mathbf{C}_{p} = \frac{E}{1 - \nu^{2}}\begin{bmatrix} 1 & \nu & \\ \nu & 1 & \\ & & \frac{1 - \nu}{2} \\ \end{bmatrix}$$

where $E,\nu$ are the Young’s modulus and Poisson’s ratio respectively.

Note that there is an out of plane strain $\varepsilon_{zz}$, which we can ignore as it plays no part in the element formulation.

The strains are defined in terms of the displacements as

$$\varepsilon_{xx} = \frac{\partial u}{\partial x}\qquad\varepsilon_{yy} = \frac{\partial v}{\partial y}\qquad\gamma_{xy} = \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}$$

The simplest elements to consider are the 4 noded and 8 noded quadrilateral elements, of which the 4 noded element can be considered a simplification of the 8 noded element. A typical 8 noded element is shown below. The element has an arbitrary local coordinate system based on the location of the nodes and the element property axis $(x, y)$, and a natural (curvilinear) coordinate system $(r, s)$ based on the topology of the nodes.

We can set up interpolation functions to interpolate the geometry as follows

$$\mathbf{x}_{s}(r,s) = \sum_{i}^{}{h_{i}\mathbf{x}_{m,i}}$$

where the $h_{i}$ are the interpolation functions defined below.

These interpolation functions are chosen so that

$$h_{i} = \begin{Bmatrix} 1 \\ 0 \\ \end{Bmatrix} \text{when node is} \begin{Bmatrix} i \\ \neq i \\ \end{Bmatrix}$$

As the elements are isoparametric we use the same interpolation function for the displacements so the displacement in the element is related to the nodal displacements by

$$u_{s}(r,s) = \sum_{i}^{}{h_{i}u_{m,i}}$$

To evaluate the stiffness matrix we need the strain-displacement transformation matrix. The element displacements are obtained in terms of derivatives of the element displacements with respect to the local coordinate system $(x,y)$. Because the elements displacements are in the natural coordinate system $(r,s)$ we need to relate the derivatives in the local coordinate system to those in the natural coordinate system. We can write an equation for the derivative with respect to $x$ in terms of derivatives with respect to $(r,s)$

$$\frac{\partial}{\partial x} = \frac{\partial}{\partial r}\frac{\partial r}{\partial x} + \frac{\partial}{\partial s}\frac{\partial s}{\partial x}$$

to establish these derivatives we use the chain rule to set up the following relationship

$$\begin{bmatrix} \frac{\partial}{\partial r} \\ \frac{\partial}{\partial s} \\ \end{bmatrix} = \begin{bmatrix} \frac{\partial x}{\partial r} & \frac{\partial y}{\partial r} \\ \frac{\partial x}{\partial s} & \frac{\partial y}{\partial s} \\ \end{bmatrix}\begin{bmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ \end{bmatrix}$$

or in matrix notation

$$\frac{\partial}{\partial\mathbf{r}} = \mathbf{J}\frac{\partial}{\partial\mathbf{x}}$$

where $\mathbf{J}$ is the Jacobian operator relating natural coordinate derivatives to local coordinate derivatives. Given that we know $x$ and $y$ in terms of the interpolation functions the Jacobian operator is easily found

$$\frac{\partial}{\partial\mathbf{x}} = \mathbf{J}^{- 1}\frac{\partial}{\partial\mathbf{r}}$$

This requires that the inverse of the Jacobian exist, which requires that there is a one to one correspondence between natural and local coordinates. This will be the case provided the element is not grossly distorted from a square and that it does not fold back on itself.

We can then evaluate

$$\frac{\partial u}{\partial x},\frac{\partial u}{\partial y},\frac{\partial v}{\partial x},\frac{\partial v}{\partial y}$$

and thus we can construct the strain-displacement transformation matrix, $\mathbf{B}$

$$\boldsymbol{\varepsilon} = \mathbf{Bu}$$

where $\mathbf{u}$ is the vector of nodal displacements. The element stiffness corresponding to the local element degrees of freedom is then

$$\mathbf{K} = \int_{V}^{}{\mathbf{B}^{T}\mathbf{CB}dV}$$

The elements of $\mathbf{B}$ are functions of the natural coordinate system, $(r,s)$. Therefore the volume integration extends over the natural coordinate volume so the volume differential needs to be written in terms of the natural coordinates

$$dV = \det\mathbf{J}drds$$

The volume integral is not normally amenable to an explicit integration so normally a numerical integration technique is used. The integral can be written

$$\mathbf{K} = \iint_{r,s}^{}{\mathbf{F}drds}$$

where

$$\mathbf{F} = \mathbf{B}^{T}\mathbf{CB}\det\mathbf{J}$$

and the integral is performed in the natural coordinate system of the element. This is convenient as the limits of the integration are then ±1. The stiffness can then be calculated

$$\mathbf{K} = \sum_{}^{}{\alpha_{ij}\mathbf{F}_{ij}}$$

where the matrix $\mathbf{F}$ is evaluated at the Gaussian integration points $\left( r_{i},s_{i} \right)$ and $\alpha_{ij}$ are Gaussian weights.

In a similar way the mass matrix and the load vectors are established.

$$\mathbf{M} = \int_{V}^{}{\rho\mathbf{H}^{T}\mathbf{H}dV}$$

$$\mathbf{r}_{b} = \int_{V}^{}{\mathbf{H}^{T}\mathbf{f}_{b}dV}$$

$$\mathbf{r}_{s} = \int_{S}^{}{\mathbf{H}_s^{T}\mathbf{f}_{s}dS}$$

$$\mathbf{r}_{i} = \int_{V}^{}{\mathbf{B}^{T}\mathbf{\tau}_{i}dV}$$

where $\mathbf{H}$ is a matrix of interpolation functions and $\mathbf{r}_{b},\mathbf{r}_{s},\mathbf{r}_{i}$, are the body force vector, surface force vector and initial stress vector respectively.

Geometric stiffness matrix of shell element

The geometric stiffness matrix is calculated from:

$$\begin{aligned}\mathbf{K}_g &= \int_A \mathbf{G}^T\mathbf{NG}dA\\ &=\iint_{r,s}\mathbf{G}^T\mathbf{NG}\det\mathbf{J}drds\end{aligned}$$

where

$$\mathbf{N} = \begin{bmatrix} N_{x} & N_{xy} \\ N_{yx} & N_{y} \\ \end{bmatrix}$$

With $\left( N_{x},N_{y},N_{xy} \right)$ are the in-plane forces of the shell element in $x$, $y$ and $xy$ shear directions respectively.

$$\mathbf{G} = \begin{bmatrix} G_{1} & G_{2} & \ldots & G_{n} \\ \end{bmatrix}$$

and $n$ is the number of nodes of the elements

$$\mathbf{G}_{i} = \begin{bmatrix} 0 & 0 & \frac{\partial h_{i}}{\partial x} & 0 & 0 & 0 \\ 0 & 0 & \frac{\partial h_{i}}{\partial y} & 0 & 0 & 0 \\ \end{bmatrix}$$

2D Plate Elements

The formulation of 2D plate and shell elements considers both in-plane and transverse (out-of-plane) displacements. Following Mindlin-Reissner plate theory, in addition to the bending strains we consider the effect of transverse shear strain in our complete expression for the element strain

$$\begin{Bmatrix} \gamma_{xz} \\ \gamma_{yz} \\ \end{Bmatrix} = \begin{Bmatrix} \frac{\partial w}{\partial x} - \beta_{x} \\ \frac{\partial w}{\partial y} - \beta_{y} \\ \end{Bmatrix}$$

where $w$ is the out-of-plane displacement and $\beta$ is introduced as an independent variable to express the section rotations (i.e. rotations of the transverse normals about the local $x$ and $y$ axes).

We can define separate material matrices $\mathbf{C}_{b}$ and $\mathbf{C}_{s}$ that relate stress and strain for the pure bending and shear strains respectively and so the pure bending moments and shear forces can be written

$$\begin{Bmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \tau_{xy} \\ \end{Bmatrix} = - z\, \mathbf{C}_{b}\begin{Bmatrix} \frac{\partial\beta_{x}}{\partial x} \\ \frac{\partial\beta_{y}}{\partial y} \\ \frac{\partial\beta_{x}}{\partial y} + \frac{\partial\beta_{y}}{\partial x} \\ \end{Bmatrix}$$

and

$$\begin{Bmatrix} \tau_{xz} \\ \tau_{yz} \\ \end{Bmatrix} = \mathbf{C}_{s}\begin{Bmatrix} \frac{\partial w}{\partial x} - \beta_{x} \\ \frac{\partial w}{\partial y} - \beta_{y} \\ \end{Bmatrix}$$

respectively.

In this way we can then express the local element stiffness matrix as a summation of the in-plane and out-of-plane stiffness’s

$$\mathbf{K} = \int_{V}^{}{\left( {\mathbf{B}_{\kappa}}^{T}\mathbf{C}_{b}\mathbf{B}_{\kappa} + {\mathbf{B}_{\gamma}}^{T}\mathbf{C}_{s}\mathbf{B}_{\gamma} \right)dV}$$

where $\mathbf{B}_{\kappa}$ and $\mathbf{B}_{\gamma}$ represent the strain-displacement transforms for the bending and shear components respectively.

While brief, this outlines the basic approach to the Mindlin-Reissner 2D element stiffness formulation. In GSA we label this concisely as the Mindlin formulation.

MITC Element Formulation

We find that the Mindlin formulation is an effective approach for 2D parabolic elements where the 8-noded element accommodates sufficient terms in the stiffness matrix to sufficiently define the behaviour of the element numerically. However the same formulation defined over a linear element becomes noticeably more problematic where the absence of available terms in our 4-noded stiffness matrix leads to numerical difficulties in expressing the same element behaviour. Specifically we find a difficulty in attempting to represent the transverse shear strain terms. Analytically as before we can express the transverse shear strain as

$$\begin{Bmatrix} \gamma_{xz} \\ \gamma_{yz} \\ \end{Bmatrix} = \begin{Bmatrix} \frac{\partial w}{\partial x} - \beta_{x} \\ \frac{\partial w}{\partial y} - \beta_{y} \\ \end{Bmatrix}$$

although numerically, the difference in the order of terms for the shear strain may lead to artificial stiffening of the element where the shear terms are numerically constrained from approaching zero. See reference 1 in the bibliography for further information. This restriction would be particularly noticeable where the thickness of the plate is small.

A widely practiced remedy is to under-integrate the shear term and while effective, its use is at the cost of reduced accuracy and stability for the element. The problem of stability alone is often of greatest concern where the phenomenon of hourglassing can become apparent in elements where the thickness to length ratio is large.

An alternative formulation was put forward by Bathe and Dvorkin and has been found to be especially effective at resolving these difficulties. The formulation is extendable to higher order elements although we find the approach is most effective when resolving the difficulties most apparent in linear elements. The formulation is based upon the theory of Mixed Interpolated Tensoral Components (MITC). For the pure displacement-based Mindlin formulation we use the independent variables

$$\begin{aligned}w &= \sum_{i = 1}^{n}{{\mathbf{h}_{i}}^{T}\mathbf{w}_{i}}\\ \beta &= \sum_{i = 1}^{n}{{\mathbf{h}_{i}}^{T}\theta_{i}}\end{aligned}$$

where Bathe now introduces a separate independent variable to represent the transverse shear term

$$\gamma = \sum_{i = 1}^{n}{{{\mathbf{h}}_{i}^{*}}^T\gamma_{Pi}}$$

We use ${\mathbf{h}_{i}}^{*}$ to represent an additional set of interpolation functions for our new variable $\gamma$ which we find by a direct evaluation of the shear strain at points $Pi$, that is

$$\gamma_{Pi} = \left( \frac{\partial w}{\partial x} - \beta \right)_{Pi}$$

For a linear 2D element we obtain direct values for the shear strain at four points A, B, C, D on the element and so we evaluate the displacement and section rotations at these points through direct interpolation.

We can then construct the interpolations

$$\begin{aligned}\gamma_{rz} &= \frac{1}{2}(1 + s)\gamma_{rz}^{A} + \frac{1}{2}(1 - s)\gamma_{rz}^{C}\\ \gamma_{sz} &= \frac{1}{2}(1 + r)\gamma_{sz}^{D} + \frac{1}{2}(1 - r)\gamma_{sz}^{B}\end{aligned}$$

using direct values for the shear strains obtained at the four points. This replaces our original expression for the shear terms and we continue to construct the local stiffness matrix as normal in a similar approach as before in 2D isoparametric elements. It is lastly worthwhile to note that the interpolation above assumes our element is in the isoparametric coordinate system. Further transformations are necessary to interpolate the shear strains directly through an arbitrary element in local space. Indeed, when this is done the element shows considerably improved predicative capabilities for distorted elements.

Element shape functions

The element shape functions for 2D elements are

$$u(r,s) = \sum_{i}^{}{h_{i}u_{i}}$$

These are for the different element shapes

Tri-3

$$\begin{aligned} h_{1} & = 1 - r - s\\ h_{2} & = r\\ h_{3} & = s \end{aligned}$$

Quad-4

$$\begin{aligned} h_{1} & = \frac{1}{4}(1 - r)(1 - s)\\ h_{2} & = \frac{1}{4}(1 + r)(1 - s)\\ h_{3} & = \frac{1}{4}(1 + r)(1 + s)\\ h_{4} & = \frac{1}{4}(1 - r)(1 + s) \end{aligned}$$

Tri-6

$$\begin{aligned} h_{1} & = 1 - r - s - \frac{1}{2}h_{4} - \frac{1}{2}h_{6}\\ h_{2} & = r - \frac{1}{2}h_{4} - \frac{1}{2}h_{5}\\ h_{3} & = s - \frac{1}{2}h_{5} - \frac{1}{2}h_{6}\\ h_{1} & = 4r(1 - r - s)\\ h_{2} & = 4rs\\ h_{3} & = 4s(1 - r - s) \end{aligned}$$

Quad-8

$$\begin{aligned} h_{1} & = \frac{1}{4}(1 - r)(1 - s) - \frac{1}{2}h_{5} - \frac{1}{2}h_{8}\\ h_{2} & = \frac{1}{4}(1 + r)(1 - s) - \frac{1}{2}h_{5} - \frac{1}{2}h_{6}\\ h_{3} & = \frac{1}{4}(1 + r)(1 + s) - \frac{1}{2}h_{6} - \frac{1}{2}h_{7}\\ h_{4} & = \frac{1}{4}(1 - r)(1 + s) - \frac{1}{2}h_{7} - \frac{1}{2}h_{8}\\ h_{5} & = \frac{1}{4}\left( 1 - r^{2} \right)(1 - s)\\ h_{6} & = \frac{1}{4}(1 + r)\left( 1 - s^{2} \right)\\ h_{7} & = \frac{1}{4}\left( 1 - r^{2} \right)(1 + s)\\ h_{8} & = \frac{1}{4}(1 - r)\left( 1 - s^{2} \right) \end{aligned}$$

2D Element Shape Checks

A number of element checks are carried out by GSA prior to a GSS analysis. Other analysis programs may have different limits but the same principles apply. For GSS the following warnings, severe warnings and errors are produced.

Triangle

Warning	Severe warning	Error
$5 < r_{\max} \leq 15$	$r_{\max} > 15$
$15 \leq \theta_{\min} < 30$	$\theta_{\min} < 15$
$150 \leq \theta_{\max}\leq 165$		$\theta_{\max} > 165$

Quad

Warning	Severe warning	Error
$5 < r_{\max}\leq 15$	$r_{\max} > 15$
$25 \leq \theta_{\min} < 45$	$\theta_{\min} < 25$
$135 \leq \theta_{\max} \leq 155$		$\theta_{\max} > 155$
$0.00001 < h_{\max} \leq 0.01$	$h_{\max} > 0.01$

where

$r_{\max}$ longest side / shortest side

$\theta_{\min}$ minimum angle

$\theta_{\max}$ maximum angle

$h_{\max}$ distance out of the plane of the element of edge 1 / longest side

Notes:

The distance out of plane of edge 1 is calculated as $h_{\max}$

$${h_{\max} = \frac{\mathbf{n} \cdot \left( \mathbf{c}_{2} - \mathbf{c}_{1} \right)}{s_{\max}}}$$

Where $\mathbf{n}$ is the element normal, $\mathbf{c}_{1}$, $\mathbf{c}_{2}$ are the coordinates of the first and second corner nodes and $s_{\max}$ is the length of the longest side of the element.

Mid-side node locations not checked but should be approximately halfway along edge.

No check on ratio thickness/shortest side.

Hourglassing

When Quad 4 elements with the Mindlin formulation are used in bending it is possible to encounter hourglassing problems. This is a problem that arises with under-integrated elements where there are insufficient stiffness terms to fully represent the stiffness of the element. The problem is noticeable in the results by an hourglass pattern in the mesh as shown below.

This problem is avoided using the MITC formulation for Quad 4 elements. This formulation uses a separate interpolation method for the transverse shear strains and provides considerably greater stability than the original Mindlin method. The original Mindlin method is kept in GSA for compatibility with previous models although for new models the MITC formulation is recommended.

Alternatively when parabolic accuracy is required Quad 8 elements are recommended. These also formulate elements that are still stiff in all modes of deformation, even when under-integrated.

3D Element

The element shape functions for 2D elements are

$$u(r,s,t) = \sum_{i}^{}{h_{i}u_{i}}$$

These are for the different element shapes

Tetra-4

$$\begin{aligned} h_{1} & = \frac{1}{8}(1 - r)(1 - s)(1 - t)\\ h_{2} & = \frac{1}{8}(1 + r)(1 - s)(1 - t)\\ h_{3} & = \frac{1}{4}(1 + s)(1 - t)\\ h_{4} & = \frac{1}{2}(1 + t) \end{aligned}$$

Pyramid-5

$$\begin{aligned} h_{1} & = \frac{1}{8}(1 - r)(1 - s)(1 - t)\\ h_{2} & = \frac{1}{8}(1 + r)(1 - s)(1 - t)\\ h_{3} & = \frac{1}{8}(1 + r)(1 + s)(1 - t)\\ h_{4} & = \frac{1}{8}(1 - r)(1 + s)(1 - t)\\ h_{5} & = \frac{1}{2}(1 + t) \end{aligned}$$

Wedge-6

$$\begin{aligned} h_{1} & = \frac{1}{8}(1 - r)(1 - s)(1 - t)\\ h_{2} & = \frac{1}{8}(1 + r)(1 - s)(1 - t)\\ h_{3} & = \frac{1}{4}(1 + s)(1 - t)\\ h_{4} & = \frac{1}{8}(1 - r)(1 - s)(1 + t)\\ h_{5} & = \frac{1}{8}(1 + r)(1 - s)(1 + t)\\ h_{6} & = \frac{1}{4}(1 + s)(1 + t) \end{aligned}$$

Brick-8

$$\begin{aligned} h_{1} & = \frac{1}{8}(1 - r)(1 - s)(1 - t)\\ h_{2} & = \frac{1}{8}(1 + r)(1 - s)(1 - t)\\ h_{3} & = \frac{1}{8}(1 + r)(1 + s)(1 - t)\\ h_{4} & = \frac{1}{8}(1 - r)(1 + s)(1 - t)\\ h_{5} & = \frac{1}{8}(1 - r)(1 - s)(1 + t)\\ h_{6} & = \frac{1}{8}(1 + r)(1 - s)(1 + t)\\ h_{7} & = \frac{1}{8}(1 + r)(1 + s)(1 + t)\\ h_{8} & = \frac{1}{8}(1 - r)(1 + s)(1 + t) \end{aligned}$$

The stiffness matrix for brick elements is obtained using selective reduced integration to alleviate volumetric locking. The $\bar{B}$ approach has been implemented in GSA where the strain-displacement matrix is split into dilatational and deviatoric parts$^7$ and then the dilatational part of matrix is replaced with improved dilatational component. The strain-displacement matrix for brick elements is given by

$$B = \begin{bmatrix} B_{1} & 0 & 0 \\ 0 & B_{2} & 0 \\ 0 & 0 & B_{3} \\ 0 & B_{3} & B_{2} \\ B_{3} & 0 & B_{1} \\ B_{2} & B_{1} & 0 \\ \end{bmatrix}$$

The dilatational and deviatoric parts of the matrix can be computed using

$$\begin{aligned}B^{dil} &= \begin{bmatrix} B_{1} & B_{2} & B_{3} \\ B_{1} & B_{2} & B_{3} \\ B_{1} & B_{2} & B_{3} \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{bmatrix} \\ B^{dev} &= B - B^{dil} \end{aligned}$$

$\bar{B}$ is computed using

$$\begin{aligned}\bar{B} &= B^{dev} + {\bar{B}}^{dil}\\ {\bar{B}}^{dil} &= \frac{1}{3}\begin{bmatrix} {\bar{B}}_{1} & {\bar{B}}_{2} & {\bar{B}}_{3} \\ {\bar{B}}_{1} & {\bar{B}}_{2} & {\bar{B}}_{3} \\ {\bar{B}}_{1} & {\bar{B}}_{2} & {\bar{B}}_{3} \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{bmatrix}\end{aligned}$$

and

$${\bar{B}}_{i} = \frac{\int_{}^{}B_{i}d\Omega}{\int_{}^{}{d\Omega}}\qquad i = 1,2\text{ and }3$$

for $i = 1,2$ and $3$, and where $\int d\Omega$ is the volume of the element.

Link Element

Link elements are different to the other element types in that they apply a constraint between a pair of nodes. The primary node is the first node specified in the element topology and this is the node that is retained in the solution. The other node (the constrained node) is related back to the primary node.

$$\mathbf{u}_{s} = \mathbf{T}\mathbf{u}_{m}$$

The degrees of freedom at the constrained node that are linked will depend on the type of link. The link allows the constrained node to be fixed (where the rotations at the constrained node depend on the rotation of the primary) or pinned (where the rotations at the constrained node are independent of the rotation of the primary). The primary node is always has the rotations linked to the rest of the structure. The links can act in all directions or be restricted to act in a plane ($xy$, $yz$ or $zx$) where the nodes are rigidly connected for motion in the plane but are independent for out of plane motion.

The constraint conditions are summarised below:

All directions

$$\begin{bmatrix} 1 & 0 & 0 & 0 & z & - y \\ 0 & 1 & 0 & - z & 0 & x \\ 0 & 0 & 1 & y & - x & 0 \\ 0 & 0 & 0 & \delta & 0 & 0 \\ 0 & 0 & 0 & 0 & \delta & 0 \\ 0 & 0 & 0 & 0 & 0 & \delta \\ \end{bmatrix} \text{ where } \delta = 1/0 \text{ if fixed / pinned}$$

Plane ($xy$ plane as example)

$$\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & - y \\ 0 & 1 & 0 & 0 & 0 & x \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \delta \\ \end{bmatrix} \text{ where }\delta = 1/0 \text{ if fixed / pinned}$$

Loads applied to link elements will be correctly transferred to the primary degree of freedom as a force + moment so no spurious moments result.

The inertia properties of a link element can be calculated from the masses at the nodes attached to the rigid element as follows. The mass is

$$m = \sum_{i}^{}m_{i}$$

and the coordinates of the centre of mass are

$$x_{c} = \frac{\sum_{i}^{}{m_{i}x_{i}}}{m}\quad y_{c} = \frac{\sum_{i}^{}{m_{i}y_{i}}}{m}\quad z_{c} = \frac{\sum_{i}^{}{m_{i}z_{i}}}{m}$$

The inertia about the global origin is then

$$\begin{aligned}{\bar{I}}_{xx} = \sum_{i}^{}{m_{i}\left( {y_{i}}^{2} + {z_{i}}^{2} \right)}\quad{\bar{I}}_{xy} = - \sum_{i}^{}{m_{i}x_{i}y_{i}}\\ {\bar{I}}_{yy} = \sum_{i}^{}{m_{i}\left( {z_{i}}^{2} + {x_{i}}^{2} \right)}\quad{\bar{I}}_{yz} = - \sum_{i}^{}{m_{i}y_{i}z_{i}}\\ {\bar{I}}_{zz} = \sum_{i}^{}{m_{i}\left( {x_{i}}^{2} + {y_{i}}^{2} \right)}\quad{\bar{I}}_{zx} = - \sum_{i}^{}{m_{i}z_{i}x_{i}}\end{aligned}$$

and relative to the centre of mass $\left( x_{c},y_{c},z_{c} \right)$ this is

$$\begin{aligned}I_{xx} = \sum_{i}^{}{m_{i}\left( \left( y_{i} - y_{c} \right)^{2} + \left( z_{i} - z_{c} \right)^{2} \right)}\quad{\bar{I}}_{xy} = - \sum_{i}^{}{m_{i}\left( x_{i} - x_{c} \right)\left( y_{i} - y_{c} \right)}\\ I_{yy} = \sum_{i}^{}{m_{i}\left( \left( z_{i} - z_{c} \right)^{2} + \left( x_{i} - x_{c} \right)^{2} \right)}\quad{\bar{I}}_{yz} = - \sum_{i}^{}{m_{i}\left( y_{i} - y_{c} \right)\left( z_{i} - z_{c} \right)}\\ I_{zz} = \sum_{i}^{}{m_{i}\left( \left( x_{i} - x_{c} \right)^{2} + \left( y_{i} - y_{c} \right)^{2} \right)}\quad{\bar{I}}_{zx} = - \sum_{i}^{}{m_{i}\left( z_{i} - z_{c} \right)\left( x_{i} - x_{c} \right)}\end{aligned}$$

The constraint equations for a link element assume small displacements. When large displacements are applied to a link element the constraint equations no longer apply and the links between constrained and primary get stretched. This effect can be noticeable in a dynamic analysis where the results are scaled to an artificially large value. When these are scaled to realistic value this error should be insignificant.

---

$^7$ Hughes T. J. R. The Finite Element Method – Linear Static and Dynamic Finite Element Analysis, Dover, 2000