結合面微觀接觸面積分布數值仿真分析方法

蘇 蕾，荊建平，2

( 1.上海交通大學 機械系統與振動國家重點實驗室，上海200240；2.上海交通大學 燃氣輪機研究院，上海200240)

機械系統通常由各個零部件通過諸多結合面之間的配合連接構成。任何零件表面在微觀下都不是絕對光滑的，存在幾何形貌的偏差，實際接觸只在存在于較高的微凸體上，因此兩粗糙面構成的結合面實際接觸面積必然遠小于配合面的面積。結合面的靜、動態特性直接影響整機的剛度、阻尼等[1]，進而影響整機振動、噪聲等特性。因此要獲得更可靠的整機特性必須加深對結合面的微觀接觸特性的研究。

結合面接觸研究由來已久，早期研究將Hertz理論用于單個接觸點研究，假設粗糙面由等高但相互獨立的微凸體構成，所有微凸體產生相同的彈性變形[2]。由于粗糙表面微凸體的等高分布與表面實際輪廓不符，Greenwood 和Williamson[3]提出基于統計分析的接觸模型，假設粗糙峰高度隨機。根據粗糙面的分形特征，Majumdar 和Bhushan 提出以分形函數為基礎的接觸模型。Komvopoulos[4]提出修正的MB 接觸模型，并指出結合面真實接觸面積只有名義接觸面積的1%或更少。Jiang[5]根據分形理論建立了結合面法向接觸剛度模型，并與試驗數據對比。李小鵬[6]等建立了引入摩擦因素的結合面剛度阻尼模型。

MB 分形接觸模型是目前國內外研究人員分析結合面接觸狀態的常用模型。本文以MB分形接觸模型為基礎，在粗糙表面形貌的有限分形嵌套的實際情況下，依據W-M函數的最大空間頻率以連續可導曲線表征分形特征。提出接觸面積分布的數值求解方法，通過仿真直接求得分形函數描述的粗糙面與平面接觸的面積分布，并與早期研究的分布函數對比。

1 粗糙面微觀形貌描述

以往學者通過對粗糙面實測輪廓的研究發現粗糙表面輪廓曲線具有統計意義上的相似，即無規則分形特性。傳統歐氏幾何的各類粗糙度系數難以全面表征粗糙表面形貌的特征。隨機化的Weierstrass函數是不同頻域下的幾何累加，Mandelbrot 以此為基礎構造了具有分形特征的函數——W-M 函數[7]，其表達式如公式(1)所示。分形幾何為基礎的W-M函數利用不同波長和振幅的余弦波相互疊加以表征粗糙表面輪廓更完整的形貌特征，因此廣泛應用于粗糙面接觸問題。該函數具有處處連續但不可導、自仿射等特性，圖1為該函數模擬出的輪廓曲線

Z——表面輪廓高度；

x——表面輪廓位置坐標；

G——特征尺度系數；

D——分形維數；

γn——輪廓空間頻率，γ為大于1的常數，高斯分布的隨機輪廓一般γ取值為1.5；輪廓最低頻率γnl≈1/L，其中L為粗糙面樣本的長度。

圖1 W-M函數模擬輪廓曲線

2 MB分形接觸模型

機械系統中結合面的真實狀態是兩不光滑的表面相互接觸，具體表現為具有分形特征粗糙面與粗糙面的接觸，如圖2所示。下表面輪廓高于上表面輪廓的陰影部分為兩粗糙面的實際接觸變性部分。

圖2 粗糙面與粗糙面接觸模擬圖

Majumdar 和Bhushan[8]提出基于分形幾何的粗糙表面接觸模型，將無摩擦情況下兩粗糙表面的靜態接觸問題簡化為理想剛性平面與粗糙平面的接觸，并且忽略材料硬度隨表面深度的變化和相鄰微凸體變形產生的相互作用，如圖3所示。該模型以W-M分形函數和全球島嶼截面積分布規律為基礎，假定微凸體的曲率半徑與表面輪廓有關，使得模型具有尺度獨立性。

圖3 粗糙面與剛性平面接觸模擬圖

MB 模型將兩相互無摩擦的粗糙平面受載情況下的靜態接觸簡化為各個獨立的微凸體與剛性平面的接觸問題。根據hertz接觸理論，微凸體發生彈性變形，實際接觸面積at為變形前微凸體與平面截面積a的一半[9]。單個微凸體接觸模型如4所示。

