侯 磊, 吳守志, 劉芳媛, 伍星光

(中國石油大學油氣管道輸送安全國家工程實驗室,石油工程教育部重點實驗室,北京 102249)

隨著能源化工行業(yè)設備多樣復雜化、生產(chǎn)技術更新和危險品數(shù)量增加,發(fā)生二次及更高次事故可能性大大增加,潛在多米諾事故危險性也急劇增加[1]。多米諾效應具有事故小概率性、擴展不確定性等特點,會造成災難性后果。歐洲從1982年開始對多米諾危害進行評估,經(jīng)過多年研究,學者們在多米諾效應頻率估計[2]、輻射引起的概率估計[3]、事故模型[4-6]等方面取得重大進展。Cozzani等[7]對危險源進行全面分析,研發(fā)了一種計算多米諾事件概率的方法,通過初始事件、升級向量和目標單元識別,計算事故擴展概率和二次事件概率,稱為解析法,該方法只考慮一級多米諾效應,忽略二次場景可能引起的進一步升級。Abdolhamidzadeh等[8]對224起多米諾事故進行分析,發(fā)現(xiàn)超過一半的事故都引發(fā)二次及更高級別事件,認為解析法忽略高級別多米諾效應是不合理的,會造成計算結(jié)果失真。因此,Abdolhamidzadeh等[9]提出蒙特卡洛模擬方法,其在處理多米諾效應不確定性和復雜性上有很大優(yōu)勢。為確定多米諾效應最可能的事故序列,Khakzad等[10]提出一種基于貝葉斯網(wǎng)絡的多米諾效應傳播和概率計算方法,該方法考慮事故擴展不確定性、設備間協(xié)同和相互作用,能根據(jù)新信息更新事件概率確定事故傳播模式和各級多米諾事件概率。目前,這三種方法得到了不同程度的應用,但對油庫多米諾效應,沒有查閱到方法優(yōu)劣性和適用性的研究報道。筆者從這3種計算方法的理論基礎、計算過程出發(fā),分析并對比其優(yōu)缺點,結(jié)合實際案例計算,確定考慮事故傳播更全面、事故概率更準確且事故擴展更符合實際的定量分析方法,用于油庫多米諾效應預防和控制。

1 油庫池火災多米諾效應與安全屏障

1.1 油庫池火災多米諾效應概述

自1907年發(fā)生第一起多米諾事故,學者們相繼提出多米諾效應的定義,目前還沒有公認的定義[11]。Reniners等[12]結(jié)合前人研究,提出一個較全面的定義:一個初始事件傳播到附近設備,觸發(fā)一個或多個二次事件,進而引發(fā)高階事件,導致比初始事件后果更嚴重的事故。多米諾效應的特征[7]為初始事件、升級向量和目標單元,升級向量包括熱輻射、沖擊波、碎片和有毒氣體擴散。

油庫池火災多米諾效應是指一個儲罐發(fā)生池火災,以熱輻射傳播到附近儲罐,觸發(fā)一個或多個儲罐發(fā)生二次事故,或觸發(fā)更高級別事故,導致比單罐池火災更嚴重的罐區(qū)事故。常壓儲罐是中國大型原油儲庫中常用儲罐,一旦單罐池火災引發(fā)多米諾效應和群罐火災,會造成災難性后果。常壓儲罐設置有安全屏障,對多米諾效應預防和控制有重要作用[13]。

1.2 油庫池火災安全屏障

根據(jù)相關研究,安全屏障分為4類[14]:①本質(zhì)安全設計;②主動防護屏障;③被動防護屏障;④應急和響應措施。因儲罐及相關設施已建成,設計、安全距離等都已確定,應急和響應措施受管理、人為影響較大,不確定性較多。本研究對主動和被動防護屏障下多米諾效應概率計算。

1.2.1 主動防護屏障

主動防護屏障需要外界提供能量或操作才能起作用,如噴淋水冷卻裝置、泡沫滅火系統(tǒng)等,用于控制和撲滅火災,通常安裝在儲罐上部。

