鋼結構中疲勞裂紋增長概率模型和Bayesian動態預測

陳夢成，方 葦，楊 超，謝 力

(1.華東交通大學 省部共建軌道交通基礎設施性能監測與保障國家重點實驗室，江西 南昌 330013；2.華東交通大學 土木建筑學院，江西 南昌 330013)

近50年來，隨著我國改革開放和經濟實力提升，基礎設施建設得到蓬勃發展。大量既有橋梁結構在使用環境和使用荷載作用下結構性能逐漸下降，產生疲勞裂紋并不斷擴大，嚴重影響了橋梁的耐久性和適用性，甚至有部分橋梁尚未達到設計使用年限就已經出現重大安全隱患。準確地預測橋梁工程中鋼結構性能退化和疲勞壽命評估是進行橋梁結構維修與決策、延長橋梁使用壽命的關鍵[1-2]。工程結構可靠度理論的研究是個長久的課題，早在20世紀70年代發達國家就已開始進行服役工程構件及結構系統的可靠性分析及預測的研究[3-4]，從20世紀90年代開始，結構時變可靠性評估已成為國內外學者的研究熱點[5-6]，但結構時變可靠性的計算方法還遠不成熟。

結構可靠性的研究通常是通過考慮影響結構性能退化的各不確定因素，推算結構抗力的時變概率特性[7-10]，從而建立結構抗力隨機過程模型[7-8,11-12]。但在具體的實踐中，因橋梁結構耐久性問題的時間周期長、影響因素復雜、不確定性大、不同結構從材料到構件再到結構設計的差異性等因素，很難用某一種模型來概括。已有的工程結構失效預測和可靠性評估模型也無法做到充分考慮實際工程中各種復雜因素對各參數的影響。因此研究結構可靠性就要求建立合理的結構性能隨機過程模型[13]。現代無損檢測技術常用于橋梁結構系統的健康監測中，受此啟發，眾多研究者提出利用健康監測數據進行橋梁結構的疲勞壽命預測[14-15]，而貝葉斯更新理論可以實現耦合服役工程結構性能退化過程(模型的不確定性)和監測數據，動態更新結構性能退化隨機過程模型。Bayesian更新理論的數值計算方法中，馬爾卡夫鏈蒙特卡羅(Markov Chain Monte Carlo，簡稱MCMC)抽樣方法最為流行，因為它利用計算機的強大計算能力和仿真能力[16-17]。

疲勞是導致橋梁結構失效的原因之一，而疲勞裂紋在既有橋梁結構中非常普遍，其增長過程是一個隨機過程，疲勞裂紋的增長會引起結構性能退化。本文擬針對橋梁結構中因疲勞裂紋損傷而引起的結構性能退化以及時變可靠性預測的科學問題，采用貝葉斯更新理論，開展相關研究。文中介紹了Bayesian更新理論在疲勞裂紋增長模型的參數更新和構件壽命預測的應用中如何實施，并進行了具體的推導論證。最后結合具體算例，采用Matlab編程，對疲勞裂紋損傷構件在未來一段時間內的裂紋增長和時變可靠性進行預測，以檢驗本文方法的可行性和有效性。

1 Bayesian更新定理

(1)

依據上面說明，Bayesian更新公式可用文字表述為：更新概率(后驗概率)=新觀測信息的條件概率(似然度×先驗概率)/新觀測信息的全概率(標準化常量)，也就是說，后驗概率與先驗概率和似然度的乘積成正比。另外，比例P(B/Ai)/P(B)有時也被稱作標準似然度，因此，Bayesian更新公式可進一步用文字表述為：后驗概率= 標準似然度×先驗概率。

2 構件疲勞裂紋損傷性能隨機退化建模

2.1 疲勞裂紋擴展確定性模型

對于一個含有表面裂紋的鋼構件，依據線性彈性斷裂力學(LEFM)理論，裂紋端部應力強度因子K在某一應力水平σ作用下為

(2)

