On the Computational Methods for Solving the Differential-Algebraic Equations of Motion of Multibody Systems

In this investigation, different computational methods for the analytical development and the computer implementation of the differential-algebraic dynamic equations of rigid multibody systems are examined. The analytical formulations considered in this paper are the Reference Point Coordinate Formulation based on Euler Parameters (RPCF-EP) and the Natural Absolute Coordinate Formulation (NACF). Moreover, the solution approaches of interest for this study are the Augmented Formulation (AF) and the Udwadia–Kalaba Equations (UKE). As shown in this paper, the combination of all the methodologies analyzed in this work leads to general, effective, and efficient multibody algorithms that can be readily implemented in a general-purpose computer code for analyzing the time evolution of mechanical systems constrained by kinematic joints. This study demonstrates that multibody algorithm based on the combination of the NACF with the UKE turned out to be the most effective and efficient computational method. The conclusions drawn in this paper are based on the numerical results obtained for a benchmark multibody system analyzed by means of dynamical simulations.


Introduction
In the last three decades, multibody system dynamics has emerged as an independent and interdisciplinary research field dedicated to the analysis and the synthesis of the motion of mechanical systems connected by kinematic pairs [1]. In the scientific literature, this particular class of dynamical system is referred to as multibody mechanical systems [2,3]. Multibody systems are mechanical systems defined by a collection of rigid and deformable continuum bodies, mechanical joints, force elements, and force fields [4][5][6][7][8][9][10]. The mathematical description of the time evolution of a multibody system is characterized by the presence of intrinsic nonlinearities that induce large reference displacements and large finite rotations [11][12][13]. Therefore, the resulting complex dynamic behavior of a mechanical system constrained by kinematic pairs can be described by a large set of nonlinear differential-algebraic dynamic equations [14,15]. In the field of multibody system dynamics, general analysis approaches are required in order to capture the dynamic behavior of a given multibody mechanical system. In general, the analytical techniques used for modelling multibody systems need to facilitate the formulation of the differential dynamic equations and lead to a consistent modelling of the mechanical joints mathematically represented by nonlinear algebraic equations [16][17][18][19][20]. The correct modelling of a multibody system is of paramount importance in several industrial applications such as, for example, in vehicle system dynamics, in aerospace engineering, and, more generally, in the problem of the engineering design of control actuators for mechanical systems formed by rigid and flexible components [21][22][23][24][25][26][27][28][29]. In particular, in the case of rigid-flexible multibody systems, advanced methods must be adopted for obtaining an estimation of the system state, for identifying the unknown or unmeasurable system parameters, and for controlling the system dynamical behavior [30][31][32][33][34][35][36][37][38][39][40][41][42].
As shown in several investigations, the multibody approach used for modelling the dynamics of mechanical systems constrained by mechanical joints represents an effective and efficient method for describing the kinematic structure of a given mechanical system and for analyzing the dynamic behavior resulting from a prescribed loading condition [43][44][45][46][47][48][49][50][51]. In this paper, therefore, a comparative study is carried out considering two general formulation strategies and two important solution procedures for solving the differential-algebraic dynamic equations of mechanical systems composed of multiple rigid bodies connected by mechanical joints. To this end, two modelling approaches based on Cartesian coordinates are considered, namely the Reference Point Coordinate Formulation based on Euler Parameters (RPCF-EP) and the Natural Absolute Coordinate Formulation (NACF) [52]. Unlike the RPCF-EP, the kinematic description of the NACF is based on the separation between variables that are space-dependent and the coordinates that are time-dependent [53]. The property of separation of variables used in the kinematic description of the NACF allows for formulating a system of equations of motion, which is characterized by a mass matrix that is independent from the generalized coordinate vector. As a result, the Coriolis and centrifugal generalized inertia terms do not appear in the mathematical derivation of the multibody dynamic equations [54]. In this investigation, on the other hand, two effective solution procedures for the calculation of the generalized acceleration vector of a multibody mechanical system are examined. The solution procedures considered in this paper are the Augmented Formulation (AF) and the Udwadia-Kalaba Equations (UKE) [55]. While the AF is a well-established multibody computational procedure, the UKE represents a new methodology recently discovered in the Lagrangian reformulation of classical mechanics. The UKE can be effectively employed for computing the generalized acceleration vector of mechanical systems constrained by holonomic and/or nonholonomic algebraic equations. Furthermore, the UKE represents one of the most general methods of classical mechanics for obtaining closed-form solutions of the fundamental problem of constrained dynamics. Therefore, the use of the UKE in the development of new, effective, and efficient computational algorithms is of interest for the multibody research community. In order to obtain a systematic comparison of the computational methodologies considered in this paper, a mechanical model of a benchmark multibody system is employed as a numerical example for performing numerical experiments.
The structure of this paper can be summarized as follows. In Section 2, the key equations of the two general multibody formulation approaches of interest for this investigation are described. In Section 3, the main features of the two general multibody solution procedures considered in this work are illustrated. In Section 4, the numerical results obtained implementing the methods analyzed in the paper in the case of a simple multibody benchmark problem are discussed. In Section 5, a discussion on the methodologies considered in this paper, the conclusions obtained in this investigation, the summary of the work, and some suggestions on future directions of research are provided.

