馬秀騰翟彥博謝守勇

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

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

多體系統(tǒng)指標2運動方程HHT方法違約校正1)

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

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

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

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

運動方程，指標，HHT方法，Moore--Penrose廣義逆，違約校正

引言

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

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

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

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

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

1 多體系統(tǒng)運動方程

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

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

2 指標2 運動方程的HHT方法

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

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

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

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

設由HHT方法求解方程(2)，數(shù)值積分到第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.

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

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

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

由數(shù)值實驗可知，本文方法和原HHT方法一樣，α在[-1/3,0]區(qū)間內，可以通過改變α控制方法數(shù)值阻尼的大小，隨著α的增大，數(shù)值阻尼是減少的.

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

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

圖8 精度Fig.8 Order of accuracy

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

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

5 結論

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

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

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

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

位置校正公式含有與約束Jacobian矩陣相關的m×m維矩陣求逆，對于復雜系統(tǒng)的仿真求解，如果每一步都進行校正，矩陣求逆將會耗費大量計算時間.可以根據(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)動力學微分/代數(shù)方程組數(shù)值方法.力學進展,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ù)值算法的研究進展.力學進展,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方程的最小二乘法與違約修正.力學學報,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)動力學剛性方程廣義α投影法.中國科學：物理學力學天文學.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 姚廷強,遲毅林,黃亞宇.柔性多體系統(tǒng)動力學新型廣義α數(shù)值分析方法.機械工程學報,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)動力學微分--代數(shù)方程廣義α投影法.工程力學,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方法的多體動力學數(shù)值算法研究.力學學報,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 劉穎,馬建敏.多體系統(tǒng)動力學方程的基于離散零空間理論的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 陳立平,張云清,任衛(wèi)群等.機械系統(tǒng)動力學分析及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 于清,洪嘉振.受約束多體系統(tǒng)一種新的違約校正方法.力學學報,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)絡版發(fā)表.

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

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

馬秀騰,翟彥博,謝守勇.多體系統(tǒng)指標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