基于比例型Paris 公式和逆高斯過程的金屬疲勞裂紋擴展隨機模型

陳 龍，黃天立，周 浩

(中南大學(xué)土木工程學(xué)院，湖南，長沙 410075)

金屬構(gòu)件的疲勞可靠性及剩余壽命預(yù)測一直是抗疲勞設(shè)計的主要任務(wù)之一，而疲勞裂紋擴展分析作為金屬構(gòu)件疲勞壽命預(yù)測的主要內(nèi)容，其相關(guān)問題一直是該領(lǐng)域的研究熱點。由于外部的不確定因素和材料內(nèi)部的非均質(zhì)等影響，疲勞裂紋擴展過程中的波動范圍通常較大，如著名的Virkler 金屬疲勞裂紋擴展試驗數(shù)據(jù)[1]等。探討這種隨機性的原因及規(guī)律，對于橋梁、飛機、船舶等安全性和可靠性要求較高的構(gòu)件設(shè)計和壽命預(yù)測是十分重要的，也是必需的。

大量研究表明[2-5]，金屬疲勞裂紋擴展過程中的隨機性主要包括兩種：1)材料內(nèi)部阻力的不恒定、裂紋面的不規(guī)則導(dǎo)致裂紋擴展路徑的不規(guī)則波動，處理此類隨機性的通常做法是，在確定性疲勞裂紋擴展速率模型的基礎(chǔ)上引入隨機因子，借助數(shù)學(xué)理論和方法對概率壽命進行估計，形成了基于隨機過程的疲勞裂紋擴展模型，蔣鳴曉等[2]、吳圣川等[3]對這類模型進行了綜述，其中包括B-model、連續(xù)的馬爾科夫過程模型、相關(guān)過程模型等，楊冰等[4]在Elber 方程的基礎(chǔ)上，假設(shè)疲勞裂紋擴展速率服從一個平穩(wěn)的對數(shù)正態(tài)分布，建立了相應(yīng)的裂紋擴展隨機模型，劉勇等[5]基于Paris公式，利用概率密度演化方法，解釋了疲勞裂紋擴展過程中的波動行為，并對構(gòu)件的疲勞壽命進行預(yù)測；2)材料微觀結(jié)構(gòu)、缺陷分布的不均勻性導(dǎo)致了樣本間裂紋擴展速率的差異，處理這種隨機性的通常方法是，在第一類模型基礎(chǔ)上，隨機化代表材料固有屬性的參數(shù)(如Paris 公式中的材料參數(shù)C和m)，從而引入隨機效應(yīng)來表征樣本之間的異質(zhì)性。

單個樣本的疲勞裂紋擴展路徑在宏觀上表現(xiàn)為循環(huán)荷載下的間接性增長，這類由一系列微小增量累積形成的損傷很適合采用獨立增量過程模擬[6-7]。近年來，維納(WP)、伽馬(GP)和逆高斯(IGP)等具有獨立增量性質(zhì)的隨機過程，由于其在退化建模方面表現(xiàn)出的優(yōu)良特性，被廣泛用于疲勞裂紋擴展的隨機建模。Wang 等[8]采用冪函數(shù)作為退化軌道函數(shù)，聯(lián)合維納過程，建立了疲勞裂紋擴展隨機效應(yīng)模型，并推導(dǎo)了剩余壽命概率密度函數(shù)的近似解析解。Ye 等[9]利用指數(shù)函數(shù)描述疲勞裂紋的平均擴展路徑，結(jié)合伽馬過程，建立了相應(yīng)的簡單模型和隨機效應(yīng)模型，并采用了半?yún)?shù)方法估計模型參數(shù)。Peng 等[10]利用逆高斯過程描述疲勞裂紋的隨機擴展，同時比較了冪函數(shù)和指數(shù)函數(shù)兩種退化軌道下的模型擬合效果，結(jié)果表明，指數(shù)型退化軌道下的隨機模型能更好地解釋裂紋擴展的隨機行為。Hao 等[11]提出了一個具有偏正態(tài)隨機效應(yīng)的擴展逆高斯過程模型，并將其應(yīng)用到疲勞裂紋的剩余壽命預(yù)測，結(jié)果表明，相對于傳統(tǒng)的逆高斯隨機效應(yīng)模型，其所提出的擴展模型擬合效果更好。

