余建星， 任 杰， 楊政龍， 陳海成

(1. 天津大學(xué) 水利工程仿真與安全國家重點實驗室， 天津 300072；2. 高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心， 上海 200240)

0 引 言

深海采油樹系統(tǒng)是深水油氣田開發(fā)生產(chǎn)中極為重要的設(shè)備。深海的極端環(huán)境對采油樹系統(tǒng)的影響非常嚴(yán)重，極易引起泄漏事故，造成海域的大面積污染[1]。因此,采油樹系統(tǒng)的風(fēng)險評估迫在眉睫。

故障樹分析法規(guī)定了許多邏輯門和事件，并以圖形的方式表示出各事件之間的邏輯關(guān)系[2]。傳統(tǒng)的靜態(tài)故障樹基于的是靜態(tài)邏輯和靜態(tài)失效機理，因此其分析方法不適用于失效原因復(fù)雜且動態(tài)特征明顯的深海采油樹系統(tǒng)。對動態(tài)故障樹的定量分析方法國內(nèi)外學(xué)者已有較為深入的研究，如基于Markov模型的分析法[3]，Monte Carlo仿真法[4-5]，以及基于離散時間貝葉斯網(wǎng)絡(luò)的分析方法[6]，等等。但上述理論多應(yīng)用于化工及軟件領(lǐng)域，在深海結(jié)構(gòu)的風(fēng)險分析中應(yīng)用較少，而針對深海采油樹系統(tǒng)的風(fēng)險評估應(yīng)用更是空白。

因此，本文將基于動態(tài)邏輯門的動態(tài)故障樹理論引入深海采油樹系統(tǒng)的風(fēng)險評估，建立采油樹系統(tǒng)的動態(tài)故障樹模型，并根據(jù)現(xiàn)有的動態(tài)故障樹Monte Carlo定量分析方法，對動態(tài)故障樹進(jìn)行定量分析，計算頂事件故障概率并分析基本事件的概率重要度和結(jié)構(gòu)重要度系數(shù)，最后根據(jù)計算結(jié)果對采油樹系統(tǒng)提出風(fēng)險規(guī)避的建議。

1 深海采油樹系統(tǒng)動態(tài)故障樹建模

深海采油樹系統(tǒng)由采油樹主設(shè)備及其控制系統(tǒng)兩部分組成，其失效可分為采油樹樹體及其控制系統(tǒng)的失效。前者失效又可以分為3大類原因，即人為影響H1、外部影響H2、防腐失效H3。

1.1 人為影響因素

人為因素對采油樹系統(tǒng)的影響主要存在于制造過程H5和施工過程H6中。

H5的影響由部件生產(chǎn)失誤H13和部件檢驗失誤H40引起。H13可能來源于分析錯誤X2或者鋼料加工失誤H25，而H25可能由操作失誤X1、分析錯誤X2以及偶然碰撞X3引起。H40可能來源于X1以及X2。

H6的影響來源于下水過程H14、安裝過程H15以及儲運過程X6中的失誤。H14來源于X2和其他誤操作H26，H26包括X1以及X3。H15來源于零件鎖緊失效X4、密封失效H27、焊接失效H28以及H26，H27包括密封處有雜物X5、密封件損壞H41、密封處配合不完全H42，H41由X1或X3引起，H42由X2或X1引起[7]。H28由焊接施工和焊縫檢驗引起，焊接施工過程中的失誤來源于X1或X3，焊縫檢驗過程中的失誤來源于X1或X2。人為影響子樹如圖 1所示。

圖1 人為影響子樹及外部影響子樹

1.2 外部影響因素

外部影響因素分為采油樹體疲勞失效H7和突發(fā)沖擊載荷H8造成的失效。

H7由兩方面因素共同引發(fā)：首先是采油樹部件的內(nèi)因H16；其次是環(huán)境外因H17對采油樹體的影響。誘使H16事件發(fā)生的因素主要有安裝過程中擰緊力矩錯誤H29、產(chǎn)生了裂紋或尖銳缺口H30以及零件的形狀尺寸不合理H31，這些都可能由X1和X2引起。誘使H17發(fā)生的因素可能由渦激振動X7、渦旋影響H32、波浪影響H33以及X6引起。地形有狹窄通道X8和海水雷諾數(shù)過大X9同時發(fā)生會引發(fā)渦旋的破壞；而X9或風(fēng)力過大X10會引發(fā)波浪的破壞。

H8失效可能由以下因素引起：地震X15、錨擊X16、風(fēng)暴H18和土壤沖垮H19。H18由X10和采油樹部件結(jié)構(gòu)質(zhì)量過小X11共同引發(fā)，即在X10與 X11同時發(fā)生時，采油樹結(jié)構(gòu)才會失效。造成H19這一事件發(fā)生的因素可能有X9、土壤沉積物過厚X12、土壤黏性降低X13以及土壤液化X14。外部影響子樹如圖 1所示。

