王藝斐 蘇 春 謝明江
(東南大學(xué)機(jī)械工程學(xué)院,南京 211189)
管道是一種高效的石油和天然氣運輸方式.然而油氣管道長期在高溫、高壓下工作,并且油氣中的H2S、CO2、Cl-等腐蝕性介質(zhì)導(dǎo)致管道內(nèi)壁腐蝕.此外,管道外壁長期經(jīng)受不同性質(zhì)土壤、雜散電流的侵蝕.在多重因素的共同作用下,管壁會出現(xiàn)腐蝕、開裂、穿孔等現(xiàn)象[1].在突發(fā)的外力作用下或當(dāng)油氣的運行壓力大于管壁承載能力時,管線會發(fā)生破裂并導(dǎo)致油氣泄漏,造成環(huán)境污染和生態(tài)破壞[2].因此,開展管道剩余壽命預(yù)測、有效減少管道失效事件發(fā)生具有重要意義.
近年來,基于可靠性理論的油氣管道剩余壽命預(yù)測研究受到關(guān)注.Mosallam等[3]利用概率模型表征腐蝕尺寸和管道特性的不確定性.Gong等[4]考慮管段泄漏、爆裂失效模式和不同腐蝕缺陷下極限狀態(tài)函數(shù)的相關(guān)性,采用一階可靠度方法評價含腐蝕缺陷管道系統(tǒng)可靠性.Bouledroua等[5]采用二階矩方法評估受腐蝕管道因內(nèi)壓過大而失效的可靠性問題.但上述方法存在不足,主要表現(xiàn)在管道壽命預(yù)測精度不高、不適合于某些極端條件等[6].
管道腐蝕常表現(xiàn)為壁厚減薄、強(qiáng)度下降,是典型的性能退化過程.一些學(xué)者采用隨機(jī)過程描述管壁性能退化的不確定性.張新生等[7]基于逆高斯(IG)過程建立油氣管道腐蝕退化的狀態(tài)空間模型,開展管道壽命預(yù)測和維修決策優(yōu)化.Bazán等[8]提出非線性隨機(jī)過程的腐蝕增長模型,通過對比線性隨機(jī)變量腐蝕模型,驗證模型的擬合效果.Wang等[9]建立埋地管道結(jié)構(gòu)隨機(jī)腐蝕增長模型,采用幾何布朗橋過程模擬動態(tài)外部點蝕過程,用于地下管線結(jié)構(gòu)的完整性管理.Ossai等[10]采用非齊次線性增長純生馬爾可夫模型預(yù)測油氣管道內(nèi)腐蝕的深度分布.Amaya-Gómez等[11]基于在線檢測數(shù)據(jù),集成應(yīng)用Lévy過程、伽馬過程和復(fù)合泊松過程構(gòu)建管道腐蝕退化模型.Zhang等[12]根據(jù)地下能源管道腐蝕缺陷的檢測數(shù)據(jù),基于IG過程建立腐蝕深度增長模型,采用馬爾可夫鏈蒙特卡洛仿真得到預(yù)測結(jié)果.
長輸石油管道退化通常包含腐蝕穿孔、局部爆裂等多種失效模式,并且有不同的性能指標(biāo)[13].為更準(zhǔn)確地描述管道健康狀態(tài),有必要采用2個或多個性能指標(biāo)表征管道健康狀態(tài),并考慮性能指標(biāo)的相關(guān)性.Copula函數(shù)適合于描述多性能指標(biāo)間的相關(guān)性.Zhou等[14]利用Gaussian Copula、隨機(jī)過程理論描述管道缺陷生長之間的依賴關(guān)系.金曉航等[15]以Copula函數(shù)定義軸承2個性能指標(biāo)的相關(guān)性,得到軸承剩余壽命的聯(lián)合概率密度函數(shù).
本文以長輸石油管道為研究對象,考慮管道剩余壁厚、剩余強(qiáng)度2個退化性能指標(biāo),提出一種基于二元IG過程的腐蝕輸油管道剩余壽命預(yù)測方法.首先采用Copula函數(shù)描述性能指標(biāo)之間的相關(guān)性,建立管道剩余壽命的聯(lián)合密度函數(shù);采用期望值最大化(EM)算法估計失效概率密度函數(shù)中的未知參數(shù);采用貝葉斯方法更新參數(shù)分布,完成石油管道剩余壽命的實時預(yù)測.最后通過數(shù)值案例驗證所提方法的可行性與有效性.
逆高斯過程是一種具有單調(diào)性的隨機(jī)過程[16].假設(shè)產(chǎn)品退化量Y(t)服從逆高斯過程,它具有如下性質(zhì).
1) 當(dāng)t=0時,Y(0)=0.
2) 對任意t4>t3>t2>t1,Y(t4)-Y(t3)與Y(t2)-Y(t1)相互獨立.
其中,均值函數(shù)Λ(t)為一個單調(diào)增函數(shù),ΔΛ=Λ(t+Δt)-Λ(t),Λ(0)=0.逆高斯分布IG(βΛ(t),ηΛ2(t))的概率密度函數(shù)為
(1)
式中,β為退化率;η為尺寸參數(shù).
采用多退化量預(yù)測設(shè)備剩余壽命時,需要考慮退化量之間的相關(guān)性.Copula函數(shù)是一種連接一元邊緣分布函數(shù)及其對應(yīng)多元聯(lián)合分布函數(shù)的特殊函數(shù),通過構(gòu)建每個邊緣退化過程的聯(lián)合分布函數(shù)來定義多性能退化量之間的相關(guān)性.
1.2.1Copula函數(shù)定義
當(dāng)函數(shù)C(u,v)滿足如下性質(zhì)時,稱為二元Copula函數(shù)[17]:
1)C(u,v)的定義域為[0, 1]×[0, 1].
2)C(u,v)具有零基面,且二維遞增.
3) 對于任意的u,v∈[0, 1],滿足C(u, 1)=u,C(1,v)=v.
當(dāng)已知各個變量的邊緣分布函數(shù)時,可根據(jù)Sklar定理得到它們的聯(lián)合分布函數(shù).
1.2.2Sklar定理
令F(x1,x2, …,xN)為具有邊緣分布F1(x1),F2(x2),…,FN(xN)的N元聯(lián)合分布函數(shù),則存在一個Copula函數(shù)C(u1,u2, …,uN)滿足[18]
F(x1,x2,…,xN)=
C[F1(x1),F2(x2), …,FN(xN)]
(2)
若F1(x1),F2(x2),…,FN(xN)為連續(xù)函數(shù),則C(u1,u2, …,uN)唯一確定;反之,若F1(x1),F2(x2), …,FN(xN)為一元分布函數(shù),C(u1,u2, …,uN)是一個Copula函數(shù),則由式(2)確定的F(x1,x2, …,xN)是具有邊緣分布F1(x1),F2(x2),…,FN(xN)的N元聯(lián)合分布函數(shù).
由Sklar定理可知,采用Copula函數(shù)描述2個退化性能指標(biāo)之間的相關(guān)特性可表示為[19]
F(x1,x2)=C[F1(x1),F2(x2)]
(3)
式中,F(xiàn)1(x1)、F2(x2)分別為2個性能指標(biāo)對應(yīng)的分布函數(shù).
采用首達(dá)時間(first hitting time)表征退化量與壽命L之間的關(guān)系,并據(jù)此開展剩余壽命預(yù)測.不失一般性,假設(shè)管道退化量隨時間逐漸增加,其失效閾值為ω.當(dāng)退化量超過ω時即認(rèn)為發(fā)生失效,將管道剩余壽命預(yù)測問題轉(zhuǎn)化為預(yù)測退化量何時會到達(dá)失效閾值ω.
給定失效閾值ω,管道的壽命L可定義為[20]
L=inf{lk:Y(tk+lk)≥ω|Y(0)<ω}
(4)
式中,lk為tk時刻的剩余壽命.
設(shè)剩余壽命的PDF函數(shù)為fLk(lk|Yk),則有
(5)
式中,Yk= {y0,y1, … ,yk}為tk時刻的退化數(shù)據(jù)集,yk為tk時刻的退化數(shù)據(jù);μβ為均值;σβ為標(biāo)準(zhǔn)差;φ(x)為標(biāo)準(zhǔn)正態(tài)分布的概率密度函數(shù);Φ(x)為標(biāo)準(zhǔn)正態(tài)分布累計分布函數(shù).
當(dāng)管道剩余壁厚低于某臨界值時,油氣管道可能會發(fā)生泄漏失效.根據(jù)管道完整性管理規(guī)定,通常將腐蝕深度超過管壁厚度的80%作為管道需要維修的判據(jù)[21].通常,管道公司會不定期開展管道檢測以獲得腐蝕深度等基礎(chǔ)數(shù)據(jù).本文以管道剩余壁厚作為退化性能指標(biāo).考慮到不同性能指標(biāo)的量綱存在差異,為便于分析進(jìn)行歸一化處理,剩余壁厚退化量定義如下:
(6)
式中,δ為剩余壁厚退化量;d0為初始壁厚;d(t)為t時刻時的管道壁厚.
文獻(xiàn)[11-13]多采用爆破壓力表征管道的承載能力,當(dāng)所計算的爆破壓力低于管道的運行壓力時管道將會失效.相關(guān)的國際標(biāo)準(zhǔn)包括B31G、DNV F101和PCORRC.根據(jù)上述標(biāo)準(zhǔn)可以計算管道爆破壓力[22].近年來,有限元方法在管道爆破壓力計算中受到重視.
但是,在某種條件下爆破壓力并不適用于作為管道剩余強(qiáng)度退化指標(biāo).在運行過程中,隨著腐蝕缺陷的擴(kuò)展,管道承壓能力不斷降低、剩余強(qiáng)度下降.本文綜合考慮等效應(yīng)力和屈服強(qiáng)度來描述剩余強(qiáng)度退化,采用有限元方法完成油氣管道建模,通過仿真得到不同腐蝕缺陷下管壁的等效應(yīng)力.當(dāng)最大等效應(yīng)力超過屈服強(qiáng)度時,管道將發(fā)生斷裂失效.管道剩余強(qiáng)度退化量的定義如下:
(7)
式中,ε為剩余強(qiáng)度退化量;G為屈服強(qiáng)度;τmax為最大等效應(yīng)力值.
建立剩余壁厚和剩余強(qiáng)度2個性能退化量的剩余壽命邊緣密度函數(shù),利用Copula函數(shù)分析2個性能退化量的相關(guān)性,可獲得管道剩余壽命的聯(lián)合概率密度函數(shù)[23],即
(8)
式中,f(x1,x2)為聯(lián)合概率密度函數(shù);f1(x1)、f2(x2)分別為2個性能指標(biāo)對應(yīng)的概率密度函數(shù);c(x)為Copula函數(shù)的密度函數(shù).
通過計算聯(lián)合概率密度函數(shù),可得到管道失效概率及其剩余壽命.Copula函數(shù)類型眾多,常用的Copula函數(shù)有Gaussian、T、Gumbel、Frank四種類型.不同Copula函數(shù)會導(dǎo)致不同的分析結(jié)果,因此需結(jié)合工程實際選擇合適的Copula函數(shù).
AIC(Akaike information criterion)信息準(zhǔn)則可用于評價模型擬合數(shù)據(jù)的效果,具有廣泛的適用性.根據(jù)AIC信息準(zhǔn)則計算得到的R值越小說明擬合效果越好.本文采用AIC信息準(zhǔn)則選擇合適的Copula函數(shù)[24]:
R=-2lnA+2m
(9)
式中,A為模型對應(yīng)的似然函數(shù);m為模型中參數(shù)的個數(shù).
假設(shè)Δyi=yi-yi-1對應(yīng)從時刻ti-1到ti的退化增量,Δti=ti-ti-1.在給定參數(shù)β的條件下,退化數(shù)據(jù)Yk的抽樣分布函數(shù)可以表示為[25]
(10)
(11)
其中
(12)
(13)
(14)
本文采用EM算法估計模型中的未知參數(shù),基本步驟如下:① 對隱變量計算對數(shù)似然函數(shù)的期望,記為E-step;② 對求過期望后的對數(shù)似然函數(shù)最大化,記為M-step;③ 重復(fù)步驟①和②,直到滿足指定的收斂準(zhǔn)則.
(15)
(16)
(17)
(18)
(19)
η的初值可由下式得到:
(20)
(21)
式中,xi為第i個觀測時刻的退化值;n為觀測次數(shù).
本節(jié)采用文獻(xiàn)[7]中某管道24 a的腐蝕深度檢測數(shù)據(jù)開展管道剩余壽命預(yù)測.所選管道材料為APILX52,管道內(nèi)壓為10 MPa,外徑為340.8 mm,管壁厚度為9.9 mm,最小極限應(yīng)力為359 MPa.表1為該管道不同檢測時刻(tc)的管道壁厚退化量.
表1 不同檢測時刻管道壁厚退化量
考慮到腐蝕管道失效與其剩余強(qiáng)度不足有關(guān),采用有限元方法仿真求解管道剩余強(qiáng)度.對腐蝕缺陷建模時需要將缺陷簡化為規(guī)則形狀,本文將腐蝕缺陷假設(shè)為長方體[26].缺陷的幾何尺寸包括最大腐蝕長度(lmax)、寬度(bmax)和深度(dmax).上述尺寸變化會影響等效應(yīng)力和管道剩余強(qiáng)度.參照文獻(xiàn)[7]設(shè)置管道的材料屬性參數(shù),分別建立9種不同的缺陷模型.通過有限元仿真模擬不同缺陷下管道最大等效應(yīng)力(τmax),結(jié)果如表2所示.
表2 不同缺陷所對應(yīng)的等效應(yīng)力
由表2可知,當(dāng)腐蝕長度不變時,等效應(yīng)力隨著腐蝕缺陷深度的增加而增加;當(dāng)腐蝕深度不變時,等效應(yīng)力將隨缺陷長度或?qū)挾鹊脑黾佣黾?因此,管道剩余強(qiáng)度受不同方向腐蝕缺陷大小的綜合作用,不能簡單地忽略某個方向的缺陷尺寸.
建立管道腐蝕缺陷有限元模型,劃分網(wǎng)格并設(shè)定邊界條件,可以計算出管道退化過程中不同時刻的等效應(yīng)力.圖1所示為缺陷模型對應(yīng)的應(yīng)力云圖.根據(jù)式(7),可計算得到不同檢測時刻的剩余強(qiáng)度退化量,如圖2所示.
圖1 管道應(yīng)力云圖
由式(8)可得到2個性能指標(biāo)剩余壽命的邊緣概率密度函數(shù),并根據(jù)AIC信息準(zhǔn)則選擇合適的Copula函數(shù)管道剩余壽命聯(lián)合概率密度函數(shù).
圖2 管道剩余強(qiáng)度性能指標(biāo)退化趨勢
按照AIC信息準(zhǔn)則得到4種常用Copula函數(shù)的計算結(jié)果如表3所示.根據(jù)AIC準(zhǔn)則計算所得R值越小,說明所對應(yīng)的Copula函數(shù)適應(yīng)性越好.由表3可知,Gaussian Copula函數(shù)所對應(yīng)的R值最小.因此,選擇Gaussian Copula函數(shù)來表征2個性能指標(biāo)之間的相關(guān)性,得到管道剩余壽命的聯(lián)合概率密度函數(shù).
表3 4種Copula函數(shù)的R值
圖3為剩余壁厚指標(biāo)所對應(yīng)的不同檢測時刻對參數(shù)μ、σ2和η的估計值.管道剩余強(qiáng)度參數(shù)所對應(yīng)的參數(shù)估計與剩余壁厚退化指標(biāo)趨勢相似,此處不再贅述.由參數(shù)估計結(jié)果可知,隨著有效退化數(shù)據(jù)量的增加,各參數(shù)估計結(jié)果逐漸逼近真實參數(shù)值.原因在于,貝葉斯更新與EM算法結(jié)合能夠充分利用新獲得的退化數(shù)據(jù),實時更新估計結(jié)果,有助于提高預(yù)測精度.上述結(jié)果驗證了本文所提出的參數(shù)估計方法的有效性.
(a) μ
(b) σ2
(c) η
表4 不同檢測時刻模型的參數(shù)估計值
圖4 不同時刻管道剩余壽命預(yù)測的概率密度函數(shù)
由圖4可知,隨著管道運行時間的增加,失效概率密度函數(shù)峰值呈增加趨勢,管道發(fā)生失效的概率不斷增大,同時管道的剩余壽命減小.上述結(jié)論與工程實際相符.
為客觀評價所提方法剩余壽命預(yù)測效果,定義相對誤差指標(biāo)(E):
(22)
式中,S為剩余壽命預(yù)測值;T為實際的剩余壽命.
管道實際壽命采用文獻(xiàn)[7]中的腐蝕油氣管道的剩余壽命值,即設(shè)計壽命減去當(dāng)前已服役的時間.剩余壽命預(yù)測誤差結(jié)果如表5所示.基于二元IG過程的剩余壽命預(yù)測方法的最大預(yù)測誤差為10.7%,最小預(yù)測誤差為2.2%,總體上,采用IG過程的預(yù)測誤差相對較小,具有較高的預(yù)測精度.
表5 剩余壽命預(yù)測值與實際值比較
1) 考慮到單一性能退化指標(biāo)難以全面反映輸油管道的健康狀態(tài),文中提出同時利用2個性能指標(biāo)評估管道健康狀態(tài)的方法,并構(gòu)建基于二元IG過程的管道性能退化模型.
2) 結(jié)合數(shù)值算例,采用Gaussian Copula函數(shù)評估2個性能指標(biāo)的相關(guān)性;通過更新IG過程確定管道失效概率密度函數(shù),分析管道健康狀態(tài),完成管道剩余壽命預(yù)測.案例分析表明,基于二元IG過程的剩余壽命預(yù)測最大預(yù)測誤差為10.7%,最小預(yù)測誤差為2.2%,所提出的預(yù)測方法具有較高的預(yù)測精度.