物理模型與高斯過程融合驅(qū)動的殘余應(yīng)力疲勞狀態(tài)評估

梁天佑， 尹愛軍， 陳 平， 方 杰

(重慶大學(xué) 機(jī)械工程學(xué)院 機(jī)械傳動國家重點實驗室, 重慶 400044)

金屬構(gòu)件受振動、沖擊等作用易導(dǎo)致疲勞失效，如支撐、連接螺栓等在長期振動作用下導(dǎo)致疲勞斷裂，繼而引發(fā)生產(chǎn)事故[1-2]。國內(nèi)外研究表明，殘余應(yīng)力可以表征金屬構(gòu)件的疲勞狀態(tài)[3]。振動引發(fā)疲勞可由殘余應(yīng)力表征，慕琴琴等[4]研究了發(fā)動機(jī)盤類零件振動-應(yīng)力-壽命關(guān)系。周磊等[5]對不同軌道結(jié)構(gòu)形式下地鐵車輛部件振動及應(yīng)力進(jìn)行分析，研究彈簧疲勞壽命。Kodama[6]通過疲勞試驗探究了載荷反復(fù)作用下殘余應(yīng)力與疲勞加載周期之間的關(guān)系提出了殘余應(yīng)力和疲勞循環(huán)周期的對數(shù)線性關(guān)系。Zhuang等[7-8]提出了基于初始?xì)堄鄳?yīng)力的殘余應(yīng)力-疲勞循環(huán)周期模型。Zhang等[9]通過車軸鋼疲勞試驗，認(rèn)為高周循環(huán)載荷作用下殘余應(yīng)力與裂紋萌生及擴(kuò)展有關(guān)。郭振坤等[10-11]則通過仿真模型研究了構(gòu)件的疲勞壽命演化。

殘余應(yīng)力松弛模型中常用的參數(shù)估計方法，主要有非線性最小平方(nonlinear least squares，NLS)法、貝葉斯法(Bayesian method，BM)、粒子濾波(particle filter，PF)法[12]等，對于包含隱變量的則有最大期望(expectation-maximization，EM)算法[13]等。NLS法是將非線性模型轉(zhuǎn)化為線性模型，再通過最小二乘法求解。BM是先假設(shè)模型未知參數(shù)符合某種先驗分布，在某一范圍選擇合適初始值，利用先驗分布采樣值的后驗概率值變化更新參數(shù)值，并迭代到收斂。PF法類似于BM，通過采樣大量帶權(quán)粒子確定物理模型參數(shù)。在殘余應(yīng)力對疲勞壽命影響研究中，常采用NLS方法進(jìn)行參數(shù)估計。

傳統(tǒng)物理模型如Kodama提出的對數(shù)線性模型等往往未充分考慮初始?xì)堄鄳?yīng)力與構(gòu)件疲勞狀態(tài)間的關(guān)系，以及材料、環(huán)境等因素的影響。本文在有限數(shù)據(jù)條件下，考慮初始?xì)堄鄳?yīng)力，試樣差異性以及環(huán)境因素等對疲勞狀態(tài)的影響，利用傳統(tǒng)物理模型設(shè)計高斯過程核函數(shù)，建立物理模型和數(shù)據(jù)驅(qū)動的高斯過程融合的殘余應(yīng)力疲勞狀態(tài)評估模型(Kodama-Gaussian process，K-GP)，并通過鋁合金進(jìn)行了對比試驗，驗證了模型有效性。

1 殘余應(yīng)力壽命模型

基于殘余應(yīng)力的壽命模型已得到相關(guān)研究。研究表明，殘余應(yīng)力越穩(wěn)定，材料疲勞壽命越高。Kodama通過研究退火碳鋼噴丸處理后單軸循環(huán)加載下表面殘余應(yīng)力演化規(guī)律，提出如式(1)的對數(shù)線性模型。

σN,rs=A+mlnN

(1)

式中：σN,rs為當(dāng)前循環(huán)周期對應(yīng)殘余應(yīng)力；A,m為待定超參數(shù)；N為當(dāng)前循環(huán)周期。然而該模型僅適用于殘余應(yīng)力的緩慢松弛階段。

殘余應(yīng)力松弛現(xiàn)象同時與應(yīng)力梯度、應(yīng)力比、材料應(yīng)力應(yīng)變、冷作硬化、粗糙度[14-16]有關(guān)。在考慮冷作硬化等因素的條件下，有如式(2)所示的壽命模型。

(2)

(3)

式(1)和式(3)模型均僅適用于構(gòu)件的緩慢松弛階段，式(3)同時考慮了初始?xì)堄鄳?yīng)力的影響。

2 基于殘余應(yīng)力的K-GP疲勞狀態(tài)評估

