基于分層貝葉斯模型的齒輪彎曲疲勞試驗(yàn)分析

毛天雨 劉懷舉 王寶賓 侯圣文 陳地發(fā)

1.重慶大學(xué)機(jī)械傳動(dòng)國家重點(diǎn)實(shí)驗(yàn)室，重慶，4000442.陜西法士特齒輪有限責(zé)任公司，西安，710077

0 引言

齒輪廣泛應(yīng)用于航空、航天、艦船、高鐵、海洋裝備等領(lǐng)域，其服役性能直接決定了整機(jī)裝備的可靠性[1]，而且齒輪彎曲疲勞失效導(dǎo)致的裝備事故和人機(jī)安全問題日益突出[2]，因此開展齒輪彎曲疲勞試驗(yàn)非常必要。以試驗(yàn)基礎(chǔ)數(shù)據(jù)為主要支撐是進(jìn)行高可靠齒輪設(shè)計(jì)的必要前提，同時(shí)也是對(duì)齒輪疲勞壽命等服役性能最為直接和有效的評(píng)估手段。

齒輪彎曲疲勞試驗(yàn)基礎(chǔ)數(shù)據(jù)一般指可靠度-應(yīng)力-壽命(P-S-N)曲線，它也是齒輪疲勞強(qiáng)度設(shè)計(jì)的基本參數(shù)。獲取P-S-N曲線通常采用成組法，在4～5個(gè)應(yīng)力級(jí)下進(jìn)行試驗(yàn)，每個(gè)應(yīng)力級(jí)下不少于5個(gè)試驗(yàn)點(diǎn)[3]。但由于齒輪彎曲疲勞試驗(yàn)結(jié)果受齒輪材料微結(jié)構(gòu)特征、加工工藝、表面質(zhì)量等多種不確定性因素的影響，試驗(yàn)壽命數(shù)據(jù)分散明顯，一般情況下難以獲得滿足傳統(tǒng)統(tǒng)計(jì)方法所需要的大樣本數(shù)據(jù)。目前小樣本疲勞試驗(yàn)數(shù)據(jù)的分析方法，一類是通過改進(jìn)大樣本統(tǒng)計(jì)方法，如傅慧民等[4]在假設(shè)疲勞壽命服從對(duì)數(shù)正態(tài)分布的前提下，采用異方差回歸分析的方法對(duì)疲勞試驗(yàn)數(shù)據(jù)進(jìn)行整體分析；XIE等[5]提出了一種基于樣本信息聚集原理的P-S-N曲線數(shù)據(jù)處理方法，通過將在不同應(yīng)力水平下測(cè)試的疲勞壽命轉(zhuǎn)換成在任意應(yīng)力水平下的等效壽命，將小樣本數(shù)據(jù)轉(zhuǎn)換為大規(guī)模樣本；GAO等[6]基于P-S-N曲線在高應(yīng)力區(qū)某一點(diǎn)相交的假設(shè)，實(shí)現(xiàn)了P-S-N曲線的快速獲取。另一類是通過樣本擴(kuò)充，使之盡可能適應(yīng)大樣本數(shù)據(jù)統(tǒng)計(jì)方法。如馬宇鵬等[7]提出了Bootstrap-支持向量回歸-二階累次量的方法框架和多階虛擬樣本容量擴(kuò)充的方法，解決了小樣本分析的問題；趙遠(yuǎn)等[8]基于Bootstrap法充分挖掘試驗(yàn)數(shù)據(jù)的壽命與可靠性總體信息，提高了小子樣數(shù)據(jù)評(píng)估的準(zhǔn)確性。然而樣本數(shù)據(jù)擴(kuò)充方法的本質(zhì)依然是進(jìn)行頻率統(tǒng)計(jì)，當(dāng)樣本量較小時(shí)，進(jìn)行樣本擴(kuò)充可能會(huì)產(chǎn)生誤差。另外，傳統(tǒng)的概率方法往往需要大量的疲勞試驗(yàn)數(shù)據(jù)作為依托，而齒輪疲勞試驗(yàn)耗時(shí)耗力，商用諧振式高頻疲勞試驗(yàn)機(jī)加載頻率一般為50～100 Hz，依據(jù)ISO 6336標(biāo)準(zhǔn)規(guī)定的3×106次循環(huán)基數(shù)，一個(gè)試驗(yàn)點(diǎn)要進(jìn)行17 h(按50 Hz計(jì)算)的試驗(yàn)，這導(dǎo)致通常可用于分析的疲勞試驗(yàn)數(shù)據(jù)往往是小樣本的，而貝葉斯統(tǒng)計(jì)方法可用于推斷疲勞測(cè)試數(shù)據(jù)并處理不確定性，可通過不斷更新和融合數(shù)據(jù)來解決樣本量較少這一問題。GUIDA等[9]借助貝葉斯理論和先驗(yàn)信息，將材料疲勞性能參數(shù)引入先驗(yàn)分布中，發(fā)現(xiàn)無論是點(diǎn)估計(jì)還是區(qū)間估計(jì)，所提出的貝葉斯方法遠(yuǎn)優(yōu)于經(jīng)典頻率理論方法；劉建中等[10]運(yùn)用模糊綜合評(píng)判方法綜合經(jīng)驗(yàn)信息來確定先驗(yàn)分布，同時(shí)基于貝葉斯理論給出了一種由小樣本試驗(yàn)數(shù)據(jù)確定疲勞壽命分布的可靠方法。然而，上述方法均是采用經(jīng)驗(yàn)貝葉斯方法，其先驗(yàn)分布的選取易受到主觀信息的影響，在缺乏先驗(yàn)信息的情況下往往存在爭(zhēng)議。而分層貝葉斯可以避免超參數(shù)的選取，能夠提供一種更具魯棒性的統(tǒng)計(jì)分析手段[11]。

