基于AK-MCS法的主軸系統振動可靠性分析

馮吉路， 孫志禮， 趙 堅， 張 婧， 梁春芳

(1. 天津城建大學 控制與機械工程學院，天津 300384；2. 東北大學 機械工程與自動化學院，沈陽 110819)

隨著數控技術的不斷發展，高速加工對數控機床主軸系統轉速、回轉精度和可靠性的要求日益提高，數控機床主軸系統的動力學特性已然成為了當前研究的熱點問題[1]。近年來，研究人員采用傳遞矩陣法[2]、有限單元法[3-7]和集中參數法[8-10]等主軸動力學建模方法對主軸系統動態特性進行了大量的研究，主要分析了相關因素對主軸固有頻率和振幅的影響。主軸軸承滾子與滾道之間非連續、非光滑接觸會使得主軸系統產生非線性振動，進而影響主軸系統的動力學特性。特別在軸承滾子與滾道之間留有間隙時，該間隙會在軸承非線性Hertz接觸力和主軸非平衡力的作用下，引起主軸系統的非線性分岔和混沌現象，進而可能導致主軸系統振幅過大或系統失穩[11-12]。

機床主軸系統動力學性能可靠性是設計工作者應該重點關注的問題，它是衡量機床整機性能的主要指標，會受到軸承和主軸材料的剛度和阻尼等諸多參數的影響。在這些設計參數變化的情況下，主軸軸端軸心的振動幅值必然會發生變化。當振動幅值超過設計指標時，就會引起設計失效，而不發生設計失效的概率，可以被定義為主軸振動幅值的可靠度。由于主軸動力學模型是強非線性系統，采用一般的可靠性分析方法必然會增加計算工作量，延長設計周期。為了解決此問題，近年來研究人員提出了各種方差縮減技術，如分層抽樣[13]、重要抽樣[14]、控制變量[15]等方法。雖然上述方法能夠減少樣本數量，但是很難滿足強非線性系統的實際工程需求。

鑒于上述問題，建立了考慮彈性主軸剛度、阻尼和滾動軸承非線性接觸力的十自由度數控車床主軸非線性振動分析模型。利用Runge-Kutta數值積分法對數控車床主軸-滾動軸承系統的振動微分方程進行了數值求解。采用Kriging模型和Monte Carlo法相結合的方法計算了主軸動力學性能可靠度。該研究為設計工作者進行機床主軸系統動力學性能可靠性設計提供了理論基礎。

1 主軸-軸承系統模型

1.1 滾動軸承接觸力模型

軸承的內、外圈滾道是與滾動體相接觸的，在每一個接觸點，由于是Hertz接觸，滾動體與滾道之間的接觸變形會產生一個具有非線性特性的恢復力[16]

fθj=kn(rθj)n

(1)

式中：n為常數，對于球軸承，n=3/2；rθj為滾動體j在角位置為θj處的接觸變形量；kn為接觸剛度。在任意時刻，第j個滾動體在接觸點上的彈性變形取決于滾動軸承質心的位移和軸承的初始游隙

rθj=xsinθj+ycosθj-γ

(2)

式中：γ為軸承的初始游隙，將式(2)代入式(1)中，由軸承滾動體接觸力的非負性可得

(3)

總的接觸力是每一個滾動體接觸力的總和，考慮到接觸阻尼cn的影響，可以得到在x和y方向上滾動軸承軸承總接觸力

(4)

由式(4)可知，軸承的接觸力與主軸的軸心軌跡有緊密的聯系，并且滾動軸承的剛度具有時變性，可能導致主軸-滾動軸承系統運轉過程中出現混沌響應。根據文獻[17]計算了機床主軸系統中的前后軸承組的等效總接觸剛度系數分別為kn1=3.145×109N/m3/2，kn1=2.178×109N/m3/2，軸承組的總接觸阻尼大小和剛度系數成比例，為cn=0.25×10-5kn；對于軸承游隙，假設成組軸承中每個軸承的游隙完全相同，前后軸承組的游隙分別為r1=53.65 μm，r2=19 μm。

1.2 主軸-軸承系統動力學模型

典型數控車床主軸系統主要由主軸、滾動軸承、卡盤、帶輪以及定位元件等組成，其中主軸系統前、后端軸承組分別采用NSK 7016A5(DBD) 組合和 7015C (DB)組合的形式，如圖1所示。NSK 7016A5和7015C的軸承結構參數如表1所示。

表1 角接觸球軸承的原始參數

圖1 典型數控車床主軸-軸承系統Fig.1 A typical spindle-bearing system of the CNC lathe

