多體系統(tǒng)碰撞動力學中接觸力模型的研究進展1)

王庚祥 馬道林 劉洋 劉才山,2)

*(西安建筑科技大學機電工程學院,西安 710055)

?(北京大學工學院,湍流與復(fù)雜系統(tǒng)國家重點實驗室,北京 100871)

**(上海交通大學船舶海洋與建筑工程學院,上海 200240)

??(埃克塞特大學工程、數(shù)學和物理科學學院,英國?？巳?EX44 QF)

引言

接觸碰撞行為在工程機械[1]、生物醫(yī)學儀器與航空航天等領(lǐng)域[2]的多體系統(tǒng)中無處不在[3-4].典型場景如飛機起落架輪胎與地面接觸、車輛碰撞測試、嫦娥五號著陸與空間飛行器之間的對接過程等.隨著工程機械系統(tǒng)不斷向輕質(zhì)、高速、重載與精密方向發(fā)展[5-6],碰撞行為在短時間內(nèi)引起的局部彈塑性變形與劇烈沖擊力,愈發(fā)容易導(dǎo)致多體系統(tǒng)的結(jié)構(gòu)損傷與工作狀態(tài)失穩(wěn);對高端工程裝備的安全使用造成了嚴重的威脅[7-8].因此,有必要研究碰撞行為引起多體系統(tǒng)沖擊與振動的物理機制,分析碰撞與系統(tǒng)整體動態(tài)性能之間的耦合關(guān)系[9],為工程機械系統(tǒng)的安全穩(wěn)定工作提供理論依據(jù)[10].然而,多體系統(tǒng)中碰撞行為的研究與接觸力模型的發(fā)展歷程密切相關(guān)[11-12],碰撞力作為衡量碰撞行為是否對多體系統(tǒng)造成結(jié)構(gòu)損傷的重要指標高度依賴于碰撞力模型的精確性[13-14].

至目前為止,對動態(tài)接觸碰撞行為計算的模型主要分為兩類(如圖1):(1)基于Hertz 接觸模型發(fā)展的考慮能量耗散的連續(xù)接觸力模型[11,15-18],該類模型又可分為阻尼因子中分母包含與非包含初始碰撞速度的接觸力模型;(2)通過靜態(tài)壓力試驗獲得的準靜態(tài)彈塑性接觸力模型[11-12,19-20],該類模型從初始不考慮彈塑性過渡階段的非連續(xù)接觸模型發(fā)展為連續(xù)的準靜態(tài)接觸力模型.本文就兩類模型與相關(guān)的關(guān)鍵碰撞參數(shù)的發(fā)展歷程進行了著重介紹.

圖1 接觸力模型的分類Fig.1 Classification of contact force models

1 等效連續(xù)接觸力模型

碰撞行為是典型的非線性與不連續(xù)的非光滑系統(tǒng)[21],其中由于碰撞引起系統(tǒng)速度的跳躍具有典型的非光滑特征[22].當前描述碰撞體接觸行為的概念分為[17]:非光滑動力學方法(nonsmooth dynamics formulation)[23]和正則化方法(regularized approach)[24].其中基于幾何約束的非光滑動力學方法認為接觸體在碰撞過程中不發(fā)生變形,因此該方法也稱之為剛體方法[25].相反地,正則化方法認為碰撞體在接觸區(qū)域允許發(fā)生變形且接觸力是變形量的函數(shù),所以該方法也稱之為罰方法或者連續(xù)分析方法[26].

非光滑動力學方法中最常用的接觸建模方法是線性補償技術(shù)[27](linear complementary approach),該方法的核心思想是接觸剛體之間單邊約束的顯式表達[28].單邊約束被用來處理接觸問題的主要特征是利用了剛體不可穿透性的概念[29],意味著接觸體在碰撞過程中其接觸點不能越過碰撞體的邊界[30];防止碰撞發(fā)生變形保證其接觸力為正值[31].然而,該方法在處理考慮摩擦行為的接觸問題時可能會出現(xiàn)多解和無解的可能性,或者可能出現(xiàn)能量不守恒的情況[17];另外,該方法不利于多體系統(tǒng)碰撞動力學代碼的通用化.相反,正則化方法很好地避免了剛體方法遇到的諸多問題,該方法認為接觸體的每個接觸區(qū)域都覆蓋著一些分散在其表面的彈簧阻尼元件[32],該元件的變形與剛度大小依賴于碰撞體的幾何與材料屬性以及相對滲透深度[21],因此該方法也稱之為柔順接觸力模型(compliant contact force model).該方法的主要缺點在于難以確定剛度系數(shù)與阻尼系數(shù)以及能量指數(shù)等接觸參數(shù),以及在多體系統(tǒng)動力學仿真過程中引入大量高頻信息影響其求解效率[33].但是該方法的數(shù)學表達形式簡單直接且是滲透深度的連續(xù)函數(shù),便于多體系統(tǒng)碰撞動力學的程序通用化,因此該方法被廣泛應(yīng)用于多體系統(tǒng)動力學軟件[22,34-35],包括ADAMS,RecurDyn,EDEM 等.

連續(xù)接觸力模型的原型為Hertz 接觸力模型[36-37],考慮到Hertz 模型不能描述碰撞過程中的能量耗散;阻礙了Hertz 接觸力模型在多體系統(tǒng)碰撞動力學中的進一步應(yīng)用[38-40].為此,Kelvin 與Voigt 人為地在Hertz 接觸力模型中引入了阻尼項描述碰撞過程中的能量耗散,其基本形式為F=(K為Hertz接觸剛度,δ為相對接觸變形量,D為阻尼系數(shù),為相對碰撞速度),其中為彈性力項,為阻尼力項.

為了更加精確地計算多體系統(tǒng)中的碰撞行為,眾多學者在推導(dǎo)阻尼系數(shù)的過程中采用了不同的近似與求解方法[21,34,41-45],相應(yīng)地獲得了不同性質(zhì)的阻尼系數(shù).阻尼[46]是振動過程中能量耗散的一種度量,也是理解系統(tǒng)動力行為的重要組成部分,對抗震結(jié)構(gòu)的設(shè)計至關(guān)重要[47].阻尼一般分為兩種:(1)與頻率相關(guān)的阻尼稱為黏性阻尼[48-50];(2)與頻率無關(guān)的阻尼稱為遲滯阻尼[45].

另外,考慮到Hertz 接觸剛度獨立于接觸變形量,為了建立Hertz 剛度與相對接觸變形量之間的關(guān)系,修改了接觸力模型中的Hertz 剛度而忽略了能量指數(shù)與剛度系數(shù)的關(guān)系[41,51],造成了一系列的數(shù)值仿真問題.下面的內(nèi)容集中討論了接觸力模型中人為阻尼的類型與推導(dǎo)過程,闡述了能量指數(shù)與剛度系數(shù)之間的關(guān)系;最后總結(jié)了等效連續(xù)接觸力模型在計算碰撞行為時所面臨的挑戰(zhàn).

1.1 等效遲滯阻尼

1881 年Hertz[52]在研究兩個完全彈性體的接觸應(yīng)力時提出了著名的Hertz 接觸理論,其接觸力F是變形量的非線性函數(shù)

其中,δ為接觸變形量;n為能量指數(shù);K為Hertz 接觸剛度,其表達式為

其中,Ei和Ej為碰撞體的楊氏模量;Ri和Rj為碰撞體的接觸半徑;vi和vj為碰撞體的泊松比.值得注意的是Hertz 接觸剛度不同于胡克定律中彈簧的剛度,該剛度獨立于接觸變形量且完全依賴于接觸體的幾何與材料屬性(如式(2))[51],且該參數(shù)的量綱為N/m1.5,而不是胡克定律中彈簧剛度系數(shù)的量綱N/m.

眾所周知,碰撞過程實際上是碰撞體之間能量相互轉(zhuǎn)化與耗散的過程,針對純彈性碰撞行為的Hertz 接觸理論并不能滿足實際的工程需求,其原因在于式(1)并不能描述碰撞規(guī)程中的能量損耗(碰撞過程中,接觸體的壓縮路徑與恢復(fù)路徑相同)[18];于是為了描述碰撞行為的動態(tài)接觸過程以及碰撞導(dǎo)致的能量耗散,人為地在原始Hertz 接觸力模型(1)的基礎(chǔ)上引入表示能量耗散的阻尼項,將Hertz 接觸力模型擴展為可表示能量耗散的動態(tài)連續(xù)接觸力模型[53].該方法不僅能描述完整的動態(tài)接觸過程與能量損耗,并且建立了接觸力、相對變形量和變形速度之間的耦合關(guān)系[54],以及能快速獲得碰撞行為隨時間的動態(tài)變化過程,避免恢復(fù)系數(shù)法不能計算接觸進程的缺點[55];另外,為了達到精確計算碰撞過程中能量耗散的目的,眾多國內(nèi)外學者利用不同的近似推導(dǎo)方法從阻尼系數(shù)入手開展了一系列的研究工作[17].

為了在Hertz 接觸模型中引入人為的阻尼因子描述碰撞過程中的能量耗散,1960年Kelvin 與Voigt[17]利用線性彈簧阻尼模型計算了接觸過程中的能量耗散

其中式(3) 右邊的第一項線彈性項(胡克定律),K1為常系數(shù)剛度(不同于Hertz 接觸剛度);D1為常系數(shù)阻尼;為相對碰撞速度.然而,該模型在接觸開始與結(jié)束階段由于接觸速度不為零[56],導(dǎo)致其接觸力模型中的阻尼分量一直存在并引起接觸力的不連續(xù),甚至在恢復(fù)階段結(jié)束時侵入深度為零而相對接觸速度為負導(dǎo)致其接觸力也為負,顯然這種情況違反了接觸力學的物理意義(碰撞體之間不可能存在拉力)[24,26];以上原因?qū)е翶elvin-Voigt 模型產(chǎn)生的遲滯環(huán)不是一個完整的封閉區(qū)域,喪失了描述完整碰撞過程中能量耗散的能力.雖然該模型有以上缺陷,但在某些碰撞場景也有成功的應(yīng)用[56-58],例如Dubowsky等[57]利用該模型計算了三維間隙關(guān)節(jié)元素之間的接觸力,并提供了相關(guān)實驗數(shù)據(jù)驗證了該模型在一維振動沖擊系統(tǒng)中的成功使用[58].類似地,Rogers等[59]同樣將該模型應(yīng)用在考慮間隙關(guān)節(jié)的多體系統(tǒng)動力學仿真中,并認為該模型在材料阻尼較小時能近似地表現(xiàn)出黏彈性的屬性.Taylor等[60]利用該模型計算車輛輪胎與地面的法向接觸力.以上學者均認為接觸力模型的改進是精確計算多體系統(tǒng)碰撞動力學響應(yīng)的基礎(chǔ).

