Sage Journals: Discover world-class research

Abstract

The dynamics for spatial manipulator arms consisting of n flexible links and n flexible joints is presented. All the transversal, longitudinal, and torsional deformation of flexible links are considered. Within the total longitudinal deformation, the nonlinear coupling term, also known as the longitudinal shortening caused by transversal deformation, also is considered here. Each flexible joint is modeled as a linearly elastic torsional spring, and the mass of joint is considered. Lagrange's equations are adopted to derive the governing equations of motion of the system. The algorithmic procedure is based on recursive formulation using 4 × 4 homogenous transformation matrices where all the kinematical expressions as well as the final equations of motion are suited for computation. A corresponding general-purpose C++ software package for dynamic simulation is developed. Several examples are simulated to illustrate the performance of the algorithm.

1. Introduction

Flexible multibody dynamics has become a key methodology for various engineering fields, such as robotic manipulator arms, large radar antennas, solar panels, transportation vehicles, and manufacturing equipment as well as flexible ligament in human musculoskeletal system [1], and so forth.

In particular, the flexible robotic manipulators are typically precise and complex operated at high speed. They are designed to be light with low inertia in order to achieve cost-reduction, energy-saving, and high-performance. Therefore, the dynamic analysis of flexible manipulator arms is complicated due to the coupling between the large rigid body motion and deformation.

A large number of the literature related to the so-called dynamic stiffening effects have been published, which is first proposed by Kane et al. [2] due to link's high-order coupling flexibility. It was observed that an industrial link under high-speed rotational motion would exhibit instability problems when its angular velocity exceeds a certain limit. The incorrect simulation solutions are attributed to neglecting the high-order coupling deformation terms in the dynamic equations. Wu and Haug [3] modeled a flexible multibody system by means of substructure synthesis formulation to account for structural geometric non-linear effects. Xi and Fenton [4] studied a manipulator consisting of one flexible link and one flexible joint based on the assumption that the link is constrained to move only in a horizontal plane, whereas the gravity and some coupling terms between the equations of motions have been dropped. It has been demonstrated by Wallrapp and Schwertassek [5], that the so-called “geometric stiffening” problem can be solved by keeping second order terms in the expression of the deformations of the material, and the amplitude of the flexible motions must remain small in those formalisms which are generally based on kinematic restrictions as regards the flexibility.

Low [6] presented vibration analysis of a rotating beam carrying a tip mass at its end by using Hamilton's principle and the associated boundary conditions. Yoo et al. [7] used a non-Cartesian variable along with two Cartesian variables to describe the elastic deformation and investigate the dynamic stiffness effect. Ryu et al. [8] put forward a criterion on inclusion of stiffening effects such that it clarifies the limit of the validity of the linear modeling method. El-Absy and Shabana [9] introduced the effect of longitudinal deformation due to bending and studied the influence of geometric stiffness on instability problem of nonlinear elastic model. Al-Bedoor and Hamdan [10], based on the condition of inextensibility to relate the axial and transverse deformations of the material point, studied a rotating flexible arm deformation undergoing large planar motion. Hong et al. [11, 12] established higher-order rigid-flexible coupling dynamic model of rotating flexible beam undergoing large overall motion. When the deformation of flexible beam or the deformation rate is larger, the first-order or zero-order approximation coupling model will appear with divergent results, nevertheless the complete high-order coupling model's results are still converged. Schwertassek et al. [13] presented the fundamental shape functions choice in floating frame of reference formulation, by separating the flexible body motion into a reference motion and deformation.

Book [14] used the 4 × 4 homogenous transformation matrices and assumed modes to describe the kinematics of the rigid-joint and flexible-link robots, whereas the dynamic model cannot deal with the torsional deformation of links. A recursive formulation for the spatial kinematic and dynamic analysis of open chain mechanical systems containing interconnected deformable bodies is given by Changizi and Shabana [15], and Kim and Haug [16]. Jain and Rodriguez [17] developed new spatially recursive dynamics algorithms for flexible multibody systems by using spatial operators, which is based on Newton-Euler factorization and innovations factorization of the system mass matrix. Hwang [18] developed a recursive formulation for the flexible dynamic manufacturing analysis of open-loop robotic systems with the generalized Newton-Euler equations. Zhang and Zhou [19, 20] did a further work based on Book's work in [14]. Both the bending and torsional flexibility of links were taken into account. Dynamic simulation of a spatial flexible manipulator arm was given as an example to validate the algorithm. However, the dynamic stiffening effect was not considered yet.

When the system's operating speed becomes high, neglecting the flexibility of joints is quite devastating, which will usually give rise to the errors of position precision. A typical dynamic model for taking into account the joint flexibility was presented by Spong [21]. In his work, the joint flexibility is modeled as a torsional spring with more emphasis on simplifying the equations of motion for control purposes. Other dynamic analyses of robots with joint flexibilities can be seen in the work of Wasfy and Noor [22], Dwivedy and Eberhard [23], and Na and Kim [24].

As mentioned above, there are a lot of research work in the dynamic modeling and simulation. But it is still very difficult for us to deal with the dynamics of the complex multibody systems, such as the spatial flexible-link and flexible-joint robotic manipulators with consideration of the axial, bending, and torsional deformation for links, and the flexibility and mass effects for joints, and the so-called “dynamic stiffening” effects. In this paper, we will present the dynamic modeling methodology to include all such terms. In the following section, the kinematics of the system are presented, in which coordinate frames are established, and 4 × 4 homogeneous transformation matrices are used to describe the kinematics of flexible links and flexible joints. The approach of assumed modes is employed to describe the deformation of the flexible links. Section 3 firstly deals with the description of the kinetic and elastic potential energy as well as the gravitational potential energy of system, and then focuses on the derivation of the recursive rigid-flexible coupling dynamic equations of the system. In the modeling, the high-order coupling terms related to the non-linear geometry deformation are retained and the recursive strategy for kinematics is adopted. In Section 4, several examples are simulated to illustrate the performance of the algorithm proposed in the paper. Finally, Section 5 summarizes the results and draws conclusions from them.

2. Kinematics of Flexible Robots

As a point of departure, the system considered here is an assembly of n flexible links connected by n rotary joints, as shown in Figure 1.

Figure 1:

Flexible robotic arms.

2.1. Simplified Model of Flexible Joint

Figure 2 shows the flexible joint model. We can model the flexibility of joint i as a linear torsional spring with stiffness Kt_i · J_ri is the moment of inertia of rotor i about its spinning axis. And τ_i is the torque exerted at joint i. For simplicity, we neglect friction or damping in the flexible joint. Let q_1i be the theoretical rotational angle of link i, q_2i be the real rotational angle of link i, ∊_i be the torsional angle of joint i, φ_i be the angular displacement of rotor i, and ϕ_i be the gear ratio, respectively. The relationships among them are as follows:

∊_{i} = q_{2 i} - q_{1 i},

(1)

φ_{i} = ϕ_{i} q_{1 i} .

(2)

Figure 2:

Flexible joint model.

2.2. Simplified Model of Flexible Link

Assume that the links are slender beams. Analysis here is based on the Euler-Bernoulli beam theory in the elastic small displacements field.

2.3. Coordinate Systems and Transformation Matrices

To express the transformation between different coordinate systems clearly, we establish four coordinate systems for link i. Fix the coordinate system (X_bY_bZ_b)_i at the proximal end of link i (oriented so that the X_b coincides with the neutral axis of link i in undeformed shape). This will be referred to as the base reference frame of link i. Fix the coordinate system (X_dY_dZ_d)_i at the distal end of link i. This is the distal frame of link i. When link i is in its undeformed state, the distal frame can be located by a pure translation of the base reference frame (X_bY_bZ_b)_i along the length L_i of link i. Let (H_xH_yH_z)_i′ and (H_xH_yH_z)_i be two Denavit-Hartenberg (D-H) frames fixed at the proximal end (at joint i) and the distal end (at joint i + 1) of link i, respectively. When joint i is motionless, (H_xH_yH_z)_i′ is coincident with (H_xH_yH_z)_{i – 1}, and matrix $H H_{i}^{i - 1}$ , that is, the transformation matrix between them, is the function of q_2i. Matrix H_i, the 4 × 4 homogeneous transformation matrix between frames (X_dY_dZ_d)_i and (H_xH_yH_z)_i, is a constant matrix. The transformation matrix Hb_i of (H_xH_yH_z)_i′ and (X_bY_bZ_b)_i is also a constant matrix. Define the joint-transformation matrix A_i of joint i to be the transformation matrix from (X_dY_dZ_d)_{i – 1} to (X_bY_bZ_b)_i. Then,

A_{i} = d H_{i - 1} H H_{i}^{i - 1} H b_{i} .

(3)

Obviously, A_i is a function of q_2i. Define E_i to be the link-transformation matrix of link i from (X_bY_bZ_b)_i to (X_dY_dZ_d)_i. According to the assumption of small deformation of the links, the small deformable angles can be added vectorially. E_i can be written as

