1 Introduction

In the modern design of the software and applications, the highly flexible bodies built with nonlinear materials are in common use. Consequently, the multibody system (MBS) simulation software should be able to include such bodies in a reliable manner. The most commonly used method in flexible body simulations is the floating frame of reference formulation [33]. This method is usually able to include only small deformations and linear material models. The large deformations and nonlinearities might be modeled with the finite element analysis (FEA) [4], but the FEA does not work perfectly with the MBS [42].

Absolute nodal coordinate formulation (ANCF) [32] is one of the methods that can be used for efficient dynamic simulation of flexible multibody systems. This formulation has some unique features allowing for modeling nonlinear material models with beam elements in a straightforward manner. In the ANCF, to describe local orientation, the slope coordinates are used instead of rotations. This allows complex shapes to be represented using a small number of elements. The elastic forces of the beam elements can be derived with not only the classical beam theories (such as the Bernoulli beam model) but also the continuum mechanics approach may be employed. This allows for use of general constitutive equations for the nonlinear (both compressible and incompressible) material models.

The absolute nodal coordinate formulation is well suited to model flexible bodies in the framework of the multibody systems. The ANCF can easily be used with the nonincremental formulation [34]. Flexible ANCF bodies may undergo large deformations, as well as large translational and rotational displacements. Additionally, it allows for exact modeling of the rigid body motions. This characteristic makes the ANCF method well suited for the dynamic simulations of the slender, flexible bodies modeled with nonlinear material models.

Efficient and reliable modeling of the nonlinear, incompressible materials is an important task in modern industrial and engineering applications. Many areas of engineering and biomechanics account for the material incompressibility. For example, the vessel wall material, heart muscle, as well as some bioelastomers, such as elastin and resilin, are treated as incompressible materials [10]. The rubber-like materials are used in many practical applications such as automotive (tires, rubber seals, boots), defense (coating for equipment, radar absorbing materials), safety (seat belts, impact absorbers) and many others. Incompressible materials have to be modeled in commercial finite element packages with shell or solid elements [2], even for the slender structures that are usually approximated with beams. The use of the ANCF in such application allows overcoming this limitation by using the fully parameterized beam elements [34].

Most of the papers devoted to the ANCF assume linear constitutive relation between stresses and strains [26]. The isotropic and hyperelastic material models (including incompressible models) were introduced into ANCF by Maqueda and Shabana [23]. Moreover, the applications of the incompressible materials to the rubber chains build with ANCF beam elements were presented [22]. In addition, Jung et al. [18] present the experimental validation of the rubber-like beam modeled with ANCF elements, although they do not highlight the importance of the volumetric locking suppression methods. Furthermore, the formulations of plasticity were examined within the framework of the ANCF [41]; however, this topic will not be further examined in the paper.

The main objective of this paper is to examine the use of nonlinear, hyperelastic, nearly incompressible material models using the absolute nodal coordinate formulation beam elements. While this objective was addressed in previous publications [23], the main contribution of this paper is to show the importance of volumetric locking elimination in the applications when bending-dominated models are considered. For a detailed overview of the locking phenomena occurring within ANCF framework as well as general review of developments in ANCF, consult [14, 26]. Moreover, exemplary techniques that may reduce the locking influence on the results are presented and validated with several numerical examples. Two basic incompressible material models are presented, incompressible Neo-Hookean and two-parameter Mooney–Rivlin [34]; however, presented methods may be easily extended to the higher-order material models. The incompressibility of the materials is ensured in the current study by the penalty method, which is chosen due to its popularity and simplicity. Furthermore, four fully parameterized elements are employed as they allow for a straightforward implementation of the constitutive nonlinear models. The first element is standard beam element with twenty-four nodal coordinates introduced in [35, 45] that is referred in the following as beam element 24 (BE24). The second element is higher-order beam with thirty nodal coordinates that suppress shear locking introduced in [13]. This element is denoted in the paper as beam element 30 with shear locking suppression (BE30S). The third element is higher-order element with thirty nodal coordinates that include trapezoidal deformation mode of the cross section that suppresses volumetric locking introduced in [24, 25]. This element is referred in the following as beam element 30 with trapezoidal mode (BE30T). The last element is higher-order beam element with forty-two nodal coordinates that include full quadratic approximation in transversal direction introduced in [36] where it is denoted as B2 and further investigated in [28]. This element is denoted in the following as beam element 42 (BE42). Similarly like BE30T, BE42 suppresses the volumetric locking phenomena using higher-order approximation in the cross section.

Within the finite element method, a wide variety of element technologies is known which main goal is to obtain a general element formulation that provides a good bending behavior, no locking for incompressible materials and very thin elements, efficiency, distortion insensitivity and good coarse mesh accuracy [44]. The notable and important methods that are used to solve the issues associated with lockings are the selective and reduced integration techniques [21] which are effective and easy to implement. Alternatively, a family of hybrid and mixed finite elements has been proposed [5]. The most commonly known mixed formulations are based on the two-field variational Hellinger–Reissner principle [30] and three-field Hu–Washizu principle [43]; however, the development of the finite elements via these principles requires the approximation of additional tensor fields, and the resulting number of scalar variables might be very large. Therefore, these principles are seldom implemented directly [5]. Noteworthy, the shear locking suppression methods based on these variational principles were employed in elastic line and mid-surface approaches with beam [31] and plate [7] elements. More recently, the simpler, pure displacement approaches could extend the selective reduced integration method to more complex cases emerged. Among them are the B-bar strain projection method [16, 37] and an enhanced assumed strain formulation [38]. An important advantage of the B-bar approach over the enhanced formulation is that it only results in a modified strain–displacement matrix with no additional degrees of freedom or parameters [5]. The enhanced assumed strain method was applied to the ANCF beam element to alleviate lockings by Gerstmayr and Shabana [12], and more recently, it was used with shell element by Yamashita et al. [46]. In the current paper, it is shown that the volumetric locking suppression is crucial for the BE24 and BE30S elements and two methods of the locking alleviation are presented. The first one is based on the selective reduced integration, and it is similar to the Poisson locking elimination method used with linear materials within ANCF framework [11]. The second method is based on the B-bar strain projection technique. Both methods are presented in details and verified with numerical examples.

