柔性梁含摩擦斜碰撞的剛體-彈簧-質點混合模型研究*

沈煜年， 顧金紅

(南京理工大學力學與工程科學系， 江蘇 南京 210094)

?

柔性梁含摩擦斜碰撞的剛體-彈簧-質點混合模型研究*

沈煜年， 顧金紅

(南京理工大學力學與工程科學系， 江蘇 南京 210094)

為精確快速地預測柔性梁因斜碰撞導致的復雜瞬態動力響應，提出了一種可同時計及局部接觸區法向接觸柔度、切向接觸柔度以及柔性梁整體結構柔度的混合分析模型(HAM)?；谟邢薅畏ǖ乃枷?，將柔性梁的整體位移場離散為彈簧-阻尼-剛體系統，局部接觸區的法向雙線性壓縮-恢復以及切向摩擦變形過程則用含有能量恢復系數的彈塑性彈簧-質點系統進行描述。依據HAM模型，導出了斜碰撞系統在不同接觸狀態下的分段連續動力學方程。給出了法向壓縮-恢復狀態和切向黏-滑運動狀態的切換準則，并運用事件驅動法對數值算例進行了求解。通過比較HAM模型解、實驗數據以及商業有限元解，表明該模型能準確預測接觸力的時間歷程并描述反復切換的黏-滑(反向)運動模式，驗證了HAM模型的可行性。

非線性振動； 斜碰撞； 切向柔度； 摩擦接觸； 黏-滑運動

引　言

斜碰撞不僅存在于宏觀大尺度宇宙天體之間[1]，而且在機械工程、航空航天甚至微尺度機械領域中也是頻繁出現的現象[2]。例如空間探測器的交會對接[3]、飛行器表面的粒子侵蝕[4]、機械手抓取貨物的沖擊作用[5]、齒輪對的相互拍擊[6]、金屬的爆炸焊接工藝[7]以及微小管道機器人的管內運動[8]等。少數情形人們會利用含摩擦斜碰撞去實現機構特定的功能目的；但更多的情況則是斜碰撞會帶來結構失效、增加噪音以及安全性降低等許多不利的一面[9-12]。

斜碰撞發生時，碰撞結構不僅在接觸面的法向存在一個單邊約束，而且在切向還有一個摩擦約束[13]。相比正碰正碰撞，斜碰撞的瞬態特征更加復雜，切向約束會導致接觸面間的黏滯-微滑動(有時會反向)[14]。Johnson[15]在研究具有初始角速度的超彈性球體斜向撞擊水平地面時，發現了一個運用剛體碰撞模型無法解釋的特殊現象，即球質心速度的水平分量以及轉動角速度在斜碰撞后均會發生反向。通過分析，其認為接觸區較大的切向柔度可能是該現象產生的根本原因。Stronge[16]和Shen[14]考慮硬質碰撞物體接觸區的切向柔度，非接觸區部分仍視為剛體，建立能分析接觸點反向滑動的集中參數模型，討論了接觸柔度對接觸點的切向速度和接觸力影響。然而，斜碰撞事件中接觸與分離的狀態切換(即接觸約束的添加和刪除)會使桿和梁等可變形體的振動模態發生突變，且接觸區產生彈塑性變形。這些特征均與結構的柔性密切相關，目前尚缺乏深入研究。

本文將前述文獻的研究進行了發展，提出了一種可同時計及局部接觸區法向接觸柔度、切向接觸柔度以及柔性梁整體結構柔度的混合分析模型(HAM)，用以研究柔性梁的含摩擦斜碰撞問題。基于有限段法思想，將柔性梁的整體位移場離散為彈簧-阻尼-剛體系統。局部接觸區的法向雙線性壓縮-恢復過程以及切向摩擦變形過程用含有能量恢復系數的彈塑性彈簧-質點系統進行描述。根據所建混合模型，推導了斜碰撞系統在不同接觸狀態下的分段連續動力學方程，給出了法向壓縮-恢復狀態和切向黏-滑運動狀態的切換準則，并運用事件驅動法對數值算例進行了求解。驗證了本文模型的可行性。

