Performance review of locking alleviation methods for continuum ANCF beam elements

The absolute nodal coordinate formulation (ANCF) is a nonlinear finite element approach proposed for the large deformation dynamics analysis of beam- and plate/shell-type structures. In the ANCF approach, elastic forces can be defined using three-dimensional elasticity-based continuum mechanics. This approach is often straightforward, and it makes it possible to use advanced material models in the ANCF framework. However, it has been pointed out in several studies that continuum ANCF-based elements with a full three-dimensional elasticity description can suffer from locking phenomena. In this study, a comparison between various combinations of locking alleviation techniques and their applicability to different ANCF beam variants is studied using numerical examples. Furthermore, the enhanced deformation gradient (EDG) technique, which has been proposed recently in finite element literature, is demonstrated for high-order ANCF beam elements. Based on the numerical tests, none of the currently available techniques are suitable for all types of ANCF elements. The paper also shows that the efficiency and accuracy of the techniques are case-dependent. For the ANCF beam element involving higher-order terms with respect to trapezoidal mode, however, the EDG-based techniques are preferable to reduce locking phenomena.


Introduction
The absolute nodal coordinate formulation (ANCF) is a nonlinear finite element approach proposed by Shabana in [28] for beam-and plate/shell-type structures undergoing large deformation dynamics. Today, the ANCF elements are often used within the multibody community to solve various two-and three-dimensional deformation problems [13]. The key idea of ANCF elements is to use absolute nodal positions and their gradients in the kinematic description. Similar approaches were proposed by Rhim and Lee [27] (termed the vectorial approach) for nonlinear static problems involving finite elements and, recently, in the work [4].
The gradient vectors, the slope vectors in the transverse directions, are used to define cross-sectional deformations. They make it possible to use threedimensional elasticity descriptions for the ANCF elements. However, ANCF elements based on full threedimensional elasticity are often unable to rightfully represent the deformation caused by the coupling of axial and transverse normal strains [31]. As a result, an inac-curate response results when Poisson's ratio is not zero. This problem is often referred to as the Poisson locking effect [7,24,27] etc.
The theoretical example of that is given in [18] and [32] as follows. Pure bending deformation leads to the trapezoidal shape of an initially rectangular cross section. However, some ANCF elements cannot reproduce the trapezoidal shape, which leads to additional stresses during the bending deformations. As a result, it produces an inaccurate response when Poisson's ratio differs from zero.
Additionally, as shown in [10], ANCF-based elements can suffer from other locking phenomena such as curvature thickness locking and shear locking. The first can be attributed to the element kinematics volume description and related to elements with the linear cross-section interpolation description. That interpolation leads to the norm reduction in the vectors, which define the cross-sectional thickness, when the ending cross sections are not parallel, and the shrinking of the intermediate cross sections. Another locking problem is shear locking, which happens with a linearly varying bending moment. In this case, some higher-order elements demonstrate quadratic shear strain distribution, which is not correct. As a result, the shear strain component will be assigned an inflated energy value, leading to smaller displacement predictions. These phenomena might influence results and require separate, careful studies. Therefore, they have not been considered here. Moreover, beam elements can also demonstrate stiffer behavior, which is attributed to geometric nonlinearities and can be wrongly interpreted as locking [8]. That aspect has been also omitted here.
The higher-order ANCF elements proposed in [15,16,29] can be used to reduce locking phenomena due to the limited cross-section deformation description of the ANCF element. The idea of a higher-order ANCF is to enrich the polynomial basis with higher-order polynomials by using the components of higher-order derivatives as additional nodal coordinates to represent more deformation modes and by describing the cross-section deformation more precisely.
Another option is to use one of the locking alleviation techniques. So far, several methods have been proposed such as the Enhanced Assumed Strain (EAS) method [30], the Enhanced Continuum Mechanics formulation (ECM) [12,17,22], and the Strain Split Method (SSM) [24]. The idea of the ECM approach lies in the concept of reduced integration. In this approach, the matrix of elastic coefficients is divided into two parts so that the first part presents Poisson coupling between the normal strains and the second part only involves thickness deformation. The strain split method is similar to the ECM, also considering the constitutive model split. However, in the case of the ECM, the same strain matrix is used for full and reduced integration parts. On the contrary, in the SSM, the different strain matrices are used. However, the main limitation of the mentioned approaches is applicable to only the Kirchhoff-Saint-Venant material model, which significantly limits their usage. The EAS approach was first proposed by Simo and Rifai [30]. The basic idea was to enrich the element strain field to improve its nonuniform strain conditions. The main advantage of the EAS is that it can be used in all material models. However, it increases the degrees of freedom of the elements.
In the study, the usability of locking alleviation techniques for lower-order, fully parameterized and higher-order continuum-based ANCF beam elements is demonstrated via various benchmark problems. Furthermore, the enhanced deformation gradient method (EDG), recently proposed by Pfefferkorn and Betsch [25], is implemented in the ANCF framework and compared to other techniques. The EDG is similar to the EAS but in the EDG, the strain field is enriched via deformation gradients instead of via direct modification of the strain field. Therefore, it possesses the same advantages as the EAS approach. All approaches with the additional requirements imposed on them for representative purpose are collected in Table 1 below.
The introduced ANCF elements and corresponding locking alleviation methods are verified against analytical or experimental results, and their convergence speeds are compared. As a result, the most suitable compositions, those that alleviate locking and demonstrate solutions closest to analytical or experimental results, should be found by testing various combinations of the locking alleviation techniques and the ANCF elements.

