王 民，王樹華，孫鐵偉，石 龍

(1.北京工業(yè)大學(xué)機(jī)械工程與應(yīng)用電子技術(shù)學(xué)院，先進(jìn)制造技術(shù)北京市重點實驗室，北京 100124；2.電火花加工技術(shù)北京市重點實驗室，北京 100191)

現(xiàn)代制造業(yè)水平的飛速提升導(dǎo)致了人們對數(shù)控機(jī)床整體性能的要求越來越高，滾珠絲杠副是數(shù)控機(jī)床的重要傳動部件，與數(shù)控機(jī)床整體的可靠性緊密關(guān)聯(lián)，在數(shù)控機(jī)床、自動化裝備、車輛交通及航空航天等領(lǐng)域常常用到[1].因此，對滾珠絲杠副的可靠性分析工作得到越來越多的科研機(jī)構(gòu)和研究人員的重視.

近年來，隨著可靠性工程的不斷發(fā)展，相關(guān)的數(shù)據(jù)樣本量較小.對于滾珠絲杠副產(chǎn)品，由于價格昂貴、可靠性水平較高，而且試驗周期長，很難進(jìn)行大樣本的統(tǒng)計試驗，甚至有的專項計劃只生產(chǎn)1至2臺產(chǎn)品.綜合考慮上述情況，本文通過引入小子樣、加速壽命試驗的方法完成對單樣本滾珠絲杠副的壽命評估.

在實際工況下，數(shù)控機(jī)床進(jìn)給系統(tǒng)的滾珠絲杠在加減速運(yùn)行過程中，潤滑會變差，加劇接觸部件的磨損，降低其精度和壽命[2]，進(jìn)而導(dǎo)致數(shù)控機(jī)床的可靠性水平下降、加工質(zhì)量變差.另外，根據(jù)Archard[3]磨損理論知道外界負(fù)載應(yīng)力是影響磨損量的主要因素之一.因此本文將綜合考慮轉(zhuǎn)速和負(fù)載雙應(yīng)力對滾珠絲杠副可靠性壽命的影響進(jìn)行加速壽命試驗，然后通過對試驗數(shù)據(jù)進(jìn)行統(tǒng)計分析得到在加速應(yīng)力條件下的壽命分布模型，最后利用加速方程外推出在正常應(yīng)力水平下的壽命特征，完成對滾珠絲杠副的壽命評估[4-5].

為了提高小子樣單樣本滾珠絲杠副的可靠性評估精度，本文使用虛擬增廣法對加速壽命試驗得到的故障數(shù)據(jù)進(jìn)行等價擴(kuò)充，根據(jù)虛擬增廣法擴(kuò)充得到的數(shù)組和原試驗數(shù)據(jù)的分布形式相同.然后對增廣后的數(shù)組進(jìn)行統(tǒng)計分析，從而達(dá)到提高試驗評估精度的目的.

常用的可靠性壽命及故障分布模型有正態(tài)分布、指數(shù)分布、威布爾(Weibull)分布[6].其中Weibull分布是一類廣泛應(yīng)用于機(jī)械、電子、航天航空等領(lǐng)域的壽命分布模型，是目前較為常用的壽命模型之一.Weibull分布有多種參數(shù)估計方法[7]：其中圖估計法易受人為因素影響，精度偏低；極大似然估計法求解煩瑣，常常需要借助計算機(jī)；線性無偏估計法不適用于小子樣情況，在大樣本條件下還需借助查閱專用表格，雖然步驟簡單，但表格量大且在無表格處無法使用此方法；最小二乘法求解簡便，具有唯一最優(yōu)解和好的解析性質(zhì)，是Weibull分布參數(shù)估計的一種常用方法.

綜合以上特點，本文針對小子樣單樣本滾珠絲杠副提出一種雙應(yīng)力三水平的加速壽命試驗設(shè)計和統(tǒng)計分析方法，利用虛擬增廣法和最小二乘法完成在3個步進(jìn)加速雙應(yīng)力水平下的參數(shù)估計和可靠性評估，然后聯(lián)立求解出加速方程，進(jìn)而計算出滾珠絲杠副在正常應(yīng)力條件下的壽命指標(biāo)，完成對單樣本滾珠絲杠副的可靠度壽命分析.

1 加速壽命試驗?zāi)Ｐ?/h2>

1.1 基本假定

滾珠絲杠副具有加速性的3個條件為：失效分布同一性、失效機(jī)理一致性和失效過程規(guī)律性[8].根據(jù)以上具有加速性的3個條件，本文提出3個基本假定.

