基于Monte-Carlo模擬的小樣本下齒輪疲勞極限計算方法及軟件開發

李 揚 劉懷舉 魏沛堂 毛天雨 陳地發

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

0 引言

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

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

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

1.1 D-M法

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

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

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

(1)

(2)

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

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

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

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

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

圖2 CQUboot法技術路線圖Fig.2 Technical roadmap of CQUboot

(3)

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

文獻[21]指出疲勞強度分布的標準差是應力階梯步長d和估計值σDM的函數。在此基礎上，本文充分考慮樣本量的影響，構建了修正標準差估計值的函數關系式：

(4)

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

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

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

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

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

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

(5)

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

2 軟件開發

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

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

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

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

圖4 程序功能實現技術路線圖Fig.4 Technical roadmap of program functions

2.2 調用樣本擴充及標準差修正模型方法

圖5 調用模型的技術路線Fig.5 Technical roadmap of calling model

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

2.3 軟件功能模塊設計

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

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

3 試驗驗證與分析

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

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

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

表3 D-M法試驗結果的參數估計Tab.3 Parameters estimation of D-M test results

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

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

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

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

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

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

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

(a)標準差

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

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

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

4 結論

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

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