本文開展了8620H鋼表面滲碳齒輪彎曲疲勞試驗(yàn)，并針對(duì)傳統(tǒng)統(tǒng)計(jì)方法處理小樣本齒輪彎曲疲勞試驗(yàn)數(shù)據(jù)擬合失真的問題，基于貝葉斯理論建立了齒輪彎曲疲勞P-S-N曲線的分層貝葉斯模型，給出了在無先驗(yàn)信息下分布參數(shù)的選取方法，通過Gibbs采樣對(duì)模型進(jìn)行求解并獲取模型參數(shù)的后驗(yàn)分布。結(jié)合齒輪彎曲疲勞試驗(yàn)數(shù)據(jù)，以50%和99%可靠度下S-N曲線相對(duì)斜率比為擬合結(jié)果評(píng)價(jià)指標(biāo)，驗(yàn)證了分層貝葉斯模型在小樣本數(shù)據(jù)下的穩(wěn)定性與優(yōu)越性。

1 齒輪彎曲疲勞試驗(yàn)

輪齒彎曲疲勞破壞是齒輪最嚴(yán)重的失效形式之一，航空齒輪的斷齒甚至還會(huì)導(dǎo)致機(jī)毀人亡的慘劇[12]，因此，將齒輪彎曲疲勞失效作為研究對(duì)象，通過試驗(yàn)獲得齒輪的彎曲疲勞壽命數(shù)據(jù)與齒輪彎曲疲勞P-S-N曲線，對(duì)指導(dǎo)高性能齒輪正向設(shè)計(jì)具有顯著的工程意義。

1.1 齒輪試樣

齒輪試樣參數(shù)如表1所示，制造工藝路線為：毛坯鍛造→粗車→半精車與精車→滾齒→滲碳淬火→磨齒。所有試驗(yàn)齒輪均采用相同的加工設(shè)備及加工工藝，且在同一爐進(jìn)行熱處理。

表1 試驗(yàn)齒輪基本參數(shù)

1.2 試驗(yàn)方法

齒輪彎曲疲勞試驗(yàn)通常采用脈動(dòng)加載方式[13]，通過疲勞試驗(yàn)機(jī)上的夾具對(duì)試驗(yàn)齒輪輪齒進(jìn)行脈動(dòng)加載，直至輪齒發(fā)生彎曲疲勞失效或越出(循環(huán)基數(shù)通常為3×106)。試驗(yàn)中，脈動(dòng)載荷僅施加在試驗(yàn)齒輪輪齒上，試驗(yàn)齒輪不做嚙合運(yùn)轉(zhuǎn)。

