# Recursive Dynamic Algorithm of Open-Chain Multibody System

Authors: Ming Lu; Wenbin Gu; Jianqing Liu; Zhenxiong Wang; Zhisheng Jing; Guowei Qin; Shanchao Tu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

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.


---

## Body

## 1. 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 N 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–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–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–19]. Furthermore, Shanghai Jiaotong University conducted a series of studies on flexible multibody dynamics simulation software and theory [20–22], and the National Defense University studied a parachute recovery system using a multibody dynamics model [23–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 open-chain 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.


## 2. Dynamic Modeling

### 2.1. Assumption

The open-chain multibody system shown in Figure1 is assumed to be composed of hinged rigid rods with their masses concentrated at their ends. There are n rigid bodies and n+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,…,n, and the rods were numbered as 1,2,3,…,n. The mass of each node is denoted by mi, the length of each rod is denoted by li, and the angle between the y-axis and each rod is denoted by θi(t). In the generalized coordinates, the first node position is defined by [x0(t),y0(t)], and the angle between each body and the y-axis is defined by θi(t).Figure 1
Simplified model of the open-chain multibody system.


### 2.2. Kinematic Analysis

#### 2.2.1. Position Analysis

It is assumed that the position of the first node in the inertial frame at timet is [x0(t),y0(t)] and that the angle between the ith rod and the y-axis is θi(t). The position of each node is thus given by
(1)xityit=xi-1t-lisinθityi-1t-licos⁡θiti=1,2,…,n.


#### 2.2.2. Velocity Analysis

The velocity of each node is obtained as a derivative of the position of the node and is defined byθ˙i(t)=ωi(t). Hence, the velocity of the first node is x˙0(t)y˙0(t)′, and those of subsequent nodes are
(2)x˙i(t)y˙i(t)=x˙i-1t-licos⁡θitω1ty˙i-1t+lisinθitω1ti=1,2,…,n.


#### 2.2.3. Partial Velocity Analysis

The partial velocity of each node relative to the generalized velocity can be obtained as follows:(3)uij=00i<j,j≤n10j=n+101j=n+2-licos⁡θitlisinθiti≥j.
The derivative of the partial velocity can also be obtained as follows:
(4)u˙ij=00i<jlisinθitωitlicos⁡θitωiti≥j.


#### 2.2.4. Acceleration Analysis

The derivative of the velocity of each node gives the acceleration. The acceleration of each node can thus be obtained using(5)ai=∑j=0n+1uijq¨j+u˙ijq˙ji=0,1,2,…,n,
where q˙j is the jth generalized velocity and q¨j is the derivative of the jth generalized velocity.


### 2.3. Dynamic Analysis

It is assumed that the active force acting on each node isfzi, i=0,1,2,…,n, where fz is the matrix formulation.


#### 2.3.1. Generalized Active Force

Using Kane’s method, the generalized active force of thejth generalized coordinates can be obtained as
(6)Flj=∑i=1nfzi·uijj=1,2,…,n+2.
The generalized active force of the entire system is given by
(7)Fl=u′fz.


#### 2.3.2. Generalized Inertial ForceFl*

Based on Kane’s method, the generalized inertial force of thejth generalized coordinates is
(8)Flj*=∑i=1n-miai·uij=∑i=1n-miuijq¨j+u˙ijq˙j·uijaaaaaiaaaaj=1,2,…,n+2.
The generalized inertial force of the entire system is
(9)Fl*=-u′Muq¨-u′Mu˙q˙,
where u, q˙, and q¨ are the matrix formulations of uij, q˙j, and q¨j, respectively; u′ is the transpose matrix of u; and M=diag⁡m0,m1,…,mn.


### 2.4. Dynamics Equation

Based on Kane’s equation,(10)Fl+Fl*=0.
If (9) is substituted into (10), we have
(11)Fl-u′Muq¨-u′Mu˙q˙=0.
By transposition,
(12)u′Muq¨=Fl-u′Mu˙q˙.
Equation (12) can be rewritten as
(13)Aq¨=fl,
where A=u′Mu and fl=Fl-u′Mu˙q˙.Equation (13) includes n+2 equations and a total of n+2 variables. The kinematic parameters of the system can thus be obtained.


## 2.1. Assumption

The open-chain multibody system shown in Figure1 is assumed to be composed of hinged rigid rods with their masses concentrated at their ends. There are n rigid bodies and n+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,…,n, and the rods were numbered as 1,2,3,…,n. The mass of each node is denoted by mi, the length of each rod is denoted by li, and the angle between the y-axis and each rod is denoted by θi(t). In the generalized coordinates, the first node position is defined by [x0(t),y0(t)], and the angle between each body and the y-axis is defined by θi(t).Figure 1
Simplified model of the open-chain multibody system.