1　斜碰撞混合分析模型

如圖1(a)所示，考慮一個端部為半圓頭的圓形截面柔性梁以初始平動速度V斜向撞擊粗糙剛性地面。梁長為L、橫截面積為A、材料的楊氏模量為E、質量密度為ρ、泊松比為ν。梁的軸線與x軸正向夾角為u。首先基于有限段思想將柔性梁整體分割成若干離散段[17](如圖1(b)所示)，建立梁整體的彈簧-阻尼-剛體系統，然后引入考慮切向接觸柔度的彈塑性彈簧-質點系統用以描述接觸效應和計算接觸力。通過聯合前述兩個子系統，可建立斜碰撞剛體-彈簧-質點混合模型。

1.1柔性梁整體的彈簧-阻尼-剛體系統

如圖1(b)所示，將柔性梁等分成n段，每段長度和質量分別為l=L/n和m=ρAl。為了同時考慮梁的軸向變形和彎曲變形，相鄰的離散段之間通過兩根關于中心軸對稱分布的、彈簧剛度為K的彈簧阻尼器聯接，彈簧距中心軸的距離為a/2。梁的軸向彈性通過選擇適當的拉壓剛度K值確定，梁的彎曲剛度由彈簧距離a來確定。梁端部與粗糙表面的接觸效應由1.2節中的局部接觸模型進行描述。該模型可以計算斜碰撞過程中的法向接觸力F3和切向摩擦力F1，以及刻畫黏滯-微滑動運動模式。

圖1(b)所示離散系統具有2n+1個自由度，用來描述其運動狀態的廣義坐標列向量可表示為

X=(x0,y0,u0,q1,u1,u2,…,qn-1,un-1)T

圖1　柔性梁斜向碰撞粗糙剛性地面Fig.1　Flexible beam oblique impact against rough rigid ground

式中x0和y0分別為接觸點在笛卡爾坐標系x-o-y中的位置坐標，ui(i=0,1，…，n-1)為第i+1個離散段與x軸的夾角。qi為第i+1個離散段離散截面圓心到第i個離散段離散截面的垂直距離。第i個離散段的質心點的位置(xic,yic)可表示為

(1)

離散系統的總動能T為

(2)

式中Jc為每小段對其質心的轉動慣量。由于梁碰撞引起的振動衰減發生在幾秒內，而碰撞發生在幾百個微秒之內，因此本文忽略梁內部阻尼對碰撞瞬態響應的影響[17]。系統總勢能V(取y=0處為勢能零點)為

(3)

作用于系統的廣義力向量Q由非保守力Q′(即接觸端的接觸力)和保守力向量Q″(即重力)疊加而成，系統的運動可由含有接觸力的第二類拉格朗日方程表示為

(4)

將公式(2)和(3)代入方程(4)，由符號運算軟件得到2n-1個方程組成的非線性動力學方程組

(5)

式中M為系統總體質量陣，B為關于廣義位移、廣義速度和外力的列向量。下面運用胡克定律和歐拉伯努利方程估算離散段之間彈簧單元的剛度K和a。

(1)離散模型的軸向變形

假設梁在自身重力的單獨作用下，該梁垂直懸掛于空間中。此時均質連續梁的質心處產生的軸向靜位移δc為

(6)

同時，對于離散梁模型，該點的位移δd又可用梁上半部分的段間彈簧變形量計算得到

(7)

由δc=δd，可解得拉壓剛度K

(8)

(2)離散模型的彎曲變形

根據連續懸臂梁的彎曲理論，可獲得在重力單獨作用下懸臂梁質心位置的轉角為

(9)

式中Iz為橫截面對z軸的慣性矩(對圓截面的梁有Iz=πd4/64)。由于是小變形問題，離散梁模型中第i段的力矩平衡方程為

(10)

