丁宇奇 葉碧濤 蘆 燁 成佳浩 盧 宏 謝 清 王學勇

（東北石油大學機械科學與工程學院）

隨著石化工業(yè)的發(fā)展，近年來能源危機不斷加大，國家原油戰(zhàn)略儲備庫項目的實施迫使我國不得不建造更多更大的儲罐。 儲罐儲存的液體或者氣體多為易燃易爆物品，容易發(fā)生爆炸事故［1～3］。 為了減小爆炸事故造成的傷害，需要對儲罐內(nèi)爆載荷與儲罐結構破壞之間的關系進行研究。

目前， 針對爆炸載荷計算方法的研究有很多。在理論計算方面，最常見的方法是以TNT炸藥為計算標準，將氣云爆熱等效為TNT當量來計算，即TNT當量法［4］。張麗運用TNT當量法計算了不同燃氣性質(zhì)、不同泄漏量以及不同泄漏口徑等情況下燃氣泄漏爆炸產(chǎn)生的沖擊波值，并預測了可能對結構產(chǎn)生的破壞［5］。 丁宇奇等依據(jù)TNT當量法推導出了儲罐內(nèi)爆超壓分布曲線，得到了不同起爆點和不同液位會影響儲罐破壞位置的結論［6］。然而， 由于爆炸時的物理和化學過程相當復雜，致使理論求解范圍極其有限，且理論計算無法觀測到破口的形貌。 而數(shù)值模擬在獲得爆炸過程中各種動態(tài)物理量的同時，還可以得到更加全面且足夠精確的計算結果，因此得到了廣泛的應用［7，8］。常用的確定容器內(nèi)爆載荷的數(shù)值模擬方法主要有兩種：計算流體動力學（CFD）方法和TNT當量法，已有不少學者對兩種方法進行了對比研究。 張秀華等建立了相同裝藥形狀的乙炔氣團和當量的TNT兩種爆炸源模型， 在分析兩種爆源產(chǎn)生的超壓曲線后， 得到了隨著距爆源距離的增加，兩種爆炸壓力波形間的誤差逐漸減小的結論［9］。但是，這種方法是單純的將乙炔濃縮為高能爆炸氣團，沒有考慮乙炔與空氣的化學反應。 為此， 蘆燁等考慮可燃氣云爆炸時的化學反應，分別采用TNT當量法和CFD方法研究了二維儲罐模型內(nèi)的壓力和速度變化，研究結果表明：CFD方法能更好地描述儲罐內(nèi)的壓力場分布［10］。 然而，該研究僅分析了罐內(nèi)流場，沒有對爆炸載荷作用下儲罐結構破壞進行研究。 實驗方法可以直接有效地獲取結構破壞研究所需的數(shù)據(jù)，姚術健等通過多次箱型結構內(nèi)爆實驗，得出了箱體在內(nèi)爆作用下會出現(xiàn)板中心沖切破口、箱體角隅處破壞及板邊緣撕裂等多種破壞形式的結論［11］。Vanessa P等開展了金屬容器內(nèi)爆實驗，分別使用金相分析和高速攝像技術對結構破壞進行了分析，指出結構破壞是由焊接熱影響區(qū)域應變過大和材料強度降低導致的［12］。 然而，進行爆炸實驗研究的成本非常高，同時實驗過程存在一定的危險性。

綜上所述，數(shù)值模擬方法是開展爆炸載荷研究的最佳方法，但是目前鮮有人針對儲罐在不同載荷等效方法下的結構破壞進行數(shù)值模擬研究。因此，筆者以立式拱頂儲罐為研究對象，考慮儲罐角鋼結構和壁板變壁厚的特性， 分別使用CFD方法和TNT當量法對儲罐結構破壞進行數(shù)值計算。 通過分析儲罐結構的應力、變形和破壞位置，得到兩種方法在不同容量儲罐和不同爆心距下的儲罐破壞形貌， 為研究儲罐內(nèi)可燃氣云爆炸時，兩種等效方法在載荷計算和結構破壞時的誤差范圍提供參考。

