王曉軍，王琪

(1.常州工學院 機電學院，常州 213002;2.北京航空航天大學 航空科學與工程學院，北京 100191)

物體間的接觸、碰撞和摩擦現象是普遍存在的.如航天器的空中對接、飛機的著陸、步行機器人在地面上行走、機械手抓取工件以及具有非理想約束鉸鏈(含間隙與摩擦)的機械系統運轉等，都存在物體間的接觸與分離、接觸點的粘滯與滑移(stick-slip)等現象.如何描述物體間接觸或碰撞時的相互作用力，如何給出粘滯(stick)與滑移判別的計算方法，是建立系統動力學方程并通過數值仿真分析其動力學行為的重要問題之一.

20世紀末和21世紀初，人們利用多種接觸力模型研究物體的接觸與碰撞問題，如Kelvin-Voigt接觸力模型、Hertz接觸力模型、Hunt-Crossley接觸力模型和Lankarani-Nikravesh(L-N)接觸力模型等［1］.目前被廣泛使用的摩擦模型有Coulomb干摩擦模型、修正的Coulomb摩擦模型、粘性摩擦模型、Stribeck摩擦模型［2］和 Dahl摩擦模型［3］等，Coulomb摩擦模型被認為是較簡單的摩擦模型［1，4］.

文獻［5］利用修正的Kelvin模型研究了鄰近的兩個建筑物的碰撞問題.文獻［6］利用L-N接觸力模型研究了含間隙轉動鉸的速回機構的動力學建模與數值算法問題，并給出了接觸點檢測的數值方法，該文獻未涉及摩擦問題.文獻［7］利用L-N接觸力模型、線彈性接觸力模型和修正的Coulomb摩擦模型研究了含間隙滑移鉸機械系統的動力學建模與數值算法問題，當滑移鉸中的物塊頂角與滑道碰撞時用L-N接觸力模型，當物塊的一個邊(面)與滑道接觸時用線彈性接觸力模型.該文獻采用修正的Coulomb摩擦模型無法反映物塊與滑道間的stick-slip現象.文獻［8-10］分別研究了具有含摩擦轉動鉸和滑移鉸平面多剛體系統動力學的建模與數值算法問題，文獻中均采用Coulomb干摩擦模型，應用試算法或線性互補方法有效地解決了stick-slip運動狀態(tài)的切換問題，但無法反映鉸鏈處的碰撞問題，當鉸鏈的間隙不能忽略時，該方法則不適用.文獻［11］研究了非對稱庫侖干摩擦模型對振動驅動機械系統動力學行為的影響，豐富了非光滑系統動力學.

本文基于非光滑動力學理論，給出了一種物體接觸點含對稱或非對稱Coulomb干摩擦的平面運動剛體動力學的建模與數值計算方法.該方法是將研究對象視為剛體，考慮接觸點的局部變形，將物體間的法向接觸力表示成嵌入量與嵌入速度的函數;摩擦模型采用對稱或非對稱的Coulomb干摩擦模型;利用事件驅動法，將由于摩擦引起的stick-slip運動狀態(tài)切換的判斷問題以及處于粘滯(stick)狀態(tài)時靜摩擦力的計算問題轉化成線性互補問題的求解.該方法與傳統方法相比，不包含與接觸-分離相關的互補量，因此線性互補方程的維數低，計算效率高，將適用于對稱Coulomb干摩擦的互補算法推廣到非對稱Coulomb干摩擦模型.最后通過數值算例分析了非光滑平面運動剛體的動力學特性，并驗證了該方法的有效性.

1 接觸力學模型

1.1 法向接觸力模型

當物體間接觸或碰撞時，若考慮物體的局部變形，則一個物體會局部嵌入到另一個物體中，其法向嵌入量用δ表示，如圖1所示.用Hertz接觸模型，法向接觸力FN可表示［7］為

式中:K為廣義剛度系數;n為指數，其通常的取值范圍為［1，1.5］.當物體是點(局部)接觸時，該指數n=1.5，將接觸部位近似為半徑為R的圓弧，其廣義剛度系數可表示［7］為

式中:

其中:Ei和νi分別為兩接觸物體的彈性模量和泊松比.

圖1 兩物體間的接觸模型Fig.1 Contact model between two bodies

由于Hertz接觸力模型中不含阻尼項，這與實際情況有一定的差距.文獻［12］給出了一種非線性粘彈性接觸力模型，其表達式為

式中:c為廣義阻尼系數，該接觸力模型中的彈性項δK為δ的線性函數.文獻［13］使用另一種非線性粘彈性接觸力模型，其表達式為

該接觸力的彈性項和阻尼項均為非線性的，能較好地刻畫接觸力的特性.文獻［1］還給出了一種隱式的法向接觸力模型:

接觸力模型的研究還在不斷的深入和完善.

