1. Introduction

The earliest method discovered for dealing with flexible body dynamic equations is the kineto-elasto-dynamic (KED) method [1-3]. This method is based on assumptions of low speed and small displacement, while the coupling terms of rigid-body motion and elastic motion in the dynamic equations are neglected. The inertia characteristics of flexible body rigid motion are loaded onto the flexible body itself in the form of inertial forces. With the constant emergence of lightweight, high-speed machinery, the disadvantages of the KED method have been gradually exposed [4]. To describe the effect of elastic deformation of the flexible body on its own large overall motion, Likins [5] proposed the floating frame of reference (FFR). The FFR method decomposes the configuration of the flexible body into two parts: the large overall rigid motions of the floating coordinate system, and the elastic deformation of the flexible body with respect to the floating coordinate system. Although the selection of the floating coordinate system does not affect the analysis results, the mass matrix obtained by this method is a nonlinear function matrix of generalized coordinates, resulting in inertial coupling between the rigid motion and elastic deformation of the flexible body.

To simplify the mass matrix, Simo [6-9] suggested suppressing the floating coordinate system and expressing the beam's dynamic equations in the global coordinate system directly. In describing the bending and shearing of the beam, Simo reserved the beam's cross-section local coordinate system containing the cross-section angle relative to the beam rigid configuration. Simo's approach simplifies the expression of the kinetic energy of the beam compared to the FFR method. The constant mass matrix can be obtained using variable separation; however, the expressions of elastic potential energy and the stiffness matrix become more complex.

At the end of the 20th century, Shabana and other scholars [10-14] proposed the absolute nodal coordinate formulation (ANCF) method. This method is similar to that proposed by Simo: the dynamic equations of the beam are described in the global coordinate system; the cross-section local coordinate system of the beam is used to describe bending, shearing, and twisting; and the constant mass matrix can be obtained, so the centrifugal force and Coriolis force are zero. The shortcomings of these two methods are identical: the stiffness matrix becomes highly complicated.

In contrast to Simo's method, in the ANCF approach, the cross-section slopes at the nodes (instead of the cross-section angles) become the beam's generalized coordinates. Therefore, the rigid-body inertia of the beam can be described accurately, and the constant mass matrix can be obtained without needing to perform variable separation on generalized coordinates.

For the flexible body dynamic simulation, however, a simple mass matrix and simple stiffness matrix cannot be obtained simultaneously. The stiffness matrix obtained with the FFR method is relatively simple, but the mass matrix is a nonlinear function matrix of the generalized coordinates. If the floating coordinate system is suppressed and the dynamic equations are expressed in the absolute coordinate system, then the mass matrix can be reduced to a constant matrix, but the stiffness matrix becomes highly nonlinear.

In this paper, based on assumptions of low speed and small deformation, the ANCF method is considered a finite element interpolation method. The interpolation matrices, mass matrix, and stiffness matrix of the two-node beam element are extended to those of the multi-node beam element, thereby improving calculation accuracy. The ANCF method is combined with the FFR method (for convenience, the combination of these two methods is termed the ANCF-FFR method). The stiffness matrix, which does not include generalized coordinates, is obtained along with the constant mass matrix. In the ANCF method, the bending moment is a nonlinear function vector of the generalized coordinates, so it is difficult to solve the system states directly. The split-iteration method is used to solve the system states and Lagrange multipliers. This algorithm guarantees a sufficiently accurate solution and improves the computational efficiency substantially.

2. Modeling and solving methods for flexible beam finite rotation

To improve simulation accuracy, the expressions of the generalized coordinates, interpolation matrices, mass matrix, stiffness matrix, and generalized forces of the two-node beam are extended to those of the multi-node beam. The FFR method is combined with the ANCF method to obtain the stiffness matrix, which does not contain the generalized coordinates. To enhance the efficiency of the simulation, the splitting iteration method is used to solve both the static and dynamic equations.

2.1. Generalized coordinates, interpolation matrix, and mass matrix

