GENG Guo-zhi LIU Jian-wen DING Jie-yu

(College of Information Engineering,Qingdao University,Qingdao Shandong 266071,China)

0 Introduction

Multibody system[1-2]is connected multiple objects (rigid body and elastomer/soft body,particle,etc.)of system in a certain way.In weapons,robots,aviation,machinery and so on national defense and national economic construction,such as,aircraft launch system,robot,vehicle,such as large civil machinery mechanical system can be attributed to body systems.As the national economy and national defense construction to improve the mechanical system dynamic performance requirements,need for large complex mechanical system dynamics analysis and forecast accurately and quickly.

Large complicated mechanical systems are often closed loop system,using the theory of multibody system dynamics to establish the dynamics equation is generally with differential algebraic equation[3]of Lagrange multiplier,can be expressed as follows：

There are two ways to solve the differential algebraic equations.One is the direct integral method,combines acceleration constraint equation and dynamic equation of integral.Another way is condensed and method,by using the equation of the matrix decomposition independent coordinate system consists of a set of appropriate generalized coordinates,to pure differential equations,differential algebraic equation can be converted to then integral.Commonly used direct integral method has the Euler method and Runge-Kutta method,there are some higher order numerical integral method,such as Newton-Romberg integral and Gauss integral[4],etc.

Higher order numerical integral method is within each time step to integral of integrand,integrand is obtained by interpolating fitting,we usually use interpolation method with Newton interpolation,Lagrange Interpolation, successive linear interpolation, equidistant node interpolation,piecewise interpolation and spline interpolation.Based on the Lagrange Interpolation of each step in the long function interpolation,respectively using Romberg and Gauss integral dynamics simulation.

1 Discrete Euler-Lagrangian equation and numerica integral method

Hamilton’s principle can be expressed as the following：

Where q,q·is generalized coordinates and generalized velocity,S is Hamilton integral,L is the Lagrange function,λ is the Lagrange multiplier,Φ is a vector corresponding constraints

Φ(q,t)=0（2）

Using the variation method,by（1） the Euler-Lagrange equation can be obtained as follows：

Equation （3） is a differential algebraic equation,we usually can use the Euler method and Runge-Kutta method to solve.

By using the discrete variation principle,the integration time divided into N time interval h=tf/N,using the variation method,discrete Euler-Lagrange equation is obtained

Where DjLd(j=1,2)is j th partial derivative of Ld.

In the time interval [ti,ti+1]using Lagrangian difference,first introduced a parameter τ∈(0,1),where

State variables q(t)and its first derivative q·(t)in this time interval can be difference

So the Hamilton function integrals on the interval [ti,ti+1]can be approximate as follows

Using the discrete Hamilton principle,it can be discrete Euler-Lagrangian equation[5]is obtained（10）

Using Gauss quadrature formula,which has maximal degree of accuracy 2n-1 for a fixed number n of quadrature points,we can get

Where Aris the weight and τ~ris the quadrature point.When n=2,l

Using Romberg quadrature formula,which has a degree of accuracy 2log2(n-1)for a fixed number n of quadrature points,we also can get（11）,here we give n=8 quadrature points and weights respect to the interval[-1,1]below in Table 1

Tab.1 The 9 quadrature points and weights of Rombergquadrature formula with respect to the interval[-1,1]

Using the Romberg,Gauss integral can be the solution of equation.

2 Numerical example

Make the horizontal axis as the X-axis，longitudinal axis as Y-axis to establish a coordinate system，to put two balls in the positive X axis horizontal position of the starting point for the double pendulum,the following figure shows the double pendulum state at a time t.The quality of pendulum ball is m1=2,m2=1,length of the rod is l1=1,l2=1,Let pendulum ball 1 location coordinates as(x1,y1),pendulum ball 2 location coordinates as(x2,y2),The state variables of double pendulum is q=[x1,y1,x2,y2].

The kinetic energy of the system：

Based on the above parameters,For equation （3） simulation using the Euler method and Runge-Kutta method,for equation （10） is simulated using Romberg integration and Gauss integration points.Where h=0.01,Euler simulation results diverge,Figure 2 shows the Euler simulation energy curve and the position of the curve at h=0.001,using Runge-Kutta method,Romberg integration,Gauss integration simulation energy curve and position curve at h=0.01 shown in Figure 3,Figure4,Figure5.

Euler integration stability is poor,so basically do not use it；there are divergent trends for Runge-Kutta method,available in small steps,radiation in big steps；Romberg integration and Gauss quadrature are used in big steps,energy error floating in a small area,Gauss integrals is better than Romberg.

Euler method,Runge-Kutta method,Romberg integration,Gauss integration at h=0.001 and h=0.01 simulation running time,the maximum energy error,the biggest constraint errors are shown in Table 2,Table 3

Tab.2 Comparison of the methods mentioned with time step h=0.001

Tab.3 Comparison of the methods mentioned with time step h=0.01

3 Conclusion

By comparing four different integration methods in h=0.01 and h=0.001 run time,the energy constraint error and error,we get the following conclusions：Euler integration stability is poor here are divergent trends for Runge-Kutta method,available in small steps,radiation in big steps；Romberg integration and Gauss quadrature are used in big steps,energy error floating in a small area,Gauss integrals is better than Romberg.

［1］Liu Yanzhu,Hong Jiazhen,Yang Haixing.Multibody system dynamics[M].Beijing：Higher Education Press,1989.

［2］Wittenburg J.Dynamics of Systems of Rigid Bodies[M].Teubner,Stuttgrt,1977.

［3］Hong Jiazhen.Computing multi-body system dynamics[M].Beijing：Higher Education Press,1999.

［4］Li Qingyang,Wang Nengchao,Yi Dayi.Numerical Analysis[M].Huazhong University of Science and Technology Press，1986.

［5］Ding Jieyu,Pan Zhenkuan.Higher Order Variational Integrators of Multibody System Dynamics with Constraints[J].Advances in Mechanical Engineering,2014.