1 儲罐內(nèi)可燃氣云爆炸載荷不同等效方法分析及模型建立

1.1 儲罐內(nèi)可燃氣云爆炸載荷CFD方法及模型建立

可燃氣云爆炸是一種帶有劇烈化學反應的過程，表現(xiàn)出來的形態(tài)為逐層燃燒，火焰從點火點逐步擴散，火焰在擴散過程中加速并形成爆轟波。 爆轟波的傳播速度為每秒數(shù)千米，燃燒過程伴有強烈的壓力升高［13］。 筆者在使用CFD方法模擬儲罐爆轟過程時不考慮溫度作用和氣體的泄漏， 此時儲罐內(nèi)的爆炸壓力可視為定容爆炸壓力。 氣體定容爆炸壓力pm與初始環(huán)境壓力p0有以下關系［14］：

其中，n0為氣體初始摩爾量；T0為反應初始溫度；nf為反應終態(tài)摩爾量；Tf為反應終態(tài)溫度。為了使儲罐內(nèi)可燃氣云燃燒發(fā)展成爆轟波，以達到儲罐爆炸破壞效果，筆者采用設置高溫點火源的方法驅動爆轟波形成。 為簡化計算，將儲罐內(nèi)部可燃混合氣體視為理想氣體，混合氣體滿足理想氣體狀態(tài)方程。 使用Fluent流體仿真軟件來模擬可燃氣體燃燒爆炸過程。 在保證較好模擬湍流擴散燃燒反應的前提下， 采用計算量較小的EDC燃燒模型，同時使用k-ε湍流模型實現(xiàn)燃燒時的湍流流動［15］。

筆者分析的儲罐容積為3 000、5 000、10 000m3，三者直徑分別為15、20、27.5m，罐壁高均為17.82m，罐頂曲率半徑為直徑的1.2倍，壁板厚度由上到下逐漸增大。 可燃氣體點火點的位置選擇在儲罐1/2高度，距離罐壁5m（N1）、3m（N2）、1m（N3）處。并選擇最先接觸爆炸波位置設置測點A， 另外在罐頂與罐壁連接處位置設置測點B、C。因為該位置結構不連續(xù)， 同時厚度也是最薄弱的。限于篇幅，本節(jié)僅以5 000m3儲罐做詳細分析。根據(jù)5 000m3罐尺寸在Fluent中設置儲罐內(nèi)可燃氣體域。在初始化設置上，混合氣體初始溫度300K，并在點火位置設置一個2 500K的高溫區(qū)域以及少許反應產(chǎn)物。 可燃氣體濃度取乙炔反應最劇烈的濃度7.4%，5 000m3儲罐流體域有限元模型如圖1所示。

圖1 CFD分析儲罐有限元模型

1.2 儲罐內(nèi)可燃氣云爆炸載荷TNT當量法及模型建立

由于可燃氣云爆炸的復雜性，通常將可燃氣云爆炸的能量等效成當量的TNT來研究爆炸對結構的破壞［16］。使用TNT當量法模擬儲罐內(nèi)爆時，巨大的爆炸能量會在瞬間全部釋放，產(chǎn)生大量的高溫高壓爆炸產(chǎn)物，急劇膨脹的爆炸產(chǎn)物會壓縮周圍的空氣并使之偏離原有位置，形成一層高壓空氣波。 隨著爆炸氣體的繼續(xù)膨脹，其體積的增大帶來的后果就是壓力的不斷減小，最終爆炸波壓力將與大氣壓相同。TNT爆炸壓力p的計算一般使用JWL狀態(tài)方程［17］：

其中，V為爆轟產(chǎn)物相對比容；E為爆轟產(chǎn)物的比內(nèi)能；A、B、R1、R2、ω均為JWL方程參數(shù)。 爆炸威力取決于可燃氣體換算的TNT當量，TNT當量計算公式如下：

