基于Monte-Carlo擬合的二維隨機(jī)載荷作用下汽車彈簧疲勞耐久性研究

藤瑞品，宋曉琳，劉國(guó)云，曾俊偉

（1.湖南大學(xué)，汽車車身先進(jìn)設(shè)計(jì)制造國(guó)家重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410082； 2.湖南獵豹汽車股份有限公司，長(zhǎng)沙 410100）

前言

汽車在工作時(shí)受到的載荷譜是一個(gè)隨機(jī)過(guò)程，隨機(jī)載荷的表征參數(shù)包含載荷幅值和載荷均值，疲勞耐久的研究應(yīng)基于載荷幅值和均值都為隨機(jī)變量的二維隨機(jī)載荷進(jìn)行。

陳欣等［1］對(duì)汽車傳動(dòng)系統(tǒng)的載荷譜進(jìn)行了研究，基于幅值和頻次的關(guān)系建立了八級(jí)程序載荷譜；高云凱等［2］通過(guò)將二維程序載荷譜簡(jiǎn)化為一維程序載荷譜，建立了車身疲勞試驗(yàn)程序載荷譜；高孝旺［3］在基于載荷幅值為威布爾分布的基礎(chǔ)上分析了三參數(shù)威布爾分布的3個(gè)參數(shù)對(duì)等效載荷和相關(guān)評(píng)價(jià)參數(shù)的影響趨勢(shì)；胡建軍等［4］以隨機(jī)載荷為輸入對(duì)齒輪的疲勞壽命進(jìn)行了試驗(yàn)研究；武瀅和謝里陽(yáng)［5］基于載荷幅值和均值都為隨機(jī)分布的情況下以其聯(lián)合概率密度為基礎(chǔ)建立了疲勞壽命的分布預(yù)測(cè)模型。金星等［6］對(duì)隨機(jī)疲勞裂紋進(jìn)行了研究，提出了一種疲勞裂紋擴(kuò)展的新模型。鄒小理［7］提出了一個(gè)隨機(jī)荷載作用下疲勞裂紋擴(kuò)展的統(tǒng)計(jì)模型。郭虎等［8］通過(guò)對(duì)汽車前橋載荷進(jìn)行采集，得出了在各種路面下等效載荷的分布符合三參數(shù)威布爾分布的結(jié)論。賀小帆等［9］以獨(dú)立多細(xì)節(jié)結(jié)構(gòu)概率失效模型為基礎(chǔ)，進(jìn)行MonteCarlo抽樣，對(duì)獨(dú)立多細(xì)節(jié)結(jié)構(gòu)的壽命分布特性進(jìn)行了檢驗(yàn)；朱穎等［10］提出一種計(jì)算平穩(wěn)高斯荷載作用下不確定結(jié)構(gòu)疲勞損傷的新方法。

目前對(duì)于汽車的疲勞耐久的研究絕大部分是基于循環(huán)加載的方式，對(duì)于一維隨機(jī)載荷作用下的疲勞耐久也有過(guò)一些研究，對(duì)于二維隨機(jī)載荷作用下的疲勞耐久的研究不夠深入，沒有成熟的方法和理論。對(duì)于考慮載荷幅值和載荷均值各自的隨機(jī)特性的疲勞耐久的研究尚未見相關(guān)文獻(xiàn)報(bào)道。

在進(jìn)行疲勞分析時(shí)，一些相關(guān)數(shù)據(jù)由于各種實(shí)際條件限制難以準(zhǔn)確得到，而參數(shù)擬合是一種可行的近似方法。如黃學(xué)偉等［11］提出了一種預(yù)測(cè)材料與結(jié)構(gòu)裂紋在高周疲勞下疲勞裂紋擴(kuò)展速率的數(shù)值模擬新方法；田秀云等［12］提出了一種金屬材料疲勞裂紋擴(kuò)展曲線的擬合方法等。

1 二維隨機(jī)載荷作用下的疲勞損傷

根據(jù)Goodman公式，當(dāng)量載荷Seq的表達(dá)式［13］為

式中：σb為材料的拉伸強(qiáng)度；Sa為載荷幅值；Sm為載荷均值。