The paper is organized as follows. Section 2 presents the basic kinematic equations of the spatial absolute nodal coordinate formulation beam elements. In addition, the standard form of the equations of motion is shown. Section 3 introduces the well-known incompressible material models, the penalty method that ensures the incompressibility and modified version of the nonlinear material models. The formulation of the issues connected with volumetric locking phenomena is depicted in Sect. 4. This section also presents the proposed methods for suppressing the volumetric locking. In Sect. 5, several numerical examples of the models that incorporate incompressible material with and without locking elimination techniques are presented. Conclusions of this study are presented in Section 6.

2 Kinematics of the deformable bodies using the ANCF

In the absolute nodal coordinate formulation, nodal coordinates consist of the translational and slopes coordinates, described with respect to the global coordinate system. There are no rotational coordinates, such as Euler parameters or Euler angles, used to describe the orientation of the element. Due to this feature, only the displacement field is interpolated, and there is no need for the independent rotation field interpolation [40]. For example, the node of the standard beam element BE24 has the following coordinates:

$$\begin{aligned} \varvec{e}^{i}=\begin{bmatrix}{\varvec{r}^{i}}^\mathrm{T}&{\varvec{r}_{,x}^{i}}^\mathrm{T}&{\varvec{r}_{,y}^{i}}^\mathrm{T}&{\varvec{r}_{,z}^{i}}^\mathrm{T}\end{bmatrix}^\mathrm{T} \end{aligned}$$
(1)

where vector \(\varvec{e}^{i}\) contains nodal parameters of the i node, vector \(\varvec{r}^{i}\) is the global position vector of the i node, while vectors \({\varvec{r}_{,k}^{i}}=\partial \varvec{r}^{i}/\partial k\) for \(k=x,y,z\) are first-order gradients of the position vector (slope coordinates) of the i node. The vector of the nodal coordinates for single BE24 element is \(\varvec{e}=\begin{bmatrix}{\varvec{e}^{A}}^\mathrm{T}&{\varvec{e}^{B}}^\mathrm{T}\end{bmatrix}^\mathrm{T}\) where A and B denote element nodes. All elements considered in the paper are beam elements with two nodes; however, they employ higher-order gradients as nodal coordinates. For detailed description of the kinematics of these elements, see [45] for BE24, [13] for BE30S, [25] for BE30T and [36] for BE42. Using the vector of nodal coordinates \(\varvec{e}\), one can get the position of the arbitrary point on the element:

$$\begin{aligned} \varvec{r}=\varvec{S}\varvec{e} \end{aligned}$$
(2)

where \(\varvec{S}\) is the matrix that contains element shape functions and depends only on the spatial coordinates. Usually, shape functions are written as a function of element dimensionless coordinates \(\varvec{S}=\varvec{S}(\xi ,\eta ,\zeta )\), where \(\xi =x/l\), \(\eta =y/l\) and \(\zeta =z/l\) for l denoting the length of the element in the undeformed state.

The kinetic energy of the single element is defined by the following formula:

$$\begin{aligned} T=\frac{1}{2}\int _{V}{\rho {\dot{\varvec{r}}}^\mathrm{T}\dot{\varvec{r}}\hbox {d}V}=\frac{1}{2}{\dot{\varvec{e}}}^\mathrm{T}\varvec{M}\dot{\varvec{e}} \end{aligned}$$
(3)

where \(\rho \) and V are, respectively, density and volume of the element and \(\varvec{M}\) is the mass matrix of the element. After differentiation of Eq. (2) with respect to time, mass matrix can be written as follows:

$$\begin{aligned} \varvec{M}=\int _{V}{\rho \varvec{S}^\mathrm{T}\varvec{S}\hbox {d}V} \end{aligned}$$
(4)

The mass matrix of ANCF beam elements is constant.

In the absolute nodal coordinate formulation, Coriolis, tangential, centrifugal, as well as other forces resulting from differentiation of the kinetic energy are equal to zero. Thus, the nonzero terms in system equations of motion come from the vectors of the elastic and external forces. The vector of the external forces, which include the gravitational forces, can be derived using the principle of virtual work as follows:

$$\begin{aligned} \delta W_{e}={\varvec{Q}_{e}}^\mathrm{T}\delta \varvec{e} \end{aligned}$$
(5)

where \(\delta W_{e}\) is the virtual work of the external forces and \(\varvec{Q}_{e}\) is the vector of the generalized external forces.

Vector of the elastic forces is given by:

