Analysis of high-order quadrilateral plate elements based on the absolute nodal coordinate formulation for three-dimensional elasticity

Abstract

The absolute nodal coordinate formulation is a computational approach to analyze the dynamic performance of flexible bodies experiencing large deformations in multibody system dynamics applications. In the absolute nodal coordinate formulation, full three-dimensional elasticity can be used in the definition of the elastic forces. This approach makes it straightforward to implement advanced material models known from general continuum mechanics in the absolute nodal coordinate formulation. As, however, pointed out in the literature, the use of full three-dimensional elasticity can lead to severe locking problems, already present in simple, static tests. To overcome these drawbacks and to get a better understanding of these behaviors in the case of absolute nodal coordinate formulation elements, this study introduces and carefully analyses several high-order three-dimensional plate elements based on the absolute nodal coordinate formulation, primarily in meaningful static scenarios. The proposed elements are put to test in various numerical experiments intended to bring forward possible locking phenomena and to evaluate the accuracy attainable with the considered element formulations. The proposed eight- and nine-node elements that incorporate polynomial approximations of second order in all three directions prove to be advantageous both with respect to the actual performance and with regard to the numerical efficiency when compared to other absolute nodal coordinate formulation plate elements. A comparison with a four-node high-order element corroborates the supposition that the usage of in-plane slopes as nodal coordinates has a negative effect on numerical convergence properties in thin-plate use cases. An additional example showcases the functioning of two of the higher-order elements in a dynamic simulation.

Introduction

The absolute nodal coordinate formulation (ANCF) is a nonlinear finite element–based approach that is designed to describe flexible bodies undergoing large displacements and rotations in multibody system dynamics applications. ¹ In the ANCF, the kinematics of the finite element is defined by employing the spatial shape functions together with the element nodal coordinates. The nodes of the elements are located on the centerline for beams or on the mid-plane for plates. Contrary to beam and plate elements in traditional finite element analysis, ANCF beam and plate elements are described using absolute nodal positions and their derivatives. Using the components of the deformation gradient instead of the rotation angles as degrees of freedom, large three-dimensional rotations can be described without the singularity problems in a total Lagrangian formulation. Additionally, this choice of nodal coordinates leads to a constant mass matrix, which is computationally beneficial especially when employing explicit time integration schemes. With a properly selected material model and a suitable strain description, the ANCF can capture nonlinear strain–displacement relationships as well as material nonlinearities. For a more detailed discussion of the definition of ANCF elements, see Shabana. ²

The elastic forces can be defined using either structural-based approaches, which handle strain components independently, or full three-dimensional elasticity, for instance in the form of the St. Venant–Kirchhoff material law. However, the elastic forces must be defined with care as the use of full three-dimensional elasticity in conjunction with improper interpolations and choice of slope vectors has been demonstrated to be vulnerable for various locking phenomena. The shear-deformable four-node quadrilateral ANCF plate element proposed in Mikkola and Shabana, ³ for instance, suffers from transverse shear locking and Poisson locking. Poisson locking can be eliminated using the structural mechanics–based approach for the calculation of the elastic forces. With respect to advanced material modeling, however, the use of full three-dimensional elasticity is superior when compared to the structural-based approach. Accordingly, by eliminating drawbacks associated with the use of full three-dimensional elasticity, ANCF offers a straightforward platform for dynamic analysis based on advanced material models in the context of multibody systems. It is for this reason that the research of approaches to eliminate locking phenomena also for full three-dimensional elasticity is justified.

Kuhl and Ramm ⁴ give an overview of locking phenomena in traditional finite element analysis and existing techniques to prevent these phenomena. Apart from other strategies, it is mentioned that a natural approach to avoid locking phenomena is to develop high-order elements. Since ANCF elements are often quite closely related to their traditional finite element counterparts, it makes sense to consider this possibility also in the realm of ANCF. In the case of ANCF elements, higher approximation orders can be achieved using more nodes per element, by employing high-order derivatives as nodal coordinates, or using a combination of the two. This high-order approach has been examined thoroughly for ANCF beam elements, but comparatively little attention has been devoted to three-dimensional quadrilateral high-order plate elements. For the purpose of this article, elements are said to be of high order if they employ high-order derivatives as nodal coordinates.

Concerning ANCF beam elements, in Gerstmayr and Shabana, ⁵ the deterioration of the speed of convergence as a result of shear locking has been alleviated by enriching the polynomial basis with respect to the longitudinal coordinate. To reduce the influence of Poisson locking in the case of ANCF beam elements, Matikainen et al.^6,7 have proposed to include the trapezoidal cross-section mode using a second-order transversal derivative as a nodal coordinate. In Ebel et al. ⁸ and Shen et al., ⁹ many high-order beam elements were studied to demonstrate the above-mentioned locking problems.

Regarding plate elements, Olshevskiy et al. ¹⁰ introduce a triangular two-dimensional ANCF plate element making use of a second-order mixed derivative as an additional nodal coordinate. The numerical performance of further high-order triangular plate elements has also been investigated by Olshevskiy et al. ¹¹ Matikainen et al. ¹² have examined various quadrilateral ANCF plate elements, including one high-order four-node element that uses a second-order displacement interpolation in the thickness direction to successfully overcome Poisson locking. However, the element still suffers from an additional shear locking phenomenon limiting the numerical performance in the case of thin plates. There is a perception that the use of in-plane slope vectors together with derivatives in transversal directions in the nodal coordinates contributes to this problem. In previous studies, an approach that does not employ in-plane slope vectors has delivered promising results for three-node ANCF beam elements that make use of the St. Venant-Kirchhoff material law. ⁸

