李 揚 劉懷舉 魏沛堂 毛天雨 陳地發(fā)

重慶大學(xué)機械傳動國家重點實驗室，重慶，400044

0 引言

疲勞破壞是重載齒輪的主要失效形式，嚴重影響裝備的壽命及可靠性[1-3]。疲勞強度極限是齒輪等關(guān)鍵零部件設(shè)計的重要指標之一，其定義為給定壽命范圍內(nèi)抵抗疲勞失效的最大許用應(yīng)力，通常通過疲勞試驗得到。Dixon-Mood(D-M)法是分析疲勞試驗數(shù)據(jù)并獲取疲勞強度極限的常用方法[4-5]。起初，該方法的樣本量N大于50的要求通常難以滿足[6]，限制了其應(yīng)用和齒輪疲勞數(shù)據(jù)建設(shè)。LITTLE[7-8]將其推廣應(yīng)用于疲勞極限均值與標準偏差的評估后，研究人員發(fā)現(xiàn)獲得較準確均值并不需要太多樣本，BROWNLEE等[9]發(fā)現(xiàn)樣本量N縮小至5～10時，D-M法的均值μ也是可靠的。SVENSSON[10]通過分析小樣本(N≤30)模擬測試結(jié)果發(fā)現(xiàn)D-M法的標準差估計值σ有偏。張?zhí)祜w等[11]通過大量Monte-Carlo模擬發(fā)現(xiàn)，D-M法估計的均值μ無偏、σ有偏。SVENSSON等[12]提出了標準差的線性矯正因子(以下稱為Svensson-Lorén公式)，并認為該方法可改進包括D-M法在內(nèi)的所有最大似然評估方法?，F(xiàn)行的標準方法[13-15]建議D-M法運用于解釋性研究時的樣本需求量不小于14，比初版D-M法的樣本需求量大大減小。D-M法應(yīng)用的日益廣泛對D-M法本身的誤差分析及控制提出了更高要求。MüLLER等[16]利用4種方法對階梯試驗結(jié)果進行分析，發(fā)現(xiàn)樣本量小于20時，應(yīng)根據(jù)經(jīng)驗或文獻對標準差進行估計，而不是采用D-M法計算。ZHAO等[17]采用最大似然原理對疲勞極限的平均值和標準差進行評估，根據(jù)試驗結(jié)果中成對失效越出樣本的局部S-N關(guān)系，并結(jié)合疲勞極限相關(guān)物理公式和統(tǒng)計理論推導(dǎo)得出疲勞極限數(shù)據(jù)。

盡管上述研究對D-M法進行了探索和修正，但仍未能解決小樣本情況下估計標準差存在偏差的問題。本文采用Monte-Carlo模擬仿真，基于標準差修正和樣本擴充，提出改進疲勞極限標準差估計新方法(Chongqing University bootstrap，CQUboot)[18]，并編寫小樣本分析軟件。該軟件對10組不同試件工藝狀態(tài)和試驗參數(shù)的D-M法數(shù)據(jù)的分析驗證了CQUboot法的適用性。

1 D-M法與小樣本分析方法

1.1 D-M法

疲勞極限可以被理解為具有給定預(yù)期疲勞壽命NC的疲勞強度[13]。以齒輪彎曲疲勞極限為例，綜合考慮適用性和安全性，NC通常定義為3×106次循環(huán)[19]。如圖1所示，D-M法[20]測試首先選擇初始應(yīng)力S0和應(yīng)力階梯步長d，若試樣1在NC之前發(fā)生破壞，則認定為失效，并在S0-d的應(yīng)力水平測試試樣2。若試樣2未發(fā)生破壞則認定為越出，并將試樣3的應(yīng)力水平設(shè)為S0+d。重復(fù)此過程，直至數(shù)據(jù)量足夠，且最后數(shù)據(jù)點推算的下一應(yīng)力水平與初始應(yīng)力S0相等，以實現(xiàn)升降圖閉合，保證整個試驗的應(yīng)力水平都集中在疲勞極限附近。通常以5%左右的齒輪預(yù)估疲勞強度為應(yīng)力階梯步長d，以50%的材料抗拉強度為初始應(yīng)力水平S0。

圖1 D-M法示意圖Fig.1 Schematic diagram of the D-M method

