Application of discrete shape function in absolute nodal coordinate formulation

: To solve integrals in the absolute nodal coordinate method and address the difficulty in applying it to an arbitrary-section beam, this paper focuses on two methods involving single integrals ： the invariant matrix method and the Gerstmayr method, with cross-section characteristics by applying the interpolation of a discrete function. Such single integrals demonstrate that the nodal coordinate method can be applied to an arbitrary-section beam. The Euler–Bernoulli beam used in engineering structures is characterised by a symmetrical cross-section, small section size, zero odd integrals and negligible high-order even integrals, which simplifies the single integrals of the two methods. Finally, the Gaussian integration is adopted to improve the solving efficiency of elastic force and force Jacobian.


Introduction
Large size, high precision, lightweight and multi-shape have become the inevitable demand of future spacecraft, and flexible multi-body systems rich in large deformation components have been widely used in the aerospace field [1]. For example, spacecraft solar panels [2], satellite expandable antennas [3], space telescopes, etc. [4]. In order to obtain the mechanical mechanism of such mechanisms with large deformation and large rotation [5][6][7], much work has been done in this field.
The dynamics of the flexible multi-body system focused on the dynamic behavior of a system composed of a deformable object and rigid body when it experiences a wide range of spatial motions [8]. Turcic et al. [9] proposed Kineto Elastio Dynamic Analysis (KED) method firstly to deal with dynamics of flexible multi-body systems. This method decomposes flexible multi-body dynamics into multi-rigid body dynamics and structural dynamics provides high analytical accuracy for multi-body systems with low flexibility due to low-speed movement [10]. With the emergence of lightweight highspeed mechanism, the dynamic equation constructed by the KED method could not explain the coupling effect on the elastic deformation and large deformation caused by high-speed movements, so the floating frame of reference formulation (FFRF) has been proposed [11,12]. This method constructs the floating coordinate system on the flexible body and establishes the dynamic model by using the rigid body coordinates of a floating coordinate system and the node coordinates (or modal coordinates) with the flexible body. The geometrically exact formulation (GEF) can reasonably describe finite rotation, but it has higher requirements for interpolation strategy and numerical integration method. Shabana [13] proposed an absolute nodal coordinate formulation (ANCF) based on the finite element and continuum mechanics theory. In this method, the generalized coordinates of element nodes were defined in the global coordinate system, and slope vectors were used to replace the angular coordinates in the traditional finite element method. The dynamic equation has the characteristics of a constant inertia matrix, the absence of Coriolis force, and centrifugal force [14,15]. Compared with the traditional method of FFRF [16], the GEF [17], ANCF can describe the flexible multi-body system more accurately and plays an important role in dealing with the dynamics of flexible multi-body systems [18,19].
The ANCF is widely used in flexible multibody simulations. Yakoub [20] applied the threedimensional beam element to simulate the wheel-rail contact dynamics problem, which applied the ANCF to practical engineering firstly. Li et al. [21] simulated the antenna reflector mesh based on the slope deficient ANCF beam element and applied it to the form-finding analysis of the reflector with the antenna. Tian et al. [22] introduced the ANCF in several types of deployable space structure and proposed a parallel computational framework to improve the computational efficiency. Lan et al. [23] used ANCF to model the multilayer beam with a circular cross-section. Wang et al. [24] studied the dynamics of flexible multi-body systems with hybrid uncertain parameters by using the ANCF, to predict they are truly dynamic behaviors. This paper mainly studied the flexible beam element, and briefly reviews the beam element based on the absolute nodal coordinate formulation. Shabana et al. [25] considered a two-dimensional beam element based on ANCF, which could efficiently solve the elastic force and force Jacobi. However, it is impossible to describe the complex beam elements subjected to torsion and shear. For this reason, Omar and Shabana and Omar [26] further proposed a two-dimensional shear deformation beam element that could describe the shear effect. García-Vallejo et al. [27] pointed out that the shear lock problem existed in the two-dimensional shear deformation beam element, and a new locking-free shear deformable finite element has been put forward. To accurately describe the force characteristics of the flexible spatial beam, Shabana and Yakoub [28] proposed a three-dimensional beam element model based on the absolute nodal coordinate formulation. García-Vallejo et al. [29] proposed an invariant matrix method to efficiently formulate the elastic force and the elastic force Jacobi, pointed out the nonlinear stiffness matrix can be simplified using the invariant matrix. Based on the principle of virtual work, Gerstmayr and Shabana [30] obtained the expression of the elastic force expressed by the first Piola-Kirchhoff stress tensor and position gradient tensor, the calculation efficiency was improved. Based on Gerstmayr's research, Liu et al. [31] derived the analytical formula that could accurately describe the elastic force Jacobian matrix, further improved the solving efficiency of the elastic force Jacobian.
However, it is difficult to apply a three-dimensional ANCF beam element, on the one hand, the elastic force is derived from the volume integral of the elastic energy using Green-Lagrange strain and Lame's constants [32]. On the other hand, the ANCF fully parametritis beam element uses linear interpolation at the transverse directions, therefore, the cross-section of the beam element must keep rectangular. Therefore, what is the research emphases of this paper are to extend the applicability of the ANCF on arbitrary section beams and its succinct solution form.
The contents of the paper are as follows. Section 2 introduces the absolute nodal coordinate theory and the discrete form of its shape function and proposes a position gradient construction method that takes initial configurations into account. Section 3 describes the single-integral form of the invariant matrix obtained by applying the discrete form of the shape function in the García-Vallejo method. The other single-integral form is described in Section 4 by applying the discrete form of the shape function in the Gerstmayr method. Section 5 proposes a high-order interpolation function for the Poisson locking problem. Some numerical simulations are given in Section 6, and the correctness of the model is verified. Section 7 provides concluding remarks.

