Dynamic Modeling and Simulation of Flexible Beam Finite Rotation with ANCF Method and FFR Method

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.


Introduction
The earliest method discovered for dealing with flexible body dynamic equations is the kineto-elasto-dynamic (KED) method [1][2][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][7][8][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][11][12][13][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.

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 twonode 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.

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, r1, r2, r3), and the local coordinate system is (A1, x, y, z).The original point of the local coordinate system is located at the beam's neutral axis endpoint, A1; the coordinate x is the arc length coordinate of the beam neutral axis; A1y is the axis along the height direction of the beam cross-section; and A1z is the axis along the width direction of the beam cross-section.
Fig. 1 The coordinate systems of the ANCF method Let the mass density of the rectangular beam be ρ, 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 A1, A2, A3, A4, which divide the beam into three elements.The length of the neutral axis of the ith element is Li (i=1, 2, 3).The generalized coordinates of each node are defined as: ; The generalized coordinates of the beam can be expressed as: and the scale functions are: , , , , .The interpolation functions of the ith element are defined as: , , ; , , ; , , ; , , ; 1 , 2 , 3. ) , , ; Define matrices: ; where: I is the 3-order identity matrix.The interpolation matrix of the beam elements between two adjacent nodes Ai, Ai+1 can be expressed as: The symbol 0 means the zero matrix, and the elements of Si (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 Ai+1.The position vector of P in the global coordinate system can be expressed as follows: ( , , ) ; 1 , 2 , 3.
Then, the kinetic energy of the beam is:   e e SS (11) and the mass matrix of the beam is: ( ) where: M is a symmetric positive definite constant matrix.

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: 1 , 2 , 3.
where: e0 is the generalized coordinate vector of the beam in its rigid configuration.The deformation gradient of P is: where: ( ) 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., where: ψ1, ψ2, ψ3 represent the Caldan angles.The Cauchy strain tensor of point P is: and its vector form is: ) The elastic matrix Q is: where: E is the elastic modulus and μ is the Poisson's ratio.
To eliminate the Poisson locking phenomenon, we set μ=0; thus, the elastic matrix Q can be expressed as: The elastic strain energy of the ith beam element is: The elastic strain energy of the beam is: .
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: 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.

Gravity
The virtual displacement by gravity at point P on the element between two adjacent nodes Ai, Ai+1 is set to δe, and then the virtual work done by gravity at point P is: (0 , , 0) (0 , , 0) , where: g denotes gravitational acceleration.From that, we can obtain: where: V is the beam volume.The gravity on the beam el- ement between Ai, Ai+1 is: and the overall gravity of the beam is:

Bending moment
As shown in Fig. 2, a bending moment MB is placed at the cross-section of point A2, its magnitude is M. Fig. 2 The bending moment at cross-section of point A2 The virtual work of the bending moment MB is: .
From [15][16], the rotation angle θ of the cross-section satisfies the following:

Boundary conditions
Let the beam rotate anticlockwise around the Or3 axis.A constant bending moment vector MB is placed at the cross-section of point A2.We set the boundary conditions in the rotating process as follows: a.At point A1, the vector g.The coordinate r3 of point A2 is zero.The static equilibrium equation of the system at the moment of t=0 can be expressed as: ( ) where: e0 is the rigid configuration of the beam, λ is the vector of the Lagrange multipliers, and Πe is defined as: For the expression of generalized force vector: 0 ( ) , 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 N1 and N2 be two sets such that:  ( ) ; = , , .

RD
( ) where: K R R1 and K R R2 satisfy: The first equation of ( 51) can be written as: and we can get the iterative algorithm: Equation ( 58) is convergent when p∈[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 FR is a nonlinear function vector of eR and λ, the right side of Equation ( 56) should be simplified before each iteration.For small deformation problems, expanding the right side of (56) to the 2 nd -order Taylor series at eD=eDr, λ=0 before each iteration can ensure sufficient accuracy of the solution, where , , 0 , , 0 + , + , 0 .

Dr Hsin Hcos A A H cos
Finally, eR converges to the 2 nd -order Taylor series of eD and λ.: After that , the constraint Eqs: (30)-( 36) and the second equation in (51) should be linearized to the 1 st -order Taylor series at eD=eDr, λ=0, and combined with Equation (57), eR, eD and λ can be solved quickly.

Solution of dynamic equations
The dynamic equations of the beam can be expressed as: Let the damping matrix C=0, and Equation (58) can be solved with the Newmark method: Eq. ( 62) can be solved using the method described in Section 2.6.

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×10 11 Pa, the mass density as ρ=7850kg/m 3 , and Poisson's ratio as μ=0.3.The boundary conditions are described in Section 2.5, and OA1 satisfies the following: The lump mass m is fixed at the endpoint of the beam.The constant bending moment MB is loaded at the cross-section of node A2, and: The beam angular displacement functionφ(t) satisfies: Figs. 4 and 5 denote the static deflection of the endpoint 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 Or1 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 2 nd -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.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 2 nd -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 1 st -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.
Ling HAN, Ying LIU, Bin YANG, Yueqin ZHANG S u m m a r y

DYNAMIC MODELING AND SIMULATION OF FLEXIBLE BEAM FINITE ROTATION WITH ANCF METHOD AND FFR METHOD
In this paper, a method of modeling and simulating flexible beam finite rotation is investigated.Based on the assumptions of low speed and small deformation, the ANCF method is regarded as a finite element interpolation method to obtain the constant mass matrix of the flexible beam; the local coordinate system in the ANCF method is considered the floating coordinate system, and the stiffness matrix independent of the generalized coordinates is obtained.The split-iteration method is used to expand the generalized coordinates that are not contained in the constraint equations to the 2 nd -order Taylor series of the generalized coordinates that are contained in the constraint equations and the Lagrange multipliers.The nonlinear constraint equations are linearized to the 1 st -order Taylor series of the generalized coordinates.Then, the generalized coordinates and Lagrange multipliers can be solved quickly.The results show that the dynamic equations can be effectively simplified by combining the ANCF method with the FFR method for the small-deformation problems.The low-order Taylor approximation of generalized coordinates in both the dynamic equations and constrained equations does not lose substantial computational accuracy but can significantly reduce computational time.The results of this investigation have important reference values for dynamic analysis of cranes, aerial work platforms, and other engineering equipment.

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

Fig. 6
Fig. 6 The endpoint dynamic deflection (m=0 kg) 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

Fig. 15
Fig. 15 Steady-state vibrating stage (m=50 kg) 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 2 nd -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 1 st -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.