試驗完成后進行數(shù)據(jù)統(tǒng)計處理，將應(yīng)力水平按升序排列。根據(jù)D-M法提供的疲勞極限分布參數(shù)估計式，計算應(yīng)力平均值μ和標準差的估計值σ：

(1)

(2)

式中，i為排序后的應(yīng)力等級，i= 0，1，2，…；fi為應(yīng)力等級i下的失效或越出頻次。

以失效數(shù)計算時，式(1)中取“-”；以越出數(shù)計算時，式(1)中取“+”。

1.2 基于標準差修正和樣本擴充的小樣本分析方法

針對小樣本條件下D-M法標準偏差的估計偏差較大等問題，提出了一種基于標準差修正和樣本擴充的評估方法(CQUboot法)，提高標準差估計值的準確性，減小統(tǒng)計結(jié)果的分散性。

CQUboot法路線如圖2所示，階梯試驗應(yīng)力等級i=3時，采用Svensson-Lorén公式[12]對標準差估計值進行修正：

圖2 CQUboot法技術(shù)路線圖Fig.2 Technical roadmap of CQUboot

(3)

式中，σSL為修正后的估計值；σDM為D-M法的標準差估計值；n為樣本量。

文獻[21]指出疲勞強度分布的標準差是應(yīng)力階梯步長d和估計值σDM的函數(shù)。在此基礎(chǔ)上，本文充分考慮樣本量的影響，構(gòu)建了修正標準差估計值的函數(shù)關(guān)系式：

(4)

式中，σCQU為修正后的標準差估計值；C、D、m為隨樣本量變化的參數(shù)，見表1。

表1 不同樣本量下的參數(shù)值Tab.1 Tab.1 Parameter values of different sample quantity

應(yīng)力等級i=4時，采用式(4)對標準差進行修正，所得平均估計值更接近真實值，但其統(tǒng)計結(jié)果更分散[18]。本文方法在上述標準差修正的基礎(chǔ)上，將Monte-Carlo模擬引入基于失效概率的樣本擴充。針對每一個試驗點，比較試驗點應(yīng)力水平的失效概率與[0，1]區(qū)間內(nèi)的隨機數(shù)，若隨機數(shù)小于失效的概率，則判定該點失效，否則為越出。對任意組試驗所有樣本點重復(fù)該過程，生成對應(yīng)的虛擬階梯測試數(shù)據(jù)，以降低標準差的統(tǒng)計分散性，減小標準差估計誤差。

應(yīng)力等級i≥5時，階梯數(shù)量較大，數(shù)據(jù)更加分散。小步長如d<1.0σ的情況下，結(jié)合標準差修正與樣本擴充的結(jié)果分散程度??；步長d>1.0σ的情況下，結(jié)合D-M法與樣本擴充的效果較好。

CQUboot法在可靠度R下的疲勞極限為

σlim，R=μlim+σ′Φ-1(1-R)

(5)

式中，μlim為疲勞極限的均值；σ′為優(yōu)化后的標準差估計值，σ′=σSL，σDM，σCQU；Φ-1(*)為標準正態(tài)分布反函數(shù)。

2 軟件開發(fā)

圖3 齒輪疲勞極限計算軟件邏輯框圖Fig.3 Logic diagram of gear fatigue strength calculation software

2.1 基于Monte-Carlo模擬的樣本擴充及標準差修正模型的建立

采用Python語言，搭建基于Monte-Carlo模擬的樣本擴充及標準差修正計算模型，并對其進行封裝(供Qt框架調(diào)用)，實現(xiàn)數(shù)據(jù)流互通。該過程需配置Numpy、Scipy和Matplotlib的相應(yīng)環(huán)境。其中，Numpy對輸入分布進行隨機抽樣，實現(xiàn)Monte-Carlo模擬階梯試驗，并生成均勻分布隨機數(shù)，進一步與試驗中每個應(yīng)力水平的失效概率對比，實現(xiàn)樣本擴充。利用Scipy中的優(yōu)化器模塊求解本文模型中的標準差修正函數(shù)等。Matplotlib為本文模型提供多樣的可視化輸出，通過調(diào)用Matplotlib繪制相關(guān)圖形來表示失效概率數(shù)據(jù)(F-S數(shù)據(jù))、標準差和期望的頻率分布等。

