李淑艷 翟友邦 王小龍 李若晨 楊子涵 宋正河

(中國農(nóng)業(yè)大學(xué) 工學(xué)院，北京 100083)

載荷譜是機械零部件壽命設(shè)計和可靠性分析的基礎(chǔ)[1]。傳動軸作為拖拉機傳動系的主要組成部分，承擔(dān)著從變速箱向驅(qū)動橋傳遞動力的功能，是影響拖拉機整機可靠性的關(guān)鍵零部件。編制傳動軸載荷譜，可以提高疲勞試驗的可信度，對提高傳動系可靠性具有重要意義。出于成本考慮，載荷測試試驗時間相對較短，需要對測試載荷進行外推獲得全壽命載荷譜。

常用的外推方法有時域外推和雨流矩陣外推[2]。其中，時域外推能夠保留載荷時間順序，適用于平穩(wěn)性載荷，當(dāng)載荷表現(xiàn)出非平穩(wěn)性時，其計算量較大。雨流矩陣外推雖不能保留載荷時間順序，但其計算量小，對平穩(wěn)和非平穩(wěn)載荷均適用，具有較強工程應(yīng)用價值。由于拖拉機作業(yè)環(huán)境復(fù)雜，導(dǎo)致傳動軸承受非平穩(wěn)性載荷，因此適合采用雨流矩陣外推的方法。

雨流矩陣外推常用的方法有參數(shù)法和非參數(shù)法。參數(shù)法通常采用正態(tài)分布擬合均值頻次，威布爾分布擬合幅值頻次。賈海波[3]基于單參數(shù)法對輪式裝載機傳動系載荷進行外推，在此基礎(chǔ)上獲得了均幅值二維載荷譜，并給出了疲勞試驗一維程序譜加載方案。已有研究[4-5]通過參數(shù)法對輪式裝載機半軸載荷進行外推，其中，翟新婷提出的混合分布的載荷譜編制方法對載荷多峰分布的擬合難題拓寬了解決思路，但擬合精度仍需提高。陳東升等[6]以軍用車輛傳動系統(tǒng)為研究對象，采用雙參數(shù)雨流計數(shù)法，對載荷的幅值和均值進行計數(shù)，并將測試載荷外推得到全壽命歷程二維設(shè)計載荷譜，建立了可用于疲勞壽命估算的八級程序載荷譜。采用基于參數(shù)法的雨流矩陣外推方法時，由于非平穩(wěn)性載荷具有單峰和多峰的特點，存在擬合不通過或分布函數(shù)難以確定的局限。非參數(shù)法是基于核密度估計(Kernel density estimation，KDE)的思想對雨流矩陣進行外推，可以直接得出任意載荷的密度值，克服分布函數(shù)擬合的缺陷，為多峰復(fù)雜載荷的精確擬合提供了解決辦法。李鶯鶯等[7]采用非參數(shù)估計方法對挖掘機液壓泵載荷進行外推，得到全生命周期內(nèi)每個載荷循環(huán)可能出現(xiàn)的頻次，而且可以保證每個遲滯回環(huán)的結(jié)構(gòu)不被破壞。張曉晨等[8]選擇基于載荷擴展的非參數(shù)外推方法分別對液壓挖掘機的4個作業(yè)段進行外推，獲得全壽命下的長期載荷譜。李研[9]提出了基于全帶寬非參數(shù)外推模型和載荷擴展的非參數(shù)雨流外推方法，針對輪式裝載機，根據(jù)實測載荷數(shù)據(jù)進行了試驗驗證與載荷譜編制工作，并指出帶寬對非參數(shù)外推精度的局限性。Johannesson.P等[10]將極值理論與非參數(shù)法結(jié)合起來，對汽車半軸載荷進行了外推處理。

雖然國內(nèi)外都對非參數(shù)外推開展了相關(guān)研究，但外推精度依然亟待提高。針對上述問題，本研究擬以拖拉機傳動軸為研究對象，搭建傳動軸無線扭矩測試系統(tǒng)并進行載荷測試，分析測試載荷分布特征，運用拇指法(Rule-of-thumb，ROT)求帶寬，提出一種改進的基于核密度估計的載荷外推方法對測試載荷進行外推，并從參數(shù)統(tǒng)計、擬合度和偽損傷3個方面對所提出的方法的精度和有效性進行驗證。

1 拖拉機傳動軸田間試驗