1.3 防腐失效

防腐失效的途徑有樹內(nèi)防腐失效H9以及樹外防腐失效H10。

H10主要來源于緩蝕劑變質(zhì)X22、陰極保護(hù)裝置H21以及內(nèi)壁防護(hù)涂層H22的失效。H22比較特殊，可以引入動態(tài)邏輯門中的優(yōu)先與門。以產(chǎn)生酸性氣體X20作為優(yōu)先事件，黏結(jié)剝離X21作為其次事件，即當(dāng)X20先發(fā)生，防護(hù)層X21后發(fā)生時，上層事件發(fā)生。H21主要來源于陰極屏蔽的發(fā)生H36或干擾/雜散電流的產(chǎn)生H37。H36由海管配重層內(nèi)加強筋與管道短接H45或保護(hù)套管與管道短接H46引起，前者來源于X15或X16的影響，后者還需加H19的影響。研究H37時同樣可以引入優(yōu)先與門，以絕緣層破壞X18作為優(yōu)先事件，其他管線的陰極保護(hù)系統(tǒng)失效X19作為其次事件，底事件先后發(fā)生時上層事件發(fā)生。

H9主要來源于H21以及防腐絕緣層H20的失效。前者的失效途徑與H9中陰極保護(hù)失效的途徑相同，后者的失效途徑主要有微生物的破壞X17、土壤的松動H35以及防腐涂層應(yīng)力開裂H34。H35可能由X12、X13以及X14引起，H34由H32或H33的破壞以及H8的破壞引起。防腐失效子樹如圖 2所示。

圖2 防腐失效子樹及控制系統(tǒng)失效子樹

1.4 控制系統(tǒng)失效

采油樹控制系統(tǒng)故障分為水上控制系統(tǒng)故障H11和水下控制系統(tǒng)故障H12。

H11包括液壓動力單元H23、電源單元X25、電子控制系統(tǒng)X26，其中任意一部分發(fā)生故障都會導(dǎo)致頂事件的發(fā)生。在研究H23系統(tǒng)故障時，引入熱備件門來描述其故障途徑。液壓泵系統(tǒng)失效X23作為其輸入事件，水上蓄能器失效X24作為其熱備件。只有二者均故障，液壓動力系統(tǒng)才會失效。

H12包括臍帶終端分配單元X27、水下控制單元H24和水下執(zhí)行器X33，其中任意一部分發(fā)生故障都會導(dǎo)致頂事件的發(fā)生。在研究H24的失效途徑時，同樣可以引入熱備件門來描述，即以水下控制單元的失效H38作為輸入事件，失效安全執(zhí)行單元H39作為其熱備件。只有工作元件和備用件均故障，水下控制系統(tǒng)才會失效。引起H38的因素可能有回油皮囊失效X28或電液換向閥X29的失效。引起失效安全執(zhí)行單元失效H39的因素可能有水下蓄能器X30、失效安全執(zhí)行器X31以及壓力補償器X32故障，三者任一故障都會引發(fā)上層事件的發(fā)生。控制系統(tǒng)失效子樹如圖 2所示。

2 Monte Carlo仿真法

目前動態(tài)故障樹定量分析的算法非常多，但是適用于深海結(jié)構(gòu)的并不多。一是由于深海結(jié)構(gòu)某些部件的故障概率隨時間變化，并不服從指數(shù)分布，因此不能使用傳統(tǒng)的Markov Chain法；二是深海結(jié)構(gòu)動態(tài)故障樹中的動態(tài)模塊可能比較龐大，會面臨著狀態(tài)空間組合爆炸的問題。此外，多重積分的計算也是難題。

針對該問題，李堂經(jīng)等[4]和楊恒占等[5]在動態(tài)故障樹的定量分析過程中引入Monte Carlo仿真法，將多重積分等價為具有特定分布隨機函數(shù)的期望值，設(shè)計隨機試驗以獲得隨機被積函數(shù)的樣本。Monte Carlo仿真法可以歸結(jié)為3個步驟：①構(gòu)造或描述概率過程；②實現(xiàn)從已知概率分布中抽樣；③建立各種估計量。

2.1 Monte Carlo仿真法動態(tài)子樹的失效分布

(1) 優(yōu)先與門

設(shè)其底事件x1和x2的發(fā)生時間為T1和T2，則其概率分布函數(shù)分別為G1(x1)和G2(x2)，那么頂事件發(fā)生時間的概率分布函數(shù)為GY(t)為

(1)

(2) 熱備件門

熱備件門主件與備用件同時運行且其失效分布函數(shù)基本相同；當(dāng)且僅當(dāng)二者都失效時頂事件發(fā)生。設(shè)主件x1和備用件x2的發(fā)生時間為T1和T2，其概率分布函數(shù)分別為G1(x1)和G2(x2)，則頂事件發(fā)生時間的概率分布函數(shù)GR(t)為