一方面，為對比各方法的差異，該模型根據(jù)齒輪彎曲疲勞極限已知的真實分布進行多次Monte-Carlo模擬；另一方面，該模型需通過D-M法對試驗數(shù)據(jù)進行分析并實現(xiàn)樣本擴充和標準差修正。為實現(xiàn)上述兩項功能，將該模型劃分為4個模塊：Monte-Carlo模擬階梯試驗、D-M法分析、樣本擴充、標準差修正，各模塊的輸入輸出如圖4所示。

圖4 程序功能實現(xiàn)技術(shù)路線圖Fig.4 Technical roadmap of program functions

2.2 調(diào)用樣本擴充及標準差修正模型方法

圖5 調(diào)用模型的技術(shù)路線Fig.5 Technical roadmap of calling model

通過上述方法調(diào)用封裝完成的模型后，再利用Qt平臺的windeployqt集成工具將應(yīng)用軟件進行打包發(fā)布，使得軟件可在未配置任何開發(fā)環(huán)境的計算機上使用。

2.3 軟件功能模塊設(shè)計

齒輪疲勞極限計算軟件由數(shù)據(jù)庫管理、常規(guī)數(shù)據(jù)統(tǒng)計、D-M法分析和小樣本分析4個功能模塊組成。數(shù)據(jù)庫管理模塊具備存儲、查詢和修改齒輪疲勞試驗樣本信息的功能。常規(guī)數(shù)據(jù)統(tǒng)計模塊統(tǒng)計分析數(shù)據(jù)庫內(nèi)選取的試驗集，完成試驗數(shù)據(jù)的預(yù)處理。D-M法分析模塊根據(jù)階梯試驗得到各應(yīng)力級下的“越出”和“失效”試驗點的分布，選擇總點數(shù)較少的“越出”或“失效”作為分析事件，計算疲勞極限均值和標準偏差的估計值，求得各可靠度下的疲勞強度極限。小樣本分析模塊在D-M法計算的均值和標準差基礎(chǔ)上，采用CQUboot法進行樣本擴充和標準差修正，最后實現(xiàn)誤差分析和小樣本齒輪疲勞極限計算。相比于D-M法，該模塊所需輸入樣本少且計算更準確。軟件主界面如圖6所示。

圖6 軟件平臺界面Fig.6 Software platform interface

3 試驗驗證與分析

為對比D-M法和CQUboot法的分析效果，并測試齒輪疲勞極限計算軟件的性能，分析了10組不同工藝、不同加載步長的齒輪彎曲疲勞極限試驗數(shù)據(jù)，驗證CQUboot小樣本分析方法在齒輪彎曲疲勞極限評估的適用性。

齒輪試件材料為18CrNiMo-7-6，表面均經(jīng)滲碳淬火處理，基本參數(shù)如表2所示。將所有試驗數(shù)據(jù)提前錄入數(shù)據(jù)庫管理模塊，采用D-M法分析模塊計算其中的10組階梯試驗數(shù)據(jù)，得到各組試驗的均值和標準差的估計值，如表3所示，10組數(shù)據(jù)試驗步長范圍是0.58σ～1.89σ，滿足D-M法對步長與標準偏差比值的要求(0.5σ

表2 試驗齒輪基本參數(shù)Tab.2 Basic parameters of test gears

表3 D-M法試驗結(jié)果的參數(shù)估計Tab.3 Parameters estimation of D-M test results

針對上述試驗，采用與試驗相同的階梯設(shè)置，基于Monte-Carlo展開小樣本模擬階梯試驗，將表3中的估計值假設(shè)為試驗的“真實”分布參數(shù)輸入軟件計算模型，分別對每組試驗?zāi)M構(gòu)造12個樣本條件下的階梯試驗1000次，其中的1組試驗結(jié)果如圖7所示。模擬結(jié)果均存在波動，其中，期望估計值μ為[756.84，820.92]MPa，采用D-M法的標準差估計值σDM為[12.43，61.87]MPa，采用CQUboot法的標準差估計值σCQU為[16.15，39.10]MPa，均值和標準差的波動都將導(dǎo)致各可靠度下的疲勞極限在一定范圍內(nèi)波動。