設(shè)隨機(jī)變量X和Y相互獨(dú)立，其概率密度函數(shù)分別為fX（x）和fY（y），則的概率密度函數(shù)fZ（z）［14］為

根據(jù)Miner定則［15-17］，累積疲勞損傷的計(jì)算方法為

式中：D為材料的累積疲勞損傷；N為材料受到的循環(huán)載荷總數(shù)量；S為材料受到的當(dāng)量載荷；f（S）為隨機(jī)當(dāng)量載荷的概率密度函數(shù)；Nf為材料在載荷S作用下的疲勞壽命。

2 二維隨機(jī)載荷的當(dāng)量載荷的概率密度函數(shù)

2.1 概率密度函數(shù)的代數(shù)法推導(dǎo)

利用載荷幅值和載荷均值的概率密度以及當(dāng)量載荷和載荷幅值、載荷均值之間的數(shù)學(xué)關(guān)系可對(duì)當(dāng)量載荷的概率密度函數(shù)進(jìn)行推導(dǎo)，如文獻(xiàn)［18］中推導(dǎo)得到了載荷幅值和載荷均值均符合正態(tài)分布的二維隨機(jī)載荷的當(dāng)量載荷的概率密度函數(shù)。

2.2 概率密度函數(shù)的Monte-Carlo擬合

2.2.1 Monte-Carlo擬合方法和分布參數(shù)選擇

采用數(shù)學(xué)方法推導(dǎo)二維隨機(jī)載荷的當(dāng)量載荷的概率密度函數(shù)在實(shí)際應(yīng)用中存在很多局限性：一方面推導(dǎo)過(guò)程一般都比較繁瑣；另外在采用Goodman公式求當(dāng)量載荷的概率密度函數(shù)時(shí)，很多情況下無(wú)法用數(shù)學(xué)方法直接推導(dǎo)得出。

蒙特卡洛法［19-23］也稱為隨機(jī)模擬法，是通過(guò)生成隨機(jī)數(shù)序列（實(shí)際上應(yīng)該為偽隨機(jī)數(shù)）來(lái)對(duì)數(shù)學(xué)問(wèn)題進(jìn)行擬合分析的一種數(shù)值方法。假設(shè)隨機(jī)載荷譜的載荷均值和載荷幅值的分布律和概率密度函數(shù)都已知，則可通過(guò)采用蒙特卡羅法隨機(jī)生成載荷幅值和載荷均值數(shù)據(jù)并得到當(dāng)量載荷的隨機(jī)序列，并分析得到當(dāng)量載荷的隨機(jī)分布特性。

在試驗(yàn)場(chǎng)對(duì)某城市SUV進(jìn)行載荷譜采集，得到40個(gè)通道的應(yīng)變時(shí)間歷程信號(hào)，采用雨流計(jì)數(shù)法［24-28］對(duì)載荷譜時(shí)間歷程曲線進(jìn)行計(jì)數(shù)處理，得到包含不同載荷幅值和載荷均值以及各自所對(duì)應(yīng)的計(jì)數(shù)頻次的載荷譜。進(jìn)行參數(shù)估計(jì)和分布檢驗(yàn)可得，載荷幅值基本上符合兩參數(shù)威布爾分布，載荷均值基本符合正態(tài)分布。

設(shè)載荷幅值Sa符合二參數(shù)威布爾分布，其形狀參數(shù)為α，尺度參數(shù)為β；載荷均值Sm符合正態(tài)分布，分布參數(shù)標(biāo)準(zhǔn)差為σ，均值為μ，根據(jù)隨機(jī)變量的函數(shù)的分布特性，載荷均值Sm的函數(shù)同樣符合正態(tài)分布，分布參數(shù)中，標(biāo)準(zhǔn)差為，均值為

選取用于當(dāng)量載荷概率密度函數(shù)數(shù)值擬合研究的載荷幅值和載荷均值的分布參數(shù)范圍，如表1所示。

表1 分布參數(shù)范圍

為研究不同分布參數(shù)的影響，共選取表中8 085組分布參數(shù)進(jìn)行了數(shù)值擬合的研究。

