外圓磨削砂輪形貌仿真與工件表面粗糙度預(yù)測

呂長飛 李郝林

1.上海理工大學(xué)，上海，200093 2.河北農(nóng)業(yè)大學(xué)，保定，071001

0 引言

磨削加工往往是機械產(chǎn)品的終極加工工序，其加工效果直接影響到產(chǎn)品的最終質(zhì)量和性能。在很多工業(yè)國家，磨削已成為至關(guān)重要的工藝過程［1］，很多研究者對磨削尤其是磨削過程的仿真和建模進(jìn)行了深入的研究并取得了顯著的成績［2］。仿真方法總體包括兩類：根據(jù)經(jīng)驗的方法［3］和理論分析的方法［4］，這兩類方法的聚焦點都在磨削砂輪形貌的仿真和磨削運動過程，即工件和磨粒相互作用過程的仿真。

Chakrabarti等［5］假設(shè)磨粒為錐體，且均勻分布在砂輪表面，在此基礎(chǔ)上實現(xiàn)了對工件表面粗糙度的預(yù)測，但均勻分布跟砂輪形貌實際不相符。Aurich等［6］在考慮了磨削偏差和熱效應(yīng)的條件下，研究了存在顫振時磨削粗糙度的預(yù)測。Chen等［7］在假設(shè)磨粒為圓球狀并隨機分布在砂輪表面的條件下，比較全面地對磨粒形貌和磨削加工運動過程進(jìn)行了仿真，但只考慮了磨粒的切削作用，忽略了磨粒的耕犁與摩擦作用。Pandit等［8］將時間序列理論引入砂輪形貌建模中，采用自回歸或自回歸滑動平均擬合砂輪磨削表面數(shù)據(jù)，利用時間序列對砂輪表面建模，建立了二維的砂輪輪廓模型，其計算量很大，非常耗時。Hou等［9］假設(shè)不同大小的磨粒分布在砂輪表面，在確定所有磨粒大小和分布的前提下，采用隨機方法得到了砂輪形貌的近似模型。Koshy等［10］采用與Hou等建立模型相類似的方法，采用統(tǒng)計方法估計磨粒突出部分高度，提出磨粒大小滿足正態(tài)分布，并采用與Chen相似的磨削加工過程分析方法建立了磨削模型，但采用隨機方法同樣存在計算量的問題。Zhou等［11］則認(rèn)為磨粒突出部分滿足隨機域的高斯分布，并給出了單磨粒的運動軌跡，將其映射到加工工件表面坐標(biāo)系，然后推導(dǎo)了兩個磨粒軌跡的交點方程，繼而得到了兩個磨粒作用下的工件表面形貌，對工件上的每一點取所有磨粒單獨磨削生成的砂輪形貌的最小值，最終得到整個砂輪磨削加工后的工件表面形貌。目前更受肯定的觀點是磨粒突出部分滿足非高斯分布。文獻(xiàn)［12－13］提出的砂輪表面的磨粒高度分布大部分屬于非高斯隨機分布，采用映射的方法，針對平面磨削建立了工件表面形貌模型。對于外圓磨削的研究，文獻(xiàn)［14－17］均通過各種傳感器和儀器測量實驗數(shù)據(jù)，采用計算機智能算法對實驗數(shù)據(jù)進(jìn)行系統(tǒng)分析，得到了磨削模型，這些都是基于經(jīng)驗的方法。Hecker等［18］假設(shè)磨粒為理想圓錐，建立了磨削碎片無變形模型，通過確定工件表面粗糙度與碎片厚度的關(guān)系得出粗糙度預(yù)測模型，并通過外圓磨削實驗數(shù)據(jù)進(jìn)行了驗證，磨粒的形貌假設(shè)為理想圓錐雖然簡化了計算，但與砂輪實際形貌不太相符。

本文針對Nguyen等［12］提出的磨粒高度為非高斯隨機表面理論，采用Johnson變換和Gabor小波變換以及相應(yīng)的反變換，實現(xiàn)高斯域和非高斯域的轉(zhuǎn)化，在隨機域內(nèi)對磨削砂輪形貌進(jìn)行了仿真。根據(jù)外圓磨削運動過程，通過砂輪和工件相互作用過程的分析，建立磨粒運動軌跡方程和工件形貌方程，在考慮磨粒切削、耕犁與摩擦作用的條件下，對外圓磨削過程進(jìn)行了仿真。建立了外圓磨削模型，實現(xiàn)對加工工件形貌的仿真和粗糙度預(yù)測。最后通過實驗進(jìn)行了驗證。

1 砂輪形貌仿真