式中：Y為應力強度修正系數，反映了構件和裂紋的幾何形狀；a為半裂紋長度。

在構件承受靜荷載(一般使用環境)時，只有其應力水平達到臨界應力σc時，即裂紋端部的應力強度因子達到臨界值Kc時，才會立即發生失穩斷裂，此時相對應的臨界裂紋長度為ac。當靜應力水平降低到σ0(σ0<σc)，則構件不會發生失穩破壞。但如果構件承受一個與靜應力σ0大小相等的往復循環應力，則初始裂紋a0在循環應力σ0作用下發生緩慢擴展。根據著名的Paris半經驗公式[18]，疲勞裂紋擴展速率為

(3)

假設疲勞裂紋擴展過程中，Δσ和Y保持不變，則可取ΔK(a)=a，積分式(3)可得

(4)

式中：a(0)為N=0時構件的初始疲勞裂紋長度；a(N)為荷載循環次數為N時的裂紋長度。

2.2 疲勞裂紋損傷性能退化概率模型

令式(4)中材料參數C=θ1、m-1=θ2，并假設這些參數均為不確定性參數，則對式(4)兩邊取自然對數，并令y=ln(a(N)/a(0))，則有

(5)

式中：yij為第i(i=1,2,…,n)個疲勞裂紋在荷載循環到tj(j=1,2,…,m)=N時對應的相對疲勞裂紋擴展長度，工程結構構件相應性能退化；tj為時間變量，表示作用在構件上的疲勞荷載循環次數，為性能變量的觀測時間；η(·)為描述性能變量變化的軌跡函數，它是t的單調非線性函數，在許多情況下它也是參數θ(θ1i,θ2i)的非線性函數；θ(θ1i,θ2i)∈Θ?Rq為性能退化函數η(·)中的未知隨機參數向量；q為參數向量位數；Θ為參數向量的取值空間；εij為反映性能退化過程中第i個疲勞裂紋在tj時刻相對擴展長度的不確定性的誤差項。

參數向量的“隨機性”應該理解為，對于同一軌跡的性能退化，模型參數的取值是確定的；對于不同軌跡的性能退化，模型參數的取值是隨機的，這種隨機性可以由某個概率分布刻畫。為此，進一步假設：①對應特定的tj，隨機誤差項序列εij，i=1,2,…,n獨立同分布，也就是說，εij在空間上是獨立的；②對于特定的疲勞裂紋i，隨機誤差項序列εij，j=1,2,…,m是相關的；③隨機誤差項序列εij服從均值為0，未知協方差矩陣為Σε的多維隨機變量正態分布，即N(0,Σε)。

定義yi=(yi1yi2…yim)T和ηi=(ηi1ηi2…ηim)T，其中T為矩陣轉置；用fMVN(x/μ,Σ)表示均值為μ，協方差為Σ的多維隨機變量正態分布概率密度函數，因此，在給定ηi的情況下，由式(4)可知

(6)

式中：Σε為隨機誤差項m×m階的協方差方陣，其中矩陣元素等于ρklσkσl，(k=1，2，…，m；l=1，2，…，m)；ρkl為荷載循環次數k次和l次時隨機誤差項之間的相關系數；σk和σl分別為荷載循環次數k次和l次時隨機誤差項的標準方差。

3 疲勞裂紋損傷性能退化模型隨機參數的Bayesian更新

3.1 先驗分布

3.2 性能退化過程似然函數和后驗分布

對于特定的疲勞裂紋i，假定在tj(j=1,2,…,m)時刻有損傷性能退化量測值yi1,yi2,…,yim。由于隨機誤差項序列εij，j=1,2,…,m是獨立同分布的隨機變量，在給定θ和σε的條件下，性能退化樣本yi1:m的條件聯合分布函數或似然函數為

(7)

在已知聯合先驗分布π(θ,σε)的基礎上，得到性能退化量測值D(=yi1,yi2,…,yim)后就可以利用Bayesian更新式(1)來推斷后驗分布p(θ,σε|D)，即有

