多體系統接觸碰撞問題的牛頓
--歐拉線性互補方法1)

富 立 胡鴻奎 富 騰

(華北理工大學理學院，河北唐山063210)

多體系統接觸碰撞問題的牛頓
--歐拉線性互補方法1)

富 立2)胡鴻奎 富 騰

(華北理工大學理學院，河北唐山063210)

基于非光滑動力學方法的多體系統接觸碰撞分析是目前多體系統動力學的研究熱點.本文采用牛頓--歐拉方法建立多體系統接觸、碰撞問題的動力學模型，給出一種牛頓--歐拉型線性互補公式.該建模方法與目前一般采用的拉格朗日建模方法的不同之處是約束條件中除了庫侖摩擦、單邊約束之外還含有光滑等式約束.在建立系統動力學模型時，首先解除摩擦約束和單邊約束得到原系統對應的基本系統.牛頓--歐拉方法采用最大數目坐標建立基本系統的動力學方程，由于坐標不相互獨立，因此基本系統中帶有等式約束，其數學模型為一組微分代數方程.借助約束雅可比矩陣,在基本系統微分代數方程中添加摩擦接觸和單邊約束對應的拉氏乘子,就可以得到系統全局運動的具有變拓撲結構特征的動力學方程，再結合非光滑約束互補條件便可構成完備的系統動力學模型.完備的動力學模型由動力學微分方程以及等式約束和不等式約束組成.線性互補公式采用分塊矩陣形式進行推導，簡化了推導過程.數值計算采用基于線性互補的時間步進算法.時間步進算法是目前流行的非光滑數值算法，其突出特點是可以免去數值積分中繁瑣的事件檢測過程,而數值積分過程中通過對線性互補問題的求解可以確定系統的觸--離狀態.通過對典型的曲柄滑塊間隙機構進行數值分析，驗證本文方法的有效性.

非光滑，多體系統，摩擦，碰撞，線性互補問題

引言

近三十年來，非光滑多體系統動力學在基礎理論以及數值算法等方面取得顯著進展并成為力學研究的熱門領域之一.20世紀80年代，Moreau[1]將凸分析理論與測度微分包含理論相結合，應用于求解含摩擦、碰撞的剛體動力學問題，奠定了現代非光滑力學的數學和力學基礎.Panagiotopoulos[2]引入非凸變分不等式，進一步發展和完善了現代非光滑力學理論.數值算法方面，非光滑多體系統問題的關鍵是對滯-滑轉換,觸-離轉換以及碰撞等事件的檢測問題.隨著約束數目的增加，事件檢測所需的計算量呈指數速率迅速增長，這是所謂的Delassus問題[3].為克服這一難題，目前流行的做法是采用線性互補方法[4].L¨otstedt[5]首次將剛體的接觸問題歸結為一個含線性互補條件的數學規劃問題來處理,此后,經Bara ff等學者們[611]的努力，將這一理論框架成功地應用于多體系統動力學理論之中.