由于材料高應(yīng)力幅值會加速殘余應(yīng)力馳豫速率，使傳統(tǒng)物理模型計算結(jié)果在低周疲勞階段失準(zhǔn)。因此在評估時，一般不考慮低周數(shù)據(jù)。而基于擬合回歸的參數(shù)估計方法，要求物理模型本身簡潔，因此評估誤差會增大。由數(shù)據(jù)驅(qū)動的高斯過程，將殘余應(yīng)力松弛模型引入高斯過程核函數(shù)，建立構(gòu)件疲勞循環(huán)周期預(yù)測模型。

2.1 高斯過程

高斯過程回歸(Gaussian process regression，GPR)假設(shè)多變量特征符合聯(lián)合高斯分布，從而通過計算聯(lián)合分布的邊緣概率密度獲得目標(biāo)條件分布。GPR廣泛用于機(jī)械性能評估，壽命預(yù)測，金融股票走勢預(yù)測等問題有諸多研究和運用[17-19]。

(4)

式中：f(X*)為關(guān)于Y*的高斯分布，其分布參數(shù)與X*有關(guān)；K為核函數(shù)，用于定義樣本空間的內(nèi)積關(guān)系；K(X,X)為以K(xi,xj)為元素的矩陣。

則f(X*)概率分布的均值μ(X*)以及協(xié)方差Σ*如式(5)和式(6)。

μ*=K(X*,X)[K(X,X)+σ2I]-1×

[Y-μ(X)]+μ(X*)

(5)

Σ*=K(X*,X*)-K(X*,X)×

[K(X,X)+σ2I]-1K(X,X*)

(6)

式中，σ2為高斯噪聲的系數(shù)[20-21]。

根據(jù)上述分析，高斯過程回歸關(guān)鍵在于設(shè)計核函數(shù)K[22]。不同的核函數(shù)有不同的超參數(shù)，通常采取優(yōu)化極大似然邊緣概率估計核函數(shù)的超參數(shù)。設(shè)超參數(shù)向量為θ，則有

L(θ)=-ln[p(y|X,θ)]=

(7)

式中，Σθ為超參數(shù)θ下的協(xié)方差矩陣。

2.2 K-GP狀態(tài)評估模型

假設(shè)構(gòu)件當(dāng)前循環(huán)周期關(guān)于殘余應(yīng)力符合高斯過程。根據(jù)式(1)物理模型，可得到如式(8)的增量表達(dá)式。

σN,rs+ΔσN,rs=A+mS+mΔS

(8)

式中：ΔσN,rs為殘余應(yīng)力增量；S=lnN為對數(shù)疲勞周期；ΔS為對數(shù)疲勞周期增量。式(8)表明表示殘余應(yīng)力改變時(增量ΔS)引起對數(shù)疲勞周期變化(ΔσN,rs增量)成線性關(guān)系。

結(jié)合式(8)的線性過程和式(9)所示的高斯過程線性核函數(shù)

(9)

有

(10)

根據(jù)式(5)和式(10)，可建立式(11)所示的高斯預(yù)測模型。

(11)

對于式(1)所示的物理模型，通常只適用緩慢松弛階段。而由式(3)，初始?xì)堄鄳?yīng)力和疲勞周期呈指數(shù)關(guān)系?？紤]構(gòu)件材料差異、環(huán)境因素以及測量誤差等因素影響，本文在線性核函數(shù)的基礎(chǔ)上，融合式(12)所示的徑向基函數(shù)核

(12)

式中，σ2和l2為超參數(shù)，形成式(13)的殘余應(yīng)力疲勞狀態(tài)評估模型核函數(shù)。

K=KlKrbf+Kl

(13)

上述模型融合了物理模型和數(shù)據(jù)驅(qū)動方法，解決了數(shù)據(jù)驅(qū)動建模過程中的模型缺乏可解釋性問題，并綜合了初始?xì)堄鄳?yīng)力對評估的影響。

3 試驗驗證

3.1 試驗試樣

根據(jù)GB/T 228.1—2010《金屬拉伸試驗標(biāo)準(zhǔn)試驗類型及尺寸》，設(shè)計了如圖1所示的2024鋁合金試樣，其試件中間噴完部位表面粗糙度0.4 μm，材料性能參數(shù)如表1所示。

圖1 試驗樣件Fig.1 Experimental samples

表1 2024鋁合金力學(xué)性能Tab.1 Mechanical properties of 2024 aluminum alloy

3.2 試驗過程

3.2.1 表面噴丸強(qiáng)化

表面強(qiáng)化技術(shù)可在基本不改變材料基體的前提下可顯著提高其材料表層的抗疲勞性能[23]，噴丸區(qū)域如圖2框選區(qū)域，噴丸丸粒直徑0.3 mm，出口壓力0.3 MPa, 噴丸時間30 s, 覆蓋率100%。