(8)

式中：Θ、Ω分別為θ、σε的取值空間。

θ1(1)從全條件分布π(θ1(0))π(θ2(0))π(θ3(0))L(y1:m|θ1(0),θ2(0),θ3(0))產生；

θ2(1)從全條件分布π(θ1(1))π(θ2(0))π(θ3(0))L(y1:m|θ1(1),θ2(0),θ3(0))產生；

θ3(1)從全條件分布π(θ1(1))π(θ2(1))π(θ3(0))L(y1:m|θ1(1),θ2(1),θ3(0))產生；

?

以上Gibbs抽樣技術的關鍵在于僅考慮了單變量的條件分布。抽樣樣本可以通過選擇一個建議分布來生成。這里選擇隨機游走鏈作為建議分布為

(9)

(10)

重復上述迭代過程直到各分量Θk+1的遍歷均值穩定后，可生成一系列所需要采集的樣本{Θ1,Θ2,Θ3,…}。利用這些后驗(更新)樣本，即可估計后驗分布和隨機參數的置信區間。

4 疲勞裂紋損傷性能退化與時變可靠性評估

4.1 單根疲勞裂紋極限狀態方程

建立構件單根疲勞裂紋長度安全余量方程，作為進行可靠性分析的極限狀態方程為

Z=acr-y(tk+Δt)≥0

(11)

式中：Z為疲勞裂紋長度的安全余量；acr為給定的疲勞裂紋臨界值；y(tk+Δt)為疲勞裂紋在tk+Δt時刻的預測值，tk為當前時刻；y(tk+Δt)=η(tk,θ1,θ2)+εk+1，其中εk+1為疲勞裂紋損傷性能退化軌跡模型在tk+Δt的隨機誤差項，假設它與ε1,ε2,…,εk獨立同分布，即：εk+1～N(0,σε)；Δt為未來時間區間長度。

4.2 疲勞裂紋損傷構件時變可靠性Bayesian預測

要了解tk+Δt時刻的疲勞裂紋損傷性能退化后的構件可靠性，就必須依據當前時刻以及過去的疲勞裂紋長度量測數據，對其可靠度進行評價。根據式(11)，構件在未來某一時刻(tk+Δt)的可靠概率為

(12)

式中：g(θ1,θ2,σε,tk,Δt,acr)=Z=acr-y(tk+Δt)，由4.1節可知服從N(0,σε)分布。

一般說來，要解析求解式(12)積分通常很困難，因此，我們用Bayesian更新式(8)、式(9)、式(10)所建立的抽樣方法，以解決預測和可靠度分析問題。在獲得未知隨機參數和隨機誤差項的后驗分布后，可將式(12)的求解轉化為

(13)

每次達到疲勞試驗的設計試驗節點，就會收集到一個新的量測數據。利用收集到的歷史量測數據，相應地更新似然函數，即可利用式(13)計算未來一段時間Δt(可以取不同值)內疲勞裂紋損傷構件的可靠性分析和預測，并繪制一條時變可靠度變化曲線。

5 算例

為了說明以上建立的Bayesian更新方法和MCMC抽樣技術的應用，以文獻[21]中的疲勞裂紋擴展量測數據為例，對疲勞裂紋損傷構件的時變可靠性進行預測。疲勞試驗試件共21個，試驗時長(疲勞荷載循環次數)為0.12百萬次，每隔0.01百萬次對每組疲勞裂紋擴展長度進行量測。所有試件在疲勞試驗開始前的初始裂紋長度均為22.86 mm，設定試件疲勞失效時裂紋長度的臨界值acr為40.64 mm。21個試件的疲勞裂紋擴展量測數據見表1[21]，與試件相應的疲勞失效時間見表2[22]。由表1可知，疲勞試件13～21在疲勞試驗終止時，裂紋長度尚未達到臨界值。利用表1量測數據、最小二乘法和統計學知識，對式(5)中的模型參數進行近似估計，可得到模型參數的估計值θ1和θ2，參見表2，表2中還列出了文獻[22]給出的疲勞失效時間預測。

