余建星, 任 杰, 楊政龍, 陳海成
(1. 天津大學 水利工程仿真與安全國家重點實驗室, 天津 300072;2. 高新船舶與深海開發裝備協同創新中心, 上海 200240)
深海采油樹系統是深水油氣田開發生產中極為重要的設備。深海的極端環境對采油樹系統的影響非常嚴重,極易引起泄漏事故,造成海域的大面積污染[1]。因此,采油樹系統的風險評估迫在眉睫。
故障樹分析法規定了許多邏輯門和事件,并以圖形的方式表示出各事件之間的邏輯關系[2]。傳統的靜態故障樹基于的是靜態邏輯和靜態失效機理,因此其分析方法不適用于失效原因復雜且動態特征明顯的深海采油樹系統。對動態故障樹的定量分析方法國內外學者已有較為深入的研究,如基于Markov模型的分析法[3],Monte Carlo仿真法[4-5],以及基于離散時間貝葉斯網絡的分析方法[6],等等。但上述理論多應用于化工及軟件領域,在深海結構的風險分析中應用較少,而針對深海采油樹系統的風險評估應用更是空白。
因此,本文將基于動態邏輯門的動態故障樹理論引入深海采油樹系統的風險評估,建立采油樹系統的動態故障樹模型,并根據現有的動態故障樹Monte Carlo定量分析方法,對動態故障樹進行定量分析,計算頂事件故障概率并分析基本事件的概率重要度和結構重要度系數,最后根據計算結果對采油樹系統提出風險規避的建議。
深海采油樹系統由采油樹主設備及其控制系統兩部分組成,其失效可分為采油樹樹體及其控制系統的失效。前者失效又可以分為3大類原因,即人為影響H1、外部影響H2、防腐失效H3。
人為因素對采油樹系統的影響主要存在于制造過程H5和施工過程H6中。
H5的影響由部件生產失誤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 人為影響子樹及外部影響子樹
外部影響因素分為采油樹體疲勞失效H7和突發沖擊載荷H8造成的失效。
H7由兩方面因素共同引發:首先是采油樹部件的內因H16;其次是環境外因H17對采油樹體的影響。誘使H16事件發生的因素主要有安裝過程中擰緊力矩錯誤H29、產生了裂紋或尖銳缺口H30以及零件的形狀尺寸不合理H31,這些都可能由X1和X2引起。誘使H17發生的因素可能由渦激振動X7、渦旋影響H32、波浪影響H33以及X6引起。地形有狹窄通道X8和海水雷諾數過大X9同時發生會引發渦旋的破壞;而X9或風力過大X10會引發波浪的破壞。
H8失效可能由以下因素引起:地震X15、錨擊X16、風暴H18和土壤沖垮H19。H18由X10和采油樹部件結構質量過小X11共同引發,即在X10與 X11同時發生時,采油樹結構才會失效。造成H19這一事件發生的因素可能有X9、土壤沉積物過厚X12、土壤黏性降低X13以及土壤液化X14。外部影響子樹如圖 1所示。
防腐失效的途徑有樹內防腐失效H9以及樹外防腐失效H10。
H10主要來源于緩蝕劑變質X22、陰極保護裝置H21以及內壁防護涂層H22的失效。H22比較特殊,可以引入動態邏輯門中的優先與門。以產生酸性氣體X20作為優先事件,黏結剝離X21作為其次事件,即當X20先發生,防護層X21后發生時,上層事件發生。H21主要來源于陰極屏蔽的發生H36或干擾/雜散電流的產生H37。H36由海管配重層內加強筋與管道短接H45或保護套管與管道短接H46引起,前者來源于X15或X16的影響,后者還需加H19的影響。研究H37時同樣可以引入優先與門,以絕緣層破壞X18作為優先事件,其他管線的陰極保護系統失效X19作為其次事件,底事件先后發生時上層事件發生。
H9主要來源于H21以及防腐絕緣層H20的失效。前者的失效途徑與H9中陰極保護失效的途徑相同,后者的失效途徑主要有微生物的破壞X17、土壤的松動H35以及防腐涂層應力開裂H34。H35可能由X12、X13以及X14引起,H34由H32或H33的破壞以及H8的破壞引起。防腐失效子樹如圖 2所示。

