Recent Advances in the Absolute Nodal Coordinate Formulation: Literature Review From 2012 to 2020

where x, y, and z are the element coordinates. Accordingly, the ANCF element leads to a constant generalized mass matrix, which is favorable for dynamic simulations of multibody systems from a computational standpoint. The nonlinear inertia forces, including the centrifugal and Coriolis forces, are identically zero. The ANCF finite elements can be implemented in general nonincremental multibody dynamics computer algorithms in a straightforward manner without a need for ad hoc updating schemes for the finite rotation and deformation, required in large rotation vector finite element formulations.

In general, ANCF elements for beams, shells, and solids are classified according to the way the nodal gradient vectors are utilized. An element using the three position vector gradients, $rx=∂r/∂x$⁠, $ry=∂r/∂y$⁠, and $rz=∂r/∂z$ at each nodal point as in Eq. (1), is called a fully parameterized ANCF element, while an element using one or two gradient vectors at each node is called a gradient deficient ANCF element [4]. In some literature, ANCF elements using a complete set of the element coordinates x, y, and z in their shape functions are labeled as fully parameterized elements, regardless of the choice of the position vector gradients. The former definition is, however, more common and employed in this paper. In addition to the gradient vectors, the nodal coordinates are augmented by higher order derivatives of the global position vector, such as $ryy=∂2r/∂y2$ or $rzz=∂2r/∂z2$⁠, to capture complex cross section deformation modes. Such an element is generally called a higher-order ANCF element, and some mixed-coordinate ANCF elements can fall into this category. A different parameterization leads to different element properties associated with deformation modes that can be captured, degrees of interelement continuity, element convergence characteristics, and computational efficiency. As pointed out by Shabana [3], ANCF elements are characterized by (1) the use of the position vector gradient as nodal coordinates; (2) the constant mass matrix; (3) exact modeling of rigid body motion, i.e., strains are invariant under arbitrary rigid body motion; and (4) the use of the nonincremental solution procedure. Thus, elements possessing the above four properties are considered in this review paper, with a few exceptions for specialized formulations based on ANCF. The use of the nodal gradient vector r_x along the space curve leads to $C1$ interelement continuity for ANCF beam elements. It is worth noting that $C1$ interelement continuity is not always ensured for some beam and plate/shell elements, known as ANCF gradient deficient elements, including the ANCF plate element using the gradient vectors r_x and r_y tangent to the midplane. Higher degrees of continuity can be enforced by introducing continuity conditions to higher order ANCF elements, and this process is equivalent to changing a knot multiplicity of the B-spline geometry with the Cox-de-Boor recursion formula [5]. That is, the B-spline and nonuniform rational B-spline (NURBS) geometry representations of a curve, area, and volume can be directly converted to the ANCF elements for the analysis. Current engineering design requires a large amount of time and effort to generate finite element mesh from computer-aided design (CAD) data because the B-spline and NURBS functions cannot be directly converted to conventional finite element mesh. Thus, a simple linear mapping from the B-spline/NURBS geometry to the ANCF elements can facilitate the integration of CAD and analysis (ICADA) in the engineering design process.

A comprehensive review article on ANCF was published in 2013 [4], summarizing the theoretical developments and engineering application of ANCF elements with 137 journal and conference papers, textbooks, and technical reports. The topics covered in the review paper [4] included the ANCF beam and plate elements, their locking alleviation techniques, material nonlinearities, joint formulations, modeling of slope discontinuities, ICADA, and engineering applications. Since then, ANCF element capabilities have been enhanced substantially, and the range of applications has grown extensively. For example, the number of publications on the ANCF approach per year is summarized in Fig. 1 from 1997 to 2020. In this figure, blue bars show the number of publications from 1997 to 2011 covered in the review paper [4], while the number of publications has increased at a much faster pace since 2013, as shown in the red bars. In particular, the ANCF element library has grown substantially for beam, plate/shell, solid elements, addressing drawbacks of ANCF elements recognized at the time of publication of the first review paper [4]. Furthermore, various attempts have then been made to improve the computational performance of ANCF models by means of high-performance computing techniques, model order reduction, arbitrary Lagrangian–Eulerian (ALE) approach, etc. Multiphysics simulations have become active research areas to solve challenging fluid–structure interaction problems. In addition, the engineering application area has extended. These recent topics were not addressed in the first review paper. Therefore, this paper aims to conduct a comprehensive literature review on the ANCF elements and their applications published from 2012 to 2020 and supplement the first review paper [4] by overviewing the state-of-the-art ANCF simulation capabilities.

In this section, ANCF beam elements and the related issues published between 2012 and 2020 are summarized. Note that beam elements proposed before 2013 were well summarized in the first review paper [4], as well as the one by Nachbagauer [6].

The fully parameterized beam element in Fig. 2(a) leads to the general motion description that allows for describing axial, shear, bending, torsional, and cross-sectional deformations with the three nodal gradient vectors r_x, r_y, and r_z, as originally proposed in 2001 [7]. A continuum mechanics approach, therefore, can be used to define the elastic forces, for which general nonlinear material models such as hyperelasticity can be incorporated into the beam model in a straightforward manner. For the fully parameterized beam element, the element quadrature for a uniform circular cross section was examined by Orzechowski [8], while appropriate quadrature points for a multilayer beam were proposed by Lan et al. [9] to model a high-voltage electrical wire. However, the continuum mechanics approach requires many quadrature points for the volume element to describe complex cross section. The classical beam theory, for which the generalized elastic forces are calculated by the line integral with an assumption of constant section properties, has been utilized. Zheng et al. [10] presented that such a structural mechanics approach led to better convergence characteristics than the continuum mechanics approach that exhibits severe element locking. They also examined the curvature definition of the beam formulation as curvatures are defined by the spatial derivative of the cross section angle in existing low-order beam elements. In contrast, curvature definition in differential geometry has no coupling with shear deformation, and it is uniquely defined by the space curve. As a special case of the fully parameterized beam element, the Euler–Bernoulli beam element can be obtained by constraining the cross section deformation by imposing kinematic constraints, as presented by Otsuka and Makihara [11].

Locking alleviation for the fully parameterized beam element with the continuum mechanics approach has been addressed for many years, as summarized in the previous review paper [4]. To alleviate the Poisson's locking, Patel and Shabana [12] and Shabana et al. [13] proposed to decouple the bending and transverse normal strains by splitting the Green–Lagrange strain tensor into the one associated with the beam centerline and the higher-order terms that describe the beam bending, torsion, and cross-sectional distortion. The Poisson's locking of the fully parameterized beam element was then effectively alleviated.

An ANCF beam element that has only the axial gradient vector r_x, as shown in Fig. 2(b), is called a cable element as the element can describe the beam bending and axial deformations only, as presented by Gerstmayr and Shabana in 2006 [14]. The pure torsion about the axial gradient vector as well as the cross-sectional deformation, including shear deformation, cannot be captured due to the use of the single gradient vector along the beam axis. To account for the Poisson's effect in the cable element, section properties can be defined as a function of the axial strain as presented by Sheng et al. [15]. Gruber et al. [16] used an axial rotation angle at each node to define directors of the beam cross section, and the nodal angles were interpolated by a linear polynomial. It is worth noting that such a parameterization was proposed by Dmitrochenko and Pogorelov [17] in 2003. However, it was not implemented for the analysis. The cross section directors defined by the nodal rotation angles were used to describe pure torsion that cannot be captured in the original ANCF cable element. The mass matrix is, however, no longer constant in this specialized element.

