Recursive Dynamic Algorithm of Open-Chain Multibody System

Open-chain multibody systems have been extensively studied because of their widespread application. Based on the structural characteristics of such a system, the relationship between its hinged bodies was transformed into recursive constraint relationships among the position, velocity, and acceleration of the bodies. The recursive relationships were used along with the Huston-Kane method to select the appropriate generalized coordinates and determine the partial velocity of each body and to develop an algorithm of the entire system. The algorithm was experimentally validated; it has concise steps and low susceptibility to error. Further, the algorithm can readily solve and analyze open-chain multibody systems.


Introduction
With the significant development of engineering and technology over the last 40 years, experts and researchers around the world have ceaselessly pursued creativity, which has included the development of diverse dynamics methods for solving multibody systems.Hooker and Wargulies [1,2] derived a general dynamics equation of a form that consists of  rigid bodies.Roberson and Wittenburg [3,4] described the structure of multibody systems based on concepts of graph theory and derived general dynamics equations of a rigid system.Ho [5] used the direct path method to derive dynamics equations of the motion of a flexible multibody spacecraft with a topological tree configuration.Kane presented what came to be known as Kane's method, which can be used to automatically eliminate the constraints force of a system and does not require the introduction of a differential scalar energy function [6][7][8].Haug et al. [9,10] developed the equations of motion of flexible dynamic systems, which are presented in this paper, using vibration and static correction elastic deformation modes.Furthermore, de jalón and Bayo [11] proposed an absolute Cartesian method.
The above works have facilitated the development of modern engineering and technology and provided a variety of effective methods for solving the problems of dynamic mechanical systems.Many institutions and researchers have used these methods to study systems.Among them are Houston et al. [12][13][14][15][16], who subsequently used Kane's method to develop a comprehensive approach and general dynamics formula for solving a multibody system.Tianjin University also developed the new multibody dynamics software and the flexible body dynamics software and successfully applied them to the analysis, modeling, and experimental study of dynamics [12,[17][18][19].Furthermore, Shanghai Jiaotong University conducted a series of studies on flexible multibody dynamics simulation software and theory [20][21][22], and the National Defense University studied a parachute recovery system using a multibody dynamics model [23][24][25][26].
However, these methods have drawbacks.For example, some of them do not have general dynamics equations that can be applied to any multibody system, and some require the user to have sufficient experience and skill to be able to use them to select the proper general velocity that would simplify the calculation process.Moreover, some of the methods involve several differential equations and constraint forces, which increase the difficulty of solving the dynamics equations of the system, just as others require extensive and laborious derivative operations.In the present study, vector analysis of a kinematic openchain multibody system was carried out based on the structural characteristics of the system.In combination with the Huston-Kane method, the vector analysis was used to select appropriate generalized coordinates, and the generalized velocity, rather than the generalized coordinates, was used as the independent variable for describing the motion of the system.The recursive dynamics equations of the system were also obtained.

Assumption.
The open-chain multibody system shown in Figure 1 is assumed to be composed of hinged rigid rods with their masses concentrated at their ends.There are  rigid bodies and  + 1 nodes in the system.Based on the low-number array method for describing multibody systems developed by Huston, one end of the system was selected as the starting point, and the nodes were numbered from that end.The nodes were numbered as 0, 1, 2, 3, . . ., , and the rods were numbered as 1, 2, 3, . . ., .The mass of each node is denoted by   , the length of each rod is denoted by   , and the angle between the -axis and each rod is denoted by   ().In the generalized coordinates, the first node position is defined by [ 0 (),  0 ()], and the angle between each body and the axis is defined by   ().

Position Analysis.
It is assumed that the position of the first node in the inertial frame at time  is [ 0 (),  0 ()] and that the angle between the th rod and the -axis is   ().The position of each node is thus given by

Velocity Analysis.
The velocity of each node is obtained as a derivative of the position of the node and is defined by θ  () =   ().Hence, the velocity of the first node is [ ẋ 0 () ẏ 0 ()]  , and those of subsequent nodes are ]  = 1, 2, . . ., . (2) 2.2.3.Partial Velocity Analysis.The partial velocity of each node relative to the generalized velocity can be obtained as follows: The derivative of the partial velocity can also be obtained as follows: ]  ≥ . (4)

Acceleration Analysis.
The derivative of the velocity of each node gives the acceleration.The acceleration of each node can thus be obtained using where q  is the th generalized velocity and q  is the derivative of the th generalized velocity.

Dynamic Analysis.
It is assumed that the active force acting on each node is   ,  = 0, 1, 2, . . ., , where  is the matrix formulation.

