基于代理模型的廣布疲勞損傷平均行為快速確定方法

潘 微，鄭旭東，李 杰，李 彥，談嘉好，陳 昶

（航空工業(yè)成都飛機(jī)工業(yè)（集團(tuán)）有限責(zé)任公司，成都 610092）

在不受外界影響的情況下，機(jī)隊(duì)中50%飛機(jī)在所研究的廣布疲勞損傷（WFD）敏感結(jié)構(gòu)發(fā)生廣布疲勞損傷的時(shí)間點(diǎn)定義為廣布疲勞損傷平均行為[1]。其中，廣布疲勞損傷是指在結(jié)構(gòu)相似、應(yīng)力相當(dāng)?shù)亩鄠€(gè)細(xì)節(jié)部位中同時(shí)存在足夠尺度和密度的裂紋，隨著裂紋的增長(zhǎng)，裂紋間的相互影響加劇，從而加快降低結(jié)構(gòu)的剩余強(qiáng)度，使結(jié)構(gòu)不再滿(mǎn)足剩余強(qiáng)度要求[2]。基于斷裂力學(xué)的確定性分析在評(píng)估工程設(shè)計(jì)和安全性、預(yù)測(cè)結(jié)構(gòu)壽命、對(duì)不同類(lèi)型缺陷的壽命評(píng)估等方面一直起著重要的作用，但是由于結(jié)構(gòu)中出現(xiàn)孔洞、位錯(cuò)、夾雜等隨機(jī)缺陷導(dǎo)致了裂紋萌生及擴(kuò)展的不確定因素，使疲勞分散性成為材料的固有屬性，因此基于疲勞統(tǒng)計(jì)學(xué)等理論建立裂紋隨機(jī)萌生和擴(kuò)展的威布爾分布模型，進(jìn)而對(duì)預(yù)測(cè)多裂紋結(jié)構(gòu)壽命有著重要的意義。另一方面，確定廣布疲勞損傷平均行為求解復(fù)雜，需要消耗大量時(shí)間，因此基于代理模型的快速求解方法對(duì)廣布疲勞損傷平均行為工程評(píng)估有重要意義。

在確定廣布疲勞損傷的概率統(tǒng)計(jì)特性研究中，閆曉中等[3]、雷曉欣等[4]針對(duì)多部位損傷（MSD）結(jié)構(gòu)，通過(guò)蒙特卡洛法模擬，提出多裂紋疲勞損傷隨機(jī)模型，實(shí)現(xiàn)了多裂紋隨機(jī)萌生和擴(kuò)展及廣布疲勞損傷發(fā)生的過(guò)程。Ai 等[5]提出基于缺陷的大小和位置來(lái)說(shuō)明制造缺陷對(duì)疲勞壽命的影響的概率模型，通過(guò)對(duì)具有威布爾分布的表面初始裂紋的損傷機(jī)理進(jìn)行建模，得出考慮缺陷表面位置影響的校正因子。

本文在已有的研究基礎(chǔ)上，建立裂紋萌生和擴(kuò)展模型，利用ABAQUS 二次開(kāi)發(fā)得到的應(yīng)力強(qiáng)度因子數(shù)值，根據(jù)蒙特卡羅法預(yù)測(cè)廣布疲勞損傷平均行為，在保證精度的情況下預(yù)測(cè)結(jié)構(gòu)發(fā)生概率疲勞時(shí)的全壽命。最后，通過(guò)機(jī)器學(xué)習(xí)建立廣布疲勞損傷仿真程序訓(xùn)練集，使用多項(xiàng)式響應(yīng)面和高斯過(guò)程回歸的代理模型進(jìn)行訓(xùn)練，可以快速確定結(jié)構(gòu)發(fā)生廣布疲勞損傷時(shí)的平均行為。

1 概率疲勞分析方法

1.1 裂紋隨機(jī)萌生模型

由于材料本身的不均勻性、不同加工方式、試件大小、試驗(yàn)誤差及試驗(yàn)環(huán)境等因素，會(huì)引起結(jié)構(gòu)疲勞壽命的分散，從而在多孔結(jié)構(gòu)中孔邊裂紋萌生位置及次序具有隨機(jī)性，在含有多個(gè)相似細(xì)節(jié)結(jié)構(gòu)中，第i個(gè)細(xì)節(jié)的裂紋萌生壽命Ni（0），i=1，2，…，M由下式給出

式中：Ci、m為材料常數(shù)。第i個(gè)危險(xiǎn)細(xì)節(jié)的隨機(jī)參數(shù)服從對(duì)數(shù)正態(tài)分布[6]，即

式中：C0i為服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)參數(shù)，下角標(biāo)0為Ci的初始值和SlgCi分別為Ci的對(duì)數(shù)均值和對(duì)數(shù)標(biāo)準(zhǔn)差。

1.2 裂紋隨機(jī)擴(kuò)展模型

Paris 裂紋擴(kuò)展公式計(jì)算過(guò)程簡(jiǎn)便，便于運(yùn)用于工程上的具體運(yùn)算，Paris 裂紋擴(kuò)展公式為

