一種高計(jì)算收斂性的粗糙曲面三維形貌模型

喬澤龍,韓炬,董威

(華北理工大學(xué)機(jī)械工程學(xué)院,063210,河北唐山)

所有的工程表面都具有一定的粗糙性,適宜的粗糙面模型應(yīng)用在相關(guān)研究中能有效降低試驗(yàn)成本并提高分析效率,對(duì)于深入理解結(jié)合表面之間的接觸[1]、磨損[2]、潤(rùn)滑[3]、熱傳導(dǎo)[4]等機(jī)理具有重要的意義。國(guó)內(nèi)外學(xué)者研究粗糙面模型的構(gòu)建已有五十幾年[5],主要涉及兩個(gè)關(guān)鍵問題,其一是微凸體高度及分布情況能客觀反映實(shí)際粗糙面形貌特征,其二是微凸體接觸模型能有效反映實(shí)際接觸變形規(guī)律。本文重點(diǎn)關(guān)注粗糙曲面三維形貌特征模型的構(gòu)建。

機(jī)械加工面的粗糙形貌與加工工藝密切相關(guān),近年來,眾多學(xué)者結(jié)合刀具尺寸與傾角、加工路徑以及振動(dòng)等因素,實(shí)現(xiàn)了對(duì)機(jī)械加工工件表面形貌的預(yù)測(cè)[6-7],但一般用于工件表面質(zhì)量評(píng)測(cè),沒能應(yīng)用于機(jī)械結(jié)合面的接觸、潤(rùn)滑、磨損等特性分析。

現(xiàn)有獲取粗糙面的三維形貌主要包括試驗(yàn)法與數(shù)值法兩種方式。目前,表面形貌的測(cè)試方法很多[8],但受到橫縱坐標(biāo)分辨率的限制,測(cè)試獲取的形貌模型在很多細(xì)節(jié)上被簡(jiǎn)化,此外測(cè)試得到的模型只是對(duì)現(xiàn)有加工面的再現(xiàn),很難對(duì)獲取改性的曲面形貌進(jìn)行對(duì)比研究。

通過分析真實(shí)粗糙表面數(shù)據(jù),發(fā)現(xiàn)粗糙表面具有統(tǒng)計(jì)自仿射的特性[9-10],進(jìn)一步研究證實(shí)粗糙表面的分形特性與尺度無關(guān)[11]。雖然應(yīng)用數(shù)值方法構(gòu)建三維粗糙曲面的方法很多,但由于分形粗糙面在進(jìn)行接觸、潤(rùn)滑等分析時(shí)有較高的適用性[12-16]。周超等利用W-M分形函數(shù)合成粗糙表面在速度、精度以及合成表面的均方根偏差等方面均有一定優(yōu)勢(shì)[17]。因此,本文基于分形理論表征粗糙表面微凸體高度及分布特性,并在粗糙面模型中綜合體現(xiàn)曲面的宏觀輪廓與微觀形貌。黃健萌等在對(duì)納觀粗糙表面的建模研究中發(fā)現(xiàn)了分形表面不連續(xù)的問題,并選擇適合通帶的高斯濾波器處理數(shù)據(jù),獲得了較為理想的仿真粗糙面,濾波處理后的表面數(shù)據(jù)在進(jìn)行數(shù)值分析時(shí),計(jì)算收斂性有顯著提升[18]。

本文綜合應(yīng)用微分幾何理論與分形理論,構(gòu)建能表征復(fù)雜形狀曲面的宏觀輪廓與微觀形貌的三維模型,通過對(duì)模型進(jìn)行平滑處理,提升了模型在諸如接觸、疲勞、潤(rùn)滑等數(shù)值分析時(shí)的收斂性。

1 機(jī)加工粗糙曲面表面形貌模型

1.1 二維粗糙分形形貌表征的W-M函數(shù)

一般應(yīng)用W-M函數(shù)為[19]對(duì)粗糙面微觀形貌進(jìn)行分形模擬。W-M分形函數(shù)具有連續(xù)、處處不可導(dǎo)和自仿射的特點(diǎn),表達(dá)式為

(1)

