王檢耀 劉鑄永 洪嘉振

(上海交通大學(xué)工程力學(xué)系,上海 200240)

1 引 言

工程中存在大量含接觸碰撞的非連續(xù)動力學(xué)過程,如航天器的交會對接[1]、高速列車受電弓接觸網(wǎng)間的碰撞振動[2]、機(jī)械系統(tǒng)中鉸鏈間隙引起的反復(fù)撞擊[3]、輪軌之間的碰撞振動[4]、有偏心轉(zhuǎn)子的打夯機(jī)與地面的碰撞[5]等.在這些非連續(xù)過程中,系統(tǒng)的拓?fù)浒l(fā)生突變,撞擊力峰值大、作用時間短,會給系統(tǒng)造成劇烈擾動,引起系統(tǒng)動力學(xué)性態(tài)的突變.因此采用高效高精度的方法對柔性體接觸碰撞問題進(jìn)行準(zhǔn)確的動力學(xué)仿真是尤為重要的.從物理上來看,接觸雙方的界面在接觸過程中需滿足互不嵌入的約束條件.對任意形狀柔性體的接觸碰撞問題,一般采用有限元方法離散,根據(jù)處理約束條件方式的不同,目前柔性體接觸碰撞問題的建模方法可以歸納為兩類:附加約束方法[6?9]和罰函數(shù)方法[10?13].

附加約束方法(ACM)通過引入接觸約束方程,與帶有Lagrange乘子的動力學(xué)方程聯(lián)立形成封閉的方程,可求得運(yùn)動學(xué)變量和接觸力.該方法優(yōu)點(diǎn)是接觸界面互不嵌入的約束條件得到精確的滿足.不足之處在于位置約束方程與動力學(xué)方程構(gòu)成指數(shù)-3微分代數(shù)方程,使得數(shù)值求解變得復(fù)雜[14?16],并且從未接觸到接觸時,速度發(fā)生突變[15,17],數(shù)值方法上也遇到困難.另一個缺點(diǎn)是接觸約束是單邊約束,方程的維數(shù)增加且需隨接觸狀態(tài)的變化而改變.

罰函數(shù)方法(PFM)將接觸作用視為彈簧阻尼力元,當(dāng)檢測到接觸發(fā)生時,在接觸域施加與嵌入量和嵌入速率相關(guān)的接觸力.與附加約束法相比較,罰函數(shù)法的優(yōu)缺點(diǎn)正好相反.它的優(yōu)點(diǎn)是不需要求解接觸約束方程,且方程維數(shù)保持不變;不足之處在于接觸約束條件只能近似地被滿足.理論上,只有罰因子(碰撞剛度系數(shù))取的足夠大時,結(jié)果才可靠.但是過大的罰因子會引起系數(shù)矩陣的病態(tài),同時降低了數(shù)值積分時間步長的臨界值.因此使用罰函數(shù)法時,罰因子的選取要權(quán)衡精度和數(shù)值穩(wěn)定性,往往需要通過多次試算,才能獲得合適的罰因子.為了克服罰函數(shù)法依賴于罰因子的困難,增廣的Lagrange方法被提出并廣泛使用.增廣的Lagrange方法同樣把接觸作為力元,通過迭代更新力元大小,直至嵌入量小于嵌入容許量.增廣的Lagrange方法減小了對罰因子的依賴,但是Seifried等[18]指出,嵌入容許量的大小同樣需要人為設(shè)定,過小的嵌入容許量會使得迭代次數(shù)大大增加,甚至有可能造成不收斂,而過大的嵌入容許量會使得精度不夠.