試驗樣機為LX2204拖拉機，所搭載的犁耕機型號為MULTI-MASTER153T。在傳動軸中心附近粘貼應(yīng)變片，同時無線扭矩節(jié)點固定在傳動軸上以確保傳感器與傳動軸的相對位置保持不變。選用北京必創(chuàng)科技有限公司生產(chǎn)的TQ201型無線扭矩采集器，其使用應(yīng)變片為BF350-3BA型，基于TQ201搭建的拖拉機傳動軸無線扭矩測試系統(tǒng)見圖1。

圖1 拖拉機傳動軸無線扭矩測試系統(tǒng)Fig.1 Tractor propeller shaft torque test system

田間試驗地點為河南省洛陽市孟津縣金村面積約13 hm2的砂土農(nóng)田，玉米秸稈留茬高度約10 cm。采集了犁耕、旋耕、運輸和聯(lián)合耕整等工況的傳動軸載荷試驗數(shù)據(jù)，采樣頻率為100 Hz。為保證測試載荷能真實反映作業(yè)工況，根據(jù)GB/T 14225—2008《鏵式犁》中的規(guī)定，試驗過程由駕駛員按其習(xí)慣進行操作。試驗中拖拉機的檔位為中三檔，行駛速度為5～10 km/h；耕深和作業(yè)寬幅分別為0.32和2.1 m；碎土率為99.2%。對田間作業(yè)環(huán)境進行測量，風(fēng)速為4.6 m/s，環(huán)境濕度和溫度分別為21%和25.6 ℃；距土壤表面5、10和15 cm處的含水率分別為6.28%、9.84%和12.2%。

2 基于核密度估計的載荷外推方法

本研究采用雨流矩陣外推中基于核密度估計理論的非參數(shù)法和按比例外推方法對測試載荷進行外推處理，并對非參數(shù)法在外推過程中帶寬的選取進行探究，載荷的外推流程見圖2。

圖2 基于核密度估計的載荷外推方法流程Fig.2 Process of load extrapolation method based on kernel density estimation

2.1 信號預(yù)處理

田間試驗共獲得8組犁耕工況傳動軸載荷數(shù)據(jù)，受測試環(huán)境和測試系統(tǒng)等原因影響，田間測試信號中不可避免地會存在干擾噪聲。通過多項式擬合去除測試信號中的趨勢項、幅值門限法和動態(tài)標(biāo)準(zhǔn)差法去除信號中的奇異點以及巴特沃斯低通濾波器濾除高頻載荷，由于高頻載荷大部分是小載荷，因此低通濾波可以將對結(jié)構(gòu)損傷貢獻很小的小載荷濾除。圖3為預(yù)處理后的部分犁耕工況傳動軸測試載荷時間歷程。

圖3 犁耕工況傳動軸扭矩時間歷程Fig.3 Torque time history of propeller shaft in ploughing conditions

2.2 核密度估計

采用核密度估計理論對雨流矩陣進行外推屬于非參數(shù)法。核密度估計能較準(zhǔn)確估計任意分布載荷數(shù)據(jù)的密度值，當(dāng)帶寬選擇恰當(dāng)時，密度值能達到較高的精度，核密度估計的具體步驟為[11]：

1)采用雨流計數(shù)法將載荷-時間歷程轉(zhuǎn)換成雨流矩陣。

2)確定核函數(shù)K(x,y)和帶寬h，構(gòu)建核密度估計模型f(x,y)。

3)根據(jù)核密度估計模型和敏感系數(shù)計算自適應(yīng)因子λi：

(1)

式中：α為敏感系數(shù)，經(jīng)驗值為0.5；n為樣本量，xi、yi分別為雨流矩陣初值(From)和終值(To)。

4)根據(jù)自適應(yīng)因子確定自適應(yīng)帶寬的核密度估計模型：

(2)

首先采用核密度估計模型對雨流矩陣進行密度估計獲得密度矩陣，然后采用Monte Carlo法對密度矩陣進行隨機數(shù)模擬，生成外推載荷。

2.3 帶寬確定

從式(2)可以看出，核函數(shù)和帶寬影響核密度估計精度。研究表明[12-13]，核函數(shù)對核密度估計精度影響較小，帶寬是影響核密度估計的主要因素。帶寬過大，會導(dǎo)致載荷特性丟失；帶寬過小，則會在密度中引入大量虛假波型，造成密度估計不穩(wěn)定。本研究結(jié)合拖拉機傳動軸載荷特征確定ROT帶寬計算方法。

