基于分形理論的平面磨削粗糙表面接觸剛度的研究

珠海城市職業技術學院 機電工程學院,廣東 珠海 519090

機械加工表面具有一定的粗糙度特征，這種粗糙度特征將直接影響機械零件之間的連接性能，對機械系統的穩定性和可靠性有著重要的影響[1,2]。隨著機械產品精度的提高，機械系統中較為重要的連接平面多數采用平面磨削的加工方式。接觸剛度作為描述粗糙表面接觸行為的一個重要參數，若能夠準確地預測平面磨削表面的接觸剛度對于解決工程實際問題具有重要意義[3,4]。

隨著有限元技術的發展，有限元接觸仿真方法被廣泛的應用于科學研究的各個領域。近幾年也被引入到粗糙表面接觸剛度的研究上來。就目前的研究成果來看，對于粗糙表面接觸剛度的仿真主要存在兩個主流方向：（1）采用高斯或非高斯表面表征粗糙表面，繼而對粗糙表面進行接觸仿真分析。WU 等[5,6]基于傅里葉變換生成高斯或非高斯表面對磨削粗糙表面進行了模擬，基于特征表面對粗糙表面接觸剛度進行了相關的研究。WANG 等[7]在后續研究工作中結合傅里葉變換與JOHNSON轉換系統[8]對具有一定特征參數的非高斯表面進行了模擬，討論了相關粗糙度參數對于接觸剛度的影響規律。（2）采用分形表面表征，繼而對粗糙表面進行接觸分析。YAN 和KOMVOPOLOS[9]首先提出采用三維的W-M 函數構造粗糙表面，并基于三維分形表面對粗糙表面接觸剛度進行了相關的研究，得到了分形維數D與尺度系數G對接觸剛度的影響規律[10]。國內西安理工大學的傅衛平團隊[11,12]對基于分形表面接觸剛度進行了較為系統的研究。

在對平面磨削表面實測形貌的研究過程中發現，上述方法在一定程度上能夠實現對于特定粗糙表面接觸剛度的預測。但是對比實測平面磨削的采集形貌來看，上述表征方法無法實現對平面磨削表面接觸剛度的針對性準確預測。因此采用上述表征方法得到的平面磨削表面接觸剛度存在一定的誤差。

針對上述問題，本文針對平面磨削表面接觸剛度的仿真，提出一種新的途徑與方法。本方法從平面磨削表面形貌的表征出發，采用分形表面模擬砂輪表面形貌，結合磨削運動學與切削理論實現對平面磨削表面的形貌仿真，繼而基于有限元接觸仿真實現了平面磨削表面接觸剛度結果的預測。

1 平面磨削表面的生成

對于平面磨削表面的生成，本文從平面磨削加工用砂輪的表面進行分析。采用采用分形表面模擬砂輪表面形貌，結合磨削運動學與切削理論即可實現對平面磨削表面的形貌仿真。首先針對砂輪表面的生成方法進行探討。

1.1 砂輪表面的生成

本節采用三維W-M 函數實現對砂輪表面形貌進行表征，通過控制函數中相關參數的變化，即可實現對不同粗糙度參數下砂輪表面的形貌表征。YAN 和KOMVOPOLOS[12]給出了直角坐標系下各同向性的三維隨機分形表面的W-M 函數，如式(1)所示。

其中，z(x,y)為仿真表面的輪廓高度；Ds為仿真表面的分形維數，并且有21；L為仿真表面的取樣長度；M為仿真表面褶皺的重疊數；m,n為隨機相位，取值范圍為[0,2π]；nl為最低頻率，通常取nl=0；nmax為最高頻率，表達式為nmax=int(log(L/Ls)/logγ)；Ls為截止長度。

采用MATLAB 軟件進行程序的編制，即可實現對三維分形表面的模擬仿真。仿真流程如下所述：(1)確定仿真區域。對仿真區域進行插值處理，以此獲得256×256 矩陣作為仿真變量X與Y的初始值。(2)設定計算參數的初始值。初始計算參數包括：M=10；γ=1.5；L=1×10-3m；Ls=5×10-9m；D取2.5；G=10-8；φm,n=π/6?；谏鲜鰠涤嬎鉵max的數值，得到n的取值范圍。(3)將m,n作為循環變量，計算模擬表面的輪廓高度。(4)顯示圖像。同時進行數據的截取，進行圖像的顯示。最后進行數據的保存，為下一步基于分形表面的平面磨削表面形貌的獲取做好準備。由上述計算參數所得到三維分形粗糙表面如圖1 所示。

圖1 D=2.5 下三維分形表面Fig.1 3D fractal surface at D=2.5

圖2 D=2.5,G=10-8 條件下的平面磨削表面形貌Fig.2 Surface morphology of surface grinding at D=2.5 and G=10-8

