Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 54, 4, pp. 1391-1404, Warsaw 2016 DOI: 10.15632/jtam-pl.54.4.1391 COMPARISON OF NATURAL COMPLEMENT FORMULATIONS FOR MULTIBODY DYNAMICS Marcin Pękal, Janusz Frączek Warsaw University of Technology, Institute of Aeronautics and Applied Mechanics, Poland e-mail: mpekal@meil.pw.edu.pl; jfraczek@meil.pw.edu.pl The main aim of this paper is to compare the effectiveness of numerical techniques used for spatial multibody dynamics simulation by applying the natural complement method. In the paper, seven numerical schemes are considered: zero eigenvalue formulation, Pseudo Upper TriangularDecomposition, Schur decomposition, Singular Value Decomposition, QR decomposition, coordinate partitioning andWang-Huston formulation. In order to illustrate the effectiveness of these schemes, two McPherson struts are considered. Simulations are performedwith four error tolerancevalues and for three stabilizationcases. Some suggestions on possible applications of the selected methods are formulated. Keywords: multibody dynamics, numerical methods, spatial systems 1. Introduction The equations of motion of multibody systemsmay be written usingmany different coordinate sets (de Jalón and Bayo, 1994; Malczyk and Frączek, 2012). However, the most frequently used are absolute, natural and relative coordinates. Absolute and natural coordinates constitute redundant sets of coordinates hence a system of equations of motion with constraints has to be considered. Such a system is the set of differential-algebraic equations (DAEs) with differential index at most equal to 3. Numerical methods used to solve DAEs are still under intensive development and the effective integration procedures are less well-known for them than for ordinary differential equations (ODEs) (Kunkel and Mehrman, 2006). Note that there exist publications about differential equations in mechanical systems, e.g. (Awrejcewicz, 2014; Eich- Soellner and Führer, 1998). Many integrationmethods canbeused to solveDAEs (de Jalón andBayo, 1994;Haug, 1989). A large group of thesemethods transforms the equations ofmotion described in redundant coor- dinates into equationswritten in theminimal set of coordinates. To performthis transformation, the null space base of the constraints manifold is built and afterwards, equations of motion are expressed in this base. Null space base vectors are represented by the orthogonal complementary matrix. Finally, the set of ODEs is obtained, which can be solved efficiently. There are several methods proposed in literature to accomplish the orthogonal complement. The main objective of this paper is to compare numerical schemes for DAEs based on or equivalent to the orthogonal complement methods using two examples. The following methods are considered: zero eigenvalue formulation, PseudoUpper Triangular Decomposition (PUTD), Schur decomposition, SingularValueDecomposition (SVD),QRdecomposition, coordinate par- titioning and Wang-Huston formulation. Hereafter, spatial rigid multibody systems described with the absolute coordinates are considered. TheEuler angles (consecutive rotations about the z, x and z axes) are used for the orientation description and consequently, the mass matrix is nonsingular. For the chosen coordinates set, the resulting equations of motion are in the form of differential-algebraic equations that consist of differential equations of motion and algebraic 1392 M. Pękal, J. Frączek equations of constraints. Note that there are also papers analyzing impact, e.g. (Awrejcewicz et al., 2003, 2004; Awrejcewicz andKudra, 2005) or contact, e.g. (Awrejcewicz andKudra, 2014), but these phenomena are not examined here. In order to compare the efficiency of the examinedmethods, two types of McPherson struts are analyzed. The first model of the strut has no redundant constrains, while in the second mechanism the redundancy is taken into consideration (Wojtyra and Frączek, 2013). For all the considered methods, simulations are performed with four values of the error tolerances and for three constraint stabilization approaches (with no stabilization of the constraints, with the Baumgarte stabilization and with the position constraints stabilization using the Newton- -Raphsonmethod). It is worth noting that similar comparisonswere already performedby other authors (Mariti et al., 2011, 2010; Pennestr̀ı and Valentini, 2007) as well as in (Pękal, 2012). Therefore, the comparison with the previous publications (Mariti et al., 2011, 2010) is also performed. 2. Spatial system dynamics Multibody dynamics description presented in this Section is based on the absolute coordinate formulation and appears in, e.g. (Frączek and Wojtyra, 2008; de Jalón and Bayo, 1994; Haug, 1989). The vector of the absolute generalized coordinates can be written as (Haug, 1989): qn×1 = [q T 1 q T 2 · · · q T N] T = [q1 q2 · · · qn] T, where N is the number of bodies and n is the number of generalized coordinates of the multibody system. Moreover, the equation of motion has the following form (de Jalón and Bayo, 1994; Haug, 1989) Mq̈+ΦTqλ=Q e (2.1) where M is the mass matrix, Φq is the Jacobian matrix of constrains (see Eq. (2.3)), λ is the Lagrange multipliers vector and Qe is the generalized force vector. Matrix equation (2.1) is the system of n equations with n+m unknowns (q̈n×1 and λm×1), where m is the number of Lagrange multipliers. To solve this system, the additional m constraint equations should be introduced as (de Jalón and Bayo, 1994; Haug, 1989) Φm×1 =Φ(q, t)=0 (2.2) where t denotes time. After differentiating (2.2) over time, we obtain (Haug, 1989) Φ̇=Φqq̇+Φt =0 (2.3) Differentiation of Eq. (2.3) once again over time yields (Haug, 1989) Φ̈=Φqq̈+(Φqq̇)qq̇+2Φqtq̇+Φtt =Φqq̈−Γ=0 ⇒ Φqq̈=Γ (2.4) Eventually, Eqs. (2.1) and (2.4) can be written as index-1 DAEs (de Jalón and Bayo, 1994; Haug, 1989) [ M ΦTq Φq 0 ][ q̈ λ ] = [ Qe Γ ] (2.5) where the coefficient matrix is called the augmented matrix (Negrut et al., 1997; de Jalón and Gutiérrez-López, 2013). It should be pointed out that Eqs. (2.2), (2.3) and (2.4) are analytically equivalent. However, the direct numerical solution of Eq. (2.5) does not often provide the fulfil- ment of position (2.2) and velocity (2.3) constraints andmay cause the solution drift. Therefore, Comparison of natural complement formulations for multibody dynamics 1393 the stabilizationmethods are often employed.One of the simplest is theBaumgarte stabilization method (Baumgarte, 1972). In this method, theΓ vector from acceleration constraints is repla- ced by the following expression:Γ=Γ−2α̂Φ̇−β̂2Φ, where α̂, β̂ are theBaumgarte stabilization parameters. The value of these parameters is often assumed as α̂ = β̂ ∈ 〈1,20〉 (de Jalón and Bayo, 1994). It is also possible to stabilize the system using the well-known Newton-Raphson method (Frączek and Wojtyra, 2008; Haug, 1989). The q vector is corrected when constraints norm (2.2) is greater than the assumed tolerance. 3. Orthogonal complement methods The equations ofmotion of themultibody systempresented previously constitute a set ofDAEs. In the present Section, those equations are transformed into ODEs bymeans of the orthogonal complement. To perform the transformation, the projection matrix P that is orthogonal to the constraint Jacobian matrix must be computed. The application of thePmatrix leads to direct removal of the Lagrange multipliers from the considered system (de Jalón and Bayo, 1994). 3.1. Orthogonal complement for rheonomic constraints The following statements are derived for systems without redundant constraints and then generalized to the case of redundant multibody systems. The procedure for finding orthogonal complementmatrices ispresented in (deJalónandBayo, 1994;Mariti et al., 2011, 2010;Pennestr̀ı and Valentini, 2007; Pękal, 2012; de Jalón and Gutiérrez-López, 2013). Assume that the independent velocity vector v̇ can be obtained by the projection of the generalized velocity vector into the rows of the constant matrix B (de Jalón and Bayo, 1994; Pennestr̀ı and Valentini, 2007) v̇s×1 =Bs×nq̇n×1 (3.1) where s is the number of independent coordinates. Combining Eqs. (3.1) and (2.3) yields (de Jalón and Bayo, 1994; Pennestr̀ı and Valentini, 2007) [ Φq B ] q̇=Xq̇= [ −Φt v̇ ] (3.2) If X is nonsingular, the inversion of this matrix exists and can be presented in the form: X−1 = [S P], where P is the projection matrix and S is the matrix which contains the re- maining columns. Hence (de Jalón and Bayo, 1994; Pennestr̀ı and Valentini, 2007) q̇=X−1 [ −Φt v̇ ] = [ Sn×m Pn×s ] [ −Φt v̇ ] =−SΦt+Pv̇ (3.3) Orthogonality of the X matrix gives two orthogonality conditions in the form (de Jalón and Gutiérrez-López, 2013; Pennestr̀ı and Valentini, 2007) (Φq)m×nPn×s =0m×s Bs×nPn×s = Is×s (3.4) The projection matrix P can be used to eliminate the Lagrange multipliers. Differentiating Eq. (3.2) and using Eq. (2.4) yields (de Jalón and Bayo, 1994; Pennestr̀ı and Valentini, 2007) Xq̈= [ Γ v̈ ] ⇒ q̈=X−1 [ Γ v̈ ] = [ S P ] [ Γ v̈ ] =SΓ+Pv̈ (3.5) 1394 M. Pękal, J. Frączek Substituting Eq. (3.5) into Eq. (2.1) gives M(SΓ+Pv̈)+ΦTqλ=Q e (3.6) In order to eliminate the Lagrange multipliers, left-multiplication of both sides of Eq. (3.6) byPT is required. Afterwards, the use of transposition of Eq. (3.4)1 yields (de Jalón and Bayo, 1994; de Jalón and Gutiérrez-López, 2013) P T M(SΓ+Pv̈)+0λ=PTQe ⇔ PTMPv̈=PTQe−PTMSΓ (3.7) Assuming v̈s×1 =0, SΓ can be obtained fromEq. (3.5) as (de Jalón and Bayo, 1994; Pennestr̀ı and Valentini, 2007) SΓ=X−1 [ Γ 0 ] (3.8) The vector v̈ is obtained from Eq. (3.7). Substituting this vector into Eq. (3.5) yields q̈. In the case of redundant systems, it is convenient to use the pseudoinversematrixX+ instead ofX−1. Therefore, the problembecomes an optimization task, whose result has the least square norm. 3.2. Derivation of the projection matrix Methods that directly apply the projection matrixP are introduced in the following. 3.2.1. Zero eigenvalue method The first method uses the zero eigenvalue technique (Mariti et al., 2011, 2010; Pennestr̀ı and Valentini, 2007; Pękal, 2012; Walton and Steeves, 1969) which is based on the eigenvalue problem. The eigenvalue problem can be written in the form (FreeMat v4.1; Hartfiel, 2001) AΨ=ΨΛ (3.9) where thematrixA is a symmetrical matrix,Λn×n = diag(Λ1,Λ2, . . . ,Λn) contains eigenvalues andΨn×n = [Ψ1 Ψ2 · · · Ψn] is the orthogonalmodalmatrix which contains the eigenvectors. In order to determine the projection matrix P, the symmetric matrix A is considered as (Walton and Steeves, 1969; Pennestr̀ı and Valentini, 2007) An×n =Φ T qΦq (3.10) Using this matrix, the following expression is obtained (Walton and Steeves, 1969; Pennestr̀ı and Valentini, 2007) Ψ T n×n(Φ T q)n×m(Φq)m×nΨn×n =Λn×n (3.11) where Λ is unique and contains non-negative eigenvalues. Assume that Λ1 ¬ Λ2 ¬ ··· ¬ Λn, then the possible zero eigenvalues have the lowest indices. Moreover, the number of positive eigenvalues is equal to the rank r of the Jacobian matrix. Therefore, s = n− r eigenvalues are equal to zero. These eigenvalues correspond to rigid bodymotion, so their number is equal to the number of degrees of freedom (DOFs) of the system. Consequently, the system has s independent and r=n−s dependent coordinates. Assuming (Walton and Steeves, 1969) (Φq)m×nΨn×n =Dm×n (3.12) Comparison of natural complement formulations for multibody dynamics 1395 Equation (3.11) can be written as (Walton and Steeves, 1969) D T D=Λ (3.13) The order of the eigenvalues is opposite to the sequence form (Walton and Steeves, 1969), thus Dm×n = [ 0m×s D̄m×(n−s) ] (3.14) hence (Pennestr̀ı and Valentini, 2007) Φqm×nΨ 1 n×s =0m×s (3.15) where Ψ1n×s is the matrix created from the first s columns of the modal matrix, i.e. columns which correspond to the zero eigenvalues. Note that equation (3.15) satisfies first orthogonality condition (3.4)1, hence: P = Ψ 1 n×s. Moreover, the modal matrix is orthogonal so condition (3.4)2 is fulfilled for:B=P T. 3.2.2. Schur decomposition method The description of the Schur decompositionmethod can be found in (Golub andLoan, 1996; Hartfiel, 2001; Kincaid and Cheney, 2002; Mariti et al., 2011, 2010; Pennestr̀ı and Valentini, 2007; Pękal, 2012). This decomposition is based on the fact that every square matrixA can be presented in the form (Golub and Loan, 1996) U H AU=T=Λ+N (3.16) whereU is aunitarymatrix,UH denotes conjugate transpositionofUandT is ablock triangular matrix, andN is strictly the upper triangular matrix. Note that the orthogonal matrix is a particular case of the unitarymatrix for real numbers, thus (Golub and Loan, 1996) U T AU=Λ+N (3.17) Assuming theA according to Eq. (3.10) and using the fact that theA is normal, it is possible to write (Golub and Loan, 1996) U T AU=Λ (3.18) Note that this is analogous to (3.11). Hence, it is possible to follow as for the zero eigenvalue method described in Section 3.2.1. Eventually, the following condition is get (Pennestr̀ı and Valentini, 2007) (Φq)m×nU 1 n×s =0m×s (3.19) where U1 is the submatrix of U (selected analogously to Ψ1). Therefore, first orthogonality condition (3.4)1 is fulfilled when: P = U 1 n×s, while second orthogonality condition (3.4)2 is satisfied for:B=PT. Since the Schur decompositionmethod is equivalent to the zero eigenvalue method, these approaches are treated as one in the following. 1396 M. Pękal, J. Frączek 3.2.3. PUTD method ThePseudoUpperTriangularDecompositionmethod has two variants and their description can be found in (Amirouche, 2006; Ider and Amirouche, 1988; Mariti et al., 2011, 2010; Ostal- lczyk, 2008; Pennestr̀ı and Valentini, 2007; Pękal, 2012). The first type uses the Householder transformation and the second uses the Gauss elimination. The transformation matrix H is defined first. Using this matrix, the ΦTq matrix can be transformed into the upper trapezoidal formZ (in particular into the upper triangular form) as (Amirouche, 2006; Ider and Amirouche, 1988) H T n×n(Φ T q)n×m =Zn×m (3.20) whereHmatrix can be obtained using theHouseholder transformation (themethod denoted as PUTD-H) or the Gauss elimination (the method denoted as PUTD-G). The matrix orthogonal to the Z is found in the following step. This matrix can be com- puted using the Gram-Schmidt orthogonalization of the upper-triangularized Jacobian matrix, which gives an identity matrixD (Amirouche, 2006). After this, the submatrix ofD is taken as (Pennestr̀ı and Valentini, 2007) D2(n−r)×n = [ 0(n−r)×r I(n−r)×(n−r) ] (3.21) thus (Pennestr̀ı and Valentini, 2007) D2(n−r)×nH T n×n(Φ T q)n×m =0(n−r)×m (3.22) which after transposing gives (D2H T Φ T q) T =Φq(D2H T)T =Φq(HD T 2 )=0m×(n−r) (3.23) Comparing Eqs. (3.23) and (3.4)1 yields: P = HD T 2 , while (3.4)2 is fulfilled when: B = P T, what is consistent with Pennestr̀ı and Valentini (2007). Considering the redundant systems, Eq. (3.21) is fulfilled when the redundant constraints are at the end of vector (2.2). It is not the general case. This problem can be solved by the Gauss-Jordan elimination with partial pivoting of the transposed Jacobian matrix. Using this elimination, the position of the independent constraints in the Jacobian matrix is obtained. Thus, it is possible to move dependent constraints at the end of the constraints set. Moreover, the problem of redundant constrains can be also solved by themanual setting of the constraints in Eq. (2.2) during the preprocessing stage. However, this simple approach can be used only for small systems. In the implementation, themethod based on theGauss-Jordan elimination is used to detect redundant constraints. The choice of the independent constraints is done only once at the be- ginning of simulation. It should be noted that in some cases the coordinate partition must be donemore frequently due to the loss of independence of the chosen coordinates. 3.2.4. SVD method Themethod using singular value decomposition is described in (Mani et al., 1985; Mariti et al., 2011, 2010; Pennestr̀ı andValentini, 2007; Pękal, 2012). It uses the decompositionwhich can be written in the following form (Pennestr̀ı and Valentini, 2007) (ΦTq)n×m =Wn×nLn×mV T m×m (3.24) where the matrices W andV are orthogonal, and L includes singular values of the transposed Jacobianmatrix on its diagonal. Assume thatm is the total number of the eigenvalues of theΦTq Comparison of natural complement formulations for multibody dynamics 1397 and the rank of the Jacobianmatrix r is equal to the number of nonzero eigenvalues. Thus, the decomposition can be written as (Pennestr̀ı and Valentini, 2007) (ΦTq)n×m = [ (Wd)n×r (Wi)n×(n−r) ] [ Λr×m 0(n−r)×m ] V T m×m =WdΛV T (3.25) Left-multiplication of Eq. (3.25) byWTi yields (Pennestr̀ı and Valentini, 2007) W T i(n−r)×n (ΦTq)n×m =W T i WdΛV T (3.26) Using orthogonality of the matrixW gives (Pennestr̀ı and Valentini, 2007) W T i Wd =0 W T i Wi = I (3.27) hence (Pennestr̀ı and Valentini, 2007) W T i(n−r)×n (ΦTq)n×m =0(n−r)×m ⇒ ΦqWi =0m×(n−r) (3.28) Eventually, first orthogonality condition (3.4)1 is met for (Pennestr̀ı and Valentini, 2007): P=Wi and second condition (3.4)2 is fulfilled when (Pennestr̀ı and Valentini, 2007): B=P T. 3.2.5. QR decomposition method TheQRmethod is described in, e.g. (Kim andVanderploeg, 1986; Mariti et al., 2011, 2010; Pennestr̀ı and Valentini, 2007; Pękal, 2012). It is based on the QR decomposition described in details by Golub and Loan (1996). This decomposition takes the form (Kim and Vanderploeg, 1986) Φ T q =QR (3.29) where thematrixQ is orthogonal andR is upper trapezoidal orupper triangular in theparticular case of the square decomposedmatrix. ThematrixR can bewritten as (Pennestr̀ı andValentini, 2007) Rn×m = [ (R1)r×m 0(n−r)×m ] (3.30) However, when the redundant systems are considered, this form of the R matrix occurs only when the redundant constraints are placed at the end of constraints vector (2.2). There are several methods to achieve this, e.g. appropriate definition of the constraints vector, use of the Gauss-Jordan elimination with partial pivoting or application of the alternative version of the QR decomposition which can be written in the form (FreeMat v4.1; MATLABr) Φ T qE=QR (3.31) whereE is a permutation matrix, which allows the matrixR to take the form from Eq. (3.30). Eventually, the QR decomposition of the transposed Jacobian matrix yields (ΦTq)n×mEm×m =Qn×nRn×m = [ Q1n×r Q2n×(n−r) ] [ R1r×m 0(n−r)×m ] =Q1R1 (3.32) Left-multiplying Eq. (3.32) by theQT2 gives Q T 2Φ T qE=Q T 2Q1R1 (3.33) 1398 M. Pękal, J. Frączek Using the orthogonality property of theQmatrix yields (Pennestr̀ı and Valentini, 2007) Q T 2Q1 =0 Q T 2Q2 = I (3.34) Transposition of Eq. (3.33) using Eq. (3.34) gives (QT2Φ T qE) T =ETΦqQ2 =0 (3.35) Left-multiplying Eq. (3.35) by the permutation matrix yields ΦqQ2 =0 (3.36) Thus, first orthogonality condition (3.4)1 is fulfilled when: P = Q2 and second orthogonality condition (3.4)2 is met for:B=P T. 3.3. Related methods The following Section describes methods that do not use directly the projection matrix P. Note that thesemethods are equivalent to the orthogonal complementmethods. Itwas described by, e.g. de Jalón andBayo (1994) in the case of the coordinate partitioningmethod andbyWang and Huston (1989) for theWang-Huston formulation. 3.3.1. Coordinate partitioning method The coordinate partitioning method is described in (de Jalón and Bayo, 1994; Mariti et al., 2011, 2010; Pennestr̀ı and Valentini, 2007; Pękal, 2012; Wehage and Haug, 1982). The coordinates are partitioned into independent v and dependent u sets using, e.g. the Gauss-Jordan elimination with partial pivoting of the Jacobian matrix. Thus, according to the coordinate partitioning, the differential-algebraic equations of motion from Eq. (2.5) take the following form (these are dependences fromPennestr̀ı andValentini (2007), written in thematrix form)    Muu Muv ΦTu Mvu Mvv ΦTv Φu Φv 0       ü v̈ λ    =    Qeu Qev Γ    (3.37) giving    ü v̈ λ    =    Φ−1u (Γ−ΦvM̆ −1Q̆e) M̆−1Q̆e (Φ−1u ) T(Qeu−MuvM̆−1Q̆e−MuuΦ−1u (Γ−ΦvM̆ −1Q̆e))    (3.38) where (Pennestr̀ı and Valentini, 2007) M̆=Mvv−MvuΦ−1u Φv−Φ T v(Φ −1 u ) T(Muv−MuuΦ−1u Φv) Q̆ e =Qev−MvuΦ−1u Γ−Φ T v(Φ −1 u ) T(Qeu−MuuΦ−1u Γ) (3.39) Note that the Φ−1u should be replaced by its pseudoinverse Φ + u for redundant systems. Moreover, all the generalized coordinates are integrated during the simulations in order to obtain comparable results to the outcomes from the other methods. Comparison of natural complement formulations for multibody dynamics 1399 3.3.2. Wang-Huston formulation Themethod, described byMariti et al. (2011, 2010), Pękal (2012),Wang andHuston (1989) is based on pseudoinversion and is equivalent to the orthogonal complementmethod (Wang and Huston, 1989). Using Eq. (2.1) yields (Mariti et al., 2010; Wang and Huston, 1989) q̈=M−1(Qe−ΦTqλ) (3.40) Substituting Eq. (3.40) into Eq. (2.4) gives (Wang and Huston, 1989) λ=(ΦqM −1 Φ T q) −1(ΦqM −1 Q e−Γ) (3.41) and substituting Eq. (3.41) into Eq. (3.40) leads to (Wang and Huston, 1989) q̈=M−1ΦTq(ΦqM −1 Φ T q) −1 Γ− (M−1ΦTq(ΦqM −1 Φ T q) −1 Φq− I)M −1 Q e (3.42) where ΦTq(ΦqM −1ΦTq) −1 is the weighted pseudoinverse of the Φq matrix (Wang and Huston, 1989). For the redundant systems, the pseudoinversematrix (ΦqM −1ΦTq) + should be taken instead of the inverse matrix (ΦqM −1ΦTq) −1. It is important to note that the Wang-Huston formulation is not suitable for multibody systemswith singularmassmatrices. Therefore, in the numerical examples presented below, the Euler angles are used in order to avoid this issue. 4. Numerical examples 4.1. McPherson suspensions with and without redundant constraints The presented methods have been implemented and tested on two elementary examples of spatial mechanisms. The simplified McPherson strut without redundant constraints shown in Fig. 1a is the first example and the second one is the overconstrained McPherson strut presented in Fig. 1b. Elimination of the redundant constraints from the second example results in the simplifiedMcPherson strut. The idea for consideration of thesemechanisms is taken from (Haug, 1989). Both mechanisms consist of 5 rigid bodies (denoted below as i = 1,2,3,4,5) and have 4 degrees of freedom. Two of the DOFs are local mobilities (for bodies 3 and 4). Dimensions of the systems are presented in Table 1 and employ symbols presented in Figs. 1a and 1b.Moreover, origins of the local coordinate frames are located in the centers ofmass of the bodies, so that certain expressions get simplified, e.g. the gravitational torques reduce to zero. It is assumed that the centers of mass are placed in the middle of the body i lengths: |AB|, |CD|, |JK|, |IL| and |LN|, respectively. The kinematics of the mechanisms is described using the absolute coordinates. The Euler angles are used for description of the orientation in order to obtain nonsingular mass matrices. Masses of all bodies are:mi =1kg andmoments of inertia are: Jxi = Jyi = Jzi =0.1kgm 2. Table 1.Dimensions of theMcPherson struts Body 1 2 3 4 5 Symbol |AB| |AAa| |AAb| |CD| |CE| |CF| |BC| |FG| |EH| |GI| |JK| |IL| |LN| Dim. [m] 0.3 0.025 0.025 0.45 0.1 0.25 0.05 0.1 0.15 0.05 0.3 0.5 0.2 Both considered models are loaded in the same manner. The force of the constant value Fz = 40N is applied to the center of mass of the first body in the z direction and the time 1400 M. Pękal, J. Frączek Fig. 1. (a) Simplified McPherson strut, (b)McPherson strut with redundant constraints varying force acting in the y direction Fy = 0.1sin(t)N is applied to body 5. Moreover, it is assumed that the gravity is acting in the negative z direction with gravity acceleration equal to g=−9.80665m/s2. The analyzed methods have been implemented using MATLABr R2012b. The obtained results (positions, velocities and accelerations) were compared with the outcomes of the simu- lations performed in AdamsTM 2013 in order to verify their correctness. To integrate the equ- ations of motion, the ode45 method based on the Runge-Kutta scheme was used in MATLAB. Simulations were performed on the personal computer equipped with Intelr CORETM i5 CPU M520@ 2.40GHz 2.40GHz processor, 4GB of RAMand 64-bit Microsoftr Windowsr 7 Home Premium operating system. Computational times were measured using tic and toc functions. Moreover, MATLAB program was run in the single thread mode in order to avoid problems with multi thread timemeasurements. Each simulation was carried out for 10s motion. Three stabilization cases were performed for each method: without constraint stabilization, with the Baumgarte stabilization (where α̂ = β̂ = 10) and with the stabilization of the position constraints using the Newton-Raphson method. Moreover, four ode45 error tolerances were considered: AbsTol = RelTol ∈ {1e-10; 1e-8;1e-6;1e-3}, whereAbsTol andRelTol are absolute and relative error tolerances respective- ly. 4.2. Numerical results – comparison Results of the computational time versus the error tolerances are presented in Figs. 2-5. Figures 2a and 3a contain the computation time for simulations without constraint stabiliza- tion, Figs. 2b and 3b present the results for computations with the Baumgarte stabilization, and Figs. 4 and 5 depict outcomes from simulations with stabilization of the position constra- ints. Note that Fig. 5b shows the results presented in Fig. 5a but for clarity the Wang-Huston formulation results are excluded. In the most examples, PUTD methods turn out to be the slowest. Note that the PUTD-H andPUTD-G give almost the same computational time in the case of the simplifiedMcPherson strut, while for the overconstrained McPherson strut, the PUTD-G seems to be faster than the PUTD-H. It may be due to inefficient implementation of the transformation of the H matrix (Eq. (3.20)). For the remaining methods, the results of most simulation cases are very close to each other. However, in the case of the redundant McPherson strut with stabilization of the position constraints, the Wang-Huston formulation is many times slower than other methods, because of the loss of accuracy of the velocity constraints. Hence, theQR decomposition can be recognized as themost reliable and efficientmethod though, its results are a bit slower than the outcomes from theWang-Huston formulation in some cases. Comparison of natural complement formulations for multibody dynamics 1401 Fig. 2. SimplifiedMcPherson strut: (a) without constraint stabilization, (b) with the Baumgarte stabilization Fig. 3. OverconstrainedMcPherson strut: (a) without constraint stabilization, (b) with the Baumgarte stabilization Fig. 4. SimplifiedMcPherson strut with stabilization of the position constraints 1402 M. Pękal, J. Frączek Fig. 5. OverconstrainedMcPherson strut with stabilization of the position constraints: (a) with and (b) without results of theWang-Huston formulation The simulationswith stabilization of the constraints take longer than the simulationswithout stabilization, which should be expected due to greater computational cost. There are other publications that compare the computation time for spatial systems, e.g. (Mariti et al., 2011, 2010), where results for the Wang-Huston formulation are not presented (when the spatial systems are examined). This is because the Euler parameters were used there for description of the orientation, which resulted in the singularmassmatrix and, consequently, theWang-Huston formulation cannot be employed. Moreover, the zero eigenvaluemethod and the Schur decompositionmethodwere considered separately by Mariti et al. (2011, 2010), Pennestr̀ı and Valentini (2007) but the matrix A (Eq. (3.10)) is symmetric, so the Schur decomposition is equivalent to the eigenvalue problem (Golub andLoan, 1996). Therefore, thesemethods can be considered together as it is done in the current paper. Furthermore, only one type of the PUTD method was examined by Mariti et al. (2011, 2010), and in our paper 2 variants of that method are analyzed. Comparing to Mariti et al. (2010), the shorter qualitative results are obtained for the co- ordinate partitioning method. This is due to the fact that in Mariti et al. (2010) coordinate partitioning was done at each time step. In our paper, only one coordinate partitioning is suffi- cient. In the general case, it should bemonitored whether the actual partition of coordinates is still valid. The computation time growswith the increasing accuracy of computationswhat is in general consistentwith the intuition. This is not the case for all results presented byMariti et al. (2011). As mentioned earlier, the results of the zero eigenvalue, Schur decomposition, SVD andQR decompositionmethods are close to each other. The similar conclusionwas obtained byMariti et al. (2010). Note also that inMariti et al. (2010), thePUTDmethod gaves results thatwere close to the mentioned methods. This is not the case in the present work as computations for both typesof thePUTDmethodusually take longer than for thezeroeigenvalue, Schurdecomposition, SVD or QR decomposition methods. This might be caused by inefficient implementation of the Hmatrix transformation (Eq. (3.20)). 5. Conclusions Computational effectiveness of the simulation algorithms strongly depends on the selected me- thod. In the current papermethods that are based on the orthogonal complement are compared: Comparison of natural complement formulations for multibody dynamics 1403 zero eigenvalue formulation, Pseudo Upper Triangular Decomposition (PUTD), Schur decom- position, Singular Value Decomposition (SVD) and QR decomposition. Moreover, two schemes equivalent to the orthogonal complement methods are also considered: coordinate partitioning and Wang-Huston formulation. The effectiveness of these methods for two elementary mecha- nisms –McPherson struts with and without redundant constraints are considered.Moreover, it is worth noting that the obtained results are comparable with the outcomes from the previo- us publications (Mariti et al., 2010; 2011). The comparison with the mentioned publications is described in details in Section 4. The most robust and one of the most efficient formulation is the method based on the QR decomposition, although Wang-Huston formulation is faster in some cases. This is because the Wang-Huston formulation proved to be the slowest method in simulations of the redundant mechanismwith stabilization of the position constraints. It is due to the loss of accuracy of the velocity constraints, which are not stabilized in that case. The other algorithms can be slower because of several reasons, e.g. theymay require the solution of the eigenvalue problem,which is numerically expensive.Moreover, the computation time is shorter for simulations without stabi- lization than for the cases with stabilization. Despite this, there is a decrease in the constraints accuracy in theunstabilized simulations,whichmay causedifficulties in longermultibodymotion analyses. Acknowledgements This work has been supported by the National Science Centre of Poland under grant No. DEC-2012/07/B/ST8/03993. References 1. Amirouche F.M.L., 2006, Fundamentals of Multibody Dynamics, Theory and Applications, Birkhäuser Boston 2. Awrejcewicz J., 2014,Ordinary Differential Equations and Mechanical Systems, Springer Inter- national Publishing 3. Awrejcewicz J.,KudraG., 2005,Modeling,numerical analysis andapplicationof triplephysical pendulum with rigid limiters of motion,Archive of Applied Mechanics, 74, 11-12, 746-753 4. Awrejcewicz J., Kudra G., 2014,Mathematical modelling and simulation of the bifurcational wobblestone dynamics,Discontinuity, Nonlinearity, and Complexity, 3, 2, 123-132 5. Awrejcewicz J., Kudra G., Lamarque C.H., 2003, Dynamics investigation of three coupled rods with a horizontal barrier,Meccanica, 38, 6, 687-698 6. Awrejcewicz J., Kudra G., Lamarque C.H., 2004, Investigation of triple pendulum with impacts using fundamental solutionmatrices, International Journal of Bifurcation and Chaos, 14, 12, 4191-4213 7. Baumgarte J.W., 1972,Stabilization of constraints and integrals ofmotion indynamical systems, Computer Methods in Applied Mechanics and Engineering, 1, 1, 1-16 8. de Jalón J.G., Bayo E., 1994,Kinematic and Dynamic Simulation of Multibody Systems, The Real-Time Challenge, Springer-Verlag, New-York 9. de Jalón J.G.,Gutiérrez-LópezM.D., 2013,Multibody dynamicswith redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraint forces,Multibody System Dynamics, 30, 3, 311-341 10. Eich-Soellner E., Führer C., 1998,Numerical Methods inMultibody Dynamics, B.G. Teubner 11. Frączek J., Wojtyra M., 2008, Kinematyka układów wieloczłonowych. Metody obliczeniowe, Wydawnictwa Naukowo-Techniczne,Warszawa 1404 M. Pękal, J. Frączek 12. FreeMat v4.1, help 13. Golub G.H., Loan C.F.V., 1996, Matrix Computations, The Johns Hopkins University Press, 3rd edition 14. Hartfiel D.J., 2001,Matrix Theory and Applications with MATLAB, CRCPress 15. HaugE.J., 1989,Computer AidedKinematics andDynamics ofMechanical Systems,Vol. I: Basic Methods, Allyn and Bacon 16. Ider S.K., Amirouche F.M.L., 1988, Coordinate reduction in the dynamics of constrainedmul- tibody systems – a new approach,ASME Journal of Applied Mechanics, 55, 4, 899-904 17. Kim S.S.,VanderploegM.J., 1986,QRdecomposition for state space representation of constra- inedmechanical dynamic systems,ASMEJournal ofMechanisms, Transmissions, andAutomation in Design, 108, 2, 183-188 18. Kincaid D., Cheney W., 2002,Numerical Analysis: Mathematics of Scientific Computing, Bro- oks/Cole, 3-rd edition 19. Kunkel P., Mehrman V., 2006, Differential-Algebraic Equations, European Mathematical So- ciety 20. Malczyk P., Frączek J., 2012, A divide and conquer algorithm for constrained multibody system dynamics based on augmented Lagrangianmethodwith projections-based error correction, Nonlinear Dynamics, 70, 1, 871-889 21. Mani N.K., Haug E.J., Atkinson K.E., 1985, Application of singular value decomposition for analysis of mechanical system dynamics, ASME Journal of Mechanisms, Transmissions, and Automation in Design, 107, 1, 82-87 22. Mariti L., Belfiore N.P., Pennestr̀ı E., Valentini P.P., 2011, Comparison of solution strategies for multibody dynamics equations, International Journal for Numerical Methods in En- gineering, 88, 7, 637-656 23. Mariti L., Pennestr̀ı E., Valentini P.P., Belfiore N.P., 2010, Review and comparison of solution strategies for multibody dynamics equations, The 1st Joint International Conference on Multibody System Dynamics, Lappeenranta, Finland 24. MATLABr, help 25. Negrut D., Serban R., Potra F.A., 1997, A topology-based approach to exploiting sparsity in multibody dynamics: joint formulation,Mechanics of Structures and Machines, 25, 2, 221-241 26. Ostalczyk P., 2008,Wybrane zagadnienia rachunku wektorowego i macierzowego dla robotyków, Wydawnictwo Politechniki Łódzkiej, Łódź 27. Pennestr̀ı , Valentini P.P., 2007, Coordinate reduction strategies in multibody dynamics: A review,Proceedings of the Conference on Multibody System Dynamics 28. Pękal M., 2012, Porównanie metod całkowania RRA w analizie dynamiki układów wieloczłono- wych (in Polish), praca dyplomowamagisterska, PolitechnikaWarszawska,Warszawa 29. Walton W.C., Steeves E.C., 1969, A New Matrix Theorem and Its Application for Establi- shing IndependentCoordinates forComplexDynamical SystemswithConstraints,NASATechnical Report: NASATRR-326, NASA, Washington D. C. 30. Wang J.T., Huston R.L., 1989, A comparison of analysis methods of redundant multibody systems,Mechanics Research Communications, 16, 3, 175-182 31. Wehage R.A., Haug E.J., 1982,Generalized coordinate partitioning for dimension reduction in analysis of constrained dynamic systems,ASME Journal of Mechanical Design, 104, 1, 247-255 32. Wojtyra M., Frączek J., 2013, Comparison of selected methods of handling redundant con- straints in multibody systems simulations, Journal of Computational and Nonlinear Dynamics, 8, 2, 021007 (1-9) Manuscript received March 5, 2015; accepted for print March 30, 2016