粗糙界面法向接觸振動響應與能量耗散特性研究

肖會芳， 邵毅敏， 徐金梧

(1.北京科技大學 國家板帶生產先進裝備工程技術研究中心，北京 100083；2.重慶大學 機械傳動國家重點實驗室，重慶 400044)

動力系統(tǒng)中，結構之間通過接觸界面進行動力傳遞。界面不僅極大地影響機械裝備振動測量信號的特征[1]，而且是影響機械零件性能，諸如動態(tài)性能、接觸疲勞強度、摩擦功耗和磨損壽命等至關重要的因素[2]。

目前，關于接觸界面振動響應與能量耗散的研究幾乎都集中在沿接觸界面切向的滑動摩擦耗散[3-5]。盡管關于無潤滑法向接觸振動能量耗散的實驗研究表明，并未觀察到明顯的法向微動能量耗散，但卻尚未從理論上予以解釋[6]。經典的Hertz球體接觸理論被廣泛地用于描述接觸界面法向的力-變形關系、接觸面積以及壓力分布特征[7]。實際中，通過機械加工處理的工程結構表面，并非絕對平滑，而是具有不同程度的粗糙度[8-9]。對粗糙接觸界面而言，其法向接觸具有怎樣的力-變形關系，接觸振動具有怎樣的動力學響應特性，能量耗散具有怎樣的特征與機理，與Hertz接觸模型相比有何異同，對描述界面動力學機理具有重要的理論和實際意義。

本文建立粗糙界面法向接觸振動的動力學模型，提出了不同形貌粗糙界面法向接觸振動的動力學響應特征量和振動能量耗散量的計算方法，并構造了粗糙接觸界面的接觸力-變形關系；建立了接觸振動系統(tǒng)的動力學方程，計算了不同形貌粗糙界面法向接觸振動每周期的振動能量耗散率和累積振動能量耗散率；分析了粗糙界面法向接觸振動的響應特征與能量耗散特性；從理論上解釋了文獻[6]的實驗結果。

1 粗糙界面法向接觸振動模型

1.1 動力學模型描述

研究粗糙界面法向接觸振動與能量耗散特性的模型示意圖，如圖1所示。金屬板與剛性平面的接觸面為具有一定表面形貌的三維粗糙表面。在法向載荷作用下，金屬板在固定的剛性平面上進行接觸振動。由于存在界面阻尼，其振動過程伴隨著能量耗散。

圖1 粗糙界面法向接觸振動模型示意圖

圖1所示模型的動力學特性可以采用單自由度模型進行描述，如圖2所示。該單自由度模型被廣泛地用于描述Hertz“球-剛性平面”接觸振動[10-12]。其假設條件是小球質量中心的位移量與局部接觸區(qū)域的位移量相同，即小球的變形僅發(fā)生在局部接觸區(qū)域[7]。其中，彈簧的回復力為Hertz彈性接觸公式的力-變形關系式，阻尼模型為線性粘彈性阻尼。

對圖1所示的粗糙界面接觸模型，其接觸變形幾乎完全發(fā)生在接觸區(qū)的粗糙體，質量中心的位移量與局部接觸區(qū)域的位移量一致[13]?；谠撎卣?，其接觸振動模型也可采用圖2所示的單自由度模型進行描述。其中，彈簧的回復力為粗糙表面與剛性平面之間的彈性接觸力-變形關系式。對無潤滑界面，阻尼特性采用線性粘性阻尼模型表征[10]。

盡管對粗糙界面，由于實際接觸面僅為一些點和很小的面，因而容易在接觸點發(fā)生應力集中，使應力值高于屈服應力而產生塑性變形[13]，純彈性接觸幾乎不存在。但對界面的法向接觸振動而言，通常的情況是接觸體在靜平衡位置附近振動。此時，塑性變形產生于接觸體達到靜平衡位置的過程中，而其微振動則仍屬于彈性接觸范圍[14]。

圖2 接觸振動的單自由度模型

1.2 動力學方程

圖2所示的單自由度模型在靜平衡位置附近的自由振動方程可以表示為:

(1)

其中:m為金屬板質量，c為阻尼系數(shù)，fk(z,zs)為彈性回復力，即為粗糙界面的彈性接觸力-變形關系式，zs為金屬板在重力作用下的靜變形量。對于式(1)，為了保持金屬板與剛性平面接觸，必須符合z≥-zs。

采用如下的無量綱位移u，無量綱時間τ，無量綱阻尼比ζ和線性接觸頻率參數(shù)ωs，對式(1)進行無量綱化

其中，k(zs)為靜平衡位置處的接觸剛度。則式(1)對應的無量綱表達式為:

(2)

(3)

其中，無量綱的激勵幅值和激勵頻率分別為:

系統(tǒng)的無量綱初始條件為u(τ0)=u0,u′(τ0)=v0。對式(2)，金屬板與剛性平面接觸的條件變?yōu)閡≥-1。

2 接觸振動響應特征量與能量耗散量計算方法

為了獲得圖1所示接觸振動模型的動力學響應特征量與能量耗散量，首先需要確定式(1)的彈性力-變形關系表達式fk(z,zs)。

2.1 粗糙界面接觸力-變形關系式計算方法

基于粗糙表面的三維分形模型描述[13,15]，通過建立三維粗糙表面模型，采用有限元方法，對圖1所示粗糙界面模型進行彈性接觸計算分析，獲得不同表面形貌接觸界面的接觸力-變形關系式[13,16]。

三維分形表面采用修正的兩參數(shù)Weierstrass-Mandelbrot函數(shù)描述，其表達式為[15]

(4)

其中，L為樣本長度，Ls為截斷長度，γ(γ>1)為縮放參數(shù)，M為生成分形表面的脊線數(shù)，n是頻率因子，nmax=int[log(L/Ls)/logγ]，Φm,n是[0,2π]內的隨機相位。表面粗糙程度由尺度獨立的分形參數(shù)D和G控制。表面分形維數(shù)D，其物理意義是粗糙表面所占據(jù)的空間程度大小，D值越大對應于越密集的表面形態(tài)(更光滑的表面形貌)。表面分形粗糙度G是高度尺度參數(shù)，G值越大對應越粗糙的表面形貌。參數(shù)D和G一般通過實驗測試確定。根據(jù)實驗結果，D的范圍為[2.3,2.7]，G的范圍為[1.36E-13,1.36E-10]m[15]。

采用式(4)生成的具有不同表面形貌的0.9 μm×0.9 μm三維分形表面，如圖3所示。其中，Ls=1.5E-7m,M=10，γ=1.5[15]，D=2.4。圖3顯示，對相同的分形維數(shù)D，增大表面分形粗糙度G的表面形貌更粗糙。

圖3 不同表面形貌的三維分形粗糙表面

2.1.1 有限元計算模型

有限元計算模型示例，如圖4所示。三維金屬板用SOLID185單元離散，SOLID185單元的每個節(jié)點具有3個方向的平動自由度。金屬板與剛性平面之間建立接觸對，接觸對通過點-面接觸單元對TARGE170 和CONTA175識別。接觸單元CONTA175覆蓋在金屬板的接觸表面，用來描述變形體的邊界條件，并與目標單元TARGE170進行接觸。剛性平面全約束，金屬板的上表面節(jié)點僅具有Z方向的自由度。在金屬板的上表面節(jié)點施加均布的位移載荷，求取剛性平面的反力，獲得金屬板與剛性平面之間的接觸力。計算過程中，位移載荷通過多個載荷子步逐漸施加，最大載荷步和最小載荷步分別設置為400和50。接觸算法采用Augmented Lagrangian算法。力的收斂準則設為0.001。

圖4 有限元計算模型(D=2.3,G=1.36E-11m)

2.1.2 有限元計算方法有效性的驗證

采用限元方法模擬Hertz接觸與Hertz理論計算結果對比的方法，檢驗本文所采用的有限元計算方法處理接觸問題的有效性。計算模型為半徑為R=8 mm的球在法向力F=5 kN的作用下與剛性平面進行無摩擦接觸。材料的屬性參數(shù)為彈性模量E=30 GPa，泊松比ν=0.25。建立Hertz接觸的有限元模型的方法、離散球體的單元類型、識別界面接觸的接觸單元類型以及接觸算法與建立粗糙界面接觸的有限元模型保持一致。通過有限元計算求解接觸半寬b和接觸變形量d，并與Hertz接觸理論計算結果對比[7]，如表1所示。從表中可以看出，有限元計算結果與Hertz理論計算結果的誤差小于3%，表明本文所采用的有限元計算方法是準確、可靠的。

表1 有限元和Hertz理論計算結果對比

2.2 振動響應特征量與能量耗散量計算方法

通過計算獲得的接觸力-變形關系，求解動力學方程式(1)，可以獲得其振動響應特征量和能量耗散量。粗糙界面法向接觸自由振動的振動響應特征通過瞬時頻率和等效阻尼比描述。能量耗散特性則通過每周期的能量耗散率和一段時間內的累積能量耗散率描述。

2.2.1 瞬時頻率