其中，WTNT為可燃氣體的TNT當量；α為可燃氣體的效率因子，表明參與爆炸的可燃氣體的體積分數(shù)，其中乙炔的效率因子為19%；Wf為可燃氣體的總質(zhì)量；Qf為可燃氣體的燃燒熱；QTNT為TNT的爆炸熱，一般取4.61×103kJ/kg。 將3 000、5 000、10 000m3儲罐內(nèi)可燃氣體換算后， 得到TNT當量分別為276、498、962kg。

與CFD方法中流體域位于結構模型內(nèi)部不同，TNT當量法采用CEL方法計算，因此建立的流體域要大于結構模型，以保證結構發(fā)生變形破壞后數(shù)據(jù)能繼續(xù)傳遞。 為了得到較好的網(wǎng)格質(zhì)量，使用VOF法指派TNT為球形裝藥并設置初始空氣體積分數(shù)， 得到TNT當量法分析儲罐有限元模型如圖2所示。

圖2 TNT當量法分析儲罐有限元模型

2 不同等效方法下的儲罐內(nèi)爆載荷計算

2.1 CFD方法計算儲罐內(nèi)可燃氣云爆炸載荷

由于筆者采用的CFD方法無法模擬儲罐破壞后罐內(nèi)壓力的泄放， 因此在使用CFD方法計算儲罐內(nèi)可燃氣云爆炸載荷時，僅提取爆炸破壞前測點的壓力。 限于篇幅，本節(jié)僅提取5 000m3儲罐中N2點火點的工況數(shù)據(jù)。爆炸波剛接觸到罐壁時的壓力分布如圖3所示，A、B、C點的壓力隨時間變化曲線如圖4所示。

圖3 CFD方法爆炸壓力分布圖

圖4 CFD方法測點壓力曲線

由圖3可以看出，使用CFD方法模擬燃燒爆炸時，可燃氣云燃燒形成了一個球形波，而且火焰接觸到儲罐罐壁時可燃氣體僅發(fā)生了部分燃燒。 由于火焰擴散是一個逐層燃燒的過程， 因此可以看到球形波外層的壓力要大于內(nèi)層。 從圖4可以看出，A點在15ms時檢測到了壓力，壓力數(shù)值在17ms時達到峰值1.53MPa， 隨后逐漸下降直到穩(wěn)定在0.40MPa左右。B、C點由于距離起爆點較遠，因此測得壓力的時間相對滯后。但是B點測得的峰值壓力（2.01MPa） 為A點的1.31倍，C點測得的峰值壓力（3.31MPa）為A點的2.16倍。 由此可以看出，使用CFD方法模擬爆炸時，距離點火點越遠，爆炸壓力越大。這是因為在燃燒初始階段，雖然爆炸波先到達近爆點A，但是此時罐內(nèi)可燃氣體并沒有完全燃燒。 隨著時間的推移，可燃氣體邊燃燒邊擴散，壓力波依次到達B點和C點， 此時罐內(nèi)可燃氣體絕大部分已經(jīng)發(fā)生反應，產(chǎn)生的壓力也就更大。因此B、C點測得的爆炸壓力也就更大，最終出現(xiàn)了測點距離起爆點越遠爆炸壓力越大的現(xiàn)象。

2.2 TNT當量法計算儲罐內(nèi)可燃氣云爆炸載荷

TNT當量模型的工況設置與CFD方法相同，計算后僅提取N2起爆點的爆炸數(shù)據(jù)，爆炸波剛接觸到罐壁時的壓力分布如圖5所示， 各測點的壓力變化曲線如圖6所示。

圖5 TNT方法爆炸壓力分布圖

圖6 TNT方法測點壓力曲線