砂輪形貌是影響磨削過程的重要因素，砂輪磨粒的分布和形狀是磨削碎片產(chǎn)生和工件表面形貌形成最關(guān)鍵的因素，因此，對砂輪形貌的仿真尤為重要，它是實現(xiàn)磨削過程仿真的關(guān)鍵。目前最受肯定的論述是砂輪表面為隨機面，且為非高斯分布。對整個砂輪表面進(jìn)行非高斯隨機域的仿真，數(shù)據(jù)量大，非常耗時，因此非高斯域的表達(dá)和高斯域的仿真計算成為必然，選取合適的函數(shù)變換，實現(xiàn)高斯域和非高斯域的相互轉(zhuǎn)換成為必須解決的問題。Nguyen等［12］提出并論證了對于任何變換函數(shù)，滿足變換和其反變換是一對一函數(shù)，并在相關(guān)域內(nèi)其Jacobian變換非奇異，就能實現(xiàn)高斯和非高斯函數(shù)在隨機域內(nèi)的轉(zhuǎn)換。

Gabor變換函數(shù)又稱為Gabor濾波器，是一種將高斯函數(shù)作為窗函數(shù)的短時傅里葉變換，最早由 Gabor［19］于1946年提出。Daugman［20－21］最早提出了二維Gabor濾波器，并詳細(xì)論述了它的數(shù)學(xué)特性［22］。二維Gabor基函數(shù)的通用表示形式為

式中，w（x，y）為復(fù)正弦載波；P為載波的相位，在實際應(yīng)用中，通常有載波相位P＝0；（F0，0）為空間頻域中心；s（x，y）為Gaussian包絡(luò)函數(shù)；K 為一常數(shù)；a、b為系數(shù)。

二維Gabor基函數(shù)g（x，y）的傅里葉變換為

具有帶通性質(zhì)的單一Gabor濾波器只能在有效響應(yīng)范圍內(nèi)（頻域響應(yīng)幅值高于某個設(shè)定值的區(qū)域）覆蓋某一個橢圓區(qū)域，要想覆蓋整個頻域面，需構(gòu)造一個Gabor濾波器庫，即構(gòu)造具有不同尺度的Gaussian包絡(luò)與中心頻率的多個Gabor帶通濾波器來達(dá)到上述目的。則有

式中，θi為Gaussian包絡(luò)的旋轉(zhuǎn)角度；i為旋轉(zhuǎn)角度索引；j為中心頻率索引；F0，ij為 中 心頻率；φi為中 心 頻 率F0，ij的幅角，且一般情況下有θi＝φi。

對比式（1）、式（3）可以發(fā)現(xiàn)，每個gij（x，y）都可由g（x，y）經(jīng)一定的尺度伸縮和旋轉(zhuǎn)變換得到：

對g（x，y）進(jìn)行直流分量補償后可以發(fā)現(xiàn)g（x，y）即為小波理論中的母小波函數(shù)。對g（x，y）進(jìn)行直流分量補償是為了使之滿足小波理論中的容許性條件［23］，經(jīng)補償后的g（x，y）為

此時g（x，y）的傅里葉變換為

此時的傅里葉變換為

這樣，gij（x，y）也被稱為Gabor小波。

參數(shù)F0，ij、λj和aj、bj的選取原則是基于具體應(yīng)用對象實現(xiàn)相關(guān)參數(shù)的最優(yōu)選?。?4］。

對砂輪表面取樣部分磨粒高度函數(shù)h′（x，y）進(jìn)行Gabor小波變換，可得仿真磨粒高度函數(shù)h（x，y）：

對應(yīng)的傅里葉運算為

式中，H（u，v）、H′（u，v）分別為h（x，y）和h′（x，y）的傅里葉變換。

Johnson變換在工業(yè)領(lǐng)域很大范圍內(nèi)可實現(xiàn)高斯表面和非高斯表面之間的轉(zhuǎn)換，其一般表達(dá)式為

式中，Q為高斯分布隨機變量；X為非高斯分布隨機變量；γ、δ為決定X的分布形狀的變量；ξ為位置因子；λ為比例因子。

根據(jù)不同的偏差Sk和偏度Ku的要求，Johnson變換有以下幾類：

其中，參數(shù)γ、δ、ξ和λ可由靜態(tài)統(tǒng)計表［25］或通過Hill等［26］給出的數(shù)值方法確定。本文采用式（13）進(jìn)行計算，其中參數(shù)由Hill等給出的數(shù)值方法確定。

2 磨削過程分析

當(dāng)認(rèn)為磨床工件和砂輪轉(zhuǎn)動為標(biāo)準(zhǔn)圓時，外圓磨削單個磨粒運動軌跡過程等同于平面磨削，如圖1所示。

o′為磨粒與工件接觸最低點，并作為坐標(biāo)系o′xyz的原點，F(xiàn)o′F′ 為切削點擺動軌跡，C 為砂輪軸中心，ds為砂輪直徑，磨削深度為ap，vs、vw分別為砂輪和工件表面線速度，θ為砂輪在單位時間內(nèi)的旋轉(zhuǎn)角，由圖1可知：