為了克服Kelvin-Voigt 模型在壓縮起始與恢復(fù)結(jié)束時刻出現(xiàn)的大于零與小于零的碰撞力,導(dǎo)致阻尼環(huán)呈現(xiàn)出不連續(xù)的狀態(tài),文獻[61,17]基于Hertz 理論,將阻尼系數(shù)描述成遲滯阻尼因子與接觸變形量之間的函數(shù)(D=χδn,χ為遲滯阻尼因子),避免了Kelvin-Voigt 模型在壓縮開始階段接觸阻尼力的存在,同時也消除了恢復(fù)階段結(jié)束時相對接觸速度不為零引起的缺點[61],保證了接觸力模型的連續(xù)性與遲滯環(huán)封閉的特點,其表達式為

其中m為接觸參數(shù).自此開始,在隨后的近40多年時間里,國內(nèi)外眾多學者在追求更加精確的遲滯阻尼因子的這條路上正式拉開了序幕[11,16,22,51,60-62];實際上,連續(xù)接觸力模型的發(fā)展歷程本質(zhì)上是阻尼因子的發(fā)展歷史[34],因為,式(4)右邊的第一項彈性力項(原始的Hertz 接觸模型)保持不變,第二項阻尼項的發(fā)展目的是為了準確描述碰撞過程中的能量耗散.

考慮到Hunt-Crossley 模型的特殊性[61],這里有必要詳細介紹該模型的推導(dǎo)過程,兩個小球之間發(fā)生碰撞過程如圖2.

圖2 碰撞體接觸過程Fig.2 Contact process between two contact bodies

首先,根據(jù)碰撞過程能量守恒原則,碰撞過程中耗散的能量可表示為

其中,T(-)為碰撞前系統(tǒng)動能;T(+)為碰撞后系統(tǒng)動能;m1與m2為碰撞體的質(zhì)量;v1i與v2i為碰撞前接觸體的初始速度;v1j與v2j為碰撞后接觸體的速度.

其次,根據(jù)碰撞前后動量守恒原則有如下關(guān)系式

該模型在推導(dǎo)阻尼因子過程中最大的特點在于:假設(shè)衡量能量損耗的恢復(fù)系數(shù)為初始相對碰撞速度的線性函數(shù)[15]

其中 α為常數(shù),一般取0.08 s/m 至0.32 s/m.

因此,將式(9)代入式(8),并考慮到 α 小于1,在忽略 α 的高階項后式(8)重新寫為

以上式(5)～式(10)是根據(jù)碰撞前后能量守恒的原則獲得了碰撞過程的能量耗散,該過程與碰撞過程中是否發(fā)生彈塑性變形無關(guān),即對彈性與彈塑性碰撞行為均適用.

另外,對式(4)中的阻尼力積分被視為碰撞過程中耗散的能量

為了獲得阻尼項耗散的能量,式(11)中碰撞中間過程對應(yīng)的相對變形速度需要表示成相對變形量的函數(shù)以確保式(11)能被積分.由此,該模型認為在壓縮結(jié)束階段碰撞體的初始動能被分為兩部分:一部分被轉(zhuǎn)化為彈性勢能,一部分為碰撞后系統(tǒng)的動能;其表達式為

其中v12為壓縮結(jié)束時碰撞體的共同接觸速度.此時,同樣根據(jù)動量守恒原則,其共同接觸速度可表示為

將式(13) 代入式(12),壓縮結(jié)束階段的應(yīng)變能可寫為

另外,壓縮結(jié)束階段的應(yīng)變能也可通過彈性力做功表示

聯(lián)立式(14)與式(15),可以獲得相對初始速度與最大相對接觸變形量之間的關(guān)系

根據(jù)該式,式(10)改寫成

然而,對于中間接觸過程(相對變形量0<δ<δmax),根據(jù)能量守恒,其接觸過程被認為是初始碰撞動能不斷向勢能轉(zhuǎn)化的過程

將式(14)與式(15)代入式(18),中間接觸過程對應(yīng)的相對變形速度可表示為

將式(19)代入式(11),由此可以通過對阻尼項積分獲得碰撞過程中的能量損耗

聯(lián)立通過碰撞過程能量守恒獲得的式(17)與通過對阻尼項積分獲得的式(20),再利用恢復(fù)系數(shù)是初始相對碰撞速度的線性函數(shù)的假設(shè)式(9),Hunt-Crossley 模型對應(yīng)的遲滯阻尼因子可表示為

遲滯阻尼因子是Hertz 接觸剛度、恢復(fù)系數(shù)與初始相對碰撞速度的函數(shù),該遲滯阻尼因子推導(dǎo)的關(guān)鍵步驟有4 點:(1) 假設(shè)恢復(fù)系數(shù)是相對初始碰撞速度的線性函數(shù)[63];(2) 在假設(shè)條件(1)的情況下,根據(jù)能量守恒原則將耗散能量表示為相對初始碰撞速度的函數(shù)[64];(3) 假設(shè)在壓縮結(jié)束時刻,系統(tǒng)初始動能被轉(zhuǎn)化為碰撞后系統(tǒng)的動能與勢能[34];(4) 將中間碰撞過程對應(yīng)的相對變形速度表示為相對變形量的函數(shù),以便對阻尼項積分時獲得碰撞過程的能量耗散[11].通過以上兩種不同思路推導(dǎo)的碰撞過程中的能量耗散完全等價,以此獲得連續(xù)接觸力模型中的遲滯阻尼因子.

以Hunt-Crossley 模型[34]的推導(dǎo)過程為基礎(chǔ),其他8 種遲滯阻尼因子——包括Lankarani-Nikravesh模型[65-66]、Zhiying-Qishao 模型[67]、Flores 等模型[68](首次提出該模型中的遲滯阻尼因子的是文獻[69])、Hu-Guo 模型[70]、Shen 等模型[71]、Safaeifar-Farshidianfar 模型[21]、Zhang 等模型[72-73]與Zhao 等模型[22]——在推導(dǎo)過程中采用了不同的假設(shè)或近似計算的方法獲得的阻尼表達式不盡相同;例如Lankarani-Nikravesh 模型[65-66]認為恢復(fù)系數(shù)為常數(shù),并未假設(shè)其為相對初始碰撞速度的函數(shù)以此獲得另外的遲滯阻尼因子.Flores 等模型[68-69]則認為在壓縮結(jié)束階段系統(tǒng)的初始動能分為三個部分:一部分為碰撞后的系統(tǒng)動能,一部分為損耗的能量,另一部分能量被轉(zhuǎn)化為勢能;以此獲得新的遲滯阻尼因子.

除此之外,其他接觸力模型在推導(dǎo)過程中將碰撞行為等效為非受迫的彈簧阻尼系統(tǒng),通過近似求解非線性振動方程的方式獲得阻尼因子[26,74-75]

其中,D為阻尼系數(shù);為相對變形加速度.由于該方程沒有解析解,一般通過近似求解的方式進行求解,并在此過程中假設(shè)恢復(fù)系數(shù)為初始碰撞速度的函數(shù),或者直接將遲滯阻尼因子假設(shè)為的形式(其中 β為恢復(fù)系數(shù)的函數(shù)).根據(jù)此方法推導(dǎo)遲滯阻尼因子的模型包括Herbert-McWhannell 模型[76]、Gonthier 等模型[75]、Carvalho-Martins 模型[15]、Gharib-Hurmuzlu 模型[77]與Hu 等模型[78].此外,還有一種頗具爭議的推導(dǎo)方法,即Lee-Wang 模型[79]將碰撞過程視為線性的彈簧阻尼系統(tǒng),通過線性振動方程的解析解來獲得阻尼因子;然而,在其接觸力模型中仍然保留了非線性的Hertz 接觸剛度,前后存在明顯矛盾[68,80].感興趣的讀者可以查閱相關(guān)文獻[79],表1 中給出了15 種連續(xù)接觸力模型的表達式[45],圖3 描述了表1 中部分連續(xù)接觸力模型中接觸力與變形量的變化趨勢.

圖3 常見連續(xù)接觸力模型的力與變形量之間的關(guān)系[51]Fig.3 Force-deformation relationship of the common contact force models[51]

表1 中的連續(xù)接觸力模型有以下特點:(1)能量指數(shù)n與接觸參數(shù)m均等于1.5[44];(2)所有遲滯阻尼因子均是接觸剛度、恢復(fù)系數(shù)與相對初始碰撞速度的函數(shù)[43];(3)所有遲滯阻尼因子的分母中均含有初始相對碰撞速度[34,81].該類連續(xù)接觸力模型在計算具有初始速度接觸體的碰撞行為時,由于其數(shù)學結(jié)構(gòu)簡潔[34],并在計算碰撞問題時只需要判斷接觸是否發(fā)生,不需要判斷接觸行為處于壓縮階段還是恢復(fù)階段[82-83],其計算流程相對簡單;因此該類模型被廣泛應(yīng)用于多體系統(tǒng)動力學軟件中(ADAMS 與RecurDyn 等)[84].然而,該類模型也存在不可調(diào)節(jié)的缺陷,由于該類模型中阻尼項分母中含有初始碰撞速度,限制了該模型對顆粒物質(zhì)中孤立波傳播特性刻畫的能力[85-88].因為顆粒物質(zhì)的初始狀態(tài)一般處于靜止狀態(tài),以至于顆粒物質(zhì)間的相對碰撞速度等于零[89];該接觸狀態(tài)導(dǎo)致表1 中連續(xù)接觸力模型的阻尼項無窮大,這顯然違反了顆粒物質(zhì)間的物理屬性[8];同時也給顆粒物質(zhì)的數(shù)值仿真引入了數(shù)值奇異問題[8,90-92].從表1 中連續(xù)接觸力模型的遲滯阻尼因子推導(dǎo)過程發(fā)現(xiàn),導(dǎo)致阻尼項分母中含有初始相對碰撞速度的原因有以下3 點:(1)利用能量守恒原則推導(dǎo)碰撞過程的能量耗散[93];(2)假設(shè)恢復(fù)系數(shù)是初始相對碰撞速度的函數(shù)[61];(3)直接將阻尼系數(shù)假設(shè)成分母中含有初始碰撞速度的形式[78].以上三類因素均直接或者間接與初始相對碰撞速度相關(guān),以至于推導(dǎo)過程中在遲滯阻尼因子分母中引入了初始碰撞速度.

1.2 等效黏性阻尼

為了避免阻尼項分母中含有相對初始碰撞速度帶來的數(shù)值奇異問題,在推導(dǎo)阻尼因子時繞開上述相關(guān)假設(shè)[94]與避免借助碰撞前后的能量守恒原則;將碰撞行為等效為單自由度的彈簧-阻尼模型(如圖4),并直接對單自由度欠阻尼非受迫振動方程求解,就能獲得分母中不含初始碰撞速度的黏性阻尼項[29,32,68].

圖4 彈簧-阻尼模型Fig.4 The spring-damper model