式中:G為特征尺度,表征分形粗糙峰谷高度;γ為粗糙曲面的空間頻率,一般取γ=1.5;n為空間頻率指數(shù);D為分形維數(shù),二維分形模擬時(shí)D值取(1

W-M函數(shù)用于表達(dá)二維工程粗糙表面時(shí),一般取式(1)的實(shí)部,Majumdar和Bhushan根據(jù)式(1)進(jìn)一步提出了M-B函數(shù)[20],表達(dá)式為

(2)

1.2 形貌各向異性和各向同性表征函數(shù)

M-B函數(shù)可以表征二維分形曲線,但機(jī)械結(jié)合面接觸、摩擦、疲勞、潤(rùn)滑等研究的重點(diǎn)集中于面接觸問題,且工程中粗糙表面呈各向異性或各向同性,因此需要構(gòu)建粗糙表面的三維形貌模型。

通過二維W-M函數(shù)表達(dá)式,把x與y方向的W-M分形疊加在一起,實(shí)現(xiàn)x、y兩向異性的三維分形粗糙面,公式如下

(3)

式中:z(x,y)為的三維粗糙表面輪廓幅值;Gx、Gy表示x、y方向的二維尺度系數(shù);Dx、Dy表示x、y方向的二維分形維數(shù);γ為空間頻率,取值1.5;nmin為最低頻率階數(shù);nmax為最高頻率階數(shù)。若分形函數(shù)曲面尺寸為L(zhǎng)×L,取樣長(zhǎng)度L內(nèi)應(yīng)至少包含最低頻率成分一個(gè)波長(zhǎng),則最低空間頻率γnmin>(1/L),則有nmin=lg(1/L)/lg(γ);若取樣間隔為L(zhǎng)s,根據(jù)奈奎斯特采樣定理,最高空間頻率為γnmax<(1/2Ls),則有nmax=lg(1/2Ls)/lg(γ)。

式(3)中,x和y方向分形函數(shù)各空間頻率信號(hào)分量初始相位默認(rèn)為0。實(shí)際應(yīng)用中,若各頻率成分初相位為0,疊加后會(huì)導(dǎo)致曲線起始段出現(xiàn)翹起的邊緣效應(yīng)。為避免邊緣效應(yīng),引入隨機(jī)相位φnx和φny,取值范圍[0,2π],其含義為對(duì)應(yīng)空間頻率γn信號(hào)分量的隨機(jī)相位,公式為

(4)

分別基于隨機(jī)相位和0相位生成分形曲線,為了作圖清晰,空間頻率階數(shù)n取1～12。如圖1所示,圖1a、1b中分形曲線各頻率成分初相位為0,曲線起始值隨n增加而增加,曲線具有邊緣效應(yīng),曲線前端呈現(xiàn)出翹起狀態(tài);圖1c、1d中分形曲線各頻率成分具有不同的初相位時(shí),曲線邊緣效應(yīng)消除。

(a)0相位時(shí)的空間頻率分量幅值

表示各向同性三維分形粗糙面的W-M函數(shù)[17]如下

(5)

式中:h(x,y)為粗糙面微凸體高度;D為分形維數(shù)(2

2 三維分形形貌粗糙表面仿真

2.1 各向異性分形粗糙表面仿真

根據(jù)式(4),可得各向異性三維粗糙仿真曲面,如圖2所示,其中x和y方向特征尺度系數(shù)和分形維數(shù)相同,特征尺度系數(shù)取Gx=Gy=1.36×10-5mm,分形維數(shù)分別取Dx=1.488,Dy=1.518。采樣長(zhǎng)度為2 mm,采樣間隔為1 μm,則nmax=33,nmin=16。圖3為各向異性粗糙面在x方向和y方向的截面曲線。

圖2 各向異性分形粗糙面Fig.2 Anisotropic fractal rough surface

圖3 各向異性粗糙面x和y方向截面曲線Fig.3 Cross-section curves of anisotropic rough surface in x- and y-direction

2.2 各向同性分形粗糙表面仿真

根據(jù)式(5),可得到圖4所示的各向同性三維粗糙仿真曲面,取樣參數(shù)與各向異性相同,G=1.36×10-5mm,D=2.8。圖5為各向同性粗糙面在x方向和y方向截面分形粗糙輪廓曲線。

圖4 各向同性分形粗糙面Fig.4 Isotropic fractal rough surface

圖5 各向同性粗糙面x和y方向截面曲線Fig.5 Cross-section curves of isotropic rough surface in x- and y-direction

2.3 分形粗糙表面的濾波處理

分形粗糙表面具有處處連續(xù)但不可微的特點(diǎn),其輪廓放大時(shí)在更小尺度上仍然呈現(xiàn)出自仿射的粗糙特征,其上任一點(diǎn)的斜率都不存在,曲面形態(tài)表現(xiàn)為尖銳毛刺。分形函數(shù)獲得的粗糙表面直接用于接觸等數(shù)值分析計(jì)算時(shí)運(yùn)算收斂困難,有時(shí)在迭代計(jì)算初期出現(xiàn)發(fā)散現(xiàn)象,需要進(jìn)行平滑濾波,才能進(jìn)行后續(xù)的分析計(jì)算。若將不連續(xù)特征視為形貌疊加噪聲的結(jié)果,采用低通濾波器對(duì)分形粗糙面數(shù)據(jù)進(jìn)行空間濾波去噪處理,可平滑形貌毛刺,從而獲得微觀平滑的形貌數(shù)據(jù)。分形函數(shù)高頻成分會(huì)產(chǎn)生在微觀層次游離出表面的孤立物質(zhì)(原子團(tuán)),采用高斯濾波平滑可以消除由分形特性引起的不合理數(shù)據(jù)[18],但仍能保留模擬數(shù)據(jù)的分形特征。

圖6 高斯濾波后各向同性粗糙面Fig.6 Isotropic rough surface after Gaussian filtering

圖7 高斯濾波后各向同性x和y方向截面曲線Fig.7 Cross-section curves of isotropic rough surface in x- and y-direction after Gaussian filtering

采用高斯濾波對(duì)分形粗糙形貌高度進(jìn)行平滑時(shí),中心位置點(diǎn)高度為包括本點(diǎn)在內(nèi)的鄰域各位置高度數(shù)據(jù)按高斯分布規(guī)律的加權(quán)和,在實(shí)現(xiàn)對(duì)數(shù)據(jù)平滑的同時(shí),能夠保留形貌高度的原始總體分布特征。二維高斯濾波器的核心為高斯二維卷積算子,本文選取高斯核為9×9,σ=1.8。

將本文生成的分形粗糙數(shù)據(jù),應(yīng)用于二維點(diǎn)接觸彈流潤(rùn)滑分析,對(duì)濾波前后的仿真計(jì)算進(jìn)行對(duì)比。仿真節(jié)點(diǎn)數(shù)N=130,載荷w=105N,接觸球體半徑為0.02 m,最大赫茲接觸壓力為0.9 GPa,初始黏度為0.05 Pa·s。為了節(jié)約計(jì)算時(shí)間,將收斂相對(duì)精度設(shè)為兩次迭代,各節(jié)點(diǎn)油膜壓力之和相差小于1.6×10-4N。仿真計(jì)算結(jié)果如圖8和圖9所示,通過濾波前后壓力兩次迭代油膜壓力差的變化情況可以看出,濾波后迭代過程加快將近一倍,原始數(shù)據(jù)迭代136次耗時(shí)11 475 s,濾波后數(shù)據(jù)迭代69次耗時(shí)6 419 s。歸一化油膜壓力分布如圖9a所示(X、Y為x、y方向歸一化長(zhǎng)度),由于原始數(shù)據(jù)中存在高頻細(xì)節(jié),油膜存在大量尖銳壓力峰值。工程實(shí)際中,工件表面在初期剛體接觸時(shí)尖銳微凸體會(huì)首先發(fā)生塑性變形[23-24],因此在彈流潤(rùn)滑工況下,圖9b所示的濾波后油膜分布更符合實(shí)際情況。

圖8 高斯濾波前后粗糙面點(diǎn)接觸彈流潤(rùn)滑仿真Fig.8 Elastohydrodynamic lubrication simulation of rough surface point contact before and after Gaussian filtering

(a)原始數(shù)據(jù)仿真

3 基于W-M函數(shù)與矢量函數(shù)的形貌合成模型

雖然機(jī)械動(dòng)態(tài)結(jié)合面的結(jié)合區(qū)很小,接觸瞬時(shí)的曲面可以近似為粗糙平面,這是眾多粗糙面模型將接觸區(qū)簡(jiǎn)化為平面接觸的重要原因。但結(jié)合面宏觀輪廓,如齒輪齒廓、凸輪輪廓等,在傳動(dòng)過程中對(duì)傳動(dòng)性能有關(guān)鍵影響,因此三維粗糙曲面的形貌模型應(yīng)該綜合體現(xiàn)曲面的宏觀輪廓與微觀形貌,即將體現(xiàn)微觀形貌的W-M函數(shù)與體現(xiàn)宏觀形貌的矢量函數(shù)疊加。

本節(jié)以RV減速器中擺線輪輪齒齒廓為例,創(chuàng)建三維粗糙曲面模型。RV減速器擺線針輪嚙合接觸屬于柱面接觸,擺線針輪接觸仿真計(jì)算需要基于擺線輪曲面構(gòu)建粗糙曲面,由光滑曲面與非均勻的微凸體疊加構(gòu)成粗糙曲面形貌模型,其中光滑曲面由矢量函數(shù)表示,微凸體的高度值由W-M函數(shù)計(jì)算,微凸體高度方向?yàn)橄鄳?yīng)曲面點(diǎn)的法矢量方向。擺線輪齒的粗糙表面可表示為擺線輪齒齒面的矢量函數(shù)和粗糙面矢量函數(shù)之和。