The beam curvature definitions were examined by Zhang and Luo [18]. In their work, the geometrical measure of curvature, $κg=|rs×rss|$⁠, and the one in differential geometry, $κe=|rss|$⁠, were compared, where $rs=∂r/∂s$ and $rss=∂2r/∂s2$ with s being the beam arc-length coordinate. The exact curvature $κe$ in differential geometry exhibited better convergence and more accurate solutions. However, it leads to more complex expressions for the generalized elastic forces. Zemljaric and Azbe [19] derived the analytical form of the elastic forces with precomputed matrices, which led to a substantial computational speedup in the elastic force calculation. To address the complexity of the elastic force calculation in ANCF, Hara and Watanabe [20] introduced the local nodal coordinates to define the elastic forces with respect to the local coordinate system. The local nodal coordinates are then transformed back to the global ANCF nodal coordinates by algebraic constraint equations. Computational time reduction by the local coordinate approach was also discussed. Liu et al. [21] proposed formulations for the curved ANCF beam element, while radial basis functions were introduced to the curved ANCF cable element by Chen et al. [22], showing better convergence characteristics for the curved beam with the radial basis functions.

Various shear deformable gradient deficient beam elements were developed after 2012 as summarized in Figs. 2(c)–2(f), which successfully alleviated element lockings exhibited in the fully parameterized beam element. These gradient-deficient beam elements have at least one of the transverse gradient vectors r_y and r_z per node to describe the cross section rotation and deformation. Nachbagauer et al. [23] proposed a three-node shear deformable gradient deficient beam element, as shown in Fig. 2(c). Two transverse gradient vectors r_y and r_z are used at every node in the element. The space curve of the beam centerline is then approximated by the quadratic polynomial with the nodal position coordinates, while transverse gradient vectors are used to describe the rotation and deformation of the beam cross section. Since there are no longitudinal gradient vector coordinates at interface nodes, the element has $C0$ interelement continuity. Both a general continuum mechanics approach and a structural mechanics approach with an assumption of rigid cross section were applied to define the elastic forces. The numerical results in the static and eigenvalue analyses were in good agreement with those obtained by ansys, the analytical Timoshenko beam formulae, and a geometrically exact beam formulation. The element performances for the buckling analyses and dynamic multibody simulations were examined by them in Ref. [24], showing good agreement in results with reference solutions.

Ren [25] proposed three-node shear deformable beam elements, as shown in Figs. 2(d) and 2(e): one with r_y and r_z at interface nodes with $C0$ interelement continuity, while another with r_x, r_y, and r_z with $C1$ interelement continuity. The unique feature of these elements is that the middle node has different coordinates than those at the interface nodes. With the structural mechanics approach, the elements did not exhibit locking problems. The Poisson's locking was eliminated by assuming zero transverse stresses. The results for benchmark problems were in good agreement with the reference solutions in static and dynamic analyses. Otsuka et al. [26] proposed a two-node asymmetrically gradient-deficient shear deformable beam element, as shown in Fig. 2(f), which uses gradient vectors along the beam axis r_x and only one of the transverse gradients, i.e., r_y. The element is suited for modeling a slender wing for aircraft, for which the wing thickness is much smaller than the width. Thus, only transverse gradient along the wing width direction is used to reduce the model dimensionality. The cross section deformation was constrained by internal constraint equations for their applications. A wing model with the proposed elements was validated against the wind tunnel test data.

Whereas various section properties can be incorporated into the elastic force calculation with the structural mechanics approach for both fully parameterized and gradient deficient beam elements as explained above, the material model is generally limited to pure elastic. The implementation of general nonlinear material models, including hyperelasticity and plasticity, into conventional structural finite elements, requires a specialized and complex procedure. Furthermore, change of the section properties due to cross section deformation would not be negligible for large deformable thick beam structures. Therefore, the continuum mechanics approach plays a vital role in applications requiring nonlinear material models. In order to describe more complex cross section deformation modes, including warping, various higher order ANCF beam elements have been proposed by introducing the higher order derivatives of the global position vector as nodal coordinates, along with the gradient vectors, as shown in Figs. 2(g) and 2(h) as examples. Shen et al. [27] compared different orders of polynomials associated with the transverse gradients, and the higher order derivatives were introduced to the ANCF nodal coordinates to capture the cross section distortion and warping. The higher order polynomials in the transverse direction allowed for alleviating locking problems of the fully parameterized beam element.

Later on, Ebel et al. [28] proposed several three-node higher order beam elements without the longitudinal gradient. A representative example is shown in Fig. 2(h). They pointed out that conventional locking alleviation techniques such as selective integration are not effective under certain load conditions. Higher order polynomials are more effective in alleviating the ANCF locking problems. Orzechowski and Shabana [29] examined the performance of the higher order ANCF elements proposed by Shen et al. [27] in Fig. 2(g), with particular emphasis on warping and Poisson's locking. It was pointed out that the initial values for the higher order derivatives are difficult to obtain for initially curved geometry. Therefore, a transformation was proposed to convert the nodal coordinates of the higher order derivatives at interface nodes to the position coordinates at new internal nodes, leading to a mixed-coordinate ANCF beam element.

Classical rotation-based beam formulations interpolate the position and rotation fields independently. Shabana [30] refers to these elements as large rotation vector formulations (LRVFs) rather than geometrically exact beam formulations (GEBFs). Ding et al. [31] pointed out that material points defined by the two independent fields occupy two different positions fields, thereby causing a coordinate redundancy problem. That is, a space curve obtained from the rotation field is inextensible, while that obtained from the position field is extensible. Furthermore, Shabana and Patel [32] pointed out that bending and shear deformation modes must be kinematically decoupled, while they are coupled in the curvature definition used in LRVF.

To address these issues, Shabana proposed a consistent rotation-based formulation (CRBF) with the ANCF kinematic description that employs a single interpolation to define a unique rotation field [30]. The idea behind this approach is to relate the displacement gradient tensor obtained by the ANCF gradient nodal coordinates to those defined by the rotation parameters (e.g., Euler angles), along with the extensibility parameter. The velocity transformation matrix is then derived, and the equations of motion of the fully parameterized ANCF elements were converted to those expressed in terms of the three rotation and extensibility parameters as the nodal coordinates. Shear deformation can be considered in the mapping, while the cross section deformation is eliminated through the velocity transformation. Because of the time-variant velocity transformation matrix, the mass matrix is no longer constant, and nonlinear inertia forces, including the Coriolis and centrifugal forces, appear in the final form of the equations of motion. With this approach, the inconsistent independent interpolation for the position and rotation fields in existing LRVF and GEBF can be bypassed, leading to a consistent rotation-based formulation with the general continuum mechanics theory. The ANCF/CRBF approach was applied to 2D beam elements [33,34], 3D beam element [35], and plate element [36]. The computational performance as well as accuracy of the ANCF/CRBF elements was compared with those of the ANCF elements [35]. It is worth noting that the transformation can be constructed for infinitesimal rotation parameters as well for small deformation problems [37,38]. Such an element can be utilized for the floating frame of reference formulation to account for initially curved geometry described by the constant gradient vectors. Furthermore, it was pointed out that the B-spline/NURBS geometry can be converted to the consistent rotation-based finite elements through the ANCF kinematic description.