假定1：樣本在任意應(yīng)力下的壽命T均服從Weibull分布，即T～W(m,η)，累積分布函數(shù)為

(1)

式中：η>0為尺度參數(shù)，表示特征壽命；m>0為形狀參數(shù).

假定2：樣本在任意應(yīng)力下的失效機(jī)理不變，由于Weibull分布中的形狀參數(shù)反映失效機(jī)理，因此有

m1=m2=…≡m

(2)

假定3：樣本特征壽命與試驗雙應(yīng)力之間的關(guān)系模型選用Nelson廣義對數(shù)線性模型，其數(shù)學(xué)表達(dá)式為

lnη=β0+β1φ1(V)+β2φ2(P)

(3)

式中：η為特征壽命；φ1、φ2常為應(yīng)力的已知函數(shù)；β0、β1、β2為待估參數(shù)；V、P分別為轉(zhuǎn)速和負(fù)載.

考慮到本文的2個應(yīng)力均為非熱應(yīng)力，而且每個應(yīng)力的加速效果均呈逆冪律形式，因此應(yīng)力Vi、Pi與相對應(yīng)的特征壽命ηi可以通過式(3)轉(zhuǎn)化為雙應(yīng)力加速的逆冪律模型[9]：

(4)

對式(4)兩邊取對數(shù)，得到加速方程的對數(shù)形式為

lnηi=lnA-αlnVi-βlnPi

(5)

式中A、α、β均為常數(shù).

1.2 加速系數(shù)

加速系數(shù)是指滾珠絲杠副在正常應(yīng)力組合(V0,P0)與加速應(yīng)力組合(Vi,Pi)下達(dá)到特征壽命需要時間T0和Ti的比值，記為τ.即

(6)

由式(4)(6)可得加速系數(shù)理論公式為

(7)

1.3 可靠性指標(biāo)

滾珠絲杠副是可修復(fù)、可更換系統(tǒng)，出現(xiàn)故障后，可對零部件進(jìn)行修復(fù)或更換，并且記錄故障時間和修復(fù)時間，然后繼續(xù)進(jìn)行試驗.

將平均故障間隔時間作為可靠性定量指標(biāo).其表達(dá)式為

(8)

式中Γ為伽馬函數(shù)，其值可以通過查閱伽馬函數(shù)標(biāo)準(zhǔn)表[10]得到.

在特定工作條件下，將滾珠絲杠副零部件出現(xiàn)故障的概率稱為壽命可靠性.由假定1可知滾珠絲杠副可靠度的數(shù)學(xué)表達(dá)式為

(9)

2 加速壽命試驗

2.1 加速壽命試驗方案設(shè)計

試驗選用國產(chǎn)雙螺母預(yù)緊式滾珠絲杠副為研究對象，其結(jié)構(gòu)參數(shù)見表1[11].本文綜合考慮轉(zhuǎn)速以及軸向負(fù)載作為加速壽命試驗的2個應(yīng)力，分別將其記為V和P.額定轉(zhuǎn)速V0和額定負(fù)載P0分別為500 r/min、1 000 N.如果每個應(yīng)力選定3個加速水平，則共有3×3=9個應(yīng)力組合，如果將這9個應(yīng)力組合全部進(jìn)行試驗，則會大大延長試驗周期，降低試驗效率.因此按照低、中、高3個水平對加速應(yīng)力進(jìn)行合理等級劃分，根據(jù)對角線原則分別選取(V1,P1)、(V2,P2)、(V3,P3)作為3個步進(jìn)加速雙應(yīng)力進(jìn)行定時截尾試驗.由于試驗過程中滾珠絲杠副的失效機(jī)理不能改變，綜合考慮試驗方案和加載能力，結(jié)合工程經(jīng)驗，本文選取軸向負(fù)載大小分別為5%Ca、10%Ca、15%Ca(Ca為額定動載荷)，具體載荷大小如表2所示.

表1 滾珠絲杠副的參數(shù)Table 1 Parameters of the ball screw

表2 加速試驗應(yīng)力組合方案Table 2 Stress combination scheme of accelerated test