Fig. 1 displays the schematic diagram of the global coordinate system and the local coordinate system based on the ANCF method. For convenience, the number of nodes on the neutral axis of the beam is set to 4. The global coordinate system is (O, [r.sub.1], [r.sub.2], [r.sub.3]), and the local coordinate system is ([A.sub.1], x, y, z). The original point of the local coordinate system is located at the beam's neutral axis endpoint, [A.sub.1]; the coordinate x is the arc length coordinate of the beam neutral axis; [A.sub.1]y is the axis along the height direction of the beam cross-section; and [A.sub.1]z is the axis along the width direction of the beam cross-section.

Let the mass density of the rectangular beam be [rho], the height of the cross-section be h, and the width of the cross-section be w. The four nodes on the beam's neutral axis are [A.sub.1], [A.sub.2], [A.sub.3], [A.sub.4], which divide the beam into three elements. The length of the neutral axis of the ith element is [L.sub.i] (i=1, 2, 3). The generalized coordinates of each node are defined as:

[mathematical expression not reproducible] (1)

r = [([r.sub.1], [r.sub.2], [r.sub.3]).sup.T], i = 1, 2, 3, 4.

The generalized coordinates of the beam can be expressed as:

[mathematical expression not reproducible] (2)

set

[mathematical expression not reproducible] (3)

and the scale functions are:

[f.sub.L](a,b,c) = 2[a.sup.3]-3[a.sup.2]+1; [f.sub.xL](a, b ,c) = a - 2[a.sup.2]+[a.sup.3]; [f.sub.yL](a,b,c) = b - b * a; [f.sub.zL](a,b,c) = c - c * a. (4)

[f.sub.R](a,b,c) = 3[a.sup.2]-2[a.sup.3]; [f.sub.xR] (a, b, c ) = [a.sup.3]-[a.sup.2]; [f.sub.yR](a,b,c) = b-a; [f.sub.zR](a,b,c) = c-a. (5)

The interpolation functions of the ith element are defined as:

S[L.sub.i] = [f.sub.L]([[xi].sub.i] ,[[eta].sub.i],[[zeta].sub.i]); SL[x.sub.i] = [L.sub.i]*[f.sub.xL]([[xi].sub.i],[[eta].sub.i],[[zeta].sub.i]); SL[y.sub.i]= h*[f.sub.yL]([[xi].sub.i],[[eta].sub.i],[[zeta].sub.i]); SL[z.sub.i]= w * [f.sub.zL]([[xi].sub.i],[[eta].sub.i],[[zeta].sub.i]); i=1, 2, 3.

S[R.sub.i +1] = [f.sub.R]([[xi].sub.i],[[eta].sub.i],[[zeta].sub.i]); SR[x.sub.i +1] =[L.sub.i]*[f.sub.xR]([[xi].sub.i],[[eta].sub.i],[[zeta].sub.i]); SR[y.sub.i +1] =h*[f.sub.yR]([[xi].sub.i],[[eta].sub.i],[[zeta].sub.i]); R[z.sub.i +1]= w*[f.sub.zR] ([[xi].sub.i],[[eta].sub.i],[[zeta].sub.i]); (7) i=1, 2, 3.

Define matrices:

S[L.sub.i] = [[SLi * I SLxi * I SLyi * I SLzi * I].sub.3 x 12]; i = 1, 2, 3; SRj = [[S[R.sub.i] * I SR[x.sub.i] * I SRyi * I SR[z.sub.i] * I].sub.3 x 12]; j = 2, 3, 4; (8)

where: I is the 3-order identity matrix. The interpolation matrix of the beam elements between two adjacent nodes [A.sub.i] , [A.sub.i +1] can be expressed as:

[S.sub.1] = [[S[L.sub.1] S[R.sub.2] [0.sub.3x12] [0.sub.3x12]].sub.3x48]; [S.sub.2] = [[[0.sub.3x12] S[L.sub.2] S[R.sub.3] [0.sub.3x12]].sub.3x48]; [S.sub.3] = [[[0.sub.3x12] [0.sub.3x12] S[L.sub.3] S[R.sub.4]].sub.3x48]; (9)

The symbol 0 means the zero matrix, and the elements of [S.sub.i] (i= 1, 2, 3) are the functions of x, y, z. Let P be an arbitrary point on the beam element between nodes Ai and [A.sub.i +1]. The position vector of P in the global coordinate system can be expressed as follows:

r[|.sub.P] =Si(x,y,z)[|.sub.P]-e ; i=1, 2, 3. (10)

Then, the kinetic energy of the beam is:

[mathematical expression not reproducible] (11)

and the mass matrix of the beam is:

[mathematical expression not reproducible] (12)

where: M is a symmetric positive definite constant matrix.

2.2. Strain energy and stiffness matrix

The generalized displacement produced by elastic deformation at point P on the ith (i=1, 2, 3) beam element can be expressed as:

u[|.sub.P]=[S.sub.i](x, y, z)[|.sub.P]-(e - [e.sub.0]); i = 1, 2, 3. (13)

where: [e.sub.0] is the generalized coordinate vector of the beam in its rigid configuration. The deformation gradient of P is:

[mathematical expression not reproducible] (14)

where:

[mathematical expression not reproducible] (15)

Given low-speed and small-deformation assumptions and the FFR method, the ANCF local coordinate system can be considered the float coordinate system in the FFR method, i.e.,

[mathematical expression not reproducible] (16)

where: [[psi].sub.1], [[psi].sub.2], [[psi].sub.3] represent the Caldan angles. The Cauchy strain tensor of point P is:

[mathematical expression not reproducible] (16)

and its vector form is:

[mathematical expression not reproducible] (17)

The elastic matrix Q is:

where: E is the elastic modulus and [mu] is the Poisson's ratio. To eliminate the Poisson locking phenomenon, we set [mu]=0; thus, the elastic matrix Q can be expressed as:

Q = diag(E E E E/2 E/2 E/2). (18)

The elastic strain energy of the ith beam element is:

[mathematical expression not reproducible] (19)

The elastic strain energy of the beam is:

[mathematical expression not reproducible] (19)

[mathematical expression not reproducible] (20)

Based on the assumptions of low speed and small deformation, geometric nonlinearity is ignored, and the tangent stiffness matrix of the beam can be approximated as:

[mathematical expression not reproducible] (21)

where: K is a symmetric positive semi-definite matrix that does not include the generalized coordinates.

In summary, for low-speed, small-deformation problems, we can combine the ANCF method and FFR method to obtain the constant mass matrix and stiffness matrix, which does not include the generalized coordinates.

2.3. Gravity

The virtual displacement by gravity at point P on the element between two adjacent nodes [A.sub.i], [A.sub.i+1] is set to [delta]e, and then the virtual work done by gravity at point P is:

[mathematical expression not reproducible] (22)

where: g denotes gravitational acceleration. From that, we can obtain:

[mathematical expression not reproducible] (23)

where: V is the beam volume. The gravity on the beam element between [A.sub.i], [A.sub.i +1] is :

[mathematical expression not reproducible] (24)

and the overall gravity of the beam is:

[mathematical expression not reproducible] (25)

2.4. Bending moment

As shown in Fig. 2, a bending moment [M.sub.B] is placed at the cross-section of point [A.sub.2], its magnitude is M.

The virtual work of the bending moment [M.sub.B] is:

[mathematical expression not reproducible] (26)

From [15-16], the rotation angle [theta] of the cross-section satisfies the following:

[mathematical expression not reproducible] (27)

so

[mathematical expression not reproducible] (28)

and the bending moment [M.sub.B] can be expressed as:

[mathematical expression not reproducible] (29)

2.4. Boundary conditions

Let the beam rotate anticlockwise around the O[r.sub.3] axis. A constant bending moment vector [M.sub.B] is placed at the cross-section of point [A.sub.2]. We set the boundary conditions in the rotating process as follows:

a. At point [A.sub.1], the vector [mathematical expression not reproducible] is parallel with the vector [([r.sub.1], [r.sub.2], [r.sub.3]).sup.T.sub.A1]:

[[PI].sub.1] (t) = [e.sub.1] (t) * [e.sub.8] (t) - [e.sub.2] (t) * [e.sub.7] (t) = 0; (30)

b. The distance between point [A.sub.1] and O(0, 0, 0) is a constant:

[[PI].sub.2] (t) = [e.sub.1.sup.2] (t) + [e.sub.2.sup.2] (t) + [e.sub.3.sup.2] (t) [equivalent to] [H.sup.2]; (31)

c. The coordinate [r.sub.3] of point A1 is zero:

[[PI].sub.3](t) = [e.sub.3](t) [equivalent to] 0; (32)