二維擺線輪廓函數(shù)矢量平面曲線公式為[25]

r(θ)=x(θ)i+y(θ)j,r(θ)∈C0

(6)

式中:i和j為坐標(biāo)軸向的單位向量;x(θ)和y(θ)為擺線輪廓函數(shù)如下

(7)

(8)

(9)

式中:K1=e/Rg;γ0=jθ,θ為滾圓滾角;e為基圓和滾圓中心距。

用曲線弧長(zhǎng)代替式(2)中的x為自變量,可得

(10)

式中:h(s)為分形粗糙形貌輪廓高度;s為沿齒廓曲線的弧長(zhǎng),其他參數(shù)同式(2)。擺線輪齒廓弧長(zhǎng)為

(11)

粗糙齒廓曲線的二維矢量方程為擺線輪輪齒齒廓矢量方程和微凸體高度矢量之和,公式為

r*=r(θ)+h(s)m(θ)

(12)

式中:r(θ)為擺線齒廓曲線矢量;m(θ)為曲線點(diǎn)法向單位矢量,m(θ)表示式為

(13)

利用式(12),在Matlab中可生成如圖10所示的二維粗糙輪齒齒廓。

圖10 粗糙齒廓矢量曲線Fig.10 Rough surfaces and vector end curves

3.1 擺線針輪齒廓二維形貌建模