成組法是齒輪彎曲疲勞試驗(yàn)中的常用方法，該方法通過在4～5個(gè)應(yīng)力級(jí)下進(jìn)行多組壽命測(cè)試來獲得齒輪在各個(gè)應(yīng)力級(jí)下的壽命分布，進(jìn)而獲得齒輪的彎曲疲勞P-S-N曲線。圖1為P-S-N曲線示意圖，分別代表了1%、50%、99%可靠度下的S-N曲線。

1.3 試驗(yàn)設(shè)備

試驗(yàn)采用Zwick 1000 kN彎曲疲勞試驗(yàn)機(jī)和 GB/T 14230—1993標(biāo)準(zhǔn)[14]中脈動(dòng)加載的形式，對(duì)試驗(yàn)齒輪進(jìn)行脈動(dòng)循環(huán)加載，以獲得某恒定載荷級(jí)下的齒輪彎曲疲勞壽命。彎曲疲勞試驗(yàn)機(jī)基本參數(shù)如表2所示，齒輪彎曲疲勞試驗(yàn)現(xiàn)場(chǎng)如圖2所示。

表2 彎曲疲勞試驗(yàn)機(jī)基本參數(shù)

圖2 齒輪彎曲疲勞試驗(yàn)Fig.2 Gear bending fatigue test

1.4 齒輪失效的檢測(cè)與判據(jù)

試驗(yàn)過程中，疲勞試驗(yàn)機(jī)通過脈動(dòng)循環(huán)加載與齒輪發(fā)生共振，其共振頻率和疲勞試驗(yàn)機(jī)與齒輪組成的總剛度和質(zhì)量有關(guān)[15]。當(dāng)齒輪出現(xiàn)裂紋后，系統(tǒng)剛度減小，從而引起頻率的下降。根據(jù)GB/T 14230—1993標(biāo)準(zhǔn)[14]的要求，當(dāng)滿足以下兩者之一的條件即可判定齒輪發(fā)生失效：①出現(xiàn)可見裂紋或斷齒；②頻率下降5%～10%。

1.5 試驗(yàn)數(shù)據(jù)的預(yù)處理

(1)

水邊植物群落和建筑周邊植物群落在各指標(biāo)上都有一定程度的下降,根據(jù)駁岸類型和與建筑的位置關(guān)系進(jìn)一步分析(圖6)。

2 試驗(yàn)數(shù)據(jù)處理方法

齒輪彎曲疲勞試驗(yàn)數(shù)據(jù)常用的統(tǒng)計(jì)方法有正態(tài)分布、對(duì)數(shù)正態(tài)分布、雙參數(shù)威布爾分布和三參數(shù)威布爾分布[19]。由于模型可及性和許多計(jì)算工具易于找到，對(duì)數(shù)正態(tài)分布在疲勞試驗(yàn)數(shù)據(jù)處理中被廣泛使用[20-21]，因此，在本文數(shù)據(jù)處理中，采用對(duì)數(shù)正態(tài)分布擬合，通過成組法數(shù)據(jù)確定疲勞曲線的有限壽命階段，即齒輪的P-S-N曲線方程。

2.1 基于頻率理論的P-S-N曲線擬合

在有限壽命階段，可靠度P下的失效循環(huán)次數(shù)NP與應(yīng)力S滿足Basquin方程[22]：

SmPNP=CP

(2)

式中，mP、CP為材料參數(shù)，mP>0，CP>0。

通過對(duì)式(2)左右兩邊取對(duì)數(shù)，可以得到線性表達(dá)式：

lgNP=lgCP-mPlgS

(3)

由式(3)可以看出，Basquin方程在雙對(duì)數(shù)坐標(biāo)下呈現(xiàn)出線性關(guān)系。