圖1 外圓磨削單一磨粒運動軌跡

圖1中工件與砂輪相向運動，此時當(dāng)磨粒往上運動時式（16）中取“＋”，往下運動時取“－”。

因θ很小，可認(rèn)為θ＝sinθ聯(lián)立式（16）、式（17）可得

考慮所有磨粒，對每一個磨粒都有

式中，xij、zij為工件表面切削點在o′xyz坐標(biāo)系中的位置；dij為砂輪中心到磨粒頂端的距離。

顯然，dij＝ds＋hij，其中hij為砂輪表面對應(yīng)磨粒高度（圖2），hij可由上述Gabor小波變換或Johnson變換獲得。

圖2 外圓磨削過程坐標(biāo)圖

3 工件表面仿真

工件形貌采用拓?fù)渚仃噂ij表示，矩陣的每個元素gmn表示在全局坐標(biāo)系oxyz中的工件oxy平面中點（m，n）的高度。全局坐標(biāo)系的坐標(biāo)原點o的選擇應(yīng)滿足max｛gmn＝0｝，即平面中工件上所有點均在全局坐標(biāo)點下方，x、y、z軸向上，如圖2所示。對砂輪進(jìn)行離散化并用矩陣hij表示，i、j表示高度為hij的砂輪周向和軸向位置。圖2中每一切削點o′在全局坐標(biāo)中的周向等效距離近似為

式中，Δx為x方向的改變量。

則每個切削點的位置是

其中，rij為在z軸方向o′到o的距離。因此，工件上的點（m，n）的拓?fù)渚仃嚨拿恳辉乜杀硎緸?/p>

其中，Δxw為工件拓?fù)渚仃噂ij沿x軸的間距。則在經(jīng)過切點hij切削后，工件點（m，n）的拓?fù)渚仃囍忻恳辉乜杀硎緸?/p>

對所有磨粒，當(dāng)其磨粒高度滿足hij≥hmaxa時，將產(chǎn)生實際磨削過程。磨削過程包括切削、耕犁或摩擦三種不同狀態(tài)，可由磨削過程產(chǎn)生側(cè)脊作用來反映，并由磨粒切入角α確定［12］：

其中，Ar、h、l為側(cè)脊面積、高度和寬度，Ar與切槽面積成比例，比例系數(shù)為有效切除率β［7，12］；r為磨粒切深；R為假設(shè)磨粒是球形時的磨粒半徑，κ為其曲率。磨削側(cè)脊圖見圖3。

圖3 磨削側(cè)脊圖

4 實驗驗證

實驗在Schleifring公司的KC33數(shù)控萬能內(nèi)外圓磨床上進(jìn)行，砂輪為氧化鋁砂輪（A60K7V），工件為45鋼，砂輪直徑ds＝424mm，磨削深度ap＝20μm，工件速度vw＝0.3123m／s，砂輪速度vs＝35m／s。實驗和仿真結(jié)果如圖4、圖5所示。

所生成砂輪表面的均值為零，方差為0.8，偏度為±0.25，峰度為4.1。對應(yīng)不同磨削加工參數(shù)，基于仿真模型得到的工件表面粗糙度值與實驗測量值如表1所示。

圖4 砂輪形貌仿真

表1 表面粗糙度預(yù)測值與實驗測量值

由表1可知，預(yù)測模型所得表面粗糙度與實驗測量值大體相近，且變化趨勢一致，這為實現(xiàn)磨削過程預(yù)測和磨削智能化提供了依據(jù)。

5 結(jié)語

本文采用Johnson變換和Gabor小波變換以及相應(yīng)的反變換，實現(xiàn)高斯域和非高斯域的轉(zhuǎn)化，通過對磨粒高度函數(shù)的分析建模，在隨機域內(nèi)對磨削砂輪形貌進(jìn)行了仿真。針對外圓磨削運動過程，通過砂輪和工件相互作用過程的分析，建立磨粒運動軌跡方程和工件形貌方程。再假設(shè)側(cè)脊為等腰三角形，通過分析磨粒切入角來考慮磨粒切削、耕犁與摩擦的作用，實現(xiàn)對外圓磨削過程的仿真。建立了外圓磨削模型，得到工件形貌與砂輪參數(shù)（砂輪直徑、磨粒高度等）和加工參數(shù)（切削深度、砂輪和工件轉(zhuǎn)速、磨削切入角）的關(guān)系方程，實現(xiàn)了對加工工件形貌的仿真和粗糙度預(yù)測，并通過實驗進(jìn)行了驗證。該預(yù)測模型為磨削參數(shù)的選擇、磨削過程的智能化提供了重要依據(jù)，具有重要的現(xiàn)實意義。