To evaluate the performance of the ANCF beam elements, several validation problems were suggested for comparison with GEBF beam elements [39–41]. Bauchau et al. [39] reported that the 2D Euler–Bernoulli ANCF beam element with Green–Lagrange axial strain and the spatial measure of curvature led to less accurate kinematic solutions (e.g., curvature and strain) than those of the GEBF beam element, as well as in the static analysis. For spatial problems, Bauchau et al. [40] conducted validation of the fully parameterized ANCF beam element based on the cross-sectional coordinate system approach, the three-node shear deformable gradient deficient ANCF beam element [23], and the three- and four-node GEBF beam elements against test data of the Princeton beam experiment. The numerical results of the GEBF and the three-node shear deformable gradient deficient ANCF element agree well with the experimental data, while the fully parameterized ANCF beam element led to less accurate solutions due to the inability to predict the torsional behavior in this validation test problem.

Four validation problems, the Princeton beam experiment, the four-bar mechanism, the lateral buckling of a beam, and the stability of a rotating shaft were proposed to evaluate beam elements implemented in various flexible multibody dynamics simulation codes [41]. The three-node shear deformable gradient deficient ANCF beam element [6] in HOTINT code produced results that agreed well with the other GEBF beam elements. Wang et al. [42] applied the three-node ANCF beam element to the analysis of high-speed rotating shaft problems, and the results were compared with those of superelements, showing good agreement with the reference solutions.

In this section, the ANCF plate/shell elements and the relevant topics published between 2012 and 2020 are summarized.

ANCF plate and shell elements are classified according to the way the gradient vector coordinates are introduced to describe different deformation modes. The use of the three nodal gradient vectors r_x, r_y, and r_z leads to a general kinematic description that allows for capturing not only the midplane deformation of the plate but also the transverse shear and thickness deformations. Such a fully parameterized plate/shell element, originally proposed by Mikkola and Shabana [43] in 2003 as shown in Fig. 3(a), is parameterized as a volume element. On the other hand, the use of only the gradient vectors tangent to the midplane, i.e., r_x and r_y, leads to a gradient-deficient thin plate element (36DoF) that is parameterized as an area element shown in Fig. 3(b) [17,44]. Hyldahl et al. [45] stated that the original 36-DoF ANCF plate element, presented by Dufva and Shabana [44], with the shape functions based on an incomplete quartic polynomial leads to only $C0$ interelement continuity. The bi-cubic Hermite polynomials were, therefore, suggested to ensure the $C1$ interelement continuity for the thin plate element, as in the 48-DoF thin element in Fig. 3(c) utilizing the higher order derivative r_xy as the additional nodal coordinates. This element was originally proposed by Dmitrochenko and Pogorelov in 2003 [17]. Whereas these ANCF thin elements produced results comparable to those from commercial software nastran for uniform mesh models, Hyldahl et al. [45] pointed out that these elements are more sensitive to nonuniform mesh, resulting in less accurate solutions and increased computational cost. Therefore, Yan et al. [46] adopted a curvilinear coordinate system to model initially curved structures with the ANCF thin plate elements. It was presented that the issue of the nonuniform mesh was eliminated by the proposed approach, and the element convergence was comparable to that of a shell element in abaqus.

Yu and Shabana [47] pointed out that the use of the higher order derivative r_xy in the 48-DoF plate element poses a challenge in defining an arbitrary reference configuration for initially curved structures. To address this issue, they developed a procedure to convert the r_xy vectors at interface nodes to the position vectors at new internal nodes as shown in Fig. 3(d). It led to a 48-DoF mixed-coordinate ANCF thin plate element parameterized by the position and gradient vectors, r_x and r_y, only. It was also demonstrated that the proposed element can be converted to a B-spline surface with a linear transformation. The element was further extended with bi-cubic rational Bezier functions by Pappalardo et al. [48] to enable the accurate description of circular and conic geometries, which cannot be precisely modeled by nonrational polynomials used in the conventional ANCF elements. Chen et al. [49] proposed a thin plate element with the radial point interpolation to accurately describe an initially curved shell geometry.

A simplified elastic force model was explored for the 36DoF plate element in Fig. 3(b) by Otsuka and Makihara [50]. With an assumption of small midplane deformation, a constant stiffness matrix for the out-of-plane plate deformation (i.e., curvatures) can be obtained. The component mode synthesis technique is utilized to reduce the model dimensionality associated with the out-of-plane plate deformation. Furthermore, the mean strain approximation was used to obtain weakly nonlinear in-plane elastic forces for computational speedup.

To account for general nonlinear material models, including laminated composites in various engineering applications, the 48-DoF fully parameterized shear-deformable plate/shell element has been used with the continuum mechanics approach, as presented in Refs. [51] and [52] for anisotropic materials and in Ref. [53] for composite materials. On the other hand, element lockings have made the original fully parameterized plate element less attractive. Matikainen et al. [54], therefore, examined locking alleviations for the ANCF shear-deformable plate element. The transverse shear locking was alleviated by linear interpolation of the transverse shear strains with the assumed natural strain technique. Furthermore, they proposed a higher order shear-deformable plate element by introducing a higher order derivative vector r_zz to the fully parameterized plate element, leading to a 60-DoF element as shown in Fig. 3(e) to alleviate the element locking. They explained the importance of eliminating the Poisson thickness locking in these elements. Ebel et al. [55] proposed eight- and nine-node elements with either r_z and its derivative r_zz as shown in Fig. 3(f) or only r_z at each node as shown in Fig. 3(g). It is worth noting that there are no in-plane gradient vectors r_x and r_y used in these elements, and the midplane deformation is parameterized by the nodal position coordinates only. The idea behind this parameterization is similar to that of the shear-deformable gradient-deficient beam element using the transverse gradient vectors r_y and r_z. The convergence of these plate elements was examined under different loading and boundary conditions. In this category, a four-node shear-deformable shell element was developed by Yamashita et al. [56] and Valkeapää et al. [57] with the continuum and structural mechanics approaches, respectively. The transverse shear and thickness lockings were alleviated by the assumed natural strain technique, while the in-plane shear locking was eliminated by the enhanced assumed strain approach. The eight-node element was presented by Olshevskiy et al. [58] as well.

Absolute nodal coordinate formulation triangular plate elements have been proposed to produce finite element mesh for complex structures. Mohamed [59] proposed a fully parameterized triangular ANCF plate element in Fig. 4(a), showing good agreement in results with those obtained by the quadrilateral thin plate element. Pappalardo et al. [60] then proposed three fully parameterized triangular plate elements with volume parameterization of the gradient vectors: three-node triangular plate elements with two different cubic polynomials in Fig. 4(b) and a four-node mixed-coordinate element with a complete cubic polynomial in Fig. 4(c). The volume parameterization facilitates the development of shape functions as well as the element assembly.

