基于NCGM的粗糙表面數(shù)值模擬與實驗對比

唐進元 廖東日 周 煒

中南大學(xué)高性能復(fù)雜制造國家重點實驗室，長沙，410083

0 引言

很多數(shù)值方法已用于粗糙表面的模擬。Patir［1］采用線性矩陣變換方法模擬具有指定相關(guān)函數(shù)的表面，計算自相關(guān)矩陣時采用牛頓迭代法求解相關(guān)系數(shù)矩陣，該方法的主要缺點是所需的儲存空間很大，計算相當(dāng)費時，而且其收斂性也一直備受爭議。文獻［2－3］采用時間序列模型來模擬粗糙表面，采用的是自回歸滑動平均模型（ARMA）。Whitehouse［4］采用自回歸模型（AR）來模擬粗糙表面，但是ARMA與MA模型只考慮了低階系統(tǒng)，而且只描述了原點附近的自相關(guān)函數(shù)。針對儲存空間和計算速度的限制，快速傅里葉變換（FFT）方法也被用于粗糙表面的模擬。文獻［5－7］采 用 FFT 方 法 模 擬 粗 糙 表 面。Wu［6－7］指出，采用FFT方法模擬粗糙表面不能同時滿足正確的自相關(guān)函數(shù)與循環(huán)自相關(guān)函數(shù)，但是可以保證平均自相關(guān)函數(shù)的正確性，F(xiàn)FT方法只能應(yīng)用于自相關(guān)長度較小的粗糙表面的模擬。基于上述文獻的工作與Partir方法，為解決儲存空間與計算時間的限制問題，本文采用非線性共軛梯度法［8］（NCGM）計算自相關(guān)矩陣，模擬非高斯型粗糙表面時采用Johnson轉(zhuǎn)換系統(tǒng)［9－11］產(chǎn)生所需的非高斯序列。

1 粗糙表面統(tǒng)計規(guī)律

自相關(guān)函數(shù)可用于描述不同位置之間數(shù)據(jù)的相關(guān)程度，該函數(shù)提供了數(shù)據(jù)最基本的空間信息和數(shù)據(jù)相互之間的關(guān)系。因此，所有空間上的表面粗糙度參數(shù)都是從自相關(guān)函數(shù)獲取的。

三維表面的自相關(guān)函數(shù)形式如下：

其中，E是數(shù)學(xué)期望值，τx、τy分別是x 方向、y方向的相關(guān)距離，R（0，0）即為σ2（σ為均方根）。

其離散表達式為

其中，n、m分別為自相關(guān)長度在x、y方向上的最大值；N、M分別為離散點沿x、y方向的總數(shù)。

相關(guān)長度τ＊定義為在輪廓方向上自相關(guān)函數(shù)衰減到臨界值時的長度，Stout等［12］定義自相關(guān)函數(shù)從原點衰減到其值的10%時的長度為相關(guān)長度，此時相關(guān)函數(shù)的表達式如下：

2 粗糙表面數(shù)值模擬

按高度分布函數(shù)的不同，粗糙表面可分為高斯型粗糙表面與非高斯型粗糙表面，兩種粗糙表面生成的理論模型均為時間序列模型中的滑動平均模型（MA）。

2.1 高斯型粗糙表面模擬

基于Patir［1］的理論，生成（M，N）的粗糙表面需先產(chǎn)生（M＋m，N＋n）的獨立分布的隨機數(shù)和（m，n）維的自相關(guān)矩陣，生成粗糙表面的線性變換方程如下：

其中，ak，l是欲生成自相關(guān)函數(shù)的自相關(guān)系數(shù)矩陣第k行第l列元素，ηi，j是獨立的單位方差隨機數(shù)，滿足以下關(guān)系式：

自相關(guān)函數(shù)的定義如下：

式（6）是求解維數(shù)為m×n的自相關(guān)系數(shù)矩陣的元素ak，l的非線性方程組，Patir［1］采用牛頓法求解該方程組，給定迭代的近似初值為