The objective of this article is to extend this recently introduced approach of the usage of a high polynomial order in all three directions without the use of in-plane slopes to plate applications. This extension makes it possible to employ full three-dimensional elasticity in the definition of the elastic forces. This, in turn, allows the implementation of advanced material models in the ANCF in a straightforward manner. To reach this objective, this article introduces eight- and nine-node high-order ANCF plate elements. The performance of the elements is systematically compared to low- and high-order four-node elements in several static tests that are designed to focus on certain locking phenomena. Furthermore, a dynamic test exemplifies the usage of two of the functioning higher-order elements in simulations over the time domain.

Definition of the considered ANCF plate elements

In the ANCF, the current position of an arbitrary particle $r\in {ℝ}^3$ at time $t\in {ℝ}_0^{+}$ can be written as

$$\mathbf{r}(\mathbf{x},t)=\mathbf{S}(\mathbf{x})\mathbf{e}(t)$$

(1)

where $\mathbf{S}(\mathbf{x})$ is the shape function matrix and $\mathbf{e}(t)$ is the vector of nodal coordinates. The vector $\mathbf{x}$ denotes the position of the considered particle in the reference configuration, and thus, the identity

$${\mathbf{r}}_0:=\mathbf{r}(\mathbf{x},0)=\mathbf{x}$$

holds.

In this article, various plate elements with eight and nine nodes will be considered. For comparison purposes, some results obtained with four-node elements will also be given. An exemplary nine-node element is illustrated in Figure 1.

Figure 1.

Illustration of a nine-node plate element with a particle p in reference and in current configuration. The nodes are denoted by I, II,…, IX. For plate elements with fewer nodes, the nodes with higher numbers are omitted.

For a plate element with b nodes, the vector of nodal coordinates can be written as

$$\mathbf{e}={[{\mathbf{e}}_{\mathrm{I}}^{\mathrm{T}}{\mathbf{e}}_{\mathrm{II}}^{\mathrm{T}}\cdots {\mathbf{e}}_b^{\mathrm{T}}]}^{\mathrm{T}}$$

(2)

With ${\mathbf{e}}_i$ , $i\in \{\mathrm{I},\mathrm{II},\dots ,b\}$ , being the nodal coordinate vector of the ith node. The nodal coordinates can contain the positions $\mathbf{r}$ of the element’s nodes in the inertial frame of reference and selected spatial derivatives of these position vectors. For the sake of readability, for the spatial derivatives, the shorthand notations

$${\mathbf{r}}_{,x}:=\frac{\partial \mathbf{r}}{\partial x},{\mathbf{r}}_{,y}:=\frac{\partial \mathbf{r}}{\partial y},{\mathbf{r}}_{,z}:=\frac{\partial \mathbf{r}}{\partial z},{\mathbf{r}}_{,\mathrm{zz}}:=\frac{{\partial }^2\mathbf{r}}{\partial z^2}$$

will be used in the formulation of the considered elements.

One possibility to derive the shape function matrix S is to employ a polynomial basis to approximate the particle positions of the considered element. In ANCF, the same basis is used for all three dimensions. ¹ Assuming that a monomial basis consisting of p monomials is used and that the vector ${\mathbf{b}}_p(\mathbf{x})\in {ℝ}^{1\times p}$ contains all of these monomials, the polynomial expansion can be written as

$$\mathbf{r}(\mathbf{x},t)=\left[\begin{array}{ccc}{\mathbf{b}}_p(\mathbf{x})\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\mathbf{b}}_p(\mathbf{x})\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & {\mathbf{b}}_p(\mathbf{x})\hfill \end{array}\right]\mathbf{c}(t)=:\mathbf{P}(\mathbf{x})\mathbf{c}(t)$$

(3)

where $\mathbf{c}$ is the vector of polynomial coefficients. Substituting the reference positions of the nodes into equation (3) and into those spatial derivatives of equation (3) that are used as nodal coordinates directly gives a linear relationship

$$\mathbf{e}(t)=\mathbf{Tc}(t)$$

(4)

between the polynomial coefficients and the nodal coordinates. Hence, if the polynomial basis has been chosen in a way that $\mathbf{T}$ is invertible, the shape function matrix reads

$$\mathbf{S}(\mathbf{x})=\mathbf{P}(\mathbf{x}){\mathbf{T}}^{-1}$$

(5)

Therefore, given a specific type of element with a specified number and arrangement of nodes, the shape function matrix is a direct consequence of the chosen polynomial basis and the chosen set of nodal coordinates. Accordingly, the subsequent section deals with the introduction of various ANCF plate elements by giving both the employed polynomial bases and nodal coordinates.

Derivation of the shape functions

To easily identify the various elements used in this article, a four-digit code similar to the one introduced in Dmitrochenko and Mikkola ¹³ will be used. For the four digits abcd, a signifies the dimension of the element, b the number of nodes, c the number of nodal coordinates per node and dimension and d will always be of the value 3 in the following, since the same polynomial basis is used for all three dimensions of the element.

In total, seven different ANCF plate elements are considered throughout this article. Table 1 depicts the nodal degrees of freedom and the polynomial bases of all seven elements, using the abbreviations