式中：G、n為裂紋體的材料特性常數(shù)；△K為應(yīng)力強(qiáng)度因子幅值，即△K=Kmax-Kmin。

裂紋擴(kuò)展的分散性表現(xiàn)為將材料特性參數(shù)G用隨機(jī)變量Gi表示，即

三參數(shù)正態(tài)分布可以更好地表示隨機(jī)參數(shù)的可變性，引入變量ui，其表達(dá)式為

ui的概率密度函數(shù)滿(mǎn)足三參數(shù)對(duì)數(shù)正態(tài)分布[7]，如下式

式中：τ 為位置參數(shù)；α 為分散參數(shù)；β 為形狀參數(shù)，通過(guò)試驗(yàn)數(shù)據(jù)可以得到。根據(jù)三參數(shù)正態(tài)分布的性質(zhì)，Gi可以通過(guò)下式得到，即

式中：G0i為服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù)，由隨機(jī)發(fā)生器產(chǎn)生。

2 廣布疲勞損傷平均行為分析實(shí)例

結(jié)構(gòu)的全壽命Nf為裂紋萌生壽命N和裂紋擴(kuò)展壽命Nc之和，即

由第1 節(jié)提出的裂紋概率萌生模型，應(yīng)用蒙特卡羅法，得到結(jié)構(gòu)的裂紋萌生壽命；通過(guò)ABAQUS 建立七孔板裂紋尖端應(yīng)力強(qiáng)度因子計(jì)算模型，七孔板的有限元模型如圖1 所示。

圖1 七孔平板FEM 模型

通過(guò)ABAQUS 二次開(kāi)發(fā)獲得裂紋擴(kuò)展每個(gè)增量步的裂紋尖端應(yīng)力強(qiáng)度因子，建立裂紋擴(kuò)展的概率模型，最后得到整個(gè)結(jié)構(gòu)的全壽命。具體步驟如下。

第一步，由蒙特卡羅法得到M個(gè)孔邊危險(xiǎn)細(xì)節(jié)的萌生壽命Ni（i=1，2，…，M-1，M）。

第二步，按照裂紋萌生順序依次給每個(gè)危險(xiǎn)細(xì)節(jié)以最小可檢裂紋尺寸為1.27 mm 的初始裂紋。

第三步，由組合法求解應(yīng)力強(qiáng)度因子，計(jì)算復(fù)雜邊界條件下每個(gè)孔邊裂紋的應(yīng)力強(qiáng)度因子。

第四步，根據(jù)裂紋擴(kuò)展隨機(jī)模型計(jì)算出裂紋擴(kuò)展增量，由塑性連通準(zhǔn)則和Kc斷裂準(zhǔn)則共同判斷結(jié)構(gòu)是否失效，若沒(méi)有失效，則繼續(xù)重復(fù)以上計(jì)算；若失效，本次模擬失效。

為驗(yàn)證方法的準(zhǔn)確性，本文以七孔板為例進(jìn)行壽命預(yù)測(cè)，根據(jù)文獻(xiàn)[8]試驗(yàn)數(shù)據(jù)，概率模擬過(guò)程中的疲勞載荷最大值分別為Pmax=82 kN，Pmax=90 kN，Pmax=98 kN，應(yīng)力比均為0.1，材料為L(zhǎng)Y12-CZ 鋁合金，材料參數(shù)見(jiàn)表1，裂紋隨機(jī)萌生、擴(kuò)展參數(shù)見(jiàn)表2 和表3。

表1 LY12-CZ 鋁合金材料參數(shù)

表2 LY12-CZ 鋁合金裂紋隨機(jī)萌生參數(shù)

表3 LY12-CZ 鋁合金裂紋隨機(jī)擴(kuò)展參數(shù)

在3 種應(yīng)力水平下，分別進(jìn)行100 次蒙特卡羅模擬抽樣，模擬結(jié)果分別如圖2—圖4 所示。

圖2 Pmax=82 kN 壽命累計(jì)分布曲線

圖3 Pmax=90 kN 壽命累計(jì)分布曲線

圖4 Pmax=98 kN 壽命累計(jì)分布曲線

從圖2—圖4 中可以看出，在不同載荷作用下七孔板結(jié)構(gòu)的WFD 失效時(shí)刻累計(jì)分布曲線，并且可以得到50%概率下的廣布疲勞損傷敏感結(jié)構(gòu)疲勞循環(huán)數(shù)，即廣布疲勞損傷平均行為分別為Nave，82=230 450 次、Nave，90=164 665 次、Nave，98=121 700 次。平均行為計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比見(jiàn)表4，其中為對(duì)數(shù)壽命的平均值為試驗(yàn)對(duì)數(shù)壽命的平均值；SlgNc為對(duì)數(shù)壽命的標(biāo)準(zhǔn)差，SlgNt為試驗(yàn)對(duì)數(shù)壽命的標(biāo)準(zhǔn)差。

表4 七孔板計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比