2.2 Monte Carlo仿真法的抽樣公式

考慮等式

(3)

對x求解即可得分布函數(shù)為F(x)的隨機變量的一個樣本觀察值。其中，ξ是由定義在[0,1]區(qū)間上的均勻分布u(x)產(chǎn)生的一個隨機數(shù)(使用計算機產(chǎn)生)[8]。

(1) 對于服從指數(shù)分布E(λ)的概率密度函數(shù)，若ξ是u(x)產(chǎn)生的一個隨機數(shù)，則t是E(λ)的一個樣本觀察值，又因為在[0,1]區(qū)間上1-ξ也服從均勻分布，所以表達(dá)式可化為

(4)

(2) 對于服從威布爾分布W(t0,m)的概率密度函數(shù)，t0與m分別為威布爾分布的尺度參數(shù)與形狀參數(shù)。式(3)的形式為1-e(-tm/t0)=ξ，則抽樣公式變?yōu)?/p>

(5)

3 深海采油樹系統(tǒng)動態(tài)故障樹的定量分析

CHATTERJEE[9]和BIRNBAUM等[10]拓寬了“模塊”這一概念的性質(zhì)并拓展了它在靜態(tài)故障樹中的應(yīng)用，然后通過基于TARJAN[11]算法的深度優(yōu)先搜索(Depth-First Left-Most，DFLM)，在動態(tài)故障樹中查找其中所有的動態(tài)子樹和靜態(tài)子樹。

通過MATLAB對圖 1和圖 2所示動態(tài)故障樹進(jìn)行兩次DFLM。據(jù)其遍歷結(jié)果可見，H4、H11、H12、H22、H23、H24、H37為動態(tài)模塊，H38、H39為靜態(tài)模塊。

設(shè)動態(tài)故障樹的底事件為X1,X2,…,Xn并進(jìn)行模塊化處理，則完整故障樹頂事件概率的解函數(shù)f可化為

f=f2(X1,X2,…,Xl,H4,H11,H12,H22,H23,H24,H37,H38,H39)

(6)

分別求解動態(tài)模塊H4、H11、H12、H22、H23、H24、H37以及靜態(tài)模塊H38、H39的頂事件概率，并將其作為新的動態(tài)故障樹f2的底事件概率；再對f2求解，即可知動態(tài)故障樹頂事件概率。

3.1 求解動態(tài)故障樹頂事件故障概率

3.1.1 求解模塊頂事件故障概率

查閱相關(guān)數(shù)據(jù)庫[12-14]，可得動態(tài)故障樹底事件故障率，如表 1所示。

表1 動態(tài)故障樹底事件故障率 106

圖3 H38、H39模塊BDD圖

由于X23～X33均為機械電子元件，因此，可設(shè)其均服從威布爾分布W(100 000,3)，工作時間為t=104。

動態(tài)故障樹中H38、H39均為靜態(tài)模塊，可用二元決策圖(Binary Decision Diagram， BDD)法求解，其BDD圖如圖 3所示。

由此可得H38模塊的最小割集為{X28}、{X29}，H39模塊的最小割集為{X30}、{X31}、{X32}。因此,可得

(7)

(8)

H22、H23、H24、H37均為動態(tài)模塊，使用Monte Carlo仿真法進(jìn)行定量分析求解模塊頂事件故障概率。其中，H23、H24模塊服從動態(tài)邏輯門中的熱備件門。由于X23、X24以及X28～X32服從威布爾分布，因此，使用式(5)進(jìn)行隨機抽樣。

將產(chǎn)生的隨機數(shù)代入式(2)，利用MATLAB可實現(xiàn)循環(huán)計算過程，可得

P{H23}=7×10-5

P{H24}=1.66×10-3

動態(tài)模塊中H22、H37服從動態(tài)邏輯門中的優(yōu)先與門。由于X18～X21均服從指數(shù)分布，因此使用式(4)進(jìn)行隨機抽樣。將產(chǎn)生的隨機數(shù)代入式(1)可得

P{H22}=2.029 8×10-4

P{H37}=3.720 9×10-4

3.1.2 求解動態(tài)故障樹最小割集

對于完整的動態(tài)故障樹，可利用下行法求解其最小割集。(以下以Hi代替PHi，Xi代替PXi)

對于H1模塊，由深海采油樹系統(tǒng)動態(tài)故障樹可得

H5=H40·H13=(X1+X2)·(H25+X2)=(X1+X2)·(X1+X2+X3+X2)

=X1+X2+X1·X2+X2·X3+X1·X3

(9)

H6=H14+H15+X6=(H26+X2)+(H26+H27+H28+X4)+X6