分析數控車床主軸系統零部件分布可知，該系統的質量主要位于主軸的兩端以及前、后軸承組之間。為了便于分析，研究中采用張偉剛等所使用的集中質量法對國產數控車床主軸系統進行了簡化，簡化系統由5個集中質量單元和4個彈性軸段組成，考慮集中質量單元在X和Y方向的平動，可將該主軸-滾動軸承系統簡化為十自由度系統，如圖2所示。其中，非線性軸承力和轉子的偏心力為系統的激勵。簡化后主軸系統的動力學微分方程為

圖2 機床主軸-軸承系統力學簡化模型Fig.2 Simplified mechanical model for machine tool spindle-bearing system

(5)

式中： 主軸簡化后的集中質量分別為m1=9.35 kg，m2=4.61 kg，m3=8.59 kg，m4=1.77 kg，m5=13.76 kg； 集中質量間軸段的剛度系數為k12=k23=k34=k45=1.41×108N/m，阻尼系數為比例阻尼，其中c1=αk1，c2=αk2，α=0.001/ω，e=10-5m。

2 Kriging模型與學習停止條件

2.1 Kriging模型

Kriging模型[18]是一種半參數化的插值模型，不需要給出狀態函數的具體形式，這樣可以使模型的預測精度不受假定函數形式的影響。另外，Kriging模型可以應用于強非線性的問題。Kriging模型表示為

g(x)=fT(x)β+z(x)

(6)

式中：fT(x)β為回歸模型，β為回歸系數向量，f(x)為隨機變量x的多項式函數，通常可以取固定值，其取值的大小并不影響模型的近似精度。

z(x)是隨機過程函數，反映局部偏差的近似，它的均值μ是零，方差是σ2，協方差表示為

cov(z(xi),z(xj))=σ2R(xi,xj)

(7)

式中：R(xi,xj)為帶有參數θ的關于樣本點xi和xj的相關函數，模型的準確性取決于隨機過程z(x)，相關函數通常選用高斯相關方程，其表達式為

(8)

(9)

給定訓練樣本集合S={x1，x2，…，xns}，計算相應的實際功能函數響應值，將其用向量形式表達為

Y=[g(x1),g(x2),…,g(xns)]T

(10)

β和σ2的估計值為

(11)

(12)

式中：N0為向量的維數。

通過Kriging模型，得到待測點x的預測響應值為

(13)

式中：f=[f(x1),f(x2),…,f(xns)]T；r(x)=[R(x,x1),R(x,x2),…,R(x,xns)]T。

Kriging方差為

(14)

式中：u(x)=fT(x)R-1r(x)-f。

2.2 學習函數與停止準則

為了使Kriging模型不斷優化，實現自動查找和增加最合適的訓練樣本，定義學習函數U(x)為

(15)

一般而言，可以選擇U(x)取值最小的樣本點作為最合適的訓練樣本點。為滿足進一步分析的需求，還需要定義停止學習的條件。由于學習函數U(x)與Cox和John提出的用于優化的具有比較小置信區間的函數具有一定聯系，根據相關研究，本文給出了相應的停止學習準則

(16)

式中：A為將要分類的樣本集。在判斷這些樣本點所在區域時發生錯誤的概率為Ф(-2)=0.023。通過建立學習函數和相應的學習停止準則，可以保證在不斷更新后的Kriging模型具有更好的精度。

3 主軸振動可靠性的計算流程

根據建立的十自由度振動微分方程求得了主軸轉速為n=0～4 000 r/min時系統軸端軸心的最大振動幅值，并使用Kriging模型構建了樣本數據與響應值之間的關系。采用AK-MCS主動學習方法進行主軸軸端振動可靠度計算，計算流程如圖3所示。具體計算步驟如下：

圖3 AK-MCS算法流程圖Fig.3 Flow chart of the AK-MCS

步驟1由本文“4.1”節可知，主軸系統軸端振動幅值的大小會受到主軸前、后軸承組的軸承接觸剛度和游隙的因素的影響。假設設計變量服從正態分布，根據軸承廠提供的軸承設計參數，結合主軸系統中軸承組受載情況，應用赫茲接觸理論確定的設計變量的均值。各變量的均值和標準差如表2所示，其中x1為主軸系統前軸承組的等效接觸剛度、x2為后軸承組的等效接觸剛度、x3為前軸承組單個軸承的徑向間隙、x4為后軸承組單個軸承的徑向間隙，其余參數均為定值。根據Monte Carlo方法，在隨機變量的設計空間中生成的n組隨機樣本，稱為候選樣本點，用S來表示，這些候選樣本點不需計算其實際功能函數響應值，它們會被應用到之后構建Kriging模型時利用學習函數自動查找和增加最佳候選樣本點的過程當中，只有被選為最佳候選樣本點時，才需要計算該最佳候選樣本點的實際功能函數響應值。

表2 主軸系統的設計變量