Mandelbrot[10]研究全球島嶼形貌時發現，設島嶼與海平面的截面積值為A，所有島嶼中A大于一定值a的島嶼數量為N，則N與a存在冪律關系。圖5為島嶼示意圖。

圖4 微凸體接觸模型

圖5 島嶼示意圖

Majumdar[8]根據島嶼最大截面積為aL且個數為1，提出的面積分布函數如公式(2)所示

其中：D為島嶼分布規律的分形維數。對其進行微分，可得截面積分布函數如下

MB 分形模型以此表示分形表面微凸體的截面分布情況，由上述分布函數求期望可得分形表面與剛性平面總的截面積A為

MB 模型中假定微凸體截面積的分布函數僅在(0,aL]內非零，K.Komvopoulos[6]提出對原始微凸體接觸截面面積離散分布函數的非零域合理近似擴展。其假設只有一個最大接觸截面積為aL，將分布函數在a=aL鄰域內替換為迪利克雷函數，定義區間擴展因子ψ，進而提出了更加精確的分形表面微凸體截面面積分布修正函數，如公式(5)所示

并提出D與ψ的關系如公式(6)所示

對n(a) 積分求原函數可得數量N與面積a的關系如公式(7)所示。

公式(7)相比與公式(2)多了系數ψ，根據公式(6)可知ψ與分形維數D一一對應且取值范圍[1.718，2.618]，因此公式(7)所表達的接觸點面積分布仍只與分形維數D和接觸點最大面積aL相關。

3 接觸截面分布的數值計算方法

3.1 分形輪廓曲線的簡化

自然界的分形與數學上的分形幾何相比具有兩個特征。其一，自然界中的分形僅在一定尺度范圍、一定層次中才表現出分形特征；其二，數學上的分形可以具有無限層嵌套結構，而自然界的分形多是統計自相似的，僅有有限次的嵌套結構[2]。因此，W-M函數中三角函數項只需做有限次的疊加，即正整數n不能取到無窮大，存在最大值nmax.輪廓最大空間頻率與儀器分辨率有關γnmax≈1/Re，即儀器能達到的最大采樣頻率。

數學上的分形幾何中，因為存在無限次嵌套，所以W-M函數具有處處連續但不可導的特性，但對于表面形貌的具體問題，限制最大嵌套次數后，函數即具有連續可導的特性。

利用Majumdar 和Tien 在1990年[7]模擬粗糙表面的一組數據：掃描長度L為4 mm、儀器分辨率Re為5 μm、特征尺度系數G為12.445 nm、分形維數D為1.882。根據樣本長度和儀器分辨率計算得到nmin=14，nmax=30，將以上各參數帶入W-M函數的公式模擬粗糙表面。因為x=0 點時余弦函數值恒為1，與n無關，不具有一般性，為避開x=0點，所以此時x的范圍取5×10-5～4×10-3m。圖6為步長分別為5×10-6m和2×10-8m 的模擬輪廓曲線，即以儀器分辨率為采樣頻率和其他更高的采樣頻率分別模擬。

因為步長為5×10-6m時采樣點個數遠小于W-M函數的最高空間頻率γ30≈191 751，所以顯示的模擬輪廓線是首尾相連的折線。根據當前條件下WM函數連續可導的特性，將步距取為2×10-8m，從而使得采樣點個數大于最高頻率，可以獲得連續輪廓曲線。實測表面輪廓受儀器采樣頻率的限制，由圖6可看出，采樣點不足時難以表征函數完整特性。根據圖7給出的輪廓高度分布規律，采樣點大量增加并不會影響輪廓高度的分布規律(高斯分布)，在統計上仍具有相同的特性，且MB 模型假設各微凸體相互獨立。因此后續研究利用大量采樣點以獲得函數完整特征，即利用連續可導的曲線模擬粗糙面輪廓。

圖6 限定空間頻率的W-M模擬輪廓曲線

圖7 輪廓高度分布規律

3.2 微凸體截面積的求解

根據MB接觸模型的假設，建立如圖8的二維簡化接觸模型，利用3.1 中的參數模擬分形曲線，剛性平面與水平線距離為h。截線長度L即對應微凸體截面的理想直徑，由此可獲得理想截面積。

圖8 分形接觸輪廓模型

求得所有輪廓峰與截線的截斷長度，需先求得峰值高于h的極大值點。通過限定疊加次數的W-M函數，生成連續可導的光滑曲線，求出曲線的極大值即為輪廓的峰值，利用極值點與截線高度求得截斷長度。由于函數過于復雜，且需要求解所有極大值點，傳統求導獲得極值的方法難以實現，只能利用數值解法。