當(dāng)用牛頓法求解式（6）時，Bakolas［13］發(fā)現(xiàn)，當(dāng)矩陣維數(shù)一定時，該方程組不收斂。除收斂性問題外，求解式（6）的計算量相當(dāng)大，還需要很大的空間來存儲雅可比矩陣。

本文采用NCGM來求解該非線性方程組，把式（6）改寫成：

該方法把求解原函數(shù)替換成求解原函數(shù)的一階導(dǎo)數(shù)。按照NCGM，現(xiàn)有兩種搜索方向選取方法，一種為Fletcher－Reeves方法，另一種為Polak-Ribiere方法，其表達式分別為

當(dāng)初值選取為非常接近目標(biāo)函數(shù)的解時，F(xiàn)letcher－Reeves方法會收斂，而Polak－Ribiere方法在極少數(shù)情況下會不收斂，但是Polak－Ribiere方法的收斂速度比前者快很多。在迭代過程中，為避免求函數(shù)導(dǎo)數(shù)，采用割線法來進行非精確線搜索，其中割線法的步長需不斷地調(diào)整來達到收斂，本文采用Polak－Ribiere方法求解方程組，求解步驟如下：

（2）d（0）＝r（0）－f′（x（0））（d為搜索方向）。

（3）x（i＋1）＝ x（i）＋α（i）d（i），找 到 令 f（x（i）＋α（i）d（i））最小的α（i）。

（4）r（i＋1）＝－f′（x（i＋1））。

（6）d（i＋1）＝r（i＋1）＋β（i＋1）d（i）。

指定自相關(guān)函數(shù)的粗糙表面模擬步驟如下：

（1）利用計算機產(chǎn)生一高斯分布的隨機序列μi，j。

（3）根據(jù)式（10），通過NCGM 求解得出自相關(guān)系數(shù)矩陣。

（4）根據(jù)式（4）計算得出模擬的粗糙表面z（x，y）。

基于以上理論和方法，圖1和圖2分別給出了相關(guān)長度βx＝βy＝3μm和βx＝βy＝10μm時的粗糙表面，粗糙表面的均值均為0，標(biāo)準(zhǔn)差為1，模擬粗糙表面點數(shù)M＝N＝60，x與y方向采樣間隔均為1μm?？煽闯霎?dāng)相關(guān)長度變大時表面形貌更加緩和。當(dāng)βx和βy相等時，為各向同性表面。

圖1 相關(guān)長度βx＝βy＝3μm時的高斯分布表面

當(dāng)βx和βy不相等時為各向異性表面，圖3給出了βx＝30μm，βy＝3μm時的各向異性粗糙表面。當(dāng)相關(guān)長度不同且數(shù)值上差別較大時，可明顯看出紋理的走向。

圖2 相關(guān)長度βx＝βy＝10μm時的高斯分布表面

圖3 相關(guān)長度βx＝30μm、βy＝3μm時的高斯分布表面

2.2 非高斯型粗糙表面模擬

實際工程中，工程表面大多為非高斯分布表面，表面上的峰相對谷來說比較容易被除去，所以很多粗糙表面呈負(fù)偏態(tài)，故相對而言負(fù)偏態(tài)表面的接觸性能較好。因此，為了更加真實地模擬粗糙表面，非高斯型表面的模擬尤為重要。粗糙表面分布的偏態(tài)為

粗糙表面的峰態(tài)為

非高斯型表面的模擬相對高斯型表面來說，其輸入隨機序列μi，j與高斯型表面模擬時不同，模擬高斯型表面時，μi，j是高斯分布隨機序列；而模擬非高斯型表面時，μi，j為非高斯分布隨機序列，所以必須產(chǎn)生一個μi，j來滿足需模擬表面的偏態(tài)與峰態(tài)。式（4）本質(zhì)上是一個MA模型，對于一個純MA模型，有