2.2.2 隨機(jī)載荷序列的生成

采用蒙特卡羅法對(duì)表1中分布參數(shù)進(jìn)行隨機(jī)生成，得到載荷幅值和載荷均值的隨機(jī)序列ui和vi，再采用Goodman公式進(jìn)行一維隨機(jī)載荷的零均值當(dāng)量載荷的轉(zhuǎn)化，得到隨機(jī)當(dāng)量載荷序列為fi，則fi的生成方法為

式中vi為符合正態(tài)分布的載荷均值的函數(shù)所生成的隨機(jī)序列。

2.2.3 參數(shù)估計(jì)

（1）圖形估計(jì)［29］

對(duì)表1中的不同分布參數(shù)的隨機(jī)當(dāng)量載荷序列fi進(jìn)行WPP（威布爾概率圖）分布檢驗(yàn)，檢驗(yàn)結(jié)果發(fā)現(xiàn)采用威布爾分布形式來(lái)對(duì)隨機(jī)載荷序列fi進(jìn)行擬合的效果較好，圖1表示一組參數(shù)仿真后的圖形法估計(jì)參數(shù)的結(jié)果。

圖1 圖形法估計(jì)參數(shù)

（2）分布參數(shù)估計(jì)

采用最大似然法對(duì)符合表1中的二維隨機(jī)載荷的零均值當(dāng)量載荷的擬合分布進(jìn)行未知分布參數(shù)估計(jì)，擬合得到符合威布爾分布的當(dāng)量載荷的分布參數(shù)，包含形狀參數(shù)α和尺度參數(shù)β。

表2為部分?jǐn)M合的分布參數(shù)數(shù)據(jù)。

表2 分布參數(shù)擬合

由表2中的數(shù)據(jù)可以看出，采用通過(guò)蒙特卡羅法生成的隨機(jī)序列對(duì)載荷幅值和載荷均值的分布函數(shù)的分布參數(shù)進(jìn)行擬合，擬合得到的分布參數(shù)值與原始定義的分布參數(shù)值基本相符，由此驗(yàn)證了蒙特卡羅擬合法的正確性。

3 載荷譜數(shù)據(jù)采集和處理

3.1 載荷譜數(shù)據(jù)的采集

采用4個(gè)車輪六分力傳感器、10個(gè)加速度傳感器和38個(gè)應(yīng)變片傳感器，在襄陽(yáng)試驗(yàn)場(chǎng)的石塊路面采集了某城市SUV車型整車主要結(jié)構(gòu)載荷譜。

3.2 載荷譜數(shù)據(jù)的處理［30-34］

原始載荷譜為包含異常點(diǎn)（毛刺）以及連接路面的信號(hào)的時(shí)間歷程曲線，須將異常數(shù)據(jù)和連接路面的數(shù)據(jù)剔除，左前螺旋彈簧經(jīng)處理后的數(shù)據(jù)信號(hào)如圖2所示。

圖2 左前彈簧載荷譜信號(hào)

對(duì)載荷譜時(shí)間歷程采用雨流計(jì)數(shù)法統(tǒng)計(jì)處理，形成包含幅值和均值以及各自的統(tǒng)計(jì)計(jì)數(shù)信息的三維直方圖，左前螺旋彈簧的載荷譜三維直方圖如圖3所示。

圖3 左前彈簧載荷三維統(tǒng)計(jì)直方圖

左前螺旋彈簧的載荷幅值和均值各自的二維統(tǒng)計(jì)直方圖如圖4和圖5所示。

3.3 統(tǒng)計(jì)模型的圖形法檢驗(yàn)［29-35］

圖4 左前彈簧載荷幅值統(tǒng)計(jì)直方圖

圖5 左前彈簧載荷均值統(tǒng)計(jì)直方圖

為了確定載荷譜數(shù)據(jù)是否服從威布爾分布和正態(tài)分布，可以用相應(yīng)的概率圖去檢驗(yàn)。