為便于數(shù)值計(jì)算,基圓和擺線矢量方程可變形為歐拉公式形式,基圓矢量方程為

(14)

擺線矢量方程為

(15)

修形后的擺線輪齒的齒廓矢量方程可轉(zhuǎn)化為

(16)

利用式(10)～(16)可實(shí)現(xiàn)擺線輪齒的齒廓建模,如圖11所示。擺線輪半徑為52 mm,針齒半徑為2 mm,偏心距為0.9 mm,擺線輪齒數(shù)為39,等距修形量、移距修形量分別為-0.01 mm和-0.012 mm。擺線輪齒分形粗糙齒廓建模如圖12所示。

圖11 矢量法生成擺線輪齒的齒廓 Fig.11 Profile of cycloid gear tooth generated with vector method

圖12 矢量函數(shù)合成擺線輪齒粗糙齒廓曲線 Fig.12 Rough profile curve of cycloid gear tooth synthesized by vector function

3.2 擺線輪輪齒齒廓的三維形貌建模

擺線輪齒齒廓三維建模的核心問題是如何將理論光滑齒輪曲面模型與粗糙面模型相結(jié)合。粗糙面建模可依據(jù)式(4)和式(5)獲得模型數(shù)據(jù),并采用二維高斯濾波對(duì)數(shù)據(jù)進(jìn)行平滑。理論齒廓三維建?？刹捎脠A柱坐標(biāo)矢量函數(shù),獲得齒廓三維曲面,再通過數(shù)值計(jì)算獲得曲面上點(diǎn)的單位法向矢量。各曲面法向矢量以粗糙曲面相應(yīng)點(diǎn)高度值為模長(zhǎng),獲得粗糙曲面矢量矩陣,理論齒面矢量與粗糙曲面矢量之和即為粗糙齒廓的矢量模型。

擺線輪三維齒廓在圓柱坐標(biāo)系中的矢量方程為

(17)

式中:z為齒廓軸向坐標(biāo),其他參數(shù)與式(16)相同。得到的擺線理論齒廓三維模型如圖13所示。

圖13 擺線理論齒廓三維模型Fig.13 3D model of cycloid tooth profile

在齒面法向單位矢量計(jì)算中,由于曲面z方向的導(dǎo)數(shù)為0,故法向矢量表達(dá)式為

(18)

求得曲面法向矢量如圖14所示。

圖14 齒廓曲面和法向矢量Fig.14 Tooth profile surface and normal vector

擺線輪齒粗糙齒廓矢量方程為

Rc(θ,z)=R(θ,z)+h(s(θ),z)mc(θ,z)

(19)

