齊心歌, 王海清, 田英帥, 陳國明

(1.中國石油大學(xué)(華東)機電工程學(xué)院,山東青島 266580; 2.深圳燃氣集團安全管理部,廣東深圳 518040)

控制室是石化行業(yè)作業(yè)區(qū)內(nèi)的指揮中樞,是保證工作人員安全以及生產(chǎn)活動正常進行的關(guān)鍵機構(gòu),對于維持生產(chǎn)設(shè)備的運轉(zhuǎn)、作業(yè)的進行以及人員安全具有至關(guān)重要的作用[1]。因此對控制室的防爆進行研究,采用氣體探測GDS系統(tǒng)(氣體探測系統(tǒng),gas detection system),將氣體泄漏擴散造成的風(fēng)險控制在控制室的爆炸載荷范圍內(nèi),對于保障作業(yè)區(qū)的人員以及設(shè)備的安全具有重要意義。目前對控制室的防爆主要從兩個方面進行研究:一是從設(shè)計載荷的角度,以爆炸產(chǎn)生的沖擊波超壓值來表征控制室能夠承受的爆炸載荷,如國內(nèi)外標準規(guī)范ASCE 41088[2]、GB50779-2012[3]等;二是從概率爆炸載荷的角度,綜合考慮氣體的泄漏場景、氣云擴散、點火概率以及點火位置等得到作用于控制室的爆炸載荷[4-5]。以上兩者可得到控制室各個方向的設(shè)計載荷以及發(fā)生爆炸事故時作用于控制室沖擊波超壓[6-7],但是在GDS系統(tǒng)的應(yīng)用中,針對控制室載荷的探測閾值尚未給出明確的界定,探測載荷閾值須滿足一定的安全裕度,降低探測不確定度的影響,并為應(yīng)急措施的實施提供充分的預(yù)留時間。筆者基于控制室的概率爆炸載荷,充分考慮各類影響因素的不確定度,得到探測中控制室爆炸載荷的安全系數(shù),增強GDS系統(tǒng)設(shè)計輸入量的保守性。

1 不確定度理論

在實際爆炸載荷設(shè)計中,由于氣體泄漏、擴散以及氣云爆炸過程的高度復(fù)雜性,在計算過程中會涉及大量不確定性參數(shù),是GDS系統(tǒng)探測有效性的關(guān)鍵。從不確定度著手,得到控制室的爆炸載荷,作為GDS系統(tǒng)設(shè)計輸入量的基礎(chǔ)。

目前被廣泛認可的不確定度分類由Apostolakis首次提出,分為兩類:隨機不確定性和知識不確定性。隨機不確定性是系統(tǒng)本身存在的,是固有的、不可降低的,也稱為A型不確定性;知識不確定性可稱為B型不確定性,是指由于知識的欠缺導(dǎo)致的潛在不準確性,隨著知識的完備可逐漸降低[8]。因此對不確定性的分析主要圍繞知識不確定性展開。

在工程應(yīng)用中,對知識不確定性的分析需要從計算模型的參數(shù)入手,將輸入?yún)?shù)的不確定性量化,并判斷其對輸出結(jié)果不確定性的影響。對不確定性分析應(yīng)用最廣泛的是經(jīng)典概率方法,且有大量基于概率的不確定性方法被提出,包括蒙特卡洛法(MC法)、微分法、響應(yīng)面法等,其中MC法由于原理易于實現(xiàn),在多個領(lǐng)域得到廣泛應(yīng)用。

由于輸入?yún)?shù)在不同模型對結(jié)果的影響不同,需要通過參數(shù)敏感性分析,判斷其對結(jié)果不確定性的影響程度。參數(shù)敏感性分析方法包括全局敏感性分析和局部敏感性分析兩類。其中局部敏感性分析方法易于實現(xiàn),但每次只能分析單個參數(shù)不確定性的變化對結(jié)果的影響且對線性關(guān)系模型較為有效;全局敏感性分析方法可定量分析非線性模型的不同輸入?yún)?shù)同時變化時產(chǎn)生的交互作用對輸出結(jié)果不確定性的影響,本文中采用Sobol指數(shù)法對全局敏感性進行分析。