經(jīng)過(guò)雨流計(jì)數(shù)法統(tǒng)計(jì)右前螺旋彈簧的載荷譜的幅值和均值，進(jìn)行威布爾和正態(tài)分布的圖形檢驗(yàn)的結(jié)果見圖6和圖7。

圖6 威布爾分布檢驗(yàn)圖

通過(guò)觀察發(fā)現(xiàn)，圖6中除了部分比較小的幅值和總體趨勢(shì)不一致外，其它基本符合威布爾分布的假設(shè)。而通過(guò)觀察圖7可知，均值非常符合正態(tài)分布的假設(shè)，這與圖5的直方圖一致。

采用極大似然對(duì)模型的分布參數(shù)進(jìn)行估計(jì)［14］，部分參數(shù)估計(jì)的結(jié)果見表3。

圖7 正態(tài)分布檢驗(yàn)圖

表3 統(tǒng)計(jì)分布參數(shù)

其它測(cè)試點(diǎn)的載荷分布特性與表中類似，即載荷幅值符合威布爾分布，載荷均值符合正態(tài)分布。

4 疲勞損傷與耐久壽命分析

4.1 載荷譜數(shù)據(jù)樣本的截除

根據(jù)參數(shù)估計(jì)對(duì)載荷譜進(jìn)行圖形法檢驗(yàn)時(shí)，發(fā)現(xiàn)低載荷區(qū)域和高載荷區(qū)域都發(fā)生了較大的偏離，由于高載荷區(qū)域?qū)φ麄€(gè)疲勞壽命的預(yù)測(cè)影響較大，因此采用這種方法外推的載荷譜用于疲勞損傷計(jì)算會(huì)產(chǎn)生較大的誤差。

由于在載荷譜采集時(shí)采集了大量的小載荷，絕大部分小載荷低于疲勞極限，因此并不會(huì)對(duì)疲勞損傷帶來(lái)影響，而大量的小載荷會(huì)對(duì)隨機(jī)分布的參數(shù)估計(jì)造成較大影響，因此采用將小載荷去除的方法將樣本進(jìn)行截?cái)?，采用截?cái)嗪蟮臉颖具M(jìn)行隨機(jī)分布的參數(shù)估計(jì)。

圖8～圖10分別為右前螺旋彈簧載荷數(shù)據(jù)在未截除狀態(tài)、截除100 MPa以下和截除400 MPa以下數(shù)據(jù)后的估計(jì)參數(shù)檢驗(yàn)圖形。

由圖可知，右前螺旋彈簧的載荷數(shù)據(jù)在截除100 MPa以下的小載荷后，估計(jì)參數(shù)在中高載荷部分可實(shí)現(xiàn)很高的擬合度。

按同樣的方法對(duì)表3中數(shù)據(jù)進(jìn)行截取后的擬合參數(shù)見表4。

圖8 原始數(shù)據(jù)參數(shù)估計(jì)

圖9 截除100 MPa以下小載荷

圖10 截除400 MPa以下載荷

對(duì)樣本截除后擬合參數(shù)進(jìn)行參數(shù)檢驗(yàn)，在中高載荷部分均可實(shí)現(xiàn)很高的擬合度。

4.2 隨機(jī)載荷譜的外推

以石塊路工況的4個(gè)螺旋彈簧零件的載荷譜數(shù)據(jù)為研究對(duì)象，對(duì)隨機(jī)載荷譜進(jìn)行外推。石塊路總長(zhǎng)度為2.7 km，耐久試驗(yàn)總目標(biāo)里程目標(biāo)值為10 000 km。4個(gè)彈簧載荷在樣本截取后的單次試驗(yàn)循環(huán)里程的載荷樣本數(shù)和10 000 km的累計(jì)載荷循環(huán)數(shù)見表5。

表4 樣本截除后的分布參數(shù)

表5 載荷循環(huán)計(jì)數(shù)

根據(jù)分布參數(shù)的估計(jì)結(jié)果，載荷幅值基本上都符合威布爾分布，但形狀參數(shù)都不等于1，用傳統(tǒng)的數(shù)學(xué)方法無(wú)法推導(dǎo)出當(dāng)量載荷的概率密度函數(shù)，所以采用Monte-Carlo法進(jìn)行擬合。