圖2 防腐失效子樹及控制系統失效子樹
采油樹控制系統故障分為水上控制系統故障H11和水下控制系統故障H12。
H11包括液壓動力單元H23、電源單元X25、電子控制系統X26,其中任意一部分發生故障都會導致頂事件的發生。在研究H23系統故障時,引入熱備件門來描述其故障途徑。液壓泵系統失效X23作為其輸入事件,水上蓄能器失效X24作為其熱備件。只有二者均故障,液壓動力系統才會失效。
H12包括臍帶終端分配單元X27、水下控制單元H24和水下執行器X33,其中任意一部分發生故障都會導致頂事件的發生。在研究H24的失效途徑時,同樣可以引入熱備件門來描述,即以水下控制單元的失效H38作為輸入事件,失效安全執行單元H39作為其熱備件。只有工作元件和備用件均故障,水下控制系統才會失效。引起H38的因素可能有回油皮囊失效X28或電液換向閥X29的失效。引起失效安全執行單元失效H39的因素可能有水下蓄能器X30、失效安全執行器X31以及壓力補償器X32故障,三者任一故障都會引發上層事件的發生??刂葡到y失效子樹如圖 2所示。
目前動態故障樹定量分析的算法非常多,但是適用于深海結構的并不多。一是由于深海結構某些部件的故障概率隨時間變化,并不服從指數分布,因此不能使用傳統的Markov Chain法;二是深海結構動態故障樹中的動態模塊可能比較龐大,會面臨著狀態空間組合爆炸的問題。此外,多重積分的計算也是難題。
針對該問題,李堂經等[4]和楊恒占等[5]在動態故障樹的定量分析過程中引入Monte Carlo仿真法,將多重積分等價為具有特定分布隨機函數的期望值,設計隨機試驗以獲得隨機被積函數的樣本。Monte Carlo仿真法可以歸結為3個步驟:①構造或描述概率過程;②實現從已知概率分布中抽樣;③建立各種估計量。
(1) 優先與門
設其底事件x1和x2的發生時間為T1和T2,則其概率分布函數分別為G1(x1)和G2(x2),那么頂事件發生時間的概率分布函數為GY(t)為
(1)
(2) 熱備件門
熱備件門主件與備用件同時運行且其失效分布函數基本相同;當且僅當二者都失效時頂事件發生。設主件x1和備用件x2的發生時間為T1和T2,其概率分布函數分別為G1(x1)和G2(x2),則頂事件發生時間的概率分布函數GR(t)為
考慮等式

(3)
對x求解即可得分布函數為F(x)的隨機變量的一個樣本觀察值。其中,ξ是由定義在[0,1]區間上的均勻分布u(x)產生的一個隨機數(使用計算機產生)[8]。
(1) 對于服從指數分布E(λ)的概率密度函數,若ξ是u(x)產生的一個隨機數,則t是E(λ)的一個樣本觀察值,又因為在[0,1]區間上1-ξ也服從均勻分布,所以表達式可化為
(4)
(2) 對于服從威布爾分布W(t0,m)的概率密度函數,t0與m分別為威布爾分布的尺度參數與形狀參數。式(3)的形式為1-e(-tm/t0)=ξ,則抽樣公式變為
(5)
CHATTERJEE[9]和BIRNBAUM等[10]拓寬了“模塊”這一概念的性質并拓展了它在靜態故障樹中的應用,然后通過基于TARJAN[11]算法的深度優先搜索(Depth-First Left-Most,DFLM),在動態故障樹中查找其中所有的動態子樹和靜態子樹。
通過MATLAB對圖 1和圖 2所示動態故障樹進行兩次DFLM。據其遍歷結果可見,H4、H11、H12、H22、H23、H24、H37為動態模塊,H38、H39為靜態模塊。
設動態故障樹的底事件為X1,X2,…,Xn并進行模塊化處理,則完整故障樹頂事件概率的解函數f可化為
f=f2(X1,X2,…,Xl,H4,H11,H12,H22,H23,H24,H37,H38,H39)
(6)
分別求解動態模塊H4、H11、H12、H22、H23、H24、H37以及靜態模塊H38、H39的頂事件概率,并將其作為新的動態故障樹f2的底事件概率;再對f2求解,即可知動態故障樹頂事件概率。
3.1.1 求解模塊頂事件故障概率
查閱相關數據庫[12-14],可得動態故障樹底事件故障率,如表 1所示。

表1 動態故障樹底事件故障率 106