該類模型分為線性和非線性兩類系統(tǒng):(1)假設(shè)系統(tǒng)的剛度系數(shù)為常數(shù),其阻尼系數(shù)通過求解單自由度線性振動方程獲得,該類具有代表性的接觸力模型為Maxwell 模型[8];該模型與Kelvin-Voigt 模型類似,但是Maxwell 模型的阻尼系數(shù)是恢復(fù)系數(shù)的函數(shù);(2)考慮Hertz 接觸模型的非線性特征,基于非線性彈性碰撞行為特征建立的單自由度非線性振動系統(tǒng)的動力學模型(如式(22)),通過近似求解該動力學方程獲得接觸力模型中的阻尼因子;或者直接利用經(jīng)驗值確定其阻尼系數(shù)[74-75,95].為了清楚地詮釋該類模型與表1 中模型的區(qū)別,有必要分別介紹線性與非線性系統(tǒng)中阻尼系數(shù)的推導(dǎo)過程,為該文后續(xù)的接觸力模型分析奠定基礎(chǔ).

(1) 線性碰撞行為對應(yīng)的振動方程

需要注意的是:參數(shù)KL是線性彈簧的剛度,而不是Hertz 接觸剛度,所以方程(23)中彈性力項KLδ的變形量指數(shù)為1,其具體原因?qū)⒃?.1 小節(jié)解釋接觸力模型中剛度系數(shù)與指數(shù)系數(shù)關(guān)系的重要性.

式(23)作為常見的欠阻尼非受迫振動系統(tǒng)的動力學方程,其解析解為

其中系統(tǒng)的幅值和相角可由振動系統(tǒng)的初始條件獲得

將式(25)代入式(24),線性系統(tǒng)的解析解為

整個碰撞過程的持續(xù)時間為

另外,根據(jù)式(26),可以獲得振動系統(tǒng)的變形速度

因此,當碰撞行為結(jié)束時變形量完全恢復(fù),將式(27)代入式(28),其相對變形速度為

最后,根據(jù)Newton 恢復(fù)系數(shù)的定義

將式(29)與式(27)代入式(30),線性振動系統(tǒng)的阻尼系數(shù)為

Maxwell 模型的數(shù)學表達式為

(2) 非線性碰撞行為對應(yīng)的振動方程

其中 α為恢復(fù)系數(shù)的函數(shù).

為了獲得相對變形速度,對上式進行無量綱化,即定義

無量綱化后,式(35)改寫為

通過對上式積分可以獲得相對碰撞速度,以此導(dǎo)出阻尼系數(shù)

式(39) 正是被廣泛應(yīng)用于顆粒物質(zhì)仿真軟件(EDEM)中的連續(xù)接觸力模型[96].顯然,該模型中阻尼項不涉及初始相對碰撞速度引起的數(shù)值奇異問題,有效避免了遲滯阻尼針對顆粒物質(zhì)碰撞行為仿真的缺陷[97].其他接觸力模型的阻尼系數(shù)推導(dǎo)方式均借鑒了上述線性或非線性動力學方程的求解策略,采用了不同的近似手段和假設(shè)發(fā)展了其他6 種類似的模型.表2 中給出了接觸力模型的阻尼項分母中不含初始碰撞速度的接觸力模型,此類模型主要應(yīng)用于顆粒物質(zhì)的數(shù)值仿真[86,89,98].

表2 黏性接觸力模型(阻尼項分母中不含初始相對碰撞速度)Table 2 Contact force models with viscous damping factors(the denominators of the damping force do not include the initial impact velocity)

1.3 等效接觸力模型中阻尼系數(shù)的討論

顯然,1.1 節(jié)中通過碰撞過程中能量與動量守恒原則獲得阻尼與頻率無關(guān),但可保證碰撞前后系統(tǒng)的能量守恒[34];1.2 節(jié)中通過求解振動方程獲得與頻率相關(guān)的黏性阻尼,但無法保證碰撞前后系統(tǒng)的能量守恒[11,101].兩類阻尼在碰撞模型中不僅推導(dǎo)方式與使用背景不同,而且其展現(xiàn)的接觸動力學行為也不盡相同[102-103].為了詳細說明兩類阻尼的不同之處,假設(shè)兩個相同的鋼球在擁有不同初始速度(分別為0與8 m/s)下發(fā)生碰撞行為,其中小球的半徑為0.02 m,楊氏模量為2.07×1011Pa,泊松比為0.3.圖5中給出了兩種不同接觸力模型(Hunt-Crossley模型中的阻尼為遲滯阻尼,EDEM 軟件中的接觸力模型[104]對應(yīng)其黏性阻尼)中阻尼力的對比情況;其中遲滯阻尼隨著接觸變形量的增大緩慢增加,當碰撞行為處于壓縮階段的結(jié)束階段時,其碰撞速度趨近于0以至于其阻尼力從最大值(8000N)急劇減小至0N;在恢復(fù)起始階段雖然相對接觸速度較小但其接觸變形量較大,以至于相對接觸速度發(fā)生較小變化其阻尼力將急劇增加;另外,在恢復(fù)的結(jié)束階段其接觸變形量為0以至于其阻尼力也為0.因此,包含遲滯阻尼的接觸力模型中描述能量損耗的阻尼環(huán)不會出現(xiàn)“拉力”區(qū)域(見圖5),其主要原因在于恢復(fù)結(jié)束階段其阻尼力快速趨近于0.相反,EDEM 軟件中的接觸力模型中黏性阻尼力雖然在整體趨勢上與遲滯阻尼力保持一致(見圖6),符合碰撞行為中阻尼項的表現(xiàn)形式;但是在細節(jié)上存在明顯的區(qū)別.在壓縮起始階段與恢復(fù)結(jié)束階段,黏性阻尼力戲劇化地增加,其原因在于黏性阻尼力中δ0.25在接觸變形量較小時其值較大.因此黏性阻尼力在壓縮起始階段與恢復(fù)結(jié)束階段明顯大于遲滯阻尼,這也是為什么黏性阻尼對應(yīng)的阻尼環(huán)會在恢復(fù)臨近結(jié)束階段出現(xiàn)“拉力”區(qū)域(見圖6).黏性阻尼力在經(jīng)過極速增加區(qū)域后,當接觸變形量逐漸增大時,阻尼力緩慢增加;即使在靠近壓縮結(jié)束階段時其阻尼力也是緩慢減小,與遲滯阻尼力急劇減小不同.另外,關(guān)于“拉力”區(qū)域另一種直觀的解釋為:該區(qū)域?qū)?yīng)的接觸變形量大于0,但其碰撞力為負值;原因在于恢復(fù)結(jié)束階段其阻尼力大于彈性力.由于該類模型主要用于顆粒物質(zhì)之間的接觸碰撞行為描述[105-107],由于該區(qū)域在整個阻尼環(huán)中只出現(xiàn)在恢復(fù)的結(jié)束階段對大量顆粒物質(zhì)之間的碰撞行為與顆粒之間孤立波的傳播基本不產(chǎn)生影響,這也正是該模型在EDEM 軟件[104]中被廣泛使用的原因.圖5 中兩類接觸力模型中力與接觸變形量之間的關(guān)系在整體趨勢上保持一致,但包含黏性阻尼的接觸力小于其包含遲滯阻尼的接觸力,其原因在于最大遲滯阻尼力大于最大黏性阻尼力(因為兩類接觸力模型中Hertz 彈性項相同).

圖5 兩類阻尼的比較Fig.5 Comparative analysis between two different damping systems

圖6 擁有不同阻尼類型的接觸力模型Fig.6 Contact force models with different damping systems

1.4 能量指數(shù)與接觸剛度之間的關(guān)系

關(guān)于能量指數(shù)與接觸剛度的關(guān)系往往被大家所忽略而造成不易察覺的錯誤.大家通常認為Hertz 模型中能量指數(shù)等于1時[51],其數(shù)學表達式就是彈性力學中的胡克定律;其實兩者之間并沒有直接的聯(lián)系,并且該觀點從量綱的角度出發(fā)就不正確.因為Hertz 剛度的單位為N/m1.5,當能量指數(shù)等于1 時,最后獲得接觸力的量綱不是N,明顯違反了力學常識.這也是為什么Hertz 接觸模型中能量指數(shù)等于1.5,說明接觸體為金屬材料且接觸應(yīng)力呈拋物線分布(這里需要注意的是:Hertz 接觸剛度本身是一個特殊的物理量,有別于傳統(tǒng)的剛度系數(shù)(變形量相關(guān)的系數(shù));而Hertz 接觸剛度完全取決于接觸體的幾何與材料屬性且獨立于接觸變形量).因此,能量指數(shù)的大小必須與接觸剛度的量綱保持一致,否則將造成不易察覺的錯誤,例如表2 中的Kuwabara-Kono 模型[99]與Lee-Wang 模型[79],就忽視了Hertz接觸剛度量綱與能量指數(shù)的關(guān)系;很明顯,Kuwabara-Kono 模型中阻尼力的量綱為N/s;Lee-Wang 模型中阻尼力的量綱為N/m0.25.以上兩種模型已經(jīng)在顆粒物質(zhì)的動態(tài)仿真中有所應(yīng)用,并取得了 “合理”的仿真結(jié)果,但依然不能避免量綱不正確的缺點.

那么在何種情況下,基于Hertz 理論擴展的連續(xù)接觸力模型中的能量指數(shù)會發(fā)生變化？一般有兩種情況能量指數(shù)會發(fā)生變化:(1) 接觸材料不再是金屬材料,而是玻璃或者高分子聚合物等材料,直接影響了接觸應(yīng)力的分布形式;(2)其接觸形式發(fā)生變化,由點接觸變化為線接觸或者面接觸.然而,至目前為止就表1 與表2 中的接觸力模型而言,這種情況還未發(fā)生[51];也就是說,只要是基于Hertz 接觸模型拓展的考慮能量耗散的接觸力模型,其能量指數(shù)必定等于1.5;除非改變接觸剛度的數(shù)學表達式,使得接觸剛度系數(shù)的單位不再是N/m1.5.

另外,為什么要一直強調(diào),接觸剛度與能量指數(shù)的關(guān)系很容易被忽視,并且導(dǎo)致不易察覺的錯誤.因為在數(shù)值仿真當中,以考慮間隙關(guān)節(jié)的多體系統(tǒng)動力學仿真為例,間隙關(guān)節(jié)元素之間的碰撞行為屬于典型的多重碰撞與多重壓縮問題,一般采用連續(xù)接觸力模型計算其碰撞力;其中間隙關(guān)節(jié)元素之間的相對接觸變形量一般在[1×10-10,1×10-5] m,其能量指數(shù)的大小決定了接觸力的大小處于不同的數(shù)量級;直接影響間隙關(guān)節(jié)元素之間的接觸碰撞判定條件.例如文獻[51,108-109]中接觸模型的剛度系數(shù)單位為N/m2,但假設(shè)其能量指數(shù)仍然等于1.5,此時仍能獲得“合理”的仿真結(jié)果,且該類錯誤很難發(fā)現(xiàn).如果為了配合接觸剛度的量綱,假設(shè)能量指數(shù)等于2;此時其接觸力將處于另一個量級甚至約等于零,這將原本處于接觸的碰撞體根據(jù)其仿真結(jié)果判定為未接觸,給考慮間隙關(guān)節(jié)的多體系統(tǒng)動力學仿真帶來了嚴重的數(shù)值仿真問題;其根本原因在于能量指數(shù)的大小直接改變了 δn的數(shù)量級,具體細節(jié)可以參考文獻[51].