Olshevskiy et al. [58] studied thin triangular ANCF plate elements with and without the second-order derivative vector r_xy in Figs. 4(d) and 4(e), respectively. It is worth noting that such a parameterization in Fig. 4(d) was proposed by Dmitrochenko and Pogorelov [17] in 2003. Olshevskiy et al. [61] also proposed three- and six-node shear deformable triangular plate elements with the gradient vector r_z along the thickness in Fig. 4(f). Ren [62] proposed thin triangular plate elements with nonvertex nodes. The elements do not require numerical quadrature to derive the elastic forces. Three types of triangular ANCF plate elements were developed by Chang et al. [63] with the cubic Bézier triangle: a thin plate element as well as the lower and higher order shear deformable elements, as shown in Figs. 4(g)–4(i). They presented a linear transformation between the proposed thin triangular plate element and the Bézier parametric representation to facilitate direct conversion from CAD geometry. It is worth noting that the position coordinates of a control point of the elements are included in the nodal vector, along with ANCF nodal coordinates for each type of element. Wang [64] proposed two thin triangular elements with the area parameterization of the gradient coordinates for initially curved structures.

In addition to extensive studies on the ANCF beam and plate/shell elements, various fully parameterized ANCF solid elements have been developed since 2014. Olshevskiy et al. [65] were the first to develop planar triangular and square ANCF elements by introducing the gradient nodal coordinates r_x and r_y to the conventional solid elements, as shown in Figs. 5(a) and 5(b). Pappalardo et al. [66] introduced the area parameterization of the gradient vectors to the planar ANCF triangular element, as shown in Figs. 5(c) and 5(d). The gradient vectors are defined by coordinate lines along the sides of the triangular element. Thus, a standard FE assembly procedure can be utilized. The procedure is similar to the one in their work for the plate elements [60]. Lan et al. [67] proposed an incomplete cubic polynomial different from that in Ref. [66] for the triangular element. The element proposed by them led to a better convergence than the one with the complete cubic polynomial in Ref. [66]. Langerholc et al. [68] exploited the general kinematic description of the planar square ANCF solid element to perform the digital image correlation. That is, the ANCF solid element was used to estimate rigid-body movement and deformation of an object in a picture by correlating the gray scale intensity values before and after the object's deformation.

Olshevskiy et al. [69] proposed a spatial 96DOF eight-node ANCF solid brick element, as shown in Fig. 5(e), while a four-node tetrahedral ANCF solid element shown in Fig. 5(f) was proposed by Olshevskiy et al. [70]. Pappalardo et al. [71] proposed tetrahedral solid elements based on the volume parameterization of the gradient vectors defined by coordinate lines along the sides of the tetrahedral element, as shown in Figs. 5(g) and 5(h). The four-node ANCF tetrahedral element was developed with an incomplete cubic polynomial, while the use of the complete polynomial requires four additional internal nodes with the position coordinates as shown in Fig. 5(h). The elements showed a good convergence performance compared to conventional tetrahedral elements without gradient coordinates. Wang et al. [72] proposed an incomplete cubic polynomial different from the one in Ref. [71] for the ANCF tetrahedral element, showing a better element convergence than the one proposed in Ref. [71]. The element performance was similar to the one with the complete cubic polynomial in Ref. [71] that requires additional degrees-of-freedom. The need for curvature continuity conditions was discussed by Shabana and Zhang [73] to ensure a higher degree of continuity at the element interfaces for the ANCF solid elements in the analysis of soft materials and fluid.

This section summarizes modeling procedures for multibody dynamics simulations with ANCF approaches, including joint constraints, contact modeling, material nonlinearities, computer-aided design and analysis, and structural optimization.

Formulations of mechanical joints require relative rotation constraints between bodies. Since no rotational parameters are used in the ANCF elements, an orthogonal joint coordinate system is defined at the constraint definition point using the nodal gradient vector coordinates. However, it leads to highly nonlinear constraint equations in terms of the gradient coordinates at the velocity and acceleration levels. Therefore, to impose clamped joints that connect large deformable bodies with other rigid or flexible bodies, Shabana and Eldeeb [74] proposed a simpler approach to define the orthogonal triad using gradient vectors for ANCF elements, eliminating the use of the tangent or cross section frames. The compliant continuum-based joint formulation was presented by Wallin et al. [75] to model track chains with the ANCF approach, and the numerical results were compared with the conventional constraint and penalty approaches with rigid links. With the ANCF approach, some mechanical joints can be modeled by simple linear constraint equations with gradient coordinates.

Detailed joint modeling considering the clearance, lubrication, and deformation has been investigated for flexible bodies modeled by the ANCF approach. Elastohydrodynamic lubricated cylindrical joint was introduced to a flexible multibody system by Tian et al. [76,77]. In their studies, the outer hollow cylinder of the joint was modeled by the standard solid elements to account for the local joint deformation, while the rigid inner shaft was connected to the flexible ANCF beam model. This approach was applied to the elastohydrodynamic lubricated spherical joint [78]. Similarly, Zhao et al. [79,80] predicted wear caused by a joint clearance of a flexible mechanism with the 2D ANCF beam elements, where the local joint deformation is modeled by standard solid elements. The effect of joint clearance on control of a multilink flexible robotic arm was studied in Ref. [81], while the joint clearance was treated as uncertainty in Ref. [82]. Tang and Liu [83] presented a modeling procedure to account for a sliding joint clearance with the ANCF cable element. The ANCF-based multibody simulation framework considering the joint clearance and lubrication was applied to a robotic manipulator [84], multilink mechanisms [85–87], and a flexible moving platform mechanism [88,89].

Because of the nonrotational parameterization of ANCF elements, it is not straightforward to apply external torques and moments to ANCF large deformable bodies. In a standard approach, the generalized external torque vector is obtained by the principle of virtual work, and the element cross section angle is defined by the gradient vectors. Based on this approach, Otsuka and Makihara [50] proposed an external torque formulation for a torsional spring and damper between ANCF thin plate elements. By defining a cross section frame with the gradient vectors, an external torque proportional to a relative rotation angle can be formulated in a straightforward manner. On the other hand, Liu et al. [90] proposed another approach to apply an external torque to the fully parameterized ANCF beam element through the outer product relationship between the external torque and gradient vectors.

Contact and friction modeling is an integral part of the analysis of multibody system dynamics in many engineering applications. Gantoni et al. [91] and Hamper et al. [92] stated that the accurate ANCF geometry representation is beneficial in solving contact problems as the higher degree of surface continuity is ensured in the contact geometry analysis. It allows for accurate prediction of contact zones and forces on large deformable bodies. Khude et al. [93] introduced a continuous contact force model with a spherical decomposition technique for the ANCF cable element. With this approach, collision detection can be performed efficiently with many spheres defined along a large deformable cable by leveraging graphics processing units (GPUs) computing. Yang et al. [94–97], Wang et al. [98], and Ma et al. [99] used the compliant contact force model for modeling underground cables and drill string systems, while Zhang et al. [100] analyzed the steel cable capture mechanism using the cable element. Unlike discrete contact models, the master–slave contact approach was introduced to the cable element by Wang et al. [101,102] to detect continuous contact zones precisely, and elastoplastic behavior of beams resulting from frictional contact was examined in Ref. [103]. Sun et al. [104] presented a segment-to-segment frictionless contact approach considering the hyperelasticity of a planar ANCF solid element. Yang and Shen [105] introduced a lumped-parameter contact model that considers reverse sliding for soft robotic fingers.

