柱形罐爆炸碎片拋射的MonteCarlo分析*-

劉振翼，黃 平，徐亞博

（北京理工大學爆炸科學與技術國家重點實驗室，北京 100081）

柱形罐爆炸碎片拋射的MonteCarlo分析*-

劉振翼，黃 平，徐亞博

（北京理工大學爆炸科學與技術國家重點實驗室，北京 100081）

在動力學分析的基礎上，將二維體系內(nèi)的碎片軌跡方程擴展到三維體系內(nèi)，將風這一影響因素考慮在內(nèi)。以墨西哥事故中水平圓柱儲罐發(fā)生沸騰液體擴展為蒸氣云爆炸（BLEVE）為例，利用Monte-Carlo法模擬出碎片的軌跡曲線，計算得到有風和無風情況下碎片拋射距離的概率分布，得到一般情況下風對碎片危害的影響較小。計算了無風情況下碎片碰撞目標容器的概率，得到碰撞概率隨罐間距的增大而呈現(xiàn)嚴格的指數(shù)衰減趨勢，通過此關系式可以計算出符合安全要求的罐間距。研究結果對于提高儲罐的安全性、緩解和控制碎片產(chǎn)生的風險具有一定的指導意義。

爆炸力學；沸騰液體擴展為蒸氣云爆炸；Monte-Carlo法；碎片；柱形罐

1 引 言

危險設施由于裝載過量、腐蝕及火災等原因產(chǎn)生的爆炸除了產(chǎn)生沖擊波，一部分能量還會轉(zhuǎn)化為碎片的動能。這種具有動能的碎片在飛行過程中有可能造成人員傷亡或目標設備的損壞，進而造成更大的災難事故，這種現(xiàn)象被稱為多米諾效應或連鎖效應。由于多米諾事故的嚴重性，歐盟的塞維索法令II對重大危險源的多米諾效應做出專門規(guī)定，要求重大危險源企業(yè)提交的安全報告中要有多米諾效應分析的內(nèi)容，可能產(chǎn)生多米諾效應的企業(yè)之間要共享信息等［1－5］?；馂臒彷椛洹⒈_擊波和碎片是造成多米諾效應的主要原因，本文中將研究碎片造成多米諾事故這部分內(nèi)容。荷蘭的應用技術研究院和美國的化學工藝安全中心給出了碎片拋射距離的經(jīng)驗計算公式。U.Hauptmanns［6－7］通過對事故碎片的隨機和不確定性進行分析，利用Monte-Carlo方法得到爆炸碎片到達不同距離的概率曲線。G.Cozzani等［8］在U.Hauptmanns的基礎上，考慮目標的影響，得到爆炸碎片碰撞目標容器的概率。文獻［9-11］采用類似文獻［7-8］的方法來計算碎片的拋射范圍。

本文中，將在前人研究的基礎上，考慮風向和風速，將目前通用的二維體系內(nèi)的碎片軌跡方程擴展到三維體系內(nèi)，利用Monte-Carlo法得到柱形罐爆炸碎片在有風和無風情況下到達不同距離的概率曲線，及撞擊不同距離處的目標容器的概率，并對2種情況下的計算結果進行比較分析。

2 三維體系內(nèi)碎片軌跡方程

將事故源、碎片和目標放入三維體系中分析，如圖1所示。由于碎片與目標相比通常較小，為簡化計算，忽略不計碎片的大小。由于碎片主要繞著質(zhì)心運動，因此可用質(zhì)心的軌跡和速度來代表整個碎片的軌跡和速度。參考系O-xyz的原點O為碎片質(zhì)心的初始位置（即事故源），其中xOy與地面平行，z軸與重力方向相反。向量vp代表碎片的初速方向，與xOy平面的夾角為φ，在xOy平面上的投影與x軸的夾角為θ。向量vw代表風向，與x軸的夾角為φ，在此僅考慮與xOy面平行的風。