2.3.1R0T帶寬計算方法

本研究選取高斯核函數(shù)作為核密度估計模型的核函數(shù)，并在核密度估計模型的基礎(chǔ)上確定帶寬計算方法。針對傳動軸測試載荷，采用四點雨流計數(shù)法將測試載荷轉(zhuǎn)換成初值(From)- 終值(To)雨流矩陣。雨流矩陣中存在單峰載荷和多峰載荷(圖4)。ROT能根據(jù)單峰載荷和多峰載荷的特性，通過載荷樣本的標(biāo)準(zhǔn)差和標(biāo)準(zhǔn)四分位距(Standardized interquartile range)確定出最優(yōu)帶寬[14-16]。

圖4 測試載荷雨流矩陣分布特征Fig.4 Test load rain flow matrix distribution characteristics

基于高斯核函數(shù)的核密度估計模型的ROT求帶寬h的公式為：

(3)

式中：A為高斯分布的最小值；σ為樣本標(biāo)準(zhǔn)差；IQR為數(shù)據(jù)的標(biāo)準(zhǔn)四分位距。

拖拉機傳動軸載荷具有大載荷較為分散，中小載荷相對集中的特點。對犁耕工況傳動軸測試載荷雨流矩陣中的載荷進行簡化(圖5)。當(dāng)用式(3)計算的帶寬應(yīng)用于a點的密度估計時，帶寬會偏小。在以帶寬為半徑的圓內(nèi)計算a點的概率密度分布，并進行Monte Carlo隨機數(shù)模擬。如圖6所示，由于帶寬偏小，產(chǎn)生的隨機數(shù)會過于集中?？梢赃m當(dāng)增大帶寬以提高隨機數(shù)的分散性，通過減小樣本量來增大帶寬。參照時域外推選取大載荷閾值的原則和實際外推效果，選擇總樣本量的80%計算帶寬，則修正樣本量N為：

a代表分散性較強的大載荷，b代表分布相對集中的中小載荷。Point a represents a large load with strong dispersion, point b represents a relatively concentrated medium and small load.圖5 犁耕工況測試載荷雨流矩陣分布點Fig.5 Rain flow matrix distribution points of test load in plowing conditions

a*代表對a點的大載荷進行密度估計。a* represents the density estimation of the large load at point a.圖6 大載荷核密度估計Fig.6 Rain flow matrix large load kernel density estimation

N=[0.8×n]+1

(4)

2.3.2核密度估計應(yīng)用范圍

應(yīng)用式(1)計算出的載荷分布集中的區(qū)域，其自適應(yīng)因子較小，即發(fā)生概率較大，針對這類載荷采用按比例外推法可以保留其原分布，且外推后仍在該區(qū)域，符合原載荷較集中的特征。同時采用線性外推后，在Monte Carlo模擬時略去此區(qū)域，計算速度可得到大幅提升。

對于雨流矩陣分散性較強(大載荷)的區(qū)域使用基于KDE的非參數(shù)法外推。在Monte Carlo模擬時會在其周圍產(chǎn)生新的分布點(圖7)，新分布點符合原載荷分散性強的特性。

a’ 代表Monte Carlo模擬產(chǎn)生的新分布點。a’ represents the new distribution point generated by Monte Carlo simulation.圖7 KDE外推結(jié)果Fig.7 The results of KDE extrapolation

樣本量也是影響核密度估計精度的一個因素，為確保有較高的外推精度，需保證非參數(shù)法有足夠多的樣本量，因此需要在載荷的分散性和樣本量之間尋求一個合適的閾值。根據(jù)拖拉機傳動軸在各工況下的載荷分布特征及載荷對疲勞損傷的貢獻度確定閾值。

1)計算雨流矩陣中載荷的幅值a：

(5)

式中：F為雨流矩陣的初值；T為雨流矩陣的終值。

2)定義非參數(shù)法和按比例外推法之間的閾值t：

t=[0.1max(a)]+1

(6)

通過上述分析和修正，用于拖拉機傳動軸載荷外推的方法為：

(7)

式中：k為外推倍數(shù)；Ns為小載荷循環(huán)數(shù)；G(x,y)為采用基于KDE的非參數(shù)法的載荷外推函數(shù)。

3 載荷外推

以圖3所示犁耕工況為例進行載荷外推計算。因測試載荷為扭矩，雨流計數(shù)是基于材料應(yīng)力遲滯回線原理形成的，故將扭矩轉(zhuǎn)換成應(yīng)力進行計數(shù)處理，轉(zhuǎn)換公式為：