式中:h(s(θ),z)mc(θ,z)為各向異性或各向同性粗糙曲面矢量。依據(jù)式(19),按經(jīng)驗(yàn)取特征尺度系數(shù)和分形維數(shù),采用Matlab模擬擺線輪三維粗糙齒面、各向異性和各向同性的擺線輪齒粗糙曲面分別見圖15～圖17。

圖15 擺線輪三維粗糙齒面Fig.15 3D rough tooth surface of cycloid wheel

圖16 各向同性三維粗糙面局部放大圖Fig.16 Local zoom-in view of 3D isotropic rough tooth surface

圖17 各向異性三維粗糙面局部放大圖Fig.17 Zoom-in view of 3D anisotropic rough tooth surface

應(yīng)用PS50型Nanovea三維非接觸式表面形貌儀獲得真實(shí)擺線輪齒面圖像,如圖18所示。圖18a為珩磨機(jī)械加工面的微觀形貌,不同方向上的微凸體分布相對(duì)平均;圖18b為研磨機(jī)械加工面的微觀形貌,不同方向上微凸體的分布差別較大,有明顯的紋向。通過圖16與圖18a、圖17與圖18b的對(duì)比可知,構(gòu)建的表面形貌模型能夠比較客觀地反映實(shí)際機(jī)械加工面的形貌。通過改變模型參數(shù),可以快速得到不同的表面形貌模型,為后續(xù)的界面接觸、潤(rùn)滑等分析提供數(shù)據(jù)。

(a)各向同性加工面

3.3 擺線針輪傳動(dòng)混合潤(rùn)滑分析算例

在擺線針輪傳動(dòng)的混合潤(rùn)滑分析中應(yīng)用本文提出的模型進(jìn)行數(shù)值計(jì)算,主要參數(shù)如表1所示。仿真分析采用朱東等開發(fā)的MixedEHL軟件[26]。

表1 擺線針輪傳動(dòng)混合潤(rùn)滑分析基本參數(shù)Table 1 Parameters of cycloid-pin gears and grease

圖19為應(yīng)用各向異性粗糙面模型計(jì)算得到的潤(rùn)滑脂膜厚度及壓力分布情況,網(wǎng)格密度為256×256。采用本文模型,能完成相應(yīng)的粗糙面潤(rùn)滑接觸計(jì)算,得到的膜厚及壓力分布與實(shí)際一致,表明了模型的有效性。

(a)壓力分布

(a)膜厚分布對(duì)比曲線

圖20為采用實(shí)際掃描粗糙面數(shù)據(jù)與采用本文模型得到的x方向的膜厚及壓力曲線,其中Hm、Hr和Pm、Pr分別為應(yīng)用本文模型和實(shí)際掃描粗糙面計(jì)算得到的無量綱膜厚和無量綱壓力。從圖中可看出,本文模型與實(shí)際掃描粗糙面的計(jì)算結(jié)果一致,最小膜厚與最大壓力誤差分別為2.06%與1.08%,本文模型的計(jì)算時(shí)間為18 629.062 5 s,實(shí)際掃描粗糙面的計(jì)算時(shí)長(zhǎng)為42 014.453 1 s,計(jì)算效率提升55.7%。此外,本文模型還可以通過修改參數(shù),快速得到不同特征粗糙面,并進(jìn)行相應(yīng)的仿真對(duì)比,在實(shí)際應(yīng)用中有明顯的優(yōu)勢(shì)。

4 結(jié) 論

為方便對(duì)機(jī)械結(jié)合面進(jìn)行相應(yīng)的接觸、疲勞、潤(rùn)滑等數(shù)值分析,本文提出了一種可以綜合體現(xiàn)機(jī)械結(jié)合面宏觀輪廓與微觀形貌的三維粗糙曲面模型。該模型能客觀表征任意輪廓的機(jī)械零件曲面,且能有效表征粗糙面的紋向、粗糙度等微觀形貌。針對(duì)分形粗糙面模型處處不可微導(dǎo)致計(jì)算不易收斂的問題,采用二維高斯濾波對(duì)模型數(shù)據(jù)進(jìn)行平滑處理,提高了粗糙面參與數(shù)值計(jì)算時(shí)的收斂性。本文提出的模型可以用于機(jī)械結(jié)合面的接觸、潤(rùn)滑以及摩擦磨損分析,模型體現(xiàn)了結(jié)合面的宏觀輪廓與微觀形貌信息,可以通過對(duì)比分析,為機(jī)械結(jié)合面的設(shè)計(jì)參數(shù)、加工工藝參數(shù)等的優(yōu)化提供指導(dǎo)。