$$\begin{aligned} \varvec{Q}_{s}=\left( \frac{\partial U_{s}}{\partial \varvec{e}}\right) ^\mathrm{T} \end{aligned}$$
(6)

where \(U_{s}\) is the energy of the elastic deformation (for linear elastic material model) or elastic energy density function (e.g., for the isotropic, hyperelastic material models). Elastic energy of the linear elastic material model may be found in the literature [34]. The present paper is mainly focused on the nonlinear, hyperelastic, isotropic and nearly incompressible material models.

Finally, one can write the equations of motion in the independent coordinates for the three-dimensional ANCF beam as follows:

$$\begin{aligned} \varvec{M}\ddot{\varvec{e}}+\varvec{Q}_{s}=\varvec{Q}_{e} \end{aligned}$$
(7)

3 Existing nonlinear constitutive models of the incompressible materials

In the current section, a polynomial hyperelastic material model is shown in the form of the Mooney–Rivlin models [34]. The incompressibility condition is satisfied by the technique of penalty method. In addition, the required modifications to the original model of the Mooney–Rivlin polynomial are presented and two special cases of this model are described, namely the incompressible Neo-Hookean and two-parameter incompressible Mooney–Rivlin material models.

3.1 Standard Mooney–Rivlin polynomial models

Fully parameterized beam and plate elements have the following expression for the matrix of the position vector gradient:

$$\begin{aligned} \varvec{F}=\frac{\partial \varvec{r}}{\partial \varvec{x}} \end{aligned}$$
(8)

where \(\varvec{r}\) is defined by Eq. (2) and \(\varvec{x}=[x\ y\ z]^\mathrm{T}\).

The strain energy density function for the isotropic materials depends only on the invariants of the right Cauchy–Green deformation tensor \(\varvec{C}_{r}=\varvec{F}^\mathrm{T}\varvec{F}\), defined as:

$$\begin{aligned} I_{1}= & {} \mathrm {tr}(\varvec{C}_{r}) \end{aligned}$$
(9)
$$\begin{aligned} I_{2}= & {} \frac{1}{2}\left[ \mathrm {tr}(\varvec{C}_{r})^{2}-\mathrm {tr}(\varvec{C}_{r}^{2})\right] \end{aligned}$$
(10)
$$\begin{aligned} I_{3}= & {} \det (\varvec{C}_{r}) \end{aligned}$$
(11)

If the material is incompressible, the constraint \(J=1\), where \(J=\det (\varvec{F})\), or identically \(I_{3}=J^{2}=1\), needs to be ensure through the material. Therefore, the strain energy density function for the incompressible materials depends only on the \(I_{1}\) and \(I_{2}\). A general form of the strain energy function for incompressible rubber-like materials in the form of the Mooney–Rivlin models is following:

$$\begin{aligned} U_{s}=\sum _{i=0}^{\infty }\sum _{j=0}^{\infty }\mu _{ij}(I_{1}-3)^{i}(I_{2}-3)^{j},\ \ \mu _{00}=0 \end{aligned}$$
(12)

where coefficients \(\mu _{ij}\) are material constants, usually determined from the experiment.

Material models from Eq. (12) assume that material is fully incompressible, that is, J is equal to one. This condition must be fulfilled by using additional techniques, such as Lagrange multipliers or penalty method [23]. In the current paper, the penalty method is used, due to simplicity and efficiency. In the penalty method, additional terms are included into the expression of the strain energy function:

$$\begin{aligned} U_{p}=\frac{1}{2}k(J-1)^{2} \end{aligned}$$
(13)

where \(U_{p}\) is volumetric energy penalty function and k is penalty coefficient, which must be selected large enough to ensure the incompressibility. However, k must not be too large because this can cause serious numerical issues. It should be pointed out that the use of the volumetric energy function \(U_{p}\) causes material to be considered as nearly, not fully, incompressible. In addition, penalty term k describes a real material property, the bulk modulus [6].

Finally, the strain energy density function of the incompressible Mooney–Rivlin models is written as:

$$\begin{aligned} U_{sic}=U_{s}+U_{p} \end{aligned}$$
(14)

In the current study, two models based on that representation are presented.

3.2 Modified Mooney–Rivlin models

The Mooney–Rivlin models expressed by Eq. (12) accounts for the volumetric (pressure) behavior, which may cause numerical problems during computations. Therefore, it is necessary to separate the volumetric and deviatoric components of the deformation [6]. To uncouple this behavior, instead of the invariants of the right Cauchy–Green deformation tensor, the invariants of the deviatoric part of that tensor are used:

$$\begin{aligned} \bar{I}_{1}=J^{-2/3}I_{1} \end{aligned}$$
(15)
$$\begin{aligned} \bar{I}_{2}=J^{-4/3}I_{2} \end{aligned}$$
(16)

The final form of the Mooney–Rivlin models for incompressible and hyperelastic materials used in that paper is given by substitution of the invariants from Eqs. (15) and (16) into Eq. (14):

$$\begin{aligned} U_{sic}= & {} U_{s}+U_{p}=\sum _{i=0}^{\infty }\sum _{j=0}^{\infty }\mu _{ij}(\bar{I}_{1}-3)^{i}(\bar{I}_{2}-3)^{j}\nonumber \\&+\,\frac{1}{2}k(J-1)^{2},\ \ \mu _{00}=0 \end{aligned}$$
(17)

The value of the elastic forces \(\varvec{Q}_{s}\) may be obtained according to Eq. (6), after integration of the above expression over the volume of the flexible body.

