A new version of the Riccati transfer matrix method for multibody systems consisting of chain and branch bodies

Computational speed and stability are two important aspects in the dynamics analysis of large-scale complex multibody systems. In order to improve both in the context of the multibody system transfer matrix method, a new version of the Riccati transfer matrix method is presented. Based on the new version of the general transfer matrix method for multibody system simulation, recursive formulae are developed which not only retain all advantages of the transfer matrix method, but also reduce the truncation error. As a result, the computational speed, accuracy and efficiency are improved. Numerical computation results obtained by the proposed method and an ordinary multibody system formulation show good agreement. The successful computation for a spatial branch system with more than 100000 degrees of freedom validates that the proposed method is also working for huge systems.


Introduction
Multibody system methods are important for design and analysis of modern complex mechanical systems [1][2][3][4], where improving computational speed and stability of algorithms has always been one of the key research directions. In order to increase computation speed and simplify modeling, one of the authors and his co-workers presented and perfected constantly the transfer matrix method for multibody systems (MSTMM) [4][5][6][7]. Especially, the Riccati transformation was introduced to deal with long body chains. Recently, the form of B X. Rui harukarui@qq.com 1 state vector has been changed from displacements and angles to translational and angular accelerations to allow for arbitrary integration algorithms [6], and the automatic deduction theorem was presented to obtain the overall transfer equation more easily [7], which now calls for an adaption of the Riccati approach to this new transfer matrix method.
In principle, the overall transfer matrix of a system is obtained by multiplying its element transfer matrices. Thereby truncation errors are accumulated which may lead to numerical failure of the algorithm. In order to improve computational stability of the classical transfer matrix method (TMM), Horner [8] introduced the Riccati transfer matrix method (RTMM) which turns the boundary conditions of two points in the original differential equation into an initial condition of one point only. As a result, the numerical stability is improved efficiently, where additionally the Riccati MSTMM [5] decreases the matrix size by half. Consequently, Riccati MSTMM keeps all the advantages of classical TMM and MSTMM, and simultaneously improves numerical stability. However, so far the RTMM is only applied to chain systems and cannot be used directly for the new version of MSTMM.
Therefore, both aspects will be addressed in this paper. The Riccati transformation is introduced on the basis of the new version of MSTMM for chain and branch systems. Compared with ordinary methods, the proposed method in this paper benefits from advantages such as no need of global system equations, low order of system matrices and high computational speed. Besides, the proposed method ensures that the overall transfer equation and transfer matrices of elements strictly satisfy system boundary conditions during the entire dynamics computation. The new version of MSTMM doesn't need linearization, and consequently the new version of RTMM belongs to the accurate methods regarding equation setup like the new version of MSTMM. Examples of small and huge branch systems moving in plane and space are given to demonstrate and validate the proposed approach.

State vectors and transfer equations
For comprehension some basics of MSTMM will be provided in the following. Although the proposed Riccati approach is totally general, later application for validation will involve only spatial bodies and ball joints, which is why only those transfer matrices are provided.

Topology graph and general concept
A multibody system typically consists of bodies and hinges as shown in Fig. 1(a). The corresponding topology in Fig. 1(b) provides the transfer relationship among state vectors, element numbers, and transfer directions [7], where body elements are represented by circles, with inside numbers indicating the bodies' sequence numbers, whereas arrows and numbers beside indicate hinge elements and hinge sequence numbers; the arrowhead shows the transfer direction of state vectors. Especially, the number 0 represents system boundaries. Body elements can be distinguished as single-input-single-output (SISO) or multiple-inputsingle-output (MISO) bodies, such as body 11 in Fig. 1(a).
At each connection point of two elements, the motion is described by position coordinates x, y, z and angles θ x , θ y , θ z w.r.t. an inertial frame. Internal cutting forces are denoted by q x , q y , q z whereas moments are denoted as m x , m y , m z . At an input point of an element, forces are considered positive in the axis direction, whereas moments about negative frame axes are considered as positive. At an output point this sign convention is inverted [4]. According to the new version of MSTMM [6], the state vector of these points is defined as The input end I is the origin of a body-fixed coordinate system I ξηζ used to describe the inertia tensor J I and angular velocity ω I = A T I Ω I , where A I is the direction cosine matrix transforming vectors from the body-fixed coordinate system I ξηζ into the inertial frame oxyz. All other quantities are described in the inertial frame like vectors r I C , r I O from input point to mass center C and output point O, respectively, external moments m C and external forces f C ; I 3 is the 3 × 3 identity matrix and O 3×3 , O 3×1 and O 1×3 are zero matrix and zero vectors, respectively. The tilde operator on a vector denotes its skew-symmetric matrix utilized for cross-products.
(2) Spatial rigid body with multi-input ends For a spatial rigid body with multiple input ends I 1 , I 2 , . . . , I L , like body 11 in Fig. 1(b), the state vector (1) of input ends has to be extended as Similar to Eq. (3), a transfer relation may be deduced from rigid body kinematics and kinetics resulting in an extended body transfer matrix where the additional columns account for forces and moments acting on inputs I 2 , . . . , I L . The abbreviations refer to Eq. (5) where I has to be substituted by I 1 as origin of body-fixed coordinate system I 1 ξηζ . Additionally, the skew-symmetric matrices E I j = −r I 1 I j (9) are used. Note that Eq. (8) reduces to transfer matrix (4) of a single-input element after removing the middle part.