如圖3所示，獲取齒輪彎曲疲勞P-S-N曲線，首先需要確定各應(yīng)力水平下的疲勞壽命分布，其次需要確定各應(yīng)力水平下可靠度P下的疲勞壽命NP，進(jìn)而通過最小二乘法(least squares method，LSE)對(duì)式(3)進(jìn)行擬合以獲取P-S-N曲線，具體計(jì)算流程可參考GB/T 14230—1993標(biāo)準(zhǔn)[14]。

圖3 P-S-N曲線擬合原理圖Fig.3 P-S-N curve fitting principle diagram

2.2 基于分層貝葉斯模型的P-S-N曲線擬合

從貝葉斯統(tǒng)計(jì)的觀點(diǎn)來看，模型參數(shù)也被看作是具有概率分布的隨機(jī)變量[23]。分層貝葉斯方法(hierarchical Bayesian model，HBM)相對(duì)于經(jīng)驗(yàn)貝葉斯方法的優(yōu)勢(shì)在于分層先驗(yàn)中包含了模型參數(shù)的結(jié)構(gòu)信息[24]，同時(shí)避免了超參數(shù)的選擇，能夠提供一種更具魯棒性的分析手段。本文將齒輪彎曲疲勞壽命數(shù)據(jù)建模分為兩個(gè)層次，如圖4所示，一個(gè)層次是指每個(gè)應(yīng)力水平下的疲勞壽命，另一個(gè)層次是不同的應(yīng)力水平下的疲勞壽命。HBM模型將各應(yīng)力級(jí)下疲勞壽命參數(shù)分布估計(jì)與S-N曲線的線性擬合同時(shí)進(jìn)行，并將計(jì)算應(yīng)力水平之外的疲勞壽命分布信息考慮在內(nèi)，以進(jìn)行數(shù)據(jù)的融合與交換，從而獲取最優(yōu)擬合結(jié)果。

圖4 分層貝葉斯模型估計(jì)S-N曲線示意圖Fig.4 Schematic diagram of S-N curve estimated byhierarchical Bayesian model

本文中齒輪彎曲疲勞S-N曲線采用目前最為常用的Basquin模型，若考慮各應(yīng)力水平下疲勞壽命的隨機(jī)性，則Basquin方程可表達(dá)為

(4)

其中，m、C為材料參數(shù)；Nj、Sj分別為第j個(gè)應(yīng)力水平下疲勞壽命與應(yīng)力(j=1,2,…，n;n為應(yīng)力水平數(shù))；隨機(jī)變量εj表示在不同應(yīng)力水平Sj下的隨機(jī)性，包含材料疲勞特性的不確定性與觀測(cè)誤差。兩邊同時(shí)取對(duì)數(shù)，可得

lgNj=lgC-mlgSj+δj

(5)

Yj=Xjβ+δj

(6)

Yj=lgNjXj=(1,lgxj)β=(β0,β1)T

圖5 分層貝葉斯層次結(jié)構(gòu)圖Fig.5 Hierarchical Bayesian hierarchy diagram

后驗(yàn)分布的計(jì)算通常采用數(shù)值積分方法求解，馬爾可夫鏈蒙特卡羅(Markov chain Monte Carlo，MCMC)抽樣技術(shù)[29]中的Gibbs抽樣算法適合分層貝葉斯模型，能有效解決高維積分問題[30]，因此本文采用Gibbs抽樣算法來求解后驗(yàn)分布。

在樣本量較少的情況下，可能存在低應(yīng)力水平下對(duì)數(shù)疲勞壽命標(biāo)準(zhǔn)差比高應(yīng)力水平下對(duì)數(shù)疲勞壽命標(biāo)準(zhǔn)差小的情形，與目前普遍認(rèn)為的異方差或同方差理論相悖[31-32]?？紤]到本試驗(yàn)中齒輪來自于相同的材料與熱處理加工工藝，將不同應(yīng)力水平下的疲勞壽命分散性用統(tǒng)一的變異系數(shù)來度量，以此來表征齒輪材料及加工工藝的一致性。統(tǒng)一的變異系數(shù)為

(7)

從而各應(yīng)力水平下壽命分布的方差可調(diào)整為

(8)

2.3 模型的比較與驗(yàn)證