針對現(xiàn)有解決接觸碰撞問題方法的不足,本文提出基于交互模式的建模方法(IMM):將整個模型分為局部靜力學(xué)模塊及主體動力學(xué)模塊,局部靜力學(xué)模塊用于求解接觸力,主體動力學(xué)模塊用于求解運(yùn)動學(xué)變量,兩個模塊在計算中進(jìn)行位移和力的交互.在每一個積分步內(nèi),通過主體動力學(xué)模塊首先計算出局部區(qū)域的位移邊界,局部靜力學(xué)模塊據(jù)此計算出接觸力大小,再反饋到主體動力學(xué)模塊進(jìn)行下一積分步的計算.局部靜力學(xué)模塊類似附加約束法,施加約束方程使得接觸局部滿足約束條件;局部靜力學(xué)模塊計算出的接觸力以力元的形式施加到主體動力學(xué)模塊,因此主體動力學(xué)模塊的數(shù)值積分類似于罰函數(shù)法.交互模式方法兼具附加約束法和罰函數(shù)法的優(yōu)點(diǎn),同時避免了二者的不足:在接觸力計算上,接觸局部滿足約束條件保證了接觸力的精度,且無需人為選取參數(shù);在數(shù)值方面,避免了求解復(fù)雜的微分代數(shù)方程.本文利用圓柱桿撞擊方板的實(shí)驗(yàn),分別使用附加約束法、罰函數(shù)法和交互模式法進(jìn)行數(shù)值仿真,并與實(shí)驗(yàn)結(jié)果進(jìn)行比較,驗(yàn)證了交互模式法的準(zhǔn)確性及高效性.此外,通過對滑塊在有間隙滑槽內(nèi)滑動的多點(diǎn)碰撞問題仿真,進(jìn)一步驗(yàn)證了所提方法在處理一般性碰撞問題中的有效性.

2 柔性體的動力學(xué)方程

在柔性多體動力學(xué)中,浮動坐標(biāo)系方法是最常用的建模方法.該方法將柔性體的運(yùn)動分解為剛體運(yùn)動和相對的彈性變形,其中剛體運(yùn)動由浮動坐標(biāo)系描述,相對變形由模態(tài)坐標(biāo)或有限元節(jié)點(diǎn)坐標(biāo)描述.接觸碰撞是典型的非線性高瞬態(tài)過程,為了精確描述接觸的作用,本文采用有限元方法對柔性體進(jìn)行離散.

如圖1所示,柔性體Bi用集中質(zhì)量有限元離散,er是慣性基,eb是柔性體浮動基.柔性體上任一點(diǎn)P的絕對位置可表示為[19]

式中,ri為浮動基基點(diǎn)的位置坐標(biāo)陣,A為eb關(guān)于er的方向余弦陣,分別為點(diǎn)P變形前后相對浮動基原點(diǎn)的位置坐標(biāo)陣,為點(diǎn)P相對變形坐標(biāo)陣,(?)′表示(?)在浮動基下的坐標(biāo)陣.

圖1 柔性體的運(yùn)動學(xué)描述Fig.1.Kinematics of a f l exible body.

(1)式對時間t求導(dǎo)一次和兩次,可得點(diǎn)P速度、加速度的表達(dá)式:

根據(jù)速度變分原理,柔性體Bi的動力學(xué)方程為

式中,mP是P點(diǎn)質(zhì)量,fP是P點(diǎn)受到的外力在慣性基上的坐標(biāo)陣,C和K分別為有限元阻尼陣和剛度陣.

將(2)和(3)式代入(5)式中得到式中mi為質(zhì)量陣,?i為與vi二次項(xiàng)相關(guān)的慣性力陣,分別為外力陣和內(nèi)力陣.

將上述矩陣均寫為分塊矩陣形式,如下:

式中Mi=TTmiT,fi=TT(miT˙q˙i+?i?fie+fiσ).

3 接觸碰撞模型

有限元方法通過將接觸域進(jìn)行細(xì)致的空間離散化,在離散后的接觸域上形成多個接觸對,而接觸力元或約束是施加在各個接觸對上的.碰撞過程中,接觸點(diǎn)對的形成與分離表現(xiàn)了接觸區(qū)域的時變過程.最廣泛使用的接觸域離散形式為點(diǎn)-面接觸對(node to segment).如圖2所示,把撞擊物體Bi稱為從物體(slave body),被撞擊物體Bj稱為主物體(master body),從物體上的一點(diǎn)S與主物體單元表面上的最近投影點(diǎn)M構(gòu)成接觸對.接觸對間距矢量有如下形式:

接觸點(diǎn)間的法向距離可以寫為

