基于蒙特卡洛法的集裝箱堆場危險(xiǎn)品爆炸威力計(jì)算模型

陳武爭， 陳大鵬， 陳 力， 張婧卿， 方 秦

(1. 上海中交水運(yùn)設(shè)計(jì)研究有限公司, 上海 200092;2. 陸軍工程大學(xué) 爆炸沖擊防災(zāi)減災(zāi)國家重點(diǎn)實(shí)驗(yàn)室, 南京 210014;3. 上海國際港務(wù)(集團(tuán))股份有限公司，上海 200092)

危險(xiǎn)品集裝箱堆場是維持工業(yè)生產(chǎn)和城市運(yùn)行的重要基礎(chǔ)設(shè)施，由于大多數(shù)危險(xiǎn)品具有顯著的易燃、易爆特性，因此，危險(xiǎn)品堆場必須進(jìn)行爆炸安全規(guī)劃與防護(hù)設(shè)計(jì)。天津港危險(xiǎn)品倉庫“8.12”特大爆炸事故發(fā)生后，危險(xiǎn)品集裝箱堆場的選址、安全防護(hù)設(shè)計(jì)也越來越引起社會(huì)各界的重視，但目前我國尚未建立完善的針對(duì)危險(xiǎn)品集裝箱堆場的規(guī)范體系，特別是設(shè)計(jì)階段的標(biāo)準(zhǔn)、規(guī)范很不全面[1]，尤其缺少危險(xiǎn)品爆炸威力的計(jì)算方法和參數(shù)?，F(xiàn)有的爆炸相關(guān)研究多集中在軍事工程方面，民用領(lǐng)域的研究十分薄弱。

通過等效TNT當(dāng)量系數(shù)來描述爆炸物的爆炸威力是一類比較常見的方法。郝新紅等[2]提出了一種新的基于AutoReaGas的煙火藥爆炸數(shù)值模擬方法，得到了煙火藥在不同比例爆距下的TNT當(dāng)量系數(shù)。陽建紅等[3]通過10 kg TNT炸藥和NEPE(Nitrate Ester Plasticized Polyether Propellant)高能推進(jìn)劑的空中爆炸試驗(yàn)，計(jì)算出NEPE高能推進(jìn)劑在超壓為0.006～0.5 MPa范圍內(nèi)的TNT當(dāng)量系數(shù)為1.26。王肇中等[4]通過彈道拋擲法，量測了多種工業(yè)炸藥的做功能力及TNT當(dāng)量系數(shù)，但由于場地限制，很難廣泛應(yīng)用。王建靈等[5]通過熱刺激和機(jī)械刺激試驗(yàn)，獲得了疊氮硝胺在不同刺激下的響應(yīng)特征，確定在0.02～0.3 MPa峰值超壓范圍內(nèi)，疊氮硝胺TNT當(dāng)量系數(shù)為0.4～0.6。胡宏偉等[6]開展了多種非理想爆炸物在混凝土介質(zhì)內(nèi)的爆炸試驗(yàn)，通過測量爆炸在混凝土內(nèi)部產(chǎn)生的腔體體積來計(jì)算不同爆炸物的TNT當(dāng)量系數(shù)。胡廣霞等[7]通過對(duì)某危險(xiǎn)貨物集裝箱堆場進(jìn)行研究，通過理論分析與數(shù)值模擬，確定了堆場內(nèi)發(fā)生爆炸時(shí)的TNT當(dāng)量，但未考慮危險(xiǎn)貨物種類及堆存的隨機(jī)性。陳虹麗[8]以某煙花爆炸倉庫為例，基于等效TNT當(dāng)量分析，從沖擊波及廢氣污染物擴(kuò)散等方面對(duì)爆炸威力進(jìn)行了分析和預(yù)測?？梢钥闯觯F(xiàn)有研究大多針對(duì)某一特定危險(xiǎn)品或已知種類和數(shù)量的幾類危險(xiǎn)品的混合物。實(shí)際上，由于港口集裝箱堆場多為中轉(zhuǎn)堆場，堆存貨物雖有類別限制和場地區(qū)域劃分，但同一類別下不同種類貨物眾多。堆場內(nèi)危險(xiǎn)品種類和堆存量很難用單一或固定的指標(biāo)去描述，因此，準(zhǔn)確預(yù)測擁有巨大吞吐量和隨機(jī)危險(xiǎn)品種類的集裝箱堆場的爆炸威力困難重重。

