多体系统指标2运动方程HHT方法违约校正1)

马秀腾翟彦博谢守勇

∗(西南大学工程技术学院，重庆400715)

†(现代汽车零部件技术湖北省重点实验室，武汉430070)

多体系统指标2运动方程HHT方法违约校正1)

马秀腾∗,†,2)翟彦博∗谢守勇∗

∗(西南大学工程技术学院，重庆400715)

†(现代汽车零部件技术湖北省重点实验室，武汉430070)

采用Cartesian绝对坐标建模方法，完整约束多体系统运动方程是指标3的微分--代数方程(dif f erentialalgebraic equations,DAEs)，数值求解指标3的DAEs属于高指标问题，通过对位置约束方程求导，可使运动方程的指标降为2.位置约束方程求导得到的是速度约束方程.直接求解指标3的运动方程，速度约束方程得不到满足，而且高指标DAEs的数值求解存在一些问题.论文首先采用HHT(Hilber--Hughes--Taylor)直接积分方法求解降指标得到的指标2运动方程，此时速度约束方程参与离散计算，从机器精度上讲速度约束自然得到满足，而位置约束方程没有参与计算，存在“违约”.针对违约问题，采用基于Moore--Penrose广义逆理论的违约校正方法，消除位置约束方程的违约.指标2运动方程HHT方法违约校正，将HHT方法和违约校正方法很好地结合，在数值求解指标2运动方程的过程中，位置约束方程和速度约束方程都不存在违约问题，而且新方法没有引入新的未知数向量，离散得到的非线性方程组的方程数量与原指标2运动方程的方程数量相同，求解规模没有扩大.新方法的实用和有效性通过算例的数值实验得到验证，数值实验也说明新方法保持了HHT方法本身具有的数值阻尼可以控制和二阶精度的特性.最后从非线性方程组的求解规模和计算速度上与其他方法进行了比较分析，说明新方法的优势所在.

运动方程，指标，HHT方法，Moore--Penrose广义逆，违约校正

引言

为了实现多体系统程式化建模与仿真求解，采用Cartesian绝对坐标建模方法，对于完整约束多体系统，此时得到的运动方程是指标3的微分--代数方程 (dif f erential-algebraic equations,DAEs)[1].由于DAEs和常微分方程(ODEs)有本质区别，且一般没有解析解，因此要研究专门针对DAEs的数值求解算法[2-7].

已有的运动方程求解方法可分为增广法和缩并法[2]，也可分为直接法和状态空间法[8].在直接法中，结构动力学中的直接积分方法，如Newmark方法、HHT(Hilber-Hughes-Taylor)方法、广义α方法、θ1方法等，有广泛的应用[8-17].也有状态空间法和广义α方法相结合的数值求解算法[11].在直接法中，HHT方法和Newmark方法已集成到MSC.ADAMS软件的求解器上[9,18].

对于指标 3运动方程，对其位置约束方程求导，整体将得到指标2 DAEs形式的运动方程.位置约束方程求导得到的方程称为速度约束方程.直接数值求解指标3运动方程，属于高指标问题[1]，求解存在一些问题，且速度约束方程存在违约问题.直接求解指标2运动方程，位置约束方程违约存在问题.如何同时满足位置约束和速度约束方程，一种求解思路是基于GGL约束稳定[19]的方法，通过引入新的“未知数向量”，将速度约束方程并到指标3运动方程中一起，求解约束稳定指标2运动方程.从机器精度上讲，位置约束和速度约束方程同时得到满足，已有如HHT-SI2方法[8]、HHT-ADD方法[8]、投影方法[10,12]、广义α-SOI2方法[15]、θ1方法[13]和Lie群积分[20-21]等.GGL约束稳定指标2运动方程的求解形式已经集成到MSC.ADAMS求解器中，可以选用Gsti ff,Wstiff,Constant_BDF,Newmark和HHT等方法求解[22].MSC.ADAMS求解器中通过引入两个“导数向量”[1,22]来代替指标3运动方程中的Lagrange乘子和GGL约束稳定引入的“未知数向量”，得到约束稳定指标1运动方程.不论是约束稳定指标2还是约束稳定指标1运动方程，由于引入了新的未知数，离散后得到的非线性方程组的求解规模将变大.