式中n為主物體表面上點(diǎn)M位置的法向矢量.當(dāng)0,表示沒有接觸發(fā)生;當(dāng)=0,表示碰撞剛剛開始;當(dāng)0,則表示接觸點(diǎn)對之間有互相穿透.假設(shè)共有m組接觸點(diǎn)對,所有接觸點(diǎn)對的法向穿透量可組集為

3.1 附加約束法

附加約束法中,接觸力虛功項(xiàng)可寫為[20]

式中,λ為代表接觸力的Lagrange乘子,Gq=?G/?q為G的雅可比矩陣.

結(jié)合(11)式,系統(tǒng)的總虛功為

式中M=diag,f=

圖2 點(diǎn)-面接觸對Fig.2.A node-to-segment contact pair.

由于(15)式對任意的δq,δλ都成立,系統(tǒng)的動力學(xué)方程可寫為

(16)式是指數(shù)-3形式的微分代數(shù)方程,無論采用何種數(shù)值積分方法,求解該類型方程都有以下困難:雅可比矩陣的條件數(shù)過大、變步長算法中的局部誤差過大及數(shù)值穩(wěn)定性等問題[14].一般在求解之前需要先對其進(jìn)行處理,如降低指數(shù)[14]、通過增廣法[15]或縮并法[16]將微分代數(shù)方程轉(zhuǎn)變?yōu)槌Ｎ⒎址匠?

附加約束法的另一個困難是,從未接觸到接觸時(碰撞發(fā)生瞬時),接觸局部區(qū)域的速度是不連續(xù)的.如果不考慮速度突變,會引起數(shù)值違約及接觸力異常振蕩[17].Dong等[15]利用連續(xù)介質(zhì)力學(xué)的間斷面理論建立了速度突變條件.

附加約束法中由于施加了接觸約束方程,約束條件是被嚴(yán)格滿足的,但是需要求解指數(shù)-3形式的微分代數(shù)方程,并且要考慮未碰到碰瞬時的速度突變處理,給數(shù)值求解帶來困難.另一方面,當(dāng)接觸狀態(tài)改變,即有新的點(diǎn)對進(jìn)入接觸或已有接觸點(diǎn)對發(fā)生分離時,整個系統(tǒng)的系數(shù)矩陣的維數(shù)是要隨著改變的,在數(shù)值計算上要進(jìn)行多次的大維數(shù)矩陣稀疏化處理,會耗費(fèi)大量的計算時間.

3.2 罰函數(shù)法

罰函數(shù)法中,接觸力的虛功項(xiàng)可寫為[20]

其中,εN為罰因子(碰撞剛度參數(shù)).

那么系統(tǒng)動力學(xué)方程為

罰函數(shù)法沒有增加方程的維數(shù),并且無需求解微分代數(shù)方程,數(shù)值求解方便,但是罰因子εN的選取大多依賴經(jīng)驗(yàn),往往需要多次試算才能獲得可靠結(jié)果.

3.3 交互模式法

根據(jù)圣維南原理,接觸作用引起的應(yīng)力集中等現(xiàn)象只局限在接觸局部區(qū)域,因此可以將接觸局部的單元取出計算接觸力的大小,而局部單元所占質(zhì)量很小,慣性力相比于碰撞力為小量,從而可忽略局部的慣性效應(yīng),局部單元的接觸可做準(zhǔn)靜態(tài)處理.根據(jù)這個思路,將整個模型分為局部靜力學(xué)模塊及主體動力學(xué)模塊,局部靜力學(xué)模塊用于求解接觸力,主體動力學(xué)模塊用于求解運(yùn)動學(xué)變量,兩個模塊在計算中進(jìn)行位移和力的交互.在每一個積分步內(nèi),通過主體動力學(xué)模塊首先計算出局部區(qū)域的位移邊界,局部靜力學(xué)模塊據(jù)此計算出接觸力大小,再反饋到主體動力學(xué)模塊進(jìn)行下一積分步的計算.交互模式法的計算流程如圖3所示.

圖3 (網(wǎng)刊彩色)交互模式法的計算流程Fig.3.(color online)Schematic diagram of calculation steps using interactive mode method.

接觸力項(xiàng)的虛功寫為