造成以上錯誤的根本原因在于對連續(xù)接觸力模型中Hertz 接觸剛度的改進,使接觸力模型中的剛度與接觸變形量之間產(chǎn)生耦合關(guān)系,從整體上改善接觸力模型的精度.另外,為了改善表示能量耗散的阻尼因子的精度,甚至通過擬合能量指數(shù)的方式;然而,在此過程中忽略了接觸剛度與能量指數(shù)之間的協(xié)調(diào)關(guān)系.由此可見,為了改善接觸力模型的精度,最好的方式是不斷改進表示能量耗散的阻尼因子的精度;而不建議在忽略接觸剛度量綱與能量指數(shù)之間關(guān)系的情況下,對接觸剛度進行改進;更不可基于Hertz 接觸剛度的情況下對能量指數(shù)進行擬合.除非同時改變接觸剛度,并保證接觸剛度與能量指數(shù)之間的協(xié)調(diào)關(guān)系以確保接觸力量綱的正確性.

1.5 關(guān)于等效連續(xù)接觸力模型的討論

當前所有的連續(xù)接觸力模型(見表1 與表2)均是人為的在Hertz 接觸力模型[16-17](F=Kδn)中引入表示能量耗散的阻尼項(D=χδn)[110],將準靜態(tài)Hertz 接觸力模型擴展為可表示能量耗散的動態(tài)連續(xù)接觸力模型(F=Kδn+χδm)[16-17].然而,通過上述連續(xù)接觸力模型描述碰撞行為存在以下四個方面的問題:(1)雖然可通過對阻尼項積分將碰撞過程中的能量耗散均勻化分布在壓縮與恢復(fù)路徑上[65],但是無法通過人為添加的阻尼因子解釋能量耗散的物理原因[111];(2)在碰撞體幾何形狀與材料屬性確定的情況下,Hertz 接觸剛度可描述低速碰撞環(huán)境下的彈性碰撞(碰撞速度低于臨界彈性碰撞速度ρ為碰撞體的密度)[66];然而當碰撞速度大于臨界彈性碰撞速度時,碰撞體將發(fā)生彈塑性變形[112];此時的Hertz 接觸剛度高估了彈塑性碰撞階段的接觸剛度,無法正確描述彈塑性碰撞行為[113];以至于目前連續(xù)接觸力模型中的接觸剛度在彈塑性接觸過程中高估了壓縮過程中的應(yīng)變能與接觸力的大小[61];(3)在碰撞過程中,一部分能量在局部接觸碰撞過程中被損耗,另一部分能量在壓縮過程中以應(yīng)變能的形式儲存在碰撞體中,并以應(yīng)力波的形式影響未受沖擊區(qū)域的變形狀態(tài)[114],但是當前的連續(xù)接觸力模型均未考慮碰撞引起的應(yīng)力波在接觸體中的傳播效應(yīng)[115];(4)為了分析碰撞引起的應(yīng)力波傳播效應(yīng),需要借助連續(xù)介質(zhì)力學對碰撞問題進行求解,研究柔性碰撞體在吸收應(yīng)變能之后的變形狀態(tài);當碰撞時間接近或者大于接觸體的基礎(chǔ)振動周期時,碰撞激發(fā)的應(yīng)力波將從邊界返回并達到碰撞位置[116],與還沒有結(jié)束的碰撞過程產(chǎn)生耦合效應(yīng)激發(fā)多重壓縮-恢復(fù)等過程[117].然而,由于連續(xù)接觸力模型關(guān)于彈塑性碰撞過程中能量損耗與接觸力的計算存在偏差,影響了柔性體吸收能量的大小,以至于在接觸碰撞持續(xù)時間與相對碰撞速度的預(yù)測上存在誤差,使得多體系統(tǒng)中碰撞行為與柔性變形之間的耦合關(guān)系變得更加復(fù)雜.因此,正確表征碰撞過程中的能量耗散對于揭示能量耗散機理至關(guān)重要[118],也是揭示多體系統(tǒng)中碰撞與整體系統(tǒng)柔性特征之間耦合關(guān)系的前提條件.

更重要的是,由于壓縮過程中系統(tǒng)能量等于初始碰撞動能減去碰撞體吸收的應(yīng)變能

該等式清楚地描述了壓縮行為本質(zhì)上是初始動能不斷向應(yīng)變能轉(zhuǎn)化的過程.然而,在能量的轉(zhuǎn)化過程中,因為Hertz 接觸剛度明顯大于彈塑性變形階段的接觸剛度,根據(jù)式(40)可知,在彈塑性碰撞的壓縮過程中碰撞初始動能轉(zhuǎn)化為應(yīng)變能的效率變快,以至于能量快速被消耗而改變了碰撞接觸時間,同時也改變了接觸速度的變化規(guī)律.另外,由于阻尼因子是接觸剛度與恢復(fù)系數(shù)的函數(shù),所以對阻尼因子積分計算碰撞過程中能量耗散時,會發(fā)現(xiàn)壓縮過程中耗散能量[61]明顯比實際情況多(由于接觸剛度的原因),且耗散速度比實際的快,無法準確計算碰撞過程經(jīng)歷的時間.圖3 中可以看出耗散能量多(遲滯環(huán)的面積越大)的接觸力模型不僅相對侵入深度較小且相對碰撞速度變化急促.所以,當前動態(tài)連續(xù)接觸力模型描述彈塑性碰撞行為能量耗散不準確的根本原因在于壓縮過程中接觸剛度與實際情況不符.更進一步,正因為Hertz 接觸剛度高估了彈塑性碰撞階段的接觸剛度,改變碰撞體在壓縮過程中吸收能量的效率,進而影響碰撞過程中的能量損耗,使多體系統(tǒng)中彈塑性碰撞行為與柔性變形之間能量傳播與耗散機制更加復(fù)雜,無法準確獲得彈塑性碰撞引起系統(tǒng)振動機理的真實原因.

2 準靜態(tài)接觸力模型的研究現(xiàn)狀

在實際碰撞行為中,彈塑性碰撞行為相比純彈性碰撞現(xiàn)象更加常見,Lankarani等[66]指出當初始碰撞速度大于或者等于(說明較小的碰撞速度就能引起接觸體的局部彈塑性變形)時碰撞表面將會出現(xiàn)永久接觸變形.魏悅廣等[119]強調(diào)碰撞體在重載荷下必須考慮接觸體的彈塑性變形.在當前高端工程裝備與航空航天等領(lǐng)域,經(jīng)常面臨高速重載輕質(zhì)的工況環(huán)境,因此,相比純彈性碰撞現(xiàn)象,彈塑性接觸碰撞現(xiàn)象更值得關(guān)注.

最早彈塑性變形現(xiàn)象的研究主要來源于人們對接觸材料硬度與工程機械中結(jié)合面粗糙度的探索[120].其原因在于現(xiàn)實中根本不存在完全平整的接觸面,物體的接觸表面實際上是由眾多幾何尺寸不同的微凸體構(gòu)成(如圖7[121]為了便于研究一般近似假設(shè)微凸體的形狀為圓柱形、球體、正弦波或者波浪形[20,52];其中球體之間的碰撞行為在實際工程應(yīng)用中普遍存在),即從宏觀角度上看似光滑平整的零件表面,但實際上其表面接觸形貌呈現(xiàn)不同程度的粗糙和凹凸不平的特征[122].兩個機械表面的真實接觸情況是許多離散微凸體相互擠壓與滑動的過程,導(dǎo)致其粗糙表面的真實接觸面積遠小于名義接觸面積[44],造成了在極小的真實接觸面積上承受較大的載荷,以至于產(chǎn)生彈塑性變形最終導(dǎo)致其局部表面被壓潰[123],通常情況下塑性變形必定伴隨著彈性變形[124].

圖7 表面粗糙度形貌Fig.7 The surface roughness topography

一般而言,接觸行為按照是否“貼合”分為協(xié)調(diào)接觸(conforming contact)和非協(xié)調(diào)接觸(nonconforming contact)[125].其中協(xié)調(diào)接觸為兩個固體表面在無變形時精確地或者相當接近地 “貼合” 在一起;而非協(xié)調(diào)接觸為具有不相似外形且無變形的兩個固體接觸時,首先在一個點或沿一條線相碰,分別稱為點接觸和線接觸.以上對協(xié)調(diào)與非協(xié)調(diào)接觸的定義是從宏觀接觸表面出發(fā),忽視了表面粗糙度對接觸行為的影響[21,47].然而,當從微觀接觸表面出發(fā)時,任意接觸表面均不可能出現(xiàn)“貼合”的接觸行為,而是均為點接觸或者線接觸;因此,當考慮接觸表面微觀粗糙度時,其協(xié)調(diào)接觸與非協(xié)調(diào)接觸并沒有明顯的界限,均可視為非協(xié)調(diào)接觸[123,126-127].

當忽略接觸過程中應(yīng)變率的變化與塑性流,并認為接觸體各向同性,一般將整個接觸過程分為完全彈性階段、彈塑性階段和完全塑性等三個階段[128].其中在完全彈性接觸階段,Hertz 理論提供了一個封閉的力學本構(gòu)關(guān)系;一旦接觸行為達到材料的屈服極限,其接觸行為將進入彈塑性階段;作為完全彈性和完全塑性變形的過渡階段,使得整個接觸行為連續(xù)的彈塑性本構(gòu)方程則相對難以獲得[12];當達到彈塑性接觸變形的極限時,接觸行為將進入完全塑性階段,該階段的本構(gòu)關(guān)系只有在簡單的假設(shè)條件下才能獲得其解析解[12,129].