雖然采用獨立增量過程描述裂紋擴展行為的研究較多，但仍存在以下兩個問題：1)描述裂紋平均擴展路徑的軌道函數(shù)僅采用兩種形式，即冪函數(shù)和指數(shù)函數(shù)，由此建立的擴展模型實質(zhì)上是一個統(tǒng)計模型，無法從退化機理層面描述疲勞裂紋的擴展過程；2)上述研究中都采用“時間間隔[t,t+δt] 內(nèi)裂紋長度增量 δa與當(dāng)前狀態(tài)a(t)無關(guān)，僅取決于當(dāng)前t和時間間隔 δt”的基本假設(shè)。實際上，疲勞裂紋擴展過程中不規(guī)則波動的本質(zhì)是，由于靠近裂紋尖端的材料性質(zhì)的局部變化。眾多關(guān)于疲勞裂紋擴展的研究表明，裂紋擴展速率da/dt是一個與當(dāng)前狀態(tài)a有關(guān)的函數(shù)，裂紋長度增量 δa本質(zhì)上是依賴當(dāng)前狀態(tài)a的。因此，為了更好地解釋金屬疲勞裂紋擴展中的隨機性，有必要建立一個更貼近裂紋擴展機理的退化模型。

本文針對金屬疲勞裂紋擴展特性，以“首次達到給定裂紋長度a的時間t(a)”作為隨機過程的描述量，采用比例型Paris 公式作為退化軌道函數(shù)描述裂紋的平均擴展軌跡，建立基于逆高斯過程的單體樣本疲勞裂紋擴展隨機模型和考慮樣本異質(zhì)性的裂紋擴展隨機效應(yīng)模型，分別采用最大似然估計(MLE)和期望最大算法(EM)推導(dǎo)了單樣本模型和隨機效應(yīng)模型的參數(shù)估計公式。最后，利用所提出的模型來擬合Virkler 金屬疲勞裂紋擴展試驗數(shù)據(jù)，并進行擬合優(yōu)度分析。

1 疲勞裂紋擴展隨機模型

1.1 比例型Paris 公式

迄今為止，在斷裂力學(xué)基礎(chǔ)上發(fā)展出來的疲勞裂紋擴展模型較多，其中Paris-Erdogan 公式[12]因其形式簡單和實用，得到了工程界的普遍應(yīng)用。其表達式如式(1)所示：

式中： da(t)/dt為裂紋擴展的速率；a為疲勞裂紋的長度；t為試件經(jīng)受疲勞荷載的加載時間，通常以循環(huán)加載的次數(shù)表示；參數(shù)C′和m′代表材料參數(shù)，其值取決試件的材料和幾何尺寸； ΔK(a)為應(yīng)力強度因子，在實驗條件和試件尺寸確定的情況下； ΔK(a) 為一個關(guān)于裂紋長度a的函數(shù)。

圖1 給出了經(jīng)典的Virkler[1]疲勞裂紋擴展數(shù)據(jù)。該試驗樣本采用68 塊應(yīng)用于飛機結(jié)構(gòu)的2024-T3 鋁合金板件，試件形狀呈矩形，長度、寬度和厚度分別為558.8 mm、152.4 mm、2.54 mm。金屬板在試驗時受到頻率20 Hz，應(yīng)力幅為24.14 MPa，應(yīng)力比為0.2 的正弦循環(huán)荷載。試驗之前在樣本中心形成長度為 2a0=18 mm的初始裂紋長度，裂紋在疲勞荷載作用下向兩邊擴展延伸，同時只記錄一邊的裂紋長度數(shù)據(jù)，當(dāng)疲勞裂紋長度達到49.8 mm時，試驗結(jié)束。

圖1 Virkler 疲勞試驗數(shù)據(jù)Fig. 1 Experimental dataset from Virkler fatigue test

圖2 l og10C′ 和 m′的相關(guān)性Fig. 2 Correlation between l og10C′ andm′

1.2 單樣本疲勞裂紋擴展隨機模型