1.2 切向接觸力模型

機械系統中常用的摩擦模型有多種［2-3］，其中Coulomb摩擦模型又可分為Coulomb干摩擦模型和修正的Coulomb摩擦模型，前者是相對速度的非連續(xù)函數(數值計算存在困難)，后者是相對速度的連續(xù)函數(不易反映摩擦的靜動態(tài)特性)［14］.本文將采用非對稱Coulomb干摩擦模型，該模型的數學表達式［11］為

式中:

其中:Fτ為作用在物體上的摩擦力在接觸面切向上的投影;μ+、μ－和 μ0+、μ0－分別為物體間的正、負向動摩擦因數以及正、負向靜摩擦因數(當μ+=μ－，μ0+=μ0－時，即為對稱的 Coulomb 干摩擦模型分別為接觸點的相對速度和相對切向加速度.由式(7)～式(9)可知，當時，接觸點處于stick狀態(tài)，摩擦力的取值是一個范圍.

2 非光滑動力學方程及其數值求解算法

2.1 非光滑動力學方程

若設物體與固定面接觸時，其接觸面的法向量平行于y軸，用Newton-Euler方法可得平面運動剛體的動力學方程:

因為嵌入量和嵌入速度可以表示狀態(tài)變量的函數，因此，從式(5)不難看出，法向接觸力FNi是狀態(tài)變量)的函數;由式(7)和式(8)可知，當時，摩擦力Fτi也是狀態(tài)變量的函數.根據狀態(tài)變量可計算出FNi和Fτi.

2.2 Coulomb摩擦定律的互補關系

非光滑動力學方程數值計算的難點在于物體間接觸點stick-slip狀態(tài)切換的判斷以及處于stick狀態(tài)時摩擦力的計算.本文通過建立非對稱Coulomb摩擦定律的互補關系，將上述計算問題轉化為線性互補問題的求解.

可將該摩擦余量的概念推廣到非對稱Coulomb干摩擦模型，對應的正向摩擦余量和負向摩擦余量可分別定義為

由式(12)可得

由式(14)可得

2.3 動力學方程的求解算法

由式(13)、式(15)、式(17)和方程式(10)聯立，可以得到線性互補方程:

式中:wi1、wi2、Hxi和 Hθi均為狀態(tài)變量的函數.應用線性互補問題(LCP)的數值算法求解式(18)，可得，再由式(13)，可求出摩擦力 Fτi，將其代入動力學方程式(10)中，然后用常微分方程的數值方法求解該方程.

3 算例

3.1 算例模型

設圖2所示箱體為均質體，其質量為m，邊長分別是2a和2b，可與水平固定面發(fā)生接觸和摩擦，圖中C為箱體的質心，其坐標為(x，y)，用本文給出的方法對其進行動力學仿真.

圖2 算例中的箱體Fig.2 Slider in this example

該箱體頂角1和頂角2的嵌入量可分別表示為

箱體頂角1和頂角2與固定面接觸時的切向(沿x軸方向)速度可分別表示為

設系統參數為:m=2.0 kg，a=0.4m，b=0.3m;接觸力模型由式(5)給出，其中 n=1.5，K=2.96 ×107N/m1.5，c=1.89 ×106N·s/m2.

3.2 數值仿真分析

算例1 箱體在重力作用下自由下落并與地面發(fā)生接觸.設摩擦模型為對稱的Coulomb干摩擦模型，其摩擦因數為 μ+=μ－=0.3，μ0+=μ0－=0.4;初始條件為 x0=0，y0=0.8m，θ0=

圖3給出了箱體頂角1和頂角2的法向接觸力FN1和FN2的時間歷程圖;圖4給出了δ1和δ2的時間歷程圖.由此可以看出，箱體的頂角2先與地面發(fā)生碰撞，其碰撞速度為4.96m/s，然后該點發(fā)生第二次碰撞，其碰撞速度為0.65m/s;隨后頂角1與地面碰撞，其碰撞速度為3.16m/s;兩個頂角經過幾次碰撞后，箱體最終靜止在地面上.由于碰撞速度不同，對應的碰撞沖擊力的峰值也就不同.

圖5給出了x、y和θ的時間歷程圖，其運動狀態(tài)與箱體受力狀態(tài)吻合.圖6給出了箱體質心速度以及箱體所受摩擦力的時間歷程圖.由于箱體頂角2前兩次與地面發(fā)生碰撞時，有向左的相對運動，因此摩擦力的方向向右，在其作用下，物體的質心有向右的速度.當箱體與地面再次發(fā)生碰撞時，由于運動速度向右，因此摩擦力向左.經過幾次碰撞后，最終物體相對地面靜止.

圖3 F N1和F N2的時間歷程圖Fig.3 Time history diagram of F N1 and F N2

圖4 δ1和δ2的時間歷程圖Fig.4 Time history diagram ofδ1 and δ2