E_{i} = H_{i} + \sum_{j = 1}^{N_{i}} δ_{i j} M_{i j},

(4)

in which

H_{i} = [\begin{bmatrix} 1 & 0 & 0 & 0 \\ L_{i} & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}],

(5)

M_{i j} = [\begin{bmatrix} 0 & 0 & 0 & 0 \\ x_{i j} & 1 & - θ_{z i j} & θ_{y i j} \\ y_{i j} & θ_{z i j} & 0 & - θ_{x i j} \\ z_{i j} & - θ_{y i j} & θ_{x i j} & 0 \end{bmatrix}] .

(6)

Here, x_ij, y_ij, and z_ij are the x_bi, y_bi, and z_bi components of the elastic linear displacement mode j of link i at the origin of the coordinate (X_dY_dZ_d)_i, respectively. θ_xij, θ_yij, and θ_zij are the x_bi, y_bi, and z_bi rotation components of the elastic angular displacement mode j of link i at the origin of the coordinate (X_dY_dZ_d)_i, respectively. δ_ij is the time-varying amplitude of mode j of link i, and N_i is the number of modes used to describe the deformation of link i.

Define ${}^{0}W_{i}$ or W_i to be the 4 × 4 homogeneous transformation matrix from the base coordinate frame (X_bY_bZ_b)₀ to (X_bY_bZ_b)_i. Then, we have

{}^{0}{W_{i}} = W_{i} = W_{i - 1} E_{i - 1} A_{i} = {\hat{W}}_{i - 1} A_{i},

(7)

where ${\hat{W}}_{i - 1}$ is the 4 × 4 homogeneous transformation matrix from the base coordinate frame (X_b Y_b Z_b)₀ to the distal coordinate system (X_d Y_d Z_d)_{i – 1} of link i – 1.

2.4. Velocity of a Point of Link i

Let ${}^{i}{h_{i}} (η)$ be the homogeneous coordinates in the system (X_bY_bZ_b)_i of a point of the deformed link i at position η with the link under an undeformed condition from the origin of (X_bY_bZ_b)_i. Then, ⁱh_i (η) can be approximated as

^{i} h_{i} (η) = {[\begin{bmatrix} 1 & η & 0 & 0 \end{bmatrix}]}^{T} + \sum_{j = 1}^{N_{i}} δ_{i j} {[\begin{bmatrix} 0 & x_{i j} (η) & y_{i j} (η) & z_{i j} (η) \end{bmatrix}]}^{T} - \frac{1}{2} \sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} δ_{i j} δ_{i k} {[\begin{matrix} 0 & {\underset{}{x}}_{i j k} (η) & 0 & 0 \end{matrix}]}^{T} .

(8)

Here, ${\underset{}{x}}_{i j k}$ is the nonlinear strain coupling term, also known as the axial shortening due to the bending deformations of the link. When the flexible links are undergoing a high speed, this term will bring the so-called dynamic stiffening effect which will have a great influence on the dynamic behavior of flexible arms. In the dynamic modeling presented here, the high-order terms related to the non-linear coupling term are retained, which are ignored in the first-order and zero-order approximation coupling modeling.

In terms of the fixed inertial coordinates of the base (X_bY_bZ_b)₀, the position ⁰h_i or h_i of the point is given as

h_{i} = W_{i}^{} {}^{i}{h_{i}} .

(9)

Taking the time derivative of the position h_i, we have the velocity of the point as

\frac{d h_{i}}{d t} = {\dot{h}}_{i} = {\dot{W}}_{i} {}^{i}{h_{i}} + W_{i} {}^{i}{\dot{h}}_{i} .

(10)

To accelerate the computation of the matrices ${\dot{W}}_{i}$ or ${\overset{‥}{W}}_{i}$ , we use the recursive kinematics method. By differentiating (7), one obtains

\begin{matrix} {\dot{W}}_{i} = {\dot{\hat{W}}}_{i - 1} A_{i} + {\hat{W}}_{i - 1} {\dot{A}}_{i}, \\ {\overset{‥}{W}}_{i} = {\overset{‥}{\hat{W}}}_{i - 1} A_{i} + 2 {\dot{\hat{W}}}_{i - 1} {\dot{A}}_{i} + {\hat{W}}_{i - 1} {\overset{‥}{A}}_{i}, \end{matrix}

(11)

where

\begin{matrix} {\dot{A}}_{i} = U_{i} {\dot{q}}_{2 i}, \\ {\overset{‥}{A}}_{i} = U_{2 i} {\dot{q}}_{2 i}^{2} + U_{i} {\overset{‥}{q}}_{2 i} . \end{matrix}

(12)

Here, U_i ≜ ∂ A_i/∂ q_2i, U_2i ≜ ∂ ²A_i/∂ q_2i², and q_2i is the joint variable of joint i. Thus, ${\dot{W}}_{i}$ and ${\overset{‥}{W}}_{i}$ can be computed recursively from ${\hat{W}}_{i - 1}$ and its derivatives. Here, one additionally needs ${\hat{W}}_{i - 1}$ and its derivatives. These can be computed recursively from W_{i – 1} and its derivatives as follows:

\begin{matrix} {\hat{W}}_{i} = W_{i} E_{i}, \\ {\dot{\hat{W}}}_{i} = {\dot{W}}_{i} E_{i} + W_{i} {\dot{E}}_{i}, \\ {\overset{‥}{\hat{W}}}_{i} = {\overset{‥}{W}}_{i} E_{i} + 2 {\dot{W}}_{i} {\dot{E}}_{i} + W_{i} {\overset{‥}{E}}_{i}, \end{matrix}

(13)

where

\begin{matrix} {\dot{E}}_{i} = \sum_{j = 1}^{N_{i}} {\dot{δ}}_{i j} M_{i j}, \\ {\overset{‥}{E}}_{i} = \sum_{j = 1}^{N_{i}} {\overset{‥}{δ}}_{i j} M_{i j} . \end{matrix}

(14)

3. Dynamics of Flexible Robots

To use Lagrange's equations, we need the kinetic and potential energy of the system.

3.1. The System Kinetic Energy

The system kinetic energy K contains two parts: the links kinetic energy $K_{b}$ and the joints kinetic energy $K_{r}$ as follows:

K = K_{}^{b} + K_{}^{r} .

(15)

Assume that the links are slender beams, so the rotary inertia and shear effects can be neglected. Therefore, the present analysis is based on the Euler-Bernoulli beam theory. Also assume that the links can undergo a large overall rigid motion, however the elastic displacements are small. The kinetic energy of the ith link is

K_{i}^{b} = \frac{1}{2} \int_{0}^{L_{i}} μ (η) Tr {{\dot{h}}_{i} {\dot{h}}_{i}^{T}} d η + \frac{1}{2} {\int_{0}^{L_{i}} J_{x i} (η) (\frac{\partial θ_{x i} (η, t)}{\partial t})}^{2} d η,

(16)

where Tr {·} is the trace operator; μ(η) and J_xi are the mass per unit length and the polar moment of inertia per unit length of the link about the neutral axis x, respectively. For slender beams with uniform cross section area along x axis, μ(η) = μ_i. The first term in (16) is the kinetic energy of link i accounting for the rigid-body motion and the lateral and longitudinal deformation due to flexibility, whereas the second term is the kinetic energy accounting for the torsional deformation. Note that the torsional angle of link i, θ_xi, is expressed as

θ_{x i} = \sum_{j = 1}^{N_{i}} δ_{i j} θ_{x i j},

(17)

where δ_ij and θ_xij are mentioned in (4) and (6), respectively. δ_ij is the function of time, whereas θ_xij is the function of position η.

Substituting (10) and (17) into (16), expanding it, and then summing over all n links, one finds the links' kinetic energy to be

K_{}^{b} = \sum_{i = 1}^{n} K_{i}^{b} = \sum_{i = 1}^{n} Tr {{\dot{W}}_{i} B_{3 i} W_{i}^{T} + 2 {\dot{W}}_{i} B_{2 i} W_{i}^{T} + W_{i} B_{1 i} W_{i}^{T}} + \sum_{i = 1}^{n} \sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} {\dot{δ}}_{i j} {\dot{δ}}_{i k} T_{i j k},

(18)

where

\begin{matrix} T_{i j k} = \frac{1}{2} \int_{0}^{L_{i}} J_{x i} θ_{x i j} θ_{x i k} d η, \\ B_{1 i} = \frac{1}{2} \int_{0}^{L_{i}} μ_{i} {}^{i}{\dot{h}}_{i} {}^{i}{\dot{h}}_{i}^{T} d η, \\ B_{2 i} = \frac{1}{2} \int_{0}^{L_{i}} μ_{i} {}^{i}{h_{i}} {}^{i}{\dot{h}}_{i}^{T} d η, \\ B_{3 i} = \frac{1}{2} \int_{0}^{L_{i}} μ_{i} {}^{i}{h_{i}} {}^{i}{h_{i}^{T}} d η . \end{matrix}