Generalized Active
Force.Using Kane's method, the generalized active force of the th generalized coordinates can be obtained as The generalized active force of the entire system is given by

Generalized Inertial
Force  * .Based on Kane's method, the generalized inertial force of the th generalized coordinates is The generalized inertial force of the entire system is where , q , and q are the matrix formulations of   , q  , and q  , respectively;   is the transpose matrix of ; and  = diag[ 0 ,  1 , . . .,   ].
By transposition, Equation ( 12) can be rewritten as where  =    and  =  −    u q .Equation ( 13) includes  + 2 equations and a total of  + 2 variables.The kinematic parameters of the system can thus be obtained.

Example 1.
Figure 2 shows a double pendulum.The system consists of two hinged pendulums, the movements of which are limited to the plane of the inertial reference system.Their lengths and masses are, respectively,  1 and  2 and  1 and  2 .As known, the system has two degrees of freedom and can therefore be described by two generalized coordinates.Two methods for solving the double pendulum are presented, namely, the Lagrange method and the method proposed in this paper.
In Lagrange's method, the absolute value of each of  1 and  2 is selected using the inertial reference system axis as the generalized coordinate of the system.The dynamics equations can thus be written as In the method proposed in this paper, the relative value of each of  1 and  2 is selected using the inertial reference system axis as the generalized coordinate of the system.Gravity is the only active force acting on the multibody system.
To simplify the calculations, it is assumed that  1 =  2 = 1,  1 =  2 = 1, and the initial angle of each pendulum relative to the inertial reference system axis is 18 ∘ .We solved this problem using Lagrange's method and the method of this paper, respectively.The dynamic responses are shown in Figure 3.
It can be seen that the dynamic responses obtained by the two methods are in good agreement, which confirms that the proposed method is valid for solving the double pendulum problem.

Example 2.
The system shown in Figure 4 is a linethrowing rocket.The coordinate system is the inertial coordinate system, and the launch point is defined as the origin of the coordinates.The -axis is in the flight direction of the rocket and is locally parallel to the "ground." The -axis is vertical to the "ground." The entire system is located on the XOY plane and its motion is planar.The rocket is simplified as a mass point and the rope is broken into  arbitrary segments in accordance with the finite segment method.The length of each rope segment is   , with the length and mass of the last segment being variable.The segments are numbered 1, 2, 3, . . ., , and labels from the rocket pull them to the ground.If the elongation and bending of the rope in the axial direction are not considered and assuming that the mass of each rope segment is mainly distributed at the end of the segment farther from the rocket and that the different segments are connected by a hinge, if the rope is pulled out and the length of the last segment changes till it fits the setting condition, a new rope segment  + 1 would be created.Let us consider the example of a rocket with the length of 1 m, diameter of 122 mm, and total weight of 20 kg and with gunpowder weight of 2.33 kg, rocket total impulse of 4770 N⋅s, and engine working time of 0.43 s.Further, the linear density of the rope is 0.043 Kg/m, and each rope section is  1 m long.The emission angle is 30 ∘ .The simulation and experimental results for this rocket were compared.The test setup was as shown in Figure 5.
The simulation and test results are summarized in Table 1.The range, maximum velocity, and flight time obtained by simulation and the test using a launch angle of 30 ∘ are given in the table.As can be seen, the simulation range was 7.2 m longer than that of the test, which represents an error of 1.13%.The maximum velocity obtained by simulation was 7.1 m/s higher than that of the test, indicating an error of 3.4%.The simulation flight time was 0.5 s shorter than that of the test.
The velocities of the rocket measured by radar and obtained by simulation for the launch angle of 30 ∘ are shown in Figure 6.From the figure, it can be seen that the maximum velocity measured during the test was 199.5 m/s, whereas the maximum simulation velocity was 206.6 m/s.When the engine was stopped, the test velocities after 1, 2, 3, and 4 s were 154.3, 99.3, 73.8, and 60.5 m/s, respectively, whereas the simulation velocities were 158.7, 98.9, 70.3, and 58.5 s, respectively.The errors are 2.8%, 0.4%, 5.0%, and 3.4%, respectively, and the average error is 2.9%.
The data in Table 1 and the velocity curves in Figure 6 show that the simulation results, including the velocity changes, are in good agreement with the test results.This verifies the precision of the dynamic model of the line-throwing rocket.

Figure 1 :
Figure 1: Simplified model of the open-chain multibody system.

Table 1 :
Comparison of the simulation and test results for a launch angle of 30 ∘ .