d. The coordinate [partial derivative][r.sub.3]/[partial derivative]y of point A1 is zero:

[[PI].sub.4](t) = [e.sub.9](t ) [equivalent to] 0; (33)

e. Let the angle displacement function [phi](t) is a known function. The coordinate [r.sub.1] of point [A.sub.2] satisfies:

[mathematical expression not reproducible] (34)

f. The coordinate [r.sub.2] of point [A.sub.2] satisfies:

[mathematical expression not reproducible] (35)

g. The coordinate [r.sub.3] of point A2 is zero.

[[PI].sub.7](t ) = [e.sub.15](t) [equivalent to] 0. (36)

Fig. 3 is the schematic diagram of finite rotation boundary conditions of the beam.

2.5. Solution of initial states

The static equilibrium equation of the system at the moment of t=0 can be expressed as:

[mathematical expression not reproducible] (37)

where: [e.sub.0] is the rigid configuration of the beam, [lambda] is the vector of the Lagrange multipliers, and [[PI].sub.e] is defined as:

[[PI].sub.e]=([partial derivative][[PI].sub.j]/[partial derivative][e.sub.i]) (38)

For the expression of generalized force vector:

[F.sub.Gen](e) = G + M + [[PI].sub.e] * [lambda] + K * [e.sub.0], (39)

contains the nonlinear terms of the generalized coordinates and the coupling terms between the generalized coordinates and Lagrange multipliers, it is difficult to solve the system states directly. The split-iteration algorithm is thus used to solve Equation (37).

For convenience, assume the following conventions: Let [N.sub.1] and [N.sub.2] be two sets such that:

[N.sub.1] ={1, 2, 3, 9, 13, 14, 15}; [N.sub.1] [union] [N.sub.2] ={1, 2, 3, ..., 48}; (40) [N.sub.1] [intersection] [N.sub.2]=0.

Set matrix A = [([a.sub.ij]).sub.s x t], vector [mathematical expression not reproducible] is the submatrix of A, v([[OMEGA].sub.1],...,[[OMEGA].sub.l]) is the subvector of v:

[mathematical expression not reproducible] (41)

[mathematical expression not reproducible] (42)

Using the above agreements, we set:

[mathematical expression not reproducible] (43)

[mathematical expression not reproducible] (44)

[mathematical expression not reproducible] (45)

[mathematical expression not reproducible] (46)

[mathematical expression not reproducible] (47)

[mathematical expression not reproducible] (48)

[mathematical expression not reproducible] (49)

[mathematical expression not reproducible] (50)

where: [K.sup.R.sub.R] and [K.sup.D.sub.D] are invertible matrices. Equation (37) can be expressed as:

[mathematical expression not reproducible] (51)

set

[mathematical expression not reproducible] (52)

where: [K.sup.R.sub.R1] and [K.sup.R.sub.R2] satisfy:

[K.sup.R.sub.R1] + [K.sup.R.sub.R2]= [K.sup.R.sub.R], (53)

so

[mathematical expression not reproducible] (54)

The first equation of (51) can be written as:

[mathematical expression not reproducible] (55)

and we can get the iterative algorithm:

[mathematical expression not reproducible] (56)

Equation (58) is convergent when p[epsilon][0,1). If the elastic modulus of the material is sufficiently large, then the fast convergence rate can be obtained when p=0.

Because the generalized force vector [F.sub.R] is a nonlinear function vector of [e.sub.R] and [lambda], the right side of Equation (56) should be simplified before each iteration. For small deformation problems, expanding the right side of (56) to the 2nd-order Taylor series at [e.sub.D]=[e.sub.Dr], [lambda]=0 before each iteration can ensure sufficient accuracy of the solution, where [e.sub.Dr] is:

Finally, [e.sub.R] converges to the 2nd-order Taylor series of [e.sub.D] and [lambda].:

[e.sub.R] = [e.sub.R] ([e.sub.D], [lambda]). (57)

After that, the constraint Eqs: (30)-(36) and the second equation in (51) should be linearized to the 1st-order Taylor series at [e.sub.D]=[e.sub.Dr], [lambda]=0, and combined with Equation (57), [e.sub.R], [e.sub.D] and [lambda] can be solved quickly.

2.6. Solution of dynamic equations

The dynamic equations of the beam can be expressed as:

[mathematical expression not reproducible] (58)

Let the damping matrix C = 0, and Equation (58) can be solved with the Newmark method:

[mathematical expression not reproducible] (59)

[mathematical expression not reproducible] (60)

[mathematical expression not reproducible] (61)

where: [beta] =0.25, [gamma]=0.5. After rearranging Eqs. (59)-(61), we can get:

[mathematical expression not reproducible] (62)

Eq. (62) can be solved using the method described in Section 2.6.

3. Numerical example

Set the length of the rectangular section beam as 7.791m. The position of each node on the beam is shown in Figure 3.the width of the cross-section as w=0.1m, the height of the cross-section as h=0.1m, the elastic modulus as E=2.1 x [10.sup.11]Pa, the mass density as [rho]=7850kg/[m.sup.3], and Poisson's ratio as [mu]=03. The boundary conditions are described in Section 2.5, and O[A.sub.1] satisfies the following:

||O[A.sub.1]|| [equivalent to] H= 0.295m. (63)

The lump mass m is fixed at the endpoint of the beam. The constant bending moment [M.sub.B] is loaded at the cross-section of node [A.sub.2], and:

||[M.sub.B]|| = 1000Nm. (64)

The beam angular displacement function[phi](t) satisfies:

[mathematical expression not reproducible] (65)

Figs. 4 and 5 denote the static deflection of the end-point of the neutral axis of the beam obtained by the ANCF-FFR method and ABAQUS simulation when the lump mass m = 0 kg and m = 50 kg, respectively. The angles between the beam rigid configuration and the O[r.sub.1] axis are set at 0 degrees, 15 degrees, 30 degrees, 45 degrees, 60 degrees, respectively. When using the ANCF-FFR method, the numbers of nodes on the neutral axis of the beam are 4, 10, 16, 22, and 28, respectively. In the ABAQUS software, the beam model is divided into 4192 2nd-order hexahedral elements. The charts and tables indicate that when the number of nodes on the neutral axis of the beam is small, the results obtained by the ANCF-FFR are quite different than those in the ABAQUS simulation. As the number of nodes increases, the ANCF-FFR results tend to be consistent with those of ABAQUS. When the number of nodes in ANCF reaches 28, the relative error of the two methods is less than 1%, a satisfactorily accurate proportion for engineering machinery such as cranes and aerial work platforms.

Let the initial state of the beam be the horizontal static equilibrium state. Figs. 6 and 11 display the dynamic deflections of the beam during the starting stage, luffing stage, braking stage, and steady-state vibrating stage obtained by the ANCF-FFR method and ABAQUS software when the lump mass m = 0 kg and m = 50 kg, respectively. Figs. 7-10 and Figs. 12-15 are the local charts of every stage when m = 0 kg and m = 50 kg, respectively. From the starting stage to the steady-state vibrating stage after braking, the results of the 4-node beam clearly deviate from other results. The results of the 10-node beam agree well with the ABAQUS simulation results during the starting stage, but subsequent phase errors increase significantly over time during the luffing stage and the steady--state vibrating stage. The 16-node results show good pre-braking performance. In the steady-state vibrating stage, however, the phase errors between the 16-node results and the ABAQUS results are slightly larger. The 22-node and 28-node results are almost completely aligned with the ABAQUS simulation results. Of course, an increase in accuracy requires an increase in computing time.

4. Conclusions

1. Based on the assumptions of low speed and small deformation, the constant mass matrix and stiffness matrix without generalized coordinates of the multi-node beam are obtained by combining the ANCF and FFR methods; the calculation process is simplified, and the solution accuracy is improved.

2. By using the split-iteration method, the generalized coordinates not included in the constraint equations are iteratively expanded as a 2nd-order Taylor series of the Lagrange multipliers and the generalized coordinates contained in the constraint equations, and the constraint equations themselves, are linearized to a 1st-order Taylor series at the rigid configuration of the generalized coordinates, after which the system states and Lagrange multipliers can be solved quickly. The simulation results demonstrate that for small-deformation problems, the low-order Taylor approximation of the generalized displacement at the rigid configuration can significantly improve the solution speed and ensure sufficient solution accuracy.

Acknowledgements