Considering structural damping in ANCF elements is crucial to correlate simulation and experiment results. It also helps improve the computational efficiency as the ANCF high-frequency modes can be damped out. Kim and Yoo [106] summarized simpler structural damping models, including the Rayleigh proportional damping and the quadratic velocity damping with a drag coefficient. Because of the simplicity in modeling, they were employed for experimental correlations in several studies [26,107,108]. Zhao et al. [109] employed the Kelvin–Voigt visco-elastic model for the 2D fully parameterized beam element, while a viscoelastic laminated composite shell was developed by Shen et al. [110] with the 36-DoF ANCF thin shell element to examine the relaxation behaviors. Grossi and Shabana [111] proposed a viscoelastic constitutive model based on the Navier–Stokes equation to damp out insignificant high-frequency modes without dissipating the kinetic energy associated with the rigid body motion.

Various nonlinear material models, including hyperelastic and incompressible rubberlike materials, have been utilized for ANCF elements to account for complex material behavior exhibited in engineering applications such as soft robotics and biomechanics problems. The material models, such as Yeoh, Moony-Rivlin, Gent, and neo-Hookean models, have been implemented and used in Refs. [108] and [112–115]. Xu and Liu [115] employed the Yeoh model to model a silicone beam with the ANCF approach and validated it against the test data. They also compared the Neo-Hookean, Gent, and Arruda-Boyce material models for a large deformable silicone beam with the higher order ANCF beam element, showing that the Arruda-Boyce model led to the best agreement with experimental results [108].

Due to substantial time and effort required to convert CAD geometry to finite element mesh for the analysis, the use of NURBS geometry representation in both CAD and the analysis has been studied extensively in the computational mechanics community as the isogeometric analysis toward integration of CAD and analysis. On the contrary, Shabana et al. [116] examined the limitations of the B-spline representation as an analysis tool. Because of the hierarchical B-spline recurrence formula, the number of degrees-of-freedom tends to increase significantly compared to conventional finite element models. Furthermore, it is not, in general, straightforward to model structural discontinuities using the B-spline/NURBS geometry representation. In contrast, B-spline geometry can always be converted to ANCF geometry through a linear transformation that relates the B-spline control points and the ANCF position and gradient coordinates, thereby allowing for seamless integration of the computer-aided design and analysis (ICADA). With this important feature, CAD geometry can be directly converted to ANCF mesh for use in flexible multibody dynamics simulations, and the engineering design process can be facilitated.

Whereas conversion from the B-spline to ANCF representations is possible, the converse is not always possible [116]. The conversion from the B-spline to ANCF parameterizations was demonstrated for the ANCF thin plate element [117], and ANCF cable elements [118–120], for which equivalence of the B-spline and ANCF elements was examined by Yamashita and Sugiyama for different degrees of continuity and orders of polynomials [121]. It is worth noting that the degree of continuity can be easily controlled by changing a knot multiplicity in the B-spline geometry representation, while linear interelement continuity conditions need to be imposed for ANCF elements. A method to integrate a localized NURBS surface geometry into fully parameterized ANCF elements was presented in Ref. [122]. As an engineering application, a wind turbine blade was modeled with the ANCF thin shell elements through the conversion from the B-spline CAD geometry [123,124].

The ANCF finite element models have been utilized to optimize the performances of mechanical systems under design constraints. Structural optimization is, in general, classified into size, shape, and topology optimizations. In size optimization, material dimensions such as thickness and cross section area are defined as design variables. In such a parametric optimization, optimum shapes being generated are limited by the selection of the design variables. On the other hand, shape optimization aims to create an optimum shape of mechanical components by shifting the finite element nodes directly. The topology optimization goes a step further to seek an optimum solution by changing topologies of mechanical components through addition and removal of finite element nodes. Although topology-optimized shapes are difficult to manufacture, recent advances in 3D printing technologies facilitate topology optimization in practical manufacturing fields. Furthermore, the use of the ANCF geometry is advantageous for structural optimizations as the B-spline geometry in CAD can always be converted to the ANCF model for the analysis, as discussed in Sec. 5.4.

Sun et al. [125] developed a computational algorithm for structural optimization with ANCF cable and fully parameterized beam elements. By employing the equivalent static load (ESL) method, the structural optimization problem for a flexible multibody system was solved as a static problem to reduce the computational cost for function evaluations of the iterative optimization process. The heights of flexible components in a mechanism were optimized with the proposed procedure as size optimization. Using the procedure to integrate the localized NURBS and ANCF surface geometry, a shape optimization was performed by He et al. [126] using the fully parameterized beam element. Sun et al. [127] performed topology optimization for flexible links of a mechanism modeled by the ANCF thin plate element with ESL. The morphing movable component (MMC) method was employed to perform the topology optimization with fewer design variables. The topology optimization procedure using MMC and ESL for flexible multibody systems was also applied to a flexible mechanism modeled with the ANCF solid elements [128] and a variable-length structure with ALE ANCF plate and solid elements [129,130]. Similarly, frequency-response optimization was conducted for rotating structures with the MMC-based topology optimization [131,132]. Pi et al. [133] discussed the first-order sensitivity analysis procedure for multibody systems with the 2D fully parameterized elements.

Furthermore, ANCF elements were integrated into stochastic analysis procedures to account for uncertainty associated with material properties, geometry, and many other design parameters in the engineering design of structural multibody systems. Polynomial chaos expansion was used to obtain stochastic responses of large deformable bodies modeled by ANCF elements under uncertain material properties [134] and uncertain lengths [135], while Chebyshev collocation method was utilized in Refs. [136] and [137].

Despite many unique features of the ANCF approach, the use of gradient vector coordinates and higher order polynomials tends to increase the model dimensionality, resulting in an increase in computational costs for large-scale flexible multibody simulations. Therefore, there were many studies aiming to address the computational intensity of ANCF models, as summarized in Table 1. Among others, model order reduction, high-performance computing techniques, recursive approaches, and ALE method tailored for the ANCF formulation are summarized in this section.

The use of model order reduction techniques was explored to reduce the model dimensionality of ANCF models. While the mass matrix is constant, the nonlinearity of ANCF elastic forces does not allow the direct use of standard modal reduction techniques for a linear system. However, the bending elastic forces for the planar Euler–Bernoulli ANCF beam element can be expressed linear in the gradient coordinates if a small axial strain is assumed. Therefore, the ANCF constant mass and stiffness matrices can be used to obtain reduced order bases through the standard Craig–Bampton component mode synthesis method or normal mode approach [138,139]. Similarly, in the thin plate element, a constant stiffness matrix can be obtained for the out-of-plane elastic forces by assuming small in-plane deformations. Thus, a reduced-order ANCF thin plate model was developed by Otsuka and Makihara [50] with the Craig-Bampton component mode method, leading to an 80% reduction in the calculation time with roughly 3% error compared to the full order model. In these approaches, however, simplification of the elastic forces is needed to obtain a constant stiffness matrix. Therefore, it was proposed that the ANCF nonlinear elastic forces were linearized successively in the dynamic simulation such that modes are updated successively at every linearization point to obtain the reduced order ANCF equations [140,141]. The linearization point changes based on the user-defined criteria.