碎片在飛行過程中受重力、阻力和升力的聯(lián)合作用。由于引入升力會造成計算過于復雜，且升力只適用于端蓋，當升力方向與端蓋軌跡的夾角超過10°時，升力可忽略不計，因此在實際計算中，為了簡化，升力不計在內(nèi)［7］。根據(jù)牛頓第二定律得到

式中：D是阻力，g是重力加速度，ap是碎片的加速度，mp是碎片的質(zhì)量。

對于空氣阻力的作用，Lees分析了3種情況：沒有阻力；阻力與碎片速度成正比；阻力與碎片速度的平方成正比。認為在亞音速范圍內(nèi)，阻力與碎片速度的平方成正比這種情況是最合理的［10，12］。若考慮風的影響，此時的速度應為相對速度，即碎片速度與風速之差，因此阻力的計算公式為

式中：ρa是空氣密度，CD是空氣阻尼系數(shù)，AD是迎風面積，vre是碎片與風的相對速度，kD是阻力系數(shù)。

在xOy平面內(nèi)僅存在阻力對碎片的影響，在豎直方向則要考慮阻力和重力的雙重影響，因此將式（1）擴展，得到三維體系內(nèi)碎片的加速度方程為

式中：x、y是水平方向坐標，z是豎直方向坐標，t是碎片飛行時間，vw是水平方向恒定風速。碎片在下降階段n＝1，在上升階段n＝2。

參考圖1，給出式（4）～（6）的初始條件，x、y和z方向的初速度和位置分別為

式中：v0是碎片的初速度。

將初始條件（7）～（9）代入方程（4）～（6），利用 MATLAB求解得到

x方向

3 不確定參數(shù)的確定

3.1 碎片初始速度

碎片的初速度

式中：Ek為碎片的初動能，mp為碎片的質(zhì)量。

要計算碎片的動能，需要計算儲罐爆炸時釋放的能量E。目前有很多模型可用來計算氣體儲罐和過熱液體儲罐能量，本文中采用 Baum 模型［6，7，13］

式中：V為儲罐體積，p1為儲罐爆炸時的壓強，p0為大氣壓，γ為比熱容。

儲罐爆炸壓強p1取決于爆炸情況：（1）超壓情況下，爆炸壓強為最大工作壓強×安全系數(shù)；（2）物理爆炸情況下，爆炸壓強為正常的操作壓強；（3）火災情況下，爆炸壓強等于泄壓壓強的1.21倍［13］。由于在實際情況下很難確定儲罐爆炸時的真實壓強，所以儲罐爆炸壓強被當作是隨機變量，通常認為其在上述數(shù)值的90%～110%之間，且服從均勻分布［6，7］，數(shù)學表達式為

式中：p是計算中使用的容器爆炸壓強，Pa；εrand是［0，1］區(qū)間上的隨機數(shù)。

由于儲罐爆炸時，還要產(chǎn)生爆炸沖擊波、撕碎工藝設施等，只有部分能量轉(zhuǎn)化為碎片的能量，因此，碎片的初動能［13］

式中：fk為動能比例因子，其值服從［0.2，0.5］的右三角分布，且均值為0.3，數(shù)學表達式為

3.2 碎片的數(shù)量、質(zhì)量、形狀、能量和拋射角

對46次相關事故的分析表明，柱形罐爆炸產(chǎn)生的碎片數(shù)量N服從對數(shù)正態(tài)分布，其均值為0.855 16，標準差為0.524 48［7，14－16］，數(shù)學表達式為

式中：εrand，1、εrand，2是［0，1］區(qū)間上的2個相互獨立的隨機數(shù)。

根據(jù)P.L.Holden等［14］的結論，碎片質(zhì)量占總質(zhì)量的比例服從貝塔分布，且參數(shù)a＝0.412 13，b＝1.392 6。貝塔分布的概率密度函數(shù)為

