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

重慶大學(xué)機(jī)械傳動(dòng)國家重點(diǎn)實(shí)驗(yàn)室，重慶，400044

0 引言

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

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

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

1.1 D-M法

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

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

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

(1)

(2)

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

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

1.2 基于標(biāo)準(zhǔn)差修正和樣本擴(kuò)充的小樣本分析方法

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

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

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

(3)

式中，σSL為修正后的估計(jì)值；σDM為D-M法的標(biāo)準(zhǔn)差估計(jì)值；n為樣本量。

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

(4)

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

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

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

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

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

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

(5)

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

2 軟件開發(fā)

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

2.1 基于Monte-Carlo模擬的樣本擴(kuò)充及標(biāo)準(zhǔn)差修正模型的建立

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

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

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

2.2 調(diào)用樣本擴(kuò)充及標(biāo)準(zhǔn)差修正模型方法

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

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

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

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

圖6 軟件平臺(tái)界面Fig.6 Software platform interface

3 試驗(yàn)驗(yàn)證與分析

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

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

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

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

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

(a)D-M法均值 (b)D-M法標(biāo)準(zhǔn)差

(c)CQUboot法標(biāo)準(zhǔn)差圖7 1000次模擬階梯試驗(yàn)統(tǒng)計(jì)結(jié)果Fig.7 Statistical results of 1000 simulated step tests

計(jì)算模型考慮期望估計(jì)值μ、標(biāo)準(zhǔn)差估計(jì)值σDM和σCQU的波動(dòng)區(qū)間，計(jì)算不同可靠度下D-M法和CQUboot法的疲勞極限估計(jì)區(qū)間，得出兩種方法相比“真實(shí)”疲勞極限的最大誤差，具體結(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àn)樣本量下不同可靠度的疲勞極限預(yù)估誤差，如表4所示，CQUboot法對(duì)σFlim的估計(jì)值最大誤差始終小于D-M法，CQUboot法的最大誤差17.14%位于第3組試驗(yàn)的99.99%可靠度下，這是由于該組試驗(yàn)步長0.58σ較小，導(dǎo)致采用CQUboot法估計(jì)的誤差偏大。其他各組試驗(yàn)在99.99%可靠度下的疲勞極限最大誤差為8.02%～12.98%，D-M法的疲勞極限最大誤差為15.76%～24.99%。樣本減少時(shí)，CQUboot法的分析效果好，各可靠度下的疲勞極限更加接近“真實(shí)”值。

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

(a)標(biāo)準(zhǔn)差

(b)均值圖8 估計(jì)標(biāo)準(zhǔn)差和均值的區(qū)間分布Fig.8 Interval distribution of estimated standard deviation and mean

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

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

4 結(jié)論

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

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