基于Simulink仿真的凸輪軸高速磨削穩定性判定*

羅德龍， 鄧朝暉， 劉 濤， 佘帥龍， 羅程耀， 彭克立

(1. 湖南科技大學， 難加工材料高效精密加工湖南省重點實驗室， 湖南 湘潭 411201) (2. 湖南科技大學 智能制造研究院， 湖南 湘潭 411201) (3. 湖大海捷(湖南)工程技術研究有限公司， 長沙 410012)

凸輪軸作為配氣機構的關鍵零部件，其加工質量、精度等直接影響相關產品的使用性能和能耗標準[1]。近年來，高速磨削技術在凸輪軸加工中得到廣泛應用[2]。作為典型的非圓輪廓零件，凸輪軸在高速磨削過程中的磨削力、磨削接觸剛度等動力學參數均呈現出與普通外圓磨削和平面磨削不同的變化規律，給機床的動態性能帶來極大挑戰。

為提高凸輪軸磨削加工質量與精度，國內外學者在磨削顫振和工藝系統穩定性方面進行了諸多理論和實驗研究，顫振分析理論體系已基本建立[3-8]。HAHN[4]首先提出再生型顫振的概念，認為加工表面的波紋再生是引起工藝系統振動的主要原因。ALTINTAS等[5]對縱磨再生顫振動力學模型進行了建模分析及表面波紋度研究。ZHANG等[6]通過切除材料的幾何特征，建立了動態磨削力模型并分析了其與穩定性的關系。XU等[7]建立了外圓切入磨削非線性動力學模型。王龍山等[8]探討了影響工件顫振頻率的因素——工件接觸剛度和砂輪的非線性剛度，通過實驗和理論分析認為增大磨削力可以非線性提高其接觸剛度；而且，增大磨削力還可以使工件系統的顫振頻率提高。蔣永翔等[9]建立了外圓切入磨工件、砂輪再生顫振動力學模型及穩定性判定條件，并驗證了該模型及穩定性理論的正確性。遲玉倫等[10]對外圓切入磨材料去除率模型的時間常數進行研究，提出了磨削接觸剛度的有效測量方法，并建立其與系統固有頻率的關系。

上述研究都停留在理論層面，離實際工程應用相距甚遠，且多集中在平面或外圓磨削上[2]，而針對凸輪軸類非圓輪廓零件的磨削穩定性研究還不夠深入，缺乏精準的動力學參數獲取方法。因此，有必要針對凸輪軸高速磨削工藝系統穩定性問題開展研究，建立相應的磨削動力學模型，分析其穩定磨削區域，并基于Matlab/Simulink搭建磨削過程顫振仿真模型，選取適當的磨削工藝參數進行數值仿真分析與驗證，以達到避免顫振、提高加工質量和效率的目的。

1 凸輪軸高速磨削系統動力學模型建立

1.1 磨削系統簡化

一般認為磨床在切線方向(Y向)的整機動剛度足夠大，引起動態切削厚度b(t)變化的主要原因是磨削力徑向分力的波動。因此，磨削系統可以簡化為如圖1所示的單自由度系統。

圖1 凸輪軸高速磨削系統動力學模型 Fig. 1 Dynamic model of camshaft high speed grinding system

該磨削系統動力學模型的運動微分方程可表示為：

(1)

其中：Fn(t)為磨削過程法向動態磨削力，m為系統的等效質量，K為系統的等效剛度，C為系統的等效阻尼，并設X(t)為砂輪架系統的瞬時振動位移。

1.2 磨削過程參數求解

根據熊超等[11]建立的凸輪非圓輪廓法向磨削力Fn數學模型：

(2)

其中：k是與工件材料有關的磨削力系數；μ表示砂輪寬度或者被加工工件寬度；γ為砂輪磨粒的圓錐半頂角，一般取γ=60°；d0為平均磨粒間隔；ε為影響系數，且ε=0.2～0.5；ap為砂輪磨削深度；vw、vs分別表示磨削點線速度和砂輪線速度；ds為砂輪直徑，dρ為凸輪上瞬時磨削點的曲率直徑。