Multibody Coordinate Formulations
The multibody formulation approaches considered in this section are the RPCF-EP and the NACF.

RPCF-EP
In the RPCF-EP, the generalized position of a body i is defined using a generalized coordinate vector given by: where R i is the global position vector of the origin of a body-fixed reference system and p i is a vector of rotational coordinates employed as orientation parameters. The orientation parameter vector p i By using the virtual work of the generalized external forces, the generalized external force vector Q i e of the body i can be calculated in the RPCF-EP as: where L i is the Jacobian matrix of the body generalized configuration and F i e denotes a general external force applied on an arbitrary point of the rigid body i. In the RPCF-EP, the Jacobian matrix L i is given by: By using a Lagrangian formulation approach, one can write the index-three form of the differential-algebraic dynamic equations of a rigid multibody system in the framework of the RPCF-EP as follows: where q is the vector containing the total set of generalized coordinates of the multibody system, M represents the system mass matrix resulting from the multibody assembly procedure, Q v denotes the system quadratic velocity vector associated with the inertia terms, Q e indicates the generalized external force vector that takes into account all the forces acting on the multibody system, v stands for the vector of all the Lagrange multipliers relative to the algebraic constraint equations, C is the complete vector of algebraic constraint equations, and C q represents the Jacobian matrix of the entire set of algebraic constraint equations. In the RPCF-EP, the complete vector of algebraic constraint equations is defined as: where Φ represents the complete vector of normalization constraints associated with the unit quaternions that serve as orientation parameters and Ψ is the vector of all the algebraic equations relative to the mechanical joints.