為了驗(yàn)證HBM模型估計(jì)P-S-N曲線的性能，將HBM模型與傳統(tǒng)LSE模型進(jìn)行穩(wěn)定性比較。本節(jié)提出了相對(duì)斜率指標(biāo)來進(jìn)行模型穩(wěn)定性的比較。假設(shè)b0.5表示可靠度為50%下S-N曲線的斜率，b0.99表示可靠度為99%下S-N曲線的斜率，則相對(duì)斜率比指標(biāo)可定義為

(9)

(a)α=1

(b)α>1

(c)α<1圖6 99%可靠度與50%可靠度下S-N曲線的相對(duì)位置關(guān)系示意圖Fig.6 Schematic diagram of the relative positionrelationship of the S-N curve under 99%reliability and 50% reliability

式(9)中的相對(duì)斜率比α描述了50%可靠度下與99%可靠度下S-N曲線的相對(duì)位置關(guān)系，如圖6所示。當(dāng)相對(duì)斜率比α=1時(shí)，99%可靠度下S-N曲線與50%可靠度下S-N曲線平行，如圖6a所示；當(dāng)相對(duì)斜率比α>1時(shí)，99%可靠度下S-N曲線與50%可靠度下S-N曲線相交于高應(yīng)力級(jí)下的一點(diǎn)，如圖6b所示，該點(diǎn)代表該應(yīng)力級(jí)下試件必定會(huì)發(fā)生疲勞失效，這與高應(yīng)力級(jí)下易發(fā)生彎曲疲勞失效的事實(shí)一致；當(dāng)相對(duì)斜率比α<1時(shí)，99%可靠度下S-N曲線與50%可靠度下S-N曲線相交于低應(yīng)力級(jí)下的一點(diǎn)，如圖6c所示，顯然，這與較低應(yīng)力下存在越出點(diǎn)的事實(shí)是相悖的。此時(shí)，P-S-N曲線擬合出現(xiàn)失真。因此，可以接受α≥1的情況而不能接受α<1的情況。

3 結(jié)果與討論

本文采用的齒輪彎曲疲勞試驗(yàn)數(shù)據(jù)，是在4個(gè)應(yīng)力水平(該齒輪彎曲疲勞強(qiáng)度極限約為665 MPa)下進(jìn)行試驗(yàn)獲得的數(shù)據(jù)，每組試驗(yàn)11次，試驗(yàn)結(jié)果見表3。每組11個(gè)試驗(yàn)點(diǎn)提供了足夠的數(shù)據(jù)量以作為評(píng)價(jià)小樣本條件下擬合效果的基準(zhǔn)。

表3 齒輪彎曲疲勞試驗(yàn)結(jié)果

對(duì)于“11-11-11-11”(即4個(gè)應(yīng)力級(jí)下11個(gè)試驗(yàn)點(diǎn))樣本，LSE模型已被證明在大樣本情況下有效，因此將其所獲得的結(jié)果作為基準(zhǔn)曲線。由圖7a可以看出，HBM模型獲得的S-N曲線與LSE模型獲得的S-N曲線幾乎重合。以相同應(yīng)力水平下LSE模型的疲勞壽命作為基準(zhǔn)，絕對(duì)誤差計(jì)算結(jié)果如圖7b所示，四個(gè)應(yīng)力級(jí)下擬合壽命最大誤差為6.10%，表明本文所提出的HBM模型在“大樣本數(shù)據(jù)”情況下具有與LSE模型相同的建模精度。

為對(duì)比本文所提出的HBM模型與傳統(tǒng)LSE模型統(tǒng)計(jì)的穩(wěn)定性，將疲勞試驗(yàn)數(shù)據(jù)分為“3-3-3-3”、“5-5-5-5”、“7-7-7-7”、“9-9-9-9”四種樣本。在每種測(cè)試方案下，用于比較的子數(shù)據(jù)集都是從表3中完整的數(shù)據(jù)中進(jìn)行隨機(jī)抽樣的，每個(gè)測(cè)試方案下隨機(jī)抽樣104次，表5給出了不同樣本量下的統(tǒng)計(jì)結(jié)果，可以看出，在不同樣本量下，本文所提出的HBM模型的相對(duì)斜率比α均能給出滿意的結(jié)果，而傳統(tǒng)的LSE模型會(huì)出現(xiàn)P-S-N曲線擬合失真的情況。這是因?yàn)楫?dāng)疲勞壽命數(shù)據(jù)不滿足隨著應(yīng)力水平逐漸降低而疲勞壽命方差逐漸增大的規(guī)律時(shí)，傳統(tǒng)的LSE模型不能給出滿意的擬合結(jié)果，而本文所提出的模型穩(wěn)定性更好而不依賴于數(shù)據(jù)樣本本身。

