馬秀騰翟彥博謝守勇

?(西南大學(xué)工程技術(shù)學(xué)院，重慶400715)

?(現(xiàn)代汽車零部件技術(shù)湖北省重點(diǎn)實(shí)驗(yàn)室，武漢430070)

多體系統(tǒng)指標(biāo)2運(yùn)動(dòng)方程HHT方法違約校正1)

馬秀騰?,?,2)翟彥博?謝守勇?

?(西南大學(xué)工程技術(shù)學(xué)院，重慶400715)

?(現(xiàn)代汽車零部件技術(shù)湖北省重點(diǎn)實(shí)驗(yàn)室，武漢430070)

采用Cartesian絕對(duì)坐標(biāo)建模方法，完整約束多體系統(tǒng)運(yùn)動(dòng)方程是指標(biāo)3的微分--代數(shù)方程(dif f erentialalgebraic equations,DAEs)，數(shù)值求解指標(biāo)3的DAEs屬于高指標(biāo)問題，通過對(duì)位置約束方程求導(dǎo)，可使運(yùn)動(dòng)方程的指標(biāo)降為2.位置約束方程求導(dǎo)得到的是速度約束方程.直接求解指標(biāo)3的運(yùn)動(dòng)方程，速度約束方程得不到滿足，而且高指標(biāo)DAEs的數(shù)值求解存在一些問題.論文首先采用HHT(Hilber--Hughes--Taylor)直接積分方法求解降指標(biāo)得到的指標(biāo)2運(yùn)動(dòng)方程，此時(shí)速度約束方程參與離散計(jì)算，從機(jī)器精度上講速度約束自然得到滿足，而位置約束方程沒有參與計(jì)算，存在“違約”.針對(duì)違約問題，采用基于Moore--Penrose廣義逆理論的違約校正方法，消除位置約束方程的違約.指標(biāo)2運(yùn)動(dòng)方程HHT方法違約校正，將HHT方法和違約校正方法很好地結(jié)合，在數(shù)值求解指標(biāo)2運(yùn)動(dòng)方程的過程中，位置約束方程和速度約束方程都不存在違約問題，而且新方法沒有引入新的未知數(shù)向量，離散得到的非線性方程組的方程數(shù)量與原指標(biāo)2運(yùn)動(dòng)方程的方程數(shù)量相同，求解規(guī)模沒有擴(kuò)大.新方法的實(shí)用和有效性通過算例的數(shù)值實(shí)驗(yàn)得到驗(yàn)證，數(shù)值實(shí)驗(yàn)也說明新方法保持了HHT方法本身具有的數(shù)值阻尼可以控制和二階精度的特性.最后從非線性方程組的求解規(guī)模和計(jì)算速度上與其他方法進(jìn)行了比較分析，說明新方法的優(yōu)勢所在.

運(yùn)動(dòng)方程，指標(biāo)，HHT方法，Moore--Penrose廣義逆，違約校正

引言

為了實(shí)現(xiàn)多體系統(tǒng)程式化建模與仿真求解，采用Cartesian絕對(duì)坐標(biāo)建模方法，對(duì)于完整約束多體系統(tǒng)，此時(shí)得到的運(yùn)動(dòng)方程是指標(biāo)3的微分--代數(shù)方程 (dif f erential-algebraic equations,DAEs)[1].由于DAEs和常微分方程(ODEs)有本質(zhì)區(qū)別，且一般沒有解析解，因此要研究專門針對(duì)DAEs的數(shù)值求解算法[2-7].

已有的運(yùn)動(dòng)方程求解方法可分為增廣法和縮并法[2]，也可分為直接法和狀態(tài)空間法[8].在直接法中，結(jié)構(gòu)動(dòng)力學(xué)中的直接積分方法，如Newmark方法、HHT(Hilber-Hughes-Taylor)方法、廣義α方法、θ1方法等，有廣泛的應(yīng)用[8-17].也有狀態(tài)空間法和廣義α方法相結(jié)合的數(shù)值求解算法[11].在直接法中，HHT方法和Newmark方法已集成到MSC.ADAMS軟件的求解器上[9,18].