Absolute nodal coordinate formulation
In this subsection, the three-dimensional beam element for simplicity is shown in Figure 1  The co-ordinates of node are represented as, where, the vector denotes node global position vector in the beam elements, described as, And [ T T T ] T is position vector gradient, derivatives of the position vector with respect to co-ordinates , and , Because the finite element is a continuous system, the global position vector of an arbitrary point in beam elements is the function of the local parameters , and . Using the shape function as follows, And is a 3×3 identity matrix and the shape functions are defined as [30], where, represents the beam element length. The shape function of the displacement field is composed of the cubic interpolation of and linear interpolation of and in the local coordinate system. According to the theory of continuum mechanics, the position gradient tensor is used to describe the position relationship between the initial configuration of the elastomer and the current configuration after movement. The local coordinate system fixed on the material is denoted as ( = , , ). where, And, represents the generalized coordinates of the beam element nodes in the reference configuration, Since and −1 are constant, and the position gradient is written as, In order to construct the subsequent invariant matrix, the position gradient tensor should be rewritten here. García-Vallejo et al. [29] proposed a position gradient construction method to introduce the initial configuration, here, a new construction method was proposed, the position vector , can denote as, ̅ ∈ ℛ 8×1 , ̅ ∈ ℛ 3×8 , the specific form is as follows.
Therefore, the position gradient is written as follows, and ∈ ℛ 8×3 contains the reference configuration information of the beam element. According to the expression, is a constant matrix, and the Green strain tensor can be expressed as, The components of symmetric strain tensor are, , represents the th ( = 1,2,3) column of the position gradient.
Hence, the six different components of the symmetric strain tensor are, While, based on the principle of virtual work, the expression of the element elastic force can be obtained.

Application of discrete shape function to invariant matrix method
García-Vallejo pointed out that the nonlinear stiffness matrix ( ) in Eq (23) can be divided into two parts, The constant stiffness matrix 1 can be calculated and stored in the pre-processing process. The nonlinear stiffness matrix 2 ( ) , it can be solved by the invariant matrix, so the elements in the 2 ( ) matrix is expressed by the invariant matrix, By substituting the discrete form of the shape function of Eq (30), 1 and 2 ( ) can be obtained. Please refer to Appendix 2 for the specific derivation process.
In practice, the beam element commonly used in engineering is the Euler-Bernoulli beam with a symmetrical section and a small section size. Considering this situation, the cross-section integral can be simplified, that is, the odd term integral is zero, even higher-order integral tends to zero. The greatest advantage of this is that the stiffness matrix 1 , 2 ( ) is simplified to the single integral form which is only related to the area and moment of inertia of the section.
A is the section area of the beam element with arbitrary section, I is the moment of inertia around the -axis of the beam section, and I is the moment of inertia around the -axis of the beam section.
The Jacobian of the elastic force can also be derived by the invariant matrix method, And the element of the Jacobian can express as, (37)