圖3 H38、H39模塊BDD圖
由于X23~X33均為機械電子元件,因此,可設其均服從威布爾分布W(100 000,3),工作時間為t=104。
動態故障樹中H38、H39均為靜態模塊,可用二元決策圖(Binary Decision Diagram, BDD)法求解,其BDD圖如圖 3所示。
由此可得H38模塊的最小割集為{X28}、{X29},H39模塊的最小割集為{X30}、{X31}、{X32}。因此,可得
(7)

(8)
H22、H23、H24、H37均為動態模塊,使用Monte Carlo仿真法進行定量分析求解模塊頂事件故障概率。其中,H23、H24模塊服從動態邏輯門中的熱備件門。由于X23、X24以及X28~X32服從威布爾分布,因此,使用式(5)進行隨機抽樣。
將產生的隨機數代入式(2),利用MATLAB可實現循環計算過程,可得
P{H23}=7×10-5
P{H24}=1.66×10-3
動態模塊中H22、H37服從動態邏輯門中的優先與門。由于X18~X21均服從指數分布,因此使用式(4)進行隨機抽樣。將產生的隨機數代入式(1)可得
P{H22}=2.029 8×10-4
P{H37}=3.720 9×10-4
3.1.2 求解動態故障樹最小割集
對于完整的動態故障樹,可利用下行法求解其最小割集。(以下以Hi代替PHi,Xi代替PXi)
對于H1模塊,由深海采油樹系統動態故障樹可得
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),根據布爾運算法則,可得
H1=H5+H6=X1+X2+X3+X4+X5+X6+X1·X2+X2·X3+X1·X3
(11)
對于H2模塊,根據布爾運算法則,同理可得
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=H1+H2+H3+H4
(15)
將式(11)~式(14)代入式(15)中可得
p(T)=0.006 3。
(16)
概率重要度系數的計算公式[15]為
(17)
式中:p(T)為頂事件發生概率;qi為第i個基本事件發生的概率。
概率重要度的重要性質:設所有基本事件發生概率均為1/2,則結構重要度系數等于此時的概率重要度系數,即
(18)
利用該性質可求解結構重要度系數。
3.2.1 求解基本事件概率重要度
根據式(17)計算一次偏導數,并對基本事件概率重要度排序,如表 2所示。

表2 基本事件概率重要度排序 10-4
由表 2可以看出:基本事件X23~X33在基本事件概率重要度排序中最高,即減小這些基本事件的故障率能迅速減小頂事件的發生概率,而它們都是機械電子元件,因此控制系統需要定期檢查、保養,及時更換,并使用高質量的電子元件;基本事件X2、X1在概率重要度排序中緊隨電子元件之后,重要度也不可忽視,因此減小人為操作失誤和分析失誤的概率對減小頂事件發生概率同樣十分重要,可以通過聘請高素質人才來設計、操作某些重要部件以減少此類失誤。
3.2.2 求解基本事件結構重要度系數
利用概率重要度的性質,可求得其結構重要度系數,如表 3所示。

表3 基本事件結構重要度系數排序 10-4
由表 3可知:基本事件X1和X2的結構重要度系數最高,與第3位差距很大,因此減小人為操作失誤和分析錯誤可以規避很多風險,減小頂事件發生概率;基本事件X9、X10、X8的結構重要度系數緊隨其后,因此盡量規避易產生較大渦旋和波浪的海域對減小頂事件發生概率有重要作用;基本事件X3、X6、X7的結構重要度系數相同,處于排序表的中間位置,因此避免偶然碰撞、渦激振動以及儲運過程中的不利因素同樣對減小頂事件發生概率有積極作用。
綜上所述,對采油樹系統部件進行定期檢查、保養,發生故障以后及時更換,或投入較大人力、財力避免人為失誤以及減少可控故障,是規避風險的最佳選擇。
建立深海采油樹系統的動態故障樹模型,根據Monte Carlo仿真法對該模型進行了定量分析,計算頂事件故障概率并分析了基本事件的概率重要度和結構重要度系數,最后根據計算結果對采油樹系統提出了風險規避的建議。
后續還需繼續研究的方面主要有: 深海采油樹系統動態故障樹模型的底事件發生概率仍然靠經驗數據得來,未來可進一步通過調研及試驗確定其概率。
現有的Monte Carlo仿真法依舊沒有解決狀態空間組合爆炸的問題,未來可開發新的算法以使計算更精確。