由式(2)可以看出：影響磨削過程中法向動態磨削力變化的主要因素有瞬時磨削深度ap、凸輪磨削點的線速度vw和磨削點曲率半徑dρ等。由于本文凸輪軸磨削采用恒線速磨削加工方式，故動態磨削力的影響因素主要考慮瞬時磨削深度ap和當量砂輪直徑de，當量砂輪直徑的計算公式如下：

(3)

以某型號汽車凸輪軸為例，該凸輪基圓半徑15.5 mm，砂輪半徑200 mm，滾子測頭半徑8 mm，最大升程值為凸輪轉角位于98°時的5.1525 mm，且凸輪兩側帶有凹面。通過Matlab編程計算，得到凸輪曲率半徑dρ與轉角之間的對應關系，如圖2所示。

圖2 凸輪曲率半徑與轉角之間的對應關系 Fig. 2 The corresponding relationship between the Radius

根據文獻[1]，可以求得任意磨削點p處工件的角速度ω，如公式(4)所示：

(4)

同時，根據轉角與時間的函數關系式(5)，將凸輪轉角視為中間變量，即可得到凸輪曲率半徑與時間的對應關系，圖3(凸輪基圓轉速為nw=90 r/min)

(5)

圖3 曲率半徑與時間的對應關系圖 Fig. 3 The corresponding relationship between the radius of curvature and the time

根據文獻[1]，對于任意時刻t，磨削系統的輸入量為理論磨削深度；磨削系統的變量為砂輪磨損量、凸輪軸的去除量、砂輪與工件之間的接觸變形以及磨床系統的結構變形(即系統振動的振幅)；磨削系統在任意時刻的理論磨削深度等于此時的砂輪磨損量和凸輪軸的去除量、砂輪與工件之間的接觸變形、磨床系統的結構變形之和；磨削系統的輸出量為動態磨削力。由于在凸輪軸高速磨削過程中，砂輪存在定期修整，如果忽略砂輪的磨損對動態磨削過程的影響，從而得到任意時刻t凸輪磨削的實際磨削深度(凸輪軸的去除量)計算公式：

(6)

其中：a0表示凸輪磨削時單圈的理論磨削深度；ai(t)表示第i圈的實際瞬時磨削深度值；x(t)表示任意時刻t砂輪與工件之間的接觸變形量；b(t)表示工件表面的動態振紋。b(t)為：

b(t)=b0-[X(t)-X(t-T)]

(7)

其中：b0表示凸輪表面的初始振紋，X(t)為工藝系統的瞬時振動位移，T為工件主軸的旋轉周期。

從式(7)可以看出：凸輪軸高速磨削工藝系統是一個時滯非線振動系統。對于凸輪軸和砂輪之間的接觸變形x，考慮工件經過了熱處理以及局部冷激工藝，具有較高的硬度，而砂輪表面的磨粒則是結合在結合劑上的，其接觸變形較凸輪的接觸變形要大得多，所以在分析中可用砂輪的接觸變形近似代替砂輪和工件之間的總變形。砂輪的接觸變形函數表達式如下[8]：

(8)

其中：k1=2.23，e=0.706，x的單位是μm。

2 凸輪軸高速磨削過程仿真與結果分析

2.1 磨削顫振仿真模型

根據前文所建立的凸輪軸高速磨削系統動力學模型，基于Matlab/Simulink搭建磨削過程顫振仿真模型，考慮將整個仿真分為2個功能模塊，一塊主要是用于求解動力學模型所對應的振動微分方程，另一塊則是用于求解動態磨削力的子模塊。

對于微分方程求解，通常采用積分方式進行子結構建模，首先選擇Continuous中的integrator模塊表示微分方程的各階積分；采用Commonly used Blocks中的Gain模塊表示微分方程的各個系數；采用Math Operations中的Add運算模塊將微分方程進行等式表達；采用Sinks中的Scope模塊顯示輸出。