对指标2运动方程中速度约束方程求导，整体将得到指标1 DAEs形式的运动方程，MATLAB软件的ode15s和ode23t函数可以求解DAEs和刚性ODEs，但仅能求解指标1的DAEs[23].求解多体系统指标1的运动方程会存在位置约束和速度约束方程存在违约问题.Yoon等[24]、于清等[25-26]从改进Baumgarte约束违约稳定法角度，Nikravesh[27]从初始条件校正角度，分别基于Moore-Penrose广义逆理论，提出新的违约校正方法.最近Marques等[28-29]在此方法上又有深入的研究.

Gear[30]证明对于力学系统，数值求解指标2运动方程是较好的选择.本文将直接积分中的HHT方法和基于Moore-Penrose广义逆的违约校正方法相结合.用HHT方法求解指标2 DAEs形式的运动方程，位置约束方程的违约通过校正方法消除.新方法仅在求解的过程中增加一步迭代，没有引入新的变量.较之已有的方法，求解速度将得到提高.

1 多体系统运动方程

完整约束多体系统运动方程形式如下

将位置约束方程Φ(q)=0求导，方程(1)的指标由3降为2，方程形式如下

2 指标2 运动方程的HHT方法

对于方程(2)，采用HHT方法离散，并对速度约束方程进行缩放[31]，得到的非线性方程组如下

非线性方程组(3)中方程的数量为n+m，采用Newton-Raphson方法求解此方程组，得

此时数值解使速度约束方程得到满足，而位置约束方程存在违约问题.也就是说从机器精度上讲

3 位置约束方程的违约校正

设由HHT方法求解方程(2)，数值积分到第n+1步后，得到的广义坐标为qn+1，此时位置约束方程存在违约问题.可以在广义坐标qn+1上加入校正项得

位置约束方程Φ(q)=0的变分可表示为

4 算例

如图1所示的曲柄滑块机构，采用Cartesian绝对坐标建模方法，设曲柄AB长为l1，质心坐标为(x1,y1)，连杆BC长为l2，质心坐标为(x2,y2)，滑块质心坐标为(x3,0).

图1 曲柄滑块机构简图Fig.1 Sketch of slider crank mechanism

广义坐标为

位置约束方程为

取曲柄质量m1=3kg，长度l1=0.6m，连杆质量m2=0.9kg，长度l2=1.2m，滑块质量m3=1.8kg.重力加速度g=9.81m/s2.弹簧刚度k=100N/m.

取阻尼系数c=1N·s/m，设求解时间20s，步长h=1ms.

图2 机构θ角时间历程Fig.2 Time evolution of θ for mechanism

图2给出本文的HHT方法和违约校正相结合(HHT-I2-Correction)方法与HHT-SI2方法[8]求解结果的比较，完全吻合，说明方法是有效的.

图3给出HHT方法直接求解运动方程(2)时速度约束违约量的曲线.

图4给出HHT方法直接求解运动方程(2)时位置约束违约量的曲线.

图3 速度约束方程违约曲线Fig.3 Evolution of velocity constraint equations violation

图5 速度约束方程违约曲线Fig.5 Evolution of velocity constraint equations violation

图6 位置约束方程违约曲线Fig.6 Evolution of position constraint equations violation

由图5和图6可知，速度约束和位置约束方程的违约数量级都在10-16～10-14之间，从机器精度上讲，不存在违约问题.

当曲柄滑块机构中阻尼系数c=0时，设求解时间20s，步长h=1ms，图7给出当α=-1/3,-1/6, -1/48时，数值求解过程中系统能量的变化曲线.

图7 求解过程中系统能量的变化曲线Fig.7 Energy of the system during integration

由数值实验可知，本文方法和原HHT方法一样，α在[-1/3,0]区间内，可以通过改变α控制方法数值阻尼的大小，随着α的增大，数值阻尼是减少的.

当α=0，由HHT方法得γ=0.5，β=0.25，这就是梯形法则公式[32].如果α=0，γ=0.5，令β=0，这就是Strmer算法[32]，这两种方法都是没有数值阻尼的.