另外一類常見的爆炸威力計(jì)算方法可以不區(qū)分炸藥種類，而是根據(jù)爆炸后果進(jìn)行分析評(píng)價(jià)，再與TNT爆炸后果進(jìn)行類比[9]。如地震波數(shù)據(jù)推算法[10]，爆炸成坑特征推算法[11]和建筑毀傷程度推算法[12-14]。此類方法針對(duì)爆炸事故進(jìn)行評(píng)估，基于事故形成的可追溯特征推算存儲(chǔ)貨物爆炸威力。但此時(shí)爆炸災(zāi)害已經(jīng)發(fā)生，破壞后果已形成，仍然難以適用于堆場的爆炸安全設(shè)計(jì)。

本文針對(duì)危險(xiǎn)品集裝箱堆場危險(xiǎn)品量大且隨機(jī)性強(qiáng)這一特點(diǎn)，結(jié)合我國東部沿海地區(qū)某一超大規(guī)模危險(xiǎn)品集裝箱堆場的設(shè)計(jì)建設(shè)需求，基于蒙特卡洛法，建立了一種概率相關(guān)的適用于集裝箱堆場的危險(xiǎn)品等效TNT當(dāng)量的計(jì)算模型，并針對(duì)1類危險(xiǎn)品討論了最大當(dāng)量系數(shù)和防護(hù)概率等關(guān)鍵參數(shù)的影響，為工程相關(guān)的安全防護(hù)設(shè)施設(shè)計(jì)、周邊重要設(shè)施保護(hù)提供了理論方法和計(jì)算途徑。

1 國內(nèi)外危險(xiǎn)品集裝箱堆場現(xiàn)狀

國內(nèi)外危險(xiǎn)品集裝箱堆場堆存方式主要有兩種：集中堆存和分散堆存(見圖1)。

國外危險(xiǎn)品集裝箱堆場大多采用分散堆存方式，便于分散風(fēng)險(xiǎn)，不同種類的集裝箱之間多采用低價(jià)值集裝箱、空箱、存儲(chǔ)物品易與危險(xiǎn)品產(chǎn)生中和反應(yīng)的集裝箱相隔離。此類堆場如發(fā)生事故，多為小規(guī)模點(diǎn)狀事故，不易引起連鎖反應(yīng)，事故等級(jí)和規(guī)模相對(duì)較小，產(chǎn)生的危害和造成的損失也較小。但同時(shí)也會(huì)帶來裝卸管理、監(jiān)管等比較困難的問題。

國內(nèi)危險(xiǎn)品集裝箱堆場主要以集中堆存為主[15]，在港口外建設(shè)專業(yè)的堆場，以滿足危險(xiǎn)品堆存、裝卸、監(jiān)管、安全等各種要求[16]。由于國內(nèi)是集中堆存，堆場一旦發(fā)生爆炸事故極易引發(fā)連鎖反應(yīng)，且港外集中堆存的安全防護(hù)設(shè)計(jì)比一般港內(nèi)危險(xiǎn)品堆場要求更高。根據(jù)《危險(xiǎn)貨物分類和品名編號(hào)：GB 6944—2012》，危險(xiǎn)品一般可劃分為1類～9類，1類危險(xiǎn)品箱存儲(chǔ)的危險(xiǎn)品均為極不穩(wěn)定的易燃易爆品，相對(duì)其他種類危險(xiǎn)品來說，其傷害更直接、迅速，傷害范圍更大[17]。一般來說，此類集裝箱宜單獨(dú)存放，并進(jìn)行抗爆防護(hù)設(shè)計(jì)，建設(shè)抗爆圍墻。因此，如何確定危險(xiǎn)品堆場的爆炸威力顯得至關(guān)重要。

2 基于概率的堆場爆炸威力計(jì)算模型

2.1 危險(xiǎn)品集裝箱堆場工況