式中fc為接觸力坐標(biāo)陣.

那么主體動力學(xué)模塊的方程為

求解(20)式得到整體的廣義坐標(biāo)q,從而可以得到局部靜力學(xué)模塊區(qū)域的廣義坐標(biāo)p=p(q).將接觸視為約束,局部靜力學(xué)模塊的方程為

式中f(p)為局部靜力學(xué)區(qū)域的內(nèi)力陣.

局部靜力學(xué)模塊的線性迭代求解格式為

式中,KT=?H/?p為切線剛度陣,p?,λ?為迭代的初始值.

求解(22)式得到代表接觸力的Lagrange乘子λ,轉(zhuǎn)換為fc=fc(λ)再反饋到(20)式進(jìn)行下一步的求解,由此形成了兩個模塊間位移和力的交互過程.

交互模式法中,局部靜力學(xué)模塊類似于附加約束法,通過約束方程和Lagrange乘子求解接觸力;而主體動力學(xué)模塊類似于罰函數(shù)法,接觸力以力元形式加入方程,積分求解運(yùn)動學(xué)變量.相比于附加約束法,無需求解微分代數(shù)方程;相比于罰函數(shù)法,無需人為選取罰因子.

4 數(shù)值仿真與實(shí)驗(yàn)結(jié)果比較

4.1 實(shí)驗(yàn)描述

圖4為桿-板正碰撞實(shí)驗(yàn)的示意圖.一個圓柱形鋼桿從一定高度落下,到最低位置時以一定速度撞擊保持靜止的鋁質(zhì)方形板.兩個物體均用細(xì)線懸掛,桿在平衡位置時,恰好接觸到方板的中心,且此時桿的軸線垂直于板面.為了測量碰撞持續(xù)的時間,將兩個物體通過導(dǎo)線連接到直流電源上,碰撞過程中電路連通將會產(chǎn)生電流信號.本次實(shí)驗(yàn)的碰撞時間約為1.7 ms,在此時間內(nèi)可認(rèn)為碰撞的兩個物體是自由的水平運(yùn)動.

碰撞實(shí)驗(yàn)中,碰撞力是難以直接測量的,只能通過測試速度和位移等間接物理量來反映碰撞過程.為了測量物體的速度響應(yīng),使用了Polytec GmbH公司生產(chǎn)的PSV-300F型激光測振儀.本次實(shí)驗(yàn)中,兩臺激光測振儀分別測試了板背面的中心點(diǎn)及桿后端的中心點(diǎn)的速度響應(yīng).桿和板的幾何與材料參數(shù)列于表1中.

圖4 實(shí)驗(yàn)原理示意圖Fig.4.Schematic diagram of the impact experiment.

表1 實(shí)驗(yàn)幾何與材料參數(shù)Table 1.Geometrical and material parameters.

4.2 結(jié)果比較

本文求解接觸碰撞問題的方法均是基于有限元離散的,因此首先須考慮網(wǎng)格劃分的影響.接觸碰撞問題對網(wǎng)格尺寸的要求非常高,既要保證碰撞引起的高頻率的彈性波在柔性體中的傳播,又要準(zhǔn)確表現(xiàn)接觸局部范圍內(nèi)的高應(yīng)力分布及接觸域的時變歷程.碰撞問題網(wǎng)格的劃分基于以下的原則[18]:對于非接觸區(qū)域,有限元網(wǎng)格應(yīng)覆蓋撞擊產(chǎn)生的最高頻率的彈性波的傳播,網(wǎng)格尺寸需滿足l≤c/20fmax,c為彈性波速,fmax為考慮的最高頻率,一般取50—100 kHz.根據(jù)該原則本例中的柔性體非接觸區(qū)域網(wǎng)格約取為5 mm;對于接觸局部區(qū)域,需要更小的網(wǎng)格,可先通過解析方法如Hertz接觸公式來粗略估計最大的接觸半徑a,再確定接觸區(qū)域內(nèi)的網(wǎng)格尺寸,以保證接觸過程中有多個接觸點(diǎn)對發(fā)生接觸.通過計算本例中a=0.49 mm,接觸域內(nèi)的網(wǎng)格尺寸取為0.2 mm.