(19)

With the consideration of ${}^{i}{h_{i}}$ expressed in (8) and its derivative, consider

{}^{i}{\dot{h}}_{i} (η) = \sum_{j = 1}^{N_{i}} {\dot{δ}}_{i j} {[\begin{bmatrix} 0 & x_{i j} (η) & y_{i j} (η) & z_{i j} (η) \end{bmatrix}]}^{T} - \sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} {\dot{δ}}_{i j} δ_{i k} {[\begin{bmatrix} 0 & {\underset{}{x}}_{i j k} (η) & 0 & 0 \end{bmatrix}]}^{T},

(20)

the matrices B_1i, B_2i, and B_3i can be written as

\begin{matrix} B_{1 i} = \sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} {\dot{δ}}_{i j} {\dot{δ}}_{i k} D_{i j k} + \underset{}{\sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} \sum_{l = 1}^{N_{i}} {\dot{δ}}_{i j} {\dot{δ}}_{i k} δ_{i l} (- F_{i j k l} - F_{i j k l}^{T})} + \underset{}{\sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} \sum_{l = 1}^{N_{i}} \sum_{s = 1}^{N_{i}} {\dot{δ}}_{i j} δ_{i k} {\dot{δ}}_{i l} δ_{i s} E_{i j k l s}}, \end{matrix}

(21)

B_{2 i} = \sum_{j = 1}^{N_{i}} {\dot{δ}}_{i j} D_{i j} + \sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} {\dot{δ}}_{i j} δ_{i k} D_{i j k}^{T} + \underset{}{\sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} {\dot{δ}}_{i j} δ_{i k} (- E_{i j k})} + \underset{}{\sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} \sum_{l = 1}^{N_{i}} {\dot{δ}}_{i j} δ_{i k} δ_{i l} (- \frac{1}{2} F_{i j k l}^{T} - F_{i l j k})} + \underset{}{\sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} \sum_{l = 1}^{N_{i}} \sum_{s = 1}^{N_{i}} {\dot{δ}}_{i j} δ_{i k} δ_{i l} δ_{i s} (\frac{1}{2} E_{i j l k s}^{T})},

(22)

B_{3 i} = D_{i} + \sum_{j = 1}^{N_{i}} δ_{i j} (D_{i j} + D_{i j}^{T}) + \sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} δ_{i j} δ_{i k} D_{i j k} + \underset{}{\sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} δ_{i j} δ_{i k} (- \frac{1}{2} E_{i j k} - \frac{1}{2} E_{i j k}^{T})} + \underset{}{\sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} \sum_{l = 1}^{N_{i}} δ_{i j} δ_{i k} δ_{i l} (- \frac{1}{2} F_{i j k l} - \frac{1}{2} F_{i j k l}^{T})} + \underset{}{\sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} \sum_{l = 1}^{N_{i}} \sum_{s = 1}^{N_{i}} δ_{i j} δ_{i k} δ_{i l} δ_{i s} (\frac{1}{4} E_{i j l k s})},

(23)

where

\begin{matrix} D_{i} = \frac{1}{2} \int_{0}^{L_{i}} μ_{i} L_{i} L_{i}^{T} d η, \\ D_{i j} = \frac{1}{2} \int_{0}^{L_{i}} μ_{i} L_{i} U_{i j}^{T} d η, \\ D_{i j k} = \frac{1}{2} \int_{0}^{L_{i}} μ_{i} U_{i j} U_{i k}^{T} d η, \\ E_{i j k} = \frac{1}{2} \int_{0}^{L_{i}} μ_{i} L_{i} S_{i j k}^{T} d η, \\ F_{i j k l} = \frac{1}{2} \int_{0}^{L_{i}} μ_{i} U_{i j} S_{i k l}^{T} d η, \\ E_{i j k l s} = \frac{1}{2} \int_{0}^{L_{i}} μ_{i} S_{i j k} S_{i l s}^{T} d η, \end{matrix}

(24)

in which

\begin{matrix} L_{i} = {[\begin{bmatrix} 1 & η & 0 & 0 \end{bmatrix}]}^{T}, \\ U_{i j} = {[\begin{bmatrix} 0 & x_{i j} (η) & y_{i j} (η) & z_{i j} (η) \end{bmatrix}]}^{T}, \\ S_{i j k} = {[\begin{bmatrix} 0 & {\underset{}{x}}_{i j k} (η) & 0 & 0 \end{bmatrix}]}^{T} . \end{matrix}

(25)

In (21)–(23), the terms with underline remained in high-order approximation coupling model.

It should be noted that the link shape mentioned above is restricted to be the slender beam type. In fact, the link shape can further be extended to the other cases.

Case 1. Link i is the rigid-body with irregular shape. In this case, N_i = 0, B_1i = B_2i = 0, and B_3i = D_i. The term D_i is the equivalent of the inertia moment of the rigid-link. It is actually the pseu-matrix of the inertia moment in [25] with its more complex form compared with that of (23). Thus, for rigid-link, the link shape can be arbitrary, and one should only input the corresponding matrix D_i.