從圖5可以看出，TNT當量法模擬出來的爆炸壓力以環(huán)形波形式擴散，爆炸發(fā)生后中心區(qū)域的壓力與壓力波未到達的區(qū)域壓力數(shù)值相同。 由于爆炸波接觸到罐壁后形成了反射波，因此罐壁附近的壓力數(shù)值較大。 與CFD方法不同，圖6所示的A點在1.8ms時就檢測到了爆炸壓力，并且壓力峰值達到了9.78MPa，是CFD方法的6.39倍，并且壓力到達峰值之后在2ms內(nèi)迅速下降為零。B點的峰值壓力為8.15MPa，較A點下降了20%。 壓力波在到達C點之前結構已經(jīng)發(fā)生破壞，因此C點未檢測到爆炸壓力。由此可以得出結論，使用TNT當量法計算儲罐內(nèi)可燃氣云爆炸載荷時，測點離起爆中心越遠爆炸壓力越小。

綜上所述，使用CFD方法時，爆炸壓力隨著距爆炸中心距離的增大而增大；使用TNT當量法時，爆炸壓力隨著距爆炸中心距離的增大而減小。 由于CFD方法和TNT當量法本質(zhì)上的不同，導致在相同工況下，二者得到了完全不同的壓力變化規(guī)律。

3 不同等效方法下的儲罐內(nèi)爆結構破壞分析

3.1 CFD方法儲罐內(nèi)爆結構破壞分析

使用CFD方法研究儲罐結構破壞時需要借助Abaqus來進行模擬。 使用Mpcci來完成Fluent與Abaqus之間的數(shù)據(jù)傳遞， 由此實現(xiàn)CFD方法模擬儲罐結構在可燃氣云爆炸載荷下的破壞。 限于篇幅，此處僅出具爆源距罐壁1m時不同容積儲罐的破壞位置（圖7）。

圖7 不同容積儲罐的破壞位置

從圖7a可以看出，整個儲罐壁板都產(chǎn)生了應力和變形， 儲罐破壞位置位于左側頂壁連接處。由2.1節(jié)可知， 爆炸壓力隨著爆心距的增大而增大，因此3 000m3儲罐左側頂壁連接處承受的爆炸壓力最大。 同時，儲罐頂壁連接處屬于結構不連續(xù)位置，容易發(fā)生應力集中。 因此，頂壁連接處較儲罐其他位置更容易發(fā)生破壞。 對比圖7可以發(fā)現(xiàn)，隨著儲罐容積的增大，破口位置逐漸由左向右轉移。 這是因為隨著儲罐容積的增大，罐內(nèi)可燃氣體的質(zhì)量也在增大，氣體燃爆產(chǎn)生爆炸壓力也就更大。 此時，燃爆產(chǎn)生的球形火焰波還未到達最左側頂壁連接處，燃爆的壓力已足以使結構發(fā)生失效破壞。 因此，隨著容積的增大，破口位置不再出現(xiàn)在最左側頂壁連接處，并逐漸向右側轉移。 同時，儲罐的變形呈現(xiàn)上大下小的葫蘆狀，這是由于本課題中儲罐模型采用了變壁厚的設置，越靠近罐頂位置、罐壁厚度越小，所以儲罐上半部的變形比下半部更大。 為了避免儲罐結構大變形引起的單元畸變，本次計算使用J-C強度模型和相應的失效準則［18］，公式如下：

其中A′、B′、C、n、m為模型參數(shù)；εeq為等效塑性應變；σeq為等效應力；εf為失效應變；為等效塑性應變率；T*為無量綱溫度；D1～D5是材料常數(shù)；σ*為材料應力三軸度。 由此可見，J-C強度模型和失效準則均與應變率相關，等效應力和失效應變是隨著應變率的變化而不斷改變的。 筆者在統(tǒng)計儲罐破壞應力時發(fā)現(xiàn)，儲罐破壞伴隨著單元的刪除，單元刪除的瞬間應力變化幅度很大。 為了分析起爆點不同時儲罐結構破壞的規(guī)律，試驗中提取的數(shù)據(jù)為破壞前0.5ms時的數(shù)據(jù)，得到儲罐的最大應力及徑向變形（表1）。

