富 立 胡鴻奎 富 騰

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

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

富 立2)胡鴻奎 富 騰

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

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

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

引言

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

近幾年來(lái)，非光滑多體系統(tǒng)動(dòng)力學(xué)有以下研究熱點(diǎn)及進(jìn)展.在基礎(chǔ)理論方面,解的存在與唯一性、碰撞本構(gòu)關(guān)系、動(dòng)力學(xué)建模等問(wèn)題的研究取得進(jìn)展.文獻(xiàn)[12]論證了存在或不存在庫(kù)侖摩擦情況下單、雙邊約束多體系統(tǒng)解的存在與唯一性問(wèn)題.文獻(xiàn)[13]對(duì)經(jīng)典的Painlev′e問(wèn)題給出新的分析方法和結(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)定化方法推廣至含二維摩擦的多體動(dòng)力學(xué)問(wèn)題.文獻(xiàn)[19-20]給出對(duì)非光滑事件進(jìn)行有效檢測(cè)的數(shù)值算法.在應(yīng)用研究方面,已有的非光滑多體系統(tǒng)動(dòng)力學(xué)理論成果被應(yīng)用于航空航天、機(jī)械、機(jī)器人、車(chē)輛、生物力學(xué)、天文學(xué)等諸多領(lǐng)域.如小行星聚集的接觸動(dòng)力學(xué)[2122]，航天器空間對(duì)接動(dòng)力學(xué)[23]，間隙機(jī)構(gòu)動(dòng)力學(xué)[2425]，顆粒物質(zhì)動(dòng)力學(xué)[26]，危巖崩塌軌跡的數(shù)值模擬[27]、機(jī)器人動(dòng)力學(xué)[28]、人體步態(tài)分析[29[31]、車(chē)輛非光滑動(dòng)力學(xué)建模與數(shù)值模擬[32]等.

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

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

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

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

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

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

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

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

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

1.2 持續(xù)接觸階段

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

其中，gN是接觸處的法向距離列陣、γT是接觸處的切向相對(duì)速度列陣，γ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)動(dòng)過(guò)程中，系統(tǒng)的接觸狀態(tài)隨時(shí)間發(fā)生變化，因此方程中的系數(shù)矩陣WN和WT以及約束反力列陣λN和λT的維數(shù)都是隨時(shí)間變化的.

1.3 碰撞階段

設(shè)碰撞始于t?終于t+，剛體動(dòng)力學(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ù)對(duì)角陣.

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

在多點(diǎn)接觸碰撞問(wèn)題中，當(dāng)一點(diǎn)接觸狀態(tài)改變時(shí)，如觸離狀態(tài)的改變、滯滑狀態(tài)的改變或發(fā)生碰撞等，都會(huì)影響到其他各點(diǎn)的接觸狀態(tài).在動(dòng)力學(xué)計(jì)算中需要對(duì)系統(tǒng)的接觸狀態(tài)進(jìn)行判斷.當(dāng)接觸點(diǎn)增多時(shí)，這個(gè)判斷過(guò)程是非常繁復(fù)的.現(xiàn)代非光滑力學(xué)對(duì)剛體系統(tǒng)接觸狀態(tài)的判斷建立了一種嚴(yán)格而有效的方法，即線性互補(bǔ)方法.該方法將系統(tǒng)動(dòng)力學(xué)方程及其單邊約束條件轉(zhuǎn)化為一個(gè)線性互補(bǔ)問(wèn)題，通過(guò)求解線性互補(bǔ)問(wèn)題來(lái)判斷系統(tǒng)的接觸狀態(tài)，具體參見(jiàn)文獻(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)方法將動(dòng)力學(xué)方程和單邊約束條件離散化，離散化的差分方程不以加速度和力為未知量，而是以速度和沖量取而代之.時(shí)間步進(jìn)算法的突出特點(diǎn)是數(shù)值積分過(guò)程中無(wú)需對(duì)非光滑事件進(jìn)行檢測(cè)，是真正的非光滑算法.該算法詳細(xì)內(nèi)容及主要步驟參見(jiàn)文獻(xiàn)[1,11].

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

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

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

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

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

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

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

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

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

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

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

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

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

將式(32)展開(kāi)

由式(33)

代入式(34)，得

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

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

經(jīng)整理簡(jiǎn)化，式(39)簡(jiǎn)寫(xiě)為

其中

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

令公式中的恢復(fù)系數(shù)eN=eT=0，即變?yōu)槌掷m(xù)接觸問(wèn)題的線性互補(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ù)：

初始條件：

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

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

其中

法向和切向約束矩陣為

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

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

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

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

圖3 Fig.3

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

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

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

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

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

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

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

圖6 間隙尺寸對(duì)滑塊角點(diǎn)1運(yùn)動(dòng)的影響(續(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)動(dòng)力學(xué)方程，由于含有等式約束，因此其數(shù)學(xué)模型是一組微分代數(shù)方程.在基本系統(tǒng)上，考慮接觸點(diǎn)處的法向和切向互補(bǔ)關(guān)系，并通過(guò)約束矩陣，將對(duì)應(yīng)的法向、切向拉氏乘子(約束反力)添加加至基本系統(tǒng)微分代數(shù)方程中，便構(gòu)成接觸狀態(tài)下的系統(tǒng)動(dòng)力學(xué)模型.同樣地，在基本系統(tǒng)上，考慮碰撞點(diǎn)處的法向和切向互補(bǔ)關(guān)系，并通過(guò)約束矩陣將對(duì)應(yīng)的法向、切向拉氏乘子(約束反力)加至基本系統(tǒng)微分代數(shù)方程中，便構(gòu)成碰撞時(shí)的系統(tǒng)動(dòng)力學(xué)模型.

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

最后，通過(guò)對(duì)含間隙滑移鉸曲柄滑塊機(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)動(dò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)動(dò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張韻，李俊峰.碎石堆小行星的散體動(dòng)力學(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)動(dòng)力學(xué).北京:高等教育出版社,1999(Hong Jiazhen.Computational Dynamics of Multibody Systems.Beijing:Higher Education Press,1999(in Chinese))

32齊朝暉.多體系統(tǒng)動(dò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范新秀,王琪.車(chē)輛縱向非光滑多體動(dòng)力學(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)動(dòng)力學(xué)及非線性動(dòng)力學(xué).E-mail：13231554976@163.com

富立,胡鴻奎,富騰.多體系統(tǒng)接觸碰撞問(wèn)題的牛頓--歐拉線性互補(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