基于AK-MCS法的主軸系統(tǒng)振動可靠性分析

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

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

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

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

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

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

1.1 滾動軸承接觸力模型

軸承的內(nèi)、外圈滾道是與滾動體相接觸的，在每一個接觸點，由于是Hertz接觸，滾動體與滾道之間的接觸變形會產(chǎn)生一個具有非線性特性的恢復(fù)力[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)中，由軸承滾動體接觸力的非負(fù)性可得

(3)

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

(4)

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

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

典型數(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)進(jìn)行了簡化，簡化系統(tǒng)由5個集中質(zhì)量單元和4個彈性軸段組成，考慮集中質(zhì)量單元在X和Y方向的平動，可將該主軸-滾動軸承系統(tǒng)簡化為十自由度系統(tǒng)，如圖2所示。其中，非線性軸承力和轉(zhuǎn)子的偏心力為系統(tǒng)的激勵。簡化后主軸系統(tǒng)的動力學(xué)微分方程為

圖2 機床主軸-軸承系統(tǒng)力學(xué)簡化模型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模型與學(xué)習(xí)停止條件

2.1 Kriging模型

Kriging模型[18]是一種半?yún)?shù)化的插值模型，不需要給出狀態(tài)函數(shù)的具體形式，這樣可以使模型的預(yù)測精度不受假定函數(shù)形式的影響。另外，Kriging模型可以應(yīng)用于強非線性的問題。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ù)θ的關(guān)于樣本點xi和xj的相關(guān)函數(shù)，模型的準(zhǔn)確性取決于隨機過程z(x)，相關(guān)函數(shù)通常選用高斯相關(guān)方程，其表達(dá)式為

(8)

(9)

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

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

(10)

β和σ2的估計值為

(11)

(12)

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

通過Kriging模型，得到待測點x的預(yù)測響應(yīng)值為

(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 學(xué)習(xí)函數(shù)與停止準(zhǔn)則

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

(15)

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

(16)

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)果與分析

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

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

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

圖4 n=7 500 r/min時主軸前端軸心隨相關(guān)參數(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模型，由學(xué)習(xí)函數(shù)進(jìn)行樣本數(shù)據(jù)更新，當(dāng)滿足停止學(xué)習(xí)條件時，得到最終的Kriging模型和系統(tǒng)的可靠度。為了驗證方法本文提出方法的有效性，應(yīng)用Monte Carlo抽樣方法抽得5×104組樣本數(shù)據(jù)，并將其帶入振動微分方程式(5)，可分別求得不同樣本點的軸端軸心的最大振動幅值。假設(shè)軸心許用振動幅值為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)進(jìn)行求解時會失效。通過對比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)進(jìn)行了簡化，建立了考慮主軸剛度、阻尼、偏心質(zhì)量和成組滾動軸承非線性接觸力的十自由度非線性振動模型。

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

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