表1 不同容積儲罐應力及徑向變形量

由表1中3 000m3儲罐破壞數(shù)據(jù)可以看出，隨著點火點越靠近罐壁，儲罐結構破壞時間逐漸增加，徑向變形也越來越大。 由于點火點越靠近罐壁，燃燒產(chǎn)生的火焰波到達左側頂壁連接處的時間越長。 同時，爆炸壓力也隨著時間的增加而增大，因此變形也就隨之增大。 由式（4）、（5）可知，結構破壞時的應力與應變率相關，而應變率的大小取決于爆炸壓力作用在結構上時間的快慢。 不同工況下儲罐破壞時的應變率不同，因此儲罐破壞時的應力也就不同。 在儲罐容積增大的過程中，隨著罐內(nèi)可燃氣體量的大幅增長，儲罐的應力、變形都沒有很大變化，這是由不同容積的儲罐壁厚不同導致的。 為了保證結構的強度，隨著儲罐容量的增加， 儲罐的壁厚也會相應增大，壁厚的增加限制了儲罐應力和變形的大幅提升。

3.2 TNT當量法儲罐內(nèi)爆結構破壞分析

綜合考慮儲罐罐體結構、罐體內(nèi)外空氣域和可燃氣體的TNT當量質(zhì)量， 模型中爆炸產(chǎn)物與結構之間采用CEL流固耦合算法， 起爆點位置與CFD方法采用相同設置， 最終得到不同起爆點時儲罐應力分布及破壞位置（圖8）。

由圖8a、b可以看出，3 000、5 000m3儲罐的破壞位置均為近爆側頂壁連接處，同時儲罐罐壁位置出現(xiàn)了較大凸起。這是由于TNT爆炸時，距離爆炸中心越近，結構承受的壓力越大，爆炸產(chǎn)生的超壓使儲罐罐壁發(fā)生局部大變形。 當爆炸波繼續(xù)傳遞到達右側頂壁連接處時， 由于該位置是結構不連續(xù)區(qū)域，因此該位置先發(fā)生破壞。 而從圖8c可以看出10 000m3儲罐爆炸破壞位置為罐壁。這是因為10 000m3儲罐換算的TNT當量質(zhì)量足夠大，雖然同樣是在距離罐壁1m位置處爆炸， 但10 000m3當量TNT爆炸產(chǎn)生的超壓足以使儲罐罐壁直接被破壞。為了得到更準確的儲罐破壞時刻，采用TNT當量法分析時設置的時間間隔為0.1ms。 提取3種不同容積的儲罐在不同起爆位置破壞時的參數(shù)，具體見表2。

圖8 不同容積儲罐應力分布及破壞位置

表2 不同容積儲罐應力及徑向變形

從表2可以看出， 當儲罐容積相同時，TNT起爆位置每靠近罐壁2m， 徑向變形增量超過80%。這是由于TNT爆炸產(chǎn)生的超壓在空氣中損失很大，爆炸沖擊波每向前傳遞1m，爆炸壓力就會大幅下降。 觀察不同容積儲罐的破壞時間可以發(fā)現(xiàn)，隨著容積的增大，爆炸破壞時間變短，這是由于TNT當量變大后，爆炸瞬間產(chǎn)生的壓力越大，破壞時間越提前。

3.3 不同等效方法下儲罐內(nèi)爆結構破壞對比分析

為了分析CFD方法和TNT當量法各測點位置結構應力變化情況，對5 000m3儲罐，在距罐壁5m起爆位置工況中， 分別提取兩種方法A、B、C點的壓力和應力時程曲線如圖9、10所示。

圖9 不同方法壓力時程曲線

圖10 不同方法應力時程曲線