噴淋水冷卻裝置通過在儲罐外表面噴淋水流為儲罐冷卻降溫,降低儲罐表面的熱輻射強度。噴淋水裝置對儲罐各部分冷卻降溫不均勻,為簡化計算,根據(jù)相關文獻[13]引入強度減弱系數(shù)為0.8。

當噴淋水冷卻裝置起作用時,儲罐表面的熱輻射強度計算為

QWDS=αQHL.

(1)

式中,QHL為儲罐未冷卻時接受到的熱輻射強度,kW/m2;QWDS為噴淋水起作用時儲罐接受到的熱輻射強度,kW/m2;α為強度減弱系數(shù),取0.8。

1.2.2 被動防護屏障

被動防護屏障安裝在儲罐外壁,能減緩熱輻射下儲罐壁溫升速度[15],延長儲罐壁因受熱而發(fā)生應力失效時間,即失效時間。

一般原油儲罐的被動防護屏障是罐壁板外側(cè)隔熱涂層,Landucci等[16]通過數(shù)值模擬求取有隔熱涂層時儲罐失效時間,并擬合了熱輻射強度和失效時間的經(jīng)驗公式:

ttf=0.0167exp(-1.13lnI-2.67×10-5V+9.9).

(2)

式中,ttf為儲罐未考慮安全屏障的失效時間,min;I為儲罐表面的熱輻射強度,kW/m2;V為儲罐容量,m3,25 m3≤V≤17 500 m3。

考慮隔熱涂層時失效時間為

tp=ttf+ttfc.

(3)

式中,tp為最終失效時間,min;ttfc為隔熱涂層增加時間,取15 min[17]。

常壓儲罐事故升級用Probit模型將升級向量轉(zhuǎn)化為事故擴展概率PEscalation。采用Probit模型計算得到概率單位變量Pr,PEscalation與Pr服從正態(tài)分布,采用概率函數(shù)法計算PEscalation,表示為

Pr=9.261-1.85lntp,

(4)

(5)

式中,u為積分變量。

2 多米諾效應定量分析方法

由油庫池火災觸發(fā)的多米諾效應會引起一系列連鎖反應,造成多個二次事件,對整個罐區(qū)風險有較大影響。計算并分析池火災引起的多米諾效應概率,為降低事故后果和風險提供依據(jù)。

2.1 解析法

Cozzani等提出的解析法[8],從識別危險源開始,考慮初始場景及頻率,確定升級向量,確定各儲罐可能的事故場景和概率。該方法只考慮一級多米諾效應,忽略不同初始事件協(xié)同效應。

由初始事件引發(fā)的多米諾事件概率為

fde=fpePd.

(6)

式中,fde為多米諾事件概率;fpe為初始事件概率;Pd為升級概率。

假設目標區(qū)域共有N+1個設備,選擇一個作為初始事件設備,則k(k≤N)個設備同時發(fā)生事故第m種多米諾場景概率為

(7)

k個發(fā)生二次事件設備同時發(fā)生第m種多米諾場景概率為

(8)

解析法計算過程見圖1。

解析法能夠計算升級事件概率,分析事故后果較全面,在分析初始事件引發(fā)的一級多米諾效應并忽略協(xié)同效應和相互作用時有較大優(yōu)勢。

2.2 蒙特卡洛模擬

針對解析法在處理多米諾效應不確定性和復雜性的局限性,Abdolhamldzadeh[9]提出蒙特卡洛模擬,輸入設備數(shù)目、失效頻率、擴展概率、迭代次數(shù)等參數(shù),對設備失效引發(fā)多米諾效應概率進行模擬,并確定概率。實現(xiàn)過程見圖2。

該方法不受模擬系統(tǒng)大小和復雜程度限制,只考慮事件間擴展概率,用隨機數(shù)理論得出多米諾效應概率。計算結(jié)果比較準確,對不確定性輸入?yún)?shù)和多個設備多米諾效應適用性較強。

圖1 解析法計算過程Fig.1 Analytical method calculation process

圖2 蒙特卡洛模擬過程Fig.2 Monte carlo simulation process

2.3 貝葉斯網(wǎng)絡