Furthermore, the model order reduction based on the proper orthogonal decomposition (POD) was applied to ANCF models by several authors. POD is a statistical method to produce reduced order bases from the time-domain response of a nonlinear system. To this end, time history data of the generalized nodal coordinates, called snapshot data, are obtained from the full order ANCF model and then used to extract reduced order bases (modes) by singular value decomposition. It was applied to the ANCF beam elements for planar [142,143] and spatial [144,145] problems. Luo et al. [144] showed that computational time was reduced by 80% while retaining accuracy. Since a full order model simulation needs to be performed to produce snapshot data prior to the model order reduction, the POD modes are generally dependent on the parameters and simulation scenarios of the full order model; thus, the modes have to be properly adjusted according to the model parameters used in prediction [144].

For computational speedup of large-scale ANCF simulations, a parallel computing scheme with GPUs was developed by Melanz et al. [146] and Serban et al. [147] using the preconditioned iterative Krylov-subspace linear solver for parallelized implicit time integration. The proposed GPU parallel computing algorithm can be applied to general flexible multibody systems expressed by the index-3 differential-algebraic equations (DAEs).

A divided-and-conquer algorithm was implemented with the 2D ANCF Euler–Bernoulli beam elements by Khan and Anderson [148]. It led to O(N) and O(log(N)) computational complexities in serial and parallel computing, respectively, when N being the number of bodies. The transfer matrix method was applied to the 2D fully parameterized beam element by Rong [149], and recursive equations are used to solve the ANCF equations for computational speedup. Hewlett et al. [150] developed a first-order time integration algorithm for the cable element toward real-time simulations and showed its fast and robust computational performance.

The ALE approach was applied to ANCF elements to enable efficient finite element modeling of axially moving media and variable length structures undergoing large deformation. In the Lagrangian motion description, finite element nodes are attached to the material and move as it deforms, while those of Eulerian finite elements are fixed in space. In the ALE approach, a reference domain, which can move independent of the material, is introduced to allow for effective mesh updates. Pechstein and Gerstmayr [151] developed a planar ALE-ANCF beam element for the analysis of a belt drive system, where the spatially fixed reference configuration was introduced and the mass is transported through the beam element as the beam moves axially. With this approach, only contact region of the flexible belt model with a pully requires fine mesh, while other parts can be modeled with coarse mesh. It led to efficient finite element modeling of a belt drive system. Escalona [152] developed a spatial ALE-ANCF cable element for modeling reeving systems, including elevator and crane cables, toward real-time simulations. It is worth noting that the mass matrix of the ALE-ANCF elements is no longer constant because of the transport of mass through the elements. Nevertheless, Fotland et al. [153] concluded that the ALE-ANCF [152] has the best performance compared to other finite element formulations for modeling cable-pully systems.

One can define the ALE reference domain such that the element nodes can move within the element. It allows for modeling variable length structures, as presented by Sun et al. [154] with the ALE-ANCF solid elements. Zhang et al. [155] proposed an adaptive mesh refinement algorithm by automatically moving nodes according to the amount of deformation along the ANCF cable model.

Fujiwara et al. [156] proposed a planer Euler–Bernoulli beam element with variable length using the variable-domain finite element method. The multiple time-scale approach was introduced with nondimensionalization of the ANCF equations of motion. They reported that solving the nondimensional equation led to a reduction of the computational time. A similar nondimensionalization approach was applied to the 2D Euler–Bernoulli beam element [157], cable element [158] and 3D fully parameterized beam element [159].

In this section, the extension of the ANCF approach to multiphysics simulations, including fluid-structure interaction and thermomechanical coupling, is discussed.

The fluid–structure interaction is one of the major multiphysics phenomena that have drawn increasing attention in the multibody dynamics community. The number of publications in this area, utilizing the ANCF approach, is increasing as summarized in Table 2.

Abbas et al. [160,161] introduced a supersonic aerodynamic model, called the piston theory based on the potential flow assumption, into the fully parameterized plate element for the supersonic flutter analysis. Bayoumy et al. [162] utilized an aerodynamic model based on the airfoil aerodynamic coefficients obtained from a database to conduct the dynamic simulation of a wind turbine blade modeled by the 36-DoF thin plate element. The Theodorsen's unsteady aerodynamic model based on the potential flow assumption was used by Otsuka et al. [163–165] to account for the aeroelastic phenomena of aircraft wings modeled by the 2D Euler–Bernoulli beam element.

In these studies, the fluid–structure coupling was considered through fluid force models for simplicity, and the mutual coupling with the flow field around the flexible structure is not considered, i.e., a one-way coupling. Therefore, various computational fluid dynamics (CFD) models have been integrated with large deformable ANCF models to account for two-way fluid–structure coupling. Among others, the smoothed particle hydrodynamics (SPH) method was used by several authors to account for coupling with ANCF elements [166–171]. SPH is a mesh-free method widely used for modeling free surface flow. Since flow fields are parameterized by a large collection of particles, mesh refinement for the moving structural boundary is not required. This unique feature in SPH is beneficial for the large deformation analysis using ANCF models. On the contrary, the immersed boundary lattice Boltzmann method (IB-LBM) was applied to ANCF beam [172,173] and plate [174] elements. In the IB-LBM, the Boltzmann equation is utilized rather than the Navier–Stokes equation commonly used in CFD simulations. The advantage of IB-LBM is that mesh refinement for the moving structural boundary is not needed. The 3D unsteady vortex lattice method (UVLM) was applied to the 3D asymmetrically gradient deficient beam element by Otsuka et al. [26] to account for the 3D aerodynamic effect of deployable wings, including induced drags and wingtip vortices. Since the structural surface and flow wake are discretized by UVLM, rather than the 3D flow field, the computational cost is much lower than the other CFD models. Due to the potential flow assumption, however, the viscous and compressible effects, such as shock, separation, viscous drag, and wave drag, are not considered in UVLM. UVLM was also used by Yamano et al. [175] to conduct the analysis of the nonlinear aeroelastic vibration of flexible sheets.