近幾年來，非光滑多體系統動力學有以下研究熱點及進展.在基礎理論方面,解的存在與唯一性、碰撞本構關系、動力學建模等問題的研究取得進展.文獻[12]論證了存在或不存在庫侖摩擦情況下單、雙邊約束多體系統解的存在與唯一性問題.文獻[13]對經典的Painlev′e問題給出新的分析方法和結果.文獻[14-15]研究碰撞定律的能量協調條件以及基于能量恢復系數的建模方法.數值算法方面：在時間步進算法、并行算法以及單邊約束違約修正算法等方面取得進展.文獻[16-17]分別研究了消除約束漂移效應的迭代投影時間步進算法和提高計算精度的半顯式時間步進算法.文獻[18]將Baumgarte穩定化方法推廣至含二維摩擦的多體動力學問題.文獻[19-20]給出對非光滑事件進行有效檢測的數值算法.在應用研究方面,已有的非光滑多體系統動力學理論成果被應用于航空航天、機械、機器人、車輛、生物力學、天文學等諸多領域.如小行星聚集的接觸動力學[2122]，航天器空間對接動力學[23]，間隙機構動力學[2425]，顆粒物質動力學[26]，危巖崩塌軌跡的數值模擬[27]、機器人動力學[28]、人體步態分析[29[31]、車輛非光滑動力學建模與數值模擬[32]等.

經典多體系統動力學的建模方法主要有兩類[4,2931].一類是拉格朗日方法或凱恩方法，該方法選取最少數目的獨立廣義坐標，系統運動由一組常微分方程來描述；另一類是牛頓--歐拉方法，該方法選取最大數目的非獨立廣義坐標，系統運動由一組微分代數方程來描述.P ff ei ff e等提出的非光滑多體系統動力學公式采用拉格朗日方法建立基本系統的動力學方程，其數學模型為一組常微分方程.利用單邊約束雅可比矩陣，在基本系統動力學方程中添加接觸力及碰撞力進一步可得到多體系統摩擦、碰撞問題的全局動力學方程.以上方程再結合單邊約束以及庫侖摩擦約束條件，可推導出基于常微分方程的拉格朗日型非光滑多體系統的線性互補公式.

本文采用牛頓--歐拉方法建立基本系統的動力學方程.與拉格朗日方法不同的是，解除系統單邊約束和摩擦接觸后，基本系統還存在有光滑等式約束，故基本系統的動力學方程為一組微分代數方程.基本系統微分代數方程結合單邊約束以及庫侖摩擦接觸的約束條件，可得到多體系統摩擦、碰撞問題的牛頓--歐拉型動力模型.

數值計算采用基于線性互補的時間步進算法.首先需建立系統摩擦、碰撞問題的牛頓--歐拉型線性互補公式.將系統動力學方程以及其約束互補條件離散化為差分形式，最終建立速度--沖量形式的線性互補公式.時間步進算法將線性互補解法嵌入至微分代數方程數值積分過程.算法特點是能夠求解等式約束與不等式約束并存的非光滑多體系統動力學問題，在數值計算過程中無需對滯--滑狀態等非光滑事件進行檢測，避免了由事件檢測導致的繁復計算.

本文方法是對P ff ei ff e-Glocker方法的進一步推廣.建模方面，以牛頓--歐拉方法取代拉格朗日方法，將基本系統的動力學模型由常微分方程擴展為微分代數方程.數值計算采用基于牛頓--歐拉型線性互補公式的時間步進方法.最后通過對典型機構的數值分析驗證方法的有效性.

1 多體系統接觸碰撞問題的拉格朗日方法

含摩擦、碰撞等因素的非光滑多體系統具有典型的變拓撲結構特征.在運動過程中，系統的約束以及自由度等都是隨時間變化的.系統的全局運動可分為自由運動(無單邊約束作用)，持續接觸(單邊約束持續作用)以及碰撞(單邊約束瞬間作用)等運動狀態.非光滑多體系統動力學理論采用統一的廣義坐標來描述系統的全局運動.目前普遍采用的是拉格朗日坐標.采用獨立的拉格朗日坐標描述基本系統自由運動的方法，本文稱之為拉格朗日方法.

1.1 自由運動階段

此階段多體系統沒有任何單邊約束作用，稱為基本系統.設描述基本系統運動的獨立拉格朗日坐標是q，并將其作為描述系統全局運動的統一坐標.自由運動的拉格朗日方程為

其中，M是質量矩陣，q是拉格朗日坐標列陣，h為廣義力列陣.

1.2 持續接觸階段

此階段系統受有持續的單邊約束作用，采用第一類拉氏方程建立系統的動力學方程并將單邊約束表述為互補形式

其中，gN是接觸處的法向距離列陣、γT是接觸處的切向相對速度列陣，γR和γL是將γT分解后的互補速度，滿足關系γT=γR?γL.λN和λT分別是接觸處的法向反力列陣和切向摩擦力列陣.λR和λL是將λT分解后的互補摩擦力列陣，滿足關系λR=2μλN?λL(μ為摩擦系數矩陣).WN和WT分別是法向和切向的約束雅可比矩陣.在運動過程中，系統的接觸狀態隨時間發生變化，因此方程中的系數矩陣WN和WT以及約束反力列陣λN和λT的維數都是隨時間變化的.

1.3 碰撞階段

設碰撞始于t?終于t+，剛體動力學方法假設碰撞在瞬間完成(?t=t+?t?→0)，速度產生突變.

碰撞方程以及碰撞互補關系為

其中，ξN和ξT分別是法向和切向的碰撞修正速度，，；ξL和 ξR是將切向修正速度ξT分解后的互補速度，滿足關系分別是接觸處的法向反力沖量和切向摩擦力沖量；ΛR和ΛL是將ΛT分解后的互補摩擦沖量列陣，滿足關系ΛR=2μΛN?ΛL；eN和eT分別是法向、切向碰撞恢復系數對角陣.

1.4 基于線性互補的時間步進算法

在多點接觸碰撞問題中，當一點接觸狀態改變時，如觸離狀態的改變、滯滑狀態的改變或發生碰撞等，都會影響到其他各點的接觸狀態.在動力學計算中需要對系統的接觸狀態進行判斷.當接觸點增多時，這個判斷過程是非常繁復的.現代非光滑力學對剛體系統接觸狀態的判斷建立了一種嚴格而有效的方法，即線性互補方法.該方法將系統動力學方程及其單邊約束條件轉化為一個線性互補問題，通過求解線性互補問題來判斷系統的接觸狀態，具體參見文獻[1,4,7-8]等.

時間步進算法是適于非光滑多體系統的數值積分方法.眾多文獻所采用的時間步進方法中，Moreau的中點算法是最經典也是使用最廣泛的.時間步進方法將動力學方程和單邊約束條件離散化，離散化的差分方程不以加速度和力為未知量，而是以速度和沖量取而代之.時間步進算法的突出特點是數值積分過程中無需對非光滑事件進行檢測，是真正的非光滑算法.該算法詳細內容及主要步驟參見文獻[1,11].

2 多體系統接觸碰撞問題的牛頓--歐拉方法

2.1 自由運動階段的牛頓--歐拉方程

本文采用絕對笛卡兒坐標描述基本系統的自由運動，稱之為牛頓--歐拉方法.自由運動階段基本系統沒有單邊約束只有光滑雙邊約束，將光滑雙邊約束的約束反力添加至系統動力學方程(1)中，可以建立用絕對笛卡兒坐標表示的牛頓--歐拉型基本系統動力學方程如下

其中，gE是基本系統的光滑雙邊約束列陣.WE=(?gE/?q)T為雙邊約束的雅可比矩陣的轉置矩陣，λE為雙邊約束的約束力列陣.

2.2 牛頓--歐拉接觸動力學方程

持續接觸階段，在方程(2)中添加雙邊約束的約束反力，并在約束條件(7)～(9)基礎上補充光滑雙邊約束條件gE=0，得到以下持續接觸階段的牛頓--歐拉型系統動力學方程及其約束條件

建立接觸階段線性互補公式時，需要將方程及約束條件(12)～(16)離散化.離散化后的速度--沖量形式的接觸動力學模型是

2.3 牛頓--歐拉碰撞動力學方程

在碰撞動力學方程(6)中添加光滑雙邊約束的約束反力沖量并在碰撞約束條件(7)～(9)中補充光滑雙邊約束條件，得以下速度--沖量形式的牛頓--歐拉型碰撞動力學模型

2.4 牛頓--歐拉線性互補公式

接觸線性互補公式與碰撞線性互補公式的推導公式類似，以下給出牛頓--歐拉型碰撞線性互補公式.

將式(24)～式(30)，寫為矩陣形式

將式(31)寫成以下分塊矩陣形式

將式(32)展開

由式(33)

代入式(34)，得

將式(36)代入式(35)，得

將式(37)代入式(38)，得

經整理簡化，式(39)簡寫為

其中

由于x和y之間滿足互補關系，因此式(22)最終可化為以下標準線性互補公式

令公式中的恢復系數eN=eT=0，即變為持續接觸問題的線性互補公式.

3 算例

含間隙滑移鉸曲柄滑塊機構如圖1所示.其中含間隙的滑移鉸如圖2所示.

圖1 含間隙滑移鉸的曲柄滑塊機構Fig.1 Slider-crank mechanism with a translational clearance joint

圖2 由滑塊與滑道軌構成的含間隙滑移鉸Fig.2 Translational joint with clearance that is,the slider and guide

圖2 由滑塊與滑道軌構成的含間隙滑移鉸(續)Fig.2 Translational joint with clearance that is,the slider and guide(continued)

機構幾何參數：

慣性參數：

接觸參數：

初始條件：

為使問題簡化，曲柄和連桿的位形仍然用最少數目坐標，即θ1和θ2，來描述，只有滑塊的位形采用最大數目坐標(x3y3θ3)來描述，因此系統的廣義坐標為

基本系統動力學方程及其約束

其中

法向和切向約束矩陣為

根據系統的接觸狀態，WN和WT分別在以下WN1～WN4和WT1～WT4之間取舍組成

法向約束力和切向摩擦力列陣為

系統全局動力學方程為

圖3為曲柄轉動兩周過程中曲柄角速度、連桿角速度、連桿的角位移和角速度相圖以及滑塊質心的軌跡圖.

圖3 Fig.3

圖3 (續)Fig.3(continued)

圖4為曲柄轉動兩周過程中四個角點的運動.從運動曲線中可以觀察到滑塊的自由運動以及接觸碰撞事件的發生.

圖4 滑塊角點的運動Fig.4 Motion of the slider corners

圖5顯示恢復系數對運動的影響，恢復系數分別取為0.1,0.4,0.7,0.9.

圖6顯示間隙尺寸對系統響應的影響.間隙尺寸c分別取為0.1，0.2，0.5和1.0mm.

圖5 恢復系數對滑塊角點1運動的影響Fig.5 In fl uence of the restitution coefficient on the motion of corner 1

圖6 間隙尺寸對滑塊角點1運動的影響Fig.6 In fl uence of the clearance size on the motion of corner 1

圖6 間隙尺寸對滑塊角點1運動的影響(續)Fig.6 In fl uence of the clearance size on the motion of corner 1(continued)

4 結論

本文將多體系統的線性互補方法應用于含有光滑等式約束的摩擦、碰撞多體系統.采用牛頓--歐拉方法 (即最大數目坐標法)建立基本系統動力學方程，由于含有等式約束，因此其數學模型是一組微分代數方程.在基本系統上，考慮接觸點處的法向和切向互補關系，并通過約束矩陣，將對應的法向、切向拉氏乘子(約束反力)添加加至基本系統微分代數方程中，便構成接觸狀態下的系統動力學模型.同樣地，在基本系統上，考慮碰撞點處的法向和切向互補關系，并通過約束矩陣將對應的法向、切向拉氏乘子(約束反力)加至基本系統微分代數方程中，便構成碰撞時的系統動力學模型.

將單邊約束、庫侖摩擦定律和牛頓碰撞定律均轉化為互補關系來表述，將動力學方程及其所有約束互補關系離散化為速度--沖量形式，采用分塊矩陣模型簡化推導過程，給出了一種新型線性互補公式.僅考慮單邊約束的多體系統動力學問題是常微分方程與線性互補的混合動力學問題，同時含等式約束及不等式約束的多體系統動力學問題是微分代數方程與線性互補的混合動力學問題.數值計算采用時間步進算法，將線性互補解法嵌入至Moreau的時間步進算法之中.數值計算過程中，系統動力學模型將在自由運動、持續接觸運動以及碰撞三種運動狀態下相互切換，持續接觸過程中接觸狀態也隨時間變化，體現了典型的變拓撲結構特征.

最后，通過對含間隙滑移鉸曲柄滑塊機構的建模與數值分析，驗證了本文方法的有效性.

1 Moreau JJ.Unilateral contact and dry friction in fi nite freedom dynamics.In:Moreau JJ,Panagiotopoulos PD.Non-smooth Mechanics and Applications//International Centre for Mechanical Sciences,Courses and Lectures,Vol.302,New York,Springer-Verlag,1988:1-82

2 Panagiotopoulos PD.Inequality Problems in Mechanics and Applications.Boston:Stuttgart,Birkh′auser,1985

3 Brogliato B.Nonsmooth Mechanics:Models,Dynamics and Control,2nd ed.London:Springer-Verlag,1999

4富立.非光滑多體系統動力學線性互補方法.北京：清華大學出版社,2016(Fu Li.LCP Method for Non-smooth Multibody System Dynamics.Beijing:Tsinghua University Press,2016(in Chinese))

5 L¨otstedt P.Mechanical systems of rigid bodies subject to unilateral constraints.Siam Journal on Applied Mathematics,1982,42(2):281-296

6 Bara ff D.Issuesincomputingcontactforcesfornonpenetratingrigid bodies.Algorithmica,1993,10:292-352

7 Pfei ff er F,Glocker C.Multibody Dynamics with Unilateral Contacts//Non-linear Dynamics.Weinheim:John Wiley&Sons,1996

8 Glocker,C,Set-Valued Force Laws-Dynamics of Non-Smooth Systems.Berlin:Springer,2001

9 Stewart DE.Rigid body dynamics with friction and impact.Siam Review,2000,42(1):3-39

10 Anitescu M,Potra FA.Formulationg dynamics multi-rigid-body contact problems with friction as solvable linear complementarity problems.Nonlinear Dynamics,1997,14:231-247

11 Acary V,Brogliato B.Numerical Methods for Nonsmooth Dynamical Systems//Applications in Mechanics and Electronics.Berlin:Springer-Verlag,2008

12 Blumentals A,Brogliato B.The contact problem in Lagrangian sys-tems subject to bilateral and unilateral constraints,with or without sliding Coulomb’s friction:A tutorial.Multibody System Dynamics,2016,1:1-34

13 Zhao Z,Liu CS,Chen B,et al.Asymptotic analysis of Painleve¨s paradox.Multibody System Dynamics，2015,35(3):299-319

14 Glocker Ch.Energetic consistency conditions for standard impacts.In:Part II:Poisson-type inequality impact laws.Multibody System Dynamics,2014,32(1):445-509

15 Dietmayer K.Modelling of unilateral constraints using power-based restriction functions within Lagrangian mechanics.Mathematical&Computer Modelling of Dynamical System,2015,21(3):1-26

16 Glocker C.Simulation of Hard Contacts with Friction.In:An Iterative Projection Method.Recent Trends in Dynamical Systems,Springer,Basel,Switzerland,2013:493-515

17 Schindler T,Rezaei S,Kursawe J,et al.Half-explicit timestepping schemes on velocity level based on time-discontinuous Galerkin methods.Computer Methods in Applied Mechanics and Engineering,2015,290(15):250-276

18 Kikuuwe R,Brogliato B.A new representation of systems with frictional unilateral constraints and its Baumgarte-like relaxation.Multibody System Dynamics,2015,34(7):1-24

19 Pournaras A,Karaoulanis F,Natsiavas S.Dynamics of mechanical systems involving impact and friction using an effi cient contact detection algorithm. http://dx.doi.org/10.1016/j.ijnonlinmec.2016.08.007,2016-8-24

20王曉軍,王琪.含摩擦與碰撞平面多剛體系統動力學線性互補算法.力學學報,2015,47(2):814-821(Wang Xiaojun,Wang Qi.A LCP method for the dynamics of planar multibody systems with impact and friction.Chinese Journal of Theoretical and Applied Mechanics,2015,47(2):814-821(in Chinese))

21 Ferrari F，Tasora A，Masarati P,et al.N-body gravitational and contact dynamics for asteroid aggregation.Multibody System Dynamics，2017,39(1):3-20

22張韻，李俊峰.碎石堆小行星的散體動力學建模與仿真方法綜述.力學學報,2015,47(1):1-7(Zhangyun,Li Junfeng.A survey of granular dynamics modeling and simulation methods for rubble-pile asteroids.Chinese Journal of Theoretical and Applied Mechanics,2015,47(1):1-7(in Chinese))

23 Zhao Zhen，Liu Caishan，Chen Tao.Docking dynamics between two spacecrafts with APDSes.Multibody System Dynamics，2016,37(3):245-270

24 Yaqubi S,Dardel M.Daniali HM.Modeling and control of crank–slider mechanism with multiple clearance joints.Multibody System Dynamics,2016,36(2):1-25

25 Abdallah MAB,Khemili I,Aifaoui N.Numerical investigation of a fl exible slider–crank mechanism with multijoints with clearance.Multibody System Dynamics,2016,38(2):173-199

26 LimKW,KrabbenhoftK,Jos′eE,etal.Acontactdynamicsapproach to the Granular Element Method.Computer Methods in Applied Mechanics and Engineering,2014,268(1):557-573

27 Leine RI,Schweizer A,Christen M,et al.Simulation of rockfall trajectories with consideration of rock shapein.Multibody System Dynamics,2014,32:241-271

28 Zobova AZ,Nicolas TH,Noot V,et al.Multi-physics modelling of a compliant humanoid robot.Multibody System Dynamics,2017,39(1):95-114

29 Shourijeh MS，Mcphee J.Foot–ground contact modeling within human gait simulations:From Kelvin–Voigt to hyper-volumetric models.Multibody System Dynamics,2015,35(1):393-407

30 Josep JR,Font-Llagunes M,Plaza A,et al.Dynamic considerations of heel-strike impact in human gait.Multibody System Dynamics,2015,35(3):215-232

31洪嘉振.計算多體系統動力學.北京:高等教育出版社,1999(Hong Jiazhen.Computational Dynamics of Multibody Systems.Beijing:Higher Education Press,1999(in Chinese))

32齊朝暉.多體系統動力學.大連:科學出版社,2008(Qi zhaohui.Multi-body System Dynamics.Dalian:Science Press,2008(in Chinese))

33 Shabana AA.Computational Dynamics 3rd ed.West Sussex:John Wiley,2010

34范新秀,王琪.車輛縱向非光滑多體動力學建模與數值算法研究.力學學報,2015,47(2):301-309(Fan Xinxiu,Wang Qi.Research on modeling and simulation of longitudinal vehicle dynamics based on non-smooth dynamics of multibody systems.Journal of Theoretical and Applied Mechanics,2015,47(2):301-309(in Chinese))

CONTACT-IMPACT ANALYSIS IN MULTI-BODY SYSTEMS BASED ON NEWTON-EULER LCP APPROACH1)