4.3 累計(jì)損傷計(jì)算

4.3.1 材料的S-N曲線擬合

在沒有真實(shí)的S-N曲線可用時(shí)，根據(jù)已有的試驗(yàn)數(shù)據(jù)和素材，對(duì)材料的S-N曲線進(jìn)行估算和擬合，可以作為工程實(shí)踐中進(jìn)行疲勞耐久設(shè)計(jì)的一種近似手段。

S-N曲線的擬合方法主要有兩種：三參數(shù)公式［13，36-38］和分段擬合法［39，43］。

彈簧鋼材料為60Si2Mn，抗拉強(qiáng)度為1 625 MPa［44］，疲勞極限為660 MPa，采用分段擬合法擬合S-N曲線見圖11。

圖11中，S1000和Sbe分別為材料在標(biāo)準(zhǔn)循環(huán)加載103和106次的疲勞強(qiáng)度，Sf為材料的抗拉強(qiáng)度。

對(duì)于鋼材料，S1000取0.72Sf。擬合得到高周疲勞區(qū)疲勞壽命Nfh與名義應(yīng)力幅Sa的關(guān)系式為

圖11 S-N曲線分段擬合

低周疲勞區(qū)疲勞壽命Nfl與名義應(yīng)力幅Sa的關(guān)系式為

4.3.2 蒙特卡羅擬合法累積損傷計(jì)算

蒙特卡羅法擬合螺旋彈簧的分布數(shù)據(jù)見表6。

表6 擬合分布參數(shù)

可得左前、右前、左后和右后4個(gè)螺旋彈簧的當(dāng)量載荷的概率密度函數(shù)f1（S1），f2（S2），f3（S3）和f4（S4）：

根據(jù)式（3）可得

當(dāng)累積疲勞損傷D＝1時(shí)，可得疲勞壽命N為

由式（7）～式（10）4個(gè)螺旋彈簧的零均值一維當(dāng)量隨機(jī)載荷譜的概率密度函數(shù)可得4個(gè)螺旋彈簧的疲勞壽命N1～N4，分別為

采用八點(diǎn)高斯-勒讓德求積公式［45］進(jìn)行積分運(yùn)算，得到4個(gè)螺旋彈簧的疲勞壽命。以N1的計(jì)算為例，計(jì)算方法如下。

將累計(jì)損傷D分為D1和D2兩個(gè)階段，則

首先計(jì)算D1，根據(jù)D1的積分范圍［1 170，1 625］計(jì)算得到8個(gè)積分計(jì)算點(diǎn)，見表7。其中Si為積分計(jì)算點(diǎn)，i＝1，2，…，8。計(jì)算過(guò)程和結(jié)果見表8。

表7 積分點(diǎn)

表8 求積計(jì)算結(jié)果

因此可求得

同樣的方法求得

D2＝4.43738×10-8

由此可算得左前螺旋彈簧的疲勞壽命N1：

同樣的方法求得右前、左后和右后螺旋彈簧的疲勞壽命N2，N3和N4分別為：N2＝9429618，N3＝1020328144，N4＝446309357。

根據(jù)表5中的每10 000 km的載荷循環(huán)數(shù)，可以得到4個(gè)螺旋彈簧的總耐久壽命里程，見表9。

表9 耐久壽命里程

4.3.3 載荷譜分級(jí)法累積損傷計(jì)算［2］

為了進(jìn)行比較，同時(shí)采用載荷譜分級(jí)法對(duì)疲勞壽命進(jìn)行計(jì)算。

載荷的分級(jí)會(huì)影響損傷計(jì)算結(jié)果，研究不同的載荷分級(jí)狀態(tài)下的疲勞損傷結(jié)果表明，當(dāng)載荷分級(jí)大于40級(jí)時(shí)，計(jì)算結(jié)果趨于穩(wěn)定。

采用50×50級(jí)載荷分級(jí)，計(jì)算4個(gè)彈簧的損傷和疲勞壽命，數(shù)據(jù)見表10。

表10 累積損傷和疲勞壽命計(jì)算