我國東部某超大規(guī)模危險(xiǎn)品集裝箱堆場規(guī)劃堆存1類～6類和8類～9類危險(xiǎn)品，各類集裝箱在周邊港口卸船后由集裝箱卡車運(yùn)至該堆場進(jìn)行堆存、拆裝箱或拼箱作業(yè)，堆存周期相對(duì)較長。危險(xiǎn)品發(fā)生意外爆炸是該堆場可能發(fā)生的重大危險(xiǎn)事故之一。根據(jù)規(guī)劃，1類危險(xiǎn)品箱集中堆放至堆場東南角區(qū)域，需進(jìn)行抗爆防護(hù)設(shè)計(jì)計(jì)算并建設(shè)相應(yīng)的抗爆圍墻，如圖2所示。

圖2 危險(xiǎn)品集裝箱堆場示意圖Fig.2 Schematic diagram of the dangerous goods container yard

在圖2的危險(xiǎn)品集裝箱堆場中，1類危險(xiǎn)品區(qū)域設(shè)計(jì)平面箱位40 TEU，堆高兩層，最大容量80 TEU，常態(tài)容量60 TEU。1 TEU集裝箱容積為長609 cm×244 cm×259 cm，內(nèi)容積為569 cm×213 cm×218 cm，配貨毛重質(zhì)量一般為17.5 t。TEU(Twenty-feet Equivalent Unit)表示長度為609 cm的國際標(biāo)準(zhǔn)集裝箱單位。由于該堆場所處區(qū)域周邊港口擁有著巨大的集裝箱吞吐量，導(dǎo)致該堆場堆放的1類危險(xiǎn)品的種類及數(shù)量具有非常大的不確定性。

2.2 隨機(jī)模擬方法

隨機(jī)模擬法又稱蒙特卡洛(Monte Carlo)法，是一種通過設(shè)定隨機(jī)過程，反復(fù)生成時(shí)間序列，計(jì)算參數(shù)估計(jì)量和統(tǒng)計(jì)量，進(jìn)而研究其分布特征的方法。具體來說，當(dāng)系統(tǒng)中各個(gè)單元的可靠性特征量已知，但系統(tǒng)的可靠性過于復(fù)雜，難以建立可靠性預(yù)計(jì)的精確數(shù)學(xué)模型或模型太復(fù)雜而不便應(yīng)用時(shí)，可用隨機(jī)模擬法近似計(jì)算出系統(tǒng)可靠性的預(yù)計(jì)值。隨著模擬次數(shù)的增多，其預(yù)計(jì)精度也逐漸增高，在不斷對(duì)某一過程進(jìn)行模擬時(shí)，可通過預(yù)計(jì)值落在某一范圍內(nèi)的次數(shù)與模擬總次數(shù)的比值來近似確定預(yù)計(jì)值落在該范圍的概率。由于涉及到隨機(jī)過程的反復(fù)多次生成，蒙特卡洛模擬法是以高容量和高速度的計(jì)算機(jī)為前提條件的，因此只是在近些年才得到廣泛推廣。相比于傳統(tǒng)的數(shù)學(xué)統(tǒng)計(jì)方法，蒙特卡洛法可以通過借助計(jì)算機(jī)來減少人為計(jì)算的困難，且可操作性強(qiáng)，直觀易懂[18]。

2.3 理論計(jì)算模型

仍然選擇通過等效TNT當(dāng)量來表征危險(xiǎn)品的爆炸威力。通過對(duì)該危險(xiǎn)品集裝箱堆場周邊港口1類危險(xiǎn)品(爆炸品)進(jìn)行統(tǒng)計(jì)，將未來擬堆存在該項(xiàng)目堆場內(nèi)的爆炸品分為A～E五小類，分別對(duì)應(yīng)《危險(xiǎn)貨物分類和品名編號(hào)》中的1.1項(xiàng)～1.5項(xiàng)。由于危險(xiǎn)品爆炸產(chǎn)生的能量等于危險(xiǎn)品的質(zhì)量與其爆熱的乘積，可通過能量相似理論計(jì)算出各種危險(xiǎn)品的等效TNT當(dāng)量。假設(shè)每一小類危險(xiǎn)品的TNT當(dāng)量系數(shù)、密度和體積分別為γi,ρi,Vi(i取1～5)，則1類危險(xiǎn)品的TNT當(dāng)量可表示為

(1)