貝葉斯網(wǎng)絡是一種結(jié)合圖論和概率論研究不確定性問題的方法。貝葉斯網(wǎng)絡中節(jié)點表示隨機變量,通過有向弧連接,弧表示連接節(jié)點間依賴關系,分配給節(jié)點的條件概率表決定依賴關系類型和強度?；∫龅墓?jié)點稱為父節(jié)點,弧指向的節(jié)點稱為子節(jié)點,一個節(jié)點可以既是子節(jié)點又是父節(jié)點。沒有父節(jié)點或子節(jié)點的節(jié)點叫做根節(jié)點或葉節(jié)點[8]。

利用鏈式法則和D-分離準則,貝葉斯網(wǎng)絡擴展了一組鏈接節(jié)點聯(lián)合概率分布,如U=(X1,X2,…,Xn),其概率計算為

(9)

式中，P(U)為變量的聯(lián)合概率分布;Pa(Xi)為變量Xi的父節(jié)點。

貝葉斯網(wǎng)絡推理是在給定證據(jù)節(jié)點時計算其他節(jié)點后驗概率。在N個節(jié)點的貝葉斯網(wǎng)絡V=(V1,V2,…,VN)中,給定證據(jù)ε=e,條件概率P(Vi=vε=e)計算公式為

(10)

貝葉斯網(wǎng)絡建模計算過程見圖3。

圖3 貝葉斯網(wǎng)絡建模計算過程Fig.3 Bayesian network modeling calculation process

貝葉斯網(wǎng)絡能夠考慮事故多級多米諾效應、事故間協(xié)同效應及多種事故場景,使對多米諾效應概率分析更準確、傳播過程更完整,并確定最可能的事故場景和序列。目前其在管道腐蝕、泄漏、火災爆炸[18]等方面有一些應用。但因不確定性和證據(jù)推理復雜性,很難集成到其他軟件中使用。

3 案例分析

解析法、蒙特卡洛模擬和貝葉斯網(wǎng)絡是多米諾效應定量分析的常用方法,為掌握其各自適用性,分別計算有無安全屏障下多米諾概率,確定更準確、更符合實際的多米諾效應概率計算方法。

3.1 案例一

某罐區(qū)有6個常壓外浮頂式原油儲罐,容量均為1×104m3,罐壁高度為15.08 m,儲罐內(nèi)徑為28.5 m,兩相鄰儲罐間距為12 m,存儲相同性質(zhì)原油,原油燃燒熱為41 800 kJ/kg,蒸發(fā)熱為400 kJ/kg,比熱容為2.4 kJ/(kg·℃),20 ℃下密度為810 kg/m3,罐區(qū)平面布置見圖4。

圖4 罐區(qū)平面布局圖Fig.4 Tank farm plot layout

假定初始事故為儲罐T1泄漏油品被點燃,引發(fā)浮頂全表面池火災,池火災概率取1×10-5[8],事故升級臨界閾值為15 kW/m2。

3.1.1 解析法計算事故概率

對T2、T4、T5儲罐,對比點源模型[19]、Shokri-Beyler模型及Mudan模型等常見池火災輻射模型[20-22],綜合考慮大型儲罐布局和熱輻射傷害范圍,采用Mudan模型[23-24]計算儲罐表面受到的熱輻射強度I。由式(3)、(4)計算儲罐失效時間,再由式(5)計算Pr,最后由式(1)計算PEscalation,結(jié)果見表1。

表1 各儲罐接受熱輻射強度與擴展概率Table 1 Heat radiation intensity and escalation probability

其余各儲罐事故概率為P2=P4=5.859×10-6,P3=3.719×10-6,P5=3.464×10-6,P6=3.145×10-6。

3.1.2 蒙特卡洛模擬事故概率

根據(jù)蒙特卡洛模擬基本思路,輸入相應的迭代矩陣,分別在1 000、1×104、10×104、100×104次迭代下模擬3次,模擬結(jié)果見表2。

隨蒙特卡洛模擬迭代和模擬次數(shù)增加,模擬結(jié)果有一定波動,總體是逐漸穩(wěn)定的。迭代次數(shù)越多,模擬時間呈倍數(shù)型增長,說明精度越高迭代時間越長。從精度上看,選取100×104次迭代3次模擬結(jié)果平均值作為最終結(jié)果。