圖5 x、y和θ的時間歷程圖Fig.5 Time history diagram of x、y and θ

圖6 和Fτ的時間歷程圖Fig.6 Time history diagram ofand Fτ

算例2 箱體放在水平地面上，其上作用一個水平力 Fx=A sin t.設摩擦因數為 μ+=0.3，μ－=0.4，μ0+=0.5，μ0－=0.6，初始條件為

算例1和算例2均用試算法進行了驗證，試算法與本文給出的算法所得數值結果完全吻合.

圖7 當A=11，12，18N時x和的時間歷程圖Fig.7 Time history diagram of x andwith A=11，12，18N

4 結論

基于非光滑動力學方法，利用線性互補理論，給出了一種含非對稱Coulomb干摩擦的平面運動剛體動力學的數值算法.

1)將物體間的法向接觸力表示成物體間嵌入量和嵌入速度的函數，無需引入互補量，更易于分析物體間的接觸與分離.

2)建立了非對稱摩擦余量互補關系，將原適用于含對稱Coulomb干摩擦的LCP算法推廣到非對稱Coulomb干摩擦系統，使其具有更廣泛的適用性.

3)數值仿真算例表明，該算法易于分析非光滑系統中物體間的接觸與分離和stick-slip現象.

References)

［1］ Flores P，Ambrósio J.Kinematics and dynamics of multibody systems with imperfect joints［M］.Berlin:Springer，2008:47-66.

［2］劉麗蘭，劉宏昭，吳子英，等.機械系統中摩擦模型的研究進展［J］.力學進展，2008，38(2):200-213.Liu L L，Liu H Z，Wu Z Y，et al.An overview of friction models in mechanical systems［J］.Advances in Mechanics，2008，38(2):200-213(in Chinese).

［3］ Pennestri E，Valentini P P，Vita L.Multibody dynamics simulation of planar linkages with Dahl friction［J］.Multibody System Dynamics，2007，17(4):321-347.

［4］ Pfeiffer F，Glocker C.Multibody dynamics with unilateral contacts［M］.New York:John Wiley ＆ Sons，Inc，1996:70-91.

［5］ Ye K，Li L，Zhu H P.A modified Kelvin impact model for pounding simulation of base-isolated building with adjacent structures［J］.Earthquake Engineering and Engineering Vibration，2009，8(3):433-446.

［6］ Flores P，Ambrósio J.On the contact detection for contact-impact analysis in multibody systems［J］.Multibody System Dynamics，2010，24(1):103-122.

［7］ Flores P，Ambrosio J.Translational joints with clearance in rigid multibody systems［J］.ASME Journal of Computational and Nonlinear Dynamics，2008，3(1):011007.

［8］莊方方，王琪.含摩擦柱鉸鏈平面多體系統動力學的建模和數值方法［J］.工程力學，2012，29(5):193-199.Zhuang FF，Wang Q.Modeling and numerical algorithm for planar multibody system with friction on revolute joints［J］.Engineering Mechanics，2012，29(5):193-199(in Chinses).

［9］ Wang Q，Peng H L，Zhuang F F.A constraint-stabilized method for multibody dynamics with friction-affected translational joints based on HLCP［J］.Discrete and Continuous Dynamical Systems Series B，2011，16(2):589-605.

［10］ Zhuang F F，Wang Q.Modeling and simulation of the nonsmooth planar rigid multibody systems with frictional translational joints［J］.Multibody System Dynamics，2013，29(4):403-423.

［11］ Fang H B，Xu J.Dynamics of a three-module vibration-driven system with non-symmetric Coulomb’s dry friction［J］.Multibody System Dynamics，2012，27(4):455-485.

［12］ Kisu L.A short note for numerical analysis of dynamic contact considering impact and a very stiff spring-damper constraint on the contact point［J］.Multibody System Dynamics，2011，26(4):425-439.

［13］劉才山，陳濱，王玉玲.考慮摩擦作用的多柔體系統點-面碰撞模型［J］.中國機械工程，2000，11(6):616-619.Liu C S，Chen B，Wang Y L.The point-surface impact model considering the friction effect in flexible multi-body system［J］.China Mechanical Engineering，2000，11(6):616-619(in Chinses).

［14］王琪，莊方方，郭易圓.非光滑多體系統動力學數值算法的研究進展［J］.力學進展，2013，43(1):101-111.Wang Q，Zhuang F F，Guo Y Y.Advances in the research on the numerical method fir non-smooth dynamics of multibody systems［J］.Advances in Mechanics，2013，43(1):101-111(in Chinese).

［15］ Pfeiffer F，Foerg M，Uibrich H.Numerical aspects of nonsmooth multibody dynamics［J］.Computer Methods in Applied Mechanics and Engineering，2006，195(50):6891-6908.