對于交互模式法,局部靜力學(xué)處理部分的區(qū)域需要單獨(dú)拿出來,該區(qū)域的尺寸同樣按最大接觸半徑確定,即把最大接觸半徑范圍內(nèi)的節(jié)點(diǎn)所在單元作為局部靜力學(xué)處理的區(qū)域,如圖5所示.

圖5 (網(wǎng)刊彩色)桿-板碰撞的局部靜力學(xué)區(qū)域Fig.5.(color online)Local static region of the rodplate impact.

分別使用三種接觸模型,對該桿-板碰撞的算例進(jìn)行仿真.使用罰函數(shù)法時,需要人為選取罰因子大小.為了研究罰因子的取值對碰撞結(jié)果的影響,分別對不同的罰因子進(jìn)行仿真.不同罰因子下的仿真結(jié)果如圖6所示,圖6(a)和圖6(b)分別是測量點(diǎn)的速度響應(yīng)和碰撞力曲線.當(dāng)罰因子εN=105N/m時,由于過小的接觸剛度導(dǎo)致碰撞時間明顯過長,并且碰撞力和速度響應(yīng)的峰值偏小,振動的高頻特性也沒有反映出來.隨著罰因子值的增加,速度和碰撞力均逐漸收斂.εN=5×106和εN=107N/m時,結(jié)果已經(jīng)幾乎一致,此時的罰因子是合適的;當(dāng)繼續(xù)增大罰因子時,有可能因?yàn)檫^大的接觸剛度發(fā)生數(shù)值不穩(wěn)定.可以看到罰函數(shù)法的仿真結(jié)果非常依賴于罰因子的選取,須通過多次試算直至結(jié)果收斂.

與罰函數(shù)法相比,附加約束法和交互模式法均無需選擇碰撞參數(shù),一次仿真即可得到結(jié)果.將附加約束法和交互模式法的結(jié)果與罰函數(shù)法的收斂解結(jié)合在一起,并與實(shí)驗(yàn)結(jié)果比較,如圖7所示.圖7(a)給出了三種仿真方法與實(shí)驗(yàn)測得速度響應(yīng)結(jié)果的比較.由于實(shí)驗(yàn)無法直接測量碰撞力,圖7(b)給出了三種仿真方法碰撞力結(jié)果的互相比較.可以看到,三種接觸模型的仿真結(jié)果均與實(shí)驗(yàn)結(jié)果符合.從局部放大曲線上看,交互模式法與附加約束法的結(jié)果幾乎一致,而罰函數(shù)法的結(jié)果與它們?nèi)杂形⑿〉钠x.這是由于附加約束法和交互模式法都是通過約束方程求解得到的接觸力,而罰函數(shù)法通過施加彈簧剛度只能近似滿足約束.

該算例的仿真程序在Matlab平臺上編寫,計算機(jī)配置參數(shù)為:Intel Core i7-4770,3.40 GHz,內(nèi)存16 GB.三種方法的仿真的CPU時間列于表2中,可以看出三種方法求解該實(shí)驗(yàn)算例的計算效率有著較大差異.這是由于罰函數(shù)法中方程維數(shù)不隨接觸狀態(tài)改變,數(shù)值求解簡單,因此計算效率最高;交互模式法中,主體動力學(xué)模塊的數(shù)值積分和罰函數(shù)法相同,但是每一步內(nèi)都增加了局部小維數(shù)的非線性代數(shù)方程組的迭代求解,因此計算效率低于罰函數(shù)法;附加約束法中,由于矩陣維數(shù)增加且隨接觸狀態(tài)改變,加上數(shù)值求解的復(fù)雜性,使得計算效率最低.雖然罰函數(shù)法的單次求解效率最高,但是考慮到其不可避免的試算過程,其計算代價仍然是很大的.因此,交互模式法相比于附加約束法和罰函數(shù)法在求解效率上都有優(yōu)勢.

圖6 (網(wǎng)刊彩色)不同罰因子下的罰函數(shù)方法仿真結(jié)果 (a)板背面P1點(diǎn)的速度響應(yīng);(b)碰撞力時間歷程Fig.6.(color online)Simulation results of penalty function method with different penalty parameters:(a)Velocity of the back point P1on the plate;(b)time history of contact force.