最早對完全塑性問題的研究源自于人們對材料硬度的探索,1900年Brinell 研究了塑性變形區(qū)域的應(yīng)力分布[52],并提出接觸體的硬度為接觸力除以接觸區(qū)域的整體面積;隨后,1908 年Meyer 認為接觸材料的硬度等于接觸力除以在壓痕方向上的接觸投影面積[130](值得注意的是,1933 年Abbott 與Firestone 在研究軸承表面特征的文章中并沒有提及接觸過程中塑性變形的工作[12]).1948 年Tabor[131]在Ishlinskii 的基礎(chǔ)上分析了一個硬性球形壓頭被壓入一個軟性金屬的表面,當外力卸載時,金屬表面發(fā)生了塑性的永久變形.1972 年,Lee等[132]利用有限元與實驗測量的方式獲得接觸材料硬度與壓痕深度沒有必然的聯(lián)系.1985 年Sindair等[133]在忽略摩擦的情況下研究了剛性球體與彈塑性半空間接觸體之間的壓痕,分析了接觸過程中的應(yīng)力與應(yīng)變之間的關(guān)系.1986 年Johnson[130]系統(tǒng)地研究了不同幾何形狀接觸體在完全彈性和彈塑性以及完全塑性接觸階段的本構(gòu)關(guān)系,獲得了一系列適用于不同接觸場合的有效彈塑性接觸力模型.

直到1987 年,Chang等[134](CEB 模型)基于統(tǒng)計模型根據(jù)G-W 模型的假設(shè)條件研究了粗糙接觸表面在大載荷下的彈塑性變形狀態(tài),首次提出了包含彈性和塑性變形階段的準靜態(tài)接觸模型.然而,CEB 模型[134]在描述完全彈性階段過渡到完全塑性階段出現(xiàn)了不連續(xù)現(xiàn)象,即塑性階段的本構(gòu)關(guān)系并不能使接觸過程從彈性階段平滑地過渡到塑性階段.為了彌補CEB 模型中的不連續(xù)特征,Zhao-Maletta-Chang(ZMC 模型)[83]利用多項式插值的方式描述彈塑性接觸階段的力與變形之間的關(guān)系;然而,該模型依然沒有很好地解決CEB 模型的不連續(xù)問題.Thornton[135]在忽略彈塑性過渡階段的情況下,基于Hertz 理論提出了可連續(xù)描述整個接觸過程(包括完全彈性和塑性變形)的接觸力模型.在此基礎(chǔ)上,Johnson[130]推導(dǎo)了剛性球壓入半空間的平均接觸壓力與下壓量的關(guān)系,獲得了簡化的彈塑性接觸力模型.Stronge等[136]基于Johnson 模型[137]通過修正壓痕與接觸半徑之間的幾何關(guān)系提出了包含彈性、彈塑性與完全塑性變形的接觸力與接觸變形量之間的函數(shù)關(guān)系式.為了進一步提高準靜態(tài)接觸力模型的精度,并考慮到有限元方法的不斷發(fā)展,Jackson等[138]基于有限元法根據(jù)Mises 屈服理論提出的JG 模型,證實了ZMC 模型并沒有很好解決CEB 模型的不連續(xù)問題.另外,Kogut 等(KE 模型)[139]采用有限元法驗證了CEB 模型在塑性變形階段低估了接觸力的大小,而ZMC 模型在彈塑性階段低估了接觸力的大小;并提出了一種無量綱的連續(xù)彈塑性本構(gòu)方程可描述整個彈塑性接觸過程.Shankar等[140]利用有限元法擴展了KE 模型與JG 模型在彈塑性接觸過中彈性階段過渡到塑性階段的臨界值.Zhang等[141]利用有限元法在同時根據(jù)Hertz 理論與Mises 屈服理論情況下,忽略碰撞過程中的摩擦效應(yīng),根據(jù)不同接觸階段變形量與載荷之間的關(guān)系基于Newton-Raphson 方法獲得彈塑性接觸階段的接觸力模型.Etsion等[142]借助有限元法在研究兩個小球之間彈塑性碰撞過程中卸載階段的力與變形之間的關(guān)系,以此為基礎(chǔ),提出了可連續(xù)描述彈塑性接觸行為的力學模型.Du等[143]分析了兩個小球在發(fā)生彈塑性法向碰撞時的能量耗散,提出了可描述低速環(huán)境下發(fā)生的彈塑性碰撞行為,并通過有限元法驗證了該模型中碰撞參數(shù)的合理性.Burgoyne等[144]在研究顆粒物質(zhì)中孤立波的傳播特征時,采用有限元分析方法導(dǎo)出了材料性能的經(jīng)驗函數(shù),提出可用于彈塑性碰撞行為仿真的準靜態(tài)彈塑性接觸力模型,并通過Hopkinson 實驗驗證了該模型的正確性.Komvopoulos等[145]采用擬合有限元結(jié)果的方式將碰撞過程分為彈性、小變形的彈塑性、中等變形的彈塑性與大變形的彈塑性等四個階段.在類似的工作中,Brake[146]做出了突出的貢獻,將彈塑性接觸行為分為彈性、彈塑性和塑性接觸階段,其中彈性階段服從Hertz 理論,塑性階段服從Johnson 理論或者利用Meyer 硬度指數(shù)描述完全塑性接觸階段;忽略接觸過程中硬度的變化,利用多項式插值的方式近似彈塑性接觸階段的本構(gòu)方程.馬道林等[113]借助Brake 模型的思想,利用合理的邊界條件改進了過渡階段的函數(shù)關(guān)系,基于Johnson彈塑性接觸理論發(fā)展了可連續(xù)從彈性變形過渡到塑性變形的準靜態(tài)彈塑性分段接觸力模型.

從準靜態(tài)彈塑性接觸力模型的發(fā)展歷程可看出該類模型的發(fā)展同樣經(jīng)歷了從不連續(xù)到連續(xù)的演化過程.建立該類模型的方法主要分為3 種:(1)以固體力學為基礎(chǔ),在簡化假設(shè)條件的情況下推導(dǎo)精度較低的彈塑性解析模型;(2)借助有限元結(jié)果進行多項式擬合獲得半解析解的彈塑性接觸力模型;(3)在半解析彈塑性接觸力模型的基礎(chǔ)上,通過合理的假設(shè)與邊界條件推導(dǎo)精確的解析解模型.表3 中給出了9 種常見的連續(xù)準靜態(tài)彈塑性接觸力模型,圖8中描述了表3 中準靜態(tài)接觸力模型中接觸力與變形量之間的關(guān)系.

表3 連續(xù)準靜態(tài)彈塑性接觸力模型Table 3 Continuous quasi-static elastoplastic contact force models

續(xù)表3

一般而言,彈塑性接觸行為在加載時分為彈性、彈塑性與完全塑性接觸階段,而在卸載時只有彈性變形(如圖8);整個彈塑性接觸行為的加載與卸載階段又稱為壓縮與恢復(fù)階段.圖8 中接觸力模型加載與卸載曲線的差異性與研究對象、接觸材料屬性、邊界條件以及曲線擬合與近似方式不一致有關(guān).為了清楚描述整個彈塑性接觸過程,圖9 給出了表3 中一般性的接觸力與彈塑性變形量之間的關(guān)系,其中彈性階段在整個彈塑性接觸行為中只會在初始接觸階段發(fā)生,并且只占整個接觸過程的極短一部分,幾乎可以忽略[113].當接觸變形量超過臨界彈性變形量 δy時,接觸行為進入彈塑性接觸階段,其接觸力與變形量之間的關(guān)系近似線性;當載荷繼續(xù)增加,接觸變形量進一步增大且超過臨界塑性變形量 δp,接觸行為進入完全塑性接觸階段,其接觸力與變形量之間的關(guān)系基本為線性.另外,準靜態(tài)彈塑性接觸力模型認為在彈性階段沒有能量耗散(其加載與卸載曲線完全重合),且接觸過程遵守Hertz 理論.當接觸變形量大于臨界彈性接觸變形量時[113],由于在壓縮過程中發(fā)生了彈塑性變形,以至于其加載曲線不再遵守Hertz 接觸理論[139](其接觸剛度也不再是Hertz 接觸剛度),而卸載過程依然符合Hertz 理論,最終準靜態(tài)彈塑性接觸力模型通過不同加卸載曲線構(gòu)成的封閉面積(如圖9)描述了接觸過程中的能量損耗,說明了接觸過程中能量耗散是由于接觸力改變了壓縮階段接觸體的本構(gòu)關(guān)系.以上的力學現(xiàn)象是表3 中所有準靜態(tài)彈塑性接觸力模型的共性特征(注意:表3 中所有模型設(shè)計的臨界接觸力與臨界變形量不盡相同,具體的表達式請參考相關(guān)文獻).

圖8 準靜態(tài)彈塑性接觸力模型[146]Fig.8 Quasi-static elastoplatsic contact force model[146]

圖9 一般彈塑性接觸力模型中力與變形量之間的關(guān)系Fig.9 Force-deformation relationship for a general elastoplastic contact force model

此外,雖然表3 中的彈塑性接觸力模型可精確捕捉彈塑性沖擊過程中碰撞力與變形量之間的關(guān)系,但是,需要在仿真計算過程中區(qū)分碰撞過程是處于壓縮階段還是恢復(fù)階段,如果處于壓縮階段,需要判定當前接觸變形量是否超過臨界彈性變形量或者臨界塑性變形量[149-150];以此為依據(jù),根據(jù)不同階段的本構(gòu)關(guān)系計算碰撞過程中的力學行為.

更重要的是,每一次彈塑性接觸的計算中均需要保存永久變形量與最大接觸變形量,為下次碰撞行為仿真做準備.如此繁瑣的計算流程尤其不適合計算多重壓縮和多重碰撞現(xiàn)象;也正因為復(fù)雜的計算流程導(dǎo)致該類模型雖然能準確描述彈塑性接觸過程,但沒有在多體系統(tǒng)動力學的商業(yè)軟件中被廣泛應(yīng)用.

3 關(guān)于恢復(fù)系數(shù)的討論

多體系統(tǒng)中碰撞行為導(dǎo)致的能量損耗不可避免[17],恢復(fù)系數(shù)作為衡量沖擊過程中能量耗散的一個全局度量指標被廣泛應(yīng)用于工程與理論研究中,它可描述不同機制導(dǎo)致的能量耗散,包括阻尼損耗、彈塑性變形以及沖擊波在接觸體內(nèi)的傳播[151].常見的恢復(fù)系數(shù)有四類[152]:Newton 恢復(fù)系數(shù),Poisson 恢復(fù)系數(shù),Stronge 恢復(fù)系數(shù)與Hedrih 恢復(fù)系數(shù).其中Newton 恢復(fù)系數(shù)根據(jù)碰撞過程中動量守恒,利用碰撞前后法向接觸相對速度之比衡量碰撞過程中的能量損耗,其表達式為

Poisson 恢復(fù)系數(shù)將碰撞過程分為壓縮與恢復(fù)兩個階段,在壓縮階段碰撞體之間的相對法向速度不斷減小直至為零,在恢復(fù)階段利用壓縮階段儲存的沖量使相對法向速度為零的接觸體分離,于是通過碰撞恢復(fù)階段與壓縮階段的沖量之比衡量碰撞過程中的能量損耗,其表達式為

其中F為接觸力;tc為壓縮結(jié)束時間;te為整個碰撞過程結(jié)束時間;ts為碰撞行為起始時間.

