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

喬澤龍,韓炬,董威

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

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

機械加工面的粗糙形貌與加工工藝密切相關,近年來,眾多學者結合刀具尺寸與傾角、加工路徑以及振動等因素,實現了對機械加工工件表面形貌的預測[6-7],但一般用于工件表面質量評測,沒能應用于機械結合面的接觸、潤滑、磨損等特性分析。

現有獲取粗糙面的三維形貌主要包括試驗法與數值法兩種方式。目前,表面形貌的測試方法很多[8],但受到橫縱坐標分辨率的限制,測試獲取的形貌模型在很多細節上被簡化,此外測試得到的模型只是對現有加工面的再現,很難對獲取改性的曲面形貌進行對比研究。

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

本文綜合應用微分幾何理論與分形理論,構建能表征復雜形狀曲面的宏觀輪廓與微觀形貌的三維模型,通過對模型進行平滑處理,提升了模型在諸如接觸、疲勞、潤滑等數值分析時的收斂性。

1 機加工粗糙曲面表面形貌模型

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

一般應用W-M函數為[19]對粗糙面微觀形貌進行分形模擬。W-M分形函數具有連續、處處不可導和自仿射的特點,表達式為

(1)

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

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

(2)

1.2 形貌各向異性和各向同性表征函數

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

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

(3)

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

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

(4)

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

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

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

(5)

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

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

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

根據式(4),可得各向異性三維粗糙仿真曲面,如圖2所示,其中x和y方向特征尺度系數和分形維數相同,特征尺度系數取Gx=Gy=1.36×10-5mm,分形維數分別取Dx=1.488,Dy=1.518。采樣長度為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 各向同性分形粗糙表面仿真

根據式(5),可得到圖4所示的各向同性三維粗糙仿真曲面,取樣參數與各向異性相同,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 分形粗糙表面的濾波處理

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

圖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

采用高斯濾波對分形粗糙形貌高度進行平滑時,中心位置點高度為包括本點在內的鄰域各位置高度數據按高斯分布規律的加權和,在實現對數據平滑的同時,能夠保留形貌高度的原始總體分布特征。二維高斯濾波器的核心為高斯二維卷積算子,本文選取高斯核為9×9,σ=1.8。

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

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

(a)原始數據仿真

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

雖然機械動態結合面的結合區很小,接觸瞬時的曲面可以近似為粗糙平面,這是眾多粗糙面模型將接觸區簡化為平面接觸的重要原因。但結合面宏觀輪廓,如齒輪齒廓、凸輪輪廓等,在傳動過程中對傳動性能有關鍵影響,因此三維粗糙曲面的形貌模型應該綜合體現曲面的宏觀輪廓與微觀形貌,即將體現微觀形貌的W-M函數與體現宏觀形貌的矢量函數疊加。

本節以RV減速器中擺線輪輪齒齒廓為例,創建三維粗糙曲面模型。RV減速器擺線針輪嚙合接觸屬于柱面接觸,擺線針輪接觸仿真計算需要基于擺線輪曲面構建粗糙曲面,由光滑曲面與非均勻的微凸體疊加構成粗糙曲面形貌模型,其中光滑曲面由矢量函數表示,微凸體的高度值由W-M函數計算,微凸體高度方向為相應曲面點的法矢量方向。擺線輪齒的粗糙表面可表示為擺線輪齒齒面的矢量函數和粗糙面矢量函數之和。

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

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

(6)

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

(7)

(8)

(9)

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

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

(10)

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

(11)

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

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

(12)

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

(13)

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

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

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

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

(14)

擺線矢量方程為

(15)

修形后的擺線輪齒的齒廓矢量方程可轉化為

(16)

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

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

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

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

擺線輪齒齒廓三維建模的核心問題是如何將理論光滑齒輪曲面模型與粗糙面模型相結合。粗糙面建模可依據式(4)和式(5)獲得模型數據,并采用二維高斯濾波對數據進行平滑。理論齒廓三維建模可采用圓柱坐標矢量函數,獲得齒廓三維曲面,再通過數值計算獲得曲面上點的單位法向矢量。各曲面法向矢量以粗糙曲面相應點高度值為模長,獲得粗糙曲面矢量矩陣,理論齒面矢量與粗糙曲面矢量之和即為粗糙齒廓的矢量模型。

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

(17)

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

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

在齒面法向單位矢量計算中,由于曲面z方向的導數為0,故法向矢量表達式為

(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)為各向異性或各向同性粗糙曲面矢量。依據式(19),按經驗取特征尺度系數和分形維數,采用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

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

(a)各向同性加工面

3.3 擺線針輪傳動混合潤滑分析算例

在擺線針輪傳動的混合潤滑分析中應用本文提出的模型進行數值計算,主要參數如表1所示。仿真分析采用朱東等開發的MixedEHL軟件[26]。

表1 擺線針輪傳動混合潤滑分析基本參數Table 1 Parameters of cycloid-pin gears and grease

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

(a)壓力分布

(a)膜厚分布對比曲線

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

4 結 論

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