［1］Doman D A，Warkentin A，Bauer R.A Survey of Recent Grinding Wheel Topography Models［J］.Machine Tools ＆ Manufacture，2006，46（3／4）：343－352.

［2］Komanduri R.Machining and Grinding：a Historical Review of the Classical Papers［J］.Applied Mechanics Reviews，1993，46（3）：80－132.

［3］Badger J A，Torrance A A.A Comparison of Two Models to Predict Grinding Forces from Wheel Surface Topography［J］.International Journal of Machine Tools and Manufacture，2000，40（8）：1099－1120.

［4］Tonshoff H K，Peters J，Inasaki I，et al.Modelling and Simulation of Grinding Processes［J］.CIRP Annals－Manufacturing Technology，1992，41（2）：677－688.

［5］Chakrabarti S，Paul S.Numerical Modelling of Surface Topography in Superabrasive Grinding［J］.Int.J.Adv.Manuf.Technol.，2008，39（1／2）：29－38.

［6］Aurich J C，Biermann D，Blum H，et al.Modelling and Simulation of Process：Machine Interaction in Grinding［J］.Prod.Eng.Res.Devel.，2009，3（1）：111－120.

［7］Chen X，Rowe W B.Analysis and Simulation of the Grinding Process.PartⅡ：Mechanics of Grinding［J］.International Journal of Machine Tools and Manufacture，1996，36（8）：883－896.

［8］Pandit S M，Wu S M.Characterization of Abrasive Tools by Continuous Time Series［J］.ASME Journal of Engineering for Industry，1973，95（3）：821－826.

［9］Hou Z B，Komanduri R.On the Mechanics of the Grinding Process－Part I：Stochastic Nature of the Grinding Process［J］.International Journal of Machine Tools and Manufacture，2003，43（15）：1579－1593.

［10］Koshy P，Jain V K，Lal G K.A Model for the Topography of Diamond Grinding Wheels［J］.Wear，1993，169（2）：237－242.

［11］Zhou X，Xi F.Modeling and Predicting Surface Roughness of the Grinding Process［J］.International Journal of Machine Tools and Manufacture，2002，42（8）：967－977.

［12］Nguyen T A，Butler D L.Simulation of Precision Grinding Process.Part I：Generation of the Grinding Wheel Surface［J］.International Journal of Machine Tools Manufacturing，2005，45（8）：1321－1328.

［13］陳東祥，田延嶺.超精密磨削加工表面形貌建模與仿真方法［J］.機械工程學(xué)報，2010，46（13）：186－191.

［14］Lee C W.A Control－oriented Model for the Cylindrical Grinding Process［J］.Int.J. Adv.Manuf.Technol.，2009，44（7／8）：657－666.

［15］Chiu N，Malkin S.Computer Simulation for Cylindrical Plunge Grinding［J］.Annals of the CIRP，1993，42（1）：383－387.

［16］Choi T J，Subrahmanya N，Li H，et al.Generalized Practical Models of Cylindrical Plunge Grinding Processes［J］.International Journal of Machine Tools ＆ Manufacture，2008，48（1）：61－72.

［17］Shih A J，Lee N L.Precision Cylindrical Face Grinding［J］.Precision Engineering，1999，23（3）：177－184.

［18］Hecker R L，Liang S Y.Predictive Modeling of Surface Roughness in Grinding［J］.International Journal of Machine Tools ＆ Manufacture，2003，43（8）：755－761.

［19］Gabor D.Theory of Communication［J］.Inst.E－lect.Eng.，1946，93（3）：429－457.

［20］Daugman J G.Complete Discrete 2－D Gabor Transforms by Neural Networks of Image Analysis and Compression［J］.IEEE Transaction on A－coustics Speech＆Signal Processing，1988，36（7）：1169－1179.

［21］Daugman J G.Uncertainty Relation for Resolution in Space，Spatial Frequency and Orientation Optimized by Two－Dimensional Visual Conical Filters［J］.Journal of the Optical Society of America A，1985，2（7）：1160－1169.

［22］Daugman J G.How Iris Recognition Works［J］.IEEE Transaction on Circuits and Systems for Video Technology，2004，14（1）：21－30.

［23］Lee T S.Image Representation Using 2DGabor Wavelets［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1996，18（10）：959－971.

［24］Dunn D F，Higgjns W F.Optimai Gabor Filters for Texture Segmentation［J］.IEEE Transactions on Image Processing，1995，4（7）：947－964.

［25］Pearson E S，Biometrika Tables for Statisticians［M］.Cambridge：Cambridge University Press，1972.

［26］Hill I D，Hill R，Holder R L.Algorithm as 99：Fitting Johnson Curves by Moments［J］.Applied Statistics，1976，25（2）：180－189.