從表4 中可以看出，MSD 結(jié)構(gòu)廣布疲勞平均行為模擬結(jié)果與試驗(yàn)誤差僅在6%之內(nèi)。

3 基于代理模型的快速分析方法

為了繼續(xù)節(jié)約計(jì)算成本，可以采用機(jī)械設(shè)計(jì)與機(jī)器學(xué)習(xí)領(lǐng)域中的代理模型方法。

代理模型的基本思想為將復(fù)雜的仿真程序看作“黑箱”，使用合適的試驗(yàn)設(shè)計(jì)方法得到抽樣矩陣，抽樣矩陣的每一組設(shè)計(jì)變量的值為一個(gè)輸入向量。將抽樣矩陣輸入到仿真程序中獲取與之對(duì)應(yīng)的一次隨機(jī)失效壽命值，建立抽樣點(diǎn)-程序響應(yīng)值間的對(duì)應(yīng)關(guān)系作為代理模型的訓(xùn)練集，對(duì)選定的代理模型進(jìn)行訓(xùn)練以得到擬合或插值關(guān)系。建立代理模型的步驟如下。

第一步，確定設(shè)計(jì)變量及變量空間。

第二步，選取合適的試驗(yàn)設(shè)計(jì)技術(shù)（DOE），獲得設(shè)計(jì)變量的抽樣矩陣。

第三步，執(zhí)行仿真程序，獲得輸入樣本點(diǎn)對(duì)應(yīng)的輸出響應(yīng)。

第四步，選取合適的代理模型，建立代理關(guān)系，通過(guò)訓(xùn)練集對(duì)代理模型進(jìn)行訓(xùn)練。

第五步，精度校核。

為避免過(guò)多地輸入維度，需對(duì)輸入?yún)?shù)進(jìn)行篩選。首先確定代理模型中的設(shè)計(jì)變量，選取應(yīng)力水平S、裂紋隨機(jī)萌生參數(shù)m、裂紋隨機(jī)擴(kuò)展參數(shù)n三個(gè)參數(shù)為設(shè)計(jì)變量。其次需要確定3 個(gè)參數(shù)的變量空間，參數(shù)m、n服從正態(tài)分布，理論上的變量空間應(yīng)該是（-∞，+∞），一般情況下選取變量空間為[μ-5σ，μ＋5σ]，其中μ 為均值，σ 為標(biāo)準(zhǔn)差。應(yīng)力水平根據(jù)工程實(shí)際選擇變量空間。見(jiàn)表5。

表5 設(shè)計(jì)變量的變量空間

本文選取拉丁超立方法獲得設(shè)計(jì)變量的抽樣矩陣，建立了50 個(gè)樣本點(diǎn)，在仿真程序中使用蒙特卡羅法抽樣50 次獲得平均行為作為輸出響應(yīng)，建立多項(xiàng)式響應(yīng)面法代理模型，由輸入-輸出響應(yīng)作為訓(xùn)練集訓(xùn)練代理模型。

基于高斯過(guò)程回歸模型和多項(xiàng)式響應(yīng)面模型（RSM）分別訓(xùn)練訓(xùn)練集，通過(guò)交叉驗(yàn)證方法進(jìn)行機(jī)器學(xué)習(xí)驗(yàn)證，驗(yàn)證精度指標(biāo)見(jiàn)表6。

表6 模型學(xué)習(xí)驗(yàn)證指標(biāo)

由表6 可以看出高斯過(guò)程回歸模型精度更高，擬合效果也很好，高斯過(guò)程回歸模型響應(yīng)結(jié)果如圖5 所示。建立測(cè)試集，由測(cè)試數(shù)據(jù)對(duì)高斯過(guò)程回歸模型進(jìn)行測(cè)試驗(yàn)證，驗(yàn)證結(jié)果顯示誤差不超過(guò)10%，且節(jié)約了大量的時(shí)間計(jì)算成本。

圖5 高斯過(guò)程回歸模型響應(yīng)圖

4 結(jié)束語(yǔ)

本文建立了裂紋隨機(jī)萌生和擴(kuò)展模型，利用ABAQUS 二次開(kāi)發(fā)得到的應(yīng)力強(qiáng)度因子數(shù)值，根據(jù)蒙特卡羅法預(yù)測(cè)廣布疲勞損傷平均行為，廣布疲勞損傷平均行為與試驗(yàn)對(duì)比誤差在6%之內(nèi)；通過(guò)建立代理模型的方式，建立訓(xùn)練集和平均行為輸出響應(yīng)對(duì)高斯過(guò)程回歸模型和多項(xiàng)式響應(yīng)面模型進(jìn)行訓(xùn)練學(xué)習(xí)，由測(cè)試數(shù)據(jù)對(duì)模型進(jìn)行驗(yàn)證顯示，高斯過(guò)程回歸模型精度更高，訓(xùn)練后的預(yù)測(cè)效果誤差不超過(guò)10%，并且可以節(jié)約90%的計(jì)算時(shí)間成本，所構(gòu)建的高斯過(guò)程回歸模型可以精確地確定廣布疲勞損傷平均行為。