图8给出通过数值实验得到的本文方法精度的阶，和原HHT方法一样都具有二阶精度.

图8 精度Fig.8 Order of accuracy

由文献[14]的研究可知，为了使求解过程中位置约束和速度约束方程同时得到满足，在求解指标2超定微分--代数方程时，步长取为h=1ms，广义α-SOI2方法、HHT-SI2方法、θ1方法和文献[14]中新广义α方法的求解速度是相当的，θ1方法稍快.前3种方法非线性方程组的求解规模为2n+2m，新广义α方法的求解规模为n+2m.

对于图1所示的模型，阻尼系数c=1N·s/m，设求解时间为20s，步长为h=1ms.相同的编程语言编写程序，在相同配置的电脑上运行，θ1方法的CPU时间为20.46s，本文HHT直接积分校正方法仅为17.38s.可见方法具有求解速度的优势，而且求解的非线性方程组的规模是最小的，仅为n+m.

5 结论

求解多体系统运动方程时，为了得到更精确的解，要求位置约束和速度约束方程都不存在违约问题.

本文用HHT方法直接求解指标2 DAEs形式的运动方程，此时速度约束方程自然满足，为了消除位置约束方程的违约，引入基于Moore-Penrose广义逆理论的违约校正方法，使得位置约束方程得到满足.

数值实验表明本文方法保持了原HHT方法可控数值阻尼和二阶精度的特性，且求解速度比广义α-SOI2方法、HHT-SI2方法、θ1方法和新广义α方法都快.

从形式上看，本文方法非线性方程组求解规模最小.对于复杂多体系统建模与仿真，广义坐标数n和约束方程数m将相应增大，求解规模小的特点将更具优势.

位置校正公式含有与约束Jacobian矩阵相关的m×m维矩阵求逆，对于复杂系统的仿真求解，如果每一步都进行校正，矩阵求逆将会耗费大量计算时间.可以根据计算精度的要求，设定当位置违约量大于某一量值后才进行校正，以提高计算速度.

1 Brenan KE,Campbell SL,Petzold LR.Numerical Solution of Initial-Value Problems in Dif f erential Algebraic Equations.2nd edn.Philadelphia：SIAM,1996

2 潘振宽,赵维加,洪嘉振等.多体系统动力学微分/代数方程组数值方法.力学进展,1996,26(1)：28-40(Pan Zhenkuan,Zhao Weijia,Hong Jiazhen,et al.On numerical algorithms for dif f erential/algebraic equations of motion of multibody systems.Advances in Mechanics,1996,26(1)：28-40(in Chinese))

3 王琪,陆启韶.多体系统Lagrange方程数值算法的研究进展.力学进展,2001,31(1)：9-17(Wang Qi,Lu Qishao.Advances in the numerical methods for Lagrange’s equations of multibody systems.Advances in Mechanics,2001,31(1)：9-17(in Chinese))

4 赵维加,潘振宽.多体系统Euler-Lagrange方程的最小二乘法与违约修正.力学学报,2002,34(4)：594-603(Zhao Weijia,Pan Zhenkuan.Least square algorithms and constraint stabilization for Euler-Lagrange equations of multibody system dynamics.Chinese Journal of Theoretical and Applied Mechanics,2002,34(4)：594-603(in Chinese))

5 Bauchau OA,Laulusa A.Review of contemporary approaches for constraint enforcement in multibody systems.ASME Journal of Computational and Nonlinear Dynamics,2008,3(1)：11005

6 Arnold M.DAE aspects of multibody system dynamics.Report No. 01.Martin-Luther-Universitt Halle-Wittenberg,2016

7 Simeon B.Computational Flexible Multibody Dynamics：A Dif f erential-Algebraic Approach.Berlin Heidelberg：Springer-Verlag,2013

8 Negrut D,Jay LO,Khude N.A discussion of low order numerical integration formulas for rigid and fl xible multibody dynamics.ASME Journal of Computational and Nonlinear Dynamics,2009, 4(2)：21008

9 Negrut D,Rampalli R,Ottarsson G,et al.On an implementation of the Hilber-Hughes-Taylor method in the context of index 3 dif f erential-algebraic equations of multibody dynamics.ASME Journal of Computational and Nonlinear Dynamics,2007,2：73-85

