基于蒙特卡洛法的集裝箱堆場危險品爆炸威力計算模型

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

(1. 上海中交水運設計研究有限公司, 上海 200092;2. 陸軍工程大學 爆炸沖擊防災減災國家重點實驗室, 南京 210014;3. 上海國際港務(集團)股份有限公司，上海 200092)

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

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

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

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

1 國內外危險品集裝箱堆場現狀

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

國外危險品集裝箱堆場大多采用分散堆存方式，便于分散風險，不同種類的集裝箱之間多采用低價值集裝箱、空箱、存儲物品易與危險品產生中和反應的集裝箱相隔離。此類堆場如發生事故，多為小規模點狀事故，不易引起連鎖反應，事故等級和規模相對較小，產生的危害和造成的損失也較小。但同時也會帶來裝卸管理、監管等比較困難的問題。

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

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

2.1 危險品集裝箱堆場工況

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

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

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

2.2 隨機模擬方法

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

2.3 理論計算模型

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

(1)

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

(2)

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

Mmax=αmaxV,Mmin=αminV

(3)

(4)

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

(5)

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

表1 不同模擬次數下的均值及標準差Tab.1 Mean value and standard deviation after different simulations

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

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

圖3 計算流程框圖Fig.3 Flow chart of calculation process

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

(6)

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

(7)

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

(8)

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

3 關鍵參數討論

3.1 有效防護概率pe

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

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

圖4 TNT當量計算值分布Fig.4 TNT equivalent value distribution

求出樣本中M的均值，標準差分別為:μ=1.460,σ=0.354并畫出均值μ=1.460，標準差σ=0.354的正態分布圖，并與TNT當量計算值分布圖進行比較，結果見圖5，比較結果也證明了通常情況下，堆場危險品TNT當量隨機值應該符合正態分布。

圖5 隨機模擬值與相應正態分布比較Fig.5 Comparison between the random simulation value and the corresponding normal distribution

由于M服從于μ=1.460,σ=0.354的正態分布，故(M-μ)/σ服從于標準正態分布。查標準正態分布表，可得堆場的有效防護概率pe=[(M-μ)/σ<(M′-μ)/σ]=0.9時，(M′-μ)/σ=1.28，即堆場的危險品等效TNT當量設計值M′=1.91。定義λ為

(9)

可以看出，當此工況的有效防護概率由100%降到90%時，等效TNT當量設計值由Mmax=2.8降到M′=1.91，降低了31.8%，相應的工程防護設計建設成本也會相應下降。

3.2 最大轉換系數

經過多組轉換系數的計算發現，當有效防護概率由100%降到90%時，TNT當量設計值下降的百分比是隨著不同工況發生變化的。由式(9)可知，在計算TNT當量設計值下降比例λ時，Mmax對計算結果產生很大的影響。由于總體積V是定值，根據式(3)可知，最大轉換系數αmax是Mmax的關鍵參數。η表示最大轉化系數αmax與其余四個轉化系數之和的比值，即

(10)

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

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

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

由表2可知，η增大時，λ也隨之增大。說明最大轉化系數與其余四個轉化系數之間差距越小，有效防護概率降到90%時，TNT當量設計值下降比例越小；最大轉化系數與其余四個轉化系數之間差距越大，TNT當量設計值下降比例越大。

4 工程計算案例

繼續以東部某超大規模危險品集裝箱堆場1類危險品為計算背景，采用提出的爆炸威力計算方法計算該堆場的1類危險品的等效TNT當量。堆場內標準集裝箱內部危險品的有效體積(除去貨架及外包裝的危險品體積)為5 m3，擬堆存的危險貨物見表3。

表3 標準箱內擬堆存的危險貨物[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%范圍內，基于“2.3”節中模擬次數的確定方法，取模擬次數n=106。通過MATLAB實現計算出106組TNT當量Mi(i=1,2,…，106)的隨機分布值，取區間(2.1,2.11],(2.11,2.12],…，(6.89,6.9],基于“3.1”節中的計算結論，認為此時TNT當量隨機分布值Mi服從于μ=4.630,σ=0.582的正態分布，故(M-μ)/σ服從于標準正態分布。查標準正態分布表，可得有效防護概率p[(M-μ)/σ<(M′-μ)/σ]=0.9時，(M′-μ)/σ=1.28。既而，得出該集裝箱的危險品等效TNT當量設計值為M′=5.37 t。

5 結 論

本文針對危險品集裝箱堆場的等效TNT當量開展研究，考慮了堆存危險品種類及數量的不確定性，建立了一種概率相關的爆炸威力計算模型，并結合工程建設算例，討論了防護概率和最大轉化系數對TNT當量設計值的影響規律，主要結論有：

(1)堆場內隨機堆放1類危險品時，等效TNT當量值基本符合正態分布。

(2)將1類危險品的有效防護概率由100%降到90%時，等效TNT當量設計值可下降30%～70%，在不影響有效防護的同時，可大大降低設計和建造成本。

(3)不同種類危險品最大轉化系數與其余轉化系數之間的差距越大，等效TNT當量設計值下降越明顯。