(3) Smooth ball-and-socket hinge
For a smooth ball-and-socket hinge i, where its outboard element (i + 1) is a spatial rigid body, the corresponding position coordinates, and thus accelerations, and the internal forces are respectively identical, whereas the internal moments are zero. The relation for the angular acceleration may be deduced from the transfer relations (3)-(4) of its outboard body [6]. The transfer equation of a smooth ball-and-socket hinge is then identical to Eq. (3) with transfer matrix where u 3,1 = E 6 E 3 + E 7 , u 3,2 = E 6 E 4 + E 8 , u 3,3 = I 3 , u 3,4 = E 6 and u 3,5 = E 6 E 5 + E 9 refer to the particular submatrices in the third block row of the partitioned transfer matrix (4) of the outboard body, respectively.
(4) Multi-hinge subsystem In a branch system, all hinges (e.g., hinges 4 and 10 in Fig. 1(b)) with the same MISO outboard body i (body 11 in Fig. 1(b)) are treated as a hinge subsystem j with multiple input and multiple output ends, named multi-hinge subsystem. Its transfer equation includes two parts, the main transfer equation and additional geometric equations. The main transfer equation describes the relation among the accelerations and the internal moments and forces, while the geometric equations describe kinematic relations among the first input end I 1 and its output ends identified with the input state (6) of its outboard body i, i.e., z j,O ≡ z i,I 1 −I L . The state vector z j,I k of each input end I k (k = 1, . . . , L) is defined according to Eq. (1).

(a) Main transfer equation
Analogously to Eqs. (3) and (10), the main transfer equation describes the relationship among the comprehensive output state vector (6) and the state vectors (1) of all input ends, where just the additional moment (m j,I k ) and force (q j,I k ) impacts from I 2 , . . . , I L have to be accounted for. This can be achieved by the relation z j,O =Û j,I 1 z j,I 1 +Û j,I 2 z j,I 2 + · · · +Û j,I L z j,I L (11) with an enlarged matrix (10), given aŝ and additional matriceŝ where (k = 2, 3, . . . , L). Similar to Eq. (10), the elements u 3,m in matrices (12) and (13) refer to the third block row of the MISO outboard body transfer matrix (8).

(b) Geometric equation
We may notice in the main transfer equation that the quantity of unknowns is greater than the number of equations. Therefore, the geometric equations have to be supplemented by additional equations resulting from rigid body kinematics. Generally, the accelerations of an arbitrary point P of the rigid body are related to those of the input end I 1 by [4]r P = r I 1 −r I 1 PΩ I 1 +Ω I 1Ω I 1 r I 1 P . Since the output ends of a multi-hinge subsystem are identical with the input ends of its rigid outboard body, and since hinge positions of corresponding output and input ends are identical, these equations may be used to relate the accelerations of the hinge input ends to its first output end as Additionally, the moments at input (and output) ends of smooth hinges vanish, i.e., m i,I k = 0 (k = 2, . . . , L). Combined with Eq. (14) and by using the extended state (6), this may be summarized as where the matrices extract the respective accelerations and the moments from the state vectors of the hinge input/output ends and are the corresponding coefficients in Eq. (14).