NACF
In the NACF, the generalized position of a body i is defined using a generalized coordinate vector given by: where R i is the global position vector of the origin of a body-fixed frame of reference and d i represents a vector of rotational coordinates employed as orientation parameters. The orientation parameter vector d i used in the NACF is formed by a set of nine rotational coordinates given by the direction cosines of the body-fixed coordinate system. The orientation vector d i can be mathematically defined as follows: where: where a i , b i , and c i are the three unit vectors based on the set of the direction cosines associated with the body reference system. In the NACF, the transformation matrix A i that describes finite rotations in the three-dimensional space can be explicitly defined as: The nine direction cosines a i 1 , and c i 3 are not independent rotational coordinates because they must form an orthonormal set of three unit vectors [56]. The normalization conditions of the direction cosines are given by: where Φ i represents a vector of algebraic constrains that contains the normalization equations of the rigid body i. In the NACF, the position vector of a generic point on the rigid body i referred to an inertial frame system can be expressed as: where: where the constant matrix S i is a rectangular matrix that defines the geometric configuration of the body i, I is a 3 × 3 identity matrix, whereasx i ,ȳ i , andz i are the Cartesian coordinates of the position vector of a generic point defined in the body-fixed coordinate frame. The mass matrix M i of a body i that appear in the equations of motion can be expressed in the framework of the NACF as follows: where I is a 3 × 3 identity matrix, m i is the total mass of the body i, are the second moments of mass of the body i. The first and second moments of mass can be obtained using the local position of the body center of mass of the body i denoted withū i G i and considering the body mass moments of inertia that appear in the rigid body inertia tensor. In the NACF, the mass matrix M i is constant, symmetric, and positive definite. Since in the NACF the mass matrix is constant, the inertia quadratic velocity vector Q i v that represents the centrifugal and Coriolis generalized inertia effects vanishes. On the other hand, in the NACF, the Jacobian matrix L i of the body i coincides with the matrix of shape functions denoted with S i . Thus, by using the virtual work of the generalized external forces, the generalized force vector Q i e acting on a rigid body i can be expressed in the framework of the NACF as: where F i e denotes a general external force applied on a generic point of the body i. The index-three form of the differential-algebraic equations of motion of a rigid multibody system can be derived considering a Lagrangian formulation approach and using the NACF multibody framework as follows: where e is the vector containing the total set of generalized coordinates of the multibody system, M represents the system mass matrix resulting from the multibody assembly procedure, Q e indicates the generalized external force vector that takes into account all the forces acting on the multibody system, v stands for the vector of all the Lagrange multipliers relative to the algebraic constraint equations, C is the complete vector of algebraic constraint equations, and C e represents the Jacobian matrix of the entire set of algebraic constraint equations. In the NACF, the complete vector of algebraic constraint equations is defined as: where Φ represents the complete vector of normalization constraints associated with the direction cosines and Ψ is the vector of all the algebraic equations relative to the mechanical joints.

Multibody Solution Methods
The multibody solution procedures considered in this section are the AF and the UKE. These two multibody solution procedures can be equally applied to the dynamic equations mathematically derived by using the RPCF-EP as well as the NACF. However, for simplicity, the index-one differential-algebraic dynamic equations devived employing the RPCF-EP are considered in this section for describing the multibody methods based on the AF and the UKE.

AF
In this subsection, the multibody solution method based on AF is discussed. To this end, consider the following system of index-one differential-algebraic dynamic equations derived employing the RPCF-EP: where Q d is a vector called constraint quadratic velocity vector that includes the terms that are quadratic in the generalized velocities. In the AF, the index-one equations of motion can be reformulated in a matrix form as follows: This matrix equation formulated using the AF can be written in a compact symbolic form as: where q a is the multibody system augmented generalized acceleration vector, M a is the multibody system augmented mass matrix, and Q a is the multibody system augmented generalized force vector that are respectively defined as: The linear system of algebraic equations formulated by using the AF can be readily solved by implementing any method for the numerical solution of a system of linear equations. By doing so, one can easily obtain the system generalized acceleration vectorq and, at the same time, the Lagrange multiplier vector v. The system generalized acceleration vectorq can be used in a standard numerical integration scheme in order to calculate the numerical solution for the dynamic state of the multibody system. The vector of Lagrange multipliers v, on the other hand, can be used for calculating the generalized force vector that mathematically models the reaction force vector of the mechanical joints.

UKE
In this subsection, the multibody solution approach based on UKE is illustrated. The UKE, also known as fundamental equations of constrained dynamics, were discovered by Udwadia and Kalaba during their important research in the field of analytical dynamics. Udwadia and Kalaba focused their research on modern linear algebra methods as well as on the foundations of classical mechanics such as the principle of least constraint originally formulated by Gauss [57]. In the general form of the UKE, an auxiliary matrixM and an auxiliary generalized force vectorQ are respectively defined as: The multibody solution procedure based on the UKE leads to closed-form analytical solutions of the generalized acceleration vectorq. For this purpose, the UKE are defined as follows: whereā represents the system acceleration vector obtained when the algebraic constrains are absent, e c denotes the error vector associated with the algebraic equations,K is referred to as the system kinetic matrix,F identifies the feedback matrix relative to the algebraic constrains generated by the action of the mechanical joints,ā c is the additional acceleration vector caused by the presence of the algebraic constraints, andq is the complete system generalized acceleration vector. In the UKE, the matrix denoted withK + represents the pseudoinverse matrix of the multibody system matrixK called kinetic matrix [58,59]. By using the UKE, one can readily find the generalized force vector Q c relative to the entire set of algebraic constraints that limit the motion of the multibody system as well as the vector of Lagrange multipliers v useful for quantifying the generalized reaction forces of the kinematic pairs. Furthermore, the generalized acceleration vectorq of the multibody system necessary for performing the progressive marching of the numerical solution of the differential-algebraic dynamic equations on the time grid can be easily calculated in a closed-form employing the approach based on the UKE.