$$\begin{array}{c}{\beta }_{\mathrm{lin}.}=\{1,x,y,z,\mathrm{xy},\mathrm{xz},\mathrm{yz},\mathrm{xyz}\}\hfill \\ {\beta }_{\mathrm{qx}}=\{x^2,x^{2}y,x^{2}z,x^2\mathrm{yz}\}\hfill \\ {\beta }_{\mathrm{qy}}=\{y^2,x{y}^2,y^{2}z,x{y}^{2}z\}\hfill \\ {\beta }_{\mathrm{qz}}=\{z^2,x{z}^2,y{z}^2,\mathrm{xy}{z}^2\}\hfill \\ {\beta }_{\mathrm{a}}=\{x^2{y}^2,x^2{z}^2,y^2{z}^2,x^2{y}^{2}z,x^{2}y{z}^2,x{y}^2{z}^2,x^2{y}^2{z}^2\}\hfill \\ {\beta }_{\mathrm{b}}=\{x^2,y^2,x^3,y^3,x^{2}y,x{y}^2,x^{3}y,x{y}^3,z^2,x{z}^2,y{z}^2,\mathrm{xy}{z}^2\}\hfill \\ {\beta }^{[3933]}={\beta }_{\mathrm{lin}.}\cup {\beta }_{\mathrm{qx}}\cup {\beta }_{\mathrm{qy}}\cup {\beta }_{\mathrm{qz}}\cup {\beta }_{\mathrm{a}}\hfill \end{array}$$

(6)

Table 1.

Nodal coordinates and bases of all considered elements.

Element [dof]	Nodal coordinates	Basis	$n_x\times n_y\times n_z$
3933 [81]	$\mathbf{r},{\mathbf{r}}_{,z}$ , ${\mathbf{r}}_{,\mathrm{zz}}$	${\beta }^{[3933]}$	$5\times 5\times 5$
3923 [54]	$\mathbf{r},{\mathbf{r}}_{,z}$	${\beta }_{\mathrm{lin}.}\cup {\beta }_{\mathrm{qx}}\cup {\beta }_{\mathrm{qy}}\cup \{x^2{y}^2,x^2{y}^{2}z\}$	$5\times 5\times 3$
3833 [72]	$\mathbf{r},{\mathbf{r}}_{,z}$ , ${\mathbf{r}}_{,\mathrm{zz}}$	${\beta }^{[3933]}\backslash \{x^2{y}^2,x^2{y}^{2}z,x^2{y}^2{z}^2\}$	$5\times 5\times 5$
3823 [48]	$\mathbf{r},{\mathbf{r}}_{,z}$	${\beta }_{\mathrm{lin}.}\cup {\beta }_{\mathrm{qx}}\cup {\beta }_{\mathrm{qy}}$	$5\times 5\times 3$
3453 [60]	$\mathbf{r},{\mathbf{r}}_{,x},{\mathbf{r}}_{,y},{\mathbf{r}}_{,z}$ , ${\mathbf{r}}_{,\mathrm{zz}}$	${\beta }^{[3453]}:={\beta }_{\mathrm{lin}.}\cup {\beta }_{\mathrm{b}}$	$7\times 7\times 5$
3443 [48]	$\mathbf{r},{\mathbf{r}}_{,x},{\mathbf{r}}_{,y},{\mathbf{r}}_{,z}$	${\beta }^{[3453]}\backslash \{z^2,x{z}^2,y{z}^2,\mathrm{xy}{z}^2\}$	$7\times 7\times 3$
3423 [24]	$\mathbf{r}$ , ${\mathbf{r}}_{,z}$	${\beta }_{\mathrm{lin}.}$	$3\times 3\times 3$

The abbreviations used are defined in equation (6). The last column contains the number of integration points used in the evaluation of the elastic forces.

Two of the elements, namely 3933 and 3923, have nine nodes. Element 3933 is characterized by the fact that it uses an approximation of order 2 in all three directions. Element 3923 employs linear polynomials for the approximation in the z direction. Furthermore, two eight-node elements 3833 and 3823 are examined. They are similar to their nine-node counterparts, but since they lack a middle node, they make use of a basis containing three monomials less than elements 3933 and 3923, respectively.

For comparison, also the three four-node elements 3453, 3443, and 3423 are considered. Again, element 3453, originally proposed in Matikainen et al., ¹² uses an approximation of order 2 in the z direction, while element 3443 is linear in the z direction. On the other hand, element 3423, which is purely considered for demonstrative purposes, is linear in all directions.

For eight-node elements, there are, in theory, many different possibilities to choose valid bases in the sense that the matrix $\mathbf{T}$ from equation (4) is invertible. In this study, the two bases were chosen in such a way that the basis of the two-dimensional serendipity element Q8 from conventional finite element analysis is a subset of the bases of both element 3833 and element 3823.

The shape function matrices for all above-mentioned elements can now be calculated using equation (5). The shape functions can be given either in physical coordinates $\mathbf{x}=[x,y,z{]}^{\mathrm{T}}$ or in local element coordinates $\mathbf{\xi }\mathbf{=}\mathbf{[}\xi \mathbf{,}\eta \mathbf{,}\zeta {\mathbf{]}}^{\mathrm{T}}$ . In this study, the latter ones are scaled so that $-1\le \xi ,\eta ,\zeta \le 1$ . For a nine-node plate element, this is illustrated in Figure 2.

Figure 2.

Illustration of a nine-node plate element in physical and in local element coordinates.