Absolute nodal coordinate formulation elements were extended for the analysis of liquid sloshing. Sloshing is a phenomenon of liquid movement occurring in a container of tanker trucks and rockets in motion. Such a sloshing simulation requires accurate prediction of free-surface fluid flow coupled with multibody dynamics models. To address the complexity of integrating the Eulerian-based CFD models into the Lagrangian-based multibody dynamics algorithm, Wei et al. [176] proposed a total Lagrangian ANCF FE approach for liquid sloshing problems. The fluid element was developed using the ANCF solid element by introducing the fluid constitutive equations, along with the incompressibility condition. By doing so, both the structure and fluid can be modeled by the ANCF approach in general multibody dynamics algorithms. This approach was applied to the sloshing analysis of the road [177] and rail [178,179] vehicles with a liquid container. The ANCF and SPH models for liquid sloshing simulations were compared by Atif et al. [180]. They showed that the ANCF fluid model can capture the complex liquid flow with 40 times fewer degrees-of-freedom, leading to 70% computational time reduction compared with the SPH model. Pan and Cao [181] summarized existing modeling methods for the analysis of liquid sloshing and used the ANCF planar solid element to predict the fluid flow.

The thermomechanical coupling is important in various engineering applications. For example, aerospace structures are subjected to a large variation of environmental temperature, making the effect of thermal deformation non-negligible in design evaluation. The temperature of the structure becomes very high due to solar radiation when it is close to the sun, while the temperature drops significantly when it is away from the sun or shaded by planets. Furthermore, aerospace structures tend to become more lightweight and flexible. Thus, accurate prediction of thermal deformation plays a critical role in controlling unwanted thermally induced deformation and vibration. Studies considering the thermal coupling effect of the ANCF elements have increased since the 2010s.

The simplest thermal deformation model used in ANCF elements is to consider a thermal strain proportional to a temperature variation in the elastic force formulation. This method was applied to the fully parameterized beam element [182] for a bimetallic strip with a prescribed temperature variation. The discretized thermodynamic equations are coupled with the ANCF equations for the analysis of beam-like structures in [183–185] with the following assumptions: (1) the temperature gradient along the beam thickness, the heat convection, and the radiation heat change were assumed negligible; (2) thermal properties are assumed to be independent of temperature; and (3) the thermodynamic equation is coupled with the ANCF equations via thermally induced strains. This approach was employed for rigid-flexible multibody systems modeled by the natural coordinate formulation and ANCF [186]. Whereas different mesh was used for the mechanical and thermal models, recent studies [187,188] adopted unified meshing for the ANCF thermomechanical models.

Absolute nodal coordinate formulation elements have been utilized to solve challenging engineering problems that require accurate prediction of large deformation of flexible multibody systems [4]. Owing to significant advances in ANCF modeling capabilities in past years as summarized in this paper, the application area has extended between 2012 and 2020. In addition to the major application areas such as cables, rotating shafts, vehicle components, and belt-drive systems, there was a significant increase in aerospace applications in the 2010s. In this section, applications of ANCF elements and the related modeling procedures published between 2012 and 2020 are summarized.

A road vehicle is one of the representative multibody systems composed of interconnected rigid and flexible bodies, as shown in Fig. 6(a). ANCF has been used to model flexible vehicle components, including pneumatic tires and leaf springs, as explained in the first review paper [4].

Contreras et al. [189] reviewed existing deformable soil models for off-road mobility simulations and explained the potential benefits of using ANCF approaches in modeling deformable terrain as well as other vehicle components including tires. Shabana proposed the ANCF reference node (ANCF-RN) [190] to assemble a flexible tire modeled by ANCF into a rigid rim efficiently [191]. The ANCF-RN is a nonelement node that has the position and gradient vector coordinates. An ANCF finite element node, called an interface node, can then be rigidly connected to the ANCF-RN with linear constraint equations, as in the rigid bar element in the finite element literature. Therefore, the dependent ANCF interface nodal coordinates can be eliminated in the preprocessing process, and the resulting equations of motion can be expressed in terms of the ANCF reference coordinates at the rim node and all other independent nodal coordinates of the tire. With this procedure, Patel et al. [192] developed a flexible tire model with the fully parameterized shell elements with a composite material. Recuero et al. [193] then developed a continuum deformable soil model with the ANCF solid elements for tire–soil interaction simulations with the multibody dynamics computer algorithm. Grossi et al. [194] discussed the use of the ANCF geometry representation to develop a reduced-order tire model with the floating frame of reference formulation (FFRF) for small deformation problems, including the frequency-domain analysis. This model can retain an accurate tire geometry as the ANCF mesh converted from the CAD geometry is used to define the reference configuration of the FFRF tire model. This approach was also applied to various tires [195].

Accurate tire modeling requires considering nonlinear tire friction characteristics and their coupling with the tire structural deformation. Yamashita et al. [196] integrated the distributed-parameter LuGre tire friction model into the planar ANCF tire model for the analysis of the longitudinal tire dynamics under severe braking conditions. With this procedure, the effect of the transient tire deformation on the history-dependent dynamic tire friction behavior was evaluated. The coupled ANCF-LuGre tire model was further extended to 3D problems with the shear deformable composite shell element to account for the combined transient braking and cornering scenarios [197]. Their tire model was used for off-road mobility simulations with the continuum FE and multiscale soil models [198–200], as well as discrete-element terrain model by leveraging high-performance computing [201]. These models were validated against test data. Lan et al. [202] proposed a fractional derivative viscosity constitutive model to describe the tire rubber material with the thin plate element.

Furthermore, modeling of leaf springs was examined with the ANCF approach [203,204] to account for the accurate leaf geometry in the analysis. The ANCF approach was also applied to the modeling of a flexible rover wheel [205].

A rail vehicle system has flexible components such as a pantograph-catenary system and railroad tracks, as shown in Fig. 6(b). A flexible catenary cable has been modeled by the ANCF cable element by many authors [206–212]. The equilibrium configuration of a catenary system subjected to the initial tension was discussed in Refs. [207] and [210], while control of the contact force between the pan head and contact cable was examined in Ref. [209]. Contact modeling plays an important role in pantograph-catenary simulations. The use of the sliding constraint and elastic contact models was discussed by Kulkarni et al. [212], and modeling of loss of contact between the pan head and the contact cable due to the aerodynamics effect was discussed with the elastic contact model. Using contact forces predicted by the full vehicle model, Daocharoenporn et al. [211] investigated the wear rate variation of the contact wire.

Track flexibility has been generally modeled by FFRF in multibody rail simulations as the rail deformation is small. On the other hand, Hamper et al. [213] took advantage of the ANCF geometry representation to define geometric properties of the track space curve, i.e., the tangent, normal, and their spatial derivatives, needed in the wheel-rail contact calculation. Recuero et al. [214] modeled a flexible track using ANCF and compared it with that of FFRF, showing a good agreement in small deformation problems. The ANCF flexible track model can potentially be used to predict large rail deformation caused by buckling.

Absolute nodal coordinate formulation was applied to solve challenging aerospace engineering problems as summarized in Table 3. The ANCF beam elements were used to account for the flexibility of robot arms in space as shown in Fig. 6(c), and these models were used to examine the optimal motion control and trajectory planning [215–220].

Space deployable structures have been used for truss antenna, solar panels, and solar sails. Figure 6(d) shows a conceptual diagram of a deployable truss antenna. The components of the recent deployable structures are flexible and lightweight. A large-size deployed antenna enables long-distance space observation, while a large size solar panel can generate more electricity. The ANCF approach is suited for modeling flexible deployable structures subjected to joint constraints. The deployment motions of a hoop-truss and mesh antennas were analyzed using the 3D fully parameterized beam elements [221–223] and the cable element [224]. The deployment motion of solar panels was examined using the 2D Euler–Bernoulli beam element [225,226] and the 3D fully parameterized beam element [227]. The ANCF plate elements were used to predict the deployment behavior of large deformable solar panels [228,229] and sails [230–232].