2 爆炸結(jié)果模型及參數(shù)分析

對于氣體泄漏引起的爆炸事故后果需要從泄漏源、擴散過程、氣云的爆炸一系列過程進行分析[9-10]。對氣體泄漏源分析主要針對儲罐等裝置的小孔泄漏以及底部管道的斷裂泄漏,因此采用伯努利方程;常用的氣體擴散模型有CFD模型、高斯擴散模型、唯象模型、相似模型等,綜合考慮計算效率與精確度,選擇高斯擴散模型;氣云爆炸的典型模型是多能法,經(jīng)過了大量實驗的驗證以及修正,充分考慮了氣體活性、局部約束、湍流加速等因素,與實際爆炸場景較為接近[11]。

2.1 氣體泄漏

利用伯努利方程計算氣體泄漏速率,表示為

(1)

式中,QL為氣體的泄漏速率,kg/s;Cd為孔流系數(shù),若泄漏孔為圓形,一般取1;A為孔徑的面積,m2;ρ為流體密度,kg/m3;p為工作壓力,kPa;p0為大氣壓力,kPa;u1為裝置中氣體流速,m/s。

根據(jù)式(1),氣體泄漏過程中對泄漏速率造成主要影響的因素包括泄漏孔的直徑和形狀、模型采用的流量系數(shù)以及運行壓力等。

2.2 氣體擴散的等價氣云

利用高斯模型計算過程中,氣候條件(如風(fēng)速、風(fēng)向、濕度)、表面粗糙度等因素對計算結(jié)果產(chǎn)生影響[12]。

高斯模型平均濃度方程為

(2)

式中,C(x,y,z)為氣體泄漏時,給定地點(x,y,z)的密度,kg·m-3;u為大氣風(fēng)速,m·s-1;σy,σz為側(cè)風(fēng)向和垂直風(fēng)向的擴散系數(shù),與大氣穩(wěn)定度以及下風(fēng)向距離x有關(guān),m;x,y,z為下風(fēng)向,側(cè)風(fēng)向和垂直風(fēng)向的距離,m;Hr為泄漏點源高度,m。

基于高斯擴散模型提出等價氣云的概念,用均相的立方體氣云來代替密度、濃度不均勻、形狀不規(guī)則的非均相氣云,以此來減少需要模擬的泄漏場景數(shù)量[13-14]。

計算濃度等值線,需要從x、y、z三個方向得到氣體擴散距離,以x為變量表示y,將高斯模型式(2)變形,表示為

(3)

式中,σy、σz均為x的函數(shù),取值依據(jù)Pasquill-Turner大氣穩(wěn)定度公式。

式(3)中,在場景范圍內(nèi)將泄漏源作為坐標原點,下風(fēng)向為x軸,側(cè)風(fēng)向為y軸,垂直風(fēng)向為z軸建立坐標系。沿x軸方向取一定步長(如1 m)并生成矩陣,在一定濃度閾值下,分別計算對應(yīng)的y值。計算過程中,若x取值過大或過小,可能會出現(xiàn)y無實數(shù)解的情況。物理意義解釋為:氣體擴散過程中形成的濃度等值線可能是封閉的,并與泄漏源有一定距離,因此在一定x范圍內(nèi)才會有對應(yīng)的y值。

SESC=πab,

(4)

其中

式中,a、b分別為濃度等值線等價橢圓的長、短半軸;SESC為高度z一定時等價橢圓的面積。

(5)

2.3 多能法計算爆炸后果

氣云發(fā)生擴散后,若云團中氣體濃度在爆炸極限內(nèi),遇點火源會引發(fā)爆炸事故[15]。對氣云爆炸后果的分析,擬采用多能法[16],得到對建筑物破壞的影響距離。具體實現(xiàn)步驟如圖1所示。