對(duì)于指標(biāo) 3運(yùn)動(dòng)方程，對(duì)其位置約束方程求導(dǎo)，整體將得到指標(biāo)2 DAEs形式的運(yùn)動(dòng)方程.位置約束方程求導(dǎo)得到的方程稱為速度約束方程.直接數(shù)值求解指標(biāo)3運(yùn)動(dòng)方程，屬于高指標(biāo)問題[1]，求解存在一些問題，且速度約束方程存在違約問題.直接求解指標(biāo)2運(yùn)動(dòng)方程，位置約束方程違約存在問題.如何同時(shí)滿足位置約束和速度約束方程，一種求解思路是基于GGL約束穩(wěn)定[19]的方法，通過引入新的“未知數(shù)向量”，將速度約束方程并到指標(biāo)3運(yùn)動(dòng)方程中一起，求解約束穩(wěn)定指標(biāo)2運(yùn)動(dòng)方程.從機(jī)器精度上講，位置約束和速度約束方程同時(shí)得到滿足，已有如HHT-SI2方法[8]、HHT-ADD方法[8]、投影方法[10,12]、廣義α-SOI2方法[15]、θ1方法[13]和Lie群積分[20-21]等.GGL約束穩(wěn)定指標(biāo)2運(yùn)動(dòng)方程的求解形式已經(jīng)集成到MSC.ADAMS求解器中，可以選用Gsti ff,Wstiff,Constant_BDF,Newmark和HHT等方法求解[22].MSC.ADAMS求解器中通過引入兩個(gè)“導(dǎo)數(shù)向量”[1,22]來代替指標(biāo)3運(yùn)動(dòng)方程中的Lagrange乘子和GGL約束穩(wěn)定引入的“未知數(shù)向量”，得到約束穩(wěn)定指標(biāo)1運(yùn)動(dòng)方程.不論是約束穩(wěn)定指標(biāo)2還是約束穩(wěn)定指標(biāo)1運(yùn)動(dòng)方程，由于引入了新的未知數(shù)，離散后得到的非線性方程組的求解規(guī)模將變大.

對(duì)指標(biāo)2運(yùn)動(dòng)方程中速度約束方程求導(dǎo)，整體將得到指標(biāo)1 DAEs形式的運(yùn)動(dòng)方程，MATLAB軟件的ode15s和ode23t函數(shù)可以求解DAEs和剛性O(shè)DEs，但僅能求解指標(biāo)1的DAEs[23].求解多體系統(tǒng)指標(biāo)1的運(yùn)動(dòng)方程會(huì)存在位置約束和速度約束方程存在違約問題.Yoon等[24]、于清等[25-26]從改進(jìn)Baumgarte約束違約穩(wěn)定法角度，Nikravesh[27]從初始條件校正角度，分別基于Moore-Penrose廣義逆理論，提出新的違約校正方法.最近Marques等[28-29]在此方法上又有深入的研究.

Gear[30]證明對(duì)于力學(xué)系統(tǒng)，數(shù)值求解指標(biāo)2運(yùn)動(dòng)方程是較好的選擇.本文將直接積分中的HHT方法和基于Moore-Penrose廣義逆的違約校正方法相結(jié)合.用HHT方法求解指標(biāo)2 DAEs形式的運(yùn)動(dòng)方程，位置約束方程的違約通過校正方法消除.新方法僅在求解的過程中增加一步迭代，沒有引入新的變量.較之已有的方法，求解速度將得到提高.

1 多體系統(tǒng)運(yùn)動(dòng)方程

完整約束多體系統(tǒng)運(yùn)動(dòng)方程形式如下

將位置約束方程Φ(q)=0求導(dǎo)，方程(1)的指標(biāo)由3降為2，方程形式如下

2 指標(biāo)2 運(yùn)動(dòng)方程的HHT方法

對(duì)于方程(2)，采用HHT方法離散，并對(duì)速度約束方程進(jìn)行縮放[31]，得到的非線性方程組如下

非線性方程組(3)中方程的數(shù)量為n+m，采用Newton-Raphson方法求解此方程組，得

此時(shí)數(shù)值解使速度約束方程得到滿足，而位置約束方程存在違約問題.也就是說從機(jī)器精度上講

3 位置約束方程的違約校正