圖7 (網(wǎng)刊彩色)三種接觸模型的數(shù)值仿真與實(shí)驗(yàn)結(jié)果 (a)板背面P1點(diǎn)的速度響應(yīng);(b)碰撞力時間歷程Fig.7.(color online)Comparison of simulation and experiment results:(a)Velocity of the back point P1 on the plate;(b)time history of contact force.

表2 三種接觸模型的計算規(guī)模Table 2.Computation scale of the three contact methods.

總體來說,基于交互模式的建模方法綜合了罰函數(shù)法和附加約束法的特點(diǎn),相比于罰函數(shù)法,其優(yōu)點(diǎn)是避免了人為選取參數(shù)及多次試算的過程;相比于附加約束法,優(yōu)點(diǎn)是避免了指數(shù)-3形式的微分代數(shù)方程和速度突變的處理,在精度上幾乎一致,并且在效率上占優(yōu).

5 多點(diǎn)碰撞算例

為了驗(yàn)證所提方法對于求解一般性碰撞問題的有效性,對滑塊在有間隙的滑槽中滑動問題進(jìn)行仿真.如圖8所示,一圓頭形滑塊置于光滑的長形滑槽中,滑塊和滑槽的幾何尺寸標(biāo)注于圖中,厚度均為1 mm,有限元離散后共有3747個節(jié)點(diǎn).滑塊與滑槽均為PVC材質(zhì),彈性模量E=3 GPa,泊松比μ=0.3,密度ρ=1800 kg/m3.滑塊的浮動坐標(biāo)系建于尾部中點(diǎn),浮動坐標(biāo)系原點(diǎn)在慣性系下的初始位置為[0 mm,?1.7 mm]T,滑塊軸線與滑槽成10°角.滑槽固定約束于地面,滑塊以沿軸線初始速度v=20 m/s彈出,滑塊滑動過程中將與滑槽發(fā)生多次多點(diǎn)的碰撞.

圖8 (網(wǎng)刊彩色)多點(diǎn)碰撞問題示意圖Fig.8.(color online)Schematic diagram of multi-point impact problem.

圖9 (網(wǎng)刊彩色)多點(diǎn)碰撞問題的仿真結(jié)果 (a)碰撞力時間歷程;(b)滑塊P點(diǎn)的位置軌跡Fig.9.(color online)Results of multi-point impact problem:(a)Time history of contact force;(b)position of point P on the slider.

采用交互模式法對該過程進(jìn)行動力學(xué)仿真,當(dāng)檢測到接觸位置發(fā)生變化時,局部靜力學(xué)模塊處理的區(qū)域也隨之改變.將交互模式法的結(jié)果與附加約束法的結(jié)果進(jìn)行比較,以驗(yàn)證其對多點(diǎn)碰撞問題仿真的準(zhǔn)確性.圖9(a)為碰撞力的時間歷程,圖9(b)為滑塊端點(diǎn)P的位置軌跡,在整個過程中共發(fā)生了四次不同位置的接觸碰撞,交互模式法的結(jié)果與附加約束法的結(jié)果均能較好地符合,說明交互模式法在多點(diǎn)多次碰撞的一般性碰撞問題中有效.在同等條件的仿真效率比較上,交互模式法在接觸階段耗時552 s,附加約束法在接觸階段耗時699 s,交互模式法的效率仍然較高,接觸階段的計算時間減少了21%.

6 結(jié) 論

針對現(xiàn)有解決柔性體接觸碰撞問題方法的不足,本文提出基于交互模式的建模方法.該方法將整個模型分為局部靜力學(xué)模塊及主體動力學(xué)模塊,局部靜力學(xué)模塊用以求解接觸力,主體動力學(xué)模塊用以求解運(yùn)動學(xué)變量,兩個模塊在計算中進(jìn)行位移和力的交互.在每一個積分步內(nèi),通過主體動力學(xué)模塊首先計算出局部區(qū)域的位移邊界,局部靜力學(xué)模塊據(jù)此計算出接觸力大小,再反饋到主體動力學(xué)模塊進(jìn)行下一積分步的計算.