Application of Riccati transformation to chain systems
The overall transfer equation of a chain system composed of n elements is [6] z n,0 = U 1−n z 1,0 where the overall transfer matrix is obtained by successive multiplication of element transfer matrices U i , i.e., The vectors z 1,0 and z n,0 represent tip and root states of the chain, respectively. As the chain length n increases, the accumulated errors caused by successive multiplication of transfer matrices U i may increase and may probably result in complete failure, see [9] for classical MSTMM. The Riccati transfer matrix method has shown to be effective for overcoming such kind of numerical problems [4]. In the following, a new version of Riccati MSTMM based on the new version of MSTMM is established, where the computing procedure is rather similar to linear or classical MSTMM. Usually, half of the state variables (2) in the boundary states of Eq. (18) are known while the other half is unknown. Let generally z a compose the known state variables while z b summarize the unknown ones, and terms associated with the artificial state "1" are separated. Then the transfer equation (3) of an element i can always be written as where T 11 , T 12 , T 21 , T 22 are the rearranged submatrices of transfer matrix U i associated with physical state quantities and f a and f b correspond to the last column of U i .
The assumption that z a is known from z b may be expressed by the Riccati transformation [4,5] based on input states. Substitution into Eq. (20) yields By rewriting it for the unknown input states, we obtain With abbreviations Eq. (23) may be rewritten as By substitution of Eqs.
This may be rearranged as which relates the known (z ai,O ) and unknown (z bi,O ) output states of element i. In a chain system, an element's output end is equal to the input end of its outboard element, i.e., z i,O ≡ z i+1,I . Due to this, Eq. (27) may be also written as similar to Eq. (21). From comparison, we may define resulting in recursion formulae for the new version of Riccati MSTMM, which according to Eq. (27) finally reads as For a chain system composed of n elements, the output point of the last element yields It should be noted that this recursion acts in the inverse transfer direction in contrast to, e.g., Eq. (3).

Algorithms using the Riccati transformation
Let us, e.g., consider the left subchain in Fig. 1(b). The input of body 1 is a boundary (denoted by 0) where half of the states are known to be zero (either accelerations for fixed node or forces/moments if free). This may be transferred into a boundary condition for S and e at the beginning. According to the definition that z a summarizes the know states, which are here z a1,I = 0, whereas z b1,I = 0 are the unknown boundary conditions, Eq. (21) may be satisfied by choosing In order to proceed, the transfer matrix of body 1 has to be rearranged according to Eq. (20). Then Eq. (24) may be applied to obtain P 1 and q 1 , which may be then substituted in (29) and (30) to obtain S 2 and e 2 . Next the hinge transfer matrix U 2 is rearranged and used in the same way to get S 3 and e 3 , and so on. The recursion ends with S n+1 and e n+1 which are the Riccati coefficients of the chain root (here output of body 3 where n = 3 in Eq. (32)). Assuming that this is already a system boundary, half of the state variables would be known and half unknown splitting the associated state vector in z an,O and z bn,O . By rearranging Eq. (32) accordingly, the unknown boundary states could be determined. By applying Eq. (33), the states z bi,I and z ai,I of any joint within the subchain may be obtained. Based on these considerations, the following general procedure is proposed for chain systems.

Procedure for chain systems
For a chain system (or subchain), the following Algorithm A may be applied: A chain system contains only a single tip while a branch system has at least two tips. Additionally, the overall transfer equation of a branch system not only includes the main transfer equation but also geometric equations (15). In order to derive the overall system equations, the branch system is split into (multiple) subchains and output chain, and a multihinge is introduced according to Fig. 1(c). By using Algorithm 4.1, the state vectors of all output ends of the subchains can be expressed in terms of input boundary conditions. These outputs may then be combined with the concept of a multi-hinge according to Sect. 2.2. Summarizing, the following procedure is proposed.

Applications
To explain the Riccati procedure step by step, a simple planar branch system according to Fig. 2(a) will be investigated in Sect. 5.1. In order to ease traceability, computations are performed symbolically which, however, may give a wrong impression. In real applications, the matrix operations in Sect. 4 are executed automatically on a numerical basis as demonstrated with a much larger second example in Sect. 5.2.