式中φi為第i-1段與第i段之間的相對轉角。方程(10)中的第1項是由重力引起的，第2項是第i+1段和第i段之間彈簧的作用力矩。此外，第i+1段的絕對轉角ui可寫成如下形式

(11)

將式(10)和(11)聯立，可獲得中心段的轉角為

(12)

由uc=ud，得到

(13)

式中a只與段數n和梁橫截面參數有關。

1.2局部接觸區的彈塑性彈簧-質點系統

本文分別采用雙線性柔性單元和線性柔性單元描述接觸區的法向彈塑性變形效應和切向變形效應[14]，建立了如圖2所示局部接觸區的彈塑性彈簧-質點系統。

圖2　平面含摩擦斜碰撞局部彈性接觸模型Fig.2　Hybrid analytical model for planar impact with friction

(14)

(15)

圖3　法向和切向柔性單元的力與變形量的關系Fig.3　Force-deflection relation for normal and tangential compliant elements

此運動過程中法向接觸力和切向接觸力滿足Coulomb摩擦定律

(16)

2　接觸狀態的判斷標準

斜碰撞問題除了會導致法向存在一個壓縮和恢復階段的切換，在切向還存在滑動和黏滯狀態的切換。在使用事件驅動法求解碰撞響應時，需要判斷接觸所處的狀態。假設切向由黏滯切換為滑動的時刻為t01，從滑動切換為黏滯的時刻為t10。各狀態的判斷標準如下：

(1)初始狀態

在碰撞初始時刻(t=0)，法向接觸處于壓縮階段。切向初始黏滯的條件為

(17)

反之，一旦入射角大于μη2，則切向的初始狀態為滑動。

(2)壓縮或恢復階段

在碰撞的初始階段，接觸為壓縮階段。壓縮階段結束和恢復階段開始于tc時刻，此時法向相對速度消失，即滿足

(18)

(3)黏滯或滑動

(19)

相反地，如果系統之前處于滑動狀態，當質點P和接觸點C′間的相對速度消失，變為零時，此時記為t10時刻，系統開始由滑動狀態轉變為黏滯狀態，即

(20)

(4) 碰撞結束

若在tf-ε時刻法向接觸力為正(即F3(tf-ε)>0)，在tf+ε時刻法向接觸力為正(即F3(tf+ε)<0)，其中ε是一個正數小量，則認為碰撞結束的時刻為tf，法向接觸力將會消失為零，即

(21)

3　算　例

若無特殊說明，算例的系統參數為：柔性梁的長度為1 m，圓形橫截面面積為10-4m2，梁的接觸端部為半圓頭，材料為鋼(即楊氏模量為210 GPa，質量密度為7 800 kg/m3，泊松比為0.3)。梁的法向初始平動速度y0=-1 m/s，切向初始平動速度為x0=3.5 m/s。梁初始時刻與剛性平面接觸，無相對轉動。梁的軸線與x軸正向夾角為3π/4，即u0=3π/4。

Johnson[16]認為在法向作用力P的準靜態作用下接觸圓面的半徑a可表示為(3PR/(4E*))1/3，其中E*和R可以由下式計算得到：

(22a)

(22b)

式中νi，Ei和Ri分別為接觸物體i(i=1, 2)泊松比、彈性模量以及接觸端的半徑。接觸區法向剛度k的計算表達式為

(23)

3.1收斂性研究

首先研究了HAM模型的數值結果收斂性，獲得結構的合理離散密度。假設恢復系數e*=1，摩擦系數為μ=2/3，獲得的法向最大壓下量值的收斂曲線如圖4所示。結果表明，當n=44，接觸點最大法向壓下量已經趨于穩定值，故本文后面計算中在柔性梁長度不變的情況下，離散段數均取n=50進行數值計算。

圖4　不同離散密度下計算得到的法向最大壓下量Fig.4　The maximum normal penetration under different discretized segments

3.2結果驗證

(1)實際實驗對比