而對于動態磨削力的求解，由以上分析可知，主要求解瞬時磨削深度和當量砂輪直徑的值。由瞬時磨削深度推導公式(6)，采用Sources中的Clock模塊來表示模型仿真時間；采用Signal Routing中的Switch模塊表示仿真的分支結構；采用Commonly used Blocks中的Constant輸入模塊表示凸輪軸磨削加工設定的理論磨削深度；采用Continuous中的Transport Delay模塊表示凸輪軸磨削加工每一圈的時滯項；采用Math Operations中的Add運算模塊對磨削深度值計算求解；采用User-Defined Functions中的Fcn模塊表示磨削力與磨削深度之間的函數關系；同時，根據前文所述的凸輪輪廓曲率半徑以及當量砂輪直徑的計算方法，為了將Matlab工作區中計算得到的當量砂輪直徑隨時間的函數關系代入到Simulink仿真模型中，采用Sources中的From Workspace模塊來讀入工作區中的計算數據，并使用Math Operations中的Product模塊表示動態磨削力的各個影響因素之間的乘積關系。基于Matlab/Simulink環境，搭建的磨削過程計算機仿真系統模型如圖4所示。

圖4 凸輪軸高速磨削加工系統仿真框圖 Fig. 4 Simulation diagram of camshaft high speed grinding

2.2 凸輪軸高速磨削過程仿真結果分析

根據前期對CNC8325B數控凸輪軸磨床進行的整機動剛度測試結果[12]：磨床整機X向(法向)動剛度為K=4.834×106N/m，系統的固有頻率ωn=187 Hz，并由頻率響應曲線用半功率法求得系統的阻尼比ξ=1.2%，由此可以求得磨削系統的等效質量、阻尼參數等。同時，磨削過程的其他仿真參數設置如下：凸輪寬度μ=25 mm，平均磨粒間隔d0=0.67 mm，影響系數ε=0.28，與工件有關的磨削力系數k=170 kg/mm，圓錐半頂角γ=60°。

仿真時，凸輪表面的初始振紋b0忽略不計，磨削過程的初始條件為X(0)=0、X′(0)=0。當砂輪線速度vs=120 m/s，工件基圓轉速nw=90 r/min，理論磨削深度a0=0.02 mm時，凸輪軸高速磨削過程振動位移時域圖如圖5所示。

圖5 凸輪軸高速磨削過程振動位移時域圖 Fig. 5 The Time domain diagram of vibration displacement of camshaft high speed grinding

由圖5可知：在上述磨削工藝參數下，工藝系統的顫振振幅最終趨于穩定，能夠實現凸輪的穩定磨削。工件基圓轉速nw=90 r/min時，通過Matlab求解得到單圈的磨削時間T=1.06 s。當凸輪從基圓上與桃尖相對的位置開始磨削前半圈，即圖6所示凸輪輪廓的0°～180°范圍時，在磨削開始階段，由于磨削系統狀態發生改變，產生劇烈的受迫振動，隨后在系統阻尼的作用下振紋逐漸減小；而當磨削進行到凸輪升程段時，系統的振幅又有所增大，之后又繼續減小，對應圖5所示磨削第一圈的前半程。凸輪輪廓示意圖如圖6所示。

圖6 凸輪輪廓

這一現象產生的原因可以從以下2個方面來分析：(1)由于凸輪升程兩側輪廓存在凹面，在磨削過程中引起動態磨削力的突變，因而振幅出現突變的現象；(2)在磨削凸輪升程段時，凸輪輪廓的曲率半徑一般是先增大后減小，而由動態磨削力的計算公式(2)可知，動態磨削力與輪廓曲率半徑之間存在正相關關系，因而系統的振幅也是先增大后減小。而當磨削過程從桃尖回到凸輪基圓的過程中，即圖6所示的180°～360°范圍內，系統的顫振振幅變化趨勢與磨削凸輪前半圈時相似，都是先增大后減小，并且在回程段與基圓的連接處，由于凸輪輪廓凹面的影響，振幅同樣存在突變現象，對應圖5所示磨削第一圈的后半程。這與前期對凸輪軸磨削力的研究結果相一致。根據文獻[13]，凸輪軸高速磨削動態磨削力的單圈采集實驗結果如圖7所示，其橫坐標為凸輪輪廓上從0°～360°之間所布置的測點數，縱坐標為相應測點處所對應的磨削力。由圖7可以看出磨削力隨凸輪輪廓測點(也即轉角)的關系基本上滿足先增大后減小，再增大再減小的變化規律。