相比于罰函數(shù)法,交互模式法通過局部模塊迭代計算接觸力,避免人為選取碰撞力參數(shù)和多次試算的過程;相比于附加約束法,交互模式法將約束方程(代數(shù)方程)和動力學(xué)方程(微分方程)分開求解,避免求解指數(shù)-3形式微分代數(shù)方程的困難,同時也無需處理從未碰到碰瞬時速度的非連續(xù)問題.

然后利用一柔性桿撞擊方板的實(shí)驗(yàn),通過激光測振儀測試撞擊點(diǎn)的響應(yīng),驗(yàn)證所提方法的準(zhǔn)確性.對該桿-板碰撞問題,分別使用罰函數(shù)法、附加約束法和交互模式法進(jìn)行仿真.結(jié)果顯示,罰函數(shù)法非常依賴罰因子的選取,當(dāng)罰因子選取合適時才能獲得可靠結(jié)果.三種方法的仿真結(jié)果都能與實(shí)驗(yàn)結(jié)果符合較好.通過比較數(shù)值仿真結(jié)果與實(shí)驗(yàn)結(jié)果,發(fā)現(xiàn)本文所提方法是行之有效的,并且所需計算代價要小于另兩種方法.此外,通過對滑塊-滑槽多點(diǎn)碰撞問題的仿真,進(jìn)一步驗(yàn)證了所提方法在處理一般性碰撞問題中的有效性.

[1]Klisch T 1998Multibody Syst.Dyn.2 335

[2]Ambrósio J,Pombo J,Rauter F,Pereira M 2009Multibody Dynamics:Computational Methods and Applications(Dordrecht:Springer Netherlands)p231

[3]Flores P,Koshy C,Lankarani H,Ambrósio J,Claro J C P 2011Nonlinear Dynam.65 383

[4]Lundberg O E,Nordborg A,Arteaga I L 2016J.Sound Vib.366 429

[5]Wang X H,Wang Q 2015Chin.J.Theor.Appl.Mech.47 814(in Chinese)[王曉軍,王琪2015力學(xué)學(xué)報47 814]

[6]Tur M,Fuenmayor F J,Wriggers P 2009Comput.Method Appl.M.198 2860

[7]Duan Y C,Zhang D G,Hong J Z 2013Appl.Math.Mech.Engl.34 1393

[8]Chen P,Liu J Y,Hong J Z 2016Acta Mech.Sin.32 1

[9]Weyler R,Oliver J,Sain T,Cante J C 2012Comput.Method Appl.M.205 68

[10]Tian Q,Zhang Y,Chen L,Flore P 2009Comput.Struct.87 913

[11]Qian Z J,Zhang D G 2015J.Vib.Eng.28 879(in Chinese)[錢震杰,章定國2015振動工程學(xué)報28 879]

[12]Zhang J,Wang Q 2016Multibody Syst.Dyn.38 367

[13]Yang Y F,Feng H B,Chen H,Wu M J 2016Acta Phys.Sin.65 240502(in Chinese)[楊永鋒,馮海波,陳虎,仵敏娟2016物理學(xué)報65 240502]

[14]Brenan K E,Campbell S L,Petzold L R 1996Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations(Philadelphia:SIAM)p45

[15]Dong F X,Hong J Z,Zhu K,Yu Z Y 2010Acta Mech.Sin.26 635

[16]Dan N,Haug E J,German H C 2003Multibody Syst.Dyn.9 121

[17]Taylor R L,Papadopoulos P 2010Int.J.Numer.Meth.Eng.36 2123

[18]Seifried R,Hu B,Eberhard P 2003Multibody Syst.Dyn.9 265

[19]Hong J Z 1999Computational Dynamics of Multibody Systems(Beijing:Higher Education Press)p53(in Chinese)[洪嘉振1999計算多體系統(tǒng)動力學(xué)(北京:高等教育出版社)第53頁]

[20]Géradin M,Cardona A 2001Flexible Multibody Dynamics:a Finite Element Approach(Chichester:John Wiley)p144