E=3.5×103VESC.

(6)

其中3.5×103kJ/m3為同化學(xué)計量濃度下烴-空氣混合物的典型燃燒熱值。

(7)

式中,R′為無量綱距離;z為爆炸中心與待分析設(shè)備的距離,m;p0為大氣壓力,通常取101.325 kPa。

ps=p0p′.

(8)

式中，ps為實際峰值側(cè)向超壓，kPa;p′為無量綱比擬最大側(cè)向超壓。

由于爆炸事故對設(shè)備以及建筑等造成損害的主要形式為沖擊波超壓,因此以沖擊波超壓值作為輸出參數(shù)。

圖1 多能法計算爆炸后果流程Fig.1 Flow chart of explosion consequence calculation by multi-energy method

3 基于MC模型的參數(shù)不確定度

基于MC模型計算不確定度步驟如下:

(1)確定不確定度來源。不確定度分析的主要目的是依據(jù)模型的不確定性與輸入的不確定性修正預(yù)測結(jié)果。對計算模型中涉及的參數(shù)進行分析,確定不確定性參數(shù)。

(2)確定不確定性參數(shù)的概率特征并生成隨機輸入樣本。依據(jù)以往的擴散場景、歷史數(shù)據(jù)等擬合得到輸入?yún)?shù)的概率密度函數(shù),基于概率密度函數(shù)用抽樣法得到離散分布的樣本[17-18]。樣本的大小以及抽樣方法的選擇決定了MC模型模擬的精度。

樣本數(shù)量可通過統(tǒng)計容許極限(a%,b%)確定。統(tǒng)計容許極限是指在b%的置信區(qū)間下,保證抽取的樣本在其分布中的比例不小于a%。樣本數(shù)量M可確定為

(9)

目前常用的樣本抽樣方法有分層抽樣法、簡單隨機抽樣法、拉丁超立方抽樣法(Latin hypercube sampling,LHS)以及漢莫斯里序列抽樣法等。其中漢莫斯里序列抽樣法與LHS樣本分布較均勻,且LHS可以在較少樣本的情況下達到更高的精度,因此選擇LHS[19-20]。假設(shè)需要在m維空間抽取n個樣本:①在參數(shù)取值范圍內(nèi),將每一維度分成n個互不重疊區(qū)間(通常使用均勻分布,使得區(qū)間長度相同);②在每一維度的每一個區(qū)間中隨機抽取一個點;③從每一維度隨機抽取步驟(2)中的點,組成矩陣,作為輸入樣本。

(3)以生成的樣本為輸入計算對應(yīng)的輸出,并利用期望與標準差的值對不確定性分析結(jié)果進行表征,得到輸出結(jié)果的不確定度,表示為

(10)

(11)

4 安全系數(shù)

4.1 不確定性參數(shù)的全局參數(shù)敏感性

在敏感性分析中,散點圖可定性獲得輸入?yún)?shù)的不確定性對結(jié)果的影響。散點圖是將輸入量作為橫坐標,輸出量作為縱坐標,利用抽樣法獲取大量樣本并得到對應(yīng)輸出結(jié)果,在坐標系中描繪,可定性觀察輸入量對輸出量的影響。雖然具有形象直觀的特點,但不能定量得到兩者的關(guān)系,因此在散點圖的基礎(chǔ)上,需采用其他全局敏感性分析法對輸入?yún)?shù)與輸出結(jié)果之間的關(guān)系進行定量分析。

目前,最常用的全局敏感度分析方法為Sobol指數(shù)法[21-22],基本思想是利用總方差表示全部變量對輸出結(jié)果的影響,利用偏方差表示單變量或多變量對輸出結(jié)果的影響,不同輸入?yún)?shù)對輸出結(jié)果的影響程度依據(jù)該參數(shù)在總方差中所占的比例來衡量。分析步驟如下:

(1)假定某一函數(shù)Y,可表示為自變量為X的一系列函數(shù)之和,即