10 丁洁玉,潘振宽.多体系统动力学刚性方程广义α投影法.中国科学：物理学力学天文学.2013,43(4)：572-578(Ding Jieyu,Pan Zhenkuan.Generalized-α projection method for stif fdynamic equations of multibody systems.Scientia Sinica Physica,Mechanica&Astronomica,2013,43(4)：572-578(in Chinese))

11 姚廷强,迟毅林,黄亚宇.柔性多体系统动力学新型广义α数值分析方法.机械工程学报,2009,45(10)：53-60(Yao Tingqiang,Chi Yilin,Huang Yayu.New generalized-α algorithms for multibody dynamics.Journal of Mechanical Engineering,2009,45(10)：53-60 (in Chinese))

12 丁洁玉,潘振宽.多体系统动力学微分--代数方程广义α投影法.工程力学,2013,30(4)：380-384(Ding Jieyu,Pan Zhenkuan. Generalized-α projection method for dif f erential-algebraic equations of multibody dynamics.Engineering Mechanics,2013,30(4)：380-384(in Chinese))

13 马秀腾,翟彦博,罗书强.基于θ1方法的多体动力学数值算法研究.力学学报,2011,43(5)：931-938(Ma Xiuteng,Zhai Yanbo, Luo Shuqiang.Numerical method of multibody dynamics based on θ1method.Chinese Journal of Theoretical and Applied Mechanics, 2011,43(5)：931-938(in Chinese))

14 马秀腾,翟彦博,罗书强.多体动力学超定运动方程广义α求解新算法.力学学报,2012,44(5)：948-952(Ma Xiuteng,Zhai Yanbo, Luo Shuqiang.New generalized-α method for over-determined motion equations in multibody dynamics.Chinese Journal of Theoretical and Applied Mechanics,2012,44(5)：948-952(in Chinese))

15 Jay LO,Negrut D.A second order extension of the generalized-α method for constrained systems in mechanics//Bottasso C L,ed. Multibody Dynamics：Computational Methods and Applications, Springer Science&Business Media B.V.2009.143-158

16 Lunk C,Simeon B.Solving constrained mechanical systems by the family of Newmark and α-methods.ZAMM,2006,86(10)：772-784

17 刘颖,马建敏.多体系统动力学方程的基于离散零空间理论的Newmark积分算法.机械工程学报,2012,48(5)：87-91(Liu Ying, Ma Jianmin.Discrete null space method for the Newmark integration of multibody dynamic systems.Journal of Mechanical Engineering,2012,48(5)：87-91(in Chinese))

18 Orlandea NV.Multibody systems history of ADAMS.ASME Journal of Computational and Nonlinear Dynamics,2016,11(6)：60301

19 Gear CW,Gupta GK,Leimkuhler B.Automatic integration of Euler-Lagrange equations with constraints.Journal of Computational and Applied Mathematics,1985,12&13：77-90

20 Arnold M,Hante S.Implementation details of a generalized-α DAE Lie group method.ASME Journal of Computational and Nonlinear Dynamics,2016,12(2)：021002

21 Arnold M,Cardona A,Brls O.A Lie algebra approach to Lie group time integration of constrained systems//Betsch P,ed.Structurepreserving Integrators in Nonlinear Structural Dynamics and Flexible Multibody Dynamics.Switzerland：Springer International Publishing,2016：91-158

22 陈立平,张云清,任卫群等.机械系统动力学分析及ADAMS应用教程.北京：清华大学出版社,2005(Chen Liping,Zhang Yunqing,Ren Weiqun,et al.Dynamic Analysis of Mechanical Systems and ADAMS Application.Beijing：Tsinghua University Press,2005 (in Chinese))

23 Shampine LF,Reichelt MW,Kierzenka JA.Solving index-1 DAEs in MATLAB and Simulink.SIAM Review,1999,42(3)：538-552

24 Yoon S,Howe RM,Greenwood DT.Geometric elimination of constraint violations in numerical simulation of Lagrangian equations.ASME Journal of Mechanical Design,1994,116(4)：1058-1064

25 Yu Q,Cheng I.A direct violation correction method in numerical simulation of constrained multibody systems.Computational Mechanics,2000,26：52-57