表2 蒙特卡洛模擬結(jié)果

3.1.3 貝葉斯網(wǎng)絡計算事故概率

運用貝葉斯網(wǎng)絡軟件GeNIe繪制6儲罐區(qū)網(wǎng)絡圖,輸入各節(jié)點概率參數(shù),考慮初始事故引發(fā)的一、二級多米諾效應。

L1、L2為或門,表示一、二級多米諾單元是否失效,其中

L1=T2∪T4∪T5,

(11)

L2=T3∪T6.

(12)

DL1、DL2為與門,表示一、二多米諾效應是否發(fā)生,其中

DL1=T1∩L1,

(13)

DL2=L1∩L2.

(14)

貝葉斯網(wǎng)絡下多米諾傳播模式見圖5,各儲罐事故概率按式(6)～(8)計算。3種方法計算結(jié)果見表3。

圖5 T1發(fā)生事故時貝葉斯法多米諾效應傳播模式Fig.5 Propagation pattern of domino effect by Bayesian method when T1 accident

表3 T1發(fā)生事故時貝葉斯法3種方法儲罐事故概率結(jié)果對比

當T2儲罐發(fā)生池火災時,此時只引發(fā)一級多米諾效應,3種方法計算過程參照前面。

解析法計算各儲罐事故概率為

P1=P3=P5=5.859×10-6,P4=P6=3.464×10-6.

蒙特卡洛模擬各儲罐事故概率為

P1=5.717×10-6,P2=9.86×10-6,P3=5.904×10-6,

P4=4.329×10-6,P5=5.329×10-6,P6=3.626×10-6.

貝葉斯網(wǎng)絡下多米諾傳播模式見圖6,各儲罐事故概率按式(6)～(8)計算。此時3種方法計算結(jié)果見表4。

圖6 T2發(fā)生事故時貝葉斯網(wǎng)絡多米諾效應傳播模式Fig.6 Propagation pattern of domino effect by Bayesian method when T2 accident

表4 T2發(fā)生事故時3種方法儲罐事故概率結(jié)果對比

3.2 案例二

本案例對安全屏障下事故概率進行分析,因蒙特卡洛模擬無法設置安全屏障,只對解析法和貝葉斯網(wǎng)絡下T1、T2計算結(jié)果進行分析。

Landucci[14]提出一種量化安全屏障性能的方法,引入可用性和有效性,安全屏障為噴淋水冷卻裝置和隔熱涂層,其相關參數(shù)見表5。

表5 安全屏障有效性和可用性

根據(jù)安全屏障相關研究[17],隔熱涂層起作用時將失效時間延長15 min,噴淋水裝置起作用時將儲罐表面接收熱輻射強度降低至80%。在不同安全屏障下相關參數(shù)計算結(jié)果見表6。

通過解析法事件樹和貝葉斯網(wǎng)絡,分別對隔熱涂層、噴淋水冷卻裝置與隔熱涂層下多米諾效應概率進行分析,對比結(jié)果及操作情況。

表6 不同安全屏障下事故擴展概率Table 6 Accident escalation probability under different safety barriers

(1)設置隔熱涂層時,事故概率計算結(jié)果見圖7和表7。

圖7 設置隔熱涂層事故概率計算Fig.7 Probability calculation affected by considering fireproofing coating

表7 設置隔熱涂層事故概率Table 7 Accident probability affected by considering fireproofing coating

(2)設置噴淋水冷卻裝置和隔熱涂層時,事故概率計算結(jié)果見表8、圖8。

表8 設置兩層安全屏障事故概率Table 8 Accident probability affected by considering two layers of safety barriers

圖8 設置兩層安全屏障事故概率計算Fig.8 Probability calculation affected by considering two layers of safety barriers

4 結(jié)果分析

4.1 案例一