4.4 疲勞壽命結(jié)果的分析比較

采用兩種方法計(jì)算的疲勞壽命比較見表11。

表11 總耐久壽命里程

根據(jù)表11的結(jié)果，采用載荷譜分級(jí)法計(jì)算得到的疲勞壽命相對(duì)于采用蒙特卡羅擬合概率密度函數(shù)外推法的計(jì)算結(jié)果要保守，但采用兩種方法的計(jì)算結(jié)果沒有數(shù)量級(jí)的差別，對(duì)于疲勞耐久分析來(lái)說(shuō)誤差在可接受范圍之內(nèi)，驗(yàn)證了總體耐久壽命分析結(jié)果的可信度。由于采用載荷譜分級(jí)法的計(jì)算結(jié)果與載荷譜的分級(jí)狀態(tài)有關(guān)，而采用概率密度函數(shù)法計(jì)算結(jié)果唯一，通過(guò)修正分布參數(shù)可以得到較為準(zhǔn)確的結(jié)果，且易實(shí)現(xiàn)程序化計(jì)算，因此概率密度函數(shù)外推法比采用載荷譜分級(jí)法具有較明顯的優(yōu)勢(shì)。

4.5 誤差原因分析

導(dǎo)致疲勞壽命計(jì)算結(jié)果產(chǎn)生誤差的主要原因有如下幾點(diǎn)：

（1）擬合的S-N曲線與實(shí)際的S-N曲線之間存在誤差；

（2）估計(jì)參數(shù)與實(shí)際分布參數(shù)之間存在誤差；

（3）實(shí)際載荷譜采集時(shí)，取樣數(shù)據(jù)為大量不連續(xù)的載荷，計(jì)算時(shí)將其近似為連續(xù)的載荷譜進(jìn)行積分運(yùn)算，導(dǎo)致計(jì)算結(jié)果產(chǎn)生誤差。

5 結(jié)論

（1）提出了采用蒙特卡羅法生成隨機(jī)序列對(duì)二維隨機(jī)載荷的當(dāng)量載荷的概率密度函數(shù)進(jìn)行擬合的方法，對(duì)不同分布參數(shù)的8 085組數(shù)據(jù)進(jìn)行了概率密度函數(shù)的擬合研究，研究表明采用該方法擬合當(dāng)量載荷的概率密度函數(shù)是可行的。

（2）采集了一款城市SUV在試驗(yàn)場(chǎng)強(qiáng)化路面的載荷譜數(shù)據(jù)，提出了采用樣本截?cái)喾椒ㄟM(jìn)行隨機(jī)載荷譜參數(shù)估計(jì)的方法。參數(shù)檢驗(yàn)結(jié)果表明，在中高載荷段，樣本截?cái)喾梢院芎玫貙?shí)現(xiàn)對(duì)所采載荷譜數(shù)據(jù)的擬合。

（3）提出了采用二維隨機(jī)載荷的等效載荷的概率密度函數(shù)進(jìn)行數(shù)值積分的疲勞壽命的計(jì)算方法。以4個(gè)螺旋彈簧為研究對(duì)象，分別采用載荷譜分級(jí)法和概率密度函數(shù)積分法兩種方法對(duì)所采集的二維隨機(jī)載荷進(jìn)行外推并計(jì)算了累積損傷和疲勞壽命，對(duì)采用兩種方法計(jì)算的結(jié)果進(jìn)行了對(duì)比分析，結(jié)果表明，采用兩種方法的計(jì)算結(jié)果沒有數(shù)量級(jí)的差別，對(duì)于疲勞耐久分析來(lái)說(shuō)誤差在可接受范圍之內(nèi)，互相證實(shí)了總體耐久壽命分析結(jié)果的可信度，也證實(shí)了通過(guò)采用蒙特卡羅法擬合二維隨機(jī)載荷的當(dāng)量載荷的概率密度函數(shù)來(lái)計(jì)算疲勞壽命的方法是可行的，該方法可以很好地解決分布情況較復(fù)雜的二維隨機(jī)載荷的概率密度函數(shù)無(wú)法推導(dǎo)的問(wèn)題。