由圖9可以看出，TNT當量法在A點和B點得到的壓力遠大于CFD方法。 TNT當量法的爆炸壓力波經(jīng)過A點的壓力峰值為9.8MPa，CFD方法的壓力波經(jīng)過A點的壓力峰值僅為1.9MPa， 二者相差超過80%。 B點的壓力峰值也相差了近74%。 C點由于在儲罐破壞前沒有檢測到壓力，因此僅有CFD方法的壓力曲線。通過觀察圖10可以發(fā)現(xiàn)，兩種方法的A點應力變化趨勢相同， 但是CFD方法的應力出現(xiàn)時間和上升速度都較小。 TNT當量法的B點出現(xiàn)了單元刪除，因此應力值在發(fā)生劇烈增大后下降為零，CFD方法的B點與A點趨勢一致，僅應力出現(xiàn)時間略有增大。 兩種方法C點位置的單元沒有發(fā)生失效，因此CFD方法的應力在結構破壞前持續(xù)上升，而TNT當量法未檢測到壓力，也就沒有產(chǎn)生應力數(shù)值。 為了分析不同爆炸模擬方法對儲罐結構破壞位置的影響，分別提取CFD方法和TNT當量法模擬后的儲罐破壞位置，具體見表3。

表3 不同容積儲罐破壞位置

由表3可以看出，相同容積的儲罐采用CFD方法時破壞位置相同， 而且隨著儲罐容積的增大，破壞位置逐漸由左向右轉移。 TNT當量法計算得到的儲罐破壞位置大部分位于右側頂壁連接處，而10 000m3儲罐的破壞位置發(fā)生了改變， 起爆點距罐壁1m時儲罐的破壞位置為近爆點罐壁。破壞位置的改變是不同容積儲罐的TNT當量不同導致的。采用CFD方法得到的儲罐破壞時間差別不大，最大時長與最小時長僅相差25%； 采用TNT當量法得到的儲罐破壞時間相差較大，最大時長與最小時長相差超過72%。

4 結論

4.1 采用CFD方法和TNT當量法研究儲罐內(nèi)爆載荷時，CFD方法的爆炸載荷傳遞速度要低于TNT當量法。 在同一測點TNT當量法檢測到壓力的時間僅為CFD方法的1/10，同時TNT當量法近爆測點A的峰值壓力是CFD方法的6.39倍。 使用CFD方法時， 爆炸壓力隨著爆炸中心的距離增大而增大；使用TNT當量法時， 爆炸壓力隨著爆炸中心的距離增大而減小。

4.2 采用CFD方法和TNT當量法研究儲罐結構破壞時，CFD方法計算出來的結構破壞位置位于遠離起爆點的左側頂壁連接處，并且隨著儲罐容積的增大逐漸由左向右轉移。 TNT當量法多數(shù)工況儲罐破壞位置位于近爆點頂壁連接處。 當TNT當量質(zhì)量足夠大或者TNT與罐壁距離足夠近時，儲罐在近爆點罐壁發(fā)生破壞。 采用CFD方法得到的儲罐破壞時間差別不大，最大時長與最小時長僅相差25%； 采用TNT當量法得到的儲罐破壞時間相差較大，最大時長與最小時長相差超過72%。

4.3 由于可燃氣體逐層燃燒的特性，CFD方法的爆炸壓力隨著爆心距的增大而增大，使得距離點火點較遠處的結構承受較大的壓力，因此遠爆點處的結構先發(fā)生破壞。 而TNT當量法計算爆炸載荷時，爆炸壓力是瞬間釋放并逐漸降低的，此時近爆點結構會承受更大載荷并發(fā)生破壞。 正是兩種等效方法本質(zhì)上的不同導致了模擬結果的規(guī)律不一樣。 CFD方法模擬的爆炸在結構與爆源較近時爆炸強度較低， 此時TNT當量法模擬可燃氣云爆炸產(chǎn)生的超壓峰值最大相差80%。 當結構與爆源相距較遠時，兩種方法模擬出的結構破壞趨于一致。