齊心歌, 王海清, 田英帥, 陳國(guó)明

(1.中國(guó)石油大學(xué)(華東)機(jī)電工程學(xué)院,山東青島 266580; 2.深圳燃?xì)饧瘓F(tuán)安全管理部,廣東深圳 518040)

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

1 不確定度理論

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

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

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

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

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

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

2.1 氣體泄漏

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

(1)

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

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

2.2 氣體擴(kuò)散的等價(jià)氣云

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

高斯模型平均濃度方程為

(2)

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

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

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

(3)

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

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

SESC=πab,

(4)

其中

式中,a、b分別為濃度等值線等價(jià)橢圓的長(zhǎng)、短半軸;SESC為高度z一定時(shí)等價(jià)橢圓的面積。

(5)

2.3 多能法計(jì)算爆炸后果

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

E=3.5×103VESC.

(6)

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

(7)

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

ps=p0p′.

(8)

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

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

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

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

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

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

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

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

(9)

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

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

(10)

(11)

4 安全系數(shù)

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

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

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

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

(12)

其中

(13)

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

(14)

總方差計(jì)算為

(15)

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

(16)

4.2 安全系數(shù)計(jì)算

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

(17)

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

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

sc=1-uc.

(18)

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

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

5 案例分析

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

5.1 輸入樣本選擇

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

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

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

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

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

爆炸超壓與輸入?yún)?shù)的散點(diǎn)圖(圖3)可直觀地展示輸入?yún)?shù)的不確定性對(duì)輸出結(jié)果的影響,各個(gè)輸入?yún)?shù)的歸一化擬合斜率比值近似為泄漏孔徑:流速:流量系數(shù):風(fēng)速=28.0∶1.0∶4.8∶21.6,可以看出輸入?yún)?shù)對(duì)結(jié)果影響依次為泄漏孔徑>風(fēng)速>流量系數(shù)>流速,其中泄漏孔徑與風(fēng)速對(duì)結(jié)果影響較大且呈現(xiàn)不明顯的線性關(guān)系,其次為介質(zhì)流速,流量系數(shù)的變化對(duì)輸出量的影響很小。為對(duì)非線性關(guān)系進(jìn)行量化,此處利用Sobol指數(shù)法對(duì)等價(jià)氣云爆炸模型作參數(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)超壓值散點(diǎn)圖Fig.3 Overpressure scatter plot in sample interval

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

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

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

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

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

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

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

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

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

6 結(jié) 論

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

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