26 于清,洪嘉振.受约束多体系统一种新的违约校正方法.力学学报,1998,30(3)：300-306(Yu Qing,Hong Jiazhen.A new violation correction method for constraint multibody systems.Chinese Journal of Theoretical and Applied Mechanics,1998,30(3)：300-306(inChinese))

27 Nikravesh PE.Initial condition correction in multibody dynamics.Multibody System Dynamics,2007,18(1)：107-115

28 Marques F,Souto AP,Flores P.On the constraints violation in forward dynamics of multibody systems.Multibody System Dynamics, online.DOI：10.1007/s11044-016-9530-y.2016

29 Flores P.A new approach to eliminate the constraints violation at the position and velocity levels in constrained mechanical multibody systems//Flores P,Viadero F,eds.New Trends in Mechanism and Machine Science,5th European Conference on Mechanism Science,Guimaraes,Portugal,September 16-20,2014.Switzerland： Springer International Publishing,2015.385-393

30 Gear CW.Dif f erential-algebraic equation index transformations.SIAM Journal on Scientifi and Statistical Computing,1988,9(1)：39-47

31 Bottasso CL,Bauchau OA,Cardona A.Time-step-size-independent conditioning and sensitivity to perturbations in the numerical solution of index three dif f erential algebraic equations.SIAM Journal on Scientifi Computing,2007,29(1)：397-414

32 Marsden JE,West M.Discrete mechanics and variational integrators.Acta Numerica,2001,10：357-514

HHT METHOD WITH CONSTRAINTS VIOLATION CORRECTION IN THE INDEX 2 EQUATIONS OF MOTION FOR MULTIBODY SYSTEMS1)

Ma Xiuteng∗,†,2)Zhai Yanbo∗Xie Shouyong∗

∗(College of Engineering&Technology，Southwest University，Chongqing400715，China)

†(Hubei Key Laboratory of Advanced Technology for Automotive Components，Wuhan430070，China)

Equations of motion for multibody system with holonomic constraints in Cartesian absolute coordinates modeling method are index 3 dif f erential-algebraic equations(DAEs).It is high index problem for numerical integration of index 3 DAEs.The index can be reduced to 2 by taking the derivative of position constraint equations,and velocity constraint equations can be obtained.During the integration of index 3 equations of motion,the velocity constraint equations are violated,and there are some problems in the integration of high index DAEs.Firstly,HHT(Hilber-Hughes-Taylor) direct integration method is used to the numerical integration of index 2 equations of motion.The velocity constraint equations involved in the integration,and they are satisfie in the view of computer precision.However,the positionconstraint equations are violated.Secondly,in order to eliminate the violation,the correction method based on Moore-Penrose generalized inverse theory is adopted.HHT method with constraints violation correction for index 2 equations of motion is the combination of HHT and correction method.There are no position and velocity constraints violation during the integration in the view of computer precision.No new unknown variables are introduced,and the quantity of equations in nonlinear equations from discretization is the same as index 2 equations of motion.The new integration method is validated by numerical experiments.In addition,some characteristics of HHT method,such as controlled numerical damping and second-order accuracy,are persisted by the new integration method.Finally,the quantity of nonlinear equations from discretization and computational efficiency are compared with some other methods.The advantages of the new method are illustrated.

equations of motion，index，HHT method，Moore-Penrose generalized inverse，violation correction

O313.7

A doi：10.6052/0459-1879-16-275

2016-09-29收稿，2016-11-14录用,2016-11-18网络版发表.

1)国家自然科学基金(51605391)和现代汽车零部件技术湖北省重点实验室开放基金(2012-04)资助项目.

2)马秀腾，副教授，主要研究方向：多体系统动力学建模与仿真算法.E-mail：maxt@swu.edu.cn

马秀腾,翟彦博,谢守勇.多体系统指标2运动方程HHT方法违约校正.力学学报,2017,49(1)：175-181

Ma Xiuteng,Zhai Yanbo,Xie Shouyong.HHT method with constraints violation correction in the index 2 equations of motion for multibody systems.Chinese Journal of Theoretical and Applied Mechanics,2017,49(1)：175-181