式中：x為［0，1］區(qū)間上的隨機變量。

柱形罐爆炸可產(chǎn)生端蓋和彈片2種碎片形式。由于很難精確確定端蓋和碎片的比例，一般采用20%的端蓋和80%的彈片［7］。

假設加載在每個端蓋上的能量相等，加載在第i個彈片上的能量占總動能的比例［7］

式中：I是總彈片數(shù)量；mi是第i個彈片的質(zhì)量。

端蓋和彈片的水平角θ均服從均勻分布［7，14－15］，如圖2所示。其中

20%在［30°，150°］內(nèi)，數(shù)學表達式為

30%在［150°，210°］內(nèi)，數(shù)學表達式為

20%在［210°，330°］內(nèi)，數(shù)學表達式為

30%在［330°，30°］內(nèi)，數(shù)學表達式為

圖2 柱形罐爆炸碎片的水平角分布Fig.2 Horizontal angle distribution of fragment from cylindrical vessel explosion

端蓋的垂直角φ服從［0°，10°］內(nèi)的均勻分布，數(shù)學表達式為

彈片的垂直角φ服從［0°，90°］內(nèi)的均勻分布，數(shù)學表達式為

3.3 阻力系數(shù)

碎片拋射阻力系數(shù)取決于碎片的幾何形狀、表面粗糙度和碎片與拋射方向的夾角，其數(shù)值一般在2.9×10－4～2.1×10－3之間。由于阻力系數(shù)較難確定，R.Pula等［17］采用5種阻力系數(shù)，分別為0，2.9×10－4，6.0×10－4，1.21×10－3和2.0×10－3，計算出5種阻力水平下碎片飛行時間對應的飛行高度及不同初始速度對應的拋射距離，得到拋射距離隨阻力系數(shù)的增加而減小，且當初速度小于35m／s時，阻力作用在碎片拋射距離上的作用非常小，幾乎可以忽略不計。2.9×10－4和6.0×10－4這2種情況對應的阻力系數(shù)較小，接近無阻力，而2.0×10－3適用于超音速的范圍。由于爆炸碎片一般在低音速范圍內(nèi)（200m／s左右），因此建議采用1.21×10－3作為阻力系數(shù)。本文中阻力系數(shù)均采用1.21×10－3。

3.4 儲罐充裝水平

儲罐爆炸能量主要來自2部分：儲罐上方的蒸氣膨脹和儲罐下方的液體汽化膨脹（蒸氣爆炸）。由于前者在儲罐內(nèi)只占很小的一部分，計算時可忽略不計。

液化石油氣汽化膨脹爆炸的能量［10］式中：i1是液化石油氣在儲罐破裂時的壓力或溫度下的焓，i2是液化石油氣在大氣壓下的焓，s1是液化石油氣在儲罐破裂時的壓力或溫度下的熵，s2是液化石油氣在大氣壓下的熵，Tb是液化石油氣在大氣壓下的沸點，ρl是液化石油氣密度，lf是充裝水平；V是儲罐容積。

由式（32）看出，充裝水平lf直接影響儲罐的爆炸能量。lf越高，儲罐爆炸能量越大，碎片的拋射距離也就越遠，反之，lf越低，爆炸能量越小，拋射距離也就越近。一般假設lf在0.1～0.8倍儲罐總容量的范圍內(nèi)，且均勻分布［8－9］，數(shù)學表達式為

4 Monte-Carlo法計算拋射碎片碰撞目標的概率

在儲罐爆炸事故中，計算拋射碎片碰撞目標的概率，需要考慮以下參數(shù)：爆炸時儲罐壓力、儲罐爆炸能量、裝載程度、初始能量轉(zhuǎn)化為碎片能量的比例、碎片數(shù)目、碎片形狀和質(zhì)量、碎片拋射角、阻力系數(shù)、風速和風向等因素，這些參數(shù)都是隨機變量，具有一定的不確定性，具體數(shù)值在一定的區(qū)間范圍內(nèi)。由于參數(shù)存在不確定性，在實際計算中，若用常規(guī)的點估計很難得到一個可信的結果，而Monte-Carlo法采用大量的抽樣對隨機變量和不確定性進行模擬，可以解決參數(shù)的不確定性問題。