步驟2采用拉丁超立法抽樣方法生成初始訓練樣本點，用S1來表示這些樣本點需要通過仿真或計算得到其實際功能函數響應值，并用Y1來表示。初始訓練樣本點中樣本點的數量一般選取的較少，會在之后的學習過程中不斷地增加最佳候選樣本點，不斷地優化模型。

步驟3利用訓練樣本點S1和實際功能函數響應值Y1， 應用MATLAB中的DACE工具箱，構建Kriging模型，Kriging模型的回歸模型部分選用一次多項式，相關函數采用高斯相關函數。

步驟4利用構建的Kriging模型，計算候選樣本點S中所有樣本點的預測值和預測方差，并計算這些候選樣本點的學習函數值，并挑選候選樣本點中學習函數值中最小的樣本點，作為最佳候選點。

步驟5如果最佳候選點的學習函數值滿足學習停止條件，則構建Kriging模型的主動學習過程結束，否則，將最佳候選點加入到訓練樣本點中，并計算其實際功能函數響應值，轉到步驟3，構建新的Kriging模型。

步驟6如果主動學習過程結束，采用當前的訓練樣本點及其實際功能函數響應值所構建的Kriging模型，此時的Kriging模型為優化后的比較精確的模型，運用Monte Carlo方法計算失效概率Pf以及可靠度Pr。

4 結果與分析

4.1 參數變化對主軸軸端振幅的影響

數控車床主軸系統是一個具有強非線性特征的滾動軸承-主軸系統。主軸系統的軸端動態響應會決定機床工作過程中的加工精度，而軸端的動態響應主要會受到前、后軸承組的接觸剛度和軸承游隙等相關參數的影響。

圖4為n=6 500 r/min時主軸前端軸心隨前、后軸承組的總接觸剛度以及前、后軸承游隙變化的穩態響應分叉圖。由圖分析可知，軸承接觸剛度和軸承游隙對主軸軸心振幅會產生影響；當其它參數保持不變時，前軸承組總接觸剛度僅僅在[1.15×109N/m1.5,1.2×109N/m1.5]內變化，會導致主軸軸端響應出現跳躍現象，說明此時系統出現了多個周期吸引子；相比之下，后軸承組的接觸剛度較低時，很容易引起系統的擬周期或者混沌運動狀態，但是當后軸承剛度增大到一定程度時，剛度幾乎不會引起主軸前端軸心振動幅值的變化。由圖4(c)和圖4(d)可知，當其它參數不發生變化時，前、后軸承游隙增加到一定程度時，將不會對主軸軸端的振動幅值產生影響。

圖4 n=7 500 r/min時主軸前端軸心隨相關參數變化的分叉圖Fig.4 Bifurcation diagram for the axis at the front of the spindle with changed parameters when rotational speed n=7 500 r/min

4.2 主軸軸端振動可靠性計算結果

采用拉丁超立方抽樣抽取60個樣本點，建立初始Kriging模型，由學習函數進行樣本數據更新，當滿足停止學習條件時，得到最終的Kriging模型和系統的可靠度。為了驗證方法本文提出方法的有效性，應用Monte Carlo抽樣方法抽得5×104組樣本數據，并將其帶入振動微分方程式(5)，可分別求得不同樣本點的軸端軸心的最大振動幅值。假設軸心許用振動幅值為15 μm，此時，可得到了極限狀態函數的樣本歷史和概率直方圖，如圖5和圖6所示。分析圖6可知，極限狀態函數值g(x)的分布圖接近正態分布概率函數曲線，比較光滑且沒有間隙，表明抽樣次數足夠多。表3為本文提出方法和Monte Carlo法的可靠度計算結果。由于主軸系統是典型的強非線性系統，而很多可靠度計算方法在針對強非線性系統進行求解時會失效。通過對比Monte Carlo法和本文所使用的方法，說明了本文提出的方法在求解強非線性系統可靠度的有效性，也論證了該方法的高效性。

表3 不同方法的計算結果

圖5 g(x)的樣本歷史Fig.5 Sample history of g(x)

圖6 極限狀態函數g(x)頻率直方圖Fig.6 Frequency histogram of limit state function g(x)

5 結 論

(1) 采用集中質量法對國產數控車床主軸系統進行了簡化，建立了考慮主軸剛度、阻尼、偏心質量和成組滾動軸承非線性接觸力的十自由度非線性振動模型。

(2) 主軸系統后軸承組接觸剛度較低時，很容易引起系統的擬周期或者混沌運動狀態，且前后軸承組軸承游隙匹配合適有利于減小主軸振幅。

(3) 通過對比Monte Carlo和AK-MCS法計算的可靠度結果可知，AK-MCS法計算效率高，收斂速度快，特別適用于解決強非線性隱式復雜問題，且有效提高數控機床主軸可靠性設計的效率。