為驗證本文模型的計算精度，本節依照表1給出的2組實驗初始條件(摩擦系數為μ=0.075，恢復系數e*=1)采用本文方法進行數值計算，并將本文結果與實驗結果[18]進行對比。

表1　兩組實驗的初始條件

Tab.1　Initial conditions of two groups of experiment

實驗L/mV(0)/(m·s-1)ω(0)/(rad·s-1)1組0.1-1.8861.0322組0.2-1.710.208

由于實驗給出的初始條件滿足初始黏滯條件(17)，所以接觸點的初始運動狀態為黏滯壓縮。圖5為分離時刻桿端的接觸點切向速度混合模型解和實驗數據。通過比較表明，數值結果與實驗結果基本吻合。二者之間誤差產生的主要原因可能是HAM模型是一維梁且在平面內運動，而實驗則是一個真實三維梁以一個初始角度自由落體運動。且與靜力學實驗相比，接觸動力學實驗結果數據本身離散性也較大。此外，經計算除了當u0(0)=π/2時，接觸點運動模式為黏滯壓縮-滑動壓縮，其余角度值下接觸點的運動模式均為黏滯壓縮-滑動壓縮-滑動恢復。

圖5　接觸端接觸結束時切向速度Fig.5　Tangential velocity of the contact point at the end of contact

(2)數值實驗對比

由于真實的實驗環境能夠獲得的實驗數據有限，為進一步地驗證本文混合分析模型的精度，將本文模型的結果與非線性商業有限元的仿真結果進行對比。在商業有限元Ls-Dyna模擬中，本著疏密有致的離散原則，為了滿足非線性迭代接觸算法以及更好地捕捉接觸響應，對接觸區實施了網格細化(如圖6所示)，顯然這樣會增加自由度，降低計算效率。

圖6　梁與剛性平面彈性碰撞的有限元離散模型Fig.6　Finite element model for bar oblique impacts on rigid half-space

圖7是不同摩擦系數μ下的法向和切向接觸力時程曲線。計算結果表明，本文的方法與FEM方法獲得的接觸力響應所得的峰值、接觸時間以及曲線的變化形式吻合較好，接觸力峰值誤差僅為12%。二者之間誤差產生的原因主要是FEM模型是一個三維梁模型，由于泊松效應梁會產生橫向變形，同時桿件也容易產生橫向運動。從圖7～9可以看出，當μ=2/3時，系統的運動模式依次為滑動壓縮-滑動恢復-分離，切向沒有經歷黏滯運動；當摩擦系數變大，即當μ=5/3和μ=8/3時，系統的運動情況由滑動壓縮-黏滯壓縮-黏滯恢復-反向滑動恢復-分離構成，在接觸面切向經歷了滑動-黏滯-反向滑動的3個運動狀態的轉換。當μ=5/3和μ=8/3時(此摩擦系數恰好處于剛體斜碰撞模型中的經典painlevé悖論區[14])，法向接觸力響應和切向接觸力響應的峰值幾乎相等，即當摩擦系數超過5/3時，接觸力的峰值保持不變。這與完全剛體模型的摩擦顯著不同，造成這個現象的原因正是由于考慮了切向接觸柔度和存在切向黏滯運動。分析還發現純滑動時接觸力的峰值明顯小于具有黏滯狀態的接觸力的峰值。

圖7　切向和法向接觸力(-V1(0)/V3(0)=3.5，u0=3π/4，η2=1.21, e*=1)Fig.7　Normal and tangential components of contact force during oblique impact of a beam against a rough half-space (-V1(0)/V3(0)=3.5，u0=3π/4，η2=1.21, e*=1)