Case 2. Link i consists of a flexible beam with the concentrated mass m_i at its proximal end. To account for the contribution of the concentrated mass to the kinetic energy of link i, the extra term $D_{i}^{'}$ should be added to the matrix D_i, where

D_{i}^{'} = \frac{1}{2} m_{i} {[\begin{bmatrix} 1 & 0 & 0 & 0 \end{bmatrix}]}^{T} [\begin{bmatrix} 1 & 0 & 0 & 0 \end{bmatrix}] .

(26)

Should the concentrated mass m_i locate at the position (x₀, y₀, z₀) near the proximal end, the extra matrix $D_{i}^{'}$ is modified to be

D_{i}^{'} = \frac{1}{2} m_{i} {[\begin{bmatrix} 1 & x_{0} & y_{0} & z_{0} \end{bmatrix}]}^{T} [\begin{bmatrix} 1 & x_{0} & y_{0} & z_{0} \end{bmatrix}] .

(27)

Case 3. Link i consists of a flexible beam with the concentrated mass m_i at its distal end. Considering the contribution of the concentrated mass to the kinetic energy of link i, the extra terms $D_{i}^{″}$ , $D_{i j}^{″}$ , $D_{i j k}^{″}$ , $E_{i j k}^{″}$ , $F_{i j k l}^{″}$ , and $E_{i j k l s}^{″}$ should be added to the corresponding matrix D_i, D_ij, D_ijk, E_ijk, F_ijkl, and E_ijkls, respectively. Here,

\begin{matrix} D_{i}^{′′} = \frac{1}{2} m_{i} {[\begin{bmatrix} 0 & L_{i} & 0 & 0 \end{bmatrix}]}^{T} [\begin{bmatrix} 0 & L_{i} & 0 & 0 \end{bmatrix}], \\ D_{i j}^{′′} = \frac{1}{2} m_{i} {[\begin{bmatrix} 0 & L_{i} & 0 & 0 \end{bmatrix}]}^{T} [\begin{bmatrix} 0 & x_{i j} (L_{i}) & y_{i j} (L_{i}) & z_{i j} (L_{i}) \end{bmatrix}], \\ D_{i j k}^{′′} = \frac{1}{2} m_{i} {[\begin{bmatrix} 0 & x_{i j} (L_{i}) & y_{i j} (L_{i}) & z_{i j} (L_{i}) \end{bmatrix}]}^{T} \times [\begin{bmatrix} 0 & x_{i k} (L_{i}) & y_{i k} (L_{i}) & z_{i k} (L_{i}) \end{bmatrix}], \\ E_{i j k}^{′′} = \frac{1}{2} m_{i} {[\begin{bmatrix} 0 & L_{i} & 0 & 0 \end{bmatrix}]}^{T} [\begin{bmatrix} 0 & {\underset{}{x}}_{i j k} (L_{i}) & 0 & 0 \end{bmatrix}], \\ F_{i j k l}^{′′} = \frac{1}{2} m_{i} {[\begin{bmatrix} 0 & x_{i j} (L_{i}) & y_{i j} (L_{i}) & z_{i j} (L_{i}) \end{bmatrix}]}^{T} \times [\begin{bmatrix} 0 & {\underset{}{x}}_{i k l} (L_{i}) & 0 & 0 \end{bmatrix}], \\ E_{i j k l s}^{′′} = \frac{1}{2} m_{i} {[\begin{bmatrix} 0 & {\underset{}{x}}_{i j k} (L_{i}) & 0 & 0 \end{bmatrix}]}^{T} [\begin{bmatrix} 0 & {\underset{}{x}}_{i l s} (L_{i}) & 0 & 0 \end{bmatrix}] . \end{matrix}

(28)

For the calculation of the kinetic energy of the joint i, we can lump its mass to link i – 1 in accordance with the assumptions made in [21] for simplicity. Thus, the kinetic energy needs to be included only in the part accounting for spinning kinetic energy of the rotor of joint i as follows:

K_{i}^{r} = \frac{1}{2} J_{r i} {\dot{φ}}_{i}^{2} .

(29)

Considering the relationship in (2). We can obtain the kinetic energy of the joints as follows:

K_{}^{r} = \sum_{i = 1}^{n} K_{i}^{r} = \sum_{i = 1}^{n} (\frac{1}{2} J_{r i} ϕ_{i}^{2} {\dot{q}}_{1 i}^{2}) .

(30)

3.2. The System Potential Energy

The potential energy of the flexible manipulators V mainly includes the elastic potential energy of flexible joints V^r, the elastic potential energy of the flexible links $V_{b}$ , and their gravitational potential energy V_g. Therefore, the whole system potential energy is

V = V_{}^{r} + V_{}^{b} + V_{g} .

(31)

Firstly, the elastic potential energy of joint i is

V_{i}^{r} = \frac{1}{2} K t_{i} {(q_{2 i} - q_{1 i})}^{2} .

(32)

Secondly, the elastic potential of the flexible links considered here have three parts: one by bending about the transverse y_i and z_i axes, one by compressing about the longitudinal x_i axis, and one by twisting about the longitudinal x_i axis. Along an incremental length dη, the elastic potential energy is

d V_{i}^{b} = \frac{1}{2} {E_{i} [I_{z i} {(\frac{\partial θ_{z i}}{\partial η})}^{2} + I_{y i} {(\frac{\partial θ_{y i}}{\partial η})}^{2} + A_{x i} {(\frac{\partial x_{i}}{\partial η})}^{2}] + G_{i} I_{x i} {(\frac{\partial θ_{x i}}{\partial η})}^{2}} d η .

(33)

Here, θ_xi, θ_yi, and θ_zi are the ith link's rotations of the neutral axis at the point η in the x_i, y_i, and z_i directions, respectively; E_i is Young's modulus of material of link i; G_i is the shear modulus of material of link i; A_xi is the ith link's area of cross section about the x_i axis; I_yi and I_zi are the area moment of inertia of ith link's cross section about the y_i and z_i axes, respectively; and I_xi is the polar area moment of inertia of the ith link's cross section about the neutral axis. Similar to the torsional angle θ_xi of (17), θ_yi and θ_zi can be expressed as

\begin{matrix} θ_{y i} = \sum_{k = 1}^{N_{i}} δ_{i k} θ_{y i k}, \\ θ_{z i} = \sum_{k = 1}^{N_{i}} δ_{i k} θ_{z i k}, \end{matrix}

(34)

where θ_yik and θ_zik are mentioned in (6). By integrating dV_i^b of (33) over the link, and summing over all n links, one can obtain $V_{b}$ as

V_{}^{b} = \sum_{i = 1}^{n} (\int_{0}^{L_{i}} d V_{i}^{b}) = \frac{1}{2} \sum_{i = 1}^{n} \sum_{j = 1}^{N_{i}} \sum_{k = 1}^{N_{i}} δ_{i j} δ_{i k} K_{i j k},

(35)

where

\begin{matrix} K_{i j k} = K_{x i j k} + K_{y i j k} + K_{z i j k} + K_{t i j k}, \\ K_{x i j k} = \int_{0}^{L_{i}} E_{i} A_{x i} (η) \frac{\partial x_{i j}}{\partial η} \frac{\partial x_{i k}}{\partial η} d η, \\ K_{y i j k} = \int_{0}^{L_{i}} E_{i} I_{y i} (η) \frac{\partial θ_{y i j}}{\partial η} \frac{\partial θ_{y i k}}{\partial η} d η, \\ K_{z i j k} = \int_{0}^{L_{i}} E_{i} I_{z i} (η) \frac{\partial θ_{z i j}}{\partial η} \frac{\partial θ_{z i k}}{\partial η} d η, \\ K_{t i j k} = \int_{0}^{L_{i}} G_{i} I_{x i} (η) \frac{\partial θ_{x i j}}{\partial η} \frac{\partial θ_{x i k}}{\partial η} d η . \end{matrix}

(36)

Finally, by using the same procedure given by Book [14], the gravity potential of the system is

V_{g} = - g^{T} \sum_{i = 1}^{n} W_{i} r_{i},

(37)

in which g is the gravity vector with respect to the inertial base, and it has the following form:

\begin{matrix} g^{T} = [\begin{bmatrix} 0 & g_{x} & g_{y} & g_{z} \end{bmatrix}]; \\ r_{i} = M_{i} r_{r i} + \sum_{k = 1}^{N_{i}} δ_{i k} ∊_{i k} - \frac{1}{2} \sum_{k =}^{1} N_{i} \sum_{l = 1}^{N_{i}} δ_{i k} δ_{i l} κ_{i k l} . \end{matrix}

(38)

Here, M_i is the total mass of link i and r_ri is the homogenous coordinates of the gravity center of link i (undeformed) in the frame (X_bY_bZ_b)_i,

\begin{matrix} r_{r i} = {[\begin{bmatrix} 1 & r_{x i} & 0 & 0 \end{bmatrix}]}^{T}, \\ ∊_{i k} = \int_{0}^{L_{i}} μ_{i} {[\begin{bmatrix} 0 & x_{i k} & y_{i k} & z_{i k} \end{bmatrix}]}^{T} d η, \\ κ_{i k l} = \int_{0}^{L_{i}} μ_{i} {[\begin{bmatrix} 0 & {\underset{}{x}}_{i k l} & 0 & 0 \end{bmatrix}]}^{T} d η . \end{matrix}

(39)

Note that ∊_ik and κ_ikl can be found in the top row of D_ik and E_ijk, respectively.

3.3. Dynamic Equations of the System

We use the Lagrange method to derive the dynamics of the system and accord with the method of [14]. Then, the form of Lagrange's equations will be as follows.

For the joint variable q_1j, the following simply joint variable equation:

\frac{d}{d t} (\frac{\partial K}{\partial {\dot{q}}_{1 j}}) - \frac{\partial K}{\partial q_{1 j}} + \frac{\partial V}{\partial q_{1 j}} = τ_{j} (j = 1,2, \dots, n) .

(40)

For the joint variable q_2j, the following simply joint variable equation:

\frac{d}{d t} (\frac{\partial K}{\partial {\dot{q}}_{2 j}}) - \frac{\partial K}{\partial q_{2 j}} + \frac{\partial V}{\partial q_{2 j}} = 0 (j = 1,2, \dots, n) .

(41)

For the deformation variable δ_jf, the following simply deformation variable equation:

\begin{matrix} \frac{d}{d t} (\frac{\partial K}{\partial {\dot{δ}}_{j f}}) - \frac{\partial K}{\partial δ_{j f}} + \frac{\partial V}{\partial δ_{j f}} = 0 \\ (j = 1,2, \dots, n; f = 1,2, \dots, N_{j}) . \end{matrix}

(42)

Upon the substitution of the system kinetic energy of (15), the system elastic potential energy of (31), and the gravity energy of (37) into (40), (41), and (42), thus becoming as follows.

For the joint variable q_1j,

J_{r j} n_{i}^{2} {\overset{‥}{q}}_{1 j} = - K t_{j} q_{1 j} + K t_{j} q_{2 j} + τ_{j} (j = 1,2, \dots, n) .

(43)

For the joint variables q_2j (j = 1, 2, …, n),

2 \sum_{i = j}^{n} Tr {\frac{\partial W_{i}}{\partial q_{2 j}} [B_{3 i} {\overset{‥}{W}}_{i}^{T} + \sum_{k = 1}^{N_{i}} H_{i k} W_{i}^{T} {\overset{‥}{δ}}_{i k} + \sum_{k = 1}^{N_{i}} \sum_{l = 1}^{N_{i}} L_{i k l} W_{i}^{T} {\dot{δ}}_{i k} {\dot{δ}}_{i l} + 2 \sum_{k = 1}^{N_{i}} H_{i k} {\dot{W}}_{i}^{T} {\dot{δ}}_{i k}]} = g^{T} \sum_{i = j}^{n} \frac{\partial W_{i}}{\partial q_{2 j}} r_{i} + K t_{j} q_{1 j} - K t_{j} q_{2 j} .

(44)

For the deformation variables δ_jf (j = 1, 2, …, n; f = 1, 2, …, N_j),

2 \sum_{i = j + 1}^{n} Tr {\frac{\partial W_{i}}{\partial δ_{j f}} [B_{3 i} {\overset{‥}{W}}_{i}^{T} + \sum_{k = 1}^{N_{i}} H_{i k} W_{i}^{T} {\overset{‥}{δ}}_{i k} + \sum_{k = 1}^{N_{i}} \sum_{l = 1}^{N_{i}} L_{i k l} W_{i}^{T} {\dot{δ}}_{i k} {\dot{δ}}_{i l} + 2 \sum_{k = 1}^{N_{i}} H_{i k} {\dot{W}}_{i}^{T} {\dot{δ}}_{i k}]} + 2 Tr {[{\overset{‥}{W}}_{j} H_{j f} + \sum_{k = 1}^{N_{j}} W_{j} N_{j f k} {\overset{‥}{δ}}_{j k} + \sum_{k = 1}^{N_{j}} \sum_{l = 1}^{N_{j}} W_{j} P_{j f k l} {\dot{δ}}_{j k} {\dot{δ}}_{j l} + 2 \sum_{k = 1}^{N_{j}} {\dot{W}}_{j} N_{j f k} {\dot{δ}}_{j k}] W_{j}^{T}} + 2 \sum_{k = 1}^{N_{j}} {\overset{‥}{δ}}_{j k} T_{j f k} = g^{T} \sum_{i = j + 1}^{n} (\frac{\partial W_{i}}{\partial δ_{j f}} r_{i}) + g^{T} W_{j} ∊_{j f} - g^{T} W_{j} \sum_{k = 1}^{N_{j}} δ_{j k} κ_{j k l} - \sum_{k = 1}^{N_{j}} δ_{j k} K_{j f k} .

(45)

Here,

\begin{matrix} H_{i j} = D_{i j} + \sum_{k = 1}^{N_{i}} δ_{i k} D_{i j k}^{T} + \sum_{k = 1}^{N_{i}} δ_{i k} (- E_{i j k}) + \sum_{k = 1}^{N_{i}} \sum_{l = 1}^{N_{i}} δ_{i k} δ_{i l} (- \frac{1}{2} F_{i j k l}^{T} - F_{i l j k}) + \sum_{k = 1}^{N_{i}} \sum_{l = 1}^{N_{i}} \sum_{s = 1}^{N_{i}} δ_{i k} δ_{i l} δ_{i s} (\frac{1}{2} E_{i j l k s}^{T}), \\ L_{i j k} = - E_{i j k} + \sum_{l = 1}^{N_{i}} δ_{i l} (- F_{i l j k}) + \sum_{l = 1}^{N_{i}} \sum_{s = 1}^{N_{i}} δ_{i l} δ_{i s} (\frac{1}{2} E_{i j k l s}^{T} + E_{i j l k s}^{T} - E_{i j l k s}), \\ N_{i j k} = D_{i j k}^{T} + \sum_{l = 1}^{N_{i}} δ_{i l} (- F_{i j k l}^{T} - F_{i k j l}) + \sum_{l = 1}^{N_{i}} \sum_{s = 1}^{N_{i}} δ_{i l} δ_{i s} E_{i j l k s}^{T}, \\ P_{i j k l} = - F_{i j k l}^{T} + \sum_{s = 1}^{N_{i}} δ_{i s} E_{i j s k l}^{T} . \end{matrix}

(46)

Finally, by complicated derivation and assembling, and using the recursive scheme to reduce the number of calculation, the dynamic equations of the system are obtained in the following formulation:

J \overset{‥}{z} = R^{I},

(47)

where the generalized coordinates column z is defined as

z = {[\begin{bmatrix} q_{11} & q_{21} & δ_{11} & δ_{12} & \cdot & δ_{1 N_{1}} & q_{12} & q_{22} & δ_{21} & \cdot & δ_{2 N_{2}} & \cdot & q_{1 n} & q_{2 n} & δ_{n 1} & \cdot & δ_{n N_{n}} \end{bmatrix}]}^{T},

(48)

where J is the inertia matrix of the system composed of coefficients of the generalized acceleration $\overset{‥}{z}$ in the system dynamic equations. It is positive-definite and symmetric. The elements of J are arranged in the order of the generalized coordinate column z, the J_i are diagonal elements, the ^jJ_i are off-diagonal elements, and these elements are given as follows:

J_{i} = (\begin{pmatrix} J_{r i} n_{i}^{2} & 0 & 0 & \cdot & 0 \\ 0 & J_{i i} & J_{i i 1} & \cdot & J_{i i N_{i}} \\ 0 & I_{i 1 i} & I_{i 1 i 1} & \cdot & I_{i 1 i N_{i}} \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & I_{i N_{i} i} & I_{i N_{i} i 1} & \cdot & I_{i N_{i} i N_{i}} \end{pmatrix}),

(49)

^{j} J_{i} = (\begin{pmatrix} 0 & 0 & 0 & \cdot & 0 \\ 0 & J_{i j} & J_{i j 1} & \cdot & J_{i j N_{j}} \\ 0 & I_{i 1 j} & I_{i 1 j 1} & \cdot & I_{i 1 j N_{j}} \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & I_{i N_{i} j} & I_{i N_{i} j 1} & \cdot & I_{i N_{i} j N_{j}} \end{pmatrix}) .

(50)

(I) In the joint variable equation j (j = 1, 2, …, n), the inertia coefficient J_jh of the second derivative ${\overset{‥}{q}}_{2 h}$ of the joint variable q_2h is

J_{j h} = 2 Tr {({\hat{W}}_{j - 1} U_{j}) {}^{j} {\tilde{F}}_{h} U_{h}^{T} {\hat{W}}_{h - 1}^{T}},

(51)

where

^{j} {\tilde{F}}_{h} = \sum_{i = \max (h, j)}^{n}^{j} {\tilde{W}}_{i} {B_{3 i}}^{h} {\tilde{W}}_{i}^{T},

(52)

^{h} {\tilde{W}}_{i} = E_{h} A_{h + 1} \cdot E_{i - 1} A_{i} .

(53)

Here, ${}^{j}{\tilde{W}}_{i}$ is the 4 × 4 homogeneous transformation matrix from the proximal reference frame (X_bY_bZ_b)_j of link j to the proximal reference frame (X_bY_bZ_b)_i of link i. From (52), we have J_jh = J_hj (symmetric).

(II) In the joint variable equation j (j = 1, 2, …, n), the inertia coefficient J_jhk of the second derivative ${\overset{‥}{δ}}_{h k}$ of the joint variable δ_hk (for 1 ≤ k ≤ N_k) is expressed as follows:

for h = n; j = 1, …, n,

J_{j h k} = 2 Tr {({\hat{W}}_{j - 1} U_{j}) {}^{j}{\tilde{W}}_{n} H_{n k} W_{n}^{T}};

(54)

for h = j, …, n – 1; j = 1, …, n – 1,

J_{j h k} = 2 Tr {{\hat{W}}_{j - 1} U_{j} [{}^{j}{F_{h}} M_{h k}^{T} + {}^{j}{\tilde{W}}_{h} H_{h k}] W_{h}^{T}};

(55)

for h = 1, …, j – 1; j = 2, …, n,

J_{j h k} = 2 Tr {({\hat{W}}_{j - 1} U_{j}) {}^{j}{F_{h}} M_{h k}^{T} W_{h}^{T}} .

(56)

Here,

\begin{matrix} ^{j} F_{h} = \sum_{i = \max (h + 1, j)}^{n}^{j} {\tilde{W}}_{i} {B_{3 i}}^{h} W_{i}^{T}, \\ ^{h} W_{i} = A_{h + 1} E_{h + 1} ˙ E_{i - 1} A_{i}, \end{matrix}

(57)

where ${}^{j}{W_{i}}$ is the 4 × 4 homogeneous transformation matrix from the distal coordinate system (X_dY_dZ_d)_j of link j to the proximal coordinate system (X_bY_bZ_b)_i of link i. Note that the inertia coefficient for the deformation variable δ_hk in the joint equation j is the same as that for the joint variable q_2h in the deformation equation h, k. This further extends the symmetry of the inertia matrix and reduces the computation necessary.

(III) In the joint variable equation j (j = 1, 2, …, n), the inertia coefficients I_jfh of the second derivative ${\overset{‥}{q}}_{2 h}$ of the joint variable q_2h are expressed as follows:

for j = n; h = 1, …, n,

I_{j f n} = 2 Tr {({\hat{W}}_{h - 1} U_{h}) {}^{h}{\tilde{W}}_{n} H_{n f} W_{n}^{T}};

(58)

for j = 1, …, n – 1; h = 1, …, j,

\begin{matrix} I_{j f h} = 2 Tr {{\hat{W}}_{h - 1} U_{h} [{}^{j}{\tilde{φ}}_{h}^{T} M_{j f}^{T} + {}^{h}{\tilde{W}}_{j} H_{j f}] W_{j}^{T}}; \end{matrix}

(59)

for j = 1, …, n – 1; h = j + 1, …, n – 1,

I_{j f h} = 2 Tr {(W_{j} M_{j f}) {}^{j}{\tilde{φ}}_{h} U_{h}^{T} {\hat{W}}_{h - 1}^{T}} .

(60)

(IV) In the deformation variable equation j, f, the inertia coefficients I_jfhk of the second derivative ${\overset{‥}{δ}}_{h k}$ of the deformation variable δ_hk (1 ≤ f ≤ N_j, 1 ≤ k ≤ N_h) are expressed as follows:

for j = h = n,

I_{n f n k} = 2 Tr {N_{n f k}} + 2 T_{n f k};

(61)

for j = n; h = 1, …, n – 1,

I_{n f h k} = 2 Tr {(W_{h} M_{h k}) {}^{h}{W_{n}} H_{n f} W_{n}^{T}};

(62)

for j = 1, …, n – 1; h = n,

I_{j f n k} = 2 Tr {(W_{j} M_{j f}) {}^{j}{W_{n}} H_{n k} W_{n}^{T}};

(63)

for j = h; h = 1, …, n – 1,

I_{j f h k} = 2 Tr {{M_{j f}}^{j} ϕ_{h} M_{j k}^{T} + N_{j f k}} + 2 T_{j f k};

(64)

for j = 2, …, n – 1; h = 1, …, j – 1,

I_{j f h k} = 2 Tr {W_{j} [{M_{j f}}^{j} ϕ_{h} + {H_{j f}^{T}}^{h} W_{j}^{T}] M_{h k}^{T} W_{h}^{T}};

(65)

for j = 1, …, n – 2; h = j + 1, …, n – 1,

I_{j f h k} = 2 Tr {W_{j} M_{j f} [{}^{j}{φ_{h}} M_{h k}^{T} + {}^{j}{W_{h}} H_{h k}] W_{h}^{T}} .

(66)

Here,

\begin{matrix} ^{j} {\tilde{ϕ}}_{h} = \sum_{i = \max (j + 1, h)}^{n}^{j} W_{i} {B_{3 i}}^{h} {\tilde{W}}_{i}^{T}, \\ ^{j} ϕ_{h} = \sum_{i = \max (j + 1, h + 1)}^{n}^{j} W_{i} {B_{3 i}}^{h} W_{i}^{T} . \end{matrix}

(67)

For the symmetry of the coefficients, we have I_jfhk = I_hkjf.

(V) R^I is generalized force, containing the remaining dynamics and external forcing terms except impact force as follows:

R^{I} = {[\begin{bmatrix} R_{r 1}^{I} & R_{J 1}^{I} & R_{11}^{I} & \cdot & R_{1 N_{1}}^{I} & R_{r 2}^{I} & R_{J 2}^{I} & R_{21}^{I} & \cdot & R_{2 N_{2}}^{I} & \cdot & R_{r n}^{I} & R_{J n}^{I} & R_{n 1}^{I} & \cdot & R_{n N_{n}}^{I} \end{bmatrix}]}^{T} .

(68)

The R_rj^I is the remains of (43) with the second derivatives removed as follows:

R_{r j}^{I} = - K t_{j} q_{1 j} + K t_{j} q_{2 j} + τ_{j} (j = 1, \dots, n) .

(69)

The R_Jj^I is the remains of (44) with the second derivatives removed as follows:

for j = 1,

R_{J 1}^{I} = - 2 Tr {U_{1} Q_{1}} + g^{T} U_{1} P_{1} + K t_{1} q_{11} - K t_{1} q_{21};

(70)

for j = 2, …, n,

R_{J j}^{I} = - 2 Tr {{\hat{W}}_{j - 1} U_{j} Q_{j}} + g^{T} {\hat{W}}_{j - 1} U_{j} P_{j} + K t_{j} q_{1 j} - K t_{j} q_{2 j} .

(71)

The R_jf^I is the remains of (45) with the second derivatives removed as follows.

For j = 1, …, n – 1,

R_{J j}^{I} = - 2 Tr {W_{j} M_{j f} A_{j + 1} Q_{j + 1} + [{\overset{‥}{W}}_{v j} H_{j f} + 2 {\dot{W}}_{j} (\sum_{k = 1}^{N_{j}} N_{j f k} {\dot{δ}}_{j k}) + W_{j} (\sum_{k = 1}^{N_{j}} \sum_{l = 1}^{N_{j}} P_{j f k l} {\dot{δ}}_{j k} {\dot{δ}}_{j l})] \times W_{j}^{T}} - \sum_{k = 1}^{N_{j}} δ_{j k} K_{j f k} + g^{T} W_{j} M_{j f} A_{j + 1} P_{j + 1} + g^{T} W_{j} ∊_{j f} - g^{T} W_{j} \sum_{k = 1}^{N_{j}} δ_{j k} κ_{j f k};

(72)

for j = n,

R_{n f}^{I} = - 2 Tr [{\overset{‥}{W}}_{v n} H_{n f} + 2 {\dot{W}}_{n} (\sum_{k = 1}^{N_{n}} N_{n f k} {\dot{δ}}_{n k}) + W_{j} (\sum_{k = 1}^{N_{n}} \sum_{l = 1}^{N_{n}} P_{n f k l} {\dot{δ}}_{n k} {\dot{δ}}_{n l}) W_{n}^{T}] - \sum_{k = 1}^{N_{n}} δ_{n k} K_{n f k} + g^{T} W_{n} ∊_{n f} - g^{T} W_{n} \sum_{k = 1}^{N_{n}} δ_{n k} κ_{n f k} .

(73)

Here,

\begin{matrix} Q_{j} = \sum_{i = j}^{n} (^{j} {\tilde{W}}_{i} B_{3 i} {\overset{‥}{W}}_{v i} + 2 \sum_{k = 1}^{N_{i}}^{j} {\tilde{W}}_{i} H_{i k} {\dot{W}}_{i}^{T} {\dot{δ}}_{i k} + \sum_{k = 1}^{N_{i}} \sum_{l = 1}^{N_{i}} {}^{j}{\tilde{W}}_{i} L_{i k l} W_{i}^{T} {\dot{δ}}_{i k} {\dot{δ}}_{i l}), \\ P_{j} = \sum_{i = j}^{n} {}^{j}{\tilde{W}}_{i} r_{i}, \end{matrix}

(74)

where the value of ${\overset{‥}{W}}_{v i}$ is the remainder of ${\overset{‥}{W}}_{i}$ by only eliminating the terms involving ${\overset{‥}{q}}_{2 h}$ and ${\overset{‥}{δ}}_{h k}$ , and it can be calculated recursively by

{\overset{‥}{W}}_{i} = \sum_{h = 1}^{i} {\hat{W}}_{h - 1} {U_{h}}^{h} {\tilde{W}}_{i} {\overset{‥}{q}}_{2 h} + \sum_{h = 1}^{i - 1} \sum_{k = 1}^{N_{h}} W_{h} {M_{h k}}^{h} W_{i} {\overset{‥}{δ}}_{h k} + {\overset{‥}{W}}_{v i} .

(75)

The value of ${\overset{‥}{\hat{W}}}_{i}$ can be obtained similarly as follows:

{\overset{‥}{\hat{W}}}_{i} = \sum_{h = 1}^{i} ({\hat{W}}_{h - 1} {U_{h}}^{h} {\bar{W}}_{i} {\overset{‥}{q}}_{2 h} + \sum_{k = 1}^{N_{h}} W_{h} {M_{h k}}^{h} {\hat{W}}_{i} {\overset{‥}{δ}}_{h k}) + {\overset{‥}{\hat{W}}}_{v i},

(76)

where

\begin{matrix} ^{h} {\bar{W}}_{i} = E_{h} A_{h + 1} \cdot A_{i} E_{i}, \\ ^{h} {\hat{W}}_{i} = A_{h + 1} E_{h + 1} \cdot A_{i} E_{i} . \end{matrix}

(77)

For accelerating the computation of the generalized inertia matrix of the system, the terms ${\overset{‥}{W}}_{v i}$ , ${\overset{‥}{\hat{W}}}_{v i}$ , Q_j, P_j, ${}^{j}{\tilde{F}}_{h}$ , $^{j} {\tilde{ϕ}}_{h}$ , ^hF_j, and ${}^{j}{ϕ_{h}}$ defined in (52)–(76), respectively, can be calculated recursively as follows:

\begin{matrix} {\overset{‥}{W}}_{v 1} = U_{21} {\dot{q}}_{21}^{2}, \\ {\overset{‥}{W}}_{v i} = {\overset{‥}{\hat{W}}}_{v i - 1} A_{i} + 2 {\dot{\hat{W}}}_{i - 1} U_{i} {\dot{q}}_{2 i} + {\hat{W}}_{i - 1} U_{2 i} {\dot{q}}_{2 i}^{2}, \\ {\overset{‥}{\hat{W}}}_{v i} = {\overset{‥}{W}}_{v i} E_{i} + 2 {\dot{W}}_{i} {\dot{E}}_{i}, \\ Q_{n} = B_{3 n} {\overset{‥}{W}}_{v n}^{T} + 2 (\sum_{k = 1}^{N_{n}} H_{n k} {\dot{δ}}_{n k}) {\dot{W}}_{n}^{T} + (\sum_{k = 1}^{N_{n}} \sum_{l = 1}^{N_{n}} L_{n k l} {\dot{δ}}_{n k} {\dot{δ}}_{n l}) W_{n}^{T}, \\ Q_{j} = B_{3 i} {\overset{‥}{W}}_{v i}^{T} + 2 (\sum_{k = 1}^{N_{j}} H_{j k} {\dot{δ}}_{j k}) {\dot{W}}_{j}^{T} + (\sum_{k = 1}^{N_{j}} \sum_{l = 1}^{N_{j}} L_{j k l} {\dot{δ}}_{j k} {\dot{δ}}_{j l}) W_{j}^{T} + {}^{j}{\tilde{W}}_{j + 1} Q_{j + 1}, \end{matrix}

(78)

P_{n} = r_{n} = M_{n} r_{r n} + \sum_{k = 1}^{N_{n}} δ_{n k} ∊_{n k} - \frac{1}{2} \sum_{k = 1}^{N_{n}} \sum_{l = 1}^{N_{n}} δ_{n k} δ_{n l} κ_{n k l},

(79)

P_{j} = r_{j} + {}^{j}{\tilde{W}}_{j + 1} P_{j + 1} .

(80)

As for ${}^{j}{\tilde{F}}_{h}$ , for j = h = n,

^{n} {\tilde{F}}_{n} = B_{3 n},

(81)

for j = n; h = 1, …, n – 1,

^{n} {\tilde{F}}_{h} = {B_{3 i}}^{h} {\tilde{W}}_{n}^{T},

(82)

for j = h; h = 1, …, n – 1,

^{h} {\tilde{F}}_{h} = B_{3 h} + E_{h} {A_{h + 1}}^{h + 1} {\tilde{F}}_{h} .

(83)

And for j ≠ h ≠ n,

\begin{matrix} ^{j} {\tilde{F}}_{h} = B_{3 j} \cdot {}^{h}{\tilde{W}}_{j}^{T} + (E_{i} A_{j + 1}) \cdot {}^{j + 1}{\tilde{F}}_{h}, \\ ^{j} {\tilde{F}}_{h} = {({}^{h}{\tilde{F}}_{j})}^{T} . \end{matrix}

(84)

As for ${}^{j}{F_{h}}$ , only if N_h > 0, calculate

{}^{j}{F_{h}} = {}^{j}{\tilde{F}}_{h + 1} \cdot A_{h + 1}^{T} (j = 1, \dots, n; h = 1, \dots, n - 1) .

(85)

As for $^{j} {\tilde{ϕ}}_{h}$ , it can be obtained from ${}^{h}{F_{j}}$ easily as follows:

^{j} {\tilde{φ}}_{h} = {({}^{h}{F_{j}})}^{T} .

(86)

As for ^jϕ_h, only if N_j > 0 and N_h > 0, calculate

^{j} φ_{h} = {}^{j}{\tilde{φ}}_{h + 1} \cdot A_{h + 1}^{T} (j = 1, \dots, n - 1; h = j, \dots, n - 1) .

(87)

4. Example of Dynamic Simulation

A general-purpose software package for dynamic simulation of flexible-links and flexible-joints robotic manipulators based on the said algorithm is written in C++. The data structure is pointers and linked lists, so that the data length can be changed based on different initial conditions of the system. The first step is setting up a given interface to input system parameters via reading files. Then, the basic module matrix can be generated automatically. With that, reducing a second order differential equation to a first order linear ordinary differential equation. Eventually, a process based on Adams predictor-corrector method for solving ordinary differential equation calculates till the end time, with an output file per step. It is universal for usage, computationally efficient, and numerically stable.

4.1. Flexible Single-Link Undergoing Low/High Rotational Speed

Ignoring the gravity and the flexible effect of the joint, here, consider the flexible single-link undergoing low/high rotational speed in order to discover the applicability of the high-order geometric nonlinearity coupling model. The properties of the flexible beam are the same as those in [3, 8], and they are given as follows. The length L₁ = 8 m, the cross-section area A_x1 = 7.2968 × 10⁻⁵ m², the area moment of inertia I = 8.2189 × 10⁻⁹ m⁴, the mass density ρ = 2.7667 × 10³ kg/m³, and the elasticity modulus E₁ = 6.8952 × 10¹⁰ N/m². For the flexible cantilever beam without large overall motion, the fundamental frequency is 0.464 Hz. The large rotating motion law of the system adopted $\dot{θ}$ in [3, 8] is given by

\dot{θ} = {\begin{cases} \frac{ω_{0}}{T} t - \frac{ω_{0}}{2 π} sin (\frac{2 π}{T} t), & 0 \leq t \leq T, \\ ω_{0}, & t > T, \end{cases}

(88)

where T = 15 s, and ω₀ is chosen here to be 0.6, 2, 3, 8, and 10 rad/s, respectively, in the following numerical simulations. The tip deformations of link u_y are shown from Figures 3, 4, and 5.

Figure 3:

ω₀ < 3 rad/s.

Figure 4:

ω₀ = 3 rad/s.

Figure 5:

ω₀ > 3 rad/s.

It can be observed that, undergoing a low rotating speed (0.6 rad/s or 2 rad/s), there is an obvious difference between zero-order and high-order coupling model, and the difference becomes larger and larger with the increasing of the rotating speed. When the high-order coupling model is at the rotating speed 3 rad/s (≈0.477 Hz), closing to the fundamental frequency (0.464 Hz), the vibrating amplitude of the link tip is aroused much higher, but it is convergent. Whereas, the results form zero-order coupling model turn to be divergente. When reaching a high rotating speed, 8 rad/s (≈1.272 Hz) or 10 rad/s (≈1.59 Hz), the zero-order coupling model will fail to get the convergent results. These phenomena are coincided with the numerical result in [2, 11].

4.2. Spatial Flexible-Link and Flexible-Joint Robotic Arm

In order to validate the algorithm and package presented in this paper, the dynamic simulation of Canadarm2 (the Space Station remote manipulator system serving on the international space station) is given as an example.

As shown in Figure 6, Canadarm2 consists of four parts: shoulder, upper arm boom, lower arm boom, and wrist, in which the shoulder and the wrist both have 3 rotational degrees of freedom (DOF), and the elbow connecting the upper arm and the lower arm has 1 rotational DOF. Hence, the total system can be regarded as a 7-DOF manipulator with seven links and seven revolute joints. The links making up the shoulder and the wrist are short and thick in shape; therefore, they can be regarded as rigid links, whereas the upper and lower arms are assumed to be flexible due to their slender beam type. Therefore, Canadarm2 system is simplified to contain two flexible links and three rigid links which are connected through flexible joints, as shown in Figure 7. In this example, there is a distance d = 0.475 m between the 2nd link and 1st link along the Z-axis, which consequently arouse the rotation of the flexible links. And the longitudinal stretching (along X-axis), transversal bending (along Y-axis and Z-axis), and torsion (around X-axis) are considered here. Correspondingly, a total of 12 eigenfunctions (the first three preorientation of deformation) of a clamped-free beam are utilized to represent the mode shapes of each link deformation. The links parameters are shown in Table 1. The material of the links is aluminum, cross-section parameter E_iI_yi = E_iI_zi = 3.8 × 10⁶ Nm², E_iA_i = 2.84 × 10⁷ N, G_i = 0.27 × 10⁵ MPa, and J_xi = 15.41 kg × m². Initial joint angle is q₁₁ = 0.4π rad, q₂₁ = 0.4π rad, q₁₂ = – 0.4π rad, and q₂₂ = – 0.4π rad, respectively. Initial joint angle velocity is zero.

Table 1:

The parameters of space robot system.

Link	Mass (kg)	Centre (m)	Transverse parameters (kg × m⁻²)	Length (m)

1	139.35	(3.1883, 0.0, 0.0282)	I_xx = 3.225, I_yy = 998.630, I_zz = 997.180, I_yx = – 1.928, I_zx = 9.185, I_yz = – 0.0302	(6.377, 0, 0.1524)
2	87.42	(3.53, 0.0, −0.0066)	I_xx = 1.59, I_yy = 632.9, I_zz = 632.51, I_yx = 0.8021, I_zx = – 3.786, I_yz = – 0.08	(7.06, 0, 0.1524)
3	8.48	(0.2286, 0.0, 0.0)	I_xx = 0.0893, I_yy = 0.402, I_zz = 0.402, I_yx = 0.0, I_zx = 0.0, I_yz = 0.0	(0.4572, 0, 0)
4	45.94	(0.381, 0.0, 0.0)	I_xx = 0.408, I_yy = 4.86,I_zz = 4.86, I_yx = 0.0, I_zx = 0.0, I_yz = 0.0	(0.762, 0, 0)
5	45.21	(0.332, 0.0, 0.0)	I_xx = 1.903, I_yy = 3.816,I_zz = 2.823, I_yx = 0.385, I_zx = 0.966, I_yz = 0.353	(0.6604, 0, 0)

Figure 6:

The space robot arm Canadarm2.

Figure 7:

A spatial flexible robot arms.

4.2.1. Dropping Motion in the Gravity Field

Figure 8 describes the motion of the high-order coupling flexible five-links system dropping in the gravity field.

Figure 8:

The flexible robot arms motion.

4.2.2. Effects of Joint Flexibility on the System

Figure 9 shows the torsional angle of the 1st joint ∊₁ in the high-order coupling model with different mass J_rin_i², 5, 10, or 15, respectively. Figures 10, 11, 12, and 13 show the X, Y, and Z deformation u_x, v_y, w_z, and torsional angle θ_x of 1st link tip with different joint mass. The larger joint's mass is, the larger joint angles of motion the system has. The same phenomena can be seen form Figures 10 to 13.

Figure 9:

The torsional angle of 1st joint, ∊₁.

Figure 10:

u_x of 1st link tip.

Figure 11:

v_y of 1st link tip.

Figure 12:

w_z of 1st link tip.

Figure 13:

θ_x of 1st link tip.

The real rotational angle q_2i of five joint in the high-order coupling model is shown in Figures 14, 15, 16, 17, and 18 with different stiffness. The joint stiffness Kt_i is chosen here to be 0.60, 1.33, 2.66, and 5.00 × 10⁶ N/mm, respectively.

Figure 14:

The 1st joint, q₂₁.

Figure 15:

The 2nd joint, q₂₂.

Figure 16:

The 3rd joint, q₂₃.

Figure 17:

The 4th joint, q₂₄.

Figure 18:

The 5th joint, q₂₅.

5. Conclusions

This paper aims at building a comprehensive model of spatial manipulator arms consisting of n flexible links and n flexible joints. The deformations of stretching, bending, and torsion of the links are all considered. And the flexibility and the mass of the joints are also included. In the modeling, the so-called dynamic stiffening effects are considered via the adoption of nonlinear deformation description. Based on the presented method, a general-purpose software package of multilink manipulator arms is developed. The results demonstrate that dynamic stiffening effects cannot be neglected. The coupling effects of flexible links and flexible joints on dynamic characteristics are significant to the system's accurate dynamic response predictions, performance evaluation, and implementation of procedure control.

Footnotes

Acknowledgments

This work is supported by the National Natural Science Foundations of China (11272155, 11132007, and 10772085), 333 Project of Jiangsu Province (BRA2011172), and the Fundamental Research Funds for Central Universities (30920130112009).

References

Shabana

A. A.

Gantoi

F. M.

, and Brown

M. A.

, “Integration of finite element and multibody system algorithms for the analysis of human body motion,” in Proceedings of the IUTAM Symposium on Human Body Dynamics: From Multibody Systems to Biomechanics, vol. 2, pp. 233–240, June 2011.

Kane

T. R.

Ryan

R. R.

, and Banerjee

A. K.

, “Dynamics of a cantilever beam attached to a moving base,” Journal of Guidance, Control, and Dynamics, vol. 10, no. 2, pp. 139–151, 1987.

S.-C.

and Haug

E. J.

, “Geometric non-linear substructuring for dynamics of flexible mechanical systems,” International Journal for Numerical Methods in Engineering, vol. 26, no. 10, pp. 2211–2226, 1988.

and Fenton

R. G.

, “Coupling effect of a flexible link and a flexible joint,” International Journal of Robotics Research, vol. 13, no. 5, pp. 443–453, 1994.

Wallrapp

and Schwertassek

, “Representation of geometric stiffening in multibody system simulation,” International Journal for Numerical Methods in Engineering, vol. 32, no. 8, pp. 1833–1850, 1991.

Low

K. H.

, “Vibration analysis of a tip-loaded beam attached to a rotating joint,” Computers and Structures, vol. 52, no. 5, pp. 955–968, 1994.

Yoo

H. H.

Ryan

R. R.

, and Scott

R. A.

, “Dynamics of flexible beams undergoing overall motions,” Journal of Sound and Vibration, vol. 181, no. 2, pp. 261–278, 1995.

Ryu

Kim

S. S.

, and Kim

S.-S.

, “A criterion on inclusion of stress stiffening effects in flexible multibody dynamic system simulation,” Computers and Structures, vol. 62, no. 6, pp. 1035–1048, 1997.

El-Absy

and Shabana

A. A.

, “Geometric stiffness and stability of rigid body modes,” Journal of Sound and Vibration, vol. 207, no. 4, pp. 465–496, 1997.

10.

Al-Bedoor

B. O.

and Hamdan

M. N.

, “Geometrically non-linear dynamic model of a rotating flexible arm,” Journal of Sound and Vibration, vol. 240, no. 1, pp. 59–72, 2001.

11.

Liu

J. Y.

and Hong

J. Z.

, “Geometric stiffening effect on rigid flexible coupling dynamics of an elastic beam,” Journal of Sound and Vibration, vol. 278, no. 4–5, pp. 1147–1162, 2004.

12.

Liu

Z.-Y.

Hong

J.-Z.

, and Liu

J.-Y.

, “Complete geometric nonlinear formulation for rigid-flexible coupling dynamics,” Journal of Central South University of Technology, vol. 16, no. 1, pp. 119–124, 2009.

13.

Schwertassek

Wallrapp

, and Shabana

A. A.

, “Flexible multibody simulation and choice of shape functions,” Nonlinear Dynamics, vol. 20, no. 4, pp. 361–380, 1999.

14.

Book

W. J.

, “Recursive Lagrangian dynamics of flexible manipulator arms,” International Journal of Robotics Research, vol. 3, no. 3, pp. 87–101, 1984.

15.

Changizi

and Shabana

A. A.

, “Recursive formulation for the dynamic analysis of open loop deformable multibody systems,” Journal of Applied Mechanics, Transactions ASME, vol. 55, no. 3, pp. 687–693, 1988.

16.

Kim

S.-S.

and Haug

E. J.

, “A recursive formulation for flexible multibody dynamics—part I: open-loop systems,” Computer Methods in Applied Mechanics and Engineering, vol. 71, no. 3, pp. 293–314, 1988.

17.

Jain

and Rodriguez

, “Recursive flexible multibody system dynamics using spatial operators,” Journal of Guidance, Control, and Dynamics, vol. 15, no. 6, pp. 1453–1466, 1992.

18.

Hwang

Y.-L.

, “Recursive Newton-Euler formulation for flexible dynamic manufacturing analysis of open-loop robotic systems,” International Journal of Advanced Manufacturing Technology, vol. 29, no. 5, pp. 598–604, 2006.

19.

Zhang

D.-G.

and Zhou

S.-F.

, “Dynamic analysis of flexible-link and flexible-joint robots,” Applied Mathematics and Mechanics, vol. 27, no. 5, pp. 695–704, 2006.

20.

Zhang

D.-G.

, “Recursive Lagrangian dynamic modeling and simulation of multi-link spatial flexible manipulator arms,” Applied Mathematics and Mechanics, vol. 30, no. 10, pp. 1283–1294, 2009.

21.

Spong

M. W.

, “Modeling and control of elastic joint robots,” Journal of Dynamic Systems, Measurement and Control, Transactions of the ASME, vol. 109, no. 4, pp. 310–319, 1987.

22.

Wasfy

T. M.

and Noor

A. K.

, “Computational strategies for flexible multibody systems,” Applied Mechanics Reviews, vol. 56, no. 6, pp. 553–613, 2003.

23.

Dwivedy

S. K.

and Eberhard

, “Dynamic analysis of flexible manipulators, a literature review,” Mechanism and Machine Theory, vol. 41, no. 7, pp. 749–777, 2006.

24.

K.-S.

and Kim

J.-H.

, “Deployment of a multi-link flexible structure,” Journal of Sound and Vibration, vol. 294, no. 1–2, pp. 298–313, 2006.

25.

K. S.

Gonzalez

R. C.

, and Lee

C. S. G.

, Robotics: Control, Sensing, Vision, and Intelligence, McGraw-Hill, New York, NY, USA, 1987.

Recursive Formulation for Dynamic Modeling and Simulation of Multilink Spatial Flexible Robotic Manipulators

Abstract

1. Introduction

2. Kinematics of Flexible Robots

2.1. Simplified Model of Flexible Joint

2.2. Simplified Model of Flexible Link

2.3. Coordinate Systems and Transformation Matrices

2.4. Velocity of a Point of Link i

3. Dynamics of Flexible Robots

3.1. The System Kinetic Energy

3.2. The System Potential Energy

3.3. Dynamic Equations of the System

4. Example of Dynamic Simulation

4.1. Flexible Single-Link Undergoing Low/High Rotational Speed

4.2. Spatial Flexible-Link and Flexible-Joint Robotic Arm

4.2.1. Dropping Motion in the Gravity Field

4.2.2. Effects of Joint Flexibility on the System

5. Conclusions

Footnotes

Acknowledgments

References