Demonstration example
The pendulum system in Fig. 2(a) consists of similar rods with mass m = 1 kg and length 2a = 1 m resulting in moment of inertia J z = 4 ma 2 /3 with respect to rod ends. This planar problem is a special case of spatial MSTMMs where the state vector (1) may be reduced to the essential coordinates z = ẍÿΩ z m z q x q y 1 T .
According to step (B1), the branch system is firstly divided into two subchains, a multihinge and multi-input body 7 as shown in Fig. 2(b). For treating the subchains in step (B2) along Algorithm A, respectively, the transfer matrices may be established in step (A1). For pendulum bodies i ∈ {1, 3, 5} (Fig. 2(c)), we may compute transfer matrix (4) from abbreviations (5) based on r I C = [ac i as i 0] T , r I O = 2r I C , Ω I = [0 0θ i ] T , f C = mg[0 − 1 0] T , m C = 0, and J I = diag{0 J z J z } where c i = cos θ i , s i = sin θ i . The resulting transfer matrix may then be reduced according to the reduced state (35) by using rows and columns 1, 2, 6, 9, 10, 11 and 13, respectively: The first subchain consists of body 1 only. Its lower end is a free boundary where forces and moments are zero, i.e., z 1 (37) A comparison with Eq. (20) yields where the body index has been omitted for simplification. Equation (29), together with S 1 = O and Eq. (24), then yields Substituting e 1 = 0 and Eq. (24) into Eq. (30) yields According to Eq. (31), the relationship between z aO,1 and z bO,1 can be written as Since the first subchain consists of a simple body only, a recursion along step (A3) is not necessary. Instead, the complete state of the first subchain tip may be written as where Now we may compute the second subchain along Algorithm A. Due to the similar boundary condition, body 3 may be treated analogously to body 1, resulting in the same matrices S 2 and e 2 as Eqs. (39) A similar order of state variables as in Eq. (37) yields Eq. (20) with submatrices Substitution into Eqs. (29) and (30) finally gives After another recursive step for body 5 with transfer matrix (36), where i = 5, we obtain S 4 and e 4 , which are too large to be shown in symbolic form. The tip state z 5,O may then be finally written analogously to Eq. (42) as The state vectors z 1,O and z 5,O are the two inputs of the multi-hinge 8 with MISO body 7 as outboard body which now have to be combined in the algorithmic step (B3). Let us start with the multi-input body 7, where, according to Fig. 2(d), r I O ≡ r I C = [ac 7 as 7 0] T and Ω I , f C , m C , and J I are identical to the analogous quantities used for computing Eq. (36). By adding intermediate columns regarding the second input given as r I 1 I 2 = 2r I C , we get from Eqs. (8) and (9) the reduced transfer matrix delivering the output state z 7,O = U 7,I 1 −I 2 z 7,I 1 −I 2 (50) according to Eq. (7), where z 7,I 1 −I 2 = ẍ I 1ÿ I 1Ω zI 1 m zI 1 q xI 1 q yI 1 m zI 2 q xI 2 q yI 2 1 T is the reduced form of Eq. (6). Its input is identical with the output of multi-hinge 8, resulting from Eq. (11) as The corresponding transfer matrices are given by Eqs. (12) and (13) By substituting Eq. (51) into Eq. (50), we obtain the main transfer equation Additionally, the geometric compliance equations (15) need to be considered, reading here in reduced form as a single equation where Eq. (51) According to algorithmic step (B3), Eqs. (53) and (54) By combining z a1,I = 0 with this result for z b1,I , the full state vector z 1,0 is obtained. With the same procedure, we can get also all state vectors in the second subchain, where the recursive formula (33) have to be applied three times in reverse transfer direction as specified in algorithmic step (B4). With the above formulas, a simulation may be performed for initial conditions θ 1 (0) = θ 3 (0) = θ 5 (0) = π/2, θ 7 (0) = 0,θ i (0) = 0 ∀i. In Fig. 3, the computational results (black line) are compared with those of ADAMS [10] (stars) showing a perfect match, which validates the proposed procedure. The numerical simulation was performed with a fourth-order Runge-Kutta integration scheme with constant step-size t = 0.001 s.

Large scale examples
The real strength of the new version of the Riccati MSTMM is the analysis of branch systems with long chains. In order to see the computational efficiency of the proposed method, numerical simulations of the branch system in Fig. 1(a) with a wide range of DOFs [11] were carried out and compared with ADAMS which uses Cartesian coordinates with constraints, resulting in differential-algebraic-equations. The parameters of the bodies are the same as in the demonstration example. Each body is connected with a 3DOF ball-and-socket hinge. The initial condition is that of Fig. 1(a) where body N hanged on the ceiling and the multi-input body 11 is horizontal, while others are vertical.
The first application example consists of 500 bodies associated with N = 999. Both methods are able to simulate the system. The simulation result of the rotation about the zaxis of two boundary bodies (1 and N ) and the multi-input body 11 is shown in Fig. 4 which agrees very well for both procedures. However, the time cost for the proposed method is by a factor of 36 lower than that with ADAMS as shown in Table 1. The CPU used is Xeon E5-2640 2.00 GHz. For higher number of bodies, ADAMS fails due to running out of computer

Conclusions
Based on the new version of MSTMM, the Riccati transformation is introduced to obtain a new version of Riccati MSTMM. This is achieved by developing a transfer equation of a multi-input subsystem for branch systems and establishing recursion formula for the transfer equation. Two algorithms are proposed, one for chain systems and one for branch systems, with clear rules how to use the recursion formula iteratively along chains and subchains. Compared with MSTMM and ordinary multibody system dynamics, the proposed method benefits from advantages regarding size of matrices and algorithmic stability. It is also ensured that the overall transfer equation satisfies strictly the boundary conditions. Application to a simple example demonstrates the easy use of the procedure, and several examples computed numerically for large branch systems demonstrate good agreement with classical strategies, validating the proposed method as well as its higher computational speed compared to classical multibody system formulations.