表1 疲勞裂紋擴展長度量測數據[21] mm

表2 疲勞裂紋損傷性能退化軌跡模型參數與疲勞失效時間預測

5.1 算例1

選擇第9條試件的疲勞裂紋損傷性能退化軌跡作為分析疲勞損傷構件時變可靠性分析及預測的研究對象，其余20條試件的退化軌跡用來提取退化軌跡模型隨機參數θ1和θ2的先驗信息。

依據第2、3、4節的理論，采用Matlab語言編制了Bayesian更新MCMC法源程序，該程序可讀性強，而且力求忠實MCMC理論框架。在數值模擬過程中，模型參數θ1、θ2和σε的Markov鏈的計算初始值分別取為8、5和0.5。θ1和θ2的建議分布采用正態分布，標準方差均為σ=0.1；σε的建議分布，按照其先驗分布形式，假定它服從對數正態分布，取其標準方差σ′=1。Markov鏈總長為10 000次。獲得的三個模型隨機參數的后驗樣本的時間序列，見圖1。由圖1可知，在抽樣過程中三個參數的后驗樣本很快進入穩定期。去掉初始抽樣(燃燒期)1 000次后的后驗分布的直方圖見圖2。由圖2 可知，模型參數θ1、θ2的后驗分布均表現出正態分布形式，σε的后驗分布表現出對數正態分布形式。需要指出的是：對于馬爾卡夫鏈燃燒期的選擇，在學術上目前還沒有一個統一的標準，本文算例燃燒期的選取根據時間序列圖中所示的馬爾卡夫鏈的收斂情況為初始抽樣的10%；在M-H算法中，由于采用了隨機游走抽樣，所以對參數的建議分布的選取要求不是很嚴格，并不要求建議分布與其后驗分布接近，而且計算結果收斂快。

圖1 參數的后驗樣本時間序列

圖2 模型參數θ1、θ2、σe的后驗分布直方圖

t∈(0.10,0.11]時間內疲勞裂紋損傷性能退化構件時變可靠度預測曲線見圖3(a)。由圖3(a)可知，t∈(0.10,0.11]時間內構件時變可靠度預測曲線隨時間增長，下降非常緩慢，可靠性變化范圍在1.0～0.999 9之間，這說明構件在t∈(0.10,0.11]內失效的可能性很小。圖3(b)給出了t∈(0.11,0.12]時間內構件時變可靠度預測曲線的變化趨勢見圖3(b)。由圖3(b)可知，該時變可靠度預測曲線與t∈(0.10,0.11]時間內的時變可靠度預測曲線的變化特征明顯不同，下降速度非常快，當t=0.115左右時變可靠度降為0。如果取時變可靠度下降10%作為構件疲勞失效的閾值，那么從圖中曲線可以看出，此時t=0.113 4，對應的可靠度為0.900 8，這與表2給出的第9條試件退化軌跡的疲勞失效時間t=0.113 0完全吻合，證實了本文方法的有效性。

圖3 t∈(0.10,0.11]和t∈(0.11,0.12]時間內時變可靠度預測曲線

5.2 算例2

選擇在疲勞試驗終止時，裂紋長度尚未達到臨界值的第18條試件的疲勞裂紋損傷性能退化軌跡作為預測疲勞損傷構件的疲勞裂紋增長的研究對象，構建疲勞裂紋增長概率模型隨機參數的更新方法，動態預測疲勞裂紋損傷性能退化軌跡。疲勞裂紋損傷性能退化軌跡與時變可靠度預測曲線見圖4。由圖4(a)可知，當t=0.11時，疲勞裂紋長度為32.58 mm，實際量測得的疲勞裂紋長度為32.51 mm，絕對誤差為0.07；當t=0.12時，疲勞裂紋長度為34.08 mm，實際量測得的疲勞裂紋長度為34.29 mm，絕對誤差為0.21。證實了本文方法的準確性。由圖4(b)可知，第18條疲勞裂紋在t=0.151 9時，疲勞裂紋長度自40.63 mm擴展至40.66 mm，達到臨界值，此時試件疲勞失效。由圖4(c)可知t∈(0.12，0.16]時間內的時變可靠度預測曲線與疲勞裂紋損傷性能退化軌跡預測結果吻合。