表4 齒輪彎曲疲勞試驗(yàn)數(shù)據(jù)預(yù)處理

(a)50%可靠度下S-N曲線

(b)擬合壽命絕對(duì)誤差圖7 11-11-11-11樣本量下50%可靠度下S-N曲線與每個(gè)應(yīng)力水平壽命的誤差Fig.7 The error between the S-N curve and the lifeof each stress level at 50% reliability under 11-11-11-11sample size

其次，LSE與HBM模型相對(duì)斜率比α波動(dòng)范圍均隨樣本量增加而減小，但本文所提出的HBM模型波動(dòng)范圍較小。隨著樣本量的變化，LSE模型最小波動(dòng)范圍為0.7869，是本文所提出的HBM模型最大波動(dòng)范圍(即0.2213，樣本量為“3-3-3-3”)的3.6倍左右。顯然，與傳統(tǒng)的LSE模型相比，本文所提出的HBM模型具有更好的穩(wěn)定性。

表5 不同樣本量下的相對(duì)斜率比α值

(a)不同樣本量下LSE模型α值

(b)不同樣本量下HBM模型α值圖8 不同樣本量下LSE與HBM模型α值Fig.8 LSE and HBM model α values under differentsample sizes

由圖8a和圖8b中的箱型圖可以看出，隨著樣本量的變化下，本文提出的HBM模型相對(duì)斜率比α均大于1，而LSE模型在樣本量為“3-3-3-3”和“5-5-5-5”時(shí)易出現(xiàn)擬合失真的情況。當(dāng)樣本量從“3-3-3-3”增加到“9-9-9-9”時(shí)，HBM模型相對(duì)斜率比α的中位值從1.1392增大到1.1527，變化率為1.19%；而LSE模型中位值從0.8789增大到1.2494，變化率為42.14%。結(jié)合目前彎曲疲勞試驗(yàn)數(shù)據(jù)可知，本文所提出的HBM模型在小樣本下擬合效果更佳，不會(huì)出現(xiàn)擬合失真的情況。

4 結(jié)論

本文開展了齒輪彎曲疲勞試驗(yàn)，并基于齒輪疲勞試驗(yàn)結(jié)果提出了用于估計(jì)齒輪彎曲疲勞P-S-N曲線的分層貝葉斯(HBM)模型，通過馬爾可夫鏈蒙特卡羅(MCMC)數(shù)值仿真得出以下結(jié)論：

(1)采用單齒加載方式進(jìn)行齒輪彎曲疲勞試驗(yàn)，獲得了8620H鋼滲碳齒輪材料在701 MPa、751 MPa、826 MPa、951 MPa四個(gè)彎曲應(yīng)力級(jí)下的彎曲疲勞壽命數(shù)據(jù)，并建立了HBM模型與LSE模型進(jìn)行疲勞試驗(yàn)數(shù)據(jù)的分析與處理。

(2)定義了用于評(píng)價(jià)P-S-N曲線擬合精度的相對(duì)斜率比指標(biāo)。發(fā)現(xiàn)在大樣本情況下，HBM模型具有與傳統(tǒng)LSE模型相同的建模精度。隨著樣本量的變化，LSE模型相對(duì)斜率比α變化率為42.14%，而本文所提出的HBM模型變化率僅為1.19%，具有更好的穩(wěn)定性與適應(yīng)性，為齒輪彎曲疲勞試驗(yàn)等小樣本數(shù)據(jù)場(chǎng)合提供了更具魯棒性的分析手段。