Stronge 恢復(fù)系數(shù)認為壓縮過程中系統(tǒng)的動能以應(yīng)變能的形式儲存在接觸體中,在恢復(fù)階段通過釋放應(yīng)變能將接觸體分離,因此利用碰撞前后吸收的應(yīng)變能來衡量碰撞過程中的能量損耗,其表達式為

其中 δe為碰撞后對應(yīng)的相對變形量; δc為碰撞過程中的最大相對變形量; δs為接觸開始的變形量.

Hedrih[153-154]在研究兩個沿著圓軌跡滾動小球的正碰撞動力學行為時,通過碰撞前后小球的相對角速度定義了恢復(fù)系數(shù),其表達式為

其中 ω(-)與 ω(+)為碰撞前后小球之間的相對角速度.

以上四種根據(jù)不同物理量定義的恢復(fù)系數(shù)均是為了衡量碰撞前后的能量耗散,在本質(zhì)上是等效的;除非在斜碰撞行為中考慮了摩擦效應(yīng)等因素.一般而言,當恢復(fù)系數(shù)等于1 時,被稱之為純彈性接觸,即整個接觸過程沒有能量耗散;當恢復(fù)系數(shù)等于0時,被稱之為完全塑性接觸,即碰撞行為完成后初始動能被完全耗散.

對于表1 與表2 中的連續(xù)接觸力模型而言,在利用該類模型計算碰撞行為時,必先確定恢復(fù)系數(shù)[3,155-156];而恢復(fù)系數(shù)的取值精度則決定了多體系統(tǒng)碰撞動力學的仿真精度.目前在利用該類模型計算碰撞行為時其恢復(fù)系數(shù)一般根據(jù)經(jīng)驗值或通過實驗確定,而為了便于計算則一般采用經(jīng)驗值.然而,恢復(fù)系數(shù)在單純表征能量耗散的情況下,未考慮碰撞引起的局部接觸變形與碰撞時間,無法描述碰撞過程中接觸力隨時間的變化歷程,在多體系統(tǒng)碰撞動力學的研究中受到限制[5];所以表1 與表2 中的連續(xù)接觸力模型采用遲滯阻尼因子來描述碰撞過程中的能量耗散(阻尼環(huán)包圍的面積如圖3).然而,由于含阻尼因子的連續(xù)接觸力模型的原型為Hertz 接觸模型,導(dǎo)致該類模型在高恢復(fù)系數(shù)(大于0.85)下有較好的計算精度,一旦恢復(fù)系數(shù)小于0.8 該類模型則不能準確描述碰撞過程中的能量耗散[4,84],降低了多體系統(tǒng)碰撞動力學的仿真精度.所以當利用該類模型計算彈塑性碰撞行為時(EDEM 軟件中的連續(xù)接觸力模型(見表2)用于顆粒物質(zhì)之間的彈塑性碰撞仿真),一般通過調(diào)節(jié)恢復(fù)系數(shù)的大小來平衡彈塑性碰撞過程中的能量耗散(大部分學者認為表1和表2 中的連續(xù)接觸力模型只能用于彈性碰撞).然而,通過阻尼因子來描述碰撞過程中能量耗散的連續(xù)接觸力模型,雖然可以通過調(diào)節(jié)恢復(fù)系數(shù)描述彈塑性碰撞過程中的能量耗散,但是Hertz 接觸剛度高估了彈塑性接觸階段的剛度系數(shù)[157];即從整體上基于連續(xù)接觸力模型依然不能準確獲得彈塑性碰撞行為的力學特征(第四節(jié)中將詳細說明該問題).

恢復(fù)系數(shù)不是關(guān)于接觸材料的常數(shù),而是描述碰撞過程中能量耗散的標志性參數(shù)[61,63,69,153,158];為了獲得準確的恢復(fù)系數(shù),眾多學者基于準靜態(tài)彈塑性本構(gòu)關(guān)系在不同的碰撞環(huán)境中推導(dǎo)了一系列的恢復(fù)系數(shù)具體表達式.其原理在于準靜態(tài)彈塑性模型可準確地獲得加載階段與卸載階段接觸力所做的功,所以可通過碰撞前后能量之比的平方根或者碰撞前后的速度之比來獲得準確的恢復(fù)系數(shù).其中包括Chang等[159]在CEB 彈塑性接觸力模型的基礎(chǔ)上利用接觸表面粗糙度的統(tǒng)計學模型推導(dǎo)了恢復(fù)系數(shù)的另一種解析解模型.Thornton 模型[135]通過線性化處理完全塑性階段的接觸力與變形之間的本構(gòu)關(guān)系,緩解了CEB 模型不連續(xù)的力學特性,并通過接觸力在碰撞前后的做功比獲得了封閉的恢復(fù)系數(shù)模型.Vu-Quoc等[147,160]在基于有限元法獲得本構(gòu)關(guān)系的基礎(chǔ)上分析了碰撞過程的恢復(fù)系數(shù),但并沒有給出恢復(fù)系數(shù)的封閉解.另外,Wu等[150]從碰撞速度的觀點基于有限元法定義了臨界速度,將恢復(fù)系數(shù)分為兩個部分描述碰撞過程中的能量損耗;Weir等[161]則推導(dǎo)了碰撞速度在低速環(huán)境下的恢復(fù)系數(shù)解析解;Jackson等[138,162-163]通過擬合有限元對恢復(fù)系數(shù)的分析結(jié)果優(yōu)化了基于碰撞速度的恢復(fù)系數(shù)精度.Ma等[113]利用其本構(gòu)方程獲得了完全塑性階段對應(yīng)的恢復(fù)系數(shù)(表4 中僅給出了6 種具有封閉解的恢復(fù)系數(shù)模型,其他形式的恢復(fù)系數(shù)模型可參考文獻[164]).從表4 中可以看出恢復(fù)系數(shù)不僅與碰撞體幾何與材料屬性有關(guān),而且與碰撞速度密切相關(guān).

表4 恢復(fù)系數(shù)模型Table 4 Coefficient of restitution(CoR) model

恢復(fù)系數(shù)作為衡量碰撞過程中能量損耗的直接參數(shù),其中連續(xù)接觸力模型在計算碰撞行為時需要依賴恢復(fù)系數(shù)的大小,而準靜態(tài)彈塑性接觸力模型可根據(jù)精確的本構(gòu)關(guān)系確定恢復(fù)系數(shù)的大小.因此,對同一種接觸行為,完全可以先通過準靜態(tài)彈塑性接觸力模型確定恢復(fù)系數(shù)的大小,并以此作為輸入至連續(xù)接觸力模型中,避免了根據(jù)經(jīng)驗值與實驗測試方式獲得恢復(fù)系數(shù)的弊端.在此基礎(chǔ)上,要將面臨一個關(guān)鍵的問題:連續(xù)接觸力模型中阻尼因子做功代表的能量耗散是否與準靜態(tài)彈塑性接觸力模型的碰撞前后接觸力做功之差描述的能量耗散具有內(nèi)在聯(lián)系,即阻尼環(huán)包圍的面積(如圖3)是否等于加卸載曲線包圍的面積(如圖9)？下面一節(jié)將以恢復(fù)系數(shù)為橋梁,詳細解釋兩類模型之間的內(nèi)在聯(lián)系.

4 兩類接觸力模型之間的聯(lián)系

考慮到表1 與表2 中連續(xù)接觸力模型在計算彈塑性接觸碰撞時,Hertz 接觸剛度高估了彈塑性接觸階段的接觸剛度;導(dǎo)致其連續(xù)接觸力模型不能準確獲得動態(tài)彈塑性碰撞行為的力學特征.為此,該部分將基于本課題組在2015 年提出的連續(xù)準靜態(tài)彈塑性接觸力模型[113](ML 模型)修正連續(xù)接觸力模型在計算彈塑性碰撞時的接觸剛度.

4.1 基于材料彈塑性本構(gòu)的等效遲滯阻尼

由于當接觸行為進入彈塑性階段時其接觸力與變形量之間的關(guān)系近似為線性,所以,借助ML 模型的本構(gòu)關(guān)系將彈塑性階段的剛度線性化為

基于該線性化的彈塑性接觸剛度根據(jù)線性彈簧-阻尼模型,假設(shè)其考慮能量耗散的彈塑性接觸力模型的形式為

其中χp為新的遲滯阻尼因子.

根據(jù)1.1 節(jié)中阻尼因子的推導(dǎo)方式,在壓縮過程的結(jié)束時刻,最大應(yīng)變能為

對阻尼項積分就可獲得整個接觸過程的能量耗散

阻尼項將接觸過程耗散的能量均勻分布在壓縮和恢復(fù)階段,根據(jù)式(48),其整個接觸過程的能量耗散

以及該時刻系統(tǒng)動量守恒式(13),將式(49)改寫為

聯(lián)立式(8)、式(49)與式(51),新的阻尼因子為[43]

因此,新的連續(xù)彈塑性接觸力模型為

很明顯,由于該阻尼因子的推導(dǎo)過程利用了系統(tǒng)的能量與動能守恒原則,以至于其阻尼項的分母中含有初始碰撞速度,即該阻尼項為遲滯阻尼;但并不妨礙計算碰撞過程中的能量耗散.

假設(shè)一個間隙球面關(guān)節(jié),其接觸參數(shù)如表5 所示,其初始碰撞速度為4 m/s,此時最大接觸變形量(111.52 μm)已超過臨界彈性變形量(1.578 8 μm),接觸過程處于彈塑性接觸階段;此時的恢復(fù)系數(shù)可由ML 模型識別為0.6909.將該恢復(fù)系數(shù)作為新連續(xù)接觸力模型的輸入量可計算彈塑性碰撞行為的力學特征如圖10.

圖10 接觸力-變形量(關(guān)于彈塑性階段)[43]Fig.10 Relationship between the force and deformation(in elastoplastic phase)[43]

表5 接觸參數(shù)Table 5 Contact parameters

通過相對誤差分析,新的連續(xù)接觸力模型中阻尼因子代表的能量耗散與ML 模型描述的能量耗散誤差不超過0.25%[43],即阻尼環(huán)的面積與三角形的面積基本一致.該分析結(jié)論表明:當恢復(fù)系數(shù)保持一致時,兩類模型在描述彈塑性碰撞行為的能量耗散時是等效的,說明了在彈塑性發(fā)生時,人為阻尼項表示的能量耗散就是彈塑變形做的功,與阻尼項的積分路徑無關(guān).由于兩類模型能保證碰撞過程中能量耗散的一致性,所以兩類模型在計算最大接觸力和最大變形量方面雖然不能完全一致(畢竟兩類模型的數(shù)學表達式和計算原則均不相同),但整體上可保持協(xié)調(diào),尤其是碰撞后接觸體的運動狀態(tài)基本一致;主要歸因于新的連續(xù)彈塑性接觸力模型式(52)利用線彈塑性接觸剛度代替了Hertz 接觸剛度,避免了Hertz 接觸剛度與彈塑性接觸剛度不符的缺陷;更重要的是由于兩類模型在描述碰撞過程中的能量耗散時完全等效.