(12)

其中

(13)

(2)參數(shù)xi的方差可表示為

(14)

總方差計算為

(15)

(3)輸入?yún)?shù)xi對輸出y的影響可通過一階敏感性指數(shù)確定,表示為

(16)

4.2 安全系數(shù)計算

得到輸入?yún)?shù)的不確定性產(chǎn)生的綜合不確定性影響為

(17)

式中,p為不確定參數(shù)數(shù)量。

進而得到安全系數(shù)為

sc=1-uc.

(18)

綜上,對控制室爆炸載荷安全系數(shù)進行分析,需綜合考慮氣體泄漏、擴散以及形成氣云的爆炸過程。通過分析不確定度來源得到不確定性參數(shù),通過全局敏感性分析得到不同參數(shù)對輸出結(jié)果的影響程度,整體分析步驟如圖2所示。

圖2 控制室載荷安全系數(shù)計算流程Fig.2 Calculation flow of safety coefficient

5 案例分析

以某LNG罐區(qū)為例,假設(shè)罐區(qū)為矩形,平面尺寸為260 m×240 m;有兩座并排LNG液態(tài)儲罐,尺寸為:底面直徑80 m,高度22 m;儲罐中心與控制室的水平距離為50 m,當(dāng)?shù)貧庀髼l件為全年風(fēng)速分布在1～5 m/s。以其中某一儲罐底部管線破裂為例,分析氣體泄漏擴散導(dǎo)致的爆炸事故形成的沖擊波超壓過程,進而分析其對控制室設(shè)計的影響并得到其安全系數(shù)。

5.1 輸入樣本選擇

利用蒙特卡洛方法對不確定度進行分析,需要利用拉丁超立方抽樣法確定輸入樣本。

(1)輸入樣本大小確定。依據(jù)式(9),選取置信區(qū)間為97%,即a=b=97,解不等式得到N>116。樣本值可取為大于116的任意實數(shù)值,但樣本過大,導(dǎo)致計算效率降低,因此案例中取樣本容量為120。

(2)不確定性參數(shù)取值區(qū)間確定。選取不確定性參數(shù)為泄漏孔徑、介質(zhì)流速、流量系數(shù)、環(huán)境風(fēng)速,概率密度函數(shù)均為均勻分布。根據(jù)實際工況,選取不確定性參數(shù)的取值區(qū)間如下:泄漏孔徑的取值區(qū)間擬選擇中孔泄漏,孔徑范圍為0.01～0.05 m;介質(zhì)流速區(qū)間依據(jù)介質(zhì)在儲罐與管線以及其他設(shè)備中的不同流速來確定,取0～4 m/s;流量系數(shù)的取值與泄漏孔形狀有關(guān),當(dāng)形狀為圓型、三角形、矩形或其他形狀時,流量系數(shù)取為0.9～1.0;環(huán)境風(fēng)速取值區(qū)間的確定依據(jù)當(dāng)?shù)氐臍庀笄闆r統(tǒng)計數(shù)據(jù),為1～5 m/s。在此區(qū)間內(nèi)生成輸入?yún)?shù)的樣本矩陣。

5.2 不確定度與參數(shù)敏感性

氣體泄漏孔徑泄漏后發(fā)生擴散,基于等價氣云理論,將樣本矩陣帶入式(1)～(5)得到擴散形成的等價氣云體積,其中式(2)～(3)中濃度閾值的確定依據(jù)LNG的氣體爆炸下限;選擇多能法最為氣體爆炸事故后果分析的方法,依據(jù)式(6)～(8)得到氣云爆炸導(dǎo)致的沖擊波超壓。利用蒙特卡洛方法,依據(jù)式(10)～(11)計算得到輸出值的均值與標準差,標準差即為不確定度。