Application of discrete shape function to Gerstmayr method
The Gerstmayr method is based on the first Piola-Kirchhoff stress tensor and position gradient tensor [31], the virtual work of the elastic force is written as, here, is the second Piola-Kirchhoff stress tensor, is the Green strain tensor, is the first Piola-Kirchhoff stress tensor, and the component of the elastic force can be expressed as, Introduce, Equation (39) can be expressed as, is Kronecker ， can be represented by and , here, represents the component of the fourth-order elastic tensor, ℐ represents the component of the second-order unit tensor, Substituting Eq (42) into Eq (41), The specific derivation process can refer to Appendix 3. When the beam element is a symmetrical cross-section and small section size, the Eq (46) can be simplified to,

Poisson locking
Since most of the interpolating polynomials into ANCF beam elements in the transverse direction only use the first-order polynomial, the cross-section is still planar, thus prone to various locking effects. What is worth mentioning is that using higher-order interpolation polynomial can avoid locking effects [ The position vector expressed in terms of the nodal coordinates is given by, ̅ ∈ ℛ 10×1 , ̅ ∈ ℛ 3×10 , the specific form is as follows. According to the research content of Sections 3 and 4, the solution forms of elastic force and force Jacobi based on the García-Vallejo method and the Gerstmayr method can be obtained by reconstructing the ℍ (as given in Appendix 4).

Numerical example
In this section, some numerical examples are used to verify the correctness of the proposed model.

Statics problem
A cantilever beam with a circular cross-section is shown in Figure 2. Geometry properties of the beam are: the length is = 1 , the radius of the circular section is = 0.01 . And its material properties are: Young's modulus E = 72 GPa, the Poisson's ratio is = 0.33.  Table 1.
According to the data in Table 1, when = 0.33, the Poisson locking leads to the wrong results from Methods I and II, while Method III is consistent with the theoretical results.

Flexible pendulum
Considering the motion trajectory of the right endpoint B of the flexible pendulum with circular cross-section under the action of gravity, the length of the pendulum is 2 , as shown in Figure 3. The beam element section is a circle with a radius of 0.02 , the elastic modulus of the material E = 720 MPa, the density is 7200 kg/m 3 , Poisson's ratio = 0 and the gravitational acceleration is 9.81 m/s 2 . In order to verify the correctness of the results, the same flexible pendulum model is established in MSC.ADAMS and the iteration steps length are set to 0.001 seconds. The simulation results are shown in Figures 3-6.       Figure 9 presents the displacement of the endpoint B of the spatial flexible double pendulum in the X-direction. Compared with the displacement results of Method I and Method II, Figure 8, Method III shows the difference. As shown in Figure 10, the velocity curve of the Modified Higher-Order Element has obvious jitter, especially in the later stage, the results show that Method III makes the beam more flexible.   Run MATLAB on a laptop with an Intel Core 2.6GHz processor and 6GB of RAM, different numbers of elements were used to divide the beam, and the CPU time necessary to solve the problem as shown in Table 2.  Table 1 shows the calculation time of the three methods. Method II has better advantages, both the pre-processing speed and the simulation speed are greatly improved. While the pre-processing speed and the simulation speed of Method III are significantly reduced due to higher-order interpolation polynomials and more nodal coordinates. What is worth mentioning is that, the invariant of the modified García-Vallejo method can be shared for elements with equal dimensions, it decreases significantly the amount of data to compute in the pre-processing stage. Figure 11 shows the spatial configuration of the flexible double pendulum.

Summary and conclusions
In this paper, the shape function describing the absolute nodal coordinate formulation is divided into three parts, which are independent of the local coordinates in the floating coordinate system, the advantage of this is that the absolute nodal coordinate method can be applied to an arbitrary-section beam via the discrete form of the shape function. Moreover, considering the Euler-Bernoulli beam used in engineering structures is characterised by a symmetrical cross-section, small section size, zero odd integrals, and negligible high-order even integrals, hence the single-integral form reduces the difficulty of applying the absolute nodal coordinate method. Furthermore, single integrals are more efficient and improve pre-processing speed compared to volume integrals. Finally, the correctness of the proposed methods is verified through simulation examples. By comparing the computational costs of different methods, modified Gerstmayr Method obvious advantages in pre-processing speed and computational efficiency. Discrete shape function, The partial derivative with respect to local coordinates, The discrete form of , The specific expression form of , , are,