3.3 Incompressible Neo-Hookean material

The simplest material model prescribed by Mooney–Rivlin models is the incompressible Neo-Hookean. In that model, only the coefficient \(\mu _{10}\) does not equal zero, so the strain energy function has the following form:

$$\begin{aligned} U_{nh}=\mu _{10}(\bar{I}_{1}-3) \end{aligned}$$
(18)

where \(\mu _{10}=0.5\mu \), for \(\mu \) initially equal to the shear modulus.

3.4 Incompressible Mooney–Rivlin material

Another simple and commonly used material model is two-parameter incompressible Mooney–Rivlin material, which depends on two elastic coefficients, \(\mu _{10}\) and \(\mu _{01}\). Expression for the strain energy function of the Mooney–Rivlin material may be written as follows:

$$\begin{aligned} U_{mr}=\mu _{10}(\bar{I}_{1}-3)+\mu _{01}(\bar{I}_{2}-3) \end{aligned}$$
(19)

One can notice for \(\mu _{01}=0\) this model reduces to previously introduced Neo-Hookean material. In addition, if small strains are considered, the shear modulus is \(\mu =2(\mu _{10}+\mu _{01})\) and the Young’s modulus is \(E=6(\mu _{10}+\mu _{01})\) [4].

4 Locking elimination methods for incompressible materials

It was reported in the literature [13] that the locking phenomena affect many ANCF elements with linear elastic material model. The locking behavior is especially noticeable when continuum mechanics approach is used and revealed that flexible bodies have an erroneously stiffer characteristic during bending. However, in the case of the incompressible material model, the locking influence can be much larger than in case of the compressible materials, as will be shown in the current section.

4.1 Problem formulation

The use of hyperelastic, nearly incompressible material models [with strain energy density function given by Eqs. (18) or (19)] may lead to results that do not converge to the correct solution. This problem is denoted in the literature as volumetric locking phenomena [6, 34]. The following example presents this behavior.

As an illustration of the problem, the highly flexible, clamped beam under gravity forces is analyzed. The beam, as shown on Fig. 1, is of 1 m long and has a square cross section of 20 mm in size. The model uses six fully parameterized ANCF three-dimensional beam elements BE24 and includes two material models: compressible Neo-Hookean [23] and incompressible Neo-Hookean [see Eq. (18)]. The values of the Lamé parameters, used for the compressible Neo-Hookean, are \(\uplambda =1.15\) MPa as the first Lamé parameter and \(\mu =0.769\) MPa as the second Lamé parameter, called the shear modulus. Selected values of the parameters allow for large deformations. For the compressible Neo-Hookean, the material constant is \(\mu _{10}=0.333\) MPa and penalty constant is \(k=10^{9}\) Pa. The material density is assumed to be 7200 kg/m\(^{3}\).

Fig. 1
figure 1

Very flexible clamped beam under gravity forces

Full size image

In the Fig. 2, one can see the beam modeled with incompressible Neo-Hookean material, which is described in Sect. 3.3, exhibits a very stiff response—the tip displacements for the compressible model are much larger. This behavior is caused by the volumetric locking phenomena, which affects two considered elements: BE24 and BE30S as they do not include a trapezoidal cross-sectional deformation mode [25].

Fig. 2
figure 2

Locking influence on the incompressible material model. (

figure a
) Compressible Neo-Hookean, (
figure b
) Incompressible Neo-Hookean

Full size image

The presented example indicates that in simulations involving bending, the locking effect should be suppressed or the relevant element that is not affected by the locking should be considered. The paper discusses two methods that can be used for locking elimination in connection with incompressible material models. One is based on the selective reduced integration [11], and the second is based on the F-bar method [8].

4.2 Selective reduced integration

The integration of the integrals presented in Sect. 2, in particular the elastic forces vector, is usually performed numerically. Moreover, in the finite element analysis, usually the Gaussian quadrature is used to integrate over the cuboid element volume [4, 34]:

$$\begin{aligned}&\int _{V}G(\xi ,\eta ,\zeta )dV\nonumber \\&=\int _{0}^{1}\int _{\eta _{min}}^{\eta _{max}}\int _{\zeta _{min}}^{\zeta _{max}}G(\xi ,\eta ,\zeta )\hbox {d}\xi \hbox {d}\eta \hbox {d}\zeta \nonumber \\&\quad \approx \sum _{i=1}^{n_{\xi }}\sum _{j=1}^{n_{\eta }}\sum _{k=1}^{n_{\zeta }}w_{i}w_{j}w_{k}G(\xi _{i},\eta _{j},\zeta _{k}) \end{aligned}$$
(20)

where G is any scalar, vector or matrix function of spatial coordinates \(\xi \), \(\eta \) and \(\zeta \), \(w_{i}\), \(w_{j}\) and \(w_{k}\) are the weights of the quadrature, \(\xi _{i}\), \(\eta _{j}\) and \(\zeta _{k}\) are the points of the quadrature and \(n_{\xi }\), \(n_{\eta }\) and \(n_{\zeta }\) are the orders of the quadrature in directions, respectively, \(\xi \), \(\eta \) and \(\zeta \). The number of integration points in each direction might be different. In general, the quadrature of the function G is only an approximation of its integral. However, when function G is given as a polynomial (which is usually a case for the ANCF elements), the orders of the quadrature \(n_{\xi }\), \(n_{\eta }\) and \(n_{\zeta }\) can be selected in the way that will guarantee the exact evaluation of the integral (that order depends on the degrees of the polynomial in appropriate directions). Detailed information about the required number of integration points, types of quadratures and dependencies between the polynomial order and the quadrature order can be found in the literature [12].