當(dāng)偏態(tài)與峰態(tài)通過一個數(shù)字濾波器時，其數(shù)值會改變，所以必須對輸入的偏態(tài)與峰態(tài)進行修正，其輸入與輸出存在以下對應(yīng)關(guān)系：

其中Sη、Kη分別為輸入序列的偏態(tài)與峰態(tài)。Sz、Kz分別為輸出序列的偏態(tài)與峰態(tài)，即為所求參數(shù)。根據(jù)式（4），得

把高斯分布的隨機序列轉(zhuǎn)換為指定偏態(tài)與峰態(tài)非高斯隨機序列需使用Johnson轉(zhuǎn)換系統(tǒng)，按照偏態(tài)與峰態(tài)取值范圍的不同，Johnson轉(zhuǎn)換系統(tǒng)可分為SU、SL、SB三種不同的形式，在不同的偏態(tài)與峰態(tài)下有以下公式對應(yīng)：

其中，SU為無界系統(tǒng)，SL為對數(shù)正態(tài)系統(tǒng)，SB為有界系統(tǒng)。x是標(biāo)準(zhǔn)正態(tài)分布的隨機序列，η是所得到的序列。γ、δ、ξ、λ是系統(tǒng)得到所需偏態(tài)與峰態(tài)的4個參數(shù)。

非高斯型粗糙表面按以下步驟模擬：

（1）求解非線性方程組（式（10））。

（2）產(chǎn)生高斯分布序列。

（3）通過式（16）和式（17）得出指定生成的偏態(tài)與峰態(tài)序列。

（4）利用Johnson轉(zhuǎn)換系統(tǒng)和對應(yīng)系統(tǒng)的公式產(chǎn)生非高斯序列。

（5）根據(jù)式（4）產(chǎn)生非高斯型粗糙表面。

模擬非高斯型表面時，為了更加接近真實表面的形貌分布，采用實際加工的表面進行模擬，故自相關(guān)矩陣由真實表面計算得出。該工件采用超聲磨削，磨削進給量為0.01mm，超聲振動振幅為4μm，頻率為20kHz，進給速度為69mm／s，砂輪轉(zhuǎn)速為1440r／min，砂輪型號為 WA40L，工件材料為45鋼，低溫回火處理。采用Veeco Wyko NT9100光學(xué)輪廓測量儀進行表面微觀形貌測量，采樣間隔為0.5μm，采樣點數(shù)為640×480。為減少計算量，截取其測量表面點數(shù)為60×60進行分析重構(gòu)，經(jīng)測量計算，x與y方向采樣間隔均為4μm。

按非高斯表面模擬方法，實驗掃描完整表面，實際測量區(qū)域表面與模擬的表面如圖4～圖6所示?？煽闯龀暷ハ骷庸け砻嬗忻黠@的加工輪廓。經(jīng)計算，相關(guān)長度在y方向非常大，而x方向上的相關(guān)長度很小，完全符合加工方向的紋理分布。

圖4 實驗測量表面

圖5 實際測量區(qū)域表面

圖6 模擬表面

3 粗糙表面模擬誤差分析

粗糙表面模擬與理論模型要求的誤差可從高度分布函數(shù)與自相關(guān)函數(shù)的對比上得出。其中高度分布函數(shù)統(tǒng)計主要進行均值、方差、偏態(tài)和峰態(tài)四個參數(shù)的對比，而相關(guān)函數(shù)則采用實測表面相關(guān)函數(shù)與數(shù)值模擬表面相關(guān)函數(shù)進行對比。

3.1 表面高度分布統(tǒng)計參數(shù)誤差分析

實際加工表面的算術(shù)平均偏差Ra為0，均方根偏差Rq為6.12μm，偏態(tài)Rsk為－0.39μm，峰態(tài)Rku為3.30μm。實際加工表面與重構(gòu)表面統(tǒng)計參數(shù)與概率密度分布曲線對比如表1、圖7和圖8所示?？煽闯鰞煞N表面在統(tǒng)計參數(shù)上可以很好地符合理論模型的要求。