由于接觸振動系統(tǒng)的非線性特性，系統(tǒng)自由振動每周期內的頻率并不相同，而隨幅值變化。假設每個振動周期內的響應可以表示為:

u=Ucos(Ωτ+φ)

(5)

其中，U,Ω和φ分別為瞬時幅值、瞬時頻率和瞬時相位，在單個周期內保持恒定，且:

(6)

瞬時頻率可以表示為

(7)

其中，Δτ為相鄰峰之間的時間間隔。

2.2.2 等效阻尼比

對自由振動系統(tǒng)，描述其響應特征的等效阻尼比可以采用對數(shù)衰減率進行估算，其表達式為:

(8)

其中，δ為對數(shù)衰減率，Ui為第i個振動周期的幅值(i=1,2,3,…)。

2.2.3 能量耗散率

自由振動每周期的能量耗散量可表示為

(9)

式中:τ1和τ2為相鄰峰值對應的時間，fd為阻尼力，其表達式為:

fd(u′)=2ζu′

(10)

將式(5)，式(10)代入式(9)，并定義φ=Ωτ+φ，每周期的能量耗散量可以表示為:

(11)

式(11)顯示，系統(tǒng)每周期的能量耗散量由瞬時頻率、瞬時振幅、瞬時相位和阻尼比決定。

進一步定義系統(tǒng)每周期的能量耗散率為每周期的能量耗散量與輸入能量的比值

(12)

其中，Ein為系統(tǒng)的輸入能量，在初始條件為u(τ0)=u0，u′(τ0)=v0時，其表達式為:

(13)

系統(tǒng)在從自由振動開始時刻τ0到任意時刻τ內的累積能量耗散率為:

(14)

3 數(shù)值計算結果分析

3.1 粗糙界面彈性接觸力-變形關系表達式

界面無摩擦，表面形貌參數(shù)不同，最大法向位移δ=0.1時的力-變形曲線，如圖5所示。其中，圖5(a)為粗糙表面的表面分形粗糙度G=1.36E-11m，而表面分形尺度參數(shù)D不同，分別為D=2.3,2.4,2.5,2.6；圖5(b)為表面分形尺度參數(shù)D=2.4，而表面分形粗糙度G不同，分別為G=1.36E-10m,1.36E-11m,1.36E-12m,1.36E-13m。

圖5顯示，界面的接觸力隨變形非線性遞增，接觸剛度呈非線性遞增(曲線的斜率逐漸增大)；不同表面粗糙度界面的非線性有所差異：粗糙度較大的界面(較小的D值或較大的G值)具有更強的非線性(曲線斜率的變化更劇烈)。相同位移載荷下，粗糙度較大界面的接觸力小于粗糙度較小的界面；隨著接觸表面的粗糙度逐漸變小，界面的接觸力-變形關系逐漸接近于絕對平滑界面。

圖5 界面無摩擦，不同形貌界面的接觸力-變形曲線

進一步考慮界面摩擦對接觸力-變形關系的影響。界面存在摩擦的計算方法與前文相同，不同的是界面處增加了摩擦因子μ=0.2。此時，不同表面形貌粗糙界面的接觸力-變形曲線，如圖6所示。與界面無摩擦曲線相比，界面存在摩擦時：力-位移曲線形狀一致，即界面的接觸剛度變化趨勢一致；在相同位移下，接觸力數(shù)值增大。

不同表面形貌界面的接觸力與變形關系可用指數(shù)函數(shù)式(15)表示：

fn(δ)=kδα

(15)

圖5和圖6的擬合結果表明：式(15)能很好地描述粗糙界面的接觸力-變形關系，表達式中的各系數(shù)數(shù)值，如表2所示。為驗證力-變形關系式(15)的有效性，對比了不同表面粗糙度值G的摩擦界面(圖6(b))力-變形關系式與數(shù)值計算結果間的相對誤差，如圖7所示。最大相對誤差為3.4%，最小相對誤差為0.05%，圖7所示結果表明式(15)可準確地表征實際計算結果。

上述計算結果基于尺寸為0.9 μm×0.9 μm的分形表面模型。由于分形表面的自相似特性，表面形貌模型與尺度無關[15]。因而，式(15)同樣適用于其他尺寸。

圖6 界面存在摩擦，不同形貌界面的力-變形曲線

3.2 粗糙界面接觸與Hertz接觸力-變形關系式對比分析

表2顯示：受表面粗糙度和界面摩擦的影響，式(15)所示非線性力-變形關系指數(shù)α的范圍為α=[1.168,2.195]；且隨著表面粗糙度減小(D值較大或G值較小，接觸界面更光滑)，系數(shù)k和指數(shù)α逐漸減小。

