基于精英保留遺傳算法孔邊疲勞裂紋擴展數據處理*

馬成龍，張 晗，孟麗君

(江漢大學智能制造學院，湖北 武漢 430056)

0 引言

金屬材料在交變應力或在應變作用下會發生疲勞[1]，疲勞損傷是機械損傷的最常見形式，占全部力學破壞的50%～90%[2]。機械裝備一旦發生損傷破壞而導致結構斷裂和故障，將帶來重大的經濟損失和嚴重的社會影響。因此，及時監測機械裝備的裂紋損傷信息并對其剩余壽命進行預測，是保證設備安全運行的重要保障[3]。為此，需要測試材料的疲勞裂紋擴展速率，穩定裂紋擴展段的疲勞裂紋擴展速率常用Paris公式描述[4]，該公式是預測疲勞裂紋擴展壽命的經典公式，常用于長裂紋穩定擴展行為的研究。常見的疲勞裂紋擴展實驗的數據處理方法有割線法和遞增多項式等[5]。本文以2019 PHM Conference Data Challenge挑戰賽公布的試件裂紋長度和疲勞循環次數的實驗數據為基礎[6]，通過python編程，利用遺傳算法對施加了恒定載荷和變載荷的實驗材料穩定階段疲勞裂紋擴展的Paris公式中材料參數進行多目標優化求解，得到了較優的Paris公式參數和a-N曲線擬合結果，依據計算結果可以實現少數據點情況下材料疲勞裂紋的壽命預測。

1 Paris公式

如圖1所示，材料的裂紋擴展速率和應力強度因子的關系可以劃分為三個階段[1-2]，分別為A、B和C三個階段，在材料裂紋擴展的第一階段A內，應力強度因子ΔK和裂紋擴展速率da/dN都相對較低。在階段B內材料應力強度因子變程大于門檻值，為穩定裂紋擴展段，該區域內裂紋擴展速率da/dN與應力強度因子變程符合Paris公式，是目前研究的主要區域。在區域C，材料裂紋快速擴展，da/dN值很大，在實驗研究中通常不考慮。

Paris公式是描述穩定擴展段的疲勞裂紋擴展速率公式，其方程如下：

da/dN=C(ΔK)m

(1)

其中，a是裂紋長度或深度，N是負載循環數；da是裂紋長度增量，dN是負載周期的增量，da/dN是裂紋擴展速率；ΔK是應力強度因子變程；C和m是材料參數。

圖1 裂紋擴展曲線

應力強度因子幅值ΔK可以表示為[1]：

(2)

其中，Y是試樣的幾何因子，對于孔邊疲勞裂紋，其值為[6]：

Y=[1+0.2(1-s)+(1-s)6](2.243-2.64s+1.352s2-0.248s3)

(3)

其中，Δσ為施加的應力范圍，取Δσ=4e7 Pa，s值為：

s=a/(R+a)

(4)

其中，R為孔邊的半徑。

本文中取R=4e-4 mm。

1.1 實驗樣品施加載荷

對于實驗樣品T1至T3，施加的疲勞載荷振幅譜是恒定的；而對于樣本T4則是變化的。兩種類型的載荷振幅圖如圖2。

由圖2可以看出，在對實驗樣品施加恒定載荷譜的情況下，應力范圍可以表示為：

Δσc=σmax-σmin=100.21-4.77=95.44 MPa

(5)

圖2 疲勞載荷振幅譜

1.2 變載荷的等效應力幅

針對變載荷作用下裂紋擴展的高載遲滯和低載加速特征[7]，已有研究人員提出多種Paris 修正公式[8]，可對變幅載荷下的疲勞裂紋進行合理解釋，但是一般的修正公式參數較多，擬合也較困難。對于變幅疲勞計算通常近似按線性疲勞累計損傷Corten-Dolan準則[9]，將變化的載荷幅折算為恒定幅。其表達式為：

(6)

其中，N是應力加載的總循環數；σl是最大應力水平值；Nl是最大交變應力作用下的總循環數；σi是第i級應力水平的應力值；ai是第i級應力的循環數占總循環數的比值；d是實驗確定的常數。

對公式(1)進行積分，可以得到：

(7)

其中，a1，a2為裂紋初始長度和實驗擴展階段結束時的裂紋長度。將Δσ表示應力幅，則有：

(8)

令

(9)

則式(8)變換為：

N=β(Δσve)-m

(10)

將式(10)帶入式(6)可得：

(11)

通過上述公式計算可以將變幅載荷作用下的疲勞裂紋擴展問題，引入等效應力Δσve來表示，以便后續進行計算。

1.3 常用的七點法求材料參數C和m

依據GB/T 6398—2000《金屬材料疲勞裂紋擴展速率實驗方法》[10]，將測得的a，N數據進行擬合求導處理，則Paris公式(1)在雙對數坐標系中可變換為：

lg da/dN=lgC+mlgΔK

(12)

從式(12)中可以得到直線的斜率為m，直線在lg(da/dN)軸上的截距為lgC，最終求得Paris公式中的材料參數C和m。

但是這種方法在擬合的過程中將C和m看做定制，沒有考慮材料本身、測量誤差、環境等因素導致的C和m離散性[11]。研究表明，C和m的分布分別呈對數正態分布和正態分布[12]。