In the case of reduced integration, the quadrature orders are selected as lower than those required for full integration. In the selective reduced integration, only some parts of the spatial function G are evaluated by means of the reduced integration.

The selective reduced integration technique is used to prevent Poisson locking in the ANCF elements with a continuum mechanics approach and linear elastic material model [11], as well as for many FEA elements [47]. This method’s idea is to split the integral of the strain energy function into two parts. The first part, which does not take into account the Poisson effect, is fully integrated. The second part, which considers the Poisson effect, uses a reduced integration scheme.

The expression for the strain energy density function of the incompressible Mooney–Rivlin models, given by Eq. (17), can also be split adequately. The first part, denoted as \(U_{s}\), does not account for the volumetric behavior; therefore, it is fully integrated. The second part, the volumetric energy penalty function \(U_{p}\), accounts for the volumetric effect and is only considered at the beam axis. The strain energy density value can be expressed by the following formula:

$$\begin{aligned} U_{sic}^{SRI}=\int _{V}U_{s}\hbox {d}V+A\int _{l}U_{p}\hbox {d}l \end{aligned}$$
(21)

where index SRI denotes the selective reduced integration scheme, while A and l are, respectively, the area and the length of the element, which indicates that function \(U_{p}\) is evaluated only along beam centerline. This also means that in the formula for quadrature given by Eq. (20), the orders \(n_{\eta }\) and \(n_{\zeta }\) are equal to one.

4.3 F-bar method

The method of selective reduced integration is suitable for the problems where it is easy to split the energy density function into deviatoric and volumetric parts. In the classical finite element method, the strain projection B-bar approach is treated as an extension of the reduced and selective integration techniques [17]. In turn, the generalization of B-bar method to finite strain formulation is so-called F-bar method, introduced by de Souza Neto et al. [39], which involves the product decomposition of the position vector gradient into volumetric and deviatoric parts [6, 8].

The matrix \(\varvec{F}\) may be multiplicative split into its deviatoric (volume preserving) and volumetric-dilatational parts:

$$\begin{aligned} \varvec{F}=\varvec{F}^{dil}\varvec{F}^{dev},\quad \varvec{F}^{dev}=J^{-1/3}\varvec{F},\quad \varvec{F}^{dil}=J^{1/3}\varvec{I}\nonumber \\ \end{aligned}$$
(22)

where \(\varvec{I}\) is identity matrix. Individual parts of this split have the following properties:

$$\begin{aligned} \det (\varvec{F})=\det (\varvec{F}^{dil})=J,\quad \det (\varvec{F}^{dev})=1 \end{aligned}$$
(23)

Now we can define a modified matrix \(\varvec{F}\) in terms of the deviatoric part and the modified dilatational part as:

$$\begin{aligned} \bar{\varvec{F}}=\bar{\varvec{F}}^{dil}\varvec{F}^{dev} \end{aligned}$$
(24)

where \(\bar{\varvec{F}}^{dil}=\overline{J^{1/3}}\varvec{I}\) and \(\overline{J^{1/3}}\) may be given by:

$$\begin{aligned} \overline{J^{1/3}}=\pi (J^{1/3}) \end{aligned}$$
(25)

where \(\pi \) is the linear projection operator described below. Finally, one can write:

$$\begin{aligned} \bar{\varvec{F}}=\frac{\overline{J^{1/3}}}{J^{1/3}}\varvec{F} \end{aligned}$$
(26)

The value of the \(\bar{\varvec{F}}\) is directly used in the expression for the volumetric penalty function. There is no need to include the factors derived from \(\varvec{F}\) in the \(U_{s}\) part of the energy density function given by Eq. (17), as that part is volume preserving.

4.3.1 Projection operator

The use of the strain projection method requires further investigations on the choice of the projection operator and the space on which the projection will be performed. In the original F-bar method, the modified dilatational part of the position vector gradient is evaluated at the element center only [39]; however, it was applied to low-order elements and does not provide satisfactory results for fully parameterized ANCF beams. As suggested by Elguedj et al. [8], the \(L^{2}\) projection of strains is appropriate in the case of the higher-order elements. Moreover, the associated space on which the projection is performed should have a lower order and then the original element approximation given in Sect. 2. It turns out that the proper space for the standard fully parameterized spatial beam element should be constant in the transversal direction. Therefore, one can define the lower-order basis using tilde notation as:

$$\begin{aligned} \tilde{\varvec{S}}=\left[ 1-3\xi ^{2}+2\xi ^{3}, l(\xi -2\xi ^{2}+\xi ^{3}),3\xi ^{2}-2\xi ^{3}, l(-\xi ^{2}+\xi ^{3})\right] \nonumber \\ \end{aligned}$$
(27)

and the Eq. (25) may be written in the new space as [8, 17]:

$$\begin{aligned} \overline{J^{1/3}}=\tilde{\varvec{S}}\tilde{\varvec{J}}^{1/3} \end{aligned}$$
(28)

where:

$$\begin{aligned} \tilde{\varvec{J}}^{1/3}=\tilde{\varvec{M}}^{-1}\int _{V}\tilde{\varvec{S}}^\mathrm{T}J^{1/3}\hbox {d}V \end{aligned}$$
(29)