Kinematics of ANCF beam elements
This section explains, the kinematics of the spatial ANCF beam elements. Let r = r (x, y, z) be the position vector field at the current configuration, and r the position vector in the initial configuration being. The  connection between these vectors is: where vector u h a displacement vector. The paper considers three different elements: gradient deficient 3333 [22], fully parametrized 3243, and high-order 3343 [7]. These four numbers in the element name abcd have the following meaning: a is the dimension of the element, b denotes the number of nodes, c is the number of vectors used in the approximations, and d is the polynomial basis used to approximate all three dimensions [7]. The ANCF 3243 element has been under extensive investigation, and the Poisson locking phenomenon is known and well-documented in the literature. [11,20,31] are examples. Moreover, this problem has been discussed earlier in the work [27], where a similar vectorial approach has been presented. The ANCF 3333 element also suffers from Poisson locking, which was the main problem with this formulation [21,22]. In the case of bending, trapezoidal deformation of rectangular cross sections is unavailable, because this deformation mode is not included as one of the subject ANCF element shape functions. This leads to stiff behavior, which is the cause of the Poisson locking. The ANCF 3343 element was introduced in [7], where it was also shown that it suffers from severe Poisson locking. Figures 1  and 2 show the kinematics of the subject elements. The nodal degrees of freedom used, and the polynomial basis are also summarized in Table 2. The directional derivatives are defined as r α = ∂r ∂α , α = {x, y, z}.
The elements are isoparametric; therefore, they can be described with the local bi-normalized coordinate system ξ = {ξ, η, ζ } presented in Fig. 1, with the ranges for the local coordinates being [−1, 1].
As an example of the three-node gradient deficient beam 3333, the vector q of nodal degrees of freedom can be written as: The shape function matrix N m for this element takes the form: where I is a 3 × 3 identity matrix and where l x , l y , and l z are the physical dimensions of the element. Therefore, the finite elements can express the body motion in the form:

Equations of motion
In this section, the equations of motion for the used elements are derived based on the principle of virtual work.
where W ext is the virtual work by external forces, W elast is the virtual work by elastic forces, and W inert is the virtual work by inertial forces. The variation of virtual work by inertial forces with respect to the nodal coordinates can be written as follows.
where ρ is the mass density, and V is the volume of the element in the reference configuration. The integral expression is usually called the mass matrix, which has a constant meaning from Eq. (2).
The variation of external forces takes the following form: where vector b is the vector of body forces.

Virtual work of elastic forces
The following paragraphs present, in detail, the elastic forces for each alleviation method.

Standard continuum-based method
Going forward, the general continuum mechanics based approach with full three-dimensional elasticity will be referred to "Cont." This description suffers from Poisson locking. Modifications made to strain energy to avoid locking will be presented in the following subsections. The variation of W elast with respect to the nodal coordinates is where S is the second Piola-Kirchhoff stress, and E is the Green-Lagrange strain, which is given as From (3), the deformation gradient F is where J is the Jacobian matrix providing the transformation between the physical and local coordinate systems. The expression for the second Piola-Kirchhoff stress tensor is of the form where λ and G are Lame elastic coefficients, which relate to Young's modulus E and Poisson's coefficient ν as The expression (8) can also be represented via Voigt notation, where stresses and strains have the following vector forms.
The stress-strain relation (11) can be rewritten as where the elasticity matrix D is defined as with shear correction factors k s 2 , k s 3 . Then, the expression (8) has the form Recalling that the element is isoparametric, An advantage of the expression (15) is that integration can be performed using the Gaussian integration; furthermore, this procedure can be applied directly and independently for the three integrations.