圖8和9分別為不同摩擦系數μ下接觸點的法向和切向相對速度響應。圖中的商業有限元解選擇接觸區中心的單元的速度。當摩擦系數相同時，二者的曲線基本一致，誤差較小(例如，接觸時間誤差僅為5.6%)，說明本文給出的混合分析模型是合理有效的。分析圖8發現，摩擦系數對分離時刻接觸點的法向速度有著重要影響。摩擦系數越大，接觸點的法向速度的最小值越小。分析還發現，無論摩擦系數多大，接觸點C的法向速度的響應曲線在接觸發生不久總存在震蕩現象。當μ=5/3和μ=8/3時，接觸點的法向速度在初始階段減小，即法向加速度在接觸開始的一段時間內出現負值，使得法向速度出現負向增加，導致系統的切向摩擦力變大，若采用完全剛性體模型進行分析，此時會出現經典的Painlevé悖論現象。分析圖9發現，當μ=2/3時，C點切向速度增加，始終大于初始速度。但當μ=5/3和μ=8/3時，C點切向速度隨時間的推移發生降低，然后反向滑動增加。

圖8　接觸點C和C′法向無量綱相對速度(-V1(0)/V3(0)=3.5, u0=3π/4,η2=1.21,e*=1)Fig.8　Normal relative non-dimensional velocity between C and C′ (-V1(0)/V3(0)=3.5, u0=3π/4, η2=1.21, e*=1)

圖9　接觸點C和C′的切向無量綱相對速度(-V1(0)/V3(0)=3.5,u0=3π/4,η2=1.21,e*=1)Fig.9　Tangential relative non-dimensional velocity between C and C′ (-V1(0)/V3(0)=3.5, u0=3π/4, η2=1.21, e*=1)

4　結　論

本文提出采用混合分析模型研究柔性梁摩擦斜碰撞問題，利用事件驅動法深入研究了整個接觸過程中的切向黏滑狀態和法向壓縮恢復階段的反復切換。將HAM模型的計算結果與實驗結果以及商業有限元方法進行了對比，研究得到了以下主要結論：

(1) 本文混合模型解與實驗數據吻合較好，驗證了該模型具有較好的收斂性和較高的精度。由于考慮了接觸區切向柔度，該模型具備了分析一系列切向黏-(反向)滑運動模式的能力。

(2) 研究表明：結構的整體柔度和接觸區的切向接觸柔度會對接觸力響應和接觸點的切向黏滑運動模式帶來較大的影響。

(3) 由于引入了能量恢復系數，因此混合模型可考慮接觸區彈塑性變形導致的碰撞能量的損失。

[1]Elkins-Tanton L T. Evolutionary dichotomy for rocky planets [J]. Nature, 2013, 497: 570—572.

[2]Popov V L. Contact Mechanics and Friction: Physical Principles and Applications [M]. Berlin University of Technology, 2010.

[3]Hayne P O, Greenhagen B T, Foote M C, et al. Diviner lunar radiometer observations of the LCROSS impact [J]. Science, 2010, 330: 477—479.

[4]Humphery J A C. Fundamentals of fluid motion in erosion by solid particle impact [J]. International Journal of Heat and Fluid Flow, 1990, 11(3): 170—195.

[5]Chapnik B V, Heppler G R, Aplevich J D. Modeling impact on a one-link flexible robotic arm [J]. IEEE Transaction on Robotic and Automation, 1991, 7(4): 551—561.

[6]張鎖懷，沈允文，董海軍，等. 齒輪拍擊系統的動力響應[J]. 振動工程學報, 2003, 16(1)：62—66.

ZHANG Suohuai, SHEN Yunwen, DONG Haijun, et al. Dynamic response of a gear rattling system [J]. Journal of Vibration Engineering, 2003, 16(1)：62—66.

[7]李曉杰，莫非，閆鴻浩，等. 爆炸焊接斜碰撞過程的數值模擬[J]. 高壓物理學報, 2011, 25(2)：173—176.

Li Xiaojie, Mo Fei, Yan Honghao, et al. Numerical simulation of the oblique collision in explosive welding [J]. Chinese Journal of High Pressure Physics, 2011, 25(2): 173—176.