由于疲勞裂紋長度在宏觀上隨時間的推移不斷增長，因此，以往的研究多將裂紋長度a(t)視為隨時間t變化的一種隨機過程，記做 {A(t),t＞0}，并假設(shè)其過程是一個獨立增量過程，由此來模擬疲勞裂紋的間接性增長。然而實際上，疲勞裂紋長度的增量 δa之間是相互依賴的，這種建模方式與裂紋擴展的本質(zhì)是相違背的。Guida 等[14]認為，研究疲勞裂紋隨機擴展問題的最佳方法是，將“達到給定裂紋長度a所需要的時間t(a)”視為一個隨機過程，記為 {T(a),a≥a0}。若隨機過程{A(t),t＞0}是一個狀態(tài)依賴過程，則相應(yīng)的“逆過程” {T(a),a≥a0} 的時間增量 δt是相互獨立的。反之，若 {T(a),a≥a0}被建模成獨立增量過程，則描述狀態(tài)隨時間變化的隨機過程 {A(t),t＞0}應(yīng)該是一個純狀態(tài)依賴的過程[15]。

本文假設(shè) {T(a),a≥a0}服從一個非齊次的逆高斯過程(IGP)，由于傳統(tǒng)的IGP[16]描述的是狀態(tài)隨時間的隨機變化的過程，因此，按照文獻[15]對伽馬過程的一般定義修改方法，本文對IGP 的一般定義進行修改，

1.3 考慮樣本差異的疲勞裂紋擴展隨機效應(yīng)模型

在實際工程中，由于工藝、設(shè)計以及試驗環(huán)境等因素的差異，同類樣本不同個體之間的裂紋擴展過程會產(chǎn)生較大的差異，這種差異性在退化數(shù)據(jù)上體現(xiàn)為裂紋擴展路徑的分散性，因此有必要對其進行描述。假設(shè)單路徑隨機擴展模型中的參數(shù)C為一個隨機變量，將個體特有的隨機效應(yīng)集成到退化過程中以表征個體差異性，使退化模型能夠同時描述裂紋擴展過程中的兩種隨機性。

為了保證參數(shù)C的非負性以及數(shù)學(xué)推導(dǎo)上可行性，Wang 等[17]建議隨機參數(shù)應(yīng)服從一個截斷正態(tài)分布TN(ω,κ-2)κ＞0，概率密度函數(shù)如下所示：

首達時間增量 ΔT(δa)的邊緣概率密度函數(shù)以及邊緣累積分布函數(shù)可分別通過下列積分計算：

式中，BvN(·)表示二元標(biāo)準正態(tài)累積分布函數(shù)，其余參數(shù)如下所示：

2 隨機模型參數(shù)估計

2.1 單樣本隨機模型參數(shù)估計

2.2 考慮樣本差異的隨機效應(yīng)模型參數(shù)估計

3 模型擬合優(yōu)度檢驗

利用Kolmogorov-Smirnov 統(tǒng)計量D 衡量經(jīng)驗分布和理論分布之間的差異，同時獲取假設(shè)檢驗的P 值，對于簡單假設(shè)(25)：

1) 若P 值小于等于α ，則在顯著水平α 下拒絕H0；

2) 若P 值大于 α，則在顯著水平 α 下接受H0。

4 疲勞裂紋隨機模型擬合優(yōu)度分析

4.1 單樣本疲勞裂紋擴展模型擬合優(yōu)度分析

圖3 材料參數(shù)的常規(guī)ML 估計值和去相關(guān)ML 估計值Fig. 3 conventional and decorrelated ML estimates of material parameters

圖4 68 條疲勞路徑的K-S 檢驗P 值Fig. 4 P-value of 68 fatigue crack growth paths by using Kolmogorov-Smirnov test

4.2 考慮樣本差異的疲勞裂紋隨機效應(yīng)模型擬合優(yōu)度分析

在實際試驗過程中，同類產(chǎn)品不同個體之間在裂紋擴展路徑上也會表現(xiàn)出較大的差異，通過隨機化單體擴展模型中的參數(shù)C，令其服從一個截斷正態(tài)分布，可以建立體現(xiàn)個體差異的隨機效應(yīng)模型。表1 給出了采用EM 算法對隨機效應(yīng)模型的參數(shù)結(jié)果。

表1 隨機效應(yīng)模型參數(shù)估計值Table 1 Parameter estimates of random effect model

