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

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

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

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

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

??(?？巳卮髮W工程、數學和物理科學學院,英國?？巳?EX44 QF)

引言

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

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

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

1 等效連續接觸力模型

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

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

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

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

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

1.1 等效遲滯阻尼

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

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

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

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

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

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

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

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

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

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

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

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

其次,根據碰撞前后動量守恒原則有如下關系式

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

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

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

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

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

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

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

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

另外,壓縮結束階段的應變能也可通過彈性力做功表示

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

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

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

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

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

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

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

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

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

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

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

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

1.2 等效黏性阻尼

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

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

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

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

需要注意的是:參數KL是線性彈簧的剛度,而不是Hertz 接觸剛度,所以方程(23)中彈性力項KLδ的變形量指數為1,其具體原因將在4.1 小節解釋接觸力模型中剛度系數與指數系數關系的重要性.

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

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

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

整個碰撞過程的持續時間為

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

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

最后,根據Newton 恢復系數的定義

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

Maxwell 模型的數學表達式為

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

其中 α為恢復系數的函數.

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

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

通過對上式積分可以獲得相對碰撞速度,以此導出阻尼系數

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

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

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

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

1.4 能量指數與接觸剛度之間的關系

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

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

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

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

1.5 關于等效連續接觸力模型的討論

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

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

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

2 準靜態接觸力模型的研究現狀

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

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

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

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

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

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

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

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

表3 連續準靜態彈塑性接觸力模型Table 3 Continuous quasi-static elastoplastic contact force models

續表3

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

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

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

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

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

3 關于恢復系數的討論

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

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

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

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

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

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

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

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

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

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

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

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

4 兩類接觸力模型之間的聯系

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

4.1 基于材料彈塑性本構的等效遲滯阻尼

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

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

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

根據1.1 節中阻尼因子的推導方式,在壓縮過程的結束時刻,最大應變能為

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

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

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

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

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

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

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

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

表5 接觸參數Table 5 Contact parameters

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

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

4.2 基于材料彈塑性本構的等效黏性阻尼

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

所以該連續彈塑性接觸力模型為

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

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

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

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

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

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

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

5 接觸力模型在多體系統動力學中的應用

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

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

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

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

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

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

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

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

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

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

6 總結與發展趨勢

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

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

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

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

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