對一級多米諾單元,解析法和貝葉斯網(wǎng)絡結(jié)果均為5.859×10-6,計算結(jié)果較準確,與解析法適用于一級多米諾效應相契合;但對二級多米諾單元T3、T6,解析法結(jié)果分別為3.719×10-6、3.145×10-6,與貝葉斯網(wǎng)絡的3.444×10-6、2.991×10-6相差較大。這是由于解析法不考慮二次事件引發(fā)的進一步升級、忽略單元間相互作用,導致結(jié)果失真。

由表4和表7看出,蒙特卡洛模擬與其余兩種方法結(jié)果都比較接近;但T5模擬值4.035×10-6與其余兩種方法3.464×10-6相差較大,依據(jù)隨機數(shù)理論,內(nèi)在擴展方向不明確。蒙特卡洛模擬結(jié)果與隨機數(shù)有關,其模擬結(jié)果未正確反映位置差異的概率變化。

貝葉斯網(wǎng)絡對一級多米諾效應計算結(jié)果準確,其考慮二次事件引發(fā)的進一步升級、單元間相互作用,使二級多米諾效應計算結(jié)果更合理。貝葉斯網(wǎng)絡計算一級、二級多米諾效應概率為8.879×10-6、5.081×10-6,這是其獨特點。貝葉斯網(wǎng)絡中各節(jié)點參數(shù)和中間概率值能循環(huán)更新,進而提高計算結(jié)果準確性,這是相比解析法最大的優(yōu)勢。另外,貝葉斯網(wǎng)絡還能計算各級多米諾效應概率,通過概率值辨識事故序列。

總體上看,貝葉斯網(wǎng)絡能考慮多級多米諾效應、事故協(xié)同效應、根據(jù)單元間相互作用更新節(jié)點概率,清晰反映位置差異導致的節(jié)點概率變化,使事故傳播過程更完整,事故序列和概率值更合理,其操作簡潔、可靠性高,更適合油庫多米諾效應定量分析。

4.2 案例二

本案例考慮安全屏障下一級多米諾效應概率,兩種方法計算結(jié)果相差5%以內(nèi),計算可信度較高。

設置隔熱涂層后多米諾效應概率為5.414×10-7,相比于無安全屏障的5.859×10-6,降低一個數(shù)量級,防護效果較明顯。再設置噴淋水裝置后多米諾效應概率降低為3.654×10-7,比隔熱涂層防護效果提高30%。安全屏障越多,多米諾效應概率越低,但隨著安全屏障數(shù)目增加,計算會更復雜,表9為兩種方法設置安全屏障對比。

表9 解析法和貝葉斯網(wǎng)絡設置安全屏障對比Table 9 Comparison between analytic method and Bayesian network considering safety barriers

從表9看出,隨著安全屏障數(shù)目增加,解析法事件樹分支呈倍數(shù)型增長,規(guī)模逐漸變大,計算難度大幅增加;貝葉斯網(wǎng)絡增加防護節(jié)點,形成新網(wǎng)絡圖,可更新網(wǎng)絡中各節(jié)點概率,適用性好。因解析法和貝葉斯網(wǎng)絡對一級多米諾效應計算結(jié)果很接近,從方法操作和簡潔性上看,貝葉斯網(wǎng)絡更適用于安全屏障下多米諾效應計算。

5 結(jié) 論

(1)解析法對一級多米諾效應計算較為準確,但對二級以上多米諾效應計算誤差較大;蒙特卡洛模擬能適用于復雜系統(tǒng),但不能反映位置差異的概率變化;貝葉斯網(wǎng)絡考慮節(jié)點協(xié)同和相互作用,節(jié)點參數(shù)和概率值都能循環(huán)更新,事故傳播模式更全面,事故序列更完整,計算結(jié)果的準確度更高。

(2)隨著安全屏障數(shù)目增加,解析法事件樹規(guī)模增長快,計算難度大幅增加;貝葉斯網(wǎng)絡防護節(jié)點數(shù)少,計算適用性好,多米諾效應概率比無防護時降低一個數(shù)量級,其對多米諾效應預防有重要作用。

(3)貝葉斯網(wǎng)絡考慮事故傳播模式更全面、事故序列更完整,計算結(jié)果更符合實際,且操作方便、計算可靠性高,更加適用于油庫多米諾效應定量分析。