為了簡化計算，分析中忽略目標的形狀，認為只要碎片落在目標的外切立方形區(qū)域即算撞擊到目標，這種假設正好可彌補一部分因忽略碎片大小而產(chǎn)生的誤差。當參數(shù)和目標數(shù)據(jù)已知后，用Monte－Carlo方法進行數(shù)值模擬，將這些參數(shù)代入碎片軌跡方程，計算出碎片從拋射到落地間的軌跡曲線，當軌跡與目標的外切立方體有交集時，n取1，當無交集時，n取0。用此方法，模擬Nsim次即可得到碰撞概率Pimp，計算公式為

式中：Nsum為1次模擬中產(chǎn)生的碎片數(shù)目；Vtar為目標的外切立方形區(qū)域；Vfrag為碎片的拋射軌跡；Nsim為模擬次數(shù)。

5 案例分析

以墨西哥城BLEVE事故中盛裝丙烷的水平柱形罐為案例，模擬計算中所需參數(shù)見表1［7］，表中r0、d0、m0、ρ0分別表示柱形罐的半徑、厚度、質(zhì)量和密度。圖3～7分別給出了計算和事故數(shù)據(jù)的比較、模擬的碎片軌跡、有風和無風情況下碎片拋射距離的概率分布、落點分布及無風情況下碎片碰撞目標容器的概率。

表1 柱形罐的參數(shù)值或分布Table 1 Cylindrical vessel’s process parameters

圖3給出了計算和事故數(shù)據(jù)之間的比較?？梢钥闯?，在當前給定的計算參數(shù)下，模擬結果跟實際數(shù)據(jù)一致性較好。圖4是模擬儲罐爆炸碎片被拋射的空間分布情況。表明通過Monte-Carlo方法可以清楚地掌握每個碎片的飛行軌跡，為以后計算奠定基礎。

圖3 計算和事故數(shù)據(jù)間的比較Fig.3 Comparison between the probabilities of fragment ranges to be reached and calculation

圖4 模擬碎片軌跡圖Fig.4 Fragment trajectory of simulation

圖5是在有風和無風情況下碎片拋射距離的概率分布。為了更有代表性，計算采用了2種風速，即±5m／s和±30m／s。“＋”表示風速與爆炸點O到目標I方向相同，“－”表示風速與爆炸點O到目標I方向相反。由圖5可知，當風速為5m／s時，風對碎片拋射距離的影響可忽略不計，是由于風速與碎片速度（200m／s左右）比較相對較小。當風速為30m／s時，對碎片拋射距離的影響比風速為5m／s時大，這表明當風速提高時，風對碎片拋射距離的影響也隨之增加。

圖5 有風和無風情況下碎片拋射距離的概率分布圖Fig.5 Comparison between the probabilities of fragment ranges to be reached with and without wind

圖6是柱形罐爆炸碎片的落點分布。從圖中可清楚地看出，水平圓柱形容器的軸方向產(chǎn)生的碎片數(shù)量相對較多。因此如果目標處在水平圓柱形的［150°，210°］和［330°，360°］這2個區(qū)域內(nèi)具有較大的危險性。在實際工廠布局時應盡量避開這2個區(qū)域，而選擇在［30°，150°］和［210°，330°］這2個區(qū)域內(nèi)會相對安全。

