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

富 立 胡鴻奎 富 騰

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

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

富 立2)胡鴻奎 富 騰

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

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

非光滑，多體系統(tǒng)，摩擦，碰撞，線性互補(bǔ)問題

引言

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

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

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

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

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

本文方法是對P ff ei ff e-Glocker方法的進(jìn)一步推廣.建模方面，以牛頓--歐拉方法取代拉格朗日方法，將基本系統(tǒng)的動力學(xué)模型由常微分方程擴(kuò)展為微分代數(shù)方程.數(shù)值計(jì)算采用基于牛頓--歐拉型線性互補(bǔ)公式的時(shí)間步進(jìn)方法.最后通過對典型機(jī)構(gòu)的數(shù)值分析驗(yàn)證方法的有效性.

1 多體系統(tǒng)接觸碰撞問題的拉格朗日方法

含摩擦、碰撞等因素的非光滑多體系統(tǒng)具有典型的變拓?fù)浣Y(jié)構(gòu)特征.在運(yùn)動過程中，系統(tǒng)的約束以及自由度等都是隨時(shí)間變化的.系統(tǒng)的全局運(yùn)動可分為自由運(yùn)動(無單邊約束作用)，持續(xù)接觸(單邊約束持續(xù)作用)以及碰撞(單邊約束瞬間作用)等運(yùn)動狀態(tài).非光滑多體系統(tǒng)動力學(xué)理論采用統(tǒng)一的廣義坐標(biāo)來描述系統(tǒng)的全局運(yùn)動.目前普遍采用的是拉格朗日坐標(biāo).采用獨(dú)立的拉格朗日坐標(biāo)描述基本系統(tǒng)自由運(yùn)動的方法，本文稱之為拉格朗日方法.

1.1 自由運(yùn)動階段

此階段多體系統(tǒng)沒有任何單邊約束作用，稱為基本系統(tǒng).設(shè)描述基本系統(tǒng)運(yùn)動的獨(dú)立拉格朗日坐標(biāo)是q，并將其作為描述系統(tǒng)全局運(yùn)動的統(tǒng)一坐標(biāo).自由運(yùn)動的拉格朗日方程為

其中，M是質(zhì)量矩陣，q是拉格朗日坐標(biāo)列陣，h為廣義力列陣.

1.2 持續(xù)接觸階段

此階段系統(tǒng)受有持續(xù)的單邊約束作用，采用第一類拉氏方程建立系統(tǒng)的動力學(xué)方程并將單邊約束表述為互補(bǔ)形式

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

1.3 碰撞階段

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

碰撞方程以及碰撞互補(bǔ)關(guān)系為

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

1.4 基于線性互補(bǔ)的時(shí)間步進(jìn)算法

在多點(diǎn)接觸碰撞問題中，當(dāng)一點(diǎn)接觸狀態(tài)改變時(shí)，如觸離狀態(tài)的改變、滯滑狀態(tài)的改變或發(fā)生碰撞等，都會影響到其他各點(diǎn)的接觸狀態(tài).在動力學(xué)計(jì)算中需要對系統(tǒng)的接觸狀態(tài)進(jìn)行判斷.當(dāng)接觸點(diǎn)增多時(shí)，這個(gè)判斷過程是非常繁復(fù)的.現(xiàn)代非光滑力學(xué)對剛體系統(tǒng)接觸狀態(tài)的判斷建立了一種嚴(yán)格而有效的方法，即線性互補(bǔ)方法.該方法將系統(tǒng)動力學(xué)方程及其單邊約束條件轉(zhuǎn)化為一個(gè)線性互補(bǔ)問題，通過求解線性互補(bǔ)問題來判斷系統(tǒng)的接觸狀態(tài)，具體參見文獻(xiàn)[1,4,7-8]等.

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

2 多體系統(tǒng)接觸碰撞問題的牛頓--歐拉方法

2.1 自由運(yùn)動階段的牛頓--歐拉方程

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

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

2.2 牛頓--歐拉接觸動力學(xué)方程

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

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

2.3 牛頓--歐拉碰撞動力學(xué)方程

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

2.4 牛頓--歐拉線性互補(bǔ)公式

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

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

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

將式(32)展開

由式(33)

代入式(34)，得

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

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

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

其中

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

令公式中的恢復(fù)系數(shù)eN=eT=0，即變?yōu)槌掷m(xù)接觸問題的線性互補(bǔ)公式.