Fu Li2)Hu Hongkui Fu Teng
(Science School，North China University of Science and Technology，Tangshan 063210，Hebei，China)

The contact-impact analysis in multibody systems based on the nonsmooth dynamics approach is a hot topic in the research of multibody system dynamics.Newton-Euler approach is adopted to develop dynamics model of contactimpact analysis in non-smooth multi-body systems,and a new LCP formula is presented in this work.Di ff erent from Lagrange methods，Newton-Euler modeling method incorporate equality constraints into dynamic models with noninterpenetration constraints and frictional constraints together.In Newton-Euler modeling method,the basic system is derived by removing the non-interpenetration constraints and frictional constraints from the original multi-body system.Newton-Euler eqution of basic system is established by using the maximum coordinates method.Because the coordinates of the basic system are not independent of each other,equality constraints are involved in modeling,the basic system dynamic equations is a set of DAE(di ff erential algebra equation).With the aid of constraint Jacobian matrix,Lagrangian multipliers corresponding to the non-interpenetration constraint forces and Coulomb friction forces are added to the basic system DAE to obtain the dynamic equations of global motion of the multi-body system with characteristics of variable topologicalstructure.ThecompletedynamicmodeliscomposedofbasicsystemDAE,equalityandinequalityconstraints.In order to simplify the derivation process of LCP,a decomposed matrix form is built.The LCP-based Time-stepping method is adopted for numerical simulation.Time-stepping algorithm is a popular non-smooth numerical algorithm,Its prominent feature is that it can avoid the tedious event-detection process in numerical integration.In the process of numerical integration,the contact-detachment state of the system can be determined by solving the LCP.Our method is carried out in slider-crank mechanism with a translational clearance joint,the simulation results indicate that this method is e ff ective.

non-smooth,multi-body system,friction,collision,LCP

O313.7

A

10.6052/0459-1879-17-023

2017–01–18收稿，2017–06–07 錄用,2017–06–20 網絡版發表.

1)河北省自然科學基金資助項目(A2013209221).

2)富立，教授，主要研究方向：多體系統動力學及非線性動力學.E-mail：13231554976@163.com

富立,胡鴻奎,富騰.多體系統接觸碰撞問題的牛頓--歐拉線性互補方法.力學學報,2017,49(5):1115-1125

Fu Li,Hu Hongkui,Fu Teng.Contact-impact analysis in multi-body systems based on Newton-Euler LCP approach.Chinese Journal of Theoretical and Applied Mechanics,2017,49(5):1115-1125