爆炸超壓與輸入?yún)?shù)的散點圖(圖3)可直觀地展示輸入?yún)?shù)的不確定性對輸出結(jié)果的影響,各個輸入?yún)?shù)的歸一化擬合斜率比值近似為泄漏孔徑:流速:流量系數(shù):風(fēng)速=28.0∶1.0∶4.8∶21.6,可以看出輸入?yún)?shù)對結(jié)果影響依次為泄漏孔徑>風(fēng)速>流量系數(shù)>流速,其中泄漏孔徑與風(fēng)速對結(jié)果影響較大且呈現(xiàn)不明顯的線性關(guān)系,其次為介質(zhì)流速,流量系數(shù)的變化對輸出量的影響很小。為對非線性關(guān)系進行量化,此處利用Sobol指數(shù)法對等價氣云爆炸模型作參數(shù)敏感性分析。式(12)～(16)描述了Sobol指數(shù)法,依據(jù)式(12)～(14),建立N×2D(其中N為樣本數(shù),D為輸入變量數(shù)目)的樣本矩陣,樣本的建立同樣利用拉丁超立方抽樣法。將矩陣的前D列設(shè)置為矩陣E,后D列設(shè)置為矩陣F,構(gòu)造N×D的矩陣EFi(i=1,2,…,D),即用矩陣E中的第i列替換矩陣A的第i列,因此得到輸入矩陣E,F,EF1,EF2,EF3,EF4,由此就有(D+2)×N(即720)組輸入數(shù)據(jù),因此將有720組超壓輸出值,如圖4所示。依據(jù)式(15)～(16)可得到參數(shù)的敏感性分析結(jié)果,如表1所示。

圖3 樣本區(qū)間內(nèi)超壓值散點圖Fig.3 Overpressure scatter plot in sample interval

圖4 參數(shù)敏感度過程矩陣超壓柱形圖Fig.4 Overpressure histogram of process matrix to calculate parameter sensitivity

表1 不確定度與敏感性分析結(jié)果

在此場景下,依據(jù)式(17)得到安全系數(shù)為sc1=1-19.8%=80.2%,從敏感性指數(shù)可以看出,泄漏孔徑與環(huán)境風(fēng)速的變化對輸出結(jié)果影響較大。

作為對比,假設(shè)儲罐中心與控制室的水平距離改變?yōu)?0 m,其他參數(shù)不變,得到該場景下的不確定性與參數(shù)敏感性分析結(jié)果,如表2所示。

表2 不確定度與敏感性分析結(jié)果(距離改變)Table 2 Results of uncertainty and sensitivity (distance changing)

該場景下的安全系數(shù)為sc2=82.8%,與參數(shù)改變前得到的安全系數(shù)基本一致,不確定度與參數(shù)敏感性的分布變化不大。

為驗證安全系數(shù)的準確度,假設(shè)設(shè)備工作壓力發(fā)生改變,得到輸入?yún)?shù)的不確定度與Sobol指數(shù),如表3所示。該場景下的安全系數(shù)為sc3=82.86%。

表3 不確定度與敏感性分析結(jié)果(工作壓力改變)Table 3 Results of uncertainty and sensitivity (working pressure changing)

由此可見,安全系數(shù)隨場景的不同而發(fā)生一定改變,但數(shù)值一般約為80%,因此對建筑物能夠承受的沖擊波超壓閾值進行分析時,需考慮輸入?yún)?shù)與模型不確定度的影響,并依據(jù)泄漏場景確定對應(yīng)的安全系數(shù)。且安全系數(shù)的應(yīng)用對于增加計算結(jié)果的保守性以及事故防控具有重要意義。

6 結(jié) 論

(1)以等價氣云爆炸事故為基礎(chǔ),綜合考慮了氣體泄漏擴散過程與氣云爆炸過程中的不確定性因素,實現(xiàn)了以控制室為代表的建筑物爆炸載荷安全系數(shù)的確定。

(2)應(yīng)用于某LNG罐區(qū),得到不同場景下的安全系數(shù),增加了控制室等建筑物沖擊波超壓爆炸載荷的保守性,為GDS系統(tǒng)探測設(shè)計輸入提供理論支持。