圖7是碎片碰撞一定距離的目標的概率，本文中采用的目標為與初始事故容器一樣的水平圓柱形容器。從圖中看出，碰撞概率隨罐間距的增加而逐漸減小，且碰撞概率Pimp與罐間距R的關系式為：Pimp＝0.035 7e－0.0548R，且相關系數(shù)為0.946 7?？梢钥闯?，碰撞概率隨罐間距的增大而呈現(xiàn)嚴格的指數(shù)衰減趨勢。當R＞70m時，Pimp＜10－3，當R＞100m時，Pimp＜10－4。如果將連鎖效應控制在10－4以內(nèi)，則要求R＞100m，此時100m可作為此安全標準下的罐間所需的安全距離。同樣，如果將連鎖效應控制在10－3以內(nèi)，則要求R＞70m，此時70m可作為此安全標準下的罐間所需的安全距離。

圖7 無風情況下碎片碰撞目標容器的概率圖Fig.7 The impact probabilities of fragment on object with no wind

圖6 落點分布圖Fig.6 Projectile distribution

6 結 論

（1）將二維體系內(nèi)的碎片軌跡方程擴展到三維體系內(nèi)，由此可將風速和風向等影響因素考慮在內(nèi)，并且計算了有風和無風情況下碎片被拋射的距離。通過分析看出，在一般情況下，風對碎片拋射距離的影響較小。

（2）通過模擬得到的碰撞概率與罐間距的關系式，可以計算出符合安全要求的罐間距。這對于提高儲罐的安全性，有效緩解和控制碎片產(chǎn)生的風險，具有一定的指導意義。

（3）只對碎片碰撞目標的概率進行了分析，而未考慮目標被碰撞之后的失效概率。為了定量評價由碎片造成的多米諾風險，目標的失效概率是一個非常重要的因素，值得進一步深入的研究。

［1］劉麗，徐亞博，劉振翼.化工事故多米諾效應定量評價研究［J］.中國安全生產(chǎn)科學技術，2008，4（2）：49-52.

LIU Li，XU Ya-bo，LIU Zhen-yi.Progress on quantitative risk assessment of domino effect in chemical accidents［J］.Journal of Safety Science and Technology，2008，4（2）：49-52.

［2］師立晨，劉驥，魏利軍，等.重大危險源多米諾效應的后果分析［J］.中國安全生產(chǎn)科學技術，2007，3（6）：44-48.

SHI Li-chen，LIU Ji，WEI Li-jun，et al.Consequence analysis for domino effect of major hazardous installations［J］.Journal of Safety Science and Technology，2007，3（6）：44-48.

［3］Cozzani V，Gubinelli G，Antonioni G，et al.The assessement of risk caused by domino effect in quantitative area risk analysis［J］.Journal of Hazardous materials，2005，127（1）：14-30.

［4］Cozzani V，Gubinelli G，Salzano E.Escalation thresholds in the assessment of domino accidental events［J］.Journal of Hazardous Materials，2006，129（1）：1-21.

［5］Khan F I，Abbasi S A.DOMIFFECT（DOMIno eFFECT）：User-friendly software for domino effect analysis［J］.Environmental Modellling ＆ Software，1998，13（2）：163-177.

［6］Hauptmanns U.A Monte-Carlo based procedure for treating the flight of missiles from tank explosions［J］.Probabilistic Engineering Mechanics，2001，16（4）：307-312.

［7］Hauptmanns U.A procedure for analyzing the flight of missiles from explosions of cylindrical vessels［J］.Journal of Loss Prevention in the Process Industries，2001，14（5）：395-402.

［8］Gubinelli G，Zanelli S，Cozzani V.A simplified model for the assessment of the impact probability of fragments［J］.Journal of Hazardous Materials，2004，116（3）：175-187.

［9］楊玉勝，吳宗之.儲罐爆炸碎片最可能拋射距離的Monte-Carlo（蒙特卡羅）數(shù)值模擬［J］.中國安全科學學報，2008，18（3）：15-21.

YANG Yu-sheng，WU Zong-zhi.Monte-Carlo numerical simulation on the most probable flight distance of fragments from explosion of storage tanks［J］.China Safety Science Journal，2008，18（3）：15-21.