Calculation of the elastic forces

In this paper, all considered elements make use of the St. Venant–Kirchhoff material law which assumes a linear relationship between the Green–Lagrange strain tensor $\mathbf{E}$ and the second Piola–Kirchhoff stress tensor ${\mathbf{S}}^{\mathrm{PK}}$ . ⁵ For implementation purposes, the integration to calculate the elastic forces is evaluated using Gaussian quadrature of a sufficient degree of exactness for all three dimensions independently. For the purpose of this article, the number of integration points is always set to be large enough such that no additional error is introduced into the calculation, see Table 1.

Numerical experiments

In this section, the above-mentioned plate elements are put to test in various numerical experiments. The goal of this analysis is to gain insight into which elements deliver accurate results while still maintaining acceptable convergence rates. Of special interest are the questions whether the introduction of second-order z-derivatives can alleviate certain locking phenomena and whether omitting in-plane slopes is actually beneficial for calculations involving thin plates. The first two tests that are conducted involve small deformations only. This allows for comparisons with approximative analytical solutions in order to eliminate the possibility of any fundamental flaws in element formulation and implementation. After that, two large-displacement tests are conducted. For these tests, reference results are calculated using finite elements provided by the commercial software ANSYS. The elements involved are the solid element SOLID95 and the beam element BEAM188. SOLID95 is a three-dimensional structural solid element with three nodal degrees of freedom for each of its 20 nodes. In general, this element has multiple reference shapes, but for these calculations, hexahedral elements were applied. BEAM188, on the other hand, is a three-dimensional shear-deformable linear beam element with two nodes and 6 degrees of freedom per node. For the calculations in this study, a four-cell cross-section was used. All numeric calculations involving ANCF elements have been performed using MATLAB.

Cantilevered plate

The first test involves a cantilevered plate of the length $L=1\mathrm{m}$ , the width $W=2\mathrm{m}$ , and the thickness $H=0.01\mathrm{m}$ . Young’s modulus amounts to $E=210\times {10}^9\mathrm{N}/{\mathrm{m}}^2$ , Poisson’s ratio is set to $\nu =0.3$ , and the shear modulus is given by $G=E/(2(1+\nu ))$ . A distributed force of $F_z=-1\mathrm{N}/\mathrm{m}$ is applied to the free edge of the mid-surface of the plate. Apart from the width W, the same cantilevered plate has been studied in Matikainen et al. ¹² This setup has the advantage that the obtained results can be compared with an analytical solution for the vertical displacement of the free end of a cantilevered plate of infinite width. ¹² This analytical solution is given as

$$u_z^{W\to \infty }=\frac{12(1-{\nu }^2)}{E{H}^3}\frac{F{L}^2}{3}\approx -1.7333\times {10}^{-5}\mathrm{m}$$

(7)

Naturally, since the plate of this experiment has a finite width, one cannot expect a result that is identical to the analytical solution. Nevertheless, the results should still be comparable. In particular, for a finitely wide plate with a non-zero Poisson’s ratio, the displacement in the middle of the free end will be larger than in the free corners of the plate. For this reason, the arithmetic means of the vertical displacements of all nodes at the end edge of the cantilever are calculated. Also, in the interest of a more thorough analysis, the nodal displacements with the smallest and largest absolute values are recorded.

The discretization has been chosen so that, as a first priority, all ANCF plate elements are always square-shaped and, as a second priority, so that the resulting number of degrees of freedom is similar for the different elements. The numeric results are depicted in Table 2.

Table 2.

Numeric results for the cantilevered plate experiment with Poisson’s ratio of $\nu =0.3$ .

