Skip to content

Advanced Element Formulations

The standard displacement-based formulation presented in Formulation of 3D Solid Elements exhibits excessive artificial stiffness known as locking (volumetric locking or shear locking) when applied to nearly incompressible materials or bending-dominated thin-walled structures. To avoid this, FrontISTR provides the B-bar and F-bar methods, which replace only the volumetric part of the B matrix or deformation gradient; incompatible elements with internal degrees of freedom; u-p mixed elements that treat pressure as an independent unknown field; and MITC shell and beam elements specialized for plate and beam structures.

This chapter summarizes the formulations of these advanced elements and structural elements on an element-by-element basis.

B-bar Method

When an 8-node linear hexahedral element is used for a nearly incompressible material, the strain field within one element can conflict with the constant-volume constraint, producing excessive stiffness known as volumetric locking. The B-bar method alleviates this overconstraint by replacing the components of the B matrix that contribute to volumetric expansion with values evaluated at the element center [Hughes1980].

At the element center \(\boldsymbol{r} = \boldsymbol{0}\), let the B matrix computed from the spatial derivatives of the shape functions be \(\bar{\boldsymbol{B}}\); at integration point \(\boldsymbol{r}\), let the ordinary B matrix be \(\boldsymbol{B}(\boldsymbol{r})\). In the displacement-strain relation for node \(\alpha\) and degree of freedom \(i\), the volumetric strain components \((\varepsilon_{11}, \varepsilon_{22}, \varepsilon_{33})\) are modified by adding

\[ \Delta B_{i\alpha} = \tfrac{1}{3}\bigl(\bar{B}_{i\alpha}(\boldsymbol{0}) - B_{i\alpha}(\boldsymbol{r})\bigr) \]

For the shear components \((\varepsilon_{12}, \varepsilon_{23}, \varepsilon_{31})\), the ordinary \(\boldsymbol{B}\) is used. The element stiffness and internal-force vector are assembled using the resulting B-bar matrix.

FrontISTR provides this formulation specifically for the 8-node linear hexahedral element (element ID 361; see Element Numbering System), and it can be applied to small-deformation, Total Lagrange, and Updated Lagrange formulations.

F-bar Method

Under finite deformation, volumetric change enters nonlinearly through the deformation gradient \(\boldsymbol{F}\). The F-bar method [deSouzaNeto1996] applies a volumetric-locking treatment analogous to the B-bar method at the deformation-gradient level.

Let the volume ratio of the deformation gradient evaluated at the element center \(\boldsymbol{r} = \boldsymbol{0}\) be \(J_0 = \det \boldsymbol{F}(\boldsymbol{0})\), and let the volume ratio at an integration point be \(J = \det \boldsymbol{F}(\boldsymbol{r})\). The deformation gradient at the integration point is replaced by

\[ \bar{\boldsymbol{F}} = \left(\frac{J_0}{J}\right)^{1/3} \boldsymbol{F} \]

This gives \(\det \bar{\boldsymbol{F}} = J_0\), so the volume ratio throughout the element is matched to the value at the element center. The replaced \(\bar{\boldsymbol{F}}\) is used for stress evaluation and construction of the strain-displacement matrix, and the element tangent stiffness includes the additional terms associated with this replacement.

FrontISTR implements the F-bar method specifically for the 8-node linear hexahedral element, and it can be applied to small-deformation as well as Total Lagrange / Updated Lagrange nonlinear geometry.

Incompatible Elements

The 8-node linear hexahedral element lacks the strain modes required for bending and exhibits bending locking in bending-dominated problems. Incompatible elements [Taylor1976] introduce additional displacement modes inside the element to compensate for this deficiency.

In addition to the element nodal displacements \(\boldsymbol{u}^e\), incompatible-mode degrees of freedom \(\boldsymbol{\alpha} \in \mathbb{R}^{9}\) that exist only inside the element are introduced as 3 directions × 3 modes per element, and the displacement field is approximated as

\[ \boldsymbol{u}(\boldsymbol{r}) = \sum_{\alpha=1}^{8} N_\alpha^e(\boldsymbol{r})\, \boldsymbol{u}^e_\alpha + \sum_{k=1}^{3} M_k(\boldsymbol{r})\, \boldsymbol{\alpha}_k \]

For natural coordinates \(\boldsymbol{r} = (\xi, \eta, \zeta)\), the incompatible shape functions are taken as \(M_1 = 1 - \xi^2\), \(M_2 = 1 - \eta^2\), and \(M_3 = 1 - \zeta^2\). They do not guarantee continuity across element boundaries, but they add an internal space capable of reproducing the strains associated with bending modes.

The element stiffness is first assembled in an external-DOF/internal-DOF block form

\[ \begin{bmatrix} \boldsymbol{K}_{dd} & \boldsymbol{K}_{d\alpha} \\ \boldsymbol{K}_{\alpha d} & \boldsymbol{K}_{\alpha\alpha} \end{bmatrix} \begin{bmatrix} d\boldsymbol{u}^e \\ d\boldsymbol{\alpha} \end{bmatrix} = \begin{bmatrix} \boldsymbol{F}^e_{\text{ext}} \\ \boldsymbol{0} \end{bmatrix} \]

After assembly in this form, static condensation eliminates the internal degrees of freedom using \(d\boldsymbol{\alpha} = -\boldsymbol{K}_{\alpha\alpha}^{-1}\boldsymbol{K}_{\alpha d}\,d\boldsymbol{u}^e\), yielding the element stiffness involving only the external degrees of freedom

\[ \boldsymbol{K}^e = \boldsymbol{K}_{dd} - \boldsymbol{K}_{d\alpha}\,\boldsymbol{K}_{\alpha\alpha}^{-1}\,\boldsymbol{K}_{\alpha d} \]