在表1 模型估計的基礎(chǔ)上，本節(jié)對隨機效應(yīng)模型一些關(guān)鍵量進行了統(tǒng)計擬合分析。圖5 給出了均值函數(shù)和方差函數(shù)ML 估計和樣本估計的對比圖。當(dāng)裂紋長度較小時(a≤20 mm)，均值和方差的ML 估計值和樣本估計值相差較大，隨機效應(yīng)模型無法較好地預(yù)測樣本值；當(dāng)疲勞裂紋長度較大(a≥20 mm )時，均值E[T(a)]的ML 估計值和樣本估計值之間的誤差為±1.4%，模型和樣本的擬合程度較好，而方差 Var[T(a)]的擬合程度相對于均值E[T(a)]較低。

圖5 首達時間 T(a)的均值和方差的ML 估計以及樣本估計Fig. 5 Empirical and ML estimates of the mean and variance of first timeT(a)

圖7(a)給出了首達時間T(a)在10%、40%、70%、90%時的百分位數(shù)ML 估計和樣本估計的比較圖；圖7(b)給出了ML 估計和樣本估計的相對誤差圖，可以看到，裂紋擴展前期，模型估計值和樣本值相差較大，擬合程度并不高；當(dāng)裂紋長度a≥20 mm，相對誤差基本上穩(wěn)定在±2.5%左右，模型的擬合效果較好。

圖6(a)給出了首次到達給定裂紋長度a=20 mm,25 mm, 30 mm, 35 mm, 40 mm 的時間T(a)的累積分布函數(shù)的ML 估計和樣本估計的對比圖；圖6(b)給出了時間t=0.8 ×104cycle,1.2 ×104cycle,1.6×104cycle,2.0 ×104cycle,2.2 ×104cycle時，疲勞裂紋的擴展長度A(t)的累積分布函數(shù)的ML 估計和經(jīng)驗估計的對比圖。從圖6 可以看出，隨機變量T(a)和A(t)的CDF 的ML 估計曲線和樣本估計值的擬合程度較好。

圖6 首達時間 T(a) 和 裂紋長度 A(t)的CDF 的ML 估計以及樣本估計Fig. 6 Sample and ML estimates of the CDF of first passage time T (a) and crack lengthA(t)

圖7 首達時間 T (a)的百分位數(shù)ML 估計與樣本估計對比圖以及相對誤差圖Fig. 7 Sample and ML estimates of percentiles of first passage time T(a) and relative percent differences

為對隨機參數(shù)C的截斷正態(tài)分布假設(shè)進行檢驗，圖8 給出了隨機參數(shù)C的Q-Q 圖。從圖8 可以看出，樣本點基本分布在對角線附近。通過K-S檢驗得到的假設(shè)檢驗P 值為0.777，大于0.05 的顯著性水平 α。因此，可以認為隨機參數(shù)C服從一個截斷正態(tài)分布。

圖8 隨機變量的C 的截斷正態(tài)分布Q-Q 圖Fig. 8 The normal Q-Q plot under normal distribution hypothesis for C

5 結(jié)論

針對金屬疲勞裂紋擴展過程中的不確定性，采用比例型Pairs 疲勞公式，建立了基于逆高斯過程的單樣本疲勞裂紋擴展隨機模型和考慮樣本差異的裂紋擴展隨機效應(yīng)模型，分別采用最大似然估計和EM 算法推導(dǎo)了單樣本模型和隨機效應(yīng)模型的參數(shù)估計公式。最后將兩種模型分別應(yīng)用于Virkler 疲勞試驗數(shù)據(jù)的擬合，并對擬合效果進行檢驗，可得出如下結(jié)論：

1)采用比例型Paris 公式的隨機模型能夠有效消除參數(shù)之間的強相關(guān)性。

2)通過概率積分變換判斷疲勞裂紋增量是否符合逆高斯分布族，衡量了單樣本疲勞裂紋擴展隨機模型的擬合效果。結(jié)果表明，沒有明顯的證據(jù)顯示單個疲勞裂紋的擴展過程不服從一個逆高斯過程。

3)利用圖形法對一些關(guān)鍵量進行擬合分析，判斷了考慮樣本差異的疲勞裂紋擴展隨機效應(yīng)模型整體的擬合效果。結(jié)果表明，所選擇關(guān)鍵量的ML 估計值都能夠較好地擬合樣本估計值。

由此表明，基于比例型Pairs 公式和逆高斯過程的隨機模型能夠有效地分析和解釋金屬疲勞裂紋擴展過程中的不確定性。