表2所示的最大和最小α值對應的力-變形曲線，以及Hertz接觸理論的力-變形關系曲線，如圖8所示。圖8顯示，對不同的α值，曲線所示的非線性度以及回復力的大小存在較大的差異。因而，與Hertz接觸理論相比，對粗糙接觸界面，其力-變形關系仍然具有指數(shù)函數(shù)關系式；不同的是，受表面粗糙度的影響，指數(shù)函數(shù)的指數(shù)范圍為α>1，而非α=3/2。

3.3 接觸振動響應與能量耗散特性

基于獲得的粗糙界面接觸力-變形關系式(15)，式(1)所示的自由振動方程可以表示為:

(16)

式(2)所示的無量綱表達式為

(17)

3.3.1 瞬時頻率

求解式(17)所示的非線性振動方程。其中，不同表面粗糙度對應的接觸剛度模型以α=1,3/2,2和5/2為例進行描述；沖擊載荷的幅值為σ=0.6，頻率為Ω0=1；計算采用的阻尼系數(shù)為ζ=0.005[10-11]。

采用式(7)計算不同接觸剛度模型自由振動響應對應的瞬時頻率，如圖9所示。圖9顯示，對不同的剛度模型(表面粗糙度)，其振動瞬時頻率存在較大差異：對剛度指數(shù)為α=1的線性系統(tǒng)，瞬時頻率保持Ω=1不變；對剛度指數(shù)為α=3/2，α=2和α=5/2的非線性系統(tǒng)，瞬時頻率隨著幅值的增大而減小，其最大值為Ω=1，對應于靜平衡位置U=0；隨著非線性剛度指數(shù)α增大，瞬時頻率隨幅值的變化更劇烈，系統(tǒng)具有更強的非線性。同時，圖9顯示，相對于線性系統(tǒng)，不同剛度指數(shù)對應的非線性系統(tǒng)的瞬時頻率隨幅值的減小量很小，表明該激勵載荷下的系統(tǒng)僅表現(xiàn)出較弱的非線性；隨著激勵載荷增大，其非線性會增強。

圖8 具有不同指數(shù)大小的非線性指數(shù)函數(shù)力-變形關系曲線fn(δ)= δα

系統(tǒng)振動的瞬時頻率隨著振幅的增大而增大，表明系統(tǒng)具有硬彈簧特征；反之，系統(tǒng)則具有軟彈簧特征[17]。因此，剛度指數(shù)為α>1的系統(tǒng)均對應于軟彈簧非線性系統(tǒng)。表2所列的計算結果顯示，對粗糙接觸界面，其力-變形關系的剛度指數(shù)均為α>1，因而，粗糙界面接觸振動的響應特性具有軟彈簧非線性系統(tǒng)的動力學特征。

3.3.2 對數(shù)衰減特性

采用式(8)計算不同剛度模型自由振動曲線的等效阻尼比隨振幅UN的變化曲線，如圖10所示。其中，UN=Ui/U1(i=1,2,3,…)。圖10顯示，僅剛度指數(shù)為α=1的線性系統(tǒng)的阻尼比保持計算值ζ=0.005不變，其他剛度指數(shù)對應的等效阻尼比均隨著振幅的減小而遞減。同時可以看出，對線性系統(tǒng)，可以采用對數(shù)衰減公式估算其阻尼比[18]；但是對非線性剛度系統(tǒng)而言，即使系統(tǒng)具有線性的阻尼特性，采用對數(shù)衰減公式估算獲得的系統(tǒng)阻尼比仍然與振幅相關。因而，對瞬時頻率隨振幅變化的非線性系統(tǒng)，對數(shù)衰減公式并不適用于估算其阻尼比。

3.3.3 能量耗散特性

基于獲得的接觸力-變形關系式(15)，式(13)所示的系統(tǒng)輸入能量可以表示為:

(18)

采用式(12)，式(14)和式(18)計算不同剛度模型自由振動響應每周期的能量耗散率和累積能量耗散率，如圖11所示。圖11顯示，對不同指數(shù)的接觸剛度模型(表面形貌)，每周期的能量耗散率均隨著時間逐漸減小為0，而累積能量耗散率則逐漸增大到1；當無量綱時間τ≈600時，ND=0，NSD=1，即τ≈600時的累積能量耗散量達到系統(tǒng)的輸入能量。圖11同時顯示，對不同指數(shù)的接觸模型，其能量耗散率隨時間變化基本完全相同，該結果表明界面法向振動的能量耗散量很小，其值基本不隨剛度模型(表面形貌)的變化而變化。該計算結果與實驗研究觀測獲得的并無明顯法向微動能量耗散的結論一致，從理論上驗證了文獻[6]的實驗結果。