對(duì)于同一小類爆炸物，γi,ρi均可視為常數(shù)，因此可以引入一個(gè)新的轉(zhuǎn)換系數(shù)物理量αi=γiρi，量綱與ρi相同，表示單位體積的危險(xiǎn)品i，可以等效為質(zhì)量為αi的TNT炸藥。則此時(shí)式(1)可以轉(zhuǎn)化為

(2)

取αmax=max{αi}，αmin=min{αi}，則

Mmax=αmaxV,Mmin=αminV

(3)

(4)

蒙特卡洛法的誤差通常可以表示為[19]

(5)

對(duì)于一般工程，隨機(jī)數(shù)n通常取3 000～5 000[20]即可滿足工程精度要求。為了提高精度，分別取n=103，n=104，n=105，n=106，計(jì)算等效TNT當(dāng)量隨機(jī)分布值M的均值與標(biāo)準(zhǔn)差，每組模擬次數(shù)n進(jìn)行三次計(jì)算，得表1。

表1 不同模擬次數(shù)下的均值及標(biāo)準(zhǔn)差Tab.1 Mean value and standard deviation after different simulations

從表1可知，n=104,n=105，n=106時(shí)，標(biāo)準(zhǔn)差σ穩(wěn)定在0.275，故取σ=0.275。由正態(tài)分布表可知，顯著水平α=0.01時(shí)，λα=2.575 8，若控制誤差ε在0.1%量級(jí)，由式(5)計(jì)算可得n=5.01×105，故n>5.01×105時(shí)才能保證誤差控制在0.1%范圍內(nèi)，故取模擬次數(shù)n=1×106。

通過MATLAB完成n=106組隨機(jī)數(shù)的產(chǎn)生，利用式(4)計(jì)算出各組隨機(jī)數(shù)下的等效TNT當(dāng)量Mi(i=1,2,…,106)。等效TNT當(dāng)量隨機(jī)分布值計(jì)算框圖,如圖3所示。

圖3 計(jì)算流程框圖Fig.3 Flow chart of calculation process

由式(3)可知，M的取值范圍為[αminV,αmaxV]。任取[αminV,αmaxV]范圍內(nèi)的兩個(gè)數(shù)φ1,φ2，不妨設(shè)φ1<φ2。TNT當(dāng)量值M落在區(qū)間[φ1,φ2]的概率，等于M落在區(qū)間[φ1,φ2]的樣本個(gè)數(shù)nφ1,φ2與總樣本數(shù)n的比值，即

(6)

堆場危險(xiǎn)品等效TNT當(dāng)量隨機(jī)分布值M為離散型隨機(jī)變量，由于樣本容量n足夠大，可近似認(rèn)為M為連續(xù)型變量。令f(M)概率密度函數(shù)，簡稱概率密度。則

(7)