where \(\tilde{\varvec{M}}\) is the “mass” matrix of the tilde basis:

$$\begin{aligned} \tilde{\varvec{M}}=\int _{V}\tilde{\varvec{S}}^\mathrm{T}\tilde{\varvec{S}}\hbox {d}V \end{aligned}$$
(30)

The presented procedure corresponds to \(L^{2}\) projection of the \(J^{1/3}\) into the \(\tilde{\varvec{S}}\) and can be used with both element types affected by the volumetric locking phenomena, namely BE24 and BE30S.

4.3.2 Implementation details

The modified position vector gradient \(\bar{\varvec{F}}\) should be used to evaluate the strain energy density function:

$$\begin{aligned} U_{sic}^{F{\text {-}}bar}=\int _{V}U_{s}\hbox {d}V+A\int _{l}\bar{U}_{p}\hbox {d}l \end{aligned}$$
(31)

where index F-bar denotes the strain projection method and the modified volumetric energy penalty function \(\bar{U}_{p}\) may be calculated by substitution of the Eq. (26) into (13) and using the properties (23) as:

$$\begin{aligned} \bar{U}_{p}=\frac{1}{2}k\left[ \left( \overline{J^{1/3}}\right) ^{3}-1\right] ^{2} \end{aligned}$$
(32)

Finally:

$$\begin{aligned} \bar{U}_{p}=\frac{1}{2}k\left[ \left( \tilde{\varvec{S}}\tilde{\varvec{M}}^{-1}\int _{V}\tilde{\varvec{S}}^\mathrm{T}J^{1/3}\hbox {d}V\right) ^{3}-1\right] ^{2} \end{aligned}$$
(33)

After the evaluation of the \(\tilde{\varvec{J}}^{1/3}\) vector, the value of the \(\overline{J^{1/3}}\) depends only on the longitudinal coordinate; therefore, in Eq. (31), the \(\bar{U}_{p}\) is evaluated only along the beam centerline.

5 Numerical examples

Presented methods of preventing volumetric locking are applied to the fully parameterized spatial ANCF beam elements BE24 and BE30S. Elements BE30T and BE42 include the trapezoidal deformation mode of the cross section that suppresses the occurrence of the volumetric locking in these elements; therefore, the methods from the previous section are not applied to them. In the following examples, the efficiency of the absolute nodal coordinate formulation for modeling the nearly incompressible material models is studied using three simple models of highly deformable beams.

5.1 Physical pendulum

A flexible pendulum under gravity forces is employed for the purpose of numerical verification of the volumetric elimination methods and to shed some light on issues associated with elements BE24 and BE30S. A similar pendulum was examined in [23]. In the undeformed state, beam has a length of \(1\,\mathrm {m}\) and a uniform square cross section of dimension \(20\,\mathrm {mm}\). The density of the pendulum is assumed to be \(7200\,\mathrm {kg/m^{2}}\). The material model is highly flexible, to allow large deformations of the body. In this example, the modulus of elasticity is equal to \(6\,\mathrm {MPa}\) and Poisson ration is equal to 0.3 (for compressible material model). In the case of the Mooney–Rivlin material model, the values of its coefficients are \(\mu _{10}=0.8\,\mathrm {MPa}\) and \(\mu _{01}=0.2\,\mathrm {MPa}\). For the incompressible materials, the penalty method is used to ensure incompressibility with penalty term assumed to be \(k=10^{9}\,\mathrm {Pa}\).

The pendulum at one end is connected to the ground by a spherical join, as shown in Fig. 3. The base of the pendulum is subjected to prescribed motion \(u_{X}=-0.02\sin (2\pi t)\) m, where \(u_{X}\) denotes the displacements of the base in the X direction, and t is the time expressed in seconds. That additional base motion increases body deformations.

Fig. 3
figure 3

Very flexible physical pendulum

Full size image

Pendulum body is built with four fully parameterized ANCF BE24 elements. It should be pointed out that the behavior of the element BE30S is the same like the BE24; therefore, for clarity only, the BE24 is presented in this section. However, the results and conclusions are valid for both element types. The performed dynamical simulations’ results are shown in Figs. 4 and 5. In the current example, the following material models are used: compressible Neo-Hookean [34], incompressible Neo-Hookean with and without F-bar projection method and incompressible Mooney–Rivlin with and without F-bar projection. In addition, to present the locking influence more clearly, the results for the models with incompressible materials without locking prevention methods are compared with result of the rigid pendulum. Once again, the results are presented only for the F-bar method for clarity and for the numerical verification. The behavior of the selective reduced integration is analogous in the following examples to the F-bar projection, and the conclusions are valid for both locking elimination techniques.

Fig. 4
figure 4

Vertical position of the beam tip. (

figure c
) Compressible Neo-Hookean, (
figure d
) Incompressible Neo-Hookean with F-bar projection, (
figure e
) Incompressible Mooney–Rivlin with F-bar projection

Full size image
Fig. 5
figure 5

Vertical position of the beam tip. (

figure f
) Incompressible Neo-Hookean with F-bar projection, (
figure g
) Incompressible Neo-Hookean, (
figure h
) Rigid

Full size image