## 2.2. Kinematic Analysis

### 2.2.1. Position Analysis

It is assumed that the position of the first node in the inertial frame at timet is [x0(t),y0(t)] and that the angle between the ith rod and the y-axis is θi(t). The position of each node is thus given by
(1)xityit=xi-1t-lisinθityi-1t-licos⁡θiti=1,2,…,n.


### 2.2.2. Velocity Analysis

The velocity of each node is obtained as a derivative of the position of the node and is defined byθ˙i(t)=ωi(t). Hence, the velocity of the first node is x˙0(t)y˙0(t)′, and those of subsequent nodes are
(2)x˙i(t)y˙i(t)=x˙i-1t-licos⁡θitω1ty˙i-1t+lisinθitω1ti=1,2,…,n.


### 2.2.3. Partial Velocity Analysis

The partial velocity of each node relative to the generalized velocity can be obtained as follows:(3)uij=00i<j,j≤n10j=n+101j=n+2-licos⁡θitlisinθiti≥j.
The derivative of the partial velocity can also be obtained as follows:
(4)u˙ij=00i<jlisinθitωitlicos⁡θitωiti≥j.


### 2.2.4. Acceleration Analysis

The derivative of the velocity of each node gives the acceleration. The acceleration of each node can thus be obtained using(5)ai=∑j=0n+1uijq¨j+u˙ijq˙ji=0,1,2,…,n,
where q˙j is the jth generalized velocity and q¨j is the derivative of the jth generalized velocity.


## 2.2.1. Position Analysis

It is assumed that the position of the first node in the inertial frame at timet is [x0(t),y0(t)] and that the angle between the ith rod and the y-axis is θi(t). The position of each node is thus given by
(1)xityit=xi-1t-lisinθityi-1t-licos⁡θiti=1,2,…,n.


## 2.2.2. Velocity Analysis

The velocity of each node is obtained as a derivative of the position of the node and is defined byθ˙i(t)=ωi(t). Hence, the velocity of the first node is x˙0(t)y˙0(t)′, and those of subsequent nodes are
(2)x˙i(t)y˙i(t)=x˙i-1t-licos⁡θitω1ty˙i-1t+lisinθitω1ti=1,2,…,n.


## 2.2.3. Partial Velocity Analysis

The partial velocity of each node relative to the generalized velocity can be obtained as follows:(3)uij=00i<j,j≤n10j=n+101j=n+2-licos⁡θitlisinθiti≥j.
The derivative of the partial velocity can also be obtained as follows:
(4)u˙ij=00i<jlisinθitωitlicos⁡θitωiti≥j.


## 2.2.4. Acceleration Analysis

The derivative of the velocity of each node gives the acceleration. The acceleration of each node can thus be obtained using(5)ai=∑j=0n+1uijq¨j+u˙ijq˙ji=0,1,2,…,n,
where q˙j is the jth generalized velocity and q¨j is the derivative of the jth generalized velocity.


## 2.3. Dynamic Analysis

It is assumed that the active force acting on each node isfzi, i=0,1,2,…,n, where fz is the matrix formulation.


### 2.3.1. Generalized Active Force

Using Kane’s method, the generalized active force of thejth generalized coordinates can be obtained as
(6)Flj=∑i=1nfzi·uijj=1,2,…,n+2.
The generalized active force of the entire system is given by
(7)Fl=u′fz.


### 2.3.2. Generalized Inertial ForceFl*

Based on Kane’s method, the generalized inertial force of thejth generalized coordinates is
(8)Flj*=∑i=1n-miai·uij=∑i=1n-miuijq¨j+u˙ijq˙j·uijaaaaaiaaaaj=1,2,…,n+2.
The generalized inertial force of the entire system is
(9)Fl*=-u′Muq¨-u′Mu˙q˙,
where u, q˙, and q¨ are the matrix formulations of uij, q˙j, and q¨j, respectively; u′ is the transpose matrix of u; and M=diag⁡m0,m1,…,mn.


## 2.3.1. Generalized Active Force

Using Kane’s method, the generalized active force of thejth generalized coordinates can be obtained as
(6)Flj=∑i=1nfzi·uijj=1,2,…,n+2.
The generalized active force of the entire system is given by
(7)Fl=u′fz.


## 2.3.2. Generalized Inertial ForceFl*

Based on Kane’s method, the generalized inertial force of thejth generalized coordinates is
(8)Flj*=∑i=1n-miai·uij=∑i=1n-miuijq¨j+u˙ijq˙j·uijaaaaaiaaaaj=1,2,…,n+2.
The generalized inertial force of the entire system is
(9)Fl*=-u′Muq¨-u′Mu˙q˙,
where u, q˙, and q¨ are the matrix formulations of uij, q˙j, and q¨j, respectively; u′ is the transpose matrix of u; and M=diag⁡m0,m1,…,mn.