一般算法都是求解函數最值，諸如遺傳算法、模擬退火算法、梯度下降法等，避免獲得局部最值。本研究為了獲得所有峰值，采用將函數分段，在每段上求解最值，當區間劃分足夠精細時，每個小區間內最多只有一個極值，從而獲得所有極大值點。為降低計算量，先判斷區間是否存在極大值點然后再求解，根據區間前提條件，若存在極大值則左端點導數大于零、右端點導數小于零，如圖9所示。

圖9 分段區間求極大值

對限定截止頻率后的W-M 函數求導獲得梯度計算公式(8)

利用梯度上升法求極大值，本文構造爬坡算法，如公式(9)，其中a為步長，為加快求解速度需設置動態步長。

為判斷所求結果是否正確，驗證極值點的梯度并將區間進一步縮小再次求解，所得極值點個數基本無差別，極大值點均已獲得。

根據所獲得的極值點，若某點的值大于h，即與截線相交。構建求解算法，如公式(10)所示，從極值點出發，若前兩步對應F值的符號相同則β等于0，否則等于1。反復迭代求得極值點左右最近的高度等于h的點，左右兩點間的橫坐標差即近似為截線長度，遍歷所有符號條件極值而獲得所有截線長度。

相鄰極值點可能解得同一條截線，如圖10所示。

圖10 相鄰峰值點得同一段截線

計算過程中需要去除相鄰峰值點計算獲得重復的截線，否則不能反映真實分布情況。

3.3 面積分布規律分析

根據前述分析，求得的所有截線長度L，即為與平面相交的各微凸體的截面理想直徑。公式(2)與公式(7)中冪運算的底數均為aL/a，且易知aL/a=(Lmax/L)2，因此可根據所得截線長度進一步獲得利用W-M 函數模擬曲面與平面接觸狀態的截面分布規律。數據處理過程中去除奇異值，獲得(Lmax/L)2以及與之對應的長度大于L的個數N。取3.1中模擬分形輪廓的參數，h取1.5×10-7m時求得接觸點面積個數分布如圖11所示。

圖11 不同函數擬合所得曲線

利用公式(2)，以D作為待求參數擬合離散點，求得D的值約為1.746；利用公式(7)擬合，求得D值約為1.7；模擬輪廓所用分形維數1.882 均接近兩公式你擬合獲得的值，誤差分別為7.8%、10.7%。

表1所示為分形維數1.882 的粗糙面輪廓曲線與平面在不同名義距離下，求得接觸點面積分布對應擬合方程參數。當距離較遠時，所得參數D與分形維數更接近，且擬合曲線的均方根誤差（RMSE）較小，即可信度高；距離較近時，擬合所得結果誤差較大，且擬合曲線的RMSE也較大。

表1 平面與曲面不同距離求得分布函數

如圖8所示，大接觸點是由小接觸點的疊加獲得，平面距離較近時，發生彈性變形的部分相對更多。多個接觸點的合并成一個大接觸點，導致局部接觸點數量減少，部分接觸點面積變大，因而求解得D的值相對較小且擬合的均方根偏大。

如圖11所示，公式(2)的擬合效果并不好，本文提出對數形式的面積分布函數，如公式(11)所示。利用此公式擬合各點，所得曲線更貼近離散點的分布。

對不同D所描述的模擬粗糙面輪廓曲線，求得接觸面積分布離散點，進而求得對應的參數k，D與k的關系如圖12所示。

圖12 參數k與分形維數關系

參數k與分形維數D呈近似線性關系。當分形維數增大時，參數k增大，即接觸點的數目越多，符合粗糙面的真實接觸情況。

4 結語

本文以粗糙表面輪廓的分形特征為基礎，根據模擬輪廓的W-M函數的最大空間頻率，獲得了連續可導的模擬輪廓曲線，并對函數模擬的粗糙輪廓與平面接觸的面積分布開展數值計算及統計分析。利用數值方法求得輪廓極值分布，進而求得接觸點面積分布情況，驗證了MB 模型的分布函數近似符合W-M 函數描述輪廓的接觸面積分布。研究結果表明存在對數形式的分布函數更符合接觸面積分布，但其具體的參數仍需后續探討。該數值方法通過求解直線與曲線相交部分截線長度獲得分布規律，亦可求解兩曲線相交部分的截線，因而該方法有望用于求解兩粗糙面都具有分形特性的結合面接觸模型的接觸面積分布規律。