經由上述分析，可以獲得分形維數D為2.5、尺度系數G=10-8條件下的三維分形表面通過改變仿真的初始參數，可以實現對不同參數下三維分形表面的模擬。

通過編制相應程序，實現了對平面磨削用砂輪表面的形貌仿真，并得到了砂輪表面形貌的數據。下節將基于砂輪表面的形貌數據，結合磨削運動學與切削理論實現對平面磨削表面形貌的仿真。

1.2 平面磨削表面的生成

將基于砂輪表面的形貌數據，結合磨削運動學與切削理論實現對平面磨削表面形貌的仿真。在砂輪表面形貌數據的基礎上，結合砂輪表面所有磨粒的運動軌跡即可實現對平面磨削表面理論形貌的仿真。對于平面磨削表面理論形貌的仿真，文獻[13]中對于砂輪形貌的仿真流程及實現方法進行了詳細的介紹說明，這里不再贅述。本文按照文獻所述方法編制相應程序，即可實現對平面磨削表面理論形貌的仿真。這里對于平面磨削的加工參數取值包括：砂輪轉速n=1460 r/min，工件的進給速度為vw=5 m/s，磨削進給量ap=5 μm。圖2 所示為分形維數D為2.5、尺度系數G=10-8條件下的平面磨削理論形貌。

基于砂輪表面形貌數據，結合磨削運動學與切削理論，實現了平面磨削表面的形貌仿真，繼而得到了平面磨削表面的形貌數據。

2 基于ABAQUS 的平面磨削表面接觸剛度仿真

將對基于表面形貌數據的平面磨削接觸剛度獲取流程進行介紹。整體仿真流程可以劃分為兩個部分：基于表面形貌數據的三維實體生成流程以及有限元接觸仿真流程。下面將針對上述兩個流程進行詳細的介紹說明。

2.1 平面磨削表面三維實體的生成

對基于表面形貌數據的三維實體生成進行詳細介紹。首先，通過第1.2 分析，獲得了平面磨削表面的形貌數據，采用逆向建模軟件生成特征曲面。本文選用MESHLAB 逆向建模軟件將導入的表面形貌點云數據形成網格數據，采用滾球法實現了特征曲面的生成。圖3 所示為MESHLAB 軟件中D=2.5，G=10-8條件下通過表面形貌數據生成的特征曲面。

圖3 D=2.5,G=10-8 條件下的特征曲面Fig.3 Characteristic surface at D=2.5,G=10-8

圖4 D=2.5,G=10-8 條件下的三維實體Fig.4 3D model at D=2.5,G=10-8

將上述特征曲面導入到逆向建模軟件Pro-E 當中，采用Pro-E 的小平面特征命令進行三維建模。通過曲面剪裁、實體化命令，即可實現對三維實體的生成。圖4 所示為D=2.5，G=10-8條件下通過特征曲面生成的三維實體模型。

2.2 平面磨削表面接觸剛度的有限元仿真

采用ABAQUS 軟件對2.1 中生成的三維實體進行有限元接觸仿真。下面將對平面磨削表面接觸剛度有限元仿真流程進行詳細介紹。將2.1 得到的三維實體導入到ABAQUS 軟件當中。導入模型之后，為保證后期分析計算的效率，需要對單個模型進行拆分，以實現針對粗糙表面的局部網格細化的作用。本文將實體分割為三個部分。之后，對實體模型的材料以及截面屬性進行定義。首先將粗糙表面的接觸類型設置成彈塑性接觸，輸入彈塑性變形區間的材料參數。本文選用45#鋼為研究對象，對應的彈性變形區間的材料參數分別為：彈性模量E=210 GPa，泊松比ν=0.269，屈服強度σy=355 MPa。對于材料塑型變形區間的參數設置如表1 所示。

表1 45#鋼真實應力與塑性應變的對應曲線Table 1 The corresponding curves of true stress and plastic strain for 45#steel materials

在完成材料屬性定義之后，將接觸分析的兩個部件進行裝配。在完成兩個部件的裝配之后，下一步則是對部件進行網格劃分。由于在導入部件環節，已經對單個部件的三部分拆分，因此對于單個部件三部分的網格大小不同。本文分別將三部分的網格大小設置為0.002 μm，0.005 μm 以及0.01 μm。網格類型選用C3D8R 六面體網格單元。圖5 所示為裝配體網格劃分示意圖。

在上述操作完成之后，對分析步進行設置。將初始增量步的大小設置為1E-30；最小增量步大小設置為1E-30；最大增量步大小設置為1E-2；最大增量步數設置為50000。在上述設置完成之后，需要對場輸出以及歷程輸出進行設置。由接觸剛度的定義可知，只需選取接觸壓力與接觸位移進行輸出即可。接觸類型選為部件1 的粗糙表面與部件2 的粗糙表面的面面接觸形式。面面接觸的相互作用屬性選切向行為，摩擦公式選擇罰函數，摩擦系數設置為0.15。給定邊界條件。對于載荷加載方式的選擇，本文采用兩端（光滑表面）位移加載的方式。待上述設置完成之后，即可創建作業，提交運算。圖6 所示為D=2.5，G=10-8條件下接觸仿真應力云圖。待分析結束，打開運行結果，選擇場輸出請求中設置的面集的反作用力與位移兩個參數，輸出參數的具體數值。經過上述的分析處理，即可得到平面磨削表面接觸剛度的有限元接觸仿真結果。