考慮到試驗對象為小子樣單樣本滾珠絲杠副，為了滿足虛擬增廣法的適用要求，提高其可靠度壽命的評估精度，提出在不同步進(jìn)加速應(yīng)力組合下的試驗截尾時間設(shè)定準(zhǔn)則：在每個步進(jìn)加速雙應(yīng)力組合下出現(xiàn)至少2～3個非關(guān)聯(lián)性故障，并將最后一個故障時間記為該步進(jìn)加速雙應(yīng)力組合下的試驗截尾時間，即當(dāng)在第1個步進(jìn)加速雙應(yīng)力下發(fā)生至少2～3個非關(guān)聯(lián)性故障時，提升應(yīng)力，繼續(xù)在第2個步進(jìn)加速雙應(yīng)力水平下進(jìn)行試驗，以此類推，直至試驗結(jié)束.

整個試驗過程中環(huán)境溫度保持在(20±5) ℃，濕度最大為90%，要定期進(jìn)行潤滑，潤滑條件包括潤滑劑種類、潤滑間隔時間、潤滑量等必須按照具體說明書執(zhí)行.首先對絲杠進(jìn)行加載跑合，排除由于設(shè)計等問題引起的早期故障，確保在各個加速應(yīng)力載荷階段無異常，然后在滾珠絲杠副可靠性試驗儀器上進(jìn)行模擬加載，對其進(jìn)行較為全面的試驗，在試驗過程中對摩擦力矩、噪聲、振動、溫升等性能參數(shù)進(jìn)行實時監(jiān)測和記錄，如果滾珠絲杠副出現(xiàn)超過允許值的噪聲和振動，則判定為故障，修復(fù)或更換受損部件后繼續(xù)進(jìn)行試驗并記錄時間作為一次壽命點.如果經(jīng)過修復(fù)后仍無法正常工作即發(fā)生不可修復(fù)故障，則視為失效，停止試驗.

對可靠性進(jìn)行評估時，關(guān)鍵的步驟之一是弄清故障的判定和分析方法，這直接影響到評估結(jié)果的準(zhǔn)確性.對于滾珠絲杠副來說，故障是指產(chǎn)品在規(guī)定的工作運(yùn)行條件下，一個或某幾個功能喪失，或是性能參數(shù)超過允許界限，也包括滾珠絲杠副螺母、滾珠、回珠器等元部件的破損和斷裂.具體故障判據(jù)如下[12]：

1) 滾珠絲杠副的任一零部件發(fā)生松動、損傷、破裂或局部斷裂；

2) 檢測系統(tǒng)采集到的摩擦力矩、振動、噪聲等數(shù)值在某一段監(jiān)測時間間隔內(nèi)偏離正常范圍；

3) 絲杠傳動不平穩(wěn)，無法正常運(yùn)行.

根據(jù)以上判據(jù)對可靠性試驗過程中滾珠絲杠副出現(xiàn)的故障進(jìn)行統(tǒng)計分析，并將其按照文獻(xiàn)[13]中表5.2滾珠絲杠副的關(guān)聯(lián)故障和非關(guān)聯(lián)故障完成分類.

2.2 加速壽命試驗統(tǒng)計分析

本文在定時截尾雙應(yīng)力加速壽命試驗的基礎(chǔ)上，利用虛擬增廣法對試驗數(shù)據(jù)進(jìn)行等價擴(kuò)充以滿足參數(shù)估計方法的要求，提高在單樣本情況下的壽命評估精度.虛擬增廣法的立論依據(jù)：增廣后子樣均值應(yīng)與原來的子樣均值相等；增廣后的子樣標(biāo)準(zhǔn)差應(yīng)與類似件的子樣標(biāo)準(zhǔn)差相等[14].為使虛擬增廣得到的故障時間數(shù)據(jù)樣本更為合理，使用近似公式

(10)

對試驗得到的故障數(shù)據(jù)進(jìn)行增廣至m.式中：T0為樣本均值；c為待定系數(shù)；σ為樣本標(biāo)準(zhǔn)差；T為增廣后的樣本值.

利用近似公式(10)將試驗樣本數(shù)據(jù)增廣至m=10個，那么增廣后的樣本數(shù)據(jù)T1～T10可以表示為

(11)

根據(jù)虛擬增廣法的2個立論依據(jù)可以列出方程：

(12)

(13)

根據(jù)式(12)(13)可求解得c=0.135 0，綜合試驗數(shù)據(jù)及式(12)可以求出增廣后的樣本數(shù)據(jù)分別為

(T0-1.831 0σ)，(T0-1.089 0σ)，
(T0-0.559 0σ)，(T0-0.241 0σ)，
(T0-0.135 0σ)，(T0+0.135 0σ)，
(T0+0.241 0σ)，(T0+0.559 0σ)，
(T0+1.089 0σ)，(T0+1.831 0σ)