3 算例

含間隙滑移鉸曲柄滑塊機(jī)構(gòu)如圖1所示.其中含間隙的滑移鉸如圖2所示.

圖1 含間隙滑移鉸的曲柄滑塊機(jī)構(gòu)Fig.1 Slider-crank mechanism with a translational clearance joint

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

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

機(jī)構(gòu)幾何參數(shù)：

慣性參數(shù)：

接觸參數(shù)：

初始條件：

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

基本系統(tǒng)動力學(xué)方程及其約束

其中

法向和切向約束矩陣為

根據(jù)系統(tǒng)的接觸狀態(tài)，WN和WT分別在以下WN1～WN4和WT1～WT4之間取舍組成

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

系統(tǒng)全局動力學(xué)方程為

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

圖3 Fig.3

圖3 (續(xù))Fig.3(continued)

圖4為曲柄轉(zhuǎn)動兩周過程中四個(gè)角點(diǎn)的運(yùn)動.從運(yùn)動曲線中可以觀察到滑塊的自由運(yùn)動以及接觸碰撞事件的發(fā)生.

圖4 滑塊角點(diǎn)的運(yùn)動Fig.4 Motion of the slider corners

圖5顯示恢復(fù)系數(shù)對運(yùn)動的影響，恢復(fù)系數(shù)分別取為0.1,0.4,0.7,0.9.

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

圖5 恢復(fù)系數(shù)對滑塊角點(diǎn)1運(yùn)動的影響Fig.5 In fl uence of the restitution coefficient on the motion of corner 1

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

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

4 結(jié)論

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

將單邊約束、庫侖摩擦定律和牛頓碰撞定律均轉(zhuǎn)化為互補(bǔ)關(guān)系來表述，將動力學(xué)方程及其所有約束互補(bǔ)關(guān)系離散化為速度--沖量形式，采用分塊矩陣模型簡化推導(dǎo)過程，給出了一種新型線性互補(bǔ)公式.僅考慮單邊約束的多體系統(tǒng)動力學(xué)問題是常微分方程與線性互補(bǔ)的混合動力學(xué)問題，同時(shí)含等式約束及不等式約束的多體系統(tǒng)動力學(xué)問題是微分代數(shù)方程與線性互補(bǔ)的混合動力學(xué)問題.數(shù)值計(jì)算采用時(shí)間步進(jìn)算法，將線性互補(bǔ)解法嵌入至Moreau的時(shí)間步進(jìn)算法之中.數(shù)值計(jì)算過程中，系統(tǒng)動力學(xué)模型將在自由運(yùn)動、持續(xù)接觸運(yùn)動以及碰撞三種運(yùn)動狀態(tài)下相互切換，持續(xù)接觸過程中接觸狀態(tài)也隨時(shí)間變化，體現(xiàn)了典型的變拓?fù)浣Y(jié)構(gòu)特征.

最后，通過對含間隙滑移鉸曲柄滑塊機(jī)構(gòu)的建模與數(shù)值分析，驗(yàn)證了本文方法的有效性.

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富立.非光滑多體系統(tǒng)動力學(xué)線性互補(bǔ)方法.北京：清華大學(xué)出版社,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王曉軍,王琪.含摩擦與碰撞平面多剛體系統(tǒng)動力學(xué)線性互補(bǔ)算法.力學(xué)學(xué)報(bào),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張韻，李俊峰.碎石堆小行星的散體動力學(xué)建模與仿真方法綜述.力學(xué)學(xué)報(bào),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洪嘉振.計(jì)算多體系統(tǒng)動力學(xué).北京:高等教育出版社,1999(Hong Jiazhen.Computational Dynamics of Multibody Systems.Beijing:Higher Education Press,1999(in Chinese))

32齊朝暉.多體系統(tǒng)動力學(xué).大連:科學(xué)出版社,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范新秀,王琪.車輛縱向非光滑多體動力學(xué)建模與數(shù)值算法研究.力學(xué)學(xué)報(bào),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 網(wǎng)絡(luò)版發(fā)表.

1)河北省自然科學(xué)基金資助項(xiàng)目(A2013209221).

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

富立,胡鴻奎,富騰.多體系統(tǒng)接觸碰撞問題的牛頓--歐拉線性互補(bǔ)方法.力學(xué)學(xué)報(bào),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