唐進元 廖東日 周 煒
中南大學(xué)高性能復(fù)雜制造國家重點實驗室,長沙,410083
很多數(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)生所需的非高斯序列。
自相關(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ù)的表達式如下:
按高度分布函數(shù)的不同,粗糙表面可分為高斯型粗糙表面與非高斯型粗糙表面,兩種粗糙表面生成的理論模型均為時間序列模型中的滑動平均模型(MA)。
基于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時的高斯分布表面
實際工程中,工程表面大多為非高斯分布表面,表面上的峰相對谷來說比較容易被除去,所以很多粗糙表面呈負(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 模擬表面
粗糙表面模擬與理論模型要求的誤差可從高度分布函數(shù)與自相關(guān)函數(shù)的對比上得出。其中高度分布函數(shù)統(tǒng)計主要進行均值、方差、偏態(tài)和峰態(tài)四個參數(shù)的對比,而相關(guān)函數(shù)則采用實測表面相關(guān)函數(shù)與數(shù)值模擬表面相關(guān)函數(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)表面概率密度分布
相關(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ù)
表面高度的均值、標(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.