## 2.4. Dynamics Equation

Based on Kane’s equation,(10)Fl+Fl*=0.
If (9) is substituted into (10), we have
(11)Fl-u′Muq¨-u′Mu˙q˙=0.
By transposition,
(12)u′Muq¨=Fl-u′Mu˙q˙.
Equation (12) can be rewritten as
(13)Aq¨=fl,
where A=u′Mu and fl=Fl-u′Mu˙q˙.Equation (13) includes n+2 equations and a total of n+2 variables. The kinematic parameters of the system can thus be obtained.


## 3. Examples

### 3.1. Example 1

Figure2 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, l1 and l2 and m1 and m2. 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.Figure 2
Double pendulum.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
(14)m1+m2l1α¨1+m2l2cos⁡α2-α1α¨2-m2l2sinα2-α1α˙22+m1+m2gsinα1=0,m2l1cos⁡α2-α1α¨1+m2l2α¨2+m2l1sinα2-α1α˙12+m2gsinα2=0.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 thatm1=m2=1, l1=l2=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.Dynamic responses.
(a)
(b)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.


### 3.2. Example 2

The system shown in Figure4 is a line-throwing rocket. The coordinate system is the inertial coordinate system, and the launch point is defined as the origin of the coordinates. The x-axis is in the flight direction of the rocket and is locally parallel to the “ground.” The y-axis is vertical to the “ground.”Figure 4
Simplified model of the line-throwing rocket.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 inton arbitrary segments in accordance with the finite segment method. The length of each rope segment is li, with the length and mass of the last segment being variable. The segments are numbered 1,2,3,…,n, 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 n+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 Figure5.Test setup.
(a)
(b)The simulation and test results are summarized in Table1. 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.Table 1
Comparison of the simulation and test results for a launch angle of 30°.
Parameters
Range (m)
Maximum velocity (m/s)
Flight time (s)
Simulation
637.1
206.6
8.1
Test
629.9
199.5
8.6
Absolute error
7.2
7.1
0.5The velocities of the rocket measured by radar and obtained by simulation for the launch angle of 30° are shown in Figure6. 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%.Figure 6
Curves of simulation and test velocities for a launch angle of 30°.The data in Table1 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.


## 3.1. Example 1

Figure2 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, l1 and l2 and m1 and m2. 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.Figure 2
Double pendulum.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
(14)m1+m2l1α¨1+m2l2cos⁡α2-α1α¨2-m2l2sinα2-α1α˙22+m1+m2gsinα1=0,m2l1cos⁡α2-α1α¨1+m2l2α¨2+m2l1sinα2-α1α˙12+m2gsinα2=0.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 thatm1=m2=1, l1=l2=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.Dynamic responses.
(a)
(b)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.


## 3.2. Example 2

The system shown in Figure4 is a line-throwing rocket. The coordinate system is the inertial coordinate system, and the launch point is defined as the origin of the coordinates. The x-axis is in the flight direction of the rocket and is locally parallel to the “ground.” The y-axis is vertical to the “ground.”Figure 4
Simplified model of the line-throwing rocket.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 inton arbitrary segments in accordance with the finite segment method. The length of each rope segment is li, with the length and mass of the last segment being variable. The segments are numbered 1,2,3,…,n, 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 n+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 Figure5.Test setup.
(a)
(b)The simulation and test results are summarized in Table1. 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.Table 1
Comparison of the simulation and test results for a launch angle of 30°.
Parameters
Range (m)
Maximum velocity (m/s)
Flight time (s)
Simulation
637.1
206.6
8.1
Test
629.9
199.5
8.6
Absolute error
7.2
7.1
0.5The velocities of the rocket measured by radar and obtained by simulation for the launch angle of 30° are shown in Figure6. 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%.Figure 6
Curves of simulation and test velocities for a launch angle of 30°.The data in Table1 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.


## 4. Conclusions

By analysis, the recursive relationships among the location, velocity, and acceleration of the different parts of an open-chain multibody system were determined and were used in combination with Houston-Kane’s method to select the appropriate generalized velocity and develop a dynamic recursive algorithm of the multibody system. The proposed method does not require the solution of the constrained reaction force in the modeling process, and it involves fewer variables and equations, which improves the computational efficiency of the dynamic model.The developed dynamic recursive algorithm was used to simulate and analyze two examples. The dynamic responses of the double pendulum obtained by simulation were found to be in good agreement with the results of test for a line-throwing rocket. The validity of the dynamic recursive algorithm for solving an open-chain multibody system was thus verified.



---

*Source: 457682-2014-12-01.xml*