設(shè)由HHT方法求解方程(2)，數(shù)值積分到第n+1步后，得到的廣義坐標(biāo)為qn+1，此時(shí)位置約束方程存在違約問題.可以在廣義坐標(biāo)qn+1上加入校正項(xiàng)得

位置約束方程Φ(q)=0的變分可表示為

4 算例

如圖1所示的曲柄滑塊機(jī)構(gòu)，采用Cartesian絕對(duì)坐標(biāo)建模方法，設(shè)曲柄AB長為l1，質(zhì)心坐標(biāo)為(x1,y1)，連桿BC長為l2，質(zhì)心坐標(biāo)為(x2,y2)，滑塊質(zhì)心坐標(biāo)為(x3,0).

圖1 曲柄滑塊機(jī)構(gòu)簡圖Fig.1 Sketch of slider crank mechanism

廣義坐標(biāo)為

位置約束方程為

取曲柄質(zhì)量m1=3kg，長度l1=0.6m，連桿質(zhì)量m2=0.9kg，長度l2=1.2m，滑塊質(zhì)量m3=1.8kg.重力加速度g=9.81m/s2.彈簧剛度k=100N/m.

取阻尼系數(shù)c=1N·s/m，設(shè)求解時(shí)間20s，步長h=1ms.

圖2 機(jī)構(gòu)θ角時(shí)間歷程Fig.2 Time evolution of θ for mechanism

圖2給出本文的HHT方法和違約校正相結(jié)合(HHT-I2-Correction)方法與HHT-SI2方法[8]求解結(jié)果的比較，完全吻合，說明方法是有效的.

圖3給出HHT方法直接求解運(yùn)動(dòng)方程(2)時(shí)速度約束違約量的曲線.

圖4給出HHT方法直接求解運(yùn)動(dòng)方程(2)時(shí)位置約束違約量的曲線.

圖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可知，速度約束和位置約束方程的違約數(shù)量級(jí)都在10-16～10-14之間，從機(jī)器精度上講，不存在違約問題.

當(dāng)曲柄滑塊機(jī)構(gòu)中阻尼系數(shù)c=0時(shí)，設(shè)求解時(shí)間20s，步長h=1ms，圖7給出當(dāng)α=-1/3,-1/6, -1/48時(shí)，數(shù)值求解過程中系統(tǒng)能量的變化曲線.

圖7 求解過程中系統(tǒng)能量的變化曲線Fig.7 Energy of the system during integration

由數(shù)值實(shí)驗(yàn)可知，本文方法和原HHT方法一樣，α在[-1/3,0]區(qū)間內(nèi)，可以通過改變?chǔ)量刂品椒〝?shù)值阻尼的大小，隨著α的增大，數(shù)值阻尼是減少的.

當(dāng)α=0，由HHT方法得γ=0.5，β=0.25，這就是梯形法則公式[32].如果α=0，γ=0.5，令β=0，這就是Strmer算法[32]，這兩種方法都是沒有數(shù)值阻尼的.

圖8給出通過數(shù)值實(shí)驗(yàn)得到的本文方法精度的階，和原HHT方法一樣都具有二階精度.

圖8 精度Fig.8 Order of accuracy

由文獻(xiàn)[14]的研究可知，為了使求解過程中位置約束和速度約束方程同時(shí)得到滿足，在求解指標(biāo)2超定微分--代數(shù)方程時(shí)，步長取為h=1ms，廣義α-SOI2方法、HHT-SI2方法、θ1方法和文獻(xiàn)[14]中新廣義α方法的求解速度是相當(dāng)?shù)模?方法稍快.前3種方法非線性方程組的求解規(guī)模為2n+2m，新廣義α方法的求解規(guī)模為n+2m.

對(duì)于圖1所示的模型，阻尼系數(shù)c=1N·s/m，設(shè)求解時(shí)間為20s，步長為h=1ms.相同的編程語言編寫程序，在相同配置的電腦上運(yùn)行，θ1方法的CPU時(shí)間為20.46s，本文HHT直接積分校正方法僅為17.38s.可見方法具有求解速度的優(yōu)勢，而且求解的非線性方程組的規(guī)模是最小的，僅為n+m.