Space debris, such as artificial satellites after the end of the use and fragments of spacecraft, is a fatal problem in space engineering as the space debris causes serious damage to other spacecraft under operation when colliding with each other. A net system to capture and remove the debris is one of the promising ways to solve the problem, as shown in Fig. 6(e). ANCF is suited for modeling such a very flexible net undergoing large deformation. Shan et al. [233–235] conducted the analysis of the deployment process of a debris removal net using the cable element. The simulation results were compared with the experimental results obtained under a microgravity condition. In the comparison, they concluded that the ANCF model is superior to the conventional mass-spring model in terms of accuracy, while the computational cost of the ANCF model is higher. An orbital debris removal system was modeled by Standnyk and Ulrich [236] with the ANCF approach and validated against the test data.

Much larger scale space structures, such as space solar power stations or space solar power systems, have been analyzed by ANCF. In such large space structures, the gravitational force is no longer constant and varies depending on altitude. Thus, Sun et al. [237] modeled a tether transportation system with a moving mass along the tether using the 2D Euler–Bernoulli beam element and considered the effect of the altitude on the gravitational force. The ALE-ANCF beam element was used by Luo et al. [238] to model a variable length tether. Li et al. [239] proposed a modeling method for space structures under altitude-dependent gravitational forces using ANCF elements. Li et al. [240] investigated the gravity-gradient-induced vibration that occurs on large-scale spacecraft. The coupled orbit-attitude-vibration of large-scale space structures was examined with the ANCF elements in Refs. [241–243].

The recent trend of commercial jets is to have a high-aspect-ratio, lightweight, and very flexible wing to enable more eco-friendly flights. Such a high-aspect-ratio wing can effectively reduce the induced drag caused by wingtip vortices. However, the wing is too large to fit in the existing airport gate. Therefore, some aeronautical organizations such as Airbus, Boeing, and NASA have developed folding wing concepts as schematically shown in Fig. 6(f). In order to model a flexible folding wing subjected to mechanical joints, the ANCF modeling approach was used by Otsuka and Makihara [163]. The torsional deformation of the wing plays an important role as it drastically changes aerodynamic forces, and it is the main source of aeroelastic instability such as flutter and divergence. They introduced an infinitesimal torsional angle into the 2D ANCF beam element to consider the aeroelastic bending-torsional coupling and examined the aeroelastic vibration [244], wing deployment motion [164], and flutter instability [165]. Although the direct use of the torsional angle enabled the model to capture the aeroelastic bending–torsional coupling, its practical use was limited because of the planar model. Therefore, Otsuka et al. [26] proposed a 3D asymmetrically gradient-deficient element for modeling a high-aspect-ratio wing, considering the aerodynamics effect with UVLM. The simulation results agreed well with the wind tunnel test data.

Flexible rotating structures such as rotor blades subjected to the centrifugal stiffening effect have been investigated with the beam [245] and plate [246] elements. A wind turbine blade was modeled by the thin plate element [162]. A belt drive system experiences large rigid body motion and large deformation, and the ANCF is one of the promising modeling procedures for such an application [247]. Furthermore, a cable-pulley system such as traveling elevator cables were modeled by the ANCF elements and validated against test data in Refs. [248–250]. The ALE-ANCF beam element was applied to these applications [151,152], as discussed in Sec. 6.3. The fluid conveying pipe exhibits unstable flexible motion, and the 2D Euler–Bernoulli beam element was employed for this application [251]. A gantry crane cable was modeled with the 2D gradient-deficient beam element to examine motion control algorithms [252]. A remotely operated underwater vehicle is used for ocean survey, and it is connected to a mother ship by a flexible cable. Such a marine cable was modeled by the 3D fully parameterized beam element [253], the cable element [254,255], and the higher order beam element [256]. The piezo-electric element is useful to reduce the vibration or harvest the vibration energy of flexible structures. A piezo-electric material model was implemented in the 36-DoF thin plate element in Ref. [257]. Other applications includes the ANCF modeling of lithium-ion battery [258,259], human assistant device [260], drilling system [261], and a simple double pendulum [262], a suspension cable for bridges [263], a filament bundle [264], textile [265], moving yarn [266], and a parallel mechanism [267]. Table 4 summarizes these applications.

In this study, a comprehensive review of theoretical developments and applications of ANCF from 2012 to 2020 was presented to overview the state-of-the-art ANCF simulation capabilities that are not covered in the first review article of ANCF published in 2013. The ANCF element library has grown substantially for the beam, plate/shell, solid elements, as summarized in this paper. In particular, a number of higher order elements and shear-deformable gradient-deficient beam and plate/shell elements were developed to eliminate drawbacks of ANCF elements developed earlier. Additionally, various types of solid elements were proposed and used for the analysis of soft materials and fluid (liquid sloshing). The B-spline geometry representation in CAD can be converted directly to the ANCF geometry for the analysis through a linear transformation. This unique property in ANCF was examined to address limitations of the B-spline/NURBS representation as an analysis tool. One can facilitate structural optimization of flexible multibody systems, such as topology optimization, by exploiting the accurate geometry description of ANCF. Furthermore, computational improvements and multiphysics simulations, especially the fluid–structure interaction, have become major research topics for ANCF in recent years. There has been an increasing number of studies in aerospace applications, including the modeling of large-scale flexible deployable space structures and folding wings of aircraft.

Despite the significant developments and improvements of the ANCF simulation capabilities, there are some open issues that should be addressed in the future. Examples include, but are not limited to, (1) accurate modeling of smart materials, such as piezo-electric materials, magnetostrictive materials, dielectric elastomers, etc. for next-generation intelligent structures with application to shape control, energy harvesting, and structural health monitoring; (2) appropriate reduced order modeling and computational speedup for ANCF for modern control design of large deformable structures; and (3) more practical engineering applications of ICADA with general-purpose software to facilitate the design process.

The first and second authors acknowledge financial supports from Japan Society for the Promotion of Science (JSPS) (20K22378, 20K21041, and 21K14341); Casio Science Promotion Foundation; and Mazda Foundation. The third author acknowledges the financial support from the Automotive Research Center (ARC) in accordance with Cooperative Agreement W56HZV-19-2-0001 U.S. Army CCDC Ground Vehicle Systems Center (GVSC).

• Japan Society for the Promotion of Science (JSPS) (Grant Nos. 20K22378, 20K21041, and 21K14341; Funder ID:10.13039/501100001691).

• Casio Science Promotion Foundation (Casio Research Funding; Funder ID:10.13039/501100005936).

• Mazda Foundation (Mazda Research Funding; Funder ID: 10.13039/100002213).

• Automotive Research Center (ARC) in Accordance with Cooperative U.S. Army CCDC Ground Vehicle Systems Center (GVSC) (Agreement No. W56HZV-19-2-0001; Funder ID: 10.13039/100008192).