圖2 噴丸件實物Fig.2 Shot peening parts

3.2.2 初始?xì)堄鄳?yīng)力檢測

采用μ-X360nX射線衍射儀檢測各試樣的初始?xì)堄鄳?yīng)力，并選取其中7個殘余應(yīng)力相近的試樣進(jìn)行后續(xù)試驗。初始?xì)堄鄳?yīng)力檢測結(jié)果如表2所示。1號試樣用來載荷測試。2號試樣是3.3節(jié)中全壽命試驗樣件，3號、4號、5號、6號、7號試樣用于階段疲勞試驗。

表2 各樣件初始?xì)堄鄳?yīng)力表Tab.2 Initial residual stress of each piece

3.2.3 拉伸疲勞試驗

應(yīng)力比為0.1，正弦加載，頻率為20 Hz，加載載荷為15 kN。2號試樣全壽命加載試驗結(jié)果如表3所示。4～7號試樣階段加載，結(jié)果如表4所示。

表3 2號試樣疲勞殘余應(yīng)力檢測結(jié)果Tab.3 Fatigue residual stress test results of No.2 sample

表4 階段性試驗疲勞殘余應(yīng)力檢測結(jié)果Tab.4 Test results of fatigue residual stress in step test

3.3 疲勞狀態(tài)評估

通過全壽命試驗以及階段性試驗數(shù)據(jù)建立模型，結(jié)合蒙特卡洛方法對初始?xì)堄鄳?yīng)力以及當(dāng)前殘余應(yīng)力進(jìn)行隨機(jī)采樣增加樣本總量，評估對象是對數(shù)疲勞周期。圖3(a)是融合模型評估結(jié)果;圖3(b)中不同初始?xì)堄鄳?yīng)力下殘余應(yīng)力-疲勞循環(huán)周期曲線，每條曲線對應(yīng)一個試樣。該曲線簇反映不同試樣具有不同初始?xì)堄鄳?yīng)力下的演化特性。圖4、圖5分別是式(1)、式(3)利NLS方法得到評估結(jié)果(以Kodama模型和Omar模型表示)。

圖3 K-GP評估結(jié)果Fig.3 Assessment results of K-GP

圖4 Kodama模型評估結(jié)果Fig.4 Assessment results of Kodama model

由圖4圖5可知，可驗證傳統(tǒng)物理模型并不適用于構(gòu)件疲勞的低周階段，圖5橫坐標(biāo)為殘余應(yīng)力比即當(dāng)前殘余應(yīng)力與初始?xì)堄鄳?yīng)力比值，說明構(gòu)件疲勞過程與初始?xì)堄鄳?yīng)力有關(guān)，圖3(c)中表明不同初始?xì)堄鄳?yīng)力則對應(yīng)于不同的構(gòu)件疲勞曲線，但是其變化趨勢是相似的，且KGP適用于低周疲勞階段。

圖5 Omar模型評估結(jié)果Fig.5 Assessment results of Omar model

誤差結(jié)果如圖6所示，以2號樣件總壽命百分比表示，誤差計算公式為

圖6 評估誤差對比Fig.6 Evaluation error comparison

(14)

Nf_t取二號試樣的總疲勞周期。由于我們評估模型輸出的是對數(shù)周期，因此需將對數(shù)周期轉(zhuǎn)換為周期，對數(shù)周期在KGP中認(rèn)為是符合正態(tài)分布，則周期符合對數(shù)正態(tài)分布，其期望計算公式為

E(Np)=eE(S)+Var(S)/2

(15)

式中，S為模型評估的對數(shù)周期。

圖6給出各個點對應(yīng)誤差并通過3階多項式擬合給出誤差趨勢，相較于其他疲勞壽命模型，本文提出的K-GP模型引入初始?xì)堄鄳?yīng)力，有效的提高了疲勞狀態(tài)評估的準(zhǔn)確性，Omar模型效果最差。

4 結(jié) 論

本文分析了殘余應(yīng)力演化傳統(tǒng)物理模型的局限性，結(jié)合初始?xì)堄鄳?yīng)力，將數(shù)據(jù)驅(qū)動的高斯過程融入到物理模型中，提高了殘余應(yīng)力疲勞狀態(tài)評估模型的準(zhǔn)確性并使模型適用于低周疲勞階段。K-GP模型目前涉及初始?xì)堄鄳?yīng)力對疲勞狀態(tài)評估的影響，后續(xù)可考慮表面粗糙度，冷作硬化等對模型的影響；并優(yōu)化模型，減小對樣本數(shù)據(jù)量的依賴。