(a)D-M法均值 (b)D-M法標準差

(c)CQUboot法標準差圖7 1000次模擬階梯試驗統(tǒng)計結(jié)果Fig.7 Statistical results of 1000 simulated step tests

計算模型考慮期望估計值μ、標準差估計值σDM和σCQU的波動區(qū)間，計算不同可靠度下D-M法和CQUboot法的疲勞極限估計區(qū)間，得出兩種方法相比“真實”疲勞極限的最大誤差，具體結(jié)果見表4。

表4 12樣本量下不同可靠度的疲勞極限預(yù)估區(qū)間及誤差Tab.4 Fatigue strength error of each group test with 12 sample size under different reliabilities with two methods

利用兩種方法得到了各組試驗樣本量下不同可靠度的疲勞極限預(yù)估誤差，如表4所示，CQUboot法對σFlim的估計值最大誤差始終小于D-M法，CQUboot法的最大誤差17.14%位于第3組試驗的99.99%可靠度下，這是由于該組試驗步長0.58σ較小，導(dǎo)致采用CQUboot法估計的誤差偏大。其他各組試驗在99.99%可靠度下的疲勞極限最大誤差為8.02%～12.98%，D-M法的疲勞極限最大誤差為15.76%～24.99%。樣本減少時，CQUboot法的分析效果好，各可靠度下的疲勞極限更加接近“真實”值。

為進一步說明相同精度下CQUboot方法的樣本需求量相比D-M法的減小程度，利用軟件中的小樣本分析和D-M法分析模塊，對上述10組階梯試驗數(shù)據(jù)開展各小樣本量條件下的模擬階梯試驗，每組試驗在各樣本量條件下均模擬1000次，其中的一組試驗數(shù)據(jù)如圖8所示，兩方法的標準差估計值均隨樣本量的增大而減小，但σCQU的區(qū)間范圍始終小于σDM，σCQU更接近“真實”標準差14.66 MPa。各樣本量下期望估計值μ的均值與“真實”期望688.61 MPa幾乎一致，但在模擬過程中存在較明顯波動，且區(qū)間范圍隨樣本量增大而減小，因此，需將期望估計值μ造成的誤差考慮在內(nèi)。

(a)標準差

(b)均值圖8 估計標準差和均值的區(qū)間分布Fig.8 Interval distribution of estimated standard deviation and mean

針對不同樣本量的模擬階梯試驗參數(shù)估計區(qū)間及疲勞極限誤差如圖9所示。CQUboot法在所有樣本試驗下計算所得的齒輪彎曲疲勞極限的誤差均小于D-M法。將試驗步長從小到大排列會發(fā)現(xiàn)試驗步長較小的CQUboot法的誤差較大，例如試驗3、10；誤差隨樣本量的增大而減小的幅度較小，其中，第10組試驗樣本8下的模擬試驗誤差最大，為12.66%。試驗2和6的步長較大，但誤差比其他試驗組的小，其中，第2組試驗樣本18下的誤差最小，為5.09%。對比10組數(shù)據(jù)可知，試驗步長d為0.71σ～1.59σ時，CQUboot法比D-M法更優(yōu)。若以12.66%為允許的最大誤差，則D-M法分析需要10～18個樣本，CQUboot方法只需8個樣本，這充分證明CQUboot法對不同工藝、步長的小樣本階梯試驗具有較好效果。

圖9 兩種方法在不同樣本量下疲勞極限的預(yù)測誤差Fig.9 Prediction errors of fatigue strengthof two methods under different sample sizes

4 結(jié)論

(2)將10組不同試件和試驗參數(shù)的彎曲疲勞數(shù)據(jù)用于驗證本文方法。以D-M法在20～22個樣本下計算的99%可靠度的彎曲疲勞極限為基準，樣本數(shù)降至12時，D-M法的最大誤差為15.76%～24.99%，而CQUboot法僅為8.02%～12.98%，因此CQUboot法更適合疲勞極限的小樣本統(tǒng)計分析。

(3)CQUboot法試驗步長較小時的誤差較大，且誤差隨樣本量增大而減小的幅度較?。蝗粢?2.66%為疲勞極限預(yù)估允許的最大誤差，則D-M法至少需要10～18個樣本，CQUboot小樣本分析方法只需8個樣本。