當凸輪磨削進行到第二圈時，由于前一圈磨削時工藝系統的振動以及砂輪的接觸變形等使得磨削過程的實際磨削深度值小于理論磨削深度，從而在工件表面留下了動態振紋，根據瞬時磨削深度計算公式(6)及磨削力公式(2)可知，此時的動態磨削力發生突變，從而使得磨削系統的振幅急劇增大，對應圖5中第一圈與第二圈之間交界部分系統振幅的突變現象。從圖5中也能看出：其后的磨削過程系統的振幅變化趨勢基本與第一圈類似，凸輪磨削各圈的振幅變化規律基本都是先增大后減小，再增大再減小。

圖7 法向磨削力單圈采集時域圖 Fig. 7 Time domain diagram of single ring acquisition of grinding force

3 動態磨削仿真判定與葉瓣圖判定

根據前期對凸輪軸高速磨削過程穩定性理論及研究結果分析[1]，得到磨削系統的穩定性葉瓣圖如圖8。

由圖8可知：穩定性葉瓣圖由一系列穩定性葉瓣構成，位于葉瓣下方的區域稱為穩定磨削區，即圖8中的藍色部分；位于葉瓣上方的區域稱為不穩定磨削區；而當磨削工藝參數剛好落在葉瓣線上時，此時的工藝系統處于臨界穩定磨削狀態，即系統本身是穩定的，但是如果給系統施加一個外部干擾，磨削系統則可能會從穩定狀態過渡到磨削顫振狀態。

圖8 工藝系統穩定性葉瓣圖 Fig. 8 The diagram of stability leaf flap of process system

在穩定性葉瓣圖上選取幾個比較有代表性的點，用這些點的磨削工藝參數組合對模型進行仿真驗證。本研究選用的工藝參數組合如表1所示，表1中的實驗號對應圖8中相應點的編號，各磨削工藝參數對應的仿真結果分別如圖9a～圖9f所示。

由圖9a、圖9b可知：參數組合3、4下發生了磨削顫振，振幅值不斷增大，但其振幅增長比較緩慢，對應葉瓣圖中穩定磨削區與不穩定磨削區的臨界區域；由圖9c、圖9d可知：參數組合5、6下振動逐漸趨于穩定，對應葉瓣圖中的穩定磨削區域；由圖9e、圖9f可知：參數組合11、12下發生了磨削顫振，振動幅值不斷增大，對應葉瓣圖中的不穩定磨削區域。由此可知，基于動態磨削仿真的磨削穩定性判定與基于穩定性葉瓣圖的磨削穩定性判定基本吻合，驗證了本研究所建立的動力學模型及仿真方法的正確性。

表1 凸輪軸高速磨削仿真工藝參數表Tab. 1 The parameter table of camshaft high speed grinding simulation

(a)3號(b) 4號(c)5號(d)6號(e)11號(f)12號圖9 不同參數下的仿真結果圖Fig. 9 The simulation result graph under different parameters

4 結論

以凸輪軸數控高速磨削過程為研究對象，分析磨削過程中的動態影響因素，建立磨削系統動力學模型，通過 Matlab/Simulink環境搭建了磨削過程的顫振仿真框圖，基于凸輪軸高速磨削的穩定性分析結果，從穩定性葉瓣圖上選取有代表性的點進行數值仿真分析，通過仿真所得振動位移時域圖判定不同工藝參數組合下工藝系統的穩定性，并將仿真結果與穩定性葉瓣圖所示的穩定磨削區域進行比較，驗證了仿真模型的有效性，為凸輪軸的實際磨削加工工藝參數選擇提供參考。