圖11 不同接觸剛度模型對應的能量耗散率隨時間變化曲線

4 結 論

本文建立了粗糙界面法向接觸振動模型，提出了計算其法向接觸振動的動力學響應特征量和能量耗散量的計算方法。針對具有不同表面形貌的粗糙界面，基于自相似和尺度獨立的粗糙表面分形模型，對接觸振動過程中，界面接觸力-變形關系、接觸振動系統(tǒng)的動力學響應特征以及振動能量耗散特性進行了研究，并與Hertz接觸模型進行了對比，主要結論如下：

(1)本文的計算結果從理論上解釋了界面法向微動能量耗散的實驗結果，并進一步驗證了提出的界面接觸振動模型和計算方法是有效的。

(2)不同表面形貌粗糙界面的彈性接觸力-變形關系可以用指數(shù)函數(shù)關系表示。指數(shù)函數(shù)指數(shù)α的范圍為α>1，與Hertz接觸理論的α=3/2有所差異；指數(shù)大小由表面形貌決定，粗糙度較大的界面具有較大的α值。

(3)不同表面形貌粗糙界面法向微動的動力學響應特征量(瞬時頻率，對數(shù)衰減)存在較大差異，但是能量耗散特性卻基本相同。

參 考 文 獻

[1]Smith J D.Transmission of Smith shocks through rolling bearings [J].Journal of Sound and Vibration,1995,181(1):1-6.

[2]溫詩鑄,黃 平,等.界面科學與技術[M].北京: 清華大學出版社,2011.

[3]Berger E J.Friction modelling for dynamic system simulation [J].ASME Applied Mechanics Review,2002,55(6): 535-577.

[4]Gaul L,Nitsche R.The role of friction in mechanical joints [J].ASME Applied Mechanics Review,2001,54(2): 93-105.

[5]Ferri A A.Friction damping and isolation systems [J].ASME Journal of Vibration and Acoustics,1995,117B: 196-206.

[6]Andrew C,Cockburnt J A,Waring A E.Metal surfaces in contact under normal forces: some dynamic stiffness and damping characteristics [C].Proceedings of the Institution of Mechanical Engineers,Conference Proceedings,1967,182: 92-100.

[7]Johnson K L.Contact mechanics [M].Cambridge University Press,1985.

[8]Nayak P R.Some aspects of surface roughness measurement [J].Wear,1973,26:165-174.

[9]Al-Kindi G,Shirinzadeh B.An evaluation of surface roughness parameters measurement using vision-based data [J].International Journal of Machine Tools and Manufacture,2007,47: 697-708.

[10]Sabot J,Krempf P,Janolin C.Non-linear vibrations of a sphere-plane contact excited by a normal load [J].Journal of Sound and Vibration,1998,214: 359-375.

[11]Rigaud E,Perret-Liaudet J.Experiments and numerical results on non-linear vibrations of an impacting Hertzian contact.Part 1: harmonic excitation [J].Journal of Sound and Vibration,2003,265: 289-307.

[12]Xiao H F,Brennan M J,Shao Y M.On the undamped free vibration of a mass interacting with a Hertzian contact stiffness [J].Mechanics Research Communications,2011,38: 560-564.

[13]邵毅敏,肖會芳.動力系統(tǒng)非連續(xù)單一疊加結構界面變形與能量損耗特性[J].振動與沖擊,2011,30(4):217-222.

SHAO Yi-min,XIAO Hui-fang.Deformation and energy loss of single layered and discontinuous structural interface of a power system [J].Journal of Shock and Vibration,2011,30(4):217-222.

[14]Sherif H A.Parameters affecting contact stiffness of nominally flat surfaces [J].Wear,1991,145: 113-121.

[15]Yan W,Komvopoulos K.Contact analysis of elastic-plastic fractal surfaces [J].Journal of Applied Physics,1998,84: 3617-3624.

[16]肖會芳,邵毅敏,周曉君.非連續(xù)粗糙多界面接觸變形和能量損耗特性研究[J].振動與沖擊,2012,31(6): 83-89.

XIAO Hui-fang,SHAO Yi-min,ZHOU Xiao-jun.Contact deformation and energy dissipation characteristics of a discontinuous rough multi-interface [J].Journal of Shock and Vibration,2012,31(6): 83-89.

[17]Nayfeh A H,Mook D T.Nonlinear Oscillations [M].Wiley,New York,1979.

[18]Rayleigh J.The theory of sound,Vol.I [M].New York,reprinted by Dover,1945.