Enhanced continuum method-ECM
The elements based on the method presented in Subsect. 4.1 suffer from Poisson locking, except for some high-order ones, where the rich polynomial basis alleviates the effect naturally [7]. The ECM offers lockingfree solutions based on splitting of the elasticity matrix D. This approach was first suggested for ANCF element in [12].
where D v is responsible for the Poisson effect only In Eq. (16), diagonal matrix D 0 is as The strain energy variation following (16) can be also split into parts: where A = l y l z is the cross-sectional area. This method can be considered as selective reduced integration [22], and from Eq. (19), it is clear that the Poisson effect is considered only on the beam axis {η, ζ } = {0, 0}.

Enhanced assumed strain method-EAS
Simo and Rifai [30] proposed the enhanced assumed strain approach as a generalization of the method of incompatible modes. The key idea behind this method is improving the performance of the element with additional variables of strains via improvement of the nonuniform strain conditions. The enhanced strain field can be written in the following form: where E com is the compatible strain field obtained from (9) and E enh is the enhanced strain field, which is defined per element. The enhanced strain field must be L 2 orthogonal to the discrete stress This eliminates the independent stress. Then, the variation of the Hu-Washizu functional type in the reference configuration, assuming that α is a vector of additional enhanced strain variables, gives the following system [24]: The solution operations are performed as follows. From the first equation of (20), the problem is solved for q. Then, once q is defined, the additional enhanced strain variables α are updated from the second equation of (20). The only part left here is the discretization of E enh field. It is assumed that M m is a shape function matrix. To alleviate the Poisson locking in this work, the linear enhanced interpolation in the element transverse directions is considered as the following.
Additionally, the second M m matrix is assumed to be: In the numerical experiments, EAS methods method solutions that performed better using (22) are referred to as E AS − 2. The solutions that performed better using (23) are referred to as E AS − 4.

Enhanced deformation gradient method-EDG
The enhanced deformation gradient method was previously introduced by Pfefferkorn and Betsch [25]. It is similar to the EAS. The main difference is that the enhanced field is not added to the strain field directly, but via enrichment of the deformation gradient F as: where F com is the compatible part of the deformation gradient obtained in (10). On the other hand, the presentation of F enh might have different forms, various transformations are considered in [25]. In the case of the ANCF elements, however, they lead to approximately the same results. In this work, the form presented in [26] with small modification will be used as: Additionally, as given in Sect. 4.3, the second version of the F com with four variables is as follows: Using the referencing system from Sect. 4.3, they are designated as EDG-2 and EDG-4, respectively. Again, as in EAS approach, the system for the derivation of δW elast has the following form:

Strain split method-SSM
The strain split method was first presented and explained in detail in [24]. One of its main advantages is that it can cure locking without introducing the additional variables described in Sects. 4.3 and 4.4. The key concept of this technique is that-the axial and transverse strains of high-order terms are split via decomposition of the Green-Lagrange strain tensor (9) and the constitutive law (13). According to [24,33], the gradient deformation tensor (10) can be separated as where with r c = r| η,ζ =0 and Then, with help of (28) and (29), Green-Lagrange strain can be written as where Using Voigt notation, the second Piola-Kirchhoff stress can be defined as . Finally, the expression for the virtual work of elastic force takes the form

Results of bending tests
The following paragraphs analyze the performance of the techniques described in Sect. 4. Several often gradient deficient, fully parameterized, and high-order ANCF elements are used. Their initial models, based on standard continuum approach, suffer from Poisson locking. They are referred to in the tables as "Cont." Several benchmark problems related to the bending, presented in Fig. 3, are considered. In these tests, the boundary condition is defined by fixing all the nodal degrees of freedom at the clamped end.
Additionally, ANCF element 3363 is considered based on the conclusions from [7]. This element demonstrated reasonable results during deformation tests and maintained usable convergence rates.

Small displacement cantilever problem
The beam structure under investigation has a rectangular cross section with width W = 0.1, height H = 0.5, and length L = 2. E = 2.07 × 10 11 N m 2 and Poisson's ratio ν = 0.3. The small deformations case is tested first. To this end, the vertical tip load is set to F y = 62.5 × 10 3 N .
The convergence rate against absolute error comparison for all the above-mentioned elements and modifications is provided below. The absolute error is calculated in relation to the reference solutions presented in [7].        The tests of the beam subjected to small bending load revealed that best results for all element types are obtained using the ECM approach. See Tables 3, 4 and 5 and Figs. 4, 5 and 6. The ECM also provides the highest convergence rate among all approaches. However, the results presented by the other approaches for the high-order element agree within reasonable limits, with the reference solution (about 1% error) and even outperform the high-ordered ANCF 3363 element. Additionally, including 2 or 4 degrees of freedoms into the kinematics-based locking alleviation techniques results in minimal differences for small cantilever bending problems regardless of element type. Moreover, SSM performance deteriorates as the number of vectors per node increases.

