基于交互模式的柔性體接觸碰撞動力學建模方法?

王檢耀 劉鑄永 洪嘉振

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

1 引 言

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

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

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

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

2 柔性體的動力學方程

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

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

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

圖1 柔性體的運動學描述Fig.1.Kinematics of a f l exible body.

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

根據速度變分原理,柔性體Bi的動力學方程為

式中,mP是P點質量,fP是P點受到的外力在慣性基上的坐標陣,C和K分別為有限元阻尼陣和剛度陣.

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

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

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

3 接觸碰撞模型

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

接觸點間的法向距離可以寫為

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

3.1 附加約束法

附加約束法中,接觸力虛功項可寫為[20]

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

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

式中M=diag,f=

圖2 點-面接觸對Fig.2.A node-to-segment contact pair.

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

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

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

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

3.2 罰函數法

罰函數法中,接觸力的虛功項可寫為[20]

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

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

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

3.3 交互模式法

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

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

接觸力項的虛功寫為

式中fc為接觸力坐標陣.

那么主體動力學模塊的方程為

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

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

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

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

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

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

4 數值仿真與實驗結果比較

4.1 實驗描述

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

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

圖4 實驗原理示意圖Fig.4.Schematic diagram of the impact experiment.

表1 實驗幾何與材料參數Table 1.Geometrical and material parameters.

4.2 結果比較

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

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

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

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

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

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

圖6 (網刊彩色)不同罰因子下的罰函數方法仿真結果 (a)板背面P1點的速度響應;(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 (網刊彩色)三種接觸模型的數值仿真與實驗結果 (a)板背面P1點的速度響應;(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.

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

5 多點碰撞算例

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

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

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

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

6 結 論

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

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

然后利用一柔性桿撞擊方板的實驗,通過激光測振儀測試撞擊點的響應,驗證所提方法的準確性.對該桿-板碰撞問題,分別使用罰函數法、附加約束法和交互模式法進行仿真.結果顯示,罰函數法非常依賴罰因子的選取,當罰因子選取合適時才能獲得可靠結果.三種方法的仿真結果都能與實驗結果符合較好.通過比較數值仿真結果與實驗結果,發(fā)現(xià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力學學報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振動工程學報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物理學報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)動力學(北京:高等教育出版社)第53頁]

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