圖4 疲勞裂紋損傷性能退化軌跡與時變可靠度預測曲線

圖5 t∈(0.05,0.12]時間內疲勞裂紋損傷性能退化軌跡預測曲線

圖5為t∈(0.05，0.12]時間內疲勞裂紋損傷性能退化軌跡預測曲線，其中第1次更新是以t=0.04時刻的已收集量測數據，相應的更新似然函數后得到的疲勞裂紋性能損傷退化軌跡曲線，第2次更新是以t=0.05時刻收集到的歷史量測數據進行更新得到的退化軌跡曲線。

由圖5可知，利用數據更新疲勞裂紋增長概率模型后得到的退化軌跡與實際量測得到的退化軌跡曲線相近，第1次更新后得到的退化軌跡與實際量測得到的退化軌跡的平均絕對誤差為0.40 mm，均方根誤差為0.44，平均絕對百分比誤差為1.32%；第2次更新的退化軌跡更接近實際量測的退化軌跡，其平均絕對誤差僅為0.30 mm，均方根誤差為0.33，平均絕對百分比誤差為0.97%。對比上述三個評價指標可知，第2次更新后得到的預測模型的平均絕對誤差、均方根誤差以及平均絕對百分比誤差均小于第1次更新后得到的預測模型，可靠性更高。同時兩次更新后得到的退化軌跡預測圖相近，更進一步證實了貝葉斯更新理論在結構疲勞失效預測中的有效性和精確性，由此建立的鋼構件性能退化損傷的貝葉斯動態預測模型更符合實際工程。

6 結論

本文給出了一個基于線彈性斷裂力學和Paris半經驗疲勞裂紋增長速率公式的概率模型，該模型可用于疲勞荷載作用下的橋梁結構中。在此基礎上，進一步采用Bayesian更新理論和MCMC模擬技術，耦合實際結構新量測數據，建立了疲勞裂紋增長模型中隨機參數的更新評價方法，進而動態預測疲勞裂紋損傷性能退化軌跡和結構時變可靠性。疲勞裂紋增長概率模型考慮了實際結構量測數據的系統誤差和隨機誤差，它們在Bayesian更新過程中被處理為超級參數。通過實際應用算例研究，可以得到以下結論。

(1)本文提出的基于MCMC算法的Bayesian更新方法在進行模型參數修正的過程中充分考慮了模型參數的先驗信息，并不斷地利用新量測到的疲勞裂紋增長數據，有效地降低了抽樣過程中不確定因素對抽樣結果的影響，而且與僅考慮一次樣本數據的傳統統計學的計算結果相比，Bayesian更新統計學的計算結果更貼近于工程實際。

(2)用本文方法計算得到的可靠性預測曲線、疲勞裂紋損傷性能退化軌跡預測曲線與疲勞試驗結果高度吻合，說明本文建議的方法能夠有效預測反映結構從疲勞裂紋擴展開始至疲勞失效的全過程，適用于長期的、依時性動態過程預測且結果可靠，從而可為結構維修與養護方案、決策等提供科學依據。

(3)本文方法的準確性建立在合理的先驗信息的基礎上，因此對量測數據的精度和數量要求較高，此外，當考慮影響因素較多、樣本數增加時，MCMC算法的計算效率較低，需要發展更高效率的抽樣算法。

(4)文中的量測數據為實驗條件所得，與實際工程中的橋梁結構疲勞裂紋增長有一定的區別，在未來的研究工作中，將更注重工程實際。通過分析實際監測數據的統計特征，建立考慮塑性影響的疲勞裂紋統計模型，結合數值方法建立有限元模型，進行鋼結構疲勞壽命預測。