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 Lagrangian, and Updated Lagrangian 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, Total Lagrangian, and Updated Lagrangian formulations.

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 Lagrangian, and Updated Lagrangian 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 Lagrangian, and Updated Lagrangian formulations. In the Updated Lagrangian formulation, 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. Beam elements represent extension, bending, and torsion using only translations and rotations of nodes on the beam axis, with deformation within the cross section represented by section constants.

Each beam node has six degrees of freedom: three translational components and three rotational components. Element stiffness is obtained by forming closed-form stiffness matrices for extension, bending in two directions, and torsion in a local coordinate system defined by the element axis and the reference axis specified in the section data, then transforming them to the global coordinate system. The cross-sectional area \(A\), second moments of area \(I_y\) and \(I_z\) about the local \(y\) and \(z\) axes, and the torsional constant are specified as section data; together with the material Young's modulus \(E\) and the shear modulus \(G\) derived from Young's modulus and Poisson's ratio, they define the respective stiffness coefficients.

For bending, two formulations are available: the Bernoulli-Euler beam, which assumes the cross section remains perpendicular to the beam axis after deformation and neglects shear deformation, and the Timoshenko beam, which accounts for shear deformation. Both can be selected for two-node beam element 611, with the Bernoulli-Euler beam as the default. For the Timoshenko beam, with element length \(L\), the coefficients representing the contribution of shear deformation,

\[ \Phi_y = \frac{12 E I_z}{G A_s L^2}, \qquad \Phi_z = \frac{12 E I_y}{G A_s L^2} \]

are included in the stiffness coefficients for the respective bending planes. For example, the stiffness coefficient for deflection in the local \(y\) direction is \(12 E I_z / \{ L^3 (1 + \Phi_y) \}\). Here \(A_s\) is the shear area. Because FrontISTR section data do not provide \(A_s\) separately, \(A_s = A\) is used. This is equivalent to taking the shear correction factor as 1. Setting \(\Phi_y = \Phi_z = 0\) gives the Bernoulli-Euler beam stiffness.

Element 641 is a mixed-DOF element that represents a beam using two nodes with translational degrees of freedom and two nodes with rotational degrees of freedom, and it is formulated as a Bernoulli-Euler beam.