Numerical Results and Discussion
In this section, a numerical analysis is carried out in order to assess the performance and the reliability of the formulation approaches and solution methods discussed in the paper. To this end, a multibody computer code developed by the authors and programmed in MATLAB (R2013a version) is used for obtaining the numerical results discussed in this section. In Figure 1, a schematic representation of the multibody approach followed in this numerical study is shown.
In particular, the triple pendulum system represented in Figure 2 is considered as an illustrative example of a simple multibody mechanical system that undergoes a complex dynamic evolution.
Excluding the ground, the triple pendulum system represented in Figure 2 is composed of three rigid bodies and three spherical joints. As shown in Figure 2, the spherical joint collocated at the point A connects the body number 1 to the ground, while the spherical joints collocated at the points B and C serve as connections between the bodies 1, 2, and 3. The triple pendulum system is formed by three pendulums having the same geometric dimensions, namely half length L = 2.0 (m), breadth H = 0.2 (m), and width W = 0.2 (m). The additional numerical data used for modelling the triple pendulum system are reported in Table 1.  The configuration of the triple pendulum system at the initial stage of the dynamical simulation is shown in Figure 2 and the initial velocities of all the rigid bodies forming the multibody system are set equal to zero. The triple pendulum system is loaded with its own weight, which is due to a uniform gravity force field. The multibody dynamic equations of the triple pendulum system form a system of differential-algebraic dynamic equations that are analytically derived by using the RPCF-EP and the NACF. Subsequently, the Robust Generalized Coordinate Partitioning (RGCP) method is used for stabilizing constraint drift of the triple pendulum system as well as to enforce the kinematic constraints at both the position and velocity levels [60,61]. The constraint tolerance used in the Newton-Raphson (NR) numerical procedure implemented in the RGCP algorithm is equal to ε = 10 −9 . Furthermore, both the AF and the UKE are alternatively used for solving the equations of motion of the triple pendulum system at the acceleration level, whereas the fourth-order Adams-Bashforth (AB) method is used for the time marching of the numerical solution of the equations of motion. To this end, a constant time step equal to ∆t = 10 −3 (s) is employed for performing the dynamic analysis and a time interval equal to T = 20 (s) is considered for the dynamical simulations. Figures 3-5 respectively show the longitudinal, lateral, and vertical displacements and velocities of the point D at the tip of the triple pendulum system.     The numerical solutions represented in Figures 3-5 refer to the dynamic equations of the triple pendulum system modelled employing the RPCF-EP and the NACF as the multibody formulation approach, whereas the system generalized acceleration vector is obtained by using the UKE as the multibody solution algorithm. The numerical solutions obtained in the comparative study shown in Figures 3-5 exhibit a very good agreement and very similar numerical results are obtained employing the AF as a multibody solution algorithm. The consistency of the numerical solutions found in this work can be also observed from the trajectory of the tip of the triple pendulum system represented in Figure 6, which exhibits a chaotic motion.   Although the multibody model of a triple pendulum system is mathematically simple to derive, it is well known that this class of dynamical systems exhibits a complex physical behavior due to the nonlinearity of the equations of motion that lead to the chaos phenomenon. From a computational point of view, this property of the triple pendulum system implies that, if a small change in the initial conditions or a perturbation is induced to the numerical solution of the equations of motion because of the numerical approximations, a large difference in the resulting trajectory will be apparent in the subsequent dynamical evolution. Therefore, the simple multibody model of the triple pendulum system can be used in order to compare the numerical solution obtained employing the NACF and the RPCF-EP with the numerical solutions derived by using the AF as well as the UKE. Since there is a good agreement between the numerical solutions computed using the combination of these four different approaches, one can assume that the equations of motion are correctly solved. Therefore, the triple pendulum system can effectively serve as a multibody benchmark example in order to test the accuracy and the performance of new multibody computational algorithms by means of numerical experiments. For this purpose, in order to evaluate the accuracy of the numerical solutions obtained in this study, the norms of the constraint violations are computed for all the combinations of the methods considered in this investigation. In particular, the residuals of the constraint equations are evaluated at the position level as well as at the velocity level. To this end, Table 2 shows the maximum norms of the constraint violations for the constraint equations and their time derivatives for all the multibody methodologies analyzed in the paper. The numerical results provided in Table 2 demonstrate that all the multibody formulation approaches and solution strategies analyzed in this study are effective since they lead to a set of numerical results that are physically accurate and numerically robust. Furthermore, the time evolution obtained by means of dynamical simulations for the multibody system considered as a benchmark example is consistent with the geometric shape of the triple pendulum system. Moreover, Table 3 shows the dimensionless computational times of all the dynamic simulations. The numerical results reported in Table 3 demonstrate that, while the AF and the UKE show similar performance in terms of the computational times, the equations of motion formulated and solved employing the new multibody method based on the NACF lead to more efficient dynamical simulations when compared with the well-established multibody algorithm that relies on the RPCF-EP. This behavior can be explained by noticing the fact that the dynamic equations formulated by using the NACF have a constant mass matrix and the corresponding generalized inertia vector that contains the terms that are quadratic in the generalized velocities is always a zero vector.