因此,相比表1 與表2 中的連續(xù)接觸力模型,新的彈塑性接觸力模型可精確地計算彈塑性碰撞過程,彌補了Hertz 接觸剛度引入的誤差;另外,相比表3 中的準靜態(tài)彈塑性接觸力模型,新的連續(xù)接觸力模型不僅在力學特征上能與準靜態(tài)彈塑性接觸力模型保持一致,并且簡化了準靜態(tài)彈塑性接觸力模型在計算彈塑性碰撞過程中的復(fù)雜流程,避免了在計算過程中區(qū)分碰撞行為處于壓縮或恢復(fù)階段,更不必保存每一次碰撞行為導(dǎo)致的永久與最大變形量;而只需要判定接觸行為是否發(fā)生.然而,該模型由于遲滯阻尼項中分母含有初始速度導(dǎo)致該模型在計算顆粒物質(zhì)之間的碰撞接觸行為時會導(dǎo)致數(shù)值奇異問題(因為大部分顆粒物質(zhì)在初始階段均保持靜止狀態(tài),顆粒之間初始相對碰撞速度為零).

4.2 基于材料彈塑性本構(gòu)的等效黏性阻尼

為了避免上述所提的連續(xù)接觸力模型分母中初始相對碰撞速度導(dǎo)致的數(shù)值奇異問題[88],該小節(jié)依然基于線性化的彈塑性接觸剛度,將彈塑性碰撞過程等效為線性的振動行為(如圖4),該模型的動力學方程為單自由度非受迫阻尼自由振動方程式(23),根據(jù)該方程的解析解與初始碰撞條件式(24)～ 式(30)可推導(dǎo)系統(tǒng)的阻尼系數(shù)為[94]

所以該連續(xù)彈塑性接觸力模型為

在該模型的推導(dǎo)過程中沒有利用碰撞前后的能量與動能守恒原則,避免了獲得的阻尼項分母中含初始相對碰撞速度的特征,即阻尼項為黏性阻尼;成功回避了連續(xù)接觸力模型式(53)的數(shù)值奇異問題.

關(guān)于該模型的仿真精度,本文以一維球鏈為研究對象(如圖11),以實驗數(shù)據(jù)為基準,對比分析了EDEM 軟件中連續(xù)接觸力模型(見表2 中Tsuji等[96]模型)與新接觸力模型式(55)之間的精確度.相關(guān)仿真參數(shù)與實驗數(shù)據(jù)詳見文獻[117].

圖11 一維球鏈的碰撞實驗裝置示意圖[117]Fig.11 Experimental setup of the one-dimension granular chain[117]

圖12 給出了新連續(xù)接觸力模型式(55)在計算顆粒物質(zhì)動態(tài)碰撞行為時的計算結(jié)果,孤立波在一維球鏈中的傳播特性被很好地復(fù)現(xiàn),其結(jié)果被實驗數(shù)據(jù)證明是合理的.當與EDEM 軟件中的連續(xù)接觸力模型對比分析時,圖13 中展現(xiàn)出兩種模型基本能反映孤立波在一維球鏈中的傳播特性,但是在孤立波波峰對應(yīng)的最大接觸力方面,兩者與實驗數(shù)據(jù)的誤差百分比見表6;可以清楚地看出,新的連續(xù)接觸力模型在計算最大接觸力方面精度明顯高于EDEM軟件中被使用的連續(xù)接觸力模型[94],再一次證明了Hertz 接觸剛度在描述彈塑性接觸行為的剛度屬性時存在誤差.

圖12 一維球鏈中孤立波的傳播特征Fig.12 Propagation of solitary wave in the one-dimension granular chain

圖13 連續(xù)接觸力模型之間的對比分析Fig.13 Comparative analysis between two continuous contact force models

表6 誤差分析對比Table 6 Comparative analysis of the error percentage

5 接觸力模型在多體系統(tǒng)動力學中的應(yīng)用

多體系統(tǒng)中最早采用剛體模型基于沖量動量法研究碰撞問題,認為應(yīng)力波在碰撞體中的傳播速度無限大,幾乎同時影響了碰撞體各質(zhì)點速度;引起了多體系統(tǒng)廣義速度的跳躍,造成在考慮切向摩擦的碰撞動力學方程求解中往往出現(xiàn)無解或者多解的情況[115];并在整個碰撞過程中,采用恢復(fù)系數(shù)描述碰撞前后速度的變化以及能量耗散.然而,該方法并不能獲得碰撞力隨時間的變化歷程;要正確處理碰撞中的摩擦問題,與分析接觸力隨時間的變化規(guī)律,需要利用連續(xù)接觸力模型考慮碰撞體之間的微運動規(guī)律,解決由于切向摩擦引起的數(shù)值不穩(wěn)定問題,并通過阻尼因子闡述碰撞過程中的能量耗散.另外,連續(xù)接觸力模型從數(shù)學的角度為碰撞問題提供了一種簡單易行的方法,并且認為碰撞不再是瞬時過程,可以精細研究碰撞過程中碰撞力與侵入量以及碰撞速度之間的變化關(guān)系,這有利于研究由于碰撞而引起的結(jié)構(gòu)破壞現(xiàn)象[166],這也是連續(xù)接觸力模型在工程領(lǐng)域被廣泛應(yīng)用的原因.然而,連續(xù)接觸力模型忽略了碰撞引起的應(yīng)力波在碰撞體中的傳播效應(yīng)[115].實際多體系統(tǒng)中的碰撞問題不僅會在碰撞位置產(chǎn)生局部的彈塑性變形,而且碰撞引起的局部應(yīng)力集中會以應(yīng)力波的形式在碰撞體中傳播;為了同時描述碰撞行為引起的局部彈塑性變形以及應(yīng)力波的傳播規(guī)律,必須借助連續(xù)介質(zhì)力學研究多體系統(tǒng)中的碰撞接觸行為.

以一維球鏈的碰撞現(xiàn)象為例[117](如圖11),可以反映利用上述單點接觸力模型在計算多體系統(tǒng)中碰撞問題時遇到的困難.圖11 中的一維球鏈,采用彈簧的方式局部柔化剛體小球之間的接觸碰撞問題,該球鏈的接觸碰撞現(xiàn)象不同于兩個小球之間的單點接觸碰撞行為,球鏈中的碰撞現(xiàn)象不僅涉及多點碰撞,而且涉及接觸過程中的多重壓縮與多重碰撞現(xiàn)象.在第一對小球碰撞結(jié)束之前,第二對小球已經(jīng)進入壓縮階段并發(fā)生多點碰撞現(xiàn)象,以此類推,第一對小球發(fā)生碰撞激發(fā)的碰撞力以孤立波的形式在球鏈中傳遞并引發(fā)多重碰撞現(xiàn)象;當最后一對小球產(chǎn)生碰撞并將碰撞力重新傳遞到開始的碰撞位置時將發(fā)生多重壓縮現(xiàn)象,并且除第一次碰撞外,因為能量耗散導(dǎo)致其他碰撞激發(fā)的動態(tài)響應(yīng)逐漸衰減[167].被撞擊的一維球鏈中孤立波在球鏈中傳播與反射引發(fā)多重碰撞現(xiàn)象,當波反射到碰撞開始位置時可能在單個接觸點上發(fā)生多重壓縮現(xiàn)象,與此同時碰撞過程中伴隨著能量耗散以至于碰撞現(xiàn)象越來越微弱.

因此,當柔性體之間發(fā)生碰撞行為時,首先在碰撞的局部區(qū)域發(fā)生彈塑性變形,并產(chǎn)生高幅值且急劇變化的碰撞力[115],同時激發(fā)應(yīng)力波在柔性體中傳播,以及會誘發(fā)應(yīng)力波與碰撞之間相互耦合作用等一系列復(fù)雜的動態(tài)行為[116].多體系統(tǒng)中的能量一部分由于局部接觸區(qū)域發(fā)生的彈塑性變形被耗散,另一部分能量被傳遞到碰撞體中以瞬態(tài)應(yīng)力波的形式在柔性體中傳播與反射,以至于影響柔性體碰撞后的運動狀態(tài);另外,在碰撞結(jié)束前,應(yīng)力波在柔性體中被反射到接觸位置與碰撞行為發(fā)生耦合現(xiàn)象,進而誘發(fā)多重壓縮與多重碰撞等現(xiàn)象,進一步影響接觸力的大小.除此之外,在此期間需要對多次微碰撞過程中的碰撞力峰值、碰撞時間及碰撞次數(shù)進行準確預(yù)測,錯綜復(fù)雜的碰撞與柔性變形耦合關(guān)系增加了多體系統(tǒng)碰撞動力學的建模仿真難度[168];更進一步,上述復(fù)雜的動態(tài)交互行為令碰撞行為與柔性變形之間的耦合關(guān)系變得更加撲朔迷離.圍繞上述復(fù)雜耦合問題國內(nèi)外眾多學者開展了系統(tǒng)的研究[41,169],按照碰撞過程中是否考慮彈塑性變形將以下研究分為兩類:(1)基于動態(tài)連續(xù)接觸力模型的多體系統(tǒng)碰撞動力學;(2)基于準靜態(tài)彈塑性接觸模型的多體系統(tǒng)碰撞動力學.

5.1 未考慮彈塑性接觸變形的碰撞動力學研究

