基于精英保留遺傳算法孔邊疲勞裂紋擴(kuò)展數(shù)據(jù)處理*

馬成龍，張 晗，孟麗君

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

0 引言

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

1 Paris公式

如圖1所示，材料的裂紋擴(kuò)展速率和應(yīng)力強(qiáng)度因子的關(guān)系可以劃分為三個(gè)階段[1-2]，分別為A、B和C三個(gè)階段，在材料裂紋擴(kuò)展的第一階段A內(nèi)，應(yīng)力強(qiáng)度因子ΔK和裂紋擴(kuò)展速率da/dN都相對(duì)較低。在階段B內(nèi)材料應(yīng)力強(qiáng)度因子變程大于門檻值，為穩(wěn)定裂紋擴(kuò)展段，該區(qū)域內(nèi)裂紋擴(kuò)展速率da/dN與應(yīng)力強(qiáng)度因子變程符合Paris公式，是目前研究的主要區(qū)域。在區(qū)域C，材料裂紋快速擴(kuò)展，da/dN值很大，在實(shí)驗(yàn)研究中通常不考慮。

Paris公式是描述穩(wěn)定擴(kuò)展段的疲勞裂紋擴(kuò)展速率公式，其方程如下：

da/dN=C(ΔK)m

(1)

其中，a是裂紋長度或深度，N是負(fù)載循環(huán)數(shù)；da是裂紋長度增量，dN是負(fù)載周期的增量，da/dN是裂紋擴(kuò)展速率；ΔK是應(yīng)力強(qiáng)度因子變程；C和m是材料參數(shù)。

圖1 裂紋擴(kuò)展曲線

應(yīng)力強(qiáng)度因子幅值ΔK可以表示為[1]：

(2)

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

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

(3)

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

s=a/(R+a)

(4)

其中，R為孔邊的半徑。

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

1.1 實(shí)驗(yàn)樣品施加載荷

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

由圖2可以看出，在對(duì)實(shí)驗(yàn)樣品施加恒定載荷譜的情況下，應(yīng)力范圍可以表示為：

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

(5)

圖2 疲勞載荷振幅譜

1.2 變載荷的等效應(yīng)力幅

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

(6)

其中，N是應(yīng)力加載的總循環(huán)數(shù)；σl是最大應(yīng)力水平值；Nl是最大交變應(yīng)力作用下的總循環(huán)數(shù)；σi是第i級(jí)應(yīng)力水平的應(yīng)力值；ai是第i級(jí)應(yīng)力的循環(huán)數(shù)占總循環(huán)數(shù)的比值；d是實(shí)驗(yàn)確定的常數(shù)。

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

(7)

其中，a1，a2為裂紋初始長度和實(shí)驗(yàn)擴(kuò)展階段結(jié)束時(shí)的裂紋長度。將Δσ表示應(yīng)力幅，則有：

(8)

令

(9)

則式(8)變換為：

N=β(Δσve)-m

(10)

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

(11)

通過上述公式計(jì)算可以將變幅載荷作用下的疲勞裂紋擴(kuò)展問題，引入等效應(yīng)力Δσve來表示，以便后續(xù)進(jìn)行計(jì)算。

1.3 常用的七點(diǎn)法求材料參數(shù)C和m

依據(jù)GB/T 6398—2000《金屬材料疲勞裂紋擴(kuò)展速率實(shí)驗(yàn)方法》[10]，將測(cè)得的a，N數(shù)據(jù)進(jìn)行擬合求導(dǎo)處理，則Paris公式(1)在雙對(duì)數(shù)坐標(biāo)系中可變換為：

lg da/dN=lgC+mlgΔK

(12)

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

但是這種方法在擬合的過程中將C和m看做定制，沒有考慮材料本身、測(cè)量誤差、環(huán)境等因素導(dǎo)致的C和m離散性[11]。研究表明，C和m的分布分別呈對(duì)數(shù)正態(tài)分布和正態(tài)分布[12]。

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

2 數(shù)據(jù)處理

2.1 實(shí)驗(yàn)數(shù)據(jù)

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

表1 實(shí)驗(yàn)樣本數(shù)據(jù)

2.2 精英保留遺傳算法原理

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

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

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

圖4 程序優(yōu)化得到的材料參數(shù)m

圖5 程序優(yōu)化得到的材料參數(shù)C的對(duì)數(shù)

2.3 疲勞裂紋擴(kuò)展曲線

然后通過編程對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行處理分析，圖6給出了計(jì)算程序流程圖。

對(duì)于表1中各個(gè)實(shí)驗(yàn)樣本數(shù)據(jù)，通過編制程序進(jìn)行多項(xiàng)式最小二乘擬合得到樣本的疲勞裂紋擴(kuò)展曲線，恒定載荷樣本結(jié)果如圖7、圖8和圖9，變載荷樣本實(shí)驗(yàn)數(shù)據(jù)結(jié)果如圖10。

圖6 Python計(jì)算程序框圖

圖7 T1疲勞裂紋擴(kuò)展曲線

圖8 T2疲勞裂紋擴(kuò)展曲線

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

圖9 T3 疲勞裂紋擴(kuò)展曲線

圖10 T4 疲勞裂紋擴(kuò)展曲線

表2 擬合可決系數(shù)

2.4 材料參數(shù)

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

表3 實(shí)驗(yàn)樣本C和m

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

下面對(duì)計(jì)算結(jié)果T1的材料參數(shù)進(jìn)行驗(yàn)證。由文獻(xiàn)[11-12]可知，若把C看成一個(gè)確定的值，則lgC也是確定的，這時(shí)m服從正態(tài)分布。圖11為得到多次結(jié)果后m和lgC的關(guān)系圖。擬合得到的參數(shù)方程為：

m=1.488lgC+15.839

(13)

取lgC的均值為-8.657928，則C值為2.198226e-9，依據(jù)方程(13)可求得m的值為2.958647。此時(shí)均值對(duì)應(yīng)的Paris公式為：

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

(14)

圖11 lgC和m關(guān)系圖

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

圖12 m分布直方圖

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

由檢驗(yàn)結(jié)果可知，漸進(jìn)顯著性水平為0.200，大于0.05(小概率事件的概率值)，所以m的分布特性為正態(tài)分布。m的均值為2.961929，此時(shí)對(duì)應(yīng)的均值Paris公式為：

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

(15)

表4 單樣本K-S檢驗(yàn)

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

3 總結(jié)

針對(duì)傳統(tǒng)的使用Paris公式預(yù)測(cè)疲勞壽命問題，本文提出了基于遺傳算法對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行分析處理，預(yù)測(cè)其疲勞裂紋擴(kuò)展規(guī)律，得到其材料參數(shù)C和m。主要有以下結(jié)論：

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

2)針對(duì)獲得的少量數(shù)據(jù)點(diǎn)，在變載荷條件下算法仍能夠得到較為精確的材料參數(shù)，且隨著數(shù)據(jù)點(diǎn)的增加結(jié)果將更為精確。