本文針對傳統方法的局限性，提出了基于精英保留遺傳算法進行材料的疲勞裂紋擴展規律預測。在整個過程C和m會在給定的邊界條件中進行編碼處理，使其不再被看成確定的值，服從一定的分布規律，符合實際條件。在算法中設定C的上下邊界為[1e-10，1e-9]，m的上下邊界為[2，4]。

2 數據處理

2.1 實驗數據

對于實驗樣品施加疲勞載荷后得到數據T1、T2、T3和T4，見表1(數據來自：2019 PHM Conference Data Challenge)[13]，其中前三組數據為施加恒定載荷得到的數據，T4為施加變載荷得到的裂紋擴展數據。

表1 實驗樣本數據

2.2 精英保留遺傳算法原理

精英保留遺傳算法(Elitist Strategy Genetic Algorithm，簡稱ESGA)，將群體在進化過程中出現最好的個體不對其進行配對交叉等操作，直接復制到下一代，進而避免了一般遺傳算法中雜交變異過程中較好基因被破壞的可能性。其流程圖見圖3。本文采用精英保留策略遺傳算法對實驗數據進行處理，通過Python的遺傳算法工具箱geatpy，調用增強精英保留多染色體算法模板進行優化求解。

圖3 精英保留遺傳算法流程圖

圖4和圖5分別是程序計算后材料參數C和m的最終值(虛線)和一代種群更新對比圖，其中因得到的材料參數C的值較小，將其轉化為對數分析。從圖中可以看出，算法采用精英保留的策略，將搜索到適應度最高的個體作為精英個體在種群進化過程中得以很好的保留。

圖4 程序優化得到的材料參數m

圖5 程序優化得到的材料參數C的對數

2.3 疲勞裂紋擴展曲線

然后通過編程對實驗數據進行處理分析，圖6給出了計算程序流程圖。

對于表1中各個實驗樣本數據，通過編制程序進行多項式最小二乘擬合得到樣本的疲勞裂紋擴展曲線，恒定載荷樣本結果如圖7、圖8和圖9，變載荷樣本實驗數據結果如圖10。

圖6 Python計算程序框圖

圖7 T1疲勞裂紋擴展曲線

圖8 T2疲勞裂紋擴展曲線

本文選擇樣本可決系數R2反映樣本值擬合情況，表2為四組數據曲線擬合的可決系數值。由表2可看出，四組數據R2值均在0.9以上，擬合效果良好。

圖9 T3 疲勞裂紋擴展曲線

圖10 T4 疲勞裂紋擴展曲線

表2 擬合可決系數

2.4 材料參數

通過Python編寫的程序，采用精英保留的多染色體遺傳算法，在導入實驗數據迭代100次后得到的實驗樣本的材料參數如表3所示。

表3 實驗樣本C和m

對于實驗樣本T4，由于施加的載荷振幅是變化的，所以在對其進行求解時將其用等效應力Δσve來代替，在程序中設定為新的一個決策變量，并設定其上下界為[4e7，9e7]，最后計算得到其載荷循環的替代值為Δσve=4.0002249×107Pa。

下面對計算結果T1的材料參數進行驗證。由文獻[11-12]可知，若把C看成一個確定的值，則lgC也是確定的，這時m服從正態分布。圖11為得到多次結果后m和lgC的關系圖。擬合得到的參數方程為：

m=1.488lgC+15.839

(13)

取lgC的均值為-8.657928，則C值為2.198226e-9，依據方程(13)可求得m的值為2.958647。此時均值對應的Paris公式為：

da/dN=3.510750×10-9(ΔK)2.956003

(14)

圖11 lgC和m關系圖

上述過程將C看作是一個定值，通過文獻[13]可知，這時m服從正態分布。由程序多次運行后得到500個m的值，作出其分布直方圖如圖12所示。從圖中可以看出m近似服從正態分布。

圖12 m分布直方圖

通過SPSS軟件對m進行單樣本柯爾莫戈洛夫-斯米諾夫(Kolmogorov-Smirnov，簡稱K-S)檢驗[14]，結果如表4所示。

由檢驗結果可知，漸進顯著性水平為0.200，大于0.05(小概率事件的概率值)，所以m的分布特性為正態分布。m的均值為2.961929，此時對應的均值Paris公式為：

da/dN=3.510750×10-9(ΔK)2.961929

(15)

表4 單樣本K-S檢驗

式(14)和式(15)對比差異非常小，說明通過精英保留遺傳算法估計的材料參數較為良好。

3 總結

針對傳統的使用Paris公式預測疲勞壽命問題，本文提出了基于遺傳算法對實驗數據進行分析處理，預測其疲勞裂紋擴展規律，得到其材料參數C和m。主要有以下結論：

1)基于精英保留策略的遺傳算法，計算過程中將材料參數在邊界條件中進行編碼處理，由于材料參數本身的離散特性和實驗數據數量的局限性，作為精英個體保留輸出后的材料參數服從一定的分布，且會更加逼近實際值。

2)針對獲得的少量數據點，在變載荷條件下算法仍能夠得到較為精確的材料參數，且隨著數據點的增加結果將更為精確。