(8)

式中：τ為應(yīng)力，Pa；T為扭矩，N·m；Wt為抗扭截面系數(shù)，Wt=πd3/16，d為軸的直徑，m。

將犁耕工況傳動軸扭矩轉(zhuǎn)換成應(yīng)力(圖8)，提取測試應(yīng)力數(shù)據(jù)的轉(zhuǎn)折點并用四點雨流計數(shù)法得到雨流載荷循環(huán)的初值、終值及其頻次信息，對載荷進行分級并編制成雨流矩陣。本研究將雨流載荷分成64級，編制成64×64雨流矩陣見圖9。雨流矩陣的原始樣本量n=1 488，根據(jù)式(4)計算出修正樣本量N=1 191。根據(jù)式(5)和式(6)進一步計算出閾值t=8.2 MPa，結(jié)合式(3)可計算出最優(yōu)帶寬為2.709和2.691。非參數(shù)法的外推參數(shù)見表1。

表1 基于核密度估計的非參數(shù)法外推參數(shù)Table 1 The extrapolation of parameters based on nonparametric method of KDE

圖8 梨耕工況傳動軸應(yīng)力時間歷程Fig.8 Stress time history of propeller shaft in ploughing conditions

圖9 基于應(yīng)力數(shù)據(jù)編制的雨流矩陣Fig.9 Rainflow matrix based on stress data

根據(jù)閾值t劃分非參數(shù)法外推載荷和按比例外推載荷。非參數(shù)法外推載荷見圖10(a)，樣本數(shù)量為318，按比例外推載荷見圖10(b)，樣本數(shù)量為1 170。

圖10 犁耕工況傳動軸載荷外推結(jié)果Fig.10 The result of different extrapolation methods of propeller shaft load in ploughing condition

取敏感系數(shù)α=0.2，采用上述核密度估計模型對雨流矩陣進行密度估計，并用Monte Carlo模擬產(chǎn)生與樣本量相同的載荷。把非參數(shù)法外推的載荷和按比例外推的載荷合并形成一個雨流矩陣見圖11,對比圖9可以看出，外推1倍的雨流矩陣載荷分布與原始雨流矩陣載荷分布具有較高的相似度，保留了原有載荷的分布特征。相比于整段載荷使用核密度估計的傳統(tǒng)非參數(shù)法，外推的運算速度得到大幅提高，當(dāng)外推至全壽命頻次(106)時，運算時間從0.502 310 s縮短至0.230 247 s，速度提高了54.16%。

圖11 犁耕工況傳動軸總體外推結(jié)果Fig.11 Total result of extrapolation of propeller shaft in ploughing condition

4 外推結(jié)果驗證

當(dāng)外推方法選擇合理時，外推載荷與測試載荷應(yīng)具有相同的特征，本研究分別從統(tǒng)計參數(shù)、擬合度和偽損傷 3 個方面檢驗外推載荷的有效性。

4.1 統(tǒng)計參數(shù)檢驗

選取統(tǒng)計參數(shù)均值、標(biāo)準(zhǔn)差和最大值作為傳動軸測試載荷與外推載荷的評判依據(jù)。對犁耕工況，采用本研究提出的基于核密度估計的載荷外推方法得到的外推載荷，其均值、標(biāo)準(zhǔn)差和最大值與測試載荷的誤差分別為2.2%、0.87%和6.3%，均在工程應(yīng)用允許誤差范圍內(nèi)。旋耕工況、運輸工況和聯(lián)合耕整工況的均值、標(biāo)準(zhǔn)差和最大值誤差這3種統(tǒng)計參數(shù)均小于10%。綜上表明外推載荷能較好地保留測試載荷的統(tǒng)計特性。各工況測試載荷與外推載荷的統(tǒng)計參數(shù)見表2。

表2 不同工況下傳動軸載荷統(tǒng)計參數(shù)Table 2 Statistical parameters of propeller shaft load under different working conditions MPa

4.2 擬合度檢驗

采用擬合度校驗外推結(jié)果。載荷的循環(huán)次數(shù)和幅值是疲勞損傷的主要來源，通過繪制幅值累計頻次曲線可以反映載荷的幅值和循環(huán)次數(shù)的統(tǒng)計結(jié)果。由外推載荷和測試載荷繪制成的雨流載荷幅值累計頻次對比曲線見圖12。

