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

藤瑞品，宋曉琳，劉國云，曾俊偉

（1.湖南大學，汽車車身先進設計制造國家重點實驗室，長沙 410082； 2.湖南獵豹汽車股份有限公司，長沙 410100）

前言

汽車在工作時受到的載荷譜是一個隨機過程，隨機載荷的表征參數包含載荷幅值和載荷均值，疲勞耐久的研究應基于載荷幅值和均值都為隨機變量的二維隨機載荷進行。

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

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

在進行疲勞分析時，一些相關數據由于各種實際條件限制難以準確得到，而參數擬合是一種可行的近似方法。如黃學偉等［11］提出了一種預測材料與結構裂紋在高周疲勞下疲勞裂紋擴展速率的數值模擬新方法；田秀云等［12］提出了一種金屬材料疲勞裂紋擴展曲線的擬合方法等。

1 二維隨機載荷作用下的疲勞損傷

根據Goodman公式，當量載荷Seq的表達式［13］為

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

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

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

式中：D為材料的累積疲勞損傷；N為材料受到的循環載荷總數量；S為材料受到的當量載荷；f（S）為隨機當量載荷的概率密度函數；Nf為材料在載荷S作用下的疲勞壽命。

2 二維隨機載荷的當量載荷的概率密度函數

2.1 概率密度函數的代數法推導

利用載荷幅值和載荷均值的概率密度以及當量載荷和載荷幅值、載荷均值之間的數學關系可對當量載荷的概率密度函數進行推導，如文獻［18］中推導得到了載荷幅值和載荷均值均符合正態分布的二維隨機載荷的當量載荷的概率密度函數。

2.2 概率密度函數的Monte-Carlo擬合

2.2.1 Monte-Carlo擬合方法和分布參數選擇

采用數學方法推導二維隨機載荷的當量載荷的概率密度函數在實際應用中存在很多局限性：一方面推導過程一般都比較繁瑣；另外在采用Goodman公式求當量載荷的概率密度函數時，很多情況下無法用數學方法直接推導得出。

蒙特卡洛法［19-23］也稱為隨機模擬法，是通過生成隨機數序列（實際上應該為偽隨機數）來對數學問題進行擬合分析的一種數值方法。假設隨機載荷譜的載荷均值和載荷幅值的分布律和概率密度函數都已知，則可通過采用蒙特卡羅法隨機生成載荷幅值和載荷均值數據并得到當量載荷的隨機序列，并分析得到當量載荷的隨機分布特性。

在試驗場對某城市SUV進行載荷譜采集，得到40個通道的應變時間歷程信號，采用雨流計數法［24-28］對載荷譜時間歷程曲線進行計數處理，得到包含不同載荷幅值和載荷均值以及各自所對應的計數頻次的載荷譜。進行參數估計和分布檢驗可得，載荷幅值基本上符合兩參數威布爾分布，載荷均值基本符合正態分布。

設載荷幅值Sa符合二參數威布爾分布，其形狀參數為α，尺度參數為β；載荷均值Sm符合正態分布，分布參數標準差為σ，均值為μ，根據隨機變量的函數的分布特性，載荷均值Sm的函數同樣符合正態分布，分布參數中，標準差為，均值為

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

表1 分布參數范圍

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

2.2.2 隨機載荷序列的生成

采用蒙特卡羅法對表1中分布參數進行隨機生成，得到載荷幅值和載荷均值的隨機序列ui和vi，再采用Goodman公式進行一維隨機載荷的零均值當量載荷的轉化，得到隨機當量載荷序列為fi，則fi的生成方法為

式中vi為符合正態分布的載荷均值的函數所生成的隨機序列。

2.2.3 參數估計

（1）圖形估計［29］

對表1中的不同分布參數的隨機當量載荷序列fi進行WPP（威布爾概率圖）分布檢驗，檢驗結果發現采用威布爾分布形式來對隨機載荷序列fi進行擬合的效果較好，圖1表示一組參數仿真后的圖形法估計參數的結果。

圖1 圖形法估計參數

（2）分布參數估計

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

表2為部分擬合的分布參數數據。

表2 分布參數擬合

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

3 載荷譜數據采集和處理

3.1 載荷譜數據的采集

采用4個車輪六分力傳感器、10個加速度傳感器和38個應變片傳感器，在襄陽試驗場的石塊路面采集了某城市SUV車型整車主要結構載荷譜。

3.2 載荷譜數據的處理［30-34］

原始載荷譜為包含異常點（毛刺）以及連接路面的信號的時間歷程曲線，須將異常數據和連接路面的數據剔除，左前螺旋彈簧經處理后的數據信號如圖2所示。

圖2 左前彈簧載荷譜信號

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

圖3 左前彈簧載荷三維統計直方圖

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

3.3 統計模型的圖形法檢驗［29-35］

圖4 左前彈簧載荷幅值統計直方圖

圖5 左前彈簧載荷均值統計直方圖

為了確定載荷譜數據是否服從威布爾分布和正態分布，可以用相應的概率圖去檢驗。