Large displacement cantilever problem
The large displacement cantilever problem is analyzed in this section. The vertical tip load is set to F y = 62.5× 10 6 N . Other geometrical and physical characteristics are the same as in Sect. 5.1. Tables 6 and 7 show that for the large bending tests, the best results for the gradient-deficient and the fully parameterized elements are again provided by the ECM approach. Figures 7 and 8 also illustrate its high convergence rate. However, for this task, it does not perform as well as element 3363. In the case of the high-order element, from Table 8 as well as from Fig. 9, the EDG-4 gives the smallest error alongside the best convergence   Fig. 9 Convergence rate for the high-order element, ANCF 3343, in the large bending test rate and performs better than that provided by the 3363 ANCF element.

The Princeton beam experiment
The following paragraphs consider the Princeton beam experiment. This problem was originally conducted in [5], presented in [6], and already used within the ANCF framework by Bauchau et al. [2] and Ebel et al. [7]. In this study, beam length L = 0.508 m. Its rectangular cross section was set to height H = 12.377×10 −3 and width W = 3.2024 × 10 −3 m. For such a cross section, the torsional correction factor was set to k t = 0.198 [9]. All shear correction factors were set to 1 [7]. The physical characteristics were set as follows. Young's modulus E = 71.7 × 10 9 N m 2 and Poisson's ratio ν = 0.31. In this example, several tip loads, equal to F pb 1 = 8.896 and F pb 2 = 13.345 N, were analyzed. An angle β was defined as the angle between the force vector and the y-axis, meaning that with β = 0 • the force acted in the negative y-direction as shown in Fig. 3. The angle varied from β ∈ [0 • , 90 • ]. The results were compared against commercial software solutions using BEAM188 and Solid185 elements.
For the calculation, 32 ANCF elements were used. The higher number of the elements leads to similar results, also following from convergence speed tests shown in Figs. 4, 5, 6, 7, 8 and 9. The reference results for the ANSYS solutions were obtained using 128 beam elements (in the case of higher element number, see Appendix Appendix A1), and 500 × 12 × 4 for the model from the solid elements.
The displacements in z-and y-directions of the beam under F pb 1 tip load were considered. Figure 10a and b reveals that none of the modifications for the gradient deficient element fully passed the test. Bending along the z-axis was well-represented with the ECM and SSM approaches, with results similar to the ANSYS beam, ANCF 3363, and solid element solutions. Although the bending along the second axis was not represented well, both these modifications resulted in poorer performance than demonstrated by the ANSYS software and ANCF 3363 element (Fig. 11). Additionally, all the kinematic-based approaches, EDG and EAS, provided similar inaccurate results without differences to the number of additional DOFs.
Again, as for the gradient deficient element, in the case of the fully parameterized element, none of the modifications could correctly capture all deformations, as shown in Fig. 11. The ECM and SSM approaches gave appropriate displacement results only along the z-axis. Furthermore, the SSM approach provided solutions similar to those provided by the solid element, and the ECM gave results similar to those obtained from the ANSYS and ANCF 3363 beam elements.
For the high-order element, ANCF 3343, the best results were provided by EDG-4 and the ECM approaches, as shown in Fig. 12. The ECM element demonstrated better agreement with the ANSYS beam element solution, but EDG-4 demonstrated better agreement with the ANSYS 3D solution. Additionally, in the case of the EAS approach, there were no significant differences between EAS-2 and EAS-4. The displacements given under F pb 2 tip load were as follows.
Similarly to the previous loading case, Figs. 13 and 14 reveal that none of the approaches cured the problems associated with bidirectional bending for low-order and fully parameterized elements.
In the case of the high-order element, ANCF 3343, the best results in the z-and y-axes came from the ECM and EDG-4 approaches. The EDG-4 approach is closer to the solid element-based solution, and the ECM is closer to the beam-based commercial solution (Fig. 15).
Interestingly, regardless of the number of additional degrees of freedom, there were no significant differences in the EAS-based solutions despite element type.