=X1+X2+X3+X4+X5+X6+X1·X2+X2·X3+X1·X3

(10)

由式(9)和式(10)，根據(jù)布爾運算法則，可得

H1=H5+H6=X1+X2+X3+X4+X5+X6+X1·X2+X2·X3+X1·X3

(11)

對于H2模塊，根據(jù)布爾運算法則，同理可得

H2=H7+H8=X9+X12+X13+X14+X15+X16+X1·X6+X1·X7+X1·X9+X1·X10+

X2·X6+X2·X7+X2·X9+X2·X10+X10·X11+X1·X8·X9+X2·X8·X9

(12)

對于H3模塊，同理有

H3=H9+H10=X9+X10+X12+X13+X14+X15+X16+X17+X22+

H22+H37+X8·X9+X10·X11

(13)

對于H4模塊，同理有

H4=H11+H12=H23+H24+X25+X26+X27+X33

(14)

3.1.3 計算頂事件故障概率

深海采油樹系統(tǒng)動態(tài)故障樹的最小割集：

T=H1+H2+H3+H4

(15)

將式(11)～式(14)代入式(15)中可得

p(T)=0.006 3。

(16)

3.2 求解動態(tài)故障樹基本事件概率重要度和結(jié)構(gòu)重要度系數(shù)

概率重要度系數(shù)的計算公式[15]為

(17)

式中：p(T)為頂事件發(fā)生概率;qi為第i個基本事件發(fā)生的概率。

概率重要度的重要性質(zhì)：設(shè)所有基本事件發(fā)生概率均為1/2，則結(jié)構(gòu)重要度系數(shù)等于此時的概率重要度系數(shù)，即

(18)

利用該性質(zhì)可求解結(jié)構(gòu)重要度系數(shù)。

3.2.1 求解基本事件概率重要度

根據(jù)式(17)計算一次偏導(dǎo)數(shù)，并對基本事件概率重要度排序，如表 2所示。

表2 基本事件概率重要度排序 10-4

由表 2可以看出：基本事件X23～X33在基本事件概率重要度排序中最高，即減小這些基本事件的故障率能迅速減小頂事件的發(fā)生概率，而它們都是機械電子元件，因此控制系統(tǒng)需要定期檢查、保養(yǎng)，及時更換，并使用高質(zhì)量的電子元件；基本事件X2、X1在概率重要度排序中緊隨電子元件之后，重要度也不可忽視，因此減小人為操作失誤和分析失誤的概率對減小頂事件發(fā)生概率同樣十分重要，可以通過聘請高素質(zhì)人才來設(shè)計、操作某些重要部件以減少此類失誤。

3.2.2 求解基本事件結(jié)構(gòu)重要度系數(shù)

利用概率重要度的性質(zhì)，可求得其結(jié)構(gòu)重要度系數(shù)，如表 3所示。

表3 基本事件結(jié)構(gòu)重要度系數(shù)排序 10-4

由表 3可知：基本事件X1和X2的結(jié)構(gòu)重要度系數(shù)最高，與第3位差距很大，因此減小人為操作失誤和分析錯誤可以規(guī)避很多風(fēng)險，減小頂事件發(fā)生概率；基本事件X9、X10、X8的結(jié)構(gòu)重要度系數(shù)緊隨其后，因此盡量規(guī)避易產(chǎn)生較大渦旋和波浪的海域?qū)p小頂事件發(fā)生概率有重要作用；基本事件X3、X6、X7的結(jié)構(gòu)重要度系數(shù)相同，處于排序表的中間位置，因此避免偶然碰撞、渦激振動以及儲運過程中的不利因素同樣對減小頂事件發(fā)生概率有積極作用。

綜上所述，對采油樹系統(tǒng)部件進(jìn)行定期檢查、保養(yǎng)，發(fā)生故障以后及時更換，或投入較大人力、財力避免人為失誤以及減少可控故障，是規(guī)避風(fēng)險的最佳選擇。

4 結(jié) 語

建立深海采油樹系統(tǒng)的動態(tài)故障樹模型，根據(jù)Monte Carlo仿真法對該模型進(jìn)行了定量分析，計算頂事件故障概率并分析了基本事件的概率重要度和結(jié)構(gòu)重要度系數(shù)，最后根據(jù)計算結(jié)果對采油樹系統(tǒng)提出了風(fēng)險規(guī)避的建議。

后續(xù)還需繼續(xù)研究的方面主要有： 深海采油樹系統(tǒng)動態(tài)故障樹模型的底事件發(fā)生概率仍然靠經(jīng)驗數(shù)據(jù)得來，未來可進(jìn)一步通過調(diào)研及試驗確定其概率。

現(xiàn)有的Monte Carlo仿真法依舊沒有解決狀態(tài)空間組合爆炸的問題，未來可開發(fā)新的算法以使計算更精確。