[8]徐從啟，謝旭輝，戴一帆. 摩擦接觸約束下的微小管道機器人管內運動穩定性分析[J]. 機械工程學報, 2010, 46(15)：36—44.

Xu Congqi, Xie Xuhui, Dai Yifan. Motion stability analysis of micro in-pipe robot with frictional contacts [J]. Journal of Mechanical Engineering, 2010, 46(15): 36—44.

[9]段潔，丁千. 三質體斜碰撞振動的動力學和減振研究[J]. 振動工程學報, 2013, 26(1)：68—74.

DUAN Jie, DING Qian. Dynamics and vibration reduction of three-mass system with oblique-imapcts [J]. Journal of Vibration Engineering, 2003, 16(1)：62—66.

[10]Han W, Jin D P, Hu H Y. Dynamics of an oblique-impact vibrating system of two degrees of freedom[J]. Journal of Sound and Vibration, 2004, 275(3-5): 795—822.

[11]Han W, Hu H Y, Jin D P. Experimental study on dynamics of an oblique-impact vibrating system of two degrees of freedom[J]. Nonlinear Dynamics, 2007, 50(3): 551—573.

[12]Xie J H, Ding W C. Hopf-Hopf bifurcation and invariant torus T2of a vibro-impact system[J]. International Journal of Non-linear Mechanics, 2005, 40: 531—543.

[13]彼得.艾伯哈特，胡斌. 現代接觸動力學[M]. 南京:東南大學出版社, 2003.

Eberhart P, Hu B. Modern Contact Dynamics [M]. Nanjing: Southeast University Press, 2003.

[14]Shen Y, Stronge W J. Painlevé Paradox during oblique impact with friction[J]. European Journal of Mechanics A/Solids, 2011, 30(4): 457—467.

[15]Johnson K L. The bounce of ‘super ball’[J]. International Journal of Mechanical Engineering,1983,111: 57—63.

[16]Stronge W J. Impact Mechanics[M]. Cambridge University Press. 2000.

[17]Stoianovici D, Hurmuzlu Y. A critical study of the applicability of rigid body collision theory[J]. ASME Journal of Applied Mechanics, 1996, 63(2): 307—316.

[18]Hurmuzlu Y. An energy based coefficient of restitution for planar impacts of slender bars with massive[J]. ASME Journal of Applied Mechanics, 1996, 65(4): 952—962.

Research on rigid body-spring-particle hybrid model for flexible beam under oblique impact with friction

SHENYu-nian,GUJin-hong

(Department of Mechanics and Engineering Science, Nanjing University of Science and Technology， Nanjing 210094， China)

In order to quickly and accurately estimate the complex transient behavior of oblique collision with friction, a hybrid analysis model (HAM) is presented in this paper. This model takes account of the normal contact compliance, the tangential contact compliance of the local contact zone and the global compliance of the beam altogether. Based on the finite segment theory, the displacement field of the non-contact zone of the beam is discretized into a spring-damper-rigid body system. The bilinear process of normal compression-restitution deformation and the tangential deformation process of the local contact zone are described by an elastic-plastic spring-particle system with energy coefficient of restitution. The dynamic equations of the HAM which are piecewisely continuous under different contact states are derived. According to the switching criterion between compression and restitution as well as stick and slip, numerical examples are solved by the event-driven method. Through comparison of the results with the experimental data and those obtained by the commercial finite element software, it is shown that the hybrid model can accurately estimate the time history of contact forces and describe the repeatedly transition between stick and slip motion, which validate the applicablity of HAM.

nonlinear vibration; oblique impact; tangential compliance; frictional contact; stick-slip motion

2014-04-25；

2015-09-21

國家自然科學基金資助項目(11302107，11572157，11372138)；高等學校博士學科點專項科研基金資助項目(20123219120042)

O322; O313.4

A

1004-4523(2016)01-0001-07

10.16385/j.cnki.issn.1004-4523.2016.01.001

沈煜年(1982—),男,副教授。電話:(025)84315590; E-mail:shenyunian@aliyun.com