圖5 裝配體網格劃分示意圖Fig.5 Schematic diagram of assembly grid division

圖6 D=2.5,G=10-8 條件下接觸仿真應力云圖Fig.6 Contact simulation stress nephogram at D=2.5,G=10-8

以D=2.5，G=10-8條件下平面磨削表面形貌為例對基于表面形貌數據的平面磨削接觸剛度獲取流程進行介紹，得到了接觸剛度的有限元仿真結果。

3 仿真結果與分析

經由上述分析，可以得到平面磨削表面接觸剛度的仿真結果。本節將通過改變砂輪表面的分形參數，實現對不同平面磨削表面接觸剛度有限元結果的對比分析。按照上述處理流程，分別對砂輪表面分形維數條件下獲得的平面磨削表面進行接觸剛度有限元仿真。圖7 所示為尺度系數G=10-8，不同分形維數D條件下，平面磨削表面接觸剛度隨接觸壓力的變化曲線。分形維數分別取2.1、2.3、2.5、2.7 與2.9 共5 組。

圖7 不同分形維數下接觸剛度變化曲線Fig.7 Change curves ofcontactstiffnessunderdifferentfractaldimensions

圖8 不同尺度系數下接觸剛度變化曲線Fig.8 Changecurves of contactstiffnessunderdifferentscale coefficients

基于上圖分析，可以得到以下兩個結論：（1）平面磨削表面法向接觸剛度隨著接觸壓力的增大而逐漸增大，兩者之間呈現非線性關系，并且在不同分形維數D的條件下，兩者具有相同的變化趨勢。（2）在形同接觸壓力的條件下，平面磨削表面接觸剛度隨著分形維數D的增大而逐漸增大。并且在相同接觸壓力的條件下，不同分形維數對應的接觸剛度的差值也隨著接觸壓力的增大而逐漸增大。在接觸壓力為140 MPa 的條件下，不同分形維數下（2.1、2.3、2.5、2.7 與2.9）接觸剛度的差值達到最大，此時不同分形維數對應的接觸剛度分別為42.566 MPa/μm、58.881 MPa/μm、68.532 MPa/μm、83.977 MPa/μm、94.356 MPa/μm。

圖8 所示為分形維數D=2.1，不同尺度系數G條件下，平面磨削表面接觸剛度隨接觸壓力的變化曲。尺度系數分別取10-8、10-9、10-10與10-11共4 組。基于對圖8 的分析同樣可以得到以下兩個結論：（1）平面磨削表面法向接觸剛度隨著接觸壓力的增大而逐漸增大，兩者之間呈現非線性關系，并且在不同分形維數G的條件下，兩者具有相同的變化趨勢。（2）在形同接觸壓力的條件下，平面磨削表面接觸剛度隨著尺度系數G的增大而逐漸減小。并且在相同接觸壓力的條件下，不同尺度系數對應的接觸剛度的差值也隨著接觸壓力的增大而逐漸增大。在接觸壓力為140 MPa 的條件下，不同分形維數下（10-8、10-9、10-10與10-11）接觸剛度的差值達到最大，此時不同分形維數對應的接觸剛度分別為58.285 MPa/μm、83.301 MPa/μm、100.363 MPa/μm、133.638 MPa/μm。

本節通過改變砂輪表面的分形參數，采用本文方法獲得了不同的平面磨削表面。針對不同的平面磨削表面進行了有限元接觸分析，進而得到了不同參數條件下接觸剛度隨接觸壓力的變化關系。最后結合接觸剛度的仿真結果，將分形參數對平面磨削表面接觸剛度的影響規律進行了概括性分析。

4 結論

（1）本文采用分形表面對砂輪的表面形貌進行表征，結合機床運動學與切削理論，實現了平面磨削表面形貌的仿真。采用有限元方法對平面磨削仿真表面進行接觸分析，得到平面磨削表面接觸剛度的仿真結果；

（2）對不同分形參數砂輪表面條件下獲得的平面磨削表面接觸剛度進行仿真分析。可以得到以下兩個結論：首先，平面磨削表面接觸剛度的數值隨接觸壓力的增大而逐漸增大，并且兩者之間呈現非線性關系。不同參數條件下，接觸剛度隨接觸壓力的變化趨勢相同。其次，通過對比相同接觸壓力條件下，平面磨削表面接觸剛度隨分形維數的增大而逐漸增大，隨尺度系數的增加而逐漸減小。