然后，將虛擬增廣后的數(shù)據(jù)利用最小二乘法進(jìn)行統(tǒng)計分析，具體做法如下.

將累積故障分布函數(shù)做如下線性變換：

(14)

其中，令

(15)

則式(14)可轉(zhuǎn)換為Weibull分布線性模型

y=bx+a

(16)

在每個步進(jìn)加速雙應(yīng)力情況下，通過試驗得到的故障時間對應(yīng)的累積故障概率采用近似中位秩公式表示，其經(jīng)驗分布函數(shù)為

(17)

式中n為樣本數(shù)據(jù)總數(shù).那么對于一組增廣后的試驗數(shù)據(jù)

(ti,Fn(ti)),i=1,2,…,n

(18)

根據(jù)式(15)可轉(zhuǎn)換為

(19)

將線性模型y=bx+a通過最小二乘法進(jìn)行擬合，可以求出待估參數(shù)b、a的表達(dá)式分別為

(20)

(21)

綜上，求出Weibull分布中兩參數(shù)的點估計為

(22)

求出模型參數(shù)后進(jìn)行模型驗證，判定故障間隔時間是否符合Weibull分布模型，常用的驗證方法有擬合優(yōu)度檢驗和假設(shè)檢驗.

1) 擬合優(yōu)度檢驗

可用相關(guān)系數(shù)r表明數(shù)據(jù)xi、yi的擬合情況，其表達(dá)式為

(23)

相關(guān)系數(shù)越趨近于1，說明壽命分布服從Weibull模型越好.顯著性水平為α?xí)r通過查表確定相關(guān)系數(shù)起碼值rα.當(dāng)r>rα?xí)r，認(rèn)為回歸效果顯著；反之，認(rèn)為回歸不顯著.可以據(jù)此來判定故障數(shù)據(jù)服從Weibull分布的好壞情況.

2) 假設(shè)檢驗

基本思想：先提出兩個對立的假設(shè)，根據(jù)試驗數(shù)據(jù)計算出檢驗統(tǒng)計量，然后根據(jù)檢驗統(tǒng)計量決定其中某一個假設(shè)成立，完成檢驗.常用的模型檢驗方法有卡方檢驗和柯爾莫戈洛夫(Kolmogorov-Smirnov,KS)檢驗[15].其中，卡方檢驗對小樣本的檢驗精度不如柯爾莫戈洛夫檢驗，因此選用KS檢驗法對試驗?zāi)Ｐ瓦M(jìn)行假設(shè)檢驗.

KS檢驗法將經(jīng)驗分布函數(shù)Fn(t)和統(tǒng)計分析得到的模型累積分布函數(shù)F0(t)進(jìn)行對比，原假設(shè)H0與備擇假設(shè)H1如下.

H0：故障間隔時間服從Weibull分布

H1：故障間隔時間不服從Weibull分布

其檢驗統(tǒng)計量

Dn=sup|F0(ti)-Fn(ti)|=max{di}

(24)

通過查KS檢驗臨界值表可以得到在顯著性水平為α?xí)rDnα的值.如果有Dn≤Dnα，則H0成立，故障間隔時間數(shù)據(jù)服從Weibull分布；反之，H1成立，Weibull壽命分布模型不適用.

3 加速壽命試驗仿真算例

對樣本進(jìn)行加速壽命仿真試驗，其結(jié)構(gòu)參數(shù)如表1所示，所施加的載荷應(yīng)力如表2所示.然后借助數(shù)學(xué)方法和計算機(jī)得到需要的模擬故障數(shù)據(jù)，其具體步驟如下.

對Weibull累積分布函數(shù)式(1)進(jìn)行逆變換得到其反函數(shù)為

(25)

其中隨機(jī)變量y服從均勻分布U(0,1)，y與1-y為同分布，所以可轉(zhuǎn)換為

(26)

可以通過計算機(jī)MATLAB軟件產(chǎn)生基于均勻分布U(0,1)的隨機(jī)數(shù)y，然后計算出F-1(y)的值.

結(jié)合上述方法以及截尾時間設(shè)定準(zhǔn)則，根據(jù)文獻(xiàn)[16]中的滾珠絲杠副的形狀參數(shù)m=1.65，特征壽命η分別取800、480、350 h得到在3個步進(jìn)加速雙應(yīng)力水平下的模擬故障間隔時間數(shù)據(jù)并按從小到大排列，如表3所示.