Eigenfrequencies of the simply supported beam
The next numerical test was an eigenfrequencies analysis, which was used for ANCF beams in multiple works [3,7,12,14,19,22,23] etc. There are multiple reasons for deciding to consider the eigenfrequency analysis. For example, Orzechowski and Shabana [23] mentioned that examining the eigenvalues might provide some explanations for locking behavior. Ebel et al. [7] stated that the eigenfrequency analysis can be used as a linearized dynamics test, because it offers the advantage of coordinate-free frequency and vibration mode comparisons. In [14], the eigenfrequency analysis was used to see if beam elements could capture the Poisson modes. The solution of the nonlinear problem depends on the set of deformation shapes assumed for the models. In this study, the example of a simply supported beam used in [3,7,12,19,22] was analyzed. A beam having a rectangular cross section with width W = 0.4, height H = 0.4, length L = 2, Young's modulus E = 1 × 10 9 N m 2 , density ρ = 7850 kg m 3 , Poisson's ratio ν = 0.3 was analyzed. Each model was made up of 128 elements, as was done for several modifications for the gradient-deficient elements in [7]. In this test, the boundary conditions of an simply supported beam was undertaken.
The analytic solutions mentioned in Tables 9, 10, and 11 were based on the solutions and derivations provided in the works [1,17].
Of the solutions obtained for gradient-deficient, fully parameterized, and high-order elements, only the ECM modification was in full agreement with analytical solutions. Tables 9, 10, and 11 show that only ECM modification is in fully agreement with analytical solutions. Additionally, for the kinematic-based locking alleviation techniques, EAS and EDG, the number of additional DOFs did not significantly influence the final results.

Conclusion
Locking alleviation techniques for ANCF elements was the subject of this research. The authors considered all known techniques for continuum-based ANCF elements including the enhanced assumed strain method, the enhanced continuum mechanics approach, the strain split method, and a newly proposed enhanced deformation gradient method. These methods were applied to low-order, fully parameterized and highorder ANCF beam elements. All possible combinations were investigated and compared. None of the modifications performed better than the others for all the presented tasks. For small bending cases, the best results were obtained via the ECM approach for all element types. However, the results presented by other approaches for the high-order element were within reasonable limits and outperform the 3363 element excepting the SSM approach. When the large bending case was considered for the low-order and fully parameterized elements, the ECM approach again provided the best solutions. How- ever, in the case of the high-order element, the results of the EDG-4 approach were closest to the reference with the highest convergence rate. During the Princeton beam experiments, where the bidirectional loading cases are considered, the appropriate results were demonstrated by the ECM and EDG-4 approaches for only the high-order ANCF element. Furthermore, the EDG-4 provided solutions in good agreement with the solid element, and the ECM showed good correlation with beam-based commercial software solutions. The EAS-2 and EAS-4 methods did not differ significantly regardless of considered bending loading cases and element types, which lead to the conclusion that between them, the EAS-2 is the preferable choice, because it is less computationally intensive. For eigenfrequency analysis, the best performance was ultimately demonstrated by the ECM approach. Among the considered modifications, the locking problem can be completely cured only for the highorder element and only using the ECM and EDG-4 approaches. However, the ECM approach is limited and does not allow working with material models other than Kirchhoff-Saint-Venant. Therefore, the 3343 element with EDG-4 modification is the most universal solution. Further comparison with the element denoted 3363, and assuming as the optimal solution for various problems in [7], shows that 3343 with the EDG-4 element provides better results with the requirement of fewer DOFs. Therefore, the newly proposed approach, namely EDG-4, is a good tool to improve the performance of the high-order ANCF continuum-based beam element in the alleviation of Poisson locking in conventional load cases. Further research with practical applications to ensure the performance of the 3343 might be interesting. Additional higher-order deformation modes could be beneficial in recreating deforma-tions caused by centrifugal forces such as rotor dynamics.
Funding Open access funding provided by LUT University (previously Lappeenranta University of Technology (LUT)). Open Access funding enabled and organized by Projekt DEAL. The work provided by the Academy of Finland (Application No. 299033 for funding of Academy Research Fellow). We would like to thank for the generous grant that made this work possible.

Data Availability
The manuscript has no associated data.

Conflict of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

Appendix A1 On ANSYS converge problems
In Sect. 5.3, the authors mentioned some convergence problems within ANSYS commercial software and put the element number as 128. As a matter of fact, with increasing element number, the solutions stopped to converge to analytical solution and started significantly stepwise increasing, as shown in Fig. 16. The authors wanted to demonstrate this convergence problem for the case of Princeton beam experiment from Sect. 5.3, for two rotational angles β = 0 • and β = 30 • in the case F pb 2 .  The problem was not limited by the beam presented in Sect. 5.3, and this feature was observed for any rod model. Such high division can be considered as unnecessary exaggeration, but for the deformation analysis of the large beams' displacements or contact beam problems, this feature should be taken into account.