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

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

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

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

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

1 主軸-軸承系統(tǒng)模型

1.1 滾動軸承接觸力模型

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

fθj=kn(rθj)n

(1)

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

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

(2)

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

(3)

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

(4)

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

1.2 主軸-軸承系統(tǒng)動力學模型

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

表1 角接觸球軸承的原始參數(shù)

圖1 典型數(shù)控車床主軸-軸承系統(tǒng)Fig.1 A typical spindle-bearing system of the CNC lathe

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

圖2 機床主軸-軸承系統(tǒng)力學簡化模型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； 集中質量間軸段的剛度系數(shù)為k12=k23=k34=k45=1.41×108N/m，阻尼系數(shù)為比例阻尼，其中c1=αk1，c2=αk2，α=0.001/ω，e=10-5m。

2 Kriging模型與學習停止條件

2.1 Kriging模型

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

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

(6)

式中：fT(x)β為回歸模型，β為回歸系數(shù)向量，f(x)為隨機變量x的多項式函數(shù)，通?？梢匀」潭ㄖ?，其取值的大小并不影響模型的近似精度。

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

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

(7)

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

(8)

(9)

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

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

(10)

β和σ2的估計值為

(11)

(12)

式中：N0為向量的維數(shù)。

通過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 學習函數(shù)與停止準則

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

(15)

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

(16)

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

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

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

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

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

表2 主軸系統(tǒng)的設計變量

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

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

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

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

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

4 結果與分析

4.1 參數(shù)變化對主軸軸端振幅的影響

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

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

圖4 n=7 500 r/min時主軸前端軸心隨相關參數(shù)變化的分叉圖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模型，由學習函數(shù)進行樣本數(shù)據(jù)更新，當滿足停止學習條件時，得到最終的Kriging模型和系統(tǒng)的可靠度。為了驗證方法本文提出方法的有效性，應用Monte Carlo抽樣方法抽得5×104組樣本數(shù)據(jù)，并將其帶入振動微分方程式(5)，可分別求得不同樣本點的軸端軸心的最大振動幅值。假設軸心許用振動幅值為15 μm，此時，可得到了極限狀態(tài)函數(shù)的樣本歷史和概率直方圖，如圖5和圖6所示。分析圖6可知，極限狀態(tài)函數(shù)值g(x)的分布圖接近正態(tài)分布概率函數(shù)曲線，比較光滑且沒有間隙，表明抽樣次數(shù)足夠多。表3為本文提出方法和Monte Carlo法的可靠度計算結果。由于主軸系統(tǒng)是典型的強非線性系統(tǒng)，而很多可靠度計算方法在針對強非線性系統(tǒng)進行求解時會失效。通過對比Monte Carlo法和本文所使用的方法，說明了本文提出的方法在求解強非線性系統(tǒng)可靠度的有效性，也論證了該方法的高效性。

表3 不同方法的計算結果

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

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

5 結 論

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

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

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