表1 實際加工表面與重構(gòu)表面統(tǒng)計參數(shù)對比 μm

圖7 實際加工表面概率密度分布

圖8 重構(gòu)表面概率密度分布

3.2 表面相關(guān)函數(shù)誤差分析

相關(guān)函數(shù)反映的是表面各點之間的相關(guān)程度，相關(guān)函數(shù)和高度分布統(tǒng)計參數(shù)能完整地表征表面形貌。由于y方向相關(guān)長度很大，是主要紋理方向，故只需對比y方向自相關(guān)函數(shù)即可。實際加工表面與重構(gòu)表面的自相關(guān)函數(shù)對比如圖9所示，可以看出重構(gòu)的表面與實際加工表面的自相關(guān)函數(shù)能較好地吻合。

圖9 實際加工表面與重構(gòu)表面的相關(guān)函數(shù)

4 結(jié)語

表面高度的均值、標(biāo)準(zhǔn)差、偏態(tài)、峰態(tài)是描述高度分布函數(shù)的統(tǒng)計參數(shù)，而表面的相關(guān)性則用自相關(guān)函數(shù)來描述。本文采用MA模型和NCGM，能很好地克服牛頓法所帶來的儲存數(shù)據(jù)大、無法收斂的困難，可模擬出高度概率分布函數(shù)和自相關(guān)函數(shù)指定要求的粗糙表面。綜合使用NCGM和Johnson轉(zhuǎn)換系統(tǒng)能有效地模擬出工程表面，從而建立高精度的粗糙表面模型。

［1］Patir N.A Numerical Procedure for Random Generation of Rough Surfaces［J］.Wear，1978，47（2）：263－277.

［2］Watson W，King T，Spedding T，et al.The Machined Surface－time Series Modelling［J］.Wear，1979，57（1）：195－205.

［3］Watson W，Spedding T.The Time Series Modelling of Non－Gaussian Engineering Processes［J］.Wear，1982，83（2）：215－231.

［4］Whitehouse D.The Generation of Two Dimensional Random Surfaces Having a Specified Function［J］.CIRP Annals－Manufacturing Technology，1983，32（1）：495－498.

［5］Hu Y，Tonder K.Simulation of 3－D Random Rough Surface by 2－D Digital Filter and Fourier Analysis［J］.International Journal of Machine Tools and Manufacture，1992，32（1／2）：83－90.

［6］Wu J J.Simulation of Non－Gaussian Surfaces with FFT［J］.Tribology international，2004，37（4）：339－346.

［7］Wu J J.Simulation of Rough Surfaces with FFT［J］.Tribology International，2000，33（1）：47－58.

［8］Manesh K，Ramamoorthy B，Singaperumal M.Numerical Generation of Anisotropic 3DNon－Gaussian Engineering Surfaces with Specified 3DSurface Roughness Parameters［J］.Wear，2010，268（11／12）：1371－1379.

［9］Hill I，Hill R，Holder R.Algorithm AS 99：Fitting Johnson Curves by Moments［J］.Journal of the Royal Statistical Society Series C（Applied Statistics），1976，25（2）：180－189.

［10］Hill I.Algorithm AS 100：Normal－Johnson and Johnson－normal Transformations［J］.Journal of the Royal Statistical Society Series C（Applied Statistics），1976，25（2）：190－192.

［11］Johnson N L.Systems of Frequency Curves Generated by Methods of Translation［J］.Biometrika，1949，36（1／2）：149－176.

［12］Stout K J，Sullivan P J，Dong W P，et al.The Development of Methods for the Characterization of Roughness on Three Dimensions：EUR，0015178 EN［P］.1994－01－01.

［13］Bakolas V.Numerical Generation of Arbitrarily Oriented Non－Gaussian Three－dimensional Rough Surfaces［J］.Wear，2003，254（5／6）：546－554.