which is passed to the global assembly.

FrontISTR implements the incompatible element specifically for the 8-node linear hexahedral element (C3D8IC), and it can be applied to small-deformation, Total Lagrange, and Updated Lagrange formulations.

U-P Mixed Elements

Whereas the B-bar and F-bar methods correct the volumetric component within a displacement-based framework, the u-p mixed (U-P) element uses a mixed formulation [Bathe1996] that introduces pressure \(\lambda\) as an unknown field independent of displacement. For nearly incompressible materials, such as rubber-like materials with Poisson's ratio extremely close to 0.5 or metals after plastic deformation, enforcing the constant-volume constraint using the displacement field alone causes volumetric locking; treating pressure as an independent variable relaxes this constraint.

The stress is separated into deviatoric and pressure components as

\[ \boldsymbol{\sigma} = \boldsymbol{\sigma}_{\mathrm{dev}} + \lambda\,\boldsymbol{I}, \qquad \boldsymbol{\sigma}_{\mathrm{dev}} = \mathbf{D}_{\mathrm{dev}}\,\boldsymbol{\varepsilon} \]

Here, \(\mathbf{D}_{\mathrm{dev}}\) is the deviatoric elastic matrix obtained by removing from the elastic matrix the volumetric part proportional to the bulk modulus \(K\). The pressure \(\lambda\) and volumetric strain \(g = \mathrm{tr}\,\boldsymbol{\varepsilon}\) are related by the constraint through the compressibility \(\alpha^{-1} = 1/K\)

\[ g - \alpha^{-1}\lambda = 0 \]

Discretization with displacement \(\boldsymbol{u}\) and pressure \(\lambda\) as unknowns yields the element coupled system

\[ \begin{bmatrix} \mathbf{K}_{uu} & \mathbf{K}_{up} \\ \mathbf{K}_{up}^{T} & \mathbf{K}_{pp} \end{bmatrix} \begin{bmatrix} d\boldsymbol{u} \\ d\lambda \end{bmatrix} = \begin{bmatrix} \boldsymbol{f}_{u} \\ \boldsymbol{f}_{p} \end{bmatrix} \]

where \(\mathbf{K}_{uu}\) contains the deviatoric elastic contribution and, for finite deformation, geometric stiffness; \(\mathbf{K}_{up}\) couples volumetric strain and pressure; and \(\mathbf{K}_{pp} = -\int \alpha^{-1}\,\boldsymbol{N}_p \boldsymbol{N}_p^{T}\,dV\) is the pressure stabilization term (\(\boldsymbol{N}_p\) is the pressure shape function). Because the pressure degrees of freedom are internal to the element,

\[ \mathbf{K}_{\mathrm{eff}} = \mathbf{K}_{uu} - \mathbf{K}_{up}\,\mathbf{K}_{pp}^{-1}\,\mathbf{K}_{up}^{T} \]

static condensation is performed in this form, and the effective stiffness involving only the external (displacement) degrees of freedom is passed to the global assembly.

FrontISTR implements the U-P element specifically for the 8-node linear hexahedral element, with one pressure degree of freedom per element (constant within the element). It can be applied to small-deformation, Total Lagrange, and Updated Lagrange formulations. In Updated Lagrange, the deviatoric stress is updated using an objective stress rate (Jaumann/Hughes-Winget type), and the pressure \(\lambda\,\boldsymbol{I}\) is then imposed using the value obtained by static condensation.

Shell Elements

Shell elements based on Reissner-Mindlin plate/shell theory are used for thin-walled plate and shell structures. As the thickness decreases, low-order displacement-based plate/shell elements overestimate transverse shear strain through a mechanism analogous to volumetric locking (shear locking), causing the stiffness for bending modes to increase excessively. The MITC (Mixed Interpolation of Tensorial Components) method [Dvorkin1984] [Bathe1986] avoids this problem by resampling only the shear-strain components at predefined tying points in the element and interpolating those sampled values back over the element.

The nodes of an MITC shell element lie on the midsurface, and each node has six degrees of freedom: three translational components and three rotational components about the midsurface normal coordinate system. Element stiffness is evaluated by three-dimensional Gaussian integration over the midsurface natural coordinates and the thickness direction, and the thickness \(h\) is supplied as an element property during constitutive evaluation.

FrontISTR provides MITC3 (element ID 731), MITC4 (741), and MITC9 (743), which are represented by a single midsurface layer, as well as layered shell elements MITC3-shell361 (761, 3\(\times\)2 nodes, 3 DOFs per node) and MITC4-shell361 (781, 4\(\times\)2 nodes, 3 DOFs per node), whose nodes are arranged in two layers through the thickness. In the layered shell elements, nodal degrees of freedom consist only of three translational components, while bending modes corresponding to rotational degrees of freedom are represented by the two-layer arrangement.

Beam Elements

Line members such as girders and frame structures are discretized with beam elements. FrontISTR uses a Timoshenko beam formulation that accounts for shear deformation and expresses both bending and shear as functions of translational and rotational degrees of freedom.

Each beam node has six degrees of freedom: three translational components and three rotational components about the beam axis and transverse axes. Element stiffness is evaluated by one-dimensional numerical integration along the beam axis. The cross-sectional area \(A\) and second moments of area \(I\) in the bending and torsional directions are supplied as beam section properties, and together with the material Young's modulus \(E\) and shear modulus \(G\) they define the axial, bending, torsional, and shear stiffnesses.

FrontISTR provides a 2-node straight beam element (element ID 611) and a 4-node tetrahedral solid-beam hybrid element represented with 3 nodes (641, for mixed degrees of freedom).