Figures 4 and 5 show the pendulum’s free end-tip displacements. All the results in Fig. 4 show large deformations of the pendulum. The results for the compressible and incompressible models are very similar, and they are barely distinguishable at this resolution. This indicates influence of the volumetric locking in incompressible models has been suppressed, since no over-stiff behavior may be noticed during bending, and that models give reasonable results.

The plot in Fig. 5 shows comparisons between results of the incompressible models with locking prevention by means of the F-bar method, incompressible models without any locking elimination technique and the rigid model. It should be pointed out that the results of the incompressible models without F-bar projection are very close to the rigid body solution. At the same time, difference in results of models with and without F-bar projection applied is significant. One should conclude volumetric locking influence on incompressible hyperelastic material model results is enormous.

Figures 6 and 7 show the value of the determinant of the deformation gradient tensor \(J=\det (\varvec{F})\), which indicates the compressibility, as a function of time and position. Figure 6 shows the value of J at a fixed point on the body (point A in Fig. 3) in time. In turn, Fig. 7 shows the value of the J across beam thickness—from point A to B (see Fig. 3) for given time (\(t=0.2\) s). One can notice that material is fully incompressible in whole volume only for the original formulation without a locking elimination technique. For the compressible material models and the incompressible material with F-bar projection, the change of the value of J is noticeable. One can also note that for incompressible material with F-bar projection, J is distributed symmetrically through beam thickness as expected (the projection lower-order basis is considered only at the beam centerline).

Fig. 6
figure 6

Value of the \(\det (\varvec{F})\) in time. (

figure i
) Compressible Neo-Hookean, (
figure j
) Incompressible Neo-Hookean with F-bar projection, (
figure k
) Incompressible Neo-Hookean

Full size image
Fig. 7
figure 7

Value of the \(\det (\varvec{F})\) along beam thickness. (

figure l
) Compressible Neo-Hookean, (
figure m
) Incompressible Neo-Hookean with F-bar projection, (
figure n
) Incompressible Neo-Hookean

Full size image

5.2 Static simulation of the clamped beam

In this example, the clamped beam under gravity forces is analyzed, as shown in Fig. 8. The flexible body has the same material properties and dimensions like in previous example, and only the cross section size is changed—the height and width of the beam are increased to 40 mm. Pendulum is rigidly attached to the ground without any prescribed motion. The only external force that is acting on the system is gravitational force. In this example, the results of the ANCF models are compared with the simulation performed with commercial finite element package [2]. It worth to note that the continuity between elements is imposed at the position and first-order gradient level (that is \(C^{1}\)continuity is imposed). Moreover, the clamped condition is also imposed at those coordinates. For a more detailed discussion on the continuity condition, see, e.g., [15, 28].

Fig. 8
figure 8

Very flexible clamped beam under gravity forces

Full size image

The main purpose of this example is to compare the results of different ANCF beam elements and locking elimination techniques presented in this paper. Simulation is performed with both incompressible material models: Neo-Hookean and Mooney–Rivlin. In addition, the results for both locking elimination methods, the selective reduced integration and F-bar projection, are presented. The solution obtained with FEA analysis is treated as the reference. Therefore, in the FEA package, a very fine mesh of 1600 higher-order solid elements (SOLID186) is employed.

Figures 9 and 10 show results of the static analysis. The error in tip positions is measured by \(\mathrm {sign}(d_{z})\left| \varvec{d}\right| \), where \(\left| \bullet \right| \) denotes quadratic norm, \(d_{z}\) is a z component of the vector \(\varvec{d}\), while \(\varvec{d}=\varvec{u}_{r}-\varvec{u}_{A}\) where \(\varvec{u}_{r}\) is reference tip displacement and \(\varvec{u}_{A}\) is the ANCF tip displacement. Solutions for elements BE24 and BE30S without locking elimination techniques are not shown as they produce an enormous error—irrespective of the element type, material model or number of elements, the error in tip position has a value of about 1.19.

Fig. 9
figure 9

Convergence analysis—Neo-Hookean material model. (

figure o
) BE42, (
figure p
) BE30T, (
figure q
) BE30S SRI, (
figure r
) BE30S F-bar, (
figure s
) BE24 SRI, (
figure t
) BE24 F-bar

Full size image
Fig. 10
figure 10

Convergence analysis—Mooney–Rivlin material model. (

figure u
) BE42, (
figure v
) BE30T, (
figure w
) BE30S SRI, (
figure x
) BE30S F-bar, (
figure y
) BE24 SRI, (
figure z
) BE24 F-bar

Full size image

Results for both material models have the same characteristics. Only models with BE42 converge to the reference solution with acceptable error, and the model with 40 BE42 elements predicts the solution correctly. However, the convergence rate of this element to the steady response is the slowest, which may be caused by the shear locking [28]. Element BE30T converges to wrong over-stiff solution with noticeable error, despite that this element includes a trapezoidal displacement mode. This is probably cause by the fact that when very large deformations occur due to bending, the cross section takes the form of trapezoid but with curved sides, which requires a quadratic transversal approximation that BE42 provides. Results also show that both locking elimination techniques, selective reduced integration and F-bar, converge to almost the same solution for both BE24 and BE30S elements. Results obtained with F-bar are slightly closer to the reference solution then those for the SRI or then results for the SRI, but the difference is very small. In addition, higher-order element BE30S shows a faster convergence rate then BE24. All this results show over-soft behavior as they predicts too large displacements. The observable diverse between the ANCF beam model and the reference results is mainly because of the limitations of the locking elimination techniques employed in this paper. It is shown in Fig. 7 that the value of the \(\det (\varvec{F})\) varies linearly across the beam thickness with exactly zero value only on the element centerline when locking elimination due to F-bar (as well as SRI) is employed. Therefore, the largest violation of the incompressibility assumption can be noticed at the beam surface, and it increases as the body deformations increase. Nevertheless, presented behavior is consistent with finite element characteristics, where incompressibility is ensured only in average sense [3]. One of the solutions of this issue might be the increase of the number of elements not only along the body length, but also along its height or width. This will result in body mesh composed of two or more element layers. Since the incompressibility might be preserved in each layer centerline, its violation should propagate slower across the body height. The preliminary studies show that the solution obtained with two layers of the elements along thickness is in much better agreement with reference solutions presented in Figs. 9 and 10. However, this approach requires the more detailed analysis.