圖12 犁耕工況幅值累計頻次曲線Fig.12 Accumulated frequency curve of amplitude of ploughing condition

為判斷外推載荷與測試載荷之間的差異，對曲線擬合度進行量化，用決定系數(shù)R2(Coefficient of determination)計算曲線之間的擬合度。對犁耕工況，外推載荷幅值的R2達到0.998，其中前100個大幅值載荷的R2為0.978，表明該方法保留了拖拉機傳動軸載荷的分布特性。采用旋耕工況、運輸工況和聯(lián)合耕整工況下傳動軸測試載荷進一步驗證本研究所提外推方法的合理性，并與傳統(tǒng)非參數(shù)法進行對比，結(jié)果表明該外推方法對傳動軸不同工況載荷具有較高的適應(yīng)性，外推的整段載荷和前100個大幅值載荷的R2均大于0.9，與測試載荷高度擬合。與傳統(tǒng)非參數(shù)法相比，由本研究所提外推方法得到的R2均有較大幅度的提高，外推效果明顯優(yōu)于傳統(tǒng)非參數(shù)法。決定系數(shù)R2見表3。

表3 不同工況下外推載荷的決定系數(shù)Table 3 Determination coefficient R2 of extrapolated load under different working conditions

4.3 偽損傷檢驗

外推載荷的最終目的是編制疲勞試驗加載譜。在機械工程領(lǐng)域，分析復(fù)雜的交變載荷時，常通過計數(shù)法和損傷累積來計算疲勞。采用S-N曲線疲勞壽命模型計算構(gòu)件的疲勞壽命，循環(huán)壽命與恒幅值應(yīng)力的關(guān)系[17]為：

(9)

式中：α為材料屬性參數(shù)；β為損傷指數(shù)；Si為恒幅值應(yīng)力；Ni為循環(huán)壽命。

工程上常用Palmgren-Miner線性損傷理論計算偽損傷，該理論是基于S-N曲線各應(yīng)力循環(huán)次數(shù)對材料造成偽損傷的線性累積結(jié)果[18]，表達式為：

(10)

式中：Di為總偽損傷；ni為材料在Si載荷作用下的循環(huán)次數(shù)；Ni為載荷Si作用下的疲勞壽命循環(huán)次數(shù)。

為便于對偽損傷進行比較，定義偽損傷比為外推載荷對材料造成的偽損傷與測試載荷對材料造成的偽損傷的比值。偽損傷比越接近1，外推載荷與測試載荷的相似度越高。基于Palmgren-Miner線性損傷疲勞累積理論在Ncode軟件中搭建偽損傷計算模型。

由本研究所提外推方法得到的犁耕工況外推載荷，其偽損傷比為0.999，與測試載荷具有高度一致性。同時計算旋耕工況、運輸工況和聯(lián)合耕整工況傳動軸載荷的偽損傷比，其數(shù)值分別為0.992、0.998 和0.996，均能較完整地保留測試載荷的損傷量。

5 結(jié) 論

本研究根據(jù)拖拉機傳動軸載荷特征，在采用ROT求得帶寬的基礎(chǔ)上，提出了一種基于核密度估計的載荷外推方法，并從統(tǒng)計參數(shù)、外推擬合度和偽損傷3個方面驗證了方法的有效性。主要結(jié)論如下：

1)與傳統(tǒng)非參數(shù)法相比，當(dāng)載荷外推至全壽命頻次(106)時，采用本研究所提外推方法可將計算時間縮短54.16%，當(dāng)測試載荷量較多時，能夠有效縮短外推時間，提高計算效率。

2)將犁耕工況的外推載荷與測試載荷進行對比分析。統(tǒng)計參數(shù)檢驗表明，均值、標(biāo)準(zhǔn)差和最大值的誤差分別為2.2%、0.87%和6.3%，外推載荷很好地保留了測試載荷的統(tǒng)計特征；擬合度檢驗表明，前100個大幅值載荷的R2為0.978，整段載荷R2達0.998，外推載荷很好地保留了測試載荷的分布特征；偽損傷檢驗表明，偽損傷比達到0.99，外推載荷很好地保留了測試載荷的損傷量。

3)用本研究所提外推方法對旋耕工況、運輸工況和聯(lián)合耕整工況傳動軸載荷進行外推，統(tǒng)計參數(shù)誤差均小于10%，R2均大于0.9，偽損傷比均大于0.99，所提方法適用于拖拉機傳動軸不同工況載荷外推。