The financial supports from the Natural Science Foundation of Jiangsu Province, China (Grant No. BK20161522), Six Talent Peaks Project in Jiangsu Province, China (Grant No. JXQC-023), the State Laboratory of Advanced Design and Manufacturing for Vehicle Body (Grant No. 31415008) and Chinese Postdoctoral Science Foundation (Grant No. 2015M572243).

References

[1.] Turcic, D. A.; Midha, Ashok. 1984. Dynamic analysis of elastic mechanism systems, Part I: Applications 106(106): 249-254. https://doi.org/10.1115/1.3140681.

[2.] Turcic, D. A.; Midha, Ashok; Bosnik, J. R. 1984. Dynamic analysis of elastic mechanism systems, Part II: Experimental Results 106(106): 255-260. https://doi.org/10.1115/1.3140682.

[3.] Turcic, D. A.; Midha, Ashok. 1984. Generalized equations of motion for the dynamic analysis of elastic mechanism systems 106(4): 243-248. https://doi.org/10.1115/1.3140680.

[4.] Qiang Tian. 2010. Advances in the absolute nodal coordinate method for the flexible multibody dynamics, 40(02): 189-202. https://doi.org/0.6052/1000-0992-2010-2-J2009-024.

[5.] Likins, P. W. 1972. Finite element appendage equations for hybrid coordinate dynamic analysis, 8(5): 709-731. https://doi.org/10.1016/0020-7683(72)90038-8.

[6.] Simo, J. C.; Vu-Quoc, L. 1986. On the dynamics of flexible beams under large overall motions--the plane case: Part I, 53(4): 849-854. https://doi.org/10.1115/1.3171870.

[7.] Simo, J. C.; Vu-Quoc, L. 1986. On the dynamics of flexible beams under large overall motions--the plane case: Part II, 53(4): 855-863. https://doi.org/10.1115/1.3171871.

[8.] Simo, J. C.; Vu-Quoc, L. 1986. A three-dimensional finite-strain rod model. part II: Computational aspects, 58:79-116. https://doi.org/10.1016/0045-7825(86)90079-4.

[9.] Simo, J. C. 1985. A finite strain beam formulation. The three-dimensional dynamic problem. Part I, 49: 55-70. https://doi.org/10.1016/0045-7825(85)90050-7.

[10.] Shabana, A. 1996. Finite element incremental approach and exact rigid body inertia, 118(2): 171-178. https://doi.org/10.1115/1.2826866.

[11.] Shabana, A. 1997. Definition of the slopes and the finite element absolute nodal coordinate formulation, 1(3): 339-348. https://doi.org/10.1023/A:1009740800463.

[12.] Shabana, A.; Hussien, H. A.; Escalona, J. L. 1998. Application of the absolute nodal coordinate formulation to large rotation and large deformation problems, 120(2): 188-195. https://doi.org/10.1115/1.2826958.

[13.] Ahmed, A.; Shabana, Refaat Y. Yakoub. 2001. Three dimensional absolute nodal coordinate formulation for beam elements: theory, 123(4): 606-613. https://doi.org/10.1115/1.1410100.

[14.] Refaat Y. Yakoub, Ahmed A. Shabana. 2001. Three dimensional absolute nodal coordinate formulation for beam elements: implementation and applications, 123(4): 614-621. https://doi.org/10.1115/1.1410099.

[15.] Ahmed A. Shabana. 2005. Dynamics of multibody systems. Cambridge University Press, 312 p.

[16.] Omar, M. A.; Shabana, A. A. 2001. A two-dimensional shear deformable beam for large rotation and deformation problems, 243(3): 565-576. https://doi.org/10.1006/jsvi.2000.3416.

Ling HAN (*), Ying LIU (**), Bin YANG (***), Yueqin ZHANG (****)

(*) Nanjing Forestry University, Nanjing 210037, China, E-mail: hanlingdx@163.com

(**) Nanjing Forestry University, Nanjing 210037, China, E-mail: liuying@njfu.edu.cn

(***) Nanjing Institute of Technology, Nanjing, 211167, China, Email: yangb123@126.com

(****) Nanjing University of Aeronautics and Astronautics, Nanjing, 210016, E-mail: zyqnuaa@163.com

http://dx.doi.org/10.5755/j01.mech.24.5.20358

Received May 01, 2018

Accepted October 18, 2018