［10］邢志祥，蔣軍成，趙曉芳.液化石油氣儲罐爆炸碎片拋射的蒙特卡羅分析［J］.火災科學，2004，13（1）：39-42.

XING Zhi-xiang，JIANG Jun-cheng，ZHAO Xiao-fang.Monte-Carlo analysis of the flight of missiles from LPG tank explosion［J］.Fire Safety Science，2004，13（1）：39-42.

［11］張新梅，陳國華.爆炸碎片拋射速度及飛行軌跡分析方法［J］.華南理工大學學報（自然科學版），2009，37（4）：106-200.

ZHANG Xin-mei，CHEN Guo-hua.Analysis methods of projection velocity and flight trajectory of explosion fragments［J］.Journal of South China University of Technology（Natural Science Edition），2009，37（4）：106-200.

［12］Lees F P.Loss prevention in the process industries［M］.2nd ed.Oxford：Butterworth／Heinemann，1996：73-102.

［13］Center for Chemical Process Safety（CCPS）.Guidelines for evaluating the characteristics of vapor cloud explosions，flash fires，and BLEVEs［M］.New York：American Institute of Chemical Engineers，1994：311-335.

［14］Holden P L，Reeves A B.Fragment hazards from failures of pressurized liquefied gas vessels［C］∥Institution of Chemical Engineers Symposium Series，Assessment and Control of Major Hazards.Manchester，England，1985，93（42）：205-220.

［15］Mébarki A，Mercier F，Nguyen Q B，et al.Structural fragments and explosions in industrial facilities.Part I：Probabilistic description of the source terms［J］.Journal of Loss Prevention in the Process Industries，2009，22：408-416.

［16］Mébarki A，Nguyen Q B，Mercier F.Structural fragments and explosions in industrial facilities.Part II：Projectile trajectory and probability of impact［J］.Journal of Loss Prevention in the Process Industries，2009，22：417-425.

［17］Pula R，Khan F I，Veitch B，et al.A model for estimating the probability of missile impact：Missiles originating from bursting horizontal cylindrical vessels［J］.Process Safety Progress，2007，26（2）：129-139.

Monte-Carlo analysis of the projectile fragments from cylindrical tank boiling liquid expanding vapor explosion accident＊

LIU Zhen-yi，HUANG Ping，XU Ya-bo
（State Key Laboratory of Explosion Science and Technology，Beijing Institute of Technology，Beijing100081，China）

The generation and flight of fragments are stochastic and can cause more severe hazards.Studying on the domino effect study caused by the fragments from explosion is hotspot and difficulty.Based on dynamic analysis，the present work extended the two-dimensional fragment trajectory equation into three-dimensional system，and took wind into account.The present analysis took a cylindrical tank in Mexico City with boiling liquid expanding vapor explosion as an example，and the trajectory of fragment was simulated by Monte-Carlo method.Seen from the probabilities of fragment ranges to be reached with／without wind，wind has less impact on the hazard caused by fragments.The impact probability of a fragment on a target with no wind was attained，impact probability with the tank showed the increase of distance between the strict exponential decay trend，and through this relationship，the tank spacing in compliance with safety requirements can be calculated.The paper provides a method which is of great significance for the risk assessment and control of fragments from the explosion of cylindrical tanks.

mechanics of explosion；boiling liquid expanding vapor explosion；Monte-Carlo method；fragment；cylindrical tank

16July 2009；Revised 14September 2009

LIU Zhen-yi，zhenyiliu＠bit.edu.cn

（責任編輯 曾月蓉）

O389；X937 國標學科代碼：130·35

A

2009-07-16；

2009-09-14

國家科技重大專項基金項目（2008ZX05054）

劉振翼（1975— ），男，博士后，講師。

1001-1455（2010）06-0569-08

Supported by the National S＆T Major Project（2008ZX05054）