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

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

隨著數(shù)控技術的不斷發(fā)展，高速加工對數(shù)控機床主軸系統(tǒng)轉(zhuǎn)速、回轉(zhuǎn)精度和可靠性的要求日益提高，數(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法相結(jié)合的方法計算了主軸動力學性能可靠度。該研究為設計工作者進行機床主軸系統(tǒng)動力學性能可靠性設計提供了理論基礎。

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

1.1 滾動軸承接觸力模型

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

fθj=kn(rθj)n

(1)

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

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

(2)

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

(3)

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

(4)

由式(4)可知，軸承的接觸力與主軸的軸心軌跡有緊密的聯(lián)系，并且滾動軸承的剛度具有時變性，可能導致主軸-滾動軸承系統(tǒng)運轉(zhuǎn)過程中出現(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的軸承結(jié)構(gòu)參數(shù)如表1所示。

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

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

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

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

(5)

式中： 主軸簡化后的集中質(zhì)量分別為m1=9.35 kg，m2=4.61 kg，m3=8.59 kg，m4=1.77 kg，m5=13.76 kg； 集中質(zhì)量間軸段的剛度系數(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ù)建立的十自由度振動微分方程求得了主軸轉(zhuǎn)速為n=0～4 000 r/min時系統(tǒng)軸端軸心的最大振動幅值，并使用Kriging模型構(gòu)建了樣本數(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ù)，結(jié)合主軸系統(tǒng)中軸承組受載情況，應用赫茲接觸理論確定的設計變量的均值。各變量的均值和標準差如表2所示，其中x1為主軸系統(tǒng)前軸承組的等效接觸剛度、x2為后軸承組的等效接觸剛度、x3為前軸承組單個軸承的徑向間隙、x4為后軸承組單個軸承的徑向間隙，其余參數(shù)均為定值。根據(jù)Monte Carlo方法，在隨機變量的設計空間中生成的n組隨機樣本，稱為候選樣本點，用S來表示，這些候選樣本點不需計算其實際功能函數(shù)響應值，它們會被應用到之后構(gòu)建Kriging模型時利用學習函數(shù)自動查找和增加最佳候選樣本點的過程當中，只有被選為最佳候選樣本點時，才需要計算該最佳候選樣本點的實際功能函數(shù)響應值。

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

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

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

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

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

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

4 結(jié)果與分析

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]內(nèi)變化，會導致主軸軸端響應出現(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 主軸軸端振動可靠性計算結(jié)果

采用拉丁超立方抽樣抽取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法的可靠度計算結(jié)果。由于主軸系統(tǒng)是典型的強非線性系統(tǒng)，而很多可靠度計算方法在針對強非線性系統(tǒng)進行求解時會失效。通過對比Monte Carlo法和本文所使用的方法，說明了本文提出的方法在求解強非線性系統(tǒng)可靠度的有效性，也論證了該方法的高效性。

表3 不同方法的計算結(jié)果

圖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 結(jié) 論

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

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

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