Differential–algebraic equations in multibody system

Numerical Algorithms19(1998)183–194183 Differential–algebraic equations in multibody systemmodeling∗D.PogorelovDepartment of Applied Mechanics,Bryansk State Technical University,b.50let Oktyabrya7,241035Bryansk,RussiaE-mail:pogorelov@bitmcnit.bryansk.suNumerical methods for the efficient integration of both stiff and nonstiff equations of motion of multibody systems having the form of differential–algebraic equations(DAE)ofindex3are discussed.Linear multi-step ABM and BDF methods are considered for thenon-iterational integration of nonstiff DAE.The Park method is proposed for integration ofstiff equations.Keywords:DAE,multibody systems,linear multi-step methods1.DAE and multibody systemsThere exist two main reasons why equations for multibody systems(MBS)are DAE.First,if an MBS has kinematic loops,it is,in general,impossible to express analytically all kinematic quantities, e.g.,body positions,as functions of indepen-dent generalized coordinates.In this case,usually redundant coordinates are used and the equations of motion must be complemented by algebraic constraint equations (the descriptor form of equations).Second,large systems consisting of hundreds or even thousands of bodies can be efficiently modeled if they are preliminarily split into several subsystems by cutting some joints.The closure conditions again lead to DAE.Numerical methods which could be used in program packages for the simulation of various MBS should take into account the potentially large size of equations and a long integration period.The main requirements for the numerical methods could be the following:•reliability in direct solving of large DAE;•non-iterational process of calculation in the case of nonstiff DAE;•possibility of solving ODE,DAE and algebraic equations;•variable step size and(eventually)variable order.∗Supported by the RFBR under grant98-01-00782.©J.C.Baltzer AG,Science Publishers184 D.Pogorelov/DAE in multibody system modelingHere,some linear multi-step methods having a predictor–corrector form which satisfy the requirements for nonstiff and stiff DAE are discussed.Euler–Lagrange equations with constraints in index3,2and1formulations are given in section2.The backward differentiation formula(BDF)methods were proposed in[4]for the direct numerical integration of general DAE.Applications of the BDF methods to modeling constrained mechanical systems in the index3formulation are studied in[8]. In[5]the reduction of the system to index2by differentiating constraint equations and introducing additional Lagrange multipliers was suggested,so that the constraints are always satisfied.In particular,a5-step predictor–corrector technique was proposed for solving implicit nonlinear equations with only one“corrector”iteration.A different predictor–corrector technique exploiting a parallel solution of index3and2equations is considered in[10].It allows avoiding one step in the technique studied in[5]and, respectively,to reduce the computational effort.The Adams–Bashforth–Moulton(ABM)method is widely used for the solution of ODE[13].Its computational implementation for constrained mechanical systems based on a generalized coordinate partitioning algorithm was suggested in[14].The stabilization of the PECE formulation of the ABM method for directly solving non-stiff Euler–Lagrange equations with constraints is considered in[10].The main idea of the stabilization technique consists in solving systems of different indices for the “corrector”and“evaluation”steps.In section3some modifications of this method are studied.A computational implementation of the implicit Park method[9]for direct solu-tion of stiff DAE is proposed in section4.An approximation of the Jacobian matrix aiming towards a reduction of the computational effort for its evaluation is studied. Some examples and tests for the discussed methods are given in section5.2.Equations of motionConsider the equations of motion of a MBS with ideal and holonomic constraints. If the system has kinematic loops,redundant coordinates q are used for the description of the system configuration.In this case,the equations have the following form (Euler–Lagrange equations)[2,15]:M(q,t)¨q=b(q,˙q,t)+G(q,t)λ,(2.1)g(q,t)=0,G T=∂g/∂q T,(2.2)where equation(2.2)yields the constraint equation for the redundant coordinates,and λis the vector of Lagrange multipliers.DAE(2.1),(2.2)have index3.Differentiation of equation(2.2)gives equations of index2:M(q,t)˙v=b(q,v,t)+G(q,t)λ,(2.3)G T v+g (q,t)=0,v=˙q,g =∂g/∂t,(2.4)D.Pogorelov/DAE in multibody system modeling185 and equations of index1:M(q,t)a=b(q,v,t)+G(q,t)λ,G T a+g (q,v,t)=0,a=¨q,g =˙G T˙q+(∂g/∂t)˙.All three systems are(in the case of consistent initial conditions)equivalent from the mechanical point of view,but they are not equivalent with respect to stability[8].The third system,with index1,can be solved as an ODE.The others should always be integrated as DAE.3.Numerical methods for nonstiff DAEThe technique described below stabilizes the ABM scheme for the direct solution of DAE in the Euler–Lagrange formulation.Consider a PECE scheme of order k for calculating the unknowns q i,v i,a i,λi, i=k,k+1,...,which are approximations of the variables q(t i),v(t i),a(t i),λ(t i), t i=ih,satisfying the equations of motion.Here,h is the constant step size.The Adams–Bashforth formula is used as a predictor[13],q p i+1=q i+hσ∗(v),v pi+1=v i+hσ∗(a),whereσ∗(ξ):= kj=1βj∗ξi+1−j.The following linear systems must be solved during the correction step:M p a qi+1=b p+G pλqi+1,G p Tδq i+1+g p=0,δq i+1=q i+1−q p i+1,(3.1)M p a v i+1=b p+G pλv i+1,G p T v i+1+g p=0.The superscript p denotes that the values q pi+1,v pi+1are substituted into the correspond-ing matrices.The parameters a qi+1and a v i+1must be excluded from equation(3.1)using theAdams–Moulton relations[13]q i+1=q i+hσ(v)+hβ0v qi+1,v qi+1=v i+hσ(a)+hβ0a qi+1,v i+1=v i+hσ(a)+hβ0a v i+1,σ(ξ):=kj=1βjξi+1−j.After these computations the unknownsδq i+1,v i+1can be obtained from the following equations:M pδq i+1=k q+h2β20G pλqi+1,G p Tδq i+1+g p=0,(3.2)M p v i+1=k v+hβ0G pλv i+1,G p T v i+1+g p=0,186 D.Pogorelov/DAE in multibody system modelingwherek q=hM p c q(v,a)+h2β20b p,k v=M p c v(v,a)+hβ0b p,c q(v,a)=σ(v)−σ∗(v)+β0v i+hβ0σ(a),c v(v,a)=v i+hσ(a).Thefinal evaluation can be done in three different ways.Thefirst is described in[10].According to this method(ABM1)the values a i+1,λi+1are calculated from the equationsM c a i+1=b c+G cλi+1,G cT a i+1+g c=0,(3.3) where the superscript c corresponds to the substitution of the variables q i+1,v i+1 into the matrices.This formulation guarantees the accuracy of order δq i+1 2ofsatisfaction of the original constraint equation g(q)=0,whereas the differentiated constraint equation error is only of order δq i+1 ,although both errors do not increase over a long integration time(see section5).Further,two modifications(ABM2and ABM3)decrease the error in satisfaction of the differentiated constraint equation.For the ABM2method equation(3.3)should be supplemented by the two last equations of(3.2)evaluated for q i+1,v i+1.In other words,in the course of evaluation equations of indices1and2are solved.The ABM3modification in comparison with ABM2 includes additionally thefirst two equations of(3.2),i.e.,equations of indices1–3 should be solved during thefinal evaluation.The convergence of the considered modifications ABM2and ABM3for a constant step size is similar to that given in[10] for the ABM1method.Remark1.Eliminating a qi+1and a v i+1from equation(3.1)and later computing theunknownsδq i+1,v i+1from equation(3.2)require exactly the same number offloating-point operations as eliminatingδq i+1,v i+1and the computation of a qi+1and a v i+1.Thisis valid both for the Range-Space Method and the Null-Space Method when they are used for solving the corresponding linear equations[1].The implicit BDF method requires the calculation of partial derivatives of the mass matrix M(q,t)and the vector of generalized inertial forces with respect to q and˙q.For MBS with many degrees of freedom this leads to considerably increased computational efforts.The PECE scheme seems to be more efficient in the case of nonstiff DAE[12].To improve the accuracy in the calculation of velocities and Lagrange multipliers,a parallel solution of equations(2.1)–(2.4)is proposed in[10]. Here,the main idea of the method is considered briefly.The predictor relations for the BDF method of order k,1 k 5,are given by the interpolation formulasq p i+1=kj=1l j q i+1−j,v pi+1=kj=1l j v i+1−j.D.Pogorelov /DAE in multibody system modeling 187The accelerations a qi +1,a v i +1are excluded from equation (3.1)using the backwarddifferentiation formulasα0q i +1=−ρ∗(q )+hv q i +1,α0v q i +1=−ρ∗(v )+ha q i +1,α0v i +1=−ρ∗(v )+ha v i +1,ρ∗(ξ):=kj =1αj ξi +1−j .This finally leads to equations similar to (3.2).The implementations of the ABM and BDF methods use a variable order and step size according to Krogh’s technique [6,13].Both methods can be used for the non-iterational solution of pure algebraic equations,which is important for inverse kinematic problems.4.Numerical method for stiff DAEThe explicit predictor–corrector solvers discussed in section 3cannot be used for stiff DAE [12].Consider for this purpose the A-stable Park method,which is a linear combination of the second and the third order BDF methods [9].The Park method has a smaller truncation error in comparison with the A-stable BDF method of the second order and is,in addition,more accurate in the case of oscillatory stiff ODE.The predictor formulasq p i +1=4q i −6q i −1+4q i −2−q i −3for the coordinate vector andv p i +1=16h10q p i +1−15q i +6q i −1−q i −2 for the velocity vector are obtained from the interpolation polynomial and the Park formula,respectively.Substitution of the Park formulas v i +1=v p i +1+δq i +1 α,α=0.6h ,a i +1=16h(10v i +1−15v i +6v i −1−v i −2)=a pi +1+δq i +1/α2into equations (2.1),(2.2)gives the following nonlinear algebraic equations to be solved for δq i +1,λi +1:M q p +δq α2a p +δq =α2 b q p +δq ,v p +δq/α +G q p +δq λ ,(4.1)g q p +δq =0.Here,the lower indices as well as the dependence of matrices on t are omitted.Equations (4.1)are usually solved using Newton–Raphson iterations M k +J k ∆q k +1=−M k α2a p +δq k +α2 −k k +Q k +G k λk ,G k T ∆q k +1=−g k ,δq k +1=δq k +∆q k +1,k =0,1,...,δq 0=0,188 D.Pogorelov /DAE in multibody system modelingwhere the Jacobian matrix J k was introduced;b =−k +Q with the vectors k ,Q being the generalized inertia and applied forces.The superscript k denotes the substitution of δq k into the corresponding matrices.The convergence proof of the Park method is similar to that given in [8]for the BDF method.The implementation of the Park method includes a variable step size ensuring the given absolute error tolerance (local).The solver starts with the implicit method of the first order.Evaluation of the Jacobian matrix J k leads to a considerably increasing computa-tional effort.It is worth replacing this matrix with an approximated one [7].Stiffness in MBS simulations is caused in many cases by applied forces (elasticity,damping,contact,etc.).In these cases the motion of the system is “slow”,i.e.,the mass ma-trix M is nearly constant for several integration steps and the inertia forces k are small.For such systems it seems to be much more efficient to calculate an approximated Ja-cobian matrix taking into account only the leading terms of the vector of generalized forces Q .Other matrices and vectors in the differential part of the equations are not differentiated.Consider a bipolar force element as an example.The force acts along the line passing through two points A and B (attachment points).The points are given in body-fixed coordinate systems by the constant vectors ρi ,ρj ,where i and j are the interacting bodies.The force model used isF j =r ij rf (r ,v ,t )=e ij f ,where r is the distance between the attachment points,v =˙r ,the vector r ij connects the attachment points.The corresponding generalized force can be found as Q j = D T i+B T i ρi −D T j −B T j ρj e ij f =C T (q )e ij (q )f r (q ),v (q ,˙q ,t ),t .The matrices D i ,B i ,D j ,B j define velocities of the centers of mass and angular velocities of bodies asv i =D i ˙q ,ωi =B i ˙q ,v j =D j ˙q ,ωj =B j ˙q.The approximated Jacobian matrix corresponding to the force is defined by taking into account only the derivatives of the force value f (r ,v ,t )with respect to δq .The derivatives of the matrix C and the vector e ij are neglected.Thus,the Jacobian matrix can be computed by the expression J j =− f r α2+f v α C T e ij e T ij C ,which can easily be obtained both numerically and in a fully symbolic form.Consider a very simple example,to estimate the influence of the Jacobian matrix approximation.A rigid rod AB of length l and mass M rotates around its end A in the plane Oxy ,φis the angle of rotation.The point A coincides with the origin of the coordinate system.A linear dissipative force element connects the other end B of the rod with a fixed point (0,a )and produces the force f =−νv .The driving torqueD.Pogorelov/DAE in multibody system modeling189Table1Absolute error in position and computational effort for the rod after one revolution.lgε−4−5−6−7−8eφ3.3×10−21.2×10−23.1×10−31.1×10−31.1×10−4NFE410726140827915691is L.Simulation results for the global error in the rotation angle eφand the number of function evaluations(NEF)versus the requested error toleranceε(local)are given in table1for the following parameter values and initial conditions:l=0.4m,a=0.5m,M=1kg,L=1Nm,ν=100Ns m−1,φ(0)=0,˙φ(0)=0.The simulation was performed for one full revolution of the rod.The approximated value of the Jacobian matrix only differs considerably from the exact one if the force element is nearly parallel to the rod,i.e.,the angleφlies in the neighborhood of0orπ.But in this case,the Jacobian is small compared with the moment of inertia of the rod and the error in its value has no influence on the convergence of the iterations.The simulation results for the approximated and the exact Jacobian matrices coincide for the given parameter values.5.Numerical examplesConsider a mechanism built by a platform and6actuators(figure1).The mecha-nism consists of25bodies connected by36joints.The position of the platform relative to the inertial reference frame Oxyz is defined by3translational and3rotational coor-dinates.The actuators have identical kinematic schemes:spherical joints C,rotational joints A,B,E,cardan joints F.One translational d.o.f.corresponds to every actua-tor.The axes of the rotational joints are perpendicular to the vertical plane OAB.The joints B lie in the plane Oxy.Their coordinates are(0,1.5),(−1.2990381,0.75),(−1.2990381,−0.75),(0,−1.5),(1.2990381,−0.75), (1.2990381,0.75).The joints A lie exactly above the joints B,AB=1m.The coordinates of the joints A,E,F in the coordinate systemfixed to body2are(0,0.231),(−0.115,−0.2),(−0.115,0.2).The axes of the coordinate systemfixed to body1are parallel to those of the reference frame in the position of the mechanism shown infigure1.The coordinates of the spherical joints C in this system are(0.3000623,0.98146,0),(−1,0.23086882,0),(−1,−0.23086882,0),(0.3000623,−0.98146,0),(0.6999377,−0.7505913,0),(0.6999377,0.7505913,0).190 D.Pogorelov/DAE in multibody system modelingFigure1.Mechanism with6degrees of freedom.Table2Inertia parameters of the mechanism.Body Mass(kg)Moments of inertia(kg m2)CommentsI x I y I z110003354216752250.2335710.687510.6875Homogeneous rod,length1.5m421 6.5Center of mass lies on BE,0.8m from point B5388.5Center of mass lies on BE,0.5m from point EInertia and some geometrical parameters are presented in table2.Moments of inertia for bodies3–5are given relative to the corresponding principal axes perpendicular to the plane OAB.The model has12closed loops and6degrees of freedom.The configuration is determined by42coordinates.The equations of motion include42second-order differential equations and36algebraic constraints.To test the ABM and BDF methods of section3,consider the simulation of the mechanism dynamics taking into account only gravity.The mechanism in this case is a conservative mechanical system.The equations of motion are nonstiff but the motion of the mechanism is very complex.Let the initial coordinates and velocity of the center of mass of body1bex10=y10=0,z10=3.4,v x10=1,v y10=v z10=0,D.Pogorelov /DAE in multibody system modeling 191Figure 2.Animation results for the mechanism.the angles of initial orientation of the body be zero (the axes of the body-fixed frame are parallel to the axes of the reference frame)and the initial angular velocity beωx 0=2,ωy 0=ωz 0=0.The simulation time is 1.5s.Some animation results are shown in figure 2.Figure 3(a)shows the maximal errors in satisfaction of original and differentiated constraints e a =max g (q ) 2 ,e d =max G T ˙q 2 versus the absolute error tolerance ε(local)for the BDF and various modifications of the ABM methods.The maximal absolute errors in position and velocity variables e q ,e v at t =1.5s are given in figure 3(b).Figure 3(c)shows the maximal relative deviation of the mechanical energy dE from its initial value.Finally,the numbers of function evaluations are compared in figure 3(d).The ABM2and ABM3methods have coinciding results for the data in figures 3(b)–(d).The errors both in original and differentiated constraint equations do not increase over the integration time.As was expected,the original constraints are “exactly”satisfied for the ABM3method and the differentiated constraints –for the ABM2method.The ABM1method shows the worst accuracy in satisfaction of constraint equations and requires more function evaluations in comparison with other modifications of the ABM method.The BDF method requires fewer function evaluations for integration with a lower accuracy.The next simulations of the mechanism dynamics were performed taking into account a controlled motion of the platform (body 1).The desirable motion of the platform is defined byx 01=0.2sin 3t ,y 01=0.3cos 6t ,z 01=3.2+0.1sin 3t ,α01=β01=γ01=0,where α1,β1,γ1are the orientation angles.The control forces in the actuators are F ai =−C x ai −x 0ai (t ) −D ˙x ai −˙x 0ai(t ) ,i =1,...,6,x ai being the translational coordinates corresponding to the actuators,x 0ai (t )their values for the desirable motion obtained from solving the inverse kinematic problem,C =3.3×107N /m,D =4.7×104Ns /m.The equations of motion are stiff and cannot be solved by explicit methods.The Park method is used for the simulation of192 D.Pogorelov/DAE in multibody system modelingparison of numerical methods ABM1(◦),ABM2( ),AMB3(♦),and BDF( ).two variants of the mechanism.Thefirst variant(V1)is the original one.Bodies4 and5are weightless for the second variant(V2).The actuators for V2can be modeled by bipolar force elements.The Jacobian matrices for the control forces are calculated according to the approach of section4.The initial conditions for body1arex10=0,y10=0.3,z10=3.16,v x10=0.6,v y10=0,v z10=0.3,α10=β10=γ10=0,ωx10=ωy10=ωz10=0.The value z10differs from z01(0)and a transient process is expected.Simulation results for the absolute errors in position and velocity variables e q, e∗v,e v as well as the number of function evaluations N,N∗vs.the absolute error toleranceε(local)are shown infigure4.The values e q,e v,N correspond to the whole interval of simulation1.9s,whereas the calculation of e∗v,N∗is started at0.1s from the beginning of the simulation.The results are quite similar for both variants of the mechanism.The last example is a squeezer mechanism used as a rule as a benchmark for the comparison of different DAE solvers[1,3,11].Table3presents errors and com-putational effort versus the requested absolute error toleranceε(local)for a30msD.Pogorelov/DAE in multibody system modeling193Figure4.Errors and computational effort for V1( )and V2( )of the mechanism simulation.Table3Relative errors and computational effort for the squeezer.lgε−3−4−5−6 BDF e q3.0×10−62.2×10−62.7×10−71.6×10−8e v3.0×10−31.6×10−33.2×10−41.7×10−5NFE55577311631718 ABM e q2.2×10−64.6×10−75.8×10−85.3×10−9e v1.0×10−32.0×10−43.4×10−55.7×10−7NFE6057599211099 simulation.Here,e q,e v are the maximal relative errors in positions and velocities, parisons with the analogous simulation results given in[3]for the ODASSL solver show a decreased number of function evaluations.The simulation results of this section were obtained with the help of the program package Universal Mechanism.6.SummaryThe linear multi-step ABM and BDF methods in the predictor–corrector form can be successfully used for solution of nonstiff DAE in the Euler–Lagrange formulation. Simultaneous solution of equations of different indices stabilizes the ABM method and decreases the computational efforts for the BDF method.The Park method is proposed for simulation of stiff constrained mechanical systems.References[1]T.Andrzejewski,H.G.Bock,E.Eich and R.von Schwerin,Recent advances in the numericalintegration of multibody systems,in:Advanced Multibody System Dynamics–Simulation and Software Tools,ed.W.Schiehlen(Kluwer Academic,Dordrecht,1993)pp.127–151.194 D.Pogorelov/DAE in multibody system modeling[2]G.B.Efimov and D.Pogorelov,Some algorithms for computer-aided generation of multibody systemequations,Preprint,Inst.Appl.Math.RAS84,Moscow(1993)(in Russian).[3] C.Fuhrer and B.Leimkuhler,Numerical solution of differential–algebraic equations for constrainedmechanical motion,Numer.Math.59(1991)55–69.[4] C.W.Gear,The simultaneous numerical solution of differential/algebraic equations,IEEE Trans.Circuit Theory18(1971)89–95.[5] C.W.Gear,G.K.Gupta and B.Leimkuhler,Automatic integration of Euler–Lagrange equations withconstraints,put.Appl.Math.12,13(1985)77–90.[6]R.D.Grigorieff,Numerik Gew¨o hnlicher Differentialgleichungen,Band2,Mehrschrittverfahren(Teubner,Stuttgart,1977).[7] E.Hairer and G.Wanner,Solving Ordinary Differential Equations II,Stiff and Differential–AlgebraicProblems(Springer,Berlin/Heidelberg/New York,1996).[8]P.L¨o tstedt and L.R.Petzold,Numerical solution of nonlinear differential equations with algebraicconstraints I:Convergence results for backward differentiation formulas,put.46(1986) 491–516.[9]K.C.Park,An improved stiffly stable method for direct integration of nonlinear structural dynamicequations,J.Appl.Mech.(June1975)464–470.[10] D.Pogorelov,Numerical modelling of the motion of systems of solids,Comput.Math.Math.Phys.35(4)(1995)501–506.[11]W.Schiehlen,ed.,Multibody Systems Handbook(Springer,Berlin,1990).[12]L.F.Shampine,Numerical Solution of Ordinary Differential Equations(Chapman&Hall,NewYork,1994).[13]L.F.Shampine and M.K.Gordon,Computer Solution of Ordinary Differential Equations,The InitialValue Problem(Freeman,San Francisco,CA,1975).[14]R.A.Wehage and E.J.Haug,Generalized coordinate partitioning for dimension reduction in analysisof constrained dynamic systems,J.Mech.Design104(1982)247–255.[15]J.Wittenburg,Dynamics of Systems of Rigid Bodies(Teubner,Stuttgart,1977).。