5.3 Clamped rubber-like beam

The next example is the clamped beam model, which is made of rubber-like material and undergoes gravitational forces, similar to one shown in the literature [18]. The beam has a length of 0.35 m, and the rectangular cross section of 7 mm width and 5 mm height. The density of the material is assumed to be 2150 kg/m\(^{3}\), Kirchhoff modulus is equal to 1.91 MPa and the penalty term is equal to \(10^{9}\,\mathrm {Pa}\). The beam is composed of the number of beam elements that provide converged solution—it reflects the convergence rates discussed in Sect. 5.2. The results of the simulation obtained within the ANCF framework are compared with the results of FEA package [2]. In the FEA simulation, 336 SOLID185 [1] elements are used. In the classic FEA package, the beam elements cannot be used with nonlinear material models. Moreover, the mesh with smaller number of elements was insufficient. In this example, only the incompressible Neo-Hookean material model is considered.

Figures 11 and 12 show the results for FEA and different ANCF element formulations (element types and locking prevention techniques). A reasonably good agreement can be observed in Fig. 11 between results, and all ANCF models converge to almost the same solution. On the other hand, Fig. 12 shows a noticeable differences for BE30T in the displacement characteristic and very good agreement for BE42 when compared with reference solution. Notably, in the given application, the ANCF model involves more than four times less parameters than the FEA solid model when elements BE24, BE30S and BE30T are considered. However, for BE42 model that employs 50 elements, the number of parameters is comparable between the ANCF and FEA models (1650 for FEA and 1512 for ANCF). Nevertheless, the numerical examples show that the elimination of the volumetric effect leads to considerable improvement and can be compared to the classical finite element formulation solution. Moreover, it can be expected that the ANCF model build with two BE24 or BE30S element layers should result is greater agreement with reference FEA solution.

Fig. 11
figure 11

Beam end-tip displacements. (

figure aa
) BE24 with F-bar—30 elements, (
figure ab
) BE24 with SRI—30 elements, (
figure ac
) BE30S with F-bar—20 elements, (
figure ad
) BE30S with SRI—20 elements, (
figure ae
) reference

Full size image
Fig. 12
figure 12

Beam end-tip displacements. (

figure af
) BE30T—20 elements, (
figure ag
) BE42—50 elements, (
figure ah
) reference

Full size image

6 Conclusions

The absolute nodal coordinate formulation is known to be well suited for large deformation analysis of the flexible bodies. The present paper shows ANCF can be used with a variety of nonlinear material models in a straightforward manner, including the incompressible ones. The original implementation of the nonlinear incompressible materials shows very stiff behavior in bending-dominated tests, which suggest a lockings occurrence.

In order to perform efficient and accurate model simulations with incompressibility, four fully parameterized beam elements and two locking elimination techniques are considered. Two higher-order beam elements, BE30T and BE42, include in their displacement field the trapezoidal deformation mode, which aim is to account for the Poisson deformation mode under bending. Therefore, BE30T and BE42 kinematic description should alleviate volumetric locking phenomena. For the elements BE24 and BE30S, the influence of the volumetric locking is enormous when nearly incompressible material models are used. Hence, those elements should be used in connection with given locking elimination techniques. Simulations results acquire within the ANCF framework are compared with a commercial FEA package, and reasonable agreement may be observed. However, it should be pointed out that further development should be carried out to increase the convergence rate of the BE42 element and to improve the prediction of results when locking elimination techniques are employed.

For slender structures made of materials with nonlinear characteristics, ANCF models can be usually built with much less parameters than classic FEA models because the ANCF beam finite elements (fully parameterized) can be used, instead of the solid elements in the FEA. However, the ANCF beam elements may currently model only simple cross-sectional shapes [27]. Despite that, when practical applications are considered, such as [29], the computation cost of the method should be considered. One may try to parallelize the ANCF method by using existing approaches that are originally developed for multi-rigid-body dynamics [9, 20]. The extensions of the methods for highly flexible body systems built upon a divide-and-conquer framework can be found in [19].

As a direction of future research, it is desirable to apply the presented locking elimination methods to different material models (e.g., Ogden hyperelastic model) and to shell and solid ANCF elements. Furthermore, an alternative methods that allow for locking eliminations should be developed, and the layered ANCF beam elements assembly procedure should be formulated and thoroughly examined. Moreover, the literature studies suggest that finite element models that use incompressible materials might provide correct solution for displacements and highly inaccurate solution for stresses [4]. Therefore, the stress distribution in such models should also be concerned as important research subject.