若不計(jì)高階無窮小，則：p(M1

(8)

由此可以根據(jù)式(8)構(gòu)造概率密度函數(shù)分布圖，根據(jù)其分布特征和概率需求，確定危險(xiǎn)品堆場的等效TNT當(dāng)量M。假設(shè)M′表示堆場等效TNT當(dāng)量的設(shè)計(jì)值。當(dāng)爆炸事故發(fā)生時(shí)，堆存危險(xiǎn)品的等效TNT當(dāng)量真實(shí)值不大于設(shè)計(jì)值M′時(shí)，認(rèn)為防護(hù)工程可以完成防護(hù)任務(wù)，即達(dá)到有效防護(hù)，TNT當(dāng)量真實(shí)值不大于設(shè)計(jì)值的概率可稱為有效防護(hù)概率pe。在根據(jù)TNT當(dāng)量做防護(hù)設(shè)計(jì)時(shí)，為求保險(xiǎn)，可以采用TNT當(dāng)量的最大值Mmax來進(jìn)行設(shè)計(jì)，即M′=Mmax，此時(shí)有效防護(hù)概率為100%。將有效防護(hù)概率適當(dāng)降低，在不影響防護(hù)效果的同時(shí)，等效TNT當(dāng)量設(shè)計(jì)值也將降低，從而達(dá)到節(jié)約防護(hù)成本的目的。

3 關(guān)鍵參數(shù)討論

3.1 有效防護(hù)概率pe

常見民用危險(xiǎn)品的等效TNT當(dāng)量系數(shù)范圍為0～1.5，密度為1～2 g/cm3[21]。由于αi=γiρi，故轉(zhuǎn)換系數(shù)α的范圍為0～3。為了更具一般性，取[0,3]中的五個(gè)隨機(jī)數(shù)(0.2,0.7,1.3,2.3,2.8)作為五種爆炸品的轉(zhuǎn)換系數(shù)，即(α1,α2,α3,α4,α5)=(0.2,0.7,1.3,2.3,2.8)。為方便計(jì)算，取堆場容量V=1，則Mmax=2.8，Mmin=0.2。

取區(qū)間(0,0.01],(0.01,0.02],…,(2.99,3]，通過式(5)分別求出各區(qū)間內(nèi)的概率，并由式(8)求出概率密度f(M)的分布,如圖4所示。由圖4可知，此時(shí)TNT當(dāng)量M近似于正態(tài)分布。

圖4 TNT當(dāng)量計(jì)算值分布Fig.4 TNT equivalent value distribution

求出樣本中M的均值，標(biāo)準(zhǔn)差分別為:μ=1.460,σ=0.354并畫出均值μ=1.460，標(biāo)準(zhǔn)差σ=0.354的正態(tài)分布圖，并與TNT當(dāng)量計(jì)算值分布圖進(jìn)行比較，結(jié)果見圖5，比較結(jié)果也證明了通常情況下，堆場危險(xiǎn)品TNT當(dāng)量隨機(jī)值應(yīng)該符合正態(tài)分布。

圖5 隨機(jī)模擬值與相應(yīng)正態(tài)分布比較Fig.5 Comparison between the random simulation value and the corresponding normal distribution

由于M服從于μ=1.460,σ=0.354的正態(tài)分布，故(M-μ)/σ服從于標(biāo)準(zhǔn)正態(tài)分布。查標(biāo)準(zhǔn)正態(tài)分布表，可得堆場的有效防護(hù)概率pe=[(M-μ)/σ<(M′-μ)/σ]=0.9時(shí)，(M′-μ)/σ=1.28，即堆場的危險(xiǎn)品等效TNT當(dāng)量設(shè)計(jì)值M′=1.91。定義λ為

(9)

可以看出，當(dāng)此工況的有效防護(hù)概率由100%降到90%時(shí)，等效TNT當(dāng)量設(shè)計(jì)值由Mmax=2.8降到M′=1.91，降低了31.8%，相應(yīng)的工程防護(hù)設(shè)計(jì)建設(shè)成本也會(huì)相應(yīng)下降。

3.2 最大轉(zhuǎn)換系數(shù)

經(jīng)過多組轉(zhuǎn)換系數(shù)的計(jì)算發(fā)現(xiàn)，當(dāng)有效防護(hù)概率由100%降到90%時(shí)，TNT當(dāng)量設(shè)計(jì)值下降的百分比是隨著不同工況發(fā)生變化的。由式(9)可知，在計(jì)算TNT當(dāng)量設(shè)計(jì)值下降比例λ時(shí)，Mmax對(duì)計(jì)算結(jié)果產(chǎn)生很大的影響。由于總體積V是定值，根據(jù)式(3)可知，最大轉(zhuǎn)換系數(shù)αmax是Mmax的關(guān)鍵參數(shù)。η表示最大轉(zhuǎn)化系數(shù)αmax與其余四個(gè)轉(zhuǎn)化系數(shù)之和的比值，即

(10)

式中：η的取值范圍為[0.25,∞)，η取到0.25的情況為五個(gè)轉(zhuǎn)換系數(shù)相等。例如：α1=α2=α3=α4=α5=1；η趨于無窮的情況為最大轉(zhuǎn)化系數(shù)為常數(shù)，而另外四個(gè)轉(zhuǎn)化系數(shù)之和趨于0時(shí)。為研究η的取值對(duì)下降比例λ的影響，分別取η=0.5,η=1,η=2,η=4,η=8進(jìn)行討論。當(dāng)α1=α2=α3=α4=0.1，α5=0.2時(shí)，η=0.5，由2.3中的方法求出概率密度f(M)的分布，得出此時(shí)的均值μ=0.12，標(biāo)準(zhǔn)差σ=0.011。等效TNT當(dāng)量計(jì)算值與μ=0.12,σ=0.011的正態(tài)分布的比較如圖6所示。

圖6 隨機(jī)模擬值與相應(yīng)正態(tài)分布比較(η=0.5)Fig.6 Comparison between the random simulation value and the corresponding normal distribution at η=0.5

表2 各工況下降比例λ計(jì)算表Tab.2 Calculation table of decreasing ratio λ for each working condition

由表2可知，η增大時(shí)，λ也隨之增大。說明最大轉(zhuǎn)化系數(shù)與其余四個(gè)轉(zhuǎn)化系數(shù)之間差距越小，有效防護(hù)概率降到90%時(shí)，TNT當(dāng)量設(shè)計(jì)值下降比例越?。蛔畲筠D(zhuǎn)化系數(shù)與其余四個(gè)轉(zhuǎn)化系數(shù)之間差距越大，TNT當(dāng)量設(shè)計(jì)值下降比例越大。

4 工程計(jì)算案例

繼續(xù)以東部某超大規(guī)模危險(xiǎn)品集裝箱堆場1類危險(xiǎn)品為計(jì)算背景，采用提出的爆炸威力計(jì)算方法計(jì)算該堆場的1類危險(xiǎn)品的等效TNT當(dāng)量。堆場內(nèi)標(biāo)準(zhǔn)集裝箱內(nèi)部危險(xiǎn)品的有效體積(除去貨架及外包裝的危險(xiǎn)品體積)為5 m3，擬堆存的危險(xiǎn)貨物見表3。

表3 標(biāo)準(zhǔn)箱內(nèi)擬堆存的危險(xiǎn)貨物[22]Tab.3 Dangerous goods to be stored in standard case

由表3可知，αmax=1.38 g/cm3，αmin=0.42 g/cm3,則Mmax=6.9 t，Mmin=2.1 t。為使誤差控制在0.1%范圍內(nèi)，基于“2.3”節(jié)中模擬次數(shù)的確定方法，取模擬次數(shù)n=106。通過MATLAB實(shí)現(xiàn)計(jì)算出106組TNT當(dāng)量Mi(i=1,2,…，106)的隨機(jī)分布值，取區(qū)間(2.1,2.11],(2.11,2.12],…，(6.89,6.9],基于“3.1”節(jié)中的計(jì)算結(jié)論，認(rèn)為此時(shí)TNT當(dāng)量隨機(jī)分布值Mi服從于μ=4.630,σ=0.582的正態(tài)分布，故(M-μ)/σ服從于標(biāo)準(zhǔn)正態(tài)分布。查標(biāo)準(zhǔn)正態(tài)分布表，可得有效防護(hù)概率p[(M-μ)/σ<(M′-μ)/σ]=0.9時(shí)，(M′-μ)/σ=1.28。既而，得出該集裝箱的危險(xiǎn)品等效TNT當(dāng)量設(shè)計(jì)值為M′=5.37 t。

5 結(jié) 論

本文針對(duì)危險(xiǎn)品集裝箱堆場的等效TNT當(dāng)量開展研究，考慮了堆存危險(xiǎn)品種類及數(shù)量的不確定性，建立了一種概率相關(guān)的爆炸威力計(jì)算模型，并結(jié)合工程建設(shè)算例，討論了防護(hù)概率和最大轉(zhuǎn)化系數(shù)對(duì)TNT當(dāng)量設(shè)計(jì)值的影響規(guī)律，主要結(jié)論有：

(1)堆場內(nèi)隨機(jī)堆放1類危險(xiǎn)品時(shí)，等效TNT當(dāng)量值基本符合正態(tài)分布。

(2)將1類危險(xiǎn)品的有效防護(hù)概率由100%降到90%時(shí)，等效TNT當(dāng)量設(shè)計(jì)值可下降30%～70%，在不影響有效防護(hù)的同時(shí)，可大大降低設(shè)計(jì)和建造成本。

(3)不同種類危險(xiǎn)品最大轉(zhuǎn)化系數(shù)與其余轉(zhuǎn)化系數(shù)之間的差距越大，等效TNT當(dāng)量設(shè)計(jì)值下降越明顯。