表3 加速應(yīng)力下的模擬數(shù)據(jù)Table 3 Simulation data under accelerated stress

對加速應(yīng)力下的故障間隔時間數(shù)據(jù)進(jìn)行分析，根據(jù)式(11)(12)(13)得到每一個步進(jìn)加速雙應(yīng)力下的虛擬增廣后的數(shù)據(jù)如表4所示.

表4 經(jīng)過虛擬增廣后的數(shù)據(jù)Table 4 Data after virtual augmentation

依據(jù)最小二乘法原理，通過MATLAB軟件分別對3個步進(jìn)加速雙應(yīng)力下增廣后的數(shù)據(jù)進(jìn)行直線擬合，擬合結(jié)果如圖1所示.

圖1 數(shù)據(jù)擬合曲線Fig.1 Data fitting curve

從圖1可以看出，在3個步進(jìn)加速雙應(yīng)力水平下擬合得到的直線基本保持平行，斜率近似相等，即形狀參數(shù)保持不變，從而驗證了假定2的正確性，即在不同的加速應(yīng)力組合下，樣本的失效機(jī)理沒有發(fā)生改變.

由表5分析可知，3條擬合直線的相關(guān)系數(shù)都很接近于1，說明線性相關(guān)度很高，擬合效果較好.

表5 擬合直線數(shù)據(jù)Table 5 Data of the fitting line

結(jié)合式(22)可求得在3個步進(jìn)加速雙應(yīng)力下的形狀參數(shù)分別為

m1=1.611 0,m2=1.611 4,m3=1.609 8

為滿足假定2，形狀參數(shù)取各加速應(yīng)力組合下所求得的形狀參數(shù)的平均值，即

將m的值代入式(22)中的第2個式子，進(jìn)而求出在3個步進(jìn)加速雙應(yīng)力下的特征壽命分別為

η1=828.92,η2=498.06,η3=361.06

再將以上特征壽命數(shù)據(jù)代入式(5)，聯(lián)立求解可得

A=6.145 3×105,α=0.469 2,β=0.460 4

因此加速方程為

lnηi=13.328 6-0.469 2lnVi-0.460 4lnPi

(27)

將(V0,P0)代入式(27)求出正常應(yīng)力條件下的壽命為η0=1 383.0 h.

根據(jù)式(9)可知滾珠絲杠副在額定載荷下的可靠度函數(shù)為

由式(8)可知滾珠絲杠副的MTBF為

4 應(yīng)用討論

滾珠絲杠副在實際使用中，一般可分為高速輕載、低速重載等情況，以下結(jié)合加速方程對不同工況下的壽命可靠度進(jìn)行分析討論.

由加速方程式(27)可以求出在高速輕載(V3,P0)和低速重載(V0,P3)工況下的特征壽命分別為

η30=721.68,η03=691.93

根據(jù)式(9)可知在這兩種工況下的可靠度函數(shù)分別為

不同工況下的可靠度函數(shù)圖形如圖2所示.

圖2 滾珠絲杠副不同工況下的可靠度函數(shù)Fig.2 Reliability function diagram of ball screw under different working conditions

從圖2中可以看出，在相同可靠度的情況下，滾珠絲杠副在高速輕載下的壽命比在低速重載下的壽命略大，均遠(yuǎn)小于在額定載荷下的可靠度壽命，與實際相符.

5 結(jié)論

1) 為了在滾珠絲杠副可靠性壽命試驗中綜合考慮速度和負(fù)載兩個應(yīng)力水平對可靠性壽命的影響，本文提出一種雙應(yīng)力三水平的加速壽命試驗設(shè)計和可靠性分析方法，通過建立雙應(yīng)力加速壽命試驗?zāi)Ｐ筒⒉捎米钚《朔▽收蠑?shù)據(jù)進(jìn)行統(tǒng)計分析，建立滾珠絲杠副的加速壽命試驗的雙應(yīng)力加速方程.為進(jìn)一步開展?jié)L珠絲杠副實際不同工況條件下壽命評估提供方法和依據(jù).

2) 本文提出利用虛擬增廣法對小子樣單樣本條件下雙應(yīng)力壽命試驗數(shù)據(jù)進(jìn)行等價擴(kuò)充后得到不同應(yīng)力水平下的威布爾分布參數(shù)，并建立加速壽命試驗的加速方程的方法.該方法可以滿足小子樣單樣本試驗條件下滾珠絲桿副可靠性壽命評估的精確性要求.并且，本文給出了仿真算例驗證了所提出方法的有效性.