經過雨流計數法統計右前螺旋彈簧的載荷譜的幅值和均值，進行威布爾和正態分布的圖形檢驗的結果見圖6和圖7。

圖6 威布爾分布檢驗圖

通過觀察發現，圖6中除了部分比較小的幅值和總體趨勢不一致外，其它基本符合威布爾分布的假設。而通過觀察圖7可知，均值非常符合正態分布的假設，這與圖5的直方圖一致。

采用極大似然對模型的分布參數進行估計［14］，部分參數估計的結果見表3。

圖7 正態分布檢驗圖

表3 統計分布參數

其它測試點的載荷分布特性與表中類似，即載荷幅值符合威布爾分布，載荷均值符合正態分布。

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

4.1 載荷譜數據樣本的截除

根據參數估計對載荷譜進行圖形法檢驗時，發現低載荷區域和高載荷區域都發生了較大的偏離，由于高載荷區域對整個疲勞壽命的預測影響較大，因此采用這種方法外推的載荷譜用于疲勞損傷計算會產生較大的誤差。

由于在載荷譜采集時采集了大量的小載荷，絕大部分小載荷低于疲勞極限，因此并不會對疲勞損傷帶來影響，而大量的小載荷會對隨機分布的參數估計造成較大影響，因此采用將小載荷去除的方法將樣本進行截斷，采用截斷后的樣本進行隨機分布的參數估計。

圖8～圖10分別為右前螺旋彈簧載荷數據在未截除狀態、截除100 MPa以下和截除400 MPa以下數據后的估計參數檢驗圖形。

由圖可知，右前螺旋彈簧的載荷數據在截除100 MPa以下的小載荷后，估計參數在中高載荷部分可實現很高的擬合度。

按同樣的方法對表3中數據進行截取后的擬合參數見表4。

圖8 原始數據參數估計

圖9 截除100 MPa以下小載荷

圖10 截除400 MPa以下載荷

對樣本截除后擬合參數進行參數檢驗，在中高載荷部分均可實現很高的擬合度。

4.2 隨機載荷譜的外推

以石塊路工況的4個螺旋彈簧零件的載荷譜數據為研究對象，對隨機載荷譜進行外推。石塊路總長度為2.7 km，耐久試驗總目標里程目標值為10 000 km。4個彈簧載荷在樣本截取后的單次試驗循環里程的載荷樣本數和10 000 km的累計載荷循環數見表5。

表4 樣本截除后的分布參數

表5 載荷循環計數

根據分布參數的估計結果，載荷幅值基本上都符合威布爾分布，但形狀參數都不等于1，用傳統的數學方法無法推導出當量載荷的概率密度函數，所以采用Monte-Carlo法進行擬合。

4.3 累計損傷計算

4.3.1 材料的S-N曲線擬合

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

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

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

圖11中，S1000和Sbe分別為材料在標準循環加載103和106次的疲勞強度，Sf為材料的抗拉強度。

對于鋼材料，S1000取0.72Sf。擬合得到高周疲勞區疲勞壽命Nfh與名義應力幅Sa的關系式為

圖11 S-N曲線分段擬合

低周疲勞區疲勞壽命Nfl與名義應力幅Sa的關系式為

4.3.2 蒙特卡羅擬合法累積損傷計算

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

表6 擬合分布參數

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

根據式（3）可得

當累積疲勞損傷D＝1時，可得疲勞壽命N為

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

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

將累計損傷D分為D1和D2兩個階段，則

首先計算D1，根據D1的積分范圍［1 170，1 625］計算得到8個積分計算點，見表7。其中Si為積分計算點，i＝1，2，…，8。計算過程和結果見表8。

表7 積分點

表8 求積計算結果

因此可求得

同樣的方法求得

D2＝4.43738×10-8

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

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

根據表5中的每10 000 km的載荷循環數，可以得到4個螺旋彈簧的總耐久壽命里程，見表9。

表9 耐久壽命里程

4.3.3 載荷譜分級法累積損傷計算［2］

為了進行比較，同時采用載荷譜分級法對疲勞壽命進行計算。

載荷的分級會影響損傷計算結果，研究不同的載荷分級狀態下的疲勞損傷結果表明，當載荷分級大于40級時，計算結果趨于穩定。

采用50×50級載荷分級，計算4個彈簧的損傷和疲勞壽命，數據見表10。

表10 累積損傷和疲勞壽命計算

4.4 疲勞壽命結果的分析比較

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

表11 總耐久壽命里程

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

4.5 誤差原因分析

導致疲勞壽命計算結果產生誤差的主要原因有如下幾點：

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

（2）估計參數與實際分布參數之間存在誤差；

（3）實際載荷譜采集時，取樣數據為大量不連續的載荷，計算時將其近似為連續的載荷譜進行積分運算，導致計算結果產生誤差。

5 結論

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

（2）采集了一款城市SUV在試驗場強化路面的載荷譜數據，提出了采用樣本截斷方法進行隨機載荷譜參數估計的方法。參數檢驗結果表明，在中高載荷段，樣本截斷法可以很好地實現對所采載荷譜數據的擬合。

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