5 結(jié)論

求解多體系統(tǒng)運(yùn)動(dòng)方程時(shí)，為了得到更精確的解，要求位置約束和速度約束方程都不存在違約問題.

本文用HHT方法直接求解指標(biāo)2 DAEs形式的運(yùn)動(dòng)方程，此時(shí)速度約束方程自然滿足，為了消除位置約束方程的違約，引入基于Moore-Penrose廣義逆理論的違約校正方法，使得位置約束方程得到滿足.

數(shù)值實(shí)驗(yàn)表明本文方法保持了原HHT方法可控?cái)?shù)值阻尼和二階精度的特性，且求解速度比廣義α-SOI2方法、HHT-SI2方法、θ1方法和新廣義α方法都快.

從形式上看，本文方法非線性方程組求解規(guī)模最小.對(duì)于復(fù)雜多體系統(tǒng)建模與仿真，廣義坐標(biāo)數(shù)n和約束方程數(shù)m將相應(yīng)增大，求解規(guī)模小的特點(diǎn)將更具優(yōu)勢.

位置校正公式含有與約束Jacobian矩陣相關(guān)的m×m維矩陣求逆，對(duì)于復(fù)雜系統(tǒng)的仿真求解，如果每一步都進(jìn)行校正，矩陣求逆將會(huì)耗費(fèi)大量計(jì)算時(shí)間.可以根據(jù)計(jì)算精度的要求，設(shè)定當(dāng)位置違約量大于某一量值后才進(jìn)行校正，以提高計(jì)算速度.

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 潘振寬,趙維加,洪嘉振等.多體系統(tǒng)動(dòng)力學(xué)微分/代數(shù)方程組數(shù)值方法.力學(xué)進(jìn)展,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 王琪,陸啟韶.多體系統(tǒng)Lagrange方程數(shù)值算法的研究進(jìn)展.力學(xué)進(jìn)展,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 趙維加,潘振寬.多體系統(tǒng)Euler-Lagrange方程的最小二乘法與違約修正.力學(xué)學(xué)報(bào),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 丁潔玉,潘振寬.多體系統(tǒng)動(dòng)力學(xué)剛性方程廣義α投影法.中國科學(xué)：物理學(xué)力學(xué)天文學(xué).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 姚廷強(qiáng),遲毅林,黃亞宇.柔性多體系統(tǒng)動(dòng)力學(xué)新型廣義α數(shù)值分析方法.機(jī)械工程學(xué)報(bào),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 丁潔玉,潘振寬.多體系統(tǒng)動(dòng)力學(xué)微分--代數(shù)方程廣義α投影法.工程力學(xué),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 馬秀騰,翟彥博,羅書強(qiáng).基于θ1方法的多體動(dòng)力學(xué)數(shù)值算法研究.力學(xué)學(xué)報(bào),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 馬秀騰,翟彥博,羅書強(qiáng).多體動(dòng)力學(xué)超定運(yùn)動(dòng)方程廣義α求解新算法.力學(xué)學(xué)報(bào),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 劉穎,馬建敏.多體系統(tǒng)動(dòng)力學(xué)方程的基于離散零空間理論的Newmark積分算法.機(jī)械工程學(xué)報(bào),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 陳立平,張?jiān)魄?任衛(wèi)群等.機(jī)械系統(tǒng)動(dòng)力學(xué)分析及ADAMS應(yīng)用教程.北京：清華大學(xué)出版社,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 于清,洪嘉振.受約束多體系統(tǒng)一種新的違約校正方法.力學(xué)學(xué)報(bào),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網(wǎng)絡(luò)版發(fā)表.

1)國家自然科學(xué)基金(51605391)和現(xiàn)代汽車零部件技術(shù)湖北省重點(diǎn)實(shí)驗(yàn)室開放基金(2012-04)資助項(xiàng)目.

2)馬秀騰，副教授，主要研究方向：多體系統(tǒng)動(dòng)力學(xué)建模與仿真算法.E-mail：maxt@swu.edu.cn

馬秀騰,翟彥博,謝守勇.多體系統(tǒng)指標(biāo)2運(yùn)動(dòng)方程HHT方法違約校正.力學(xué)學(xué)報(bào),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