Conclusions
The main topics of interest for the research of the authors are multibody dynamics, system identification, and nonlinear control [62][63][64][65][66][67][68][69][70]. This investigation represents a comparative analysis of the principal computational methodologies suitable for the analytical derivation and the numerical implementation of the dynamic equations of mechanical systems composed of rigid bodies constrained by kinematic pairs. The coordinate formulations considered in this paper are the RPCF-EP and the NACF, whereas the solution methods considered in this work are the AF and the UKE. The RPCF-EP is a well-established coordinate formulation that is used in commercial and research multibody computer codes for analyzing the dynamic behavior of mechanical systems constrained by kinematic joints. The NACF, on the other hand, is a new method recently developed by the authors by combining the key features of the well-known Reference Point Coordinate Formulation with Euler Angles (RPCF-EA) and the conventional Natural Coordinate Formulation (NCF). Unlike both the RPCF-EA and the RPCF-EP, the NACF leads to a system of differential-algebraic equations of motion in which the mass matrix is constant, while the centrifugal and Coriolis generalized inertia effects are identically equal to zero. More importantly, a straightforward formulation of the nonlinear equations that model the kinematic joints by means of algebraic constraints can be systematically obtained by using the multibody framework based on the NACF. Furthermore, the AF is a well-known multibody algorithm used for computing the generalized acceleration vector and the vector of Lagrange multipliers for a mechanical system formed by rigid bodies and mechanical joints. On the other hand, the UKE represents a new set of equations of analytical dynamics discovered by Udwadia and Kalaba. The UKE have several applications that span far beyond their original interpretation. In particular, the UKE can handle a general form of the algebraic constraint equations of holonomic and/or nonholonomic nature and can be effectively used as an efficient computational method for the numerical solution of multibody system problems. In the paper, a comparative analysis is carried out considering the combination of all the analytical formulation strategies and numerical solution approaches mentioned before. The numerical results arising from the numerical analysis carried out in this study showed that the use of the NACF as a formulation approach in conjunction with the UKE as a solution procedure represents a general, effective, and efficient computational algorithm suitable for analyzing the dynamic behavior of complex mechanical systems connected by kinematic constraints. This work represents a preliminary investigation oriented towards the future development of a more detailed comparative analysis. To this end, an array of experimental benchmark examples will be used in future investigations in order to compare the effectiveness and the efficiency of the approach proposed in this paper with the computational methods already available in the literature.
Author Contributions: This research paper was principally developed by the first author (Carmine Maria Pappalardo). The detailed review carried out by the second author (Domenico Guida) considerably improved the quality of the work.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.