在忽略碰撞行為中彈塑性變形的情況下,眾多學者對多體系統(tǒng)中碰撞現(xiàn)象與柔性變形的相互作用進行了系統(tǒng)的研究.文獻[170-171]利用浮動坐標法建立曲柄滑塊機構(gòu)中柔性連桿的動力學模型,采用動態(tài)接觸力模型計算剛性滑塊之間的碰撞力,分析了剛?cè)狁詈吓c碰撞行為之間的動態(tài)關(guān)系.文獻[114,172]與胡海巖等采用絕對節(jié)點坐標法分析了多體系統(tǒng)中間隙關(guān)節(jié)碰撞現(xiàn)象與構(gòu)件變形的特征,研究了柔性機器人在抓取過程中的接觸碰撞問題.考慮到絕對節(jié)點坐標在描述變形體時引入了大量的彈性坐標,文獻[173]采用Craig-Bampton 子結(jié)構(gòu)模態(tài)綜合法縮減了彈性坐標數(shù)量,分析了間隙關(guān)節(jié)元素之間碰撞對柔性系統(tǒng)的影響.張云清等[174]利用浮動坐標法分析了柔性構(gòu)件小變形與碰撞行為的耦合現(xiàn)象,建立了曲柄滑塊機構(gòu)柔性連桿的動力學模型,討論了關(guān)節(jié)剛度與碰撞現(xiàn)象對系統(tǒng)振動頻率的影響.文獻[175-176]采用絕對節(jié)點坐標法建立柔性桿與剛性孔之間的動力學模型,考慮了柔性桿與剛性孔之間的碰撞與摩擦效應(yīng)對系統(tǒng)的動態(tài)影響;或者采用假設(shè)模態(tài)法建立柔性體的動力學模型,在彈性范圍內(nèi)采用動態(tài)接觸力模型計算柔性桿與剛體之間的碰撞力,揭示了沖擊波的傳播效應(yīng).劉昊等[177]采用絕對節(jié)點坐標法建立了空間充氣展開繩網(wǎng)捕獲系統(tǒng)的動力學模型,基于Hertz 接觸理論分析了被捕獲物體與柔性繩網(wǎng)之間的碰撞現(xiàn)象.方建士等[178]利用子系統(tǒng)法考慮柔性梁的動力剛化效應(yīng),基于考慮能量耗散的Hertz 接觸力模型研究了大范圍運動柔性梁與剛性地面之間的接觸碰撞現(xiàn)象.彭海軍等[179]針對柔性多體系統(tǒng)中接觸體之間的碰撞與摩擦問題,基于辛離散化提出了一種非光滑接觸方法,避免了因互補變量過多而導(dǎo)致求解效率較低的困難.虞磊等[180]基于絕對節(jié)點坐標法建立了柔性梁與柔性板的動力學模型,將柔性體之間的碰撞檢測轉(zhuǎn)化成基本幾何體的碰撞檢測,并利用Hertz 接觸理論計算了柔性體與剛體以及柔性體之間的碰撞力,仿真結(jié)果與RecurDyn 對比驗證了該方法的正確性.

上述工作主要集中研究碰撞行為對多體系統(tǒng)中變形構(gòu)件的動態(tài)影響[181],在彈性范圍內(nèi)分析了碰撞行為與柔性變形之間的耦合關(guān)系,對柔性體與碰撞行為之間的精細化建模與求解方面做出了貢獻.但是較少關(guān)注多體系統(tǒng)中能量在經(jīng)歷碰撞與柔性變形之后的重新分配與耗散問題,實際上多體系統(tǒng)中碰撞與柔性變形的耦合關(guān)系本質(zhì)就是能量相互轉(zhuǎn)換與耗散的過程[182-183];同時,現(xiàn)有工作較少關(guān)注碰撞與應(yīng)力波之間引發(fā)的多重壓縮與多重碰撞問題.

5.2 考慮彈塑性接觸變形的碰撞動力學研究

眾多國內(nèi)外學者發(fā)現(xiàn)多體系統(tǒng)在高速碰撞下,如果忽略碰撞行為中的彈塑性變形將導(dǎo)致其系統(tǒng)動態(tài)響應(yīng)與實際情況存在偏差[184-186].Yigit[187]研究了柔性梁與剛性地面的碰撞動力學,利用完全彈性接觸力模型計算了柔性梁與剛體之間碰撞力,仿真結(jié)果顯示完全彈性接觸模型并不能獲得精確的結(jié)果,通過實驗證明彈塑性接觸力模型更加符合實際情況.驀朋波等[116]利用動態(tài)子結(jié)構(gòu)模態(tài)綜合法推導(dǎo)了柔性構(gòu)件之間碰撞激發(fā)瞬態(tài)波傳播的動態(tài)子結(jié)構(gòu)模型,在保持接觸剛度不變的情況下,在碰撞過程中考慮了永久的彈塑性變形,基于單軸壓縮(UC)模型計算了接觸體之間的碰撞力,分析了彈塑性波的傳播特性.文獻[111,188]圍繞柔性體與剛體碰撞的動力學做了系統(tǒng)的研究,主要采用浮動坐標法建立柔性構(gòu)件的動力學模型,分別利用動態(tài)連續(xù)接觸力模型與單軸壓縮(UC)模型計算柔性體與剛體之間的接觸力,考慮了正碰撞與斜碰撞過程中的摩擦現(xiàn)象,分析了接觸在發(fā)生彈性變形與塑性變形時系統(tǒng)的動態(tài)特性,指出了動態(tài)連續(xù)接觸力模型與靜態(tài)彈塑性模型不僅對能量耗散的描述不同,而且通過兩種模型計算的接觸力大小也大相徑庭.Chen等[112]采用模態(tài)綜合法建立柔性體的動力學模型,提出接觸區(qū)域判定方法,通過罰函數(shù)的懲罰因子計算柔性體之間的彈塑性接觸力,對比分析了彈性和彈塑性碰撞后系統(tǒng)的動能相差大約30%,說明了考慮碰撞中彈塑性變形的必要性與能量耗散的主要來源.王檢耀等[168]與洪嘉振基于有限元法同時采用罰函數(shù)法和附加約束法的接觸模型計算了柔性桿之間的接觸碰撞問題,發(fā)現(xiàn)柔性梁與碰撞之間的相互作用會引發(fā)多次微碰撞問題.

上述研究考慮了碰撞行為中發(fā)生的彈塑性變形,證實了彈塑性變形在碰撞行為中較為常見且不可忽視,主要分析了彈塑性碰撞行為引起的應(yīng)力波在柔性體中的傳播規(guī)律,研究了碰撞與柔性變形之間耦合作用誘發(fā)的多重碰撞問題.然而,在碰撞力的計算過程中對于彈塑性變形引起的接觸剛度變化缺乏精細化的建模與分析(實際的接觸剛度分為彈性、彈塑性與塑性接觸剛度)[189],導(dǎo)致其碰撞中能量耗散的預(yù)測出現(xiàn)偏差,影響了碰撞的接觸時間與接觸力的大小,干擾了柔性體中應(yīng)力波與碰撞行為之間的耦合作用.另外,利用準靜態(tài)彈塑接觸力模型或者懲罰因子計算碰撞體之間的接觸力,其中懲罰因子難以表征碰撞體之間的真實接觸性質(zhì);其次,由于準靜態(tài)彈塑性接觸模型對接觸行為的描述完全取決于相對接觸變形量的大小;因此,基于該類模型計算接觸行為容易受到積分誤差的干擾.相反,連續(xù)接觸力模型從相對接觸變形量與相對變形速度上同時控制能量耗散與接觸力的大小[190-191],對積分誤差不敏感,這也是連續(xù)動態(tài)彈塑性接觸力模型的另一種優(yōu)勢.

多體系統(tǒng)碰撞動力學中碰撞與柔性體耦合效應(yīng)容易引起多重壓縮與多重碰撞現(xiàn)象,其中柔性體從碰撞行為中吸收的能量與耗散能量直接相關(guān)[192];而應(yīng)力波的傳播規(guī)律取決于柔性體吸收的能量[193],其碰撞時間也依賴于能量耗散的效率.所以,正確表征碰撞過程中的能量耗散與揭示能量耗散的本質(zhì)原因是多體系統(tǒng)碰撞動力學的核心內(nèi)容,只有保證精確計算碰撞過程中能量耗散在時間和空間上的傳播規(guī)律,才能從本質(zhì)上揭示碰撞與柔性變形之間的耦合關(guān)系.

6 總結(jié)與發(fā)展趨勢

(1) 本文主要闡述了不考慮碰撞過程中應(yīng)變率和硬度發(fā)生變化的正向接觸力模型的研究進展,分析了兩類不同接觸力模型的發(fā)展過程,指出了連續(xù)接觸力模型在彈塑性碰撞動力學計算中存在的問題,討論了準靜態(tài)彈塑性接觸力模型在計算彈塑性碰撞行為時復(fù)雜的計算流程.明確了Hertz 接觸剛度是導(dǎo)致連續(xù)接觸力模型在計算彈塑性碰撞行為的主要誤差來源,以此為依據(jù),以線性的彈塑性接觸剛度為基礎(chǔ),通過碰撞過程中的能量守恒原則推導(dǎo)了新的遲滯阻尼因子,證實了人為阻尼因子代表的能量損耗與彈塑性碰撞中彈塑性變形做的功等效,而與阻尼項的積分路徑無關(guān).

(2) 當前的連續(xù)接觸力模型與準靜態(tài)彈塑性接觸力模型均是由法向碰撞與接觸行為發(fā)展而來,忽略了碰撞行為中的切向摩擦效應(yīng),以至于在分析斜碰撞行為時,需要先根據(jù)接觸力模型計算法向接觸力,再根據(jù)摩擦模型計算接觸體之間的切向接觸力,整個過程忽略了法向接觸力與摩擦行為之間的耦合關(guān)系;因此,未來的接觸力模型的發(fā)展要從斜碰撞入手,發(fā)展考慮法向和切向耦合作用的連續(xù)接觸力模型.

(3)載人返回艙著水/著陸以及航空發(fā)動機軸承的失效問題均涉及碰撞體與流體之間的流固耦合關(guān)系,為滿足國家重大航天發(fā)展需求,結(jié)合Reynolds 方程將接觸體之間的流體引入至連續(xù)接觸力模型的阻尼因子中;因此,充分考慮流體動力學對碰撞行為的影響,建立可描述流體與法向碰撞力相耦合的連續(xù)接觸力模型是未來的發(fā)展趨勢.

(4) 考慮到高速重載與高精度是未來多體系統(tǒng)的發(fā)展趨勢,間隙關(guān)節(jié)導(dǎo)致的接觸碰撞與磨損現(xiàn)象將更加突出,然而當前對于間隙關(guān)節(jié)碰撞行為的研究均忽略了高速重載下導(dǎo)致接觸體局部的彈塑性變形,給間隙關(guān)節(jié)元素之間的碰撞力計算引入了較大的誤差;另一方面,在基于Archard 磨損模型預(yù)測間隙關(guān)節(jié)元素之間磨損深度時,也將引起接觸應(yīng)力的計算誤差,最終影響間隙關(guān)節(jié)磨損表面的幾何屬性重構(gòu).因此,考慮間隙關(guān)節(jié)局部彈塑性變形的多體系動力學性能分析與壽命預(yù)測是未來的發(fā)展趨勢.

(5) 綜合考慮多體系統(tǒng)碰撞動力學中的多種非理想因素,包括多體系統(tǒng)中構(gòu)件的柔性變形、間隙關(guān)節(jié)的碰撞行為、關(guān)節(jié)的潤滑與摩擦效應(yīng),以及碰撞行為引起的局部彈塑性變形;建立保真度更高的多體系統(tǒng)碰撞動力學模型,分析局部彈塑性變形與系統(tǒng)柔性變形之間的耦合關(guān)系,探索應(yīng)力波在柔性體中的傳播特征;揭示多體系統(tǒng)關(guān)節(jié)中流固耦合效應(yīng)與系統(tǒng)大范圍非線性運動之間的內(nèi)在聯(lián)系,開發(fā)具有通用性多體系統(tǒng)碰撞動力學軟件將是未來面臨的挑戰(zhàn)之一.