Element [#, dof]	$u_z^{\mathrm{absmin}}$ (m)	$u_z^{\mathrm{av}}$ (m)	$u_z^{\mathrm{absmax}}$ (m)
3933 [ $9\times 18$ , 5994]	$-1.6840\times {10}^{-5}$	$-1.7415\times {10}^{-5}$	$-1.7679\times {10}^{-5}$
3923 [ $11\times 22$ , 5940]	$-1.3902\times {10}^{-5}$	$-1.4638\times {10}^{-5}$	$-1.4967\times {10}^{-5}$
3833 [ $10\times 20$ , 5580]	$-1.6862\times {10}^{-5}$	$-1.7439\times {10}^{-5}$	$-1.7700\times {10}^{-5}$
3823 [ $12\times 24$ , 5328]	$-1.3906\times {10}^{-5}$	$-1.4644\times {10}^{-5}$	$-1.4972\times {10}^{-5}$
3453 [ $13\times 26$ , 5346]	$-1.6772\times {10}^{-5}$	$-1.7327\times {10}^{-5}$	$-1.7589\times {10}^{-5}$
3443 [ $15\times 30$ , 5673]	$-1.3904\times {10}^{-5}$	$-1.4625\times {10}^{-5}$	$-1.4959\times {10}^{-5}$
3423 [ $22\times 44$ , 5940]	$-2.0257\times {10}^{-6}$	$-2.0513\times {10}^{-6}$	$-2.0667\times {10}^{-6}$
Reference $(W\to \infty )$	$-1.7333\times {10}^{-5}$	$-1.7333\times {10}^{-5}$	$-1.7333\times {10}^{-5}$

The element mesh and numbers of degrees of freedom used are given in square brackets. The values $u_z^{\mathrm{av}}$ represent the arithmetic means of the displacement values for all nodes at the free end of the plate. On the other hand, $u_z^{\mathrm{absmin}}$ and $u_z^{\mathrm{absmax}}$ are the nodal displacements that have the smallest or largest absolute value.

As can be seen from the results, for elements 3933, 3833, and 3453, the analytical result is well in the range of the nodal displacements. These elements have in common that they employ a second-order approximation in the z direction. In contrast, all elements that are merely linear in the z direction deliver inaccurate results. However, since elements 3923, 3823, and 3443 still give similar values when compared with each other, it seems that a common phenomenon must be responsible for the fact that they converge to the wrong value. Indeed, when repeating the experiment with Poisson’s ratio set to be equal to zero, also elements 3923, 3823, and 3443 deliver correct results. This can be seen in Table 3.

Table 3.

Numeric results for the cantilevered plate experiment with Poisson’s ratio of $\nu =0$ .

Element [#, dof]	$u_z^{\mathrm{absmin}}$ (m)	$u_z^{\mathrm{av}}$ (m)	$u_z^{\mathrm{absmax}}$ (m)
3933 [ $9\times 18$ , 5994]	$-1.8994\times {10}^{-5}$	$-1.8994\times {10}^{-5}$	$-1.8994\times {10}^{-5}$
3923 [ $11\times 22$ , 5940]	$-1.9013\times {10}^{-5}$	$-1.9013\times {10}^{-5}$	$-1.9013\times {10}^{-5}$
3833 [ $10\times 20$ , 5580]	$-1.9005\times {10}^{-5}$	$-1.9005\times {10}^{-5}$	$-1.9005\times {10}^{-5}$
3823 [ $12\times 24$ , 5328]	$-1.9020\times {10}^{-5}$	$-1.9020\times {10}^{-5}$	$-1.9020\times {10}^{-5}$
3453 [ $13\times 26$ , 5346]	$-1.9020\times {10}^{-5}$	$-1.9020\times {10}^{-5}$	$-1.9020\times {10}^{-5}$
3443 [ $15\times 30$ , 5673]	$-1.9027\times {10}^{-5}$	$-1.9027\times {10}^{-5}$	$-1.9027\times {10}^{-5}$
3423 [ $22\times 44$ , 5940]	$-1.6812\times {10}^{-6}$	$-1.6812\times {10}^{-6}$	$-1.6812\times {10}^{-6}$
Reference $(W\to \infty )$	$-1.9048\times {10}^{-5}$	$-1.9048\times {10}^{-5}$	$-1.9048\times {10}^{-5}$

In conclusion, these elements suffer from Poisson locking, and similar to the findings for beam elements in previous studies, approximations of at least second order in all directions seem to be a functional remedy for Poisson locking. ⁸ The linear element 3423, however, suffers from a variety of severe locking phenomena and hence delivers results with absolute values that are approximately an order of magnitude smaller than the reference value.

Simply supported plate under body load

The second numeric experiment is a square, simply supported plate subjected to a body load as previously studied in Matikainen et al. ¹² The plate is of the length and width $L=W=1\mathrm{m}$ , Young’s modulus is set to $E=2.1\times {10}^{11}\mathrm{N}/{\mathrm{m}}^2$ , and Poisson’s ratio is $\nu =0.3$ . The calculations are performed for three different thickness values, namely $H_1=0.1\mathrm{m}$ , $H_2=0.01\mathrm{m}$ , and $H_3=0.005\mathrm{m}$ . The uniform body load amounts to

$$\mathbf{b}={\left[\begin{array}{ccc}0\hfill & 0\hfill & -5\times {10}^6·H_i^2\hfill \end{array}\right]}^{\mathrm{T}}$$

for the respective thickness values. The setup is illustrated in Figure 3.

Figure 3.

Illustration of the setup for the simply supported plate.

The boundary conditions are given by

$$\begin{array}{c}(\mathbf{r}{)}_{1,3}{|}_{{\Gamma }_{\mathrm{le}}}=0,(\mathbf{r}{)}_{1,3}{|}_{{\Gamma }_{\mathrm{ri}}}=0\hfill \\ (\mathbf{r}{)}_{2,3}{|}_{{\Gamma }_{\mathrm{lo}}}=0,(\mathbf{r}{)}_{2,3}{|}_{{\Gamma }_{\mathrm{up}}}=0\hfill \end{array}$$

(8)

The displacements are recorded in the middle of the plate. Once again, there exists a thin-plate analytical solution for this problem. Using Reissner–Mindlin plate theory, ¹² this analytical solution of the z displacement in the middle of the plate is of the form

$$\begin{array}{c}\multicolumn{1}{c}{u_z^{\mathrm{RM}}=\frac{12(1-{\nu }^2)}{{\pi }^4{\mathrm{EH}}_i^3}\sum _{p=1}^{\infty }\sum _{q=1}^{\infty }\left[\frac{\sin (p\pi /2)\sin (q\pi /2)}{{(p^2/L^2+q^2/W^2)}^2}\frac{4b_3{H}_i}{{\pi }^2\mathrm{pq}}(1-\cos (\pi p))(1-\cos (\pi q))\right]}\\ \multicolumn{1}{c}{+\frac{1}{G{H}_i{\pi }^2}\sum _{p=1}^{\infty }\sum _{q=1}^{\infty }\left[\frac{\sin (p\pi /2)\sin (q\pi /2)}{p^2/L^2+q^2/W^2}\frac{4b_3{H}_i}{{\pi }^2\mathrm{pq}}(1-\cos (\pi p))(1-\cos (\pi q))\right]}\end{array}$$

(9)

The results are given in Table 4. The reference values given have been obtained by evaluating the sums in equation (9) up to $p=q=1024$ , which gives results that can be considered as exact with respect to the given accuracy.

Table 4.

Results for the center displacements of the simply supported plate under body load for various thickness values.

Element [#, dof]	$u_z$ , $H_1=0.1\mathrm{m}$	$u_z$ , $H_2=0.01\mathrm{m}$	$u_z$ , $H_3=0.005\mathrm{m}$
3933 [ $16\times 16,9541$ ]	1.0782	1.0019	0.9995
3923 [ $20\times 20,9762$ ]	0.8722	0.8183	0.8165
3833 [ $18\times 18,9113$ ]	1.0787	1.0025	0.9999
3823 [ $22\times 22,8890$ ]	0.8723	0.8186	0.8166
3453 [ $24\times 24,9179$ ]	1.0717	1.0002	0.9963
3443 [ $28\times 28,9864$ ]	0.8693	0.8171	0.8151
3423 [ $40\times 40,9762$ ]	0.8631	0.4344	0.1801

The results are normalized with respect to the reference values of $-1.1018\times {10}^{-6}\mathrm{m}$ , $-1.0567\times {10}^{-6}\mathrm{m}$ , and $-1.0563\times {10}^{-6}\mathrm{m}$ for $H_1$ , $H_2$ , and $H_3$ .

As can be gathered from the results for all three thickness values, the simply supported plate experiment reinforces the impressions from the previous test. Once more, the elements 3933, 3833, and 3453 deliver correct results. The elements 3923, 3823, 3443, and 3423, on the other hand, are once again incapable of delivering correct displacements. Thus, in the interest of clarity, the remainder of this article will show only results of the three better-performing ANCF elements.

In practical applications, it is not only important that the numerical solution converges to a correct value for increasing numbers of elements, but rather that a satisfactory solution can already be obtained for a modest number of degrees of freedom. Therefore, in the following, the convergence rates for elements 3933, 3833, and 3453 will be investigated. Figure 4(a) and (b) show the relative errors for the thicknesses $H_1$ and $H_2$ .

Figure 4.

In this figure, the relative errors of the displacement results of the simply supported plate under body load are plotted for different numbers of elements. For the calculation of the relative error, the values given in Table 4 were used as the respective reference solutions: (a) H₁ = 0.1 m, (b) H₂ = 0.01 m, and (c) H₃ = 0.005 m.

While for $H_1=0.1\mathrm{m}$ , the three elements exhibit quite similar convergence behaviors, for the thinner plate of thickness $H_2=0.01\mathrm{m}$ , element 3453 seems to suffer from a locking phenomenon, since for the first few hundred degrees of freedom, the convergence rate is noticeably decreased. This can already be of practical relevance. While element 3933 delivers quite reasonable results for less than 2000 degrees of freedom with a relative error noticeably below 1%, element 3453 needs approximately twice as many degrees of freedom to achieve a relative error as small as that.

Even more so, when using exact Gaussian quadrature for all elements, the computation effort per degree of freedom for calculations with element 3453 is larger than for elements 3933 and 3833. The reason for this lies in the fact that the basis used for element 3453 contains monomials up to the third order of x while the other two elements only use monomials up to second order. Hence, for an exact quadrature, element 3453 needs to use more integration points in the x direction.

Furthermore, as demonstrated by Figure 4(c), this locking phenomenon is even more severe for the thinnest considered plate of thickness $H_3=0.005\mathrm{m}$ . This phenomenon is usually referred to as shear locking or transverse shear locking in the literature since it involves artificial shear strains that render the element to be overly stiff. In conventional finite element analysis, this phenomenon appears for instance in linear plate elements since then linear shape functions are used for both displacements and rotations and thus, pure bending cannot be replicated without artificial shear strains. ¹⁴

However, as the presented results suggest, in the case of ANCF plate elements, it seems that this phenomenon is also connected to the usage of the in-plane slope vectors ${\mathbf{r}}_x$ and ${\mathbf{r}}_y$ as nodal coordinates and hence the enforcement of the continuity of said slopes between the elements. This fits to the findings of Matikainen et al. ¹² where the element 3453 was also put to the test. The absence of noticeable shear locking in the results of the elements 3833 and 3933 shows that the usage of approximation orders of at least two in all directions without the employment of in-plane slopes as nodal coordinates seems to be a good practice to yield ANCF plate elements that suffer neither from shear nor from Poisson locking.

Large-displacement cantilever

All previous experiments have dealt with load cases that lead to small displacements. For this reason, the intention of the following experiment is to clarify whether the elements 3933, 3833, and 3453 still deliver reasonable results in the case of large deformations. To this end, a cantilevered plate of the length $L=0.5\mathrm{m}$ , the width $W=0.15\mathrm{m}$ , and the thickness $H=0.001\mathrm{m}$ is considered. The same plate was examined in Dufva and Shabana. ¹⁵ Young’s modulus is set to $E=2.07\times {10}^{11}\mathrm{N}/{\mathrm{m}}^2$ with Poisson’s ratio being $\nu =0.3$ . The total applied force is $F_z=30\mathrm{N}$ . As illustrated in Figure 5, this force is split equally on the two free corners of the plate.

Figure 5.

Illustration of the setup for the large-displacement cantilever experiment for a total applied force of $F_z$ .

For the ANCF elements, the mesh was set up in a way that always one element was used in the y direction while the number of elements in the x direction was increased step by step. The displayed reference result was obtained with the ANSYS solid element SOLID95 and a mesh consisting of 334 elements in the x direction, 100 elements in the y direction, and 1 element in the z direction.

The results of the numerical calculations for different numbers of degrees of freedom are depicted in Figure 6. For all elements, the given displacement value is the arithmetic mean of the vertical displacements of all nodes in the end cross-section.

Figure 6.

Results of the large-displacement cantilever test for various elements.

The results show that all of the high-order elements deliver reasonable results. None of the elements shows any sign of noteworthy convergence problems.

Princeton beam experiment

The previous examples have dealt with either small-displacement or simple one-dimensional bending problems. Since the proposed elements are based on full three-dimensional elasticity, they should in theory be usable in a wide range of problems. Hence, in this section, the plate elements are put to the test in the so-called Princeton beam experiment which is a demanding example involving bidirectional bending. First conceived by researchers from Princeton University, ¹⁶ this experiment has previously been used to evaluate different beam formulations, including ANCF beam elements. ¹⁷

The cantilevered beam is of the length $L=0.508\mathrm{m}$ and has a rectangular cross-section of the width $W=12.377\times {10}^{-3}\mathrm{m}$ and the height $H=3.2024\times {10}^{-3}\mathrm{m}$ . Young’s modulus and Poisson’s ratio amount to $E=71.7\times {10}^9\mathrm{N}/{\mathrm{m}}^2$ and $\nu =0.31$ , respectively. The beam is subjected to a concentrated tip load acting on the center of the end cross-section. The force vector is parallel to the y-z-plane. Its direction, as defined by the angle $\alpha$ that is included between the force vector and the direction of the y-axis, is varied in $15\circ$ steps from $0\circ$ to $90\circ$ . The experiment is performed for three different absolute values of the force, namely, $F_1^{\mathrm{pb}}=4.448\mathrm{N}$ , $F_2^{\mathrm{pb}}=8.896\mathrm{N}$ , and $F_3^{\mathrm{pb}}=13.345\mathrm{N}$ . The setup is illustrated in Figure 7.

Figure 7.

Illustration of the setup for the Princeton beam experiment for a tip load of $F^{\mathrm{pb}}$ .

For all ANCF elements involved, two elements were used in the y direction so that for all element types, a node exists that is precisely in the middle of the end cross-section. The recorded displacements are those of that middle node. For the elements 3933, 3833, and 3453, respectively, 20, 25, and 40 elements were used in the x direction. As a result, for all three elements, the system has approximately 1800 degrees of freedom. The results for the different load cases are depicted in Figures 8 and 9. For the reference results, 64 elements of the ANSYS beam element BEAM188 and a mesh of $160\times 8\times 2$ elements of the ANSYS solid element SOLID95 were used.

Figure 8.

Results for the chordwise displacements in the Princeton beam experiment.

Figure 9.

Results for the flapwise displacements in the Princeton beam experiment.

Both the results for the chordwise and flapwise displacements show that all of the elements that have passed the previous tests also deliver reasonable results in all load cases of the Princeton beam experiment. This is of particular significance since the Princeton beam experiment is a rather elaborate test case involving nonlinear effects based on the coupling of axial and transverse normal strains. This becomes manifest in the fact that as the angle $\alpha$ increases from $0\circ$ to $30\circ$ , the absolute value of the displacement in the y direction also increases while the force in that direction decreases. However, all considered elements that involve approximation orders of at least two in each direction are able to reproduce this effect. Only upon closer inspection of the results, slight differences between the elements become apparent.

In Figure 10, the y displacements are plotted over different numbers of degrees of freedom. The results have been obtained for the load case with an angle of $\alpha =30\circ$ and the largest force of magnitude $F_3^{\mathrm{pb}}$ , because for these parameters, the nonlinear contributions are most significant and the differences between the elements are most apparent.

Figure 10.

Results for the flapwise displacement of the Princeton beam experiment with an angle of $\alpha =30\circ$ and a force of absolute value $F_3^{\mathrm{pb}}$ . The displayed reference result has been obtained with the ANSYS solid element SOLID95 and a mesh of $160\times 8\times 2$ elements.

As the results show, the elements 3933 and 3833 actually seem to be able to deliver results that are in better concordance with those obtained with the ANSYS solid element SOLID95. The displacement for element 3453, on the other hand, seems to converge to an absolute value that is slightly smaller.

Dynamic test

The example case selected for the dynamic test is the moderately thick free-flying plate problem presented in Figure 11, originally developed by Matikainen et al. ¹² The example was solved using the higher-order ANCF plate elements 3833 and 3933, which have been chosen due to their good performance in the previously presented static examples. The material parameters used in the test are as follows: Young’s modulus is $E={10}^6\mathrm{N}/{\mathrm{m}}^2$ , Poisson’s ratio amounts to $\nu =0.3$ , the density is $\rho =1500\mathrm{kg}/{\mathrm{m}}^3$ , and shear modulus and correction factor are set to $G=E/(2(1+v))$ and $k_s=5/6$ , respectively. The studied plate is of the length $L=1\mathrm{m}$ , width $W=0.5\mathrm{m}$ , and thickness $H=0.1\mathrm{m}$ . The force components $F_x$ and $F_y$ only affect the node at point a, while the moment $M_z$ is distributed evenly among the nodes on the line ab. As given by

$$F_x(t)=\{\begin{array}{cc}0\hfill & \mathrm{if}0\le t\le 1.5\hfill \\ \multicolumn{1}{c}{13.5}& \mathrm{if}1.5<t\le 2.5\hfill \end{array}(\mathrm{N})$$

(10)

$$F_y(t)=\{\begin{array}{cc}0\hfill & \mathrm{if}0\le t\le 1.5\hfill \\ \multicolumn{1}{c}{18(t-1.5)}& \mathrm{if}1.5<t\le 2.5\hfill \end{array}(\mathrm{N})$$

(11)

$$M_z(t)=\{\begin{array}{cc}-47.5t\hfill & \mathrm{if}0\le t\le 1.5\hfill \\ \multicolumn{1}{c}{0}& \mathrm{if}1.5<t\le 2.5\hfill \end{array}(\mathrm{N}\mathrm{m})$$

(12)

the loading of the plate varies as a function of time. The total duration of the simulation is 2.5 s. The results are computed using MATLAB’s integrated solver ode15s, which is a variable order integrator for stiff problems, with an output step length of ${10}^{-3}\mathrm{s}$ . The integration of the elastic forces of each element is done via Gaussian quadratures of order 5. The absolute and relative error tolerances for the solution of the equations of motion are both set to $5\times {10}^{-4}$ .

Figure 11.

Free flying flexible plate.

The reference results are calculated using a $256\times 128$ mesh of SHELL181 elements in ANSYS Mechanical APDL version 17.2. Geometrical nonlinearities are taken into account in the model. The Newmark algorithm is employed in the solution process, with the integration parameters $\delta =1/2$ and $\alpha =1/4$ . The time steps are of variable length with a starting and maximum step length of ${10}^{-2}\mathrm{s}$ and a minimum step length of ${10}^{-4}\mathrm{s}$ . The full Newton–Raphson method is employed in the solution, with convergence criteria settings left at the ANSYS default values, that is, employing the Euclidean norm with tolerance ${10}^{-3}$ for force and moment.

The absolute error of the displacement of point a at time $t=0.75\mathrm{s}$ for various mesh sizes is presented in Figure 12. The error norms at the later time $t=2.5\mathrm{s}$ are depicted in Figure 13, while the computational effort required for the computations is compared in Table 5. It should be noted that the computation time values presented in the table are only suggestive, since it could not be ensured that the computation was not hindered by other processes on the machine.

Figure 12.

Absolute error of the displacement of point a at $t=0.75\mathrm{s}$ . The reference results were computed using a $256\times 128$ mesh of SHELL181 elements in ANSYS.

Figure 13.

Norm of the displacement error of point a at $t=2.5\mathrm{s}$ . The reference results were computed using a $256\times 128$ mesh of SHELL181 elements in ANSYS.

Table 5.

Comparison of the computational effort for the simulation over the time interval $t\in [0,2.5]$ (s).

Mesh	Element 3833			Element 3933
$x\times y$	DOF	Eval.	Time	DOF	Eval.	Time
$2\times 1$	117	33,870	1.00	135	36,229	1.17
$4\times 2$	333	41,665	4.98	405	38,591	5.97
$6\times 3$	657	51,576	19.6	819	52,371	27.9
$8\times 4$	1089	60,397	57.9	1377	56,681	72.7

DOF: degree of freedom.

“Eval.” denotes the number of function evaluations of the equations of motion and “Time” denotes total computation time, scaled with respect to the smallest value.

While Figure 12 shows a certain degree of difference between the ANSYS model and the higher-order plate element models, particularly in the y-directional displacement, it can be concluded from the simulation that generally the higher-order plate models behave very similarly to the ANSYS model. In particular, the plate’s rotation is reproduced in the proper way without incorrect strains. Predictably, Table 5 shows that the more complex element 3933 is more computationally demanding, while, according to Figures 12 and 13, element 3833 converges to roughly the same error value with fewer degrees of freedom.

Conclusion

This study introduces high-order three-dimensional plate elements based on the ANCF for large deformation multibody applications. The introduced elements employ eight or nine nodes.

Contrary to quadrilateral ANCF plate elements considered in previous studies, the two proposed eight- and nine-node high-order plate elements deliver reasonable results in all considered load cases while exhibiting satisfactory convergence rates. Specifically, both Poisson and shear locking are resolved through the setup of the two elements. In this study, Poisson locking disappears for all elements that have at least an approximation order of two in all directions. In addition, the absence of in-plane slope vectors in the nodal coordinates of the two elements seems to lead to elements that do not suffer from severely deteriorated convergence rates in case of very thin plates.

The considered load cases include two small-displacement tests with one being a cantilevered plate and the other one being a simply supported plate under a uniform body load. For the latter one, an examination of the convergence rates for different thicknesses revealed the advantages of the considered high-order eight- and nine-node elements when compared to the four-node fully parametrized high-order ANCF plate element that has been introduced in previous studies. ¹²

Furthermore, two large-displacement load cases, namely one large-displacement cantilevered plate in one-dimensional bending and the so-called Princeton beam experiment involving bidirectional bending, have shown that the proposed elements also behave favorably in more elaborate situations. Finally, the two proposed eight- and nine-node elements that behaved well in the static tests have been applied in an exemplary dynamic benchmark test, showing that they are also functional in simulations over the time domain.