蘇 蕾，荊建平，2

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

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

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

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

1 粗糙面微觀形貌描述

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

Z——表面輪廓高度；

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

G——特征尺度系數(shù)；

D——分形維數(shù)；

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

圖1 W-M函數(shù)模擬輪廓曲線

2 MB分形接觸模型

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

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

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

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

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

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

圖4 微凸體接觸模型

圖5 島嶼示意圖

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

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

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

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

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

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

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

3 接觸截面分布的數(shù)值計算方法

3.1 分形輪廓曲線的簡化

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

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

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

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

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

圖7 輪廓高度分布規(guī)律

3.2 微凸體截面積的求解

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

圖8 分形接觸輪廓模型

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

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

圖9 分段區(qū)間求極大值

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

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

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

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

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

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

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

3.3 面積分布規(guī)律分析

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

圖11 不同函數(shù)擬合所得曲線

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

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

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

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

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

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

圖12 參數(shù)k與分形維數(shù)關系

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

4 結(jié)語

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