丁宇奇，盧 宏，蘆 燁，謝 清，楊 明，張佳賀，劉凱

(東北石油大學(xué) 機(jī)械科學(xué)與工程學(xué)院，黑龍江大慶 163318)

0 引言

大型儲(chǔ)罐通常存儲(chǔ)易揮發(fā)、易燃、易爆物料，揮發(fā)出的氣體與空氣混合容易發(fā)生燃燒、爆炸，其爆炸后產(chǎn)生的巨大能量會(huì)導(dǎo)致儲(chǔ)罐自身結(jié)構(gòu)發(fā)生破壞，若罐體破壞后產(chǎn)生碎片將可能引發(fā)多米諾效應(yīng)[1]，對(duì)鄰近儲(chǔ)罐及建筑物造成破壞。因此，對(duì)儲(chǔ)罐內(nèi)爆破碎及碎片形成后的運(yùn)行軌跡研究就顯得尤為重要。

儲(chǔ)罐內(nèi)爆破碎問(wèn)題主要涉及爆炸載荷、結(jié)構(gòu)響應(yīng)、碎片的形成以及拋射軌跡等方面。在爆炸載荷模擬方面，管義鋒等[2]采用TNT當(dāng)量法對(duì)燃?xì)庑孤┮l(fā)的爆炸事故進(jìn)行分析，通過(guò)理論計(jì)算得到不同位置的沖擊波超壓以及危害后果，與實(shí)際爆炸后果一致；胡可等[3]采用CFD方法對(duì)大型儲(chǔ)罐內(nèi)甲烷/空氣預(yù)混氣體的爆炸過(guò)程進(jìn)行模擬，得到了不同容積、可燃?xì)怏w種類與濃度等因素對(duì)儲(chǔ)罐結(jié)構(gòu)爆炸載荷的影響。在爆炸載荷作用下的結(jié)構(gòu)響應(yīng)方面，金澤宇等[4]基于ALE流固耦合算法，對(duì)球殼在TNT炸藥爆炸載荷下的結(jié)構(gòu)響應(yīng)開(kāi)展研究，得到了流體的響應(yīng)以及球殼的抗沖擊性能；胡陳[5]考慮多相耦合作用，建立了拱頂罐爆炸數(shù)值分析模型，得到了炸藥當(dāng)量、炸藥位置及液位對(duì)儲(chǔ)罐結(jié)構(gòu)動(dòng)態(tài)響應(yīng)的影響規(guī)律；DING等[6]以立式拱頂儲(chǔ)罐為研究對(duì)象，采用有限元計(jì)算方法計(jì)算了儲(chǔ)罐爆炸過(guò)程中流體的壓力和溫度，得到了儲(chǔ)罐的應(yīng)力和變形分布。在爆炸碎片形成的數(shù)值模擬方面，張志彪[7]以TNT炸藥為爆源，通過(guò)以斷裂應(yīng)變控制的單元?jiǎng)h除法，得到了殼體碎片分布，不過(guò)單元?jiǎng)h除法會(huì)造成罐體材料的質(zhì)量損失，而且破片數(shù)量也偏少;董新剛等[8]通過(guò)損傷模型控制材料的破壞，采用光滑粒子動(dòng)力學(xué)(SPH)方法對(duì)TNT爆炸沖擊作用下絕熱擋環(huán)的破碎性能開(kāi)展研究，獲得了有無(wú)削弱槽結(jié)構(gòu)的破碎開(kāi)裂形式和碎片尺寸，但是光滑粒子法中產(chǎn)生的隨機(jī)破片形態(tài)并不理想，無(wú)法確定大碎片的邊界。在爆炸碎片拋射軌跡求解方面，錢新明等[9-10]采用蒙特卡洛方法對(duì)儲(chǔ)罐爆炸碎片的拋射開(kāi)展研究，得到了碎片的軌跡曲線、拋射距離的概率分布等。雖然蒙特卡洛方法模擬的結(jié)果與實(shí)際數(shù)據(jù)一致性較好，但該方法基于很多理想化的條件，如爆炸儲(chǔ)罐被視為一個(gè)質(zhì)點(diǎn)、忽略儲(chǔ)罐的體積等，而且也無(wú)法定量揭示爆炸沖擊對(duì)結(jié)構(gòu)的動(dòng)力響應(yīng)以及碎片的產(chǎn)生過(guò)程；王文紅[11]考慮流固耦合作用，通過(guò)理論估算方法得到碎片的拋射距離，定量分析了爆炸碎片的影響范圍及在影響區(qū)域內(nèi)對(duì)設(shè)備、人員造成的傷害程度。

綜上，現(xiàn)有關(guān)于儲(chǔ)罐內(nèi)爆破碎研究中，還未充分考慮爆炸載荷與罐體之間的多相耦合問(wèn)題，且不能準(zhǔn)確刻畫(huà)碎片產(chǎn)生后的形貌與質(zhì)量。為此，以球形儲(chǔ)罐為研究對(duì)象，考慮內(nèi)爆流體-球罐(碎片)-罐外氣體的多相耦合以及爆炸載荷的傳遞，建立球罐內(nèi)爆三維有限元模型；基于指數(shù)內(nèi)聚力失效準(zhǔn)則建立儲(chǔ)罐材料粘結(jié)單元的失效條件，既可準(zhǔn)確刻畫(huà)球罐內(nèi)爆碎片的形成過(guò)程，又可保證內(nèi)爆前后的質(zhì)量守恒；通過(guò)與蒙特卡洛方法的對(duì)比分析，驗(yàn)證本文數(shù)值模擬計(jì)算結(jié)果的準(zhǔn)確性，該研究成果可為儲(chǔ)罐內(nèi)爆事故預(yù)防與事故救援提供理論依據(jù)。

1 基于多相耦合的儲(chǔ)罐內(nèi)爆載荷傳遞與碎片軌跡計(jì)算方法

1.1 儲(chǔ)罐內(nèi)爆載荷傳遞方法

GUBINELLI等[12]綜合分析了143宗儲(chǔ)罐爆炸碎片拋射事故，得出臥式儲(chǔ)罐事故數(shù)量最多，球形儲(chǔ)罐次之。為此，本文以球形儲(chǔ)罐為例開(kāi)展研究，其內(nèi)爆并產(chǎn)生爆炸碎片的過(guò)程如圖1所示。罐內(nèi)氣體、罐內(nèi)液體、球罐和罐外空氣4個(gè)區(qū)域之間復(fù)雜的載荷數(shù)據(jù)需要通過(guò)耦合面實(shí)現(xiàn)傳遞，在球罐內(nèi)爆過(guò)程中涉及了罐內(nèi)氣體-球罐、球罐-罐外空氣、罐內(nèi)液體-球罐等多個(gè)耦合界面。通過(guò)CEL[13](耦合歐拉-拉格朗日)方法對(duì)不同相之間的耦合關(guān)系進(jìn)行描述，其中流體和結(jié)構(gòu)之間采用通用接觸算法，以達(dá)到計(jì)算中自動(dòng)追蹤結(jié)構(gòu)與流體之間的耦合接觸面的目的。

圖1 球罐內(nèi)爆過(guò)程示意

球罐內(nèi)部可燃?xì)怏w發(fā)生爆炸后產(chǎn)生的沖擊波在罐內(nèi)流體中傳播，在圖1中區(qū)域①位置，沖擊波通過(guò)罐內(nèi)氣體-球罐耦合界面后作用到球罐結(jié)構(gòu)上，球罐便發(fā)生變形破壞。球罐的響應(yīng)一方面通過(guò)罐內(nèi)氣體-球罐耦合界面將變形等數(shù)據(jù)反饋給罐內(nèi)氣體，另一方面通過(guò)球罐-罐外空氣耦合界面?zhèn)飨蚬尥饪諝?。每個(gè)時(shí)間步流體與結(jié)構(gòu)之間通過(guò)自動(dòng)追蹤的接觸面?zhèn)鬟f數(shù)據(jù)，區(qū)域①中耦合界面的載荷傳遞如圖2所示，(圖中，虛線表示耦合界面變形前的位置，實(shí)線表示耦合界面變形后的位置)。

圖2 球罐局部載荷傳遞示意

爆炸沖擊波使罐內(nèi)氣體產(chǎn)生壓力，該壓力通過(guò)罐內(nèi)氣體作用在球罐結(jié)構(gòu)上，如式(1)所示。

Pf=HfsPs1

(1)

式中,Pf為罐內(nèi)氣體壓力向量，MPa；Hfs為罐內(nèi)氣體域到球罐結(jié)構(gòu)域的傳遞矩陣；Ps1為球罐結(jié)構(gòu)內(nèi)表面的壓力向量，MPa。

在球罐壓力Ps1的作用下，球罐結(jié)構(gòu)內(nèi)表面產(chǎn)生變形并反饋給罐內(nèi)氣體，如式(2)所示。

δs1=Hfs1δf

(2)

式中,δs1為球罐結(jié)構(gòu)內(nèi)表面的位移向量,mm；Hfs1為球罐結(jié)構(gòu)域到罐內(nèi)流體域的傳遞矩陣；δf為球罐內(nèi)表面變形后罐內(nèi)流體邊界處的位移向量[14],mm。

在球罐結(jié)構(gòu)內(nèi)表面變形δs1作用下，球罐外表面產(chǎn)生變形并傳遞到罐外流體域，如式(3)所示。

δs2=Hf ′s2δf ′

(3)

式中,δs2為球罐結(jié)構(gòu)外表面的位移向量,mm；Hf ′s2為球罐結(jié)構(gòu)域到罐外流體域的傳遞矩陣；δf ′為球罐外表面變形后罐外流體邊界處的位移向量,mm。

罐外氣體域的壓力Pf′為大氣壓，直接作用到球罐結(jié)構(gòu)外表面上，如式(4)所示。

Pf ′=Hf ′s2Ps2

(4)

式中,Pf ′為罐外氣體壓力向量，MPa；Hf ′s2為罐外氣體域到球罐結(jié)構(gòu)域的傳遞矩陣；Ps2為球罐結(jié)構(gòu)外表面的壓力向量，MPa。

1.2 爆炸碎片形成后的拋射軌跡研究

1.2.1 爆炸碎片運(yùn)動(dòng)軌跡的坐標(biāo)變換

爆炸碎片產(chǎn)生后，隨著爆炸能量的消耗，碎片的速度會(huì)達(dá)到一個(gè)最大值，將該時(shí)刻作為碎片拋射的初始時(shí)刻，此時(shí)碎片在空氣中拋射的初始速度為v0。若僅考慮風(fēng)阻對(duì)碎片軌跡的影響，碎片將在一個(gè)平面內(nèi)運(yùn)動(dòng)，可將受阻力作用下的三維空間碎片拋射問(wèn)題簡(jiǎn)化為二維平面問(wèn)題，如圖3所示。

圖3 碎片二維坐標(biāo)系的轉(zhuǎn)化過(guò)程

在三維空間坐標(biāo)系中，vx′,vy′,vz′分別為v0在x′,y′,z′方向的分量;vx′y′為v0在x′o′y′面的投影。碎片在vz′,vx′y′所在平面內(nèi)運(yùn)動(dòng)，將空間坐標(biāo)系轉(zhuǎn)換為該平面內(nèi)的二維坐標(biāo)系，其中vz′=vy，vx′y′=vx，以下在碎片二維平面坐標(biāo)系中對(duì)碎片拋射軌跡進(jìn)行求解。

1.2.2 局部坐標(biāo)系下的碎片拋射

碎片在二維坐標(biāo)系中的拋射情況如圖4所示。在整體坐標(biāo)系XOY中，初始時(shí)刻碎片坐標(biāo)為(X0，Y0)，為了便于計(jì)算，在初始時(shí)刻碎片位置處建立局部坐標(biāo)系xoy；局部坐標(biāo)系中，碎片初始速度v0，拋射角度為θ(θ在第一象限為正值，在第二象限為負(fù)值)。

圖4 二維坐標(biāo)系下的碎片拋射

局部坐標(biāo)系下碎片初始條件：

(5)

(6)

碎片在空氣中運(yùn)動(dòng)總是會(huì)受到與速度成正相關(guān)的空氣阻力影響，在考慮拋體速率較低時(shí)，可將空氣阻力的大小近似為與速率成正比，即：

(7)

其中：

(8)

式中,k為比例系數(shù);ρa(bǔ)為空氣密度；S為迎風(fēng)面積；C為空氣阻力系數(shù)。

于是在二維坐標(biāo)系中，碎片在空氣阻力作用下的運(yùn)動(dòng)學(xué)方程可寫(xiě)為：

(9)

(10)

將初始條件式(5)(6)代入式(9)(10)，解出碎片的運(yùn)動(dòng)學(xué)方程表達(dá)式為：

(11)

(12)

對(duì)式(11)(12)中時(shí)間t求導(dǎo)，得到碎片運(yùn)動(dòng)的速度表達(dá)式：

(13)

(14)

令式(14)等于零，得到碎片運(yùn)動(dòng)到最高點(diǎn)的時(shí)間為：

(15)

再將式(15)代入式(12)，得到碎片的拋射高度為：

(16)

同時(shí)，由式(11)(12)可得碎片的軌跡方程：

(17)

1.2.3 球罐整體坐標(biāo)系下的碎片拋射

在球罐整體坐標(biāo)系XOY與碎片局部坐標(biāo)系xoy之間存在如下關(guān)系：

X=X0+x，Y=Y0+y

(18)

將式(18)代入軌跡方程(17)，得到球罐整體坐標(biāo)系下的軌跡方程為：

(19)

在球罐整體坐標(biāo)系下，令式(19)中Y=0，可以得到碎片的拋射距離。同時(shí)，通過(guò)式(16)、(18)可以得到碎片的拋射高度為：

H=Y0+h

(20)

2 儲(chǔ)罐內(nèi)爆碎片破壞條件的建立

內(nèi)爆引發(fā)的儲(chǔ)罐破壞包含多種表現(xiàn)形式，其中剪切失效是材料在沖擊載荷下與動(dòng)態(tài)損傷演化-破壞相關(guān)的一個(gè)重要的失效形式，該現(xiàn)象普遍存在于高速撞擊、高速成型等涉及爆炸與沖擊的高速變形過(guò)程中，采用CFD計(jì)算方法可描述這種結(jié)構(gòu)大變形、剪切效應(yīng)明顯的爆炸失效[15-16]。而儲(chǔ)罐爆炸碎片的產(chǎn)生更傾向于小變形、極高應(yīng)變率的材料失效，采用TNT當(dāng)量法則可描述傳播速度更快速的爆破過(guò)程。材料的破壞不是瞬間完成的，而是一個(gè)損傷積累的過(guò)程。本文選用指數(shù)內(nèi)聚力模型對(duì)球罐材料的初始線彈性階段和損傷后的應(yīng)力軟化階段進(jìn)行描述。有研究表明[17],指數(shù)內(nèi)聚力模型能夠較好地描述材料破壞過(guò)程中的能量耗散行為。

球罐受爆炸沖擊波作用，當(dāng)應(yīng)力超過(guò)內(nèi)聚強(qiáng)度σmax，材料就發(fā)生了損傷。在損傷階段，材料基于指數(shù)內(nèi)聚力模型的失效準(zhǔn)則定義如下：

(21)

隨著載荷的持續(xù)作用，損傷不斷發(fā)展，當(dāng)D=1時(shí)，材料完全損傷，這時(shí)球罐發(fā)生破壞。

為保證球罐內(nèi)爆發(fā)生后罐體及產(chǎn)生碎片的質(zhì)量守恒，采用在罐體材料中引入粘結(jié)單元的方法對(duì)內(nèi)爆前后儲(chǔ)罐碎片狀態(tài)進(jìn)行模擬。通過(guò)在罐體單元之間插入一層厚度很小的粘結(jié)單元，將罐體單元之間的直接聯(lián)系轉(zhuǎn)換為通過(guò)粘結(jié)單元間接聯(lián)系，利用粘結(jié)單元的失效模擬罐體單元之間的分裂，爆炸碎片形成的過(guò)程如圖5所示。

圖5 儲(chǔ)罐內(nèi)爆碎片形成過(guò)程示意

在爆炸載荷作用下，粘結(jié)單元首先進(jìn)入彈性階段，應(yīng)力與變形呈線性關(guān)系。應(yīng)力不斷增大，超過(guò)內(nèi)聚強(qiáng)度σmax后，粘結(jié)單元進(jìn)入損傷階段，其承載能力不斷降低，直至損傷參數(shù)D=1時(shí)，損傷完全，粘結(jié)單元因失去承載能力而被刪除。少量的粘結(jié)單元被刪除后，宏觀上體現(xiàn)出來(lái)的就是球罐裂紋的產(chǎn)生，當(dāng)大量粘結(jié)單元被刪除且形成一個(gè)閉合的回路時(shí)，宏觀上體現(xiàn)出來(lái)的就是球罐碎片的產(chǎn)生。這種間接的方法保證了罐體材料的質(zhì)量守恒，更符合實(shí)際情況。

3 基于多相耦合的球罐內(nèi)爆破碎過(guò)程數(shù)值模擬分析

3.1 球罐內(nèi)爆數(shù)值模型的建立

以球形儲(chǔ)罐為例進(jìn)行分析，其結(jié)構(gòu)如圖6所示，其中球罐本體是球罐結(jié)構(gòu)的主體，它是球罐儲(chǔ)存物料承受物料工作壓力和液體靜壓力的構(gòu)件；球罐支柱是用于支承球罐本體重量和儲(chǔ)存物料重量的結(jié)構(gòu)部件，赤道正切柱式支座是使用最多的一種形式。本文模擬的球罐結(jié)構(gòu)參數(shù)見(jiàn)表1。

圖6 球罐結(jié)構(gòu)示意

表1 球罐結(jié)構(gòu)參數(shù)

采用Abaqus模擬平臺(tái)進(jìn)行建模以及計(jì)算，為了便于建模，將球罐簡(jiǎn)化為等壁厚的球殼，將支座對(duì)球罐的約束作用簡(jiǎn)化為罐壁上的約束條件。為了保證儲(chǔ)罐結(jié)構(gòu)響應(yīng)、碎片飛射都在流體域內(nèi)進(jìn)行，建立包含球罐外部輪廓的空氣域，罐外空氣域尺寸為20 000 mm×20 000 mm×19 000 mm。由于模型的對(duì)稱性，建立1/4球罐及空氣域模型，如圖7所示。

圖7 球罐內(nèi)爆幾何模型

考慮到球罐內(nèi)爆產(chǎn)生拋射碎片是個(gè)三維空間問(wèn)題，所有模型采用實(shí)體單元進(jìn)行離散，球罐用拉格朗日網(wǎng)格進(jìn)行描述，TNT炸藥與流體采用歐拉網(wǎng)格進(jìn)行描述。球罐的三維有限元模型如圖8所示，流體域的有限元模型如圖9所示。對(duì)于球罐結(jié)構(gòu)，在內(nèi)爆過(guò)程中，只考慮內(nèi)爆載荷和重力的影響，并將赤道正切柱式支座簡(jiǎn)化為支座與球罐接觸位置處的完全固定約束。對(duì)于球罐內(nèi)部流體域，在對(duì)稱面(XZ平面、YZ平面)施加對(duì)稱邊界條件，其余表面施加無(wú)反射邊界條件。儲(chǔ)罐內(nèi)氣體燃爆載荷采用TNT當(dāng)量法進(jìn)行模擬[18]，本文將TNT設(shè)置在球罐正中心處，并于計(jì)算開(kāi)始時(shí)引爆。

圖8 球罐有限元模型

圖9 流體域有限元模型

TNT炸藥采用JWL狀態(tài)方程來(lái)描述爆炸過(guò)程，其狀態(tài)方程表達(dá)式為：

(22)

式中,PT為TNT爆炸壓力,MPa；A,B,R1,R2,ω為TNT材料常數(shù)；E0為TNT初始比內(nèi)能,J/kg；η為無(wú)量綱量，η=ρ/ρ0；ρ為爆轟產(chǎn)物密度,kg/m3；ρ0為炸藥密度,kg/m3。

各參數(shù)具體數(shù)值見(jiàn)表2。

表2 TNT材料參數(shù)

在高速劇烈載荷下，球罐材料表現(xiàn)出應(yīng)變率強(qiáng)化效應(yīng)，故采用Johnson-Cook[19-20]強(qiáng)化模型來(lái)描述，其表達(dá)式為：

(23)

球罐材料采用Q345鋼，其材料參數(shù)如表3所示。

表3 球罐材料(Q345)參數(shù)

為了通過(guò)粘結(jié)單元的失效模擬球罐的破壞，粘結(jié)單元應(yīng)采用與球罐材料相同的力學(xué)性能。由于斷裂能、內(nèi)聚強(qiáng)度以及失效位移等參數(shù)難以獲得，因此參考了文獻(xiàn)[21]中鋼的部分材料參數(shù)，得到的粘結(jié)單元材料參數(shù)如表4所示。

表4 粘結(jié)單元材料參數(shù)

罐內(nèi)氣體與罐外空氣均可采用理想氣體狀態(tài)方程來(lái)描述：

Pg=(γ-1)ρgE1

(24)

式中,Pg為氣體壓力,MPa；γ為氣體絕熱系數(shù)；ρg為氣體密度,kg/m3；E1為氣體內(nèi)能,J/kg。

3.2 球罐內(nèi)爆數(shù)值模擬計(jì)算結(jié)果分析

對(duì)于容積為1 600 m3的LPG球罐，在空罐狀態(tài)時(shí)，球罐下部殘留少量的儲(chǔ)液，球罐上部存在著大量的石油氣和空氣，其中石油氣的主要成分為丙烷和丁烷。若忽略殘余儲(chǔ)液體積，在烷類氣體濃度為8%時(shí)，罐內(nèi)可燃?xì)怏wTNT當(dāng)量為528 kg。

3.2.1 球罐結(jié)構(gòu)應(yīng)力分布

當(dāng)儲(chǔ)罐材料應(yīng)力達(dá)到內(nèi)聚強(qiáng)度517 MPa后，罐體發(fā)生損傷并導(dǎo)致斷裂失效。為了清晰地觀測(cè)到球罐結(jié)構(gòu)斷裂過(guò)程，在爆炸發(fā)生后第5 ms(破壞前)、第6 ms(破壞后)以及第23 ms(出現(xiàn)明顯裂紋)時(shí)的粘結(jié)單元狀態(tài)以及球罐應(yīng)力分布如圖10所示。

圖10 不同時(shí)刻球罐粘結(jié)單元狀態(tài)及應(yīng)力分布

在5 ms時(shí)，球罐結(jié)構(gòu)完整且應(yīng)力低于內(nèi)聚強(qiáng)度；在6 ms時(shí)，應(yīng)力超過(guò)內(nèi)聚強(qiáng)度517 MPa，開(kāi)始有粘結(jié)單元失效被刪除，結(jié)合罐體的應(yīng)力分布可以看出，粘結(jié)單元失效的位置罐體材料應(yīng)力很小，球罐的這些位置已經(jīng)發(fā)生斷裂破壞，不過(guò)宏觀上不是很明顯；在23 ms時(shí)，罐體材料已經(jīng)出現(xiàn)較多且明顯的裂紋。在發(fā)生破壞后(第6 ms)，不斷有新的粘結(jié)單元失效被刪除，因此球罐的應(yīng)力穩(wěn)定在略高于內(nèi)聚強(qiáng)度的水平。

3.2.2 球罐結(jié)構(gòu)速度分布

在爆炸發(fā)生23 ms后，隨著球罐裂紋越來(lái)越多，第一個(gè)碎片形成。在爆炸載荷持續(xù)作用下，碎片數(shù)量越來(lái)越多，速度也越來(lái)越大。在82 ms時(shí)，已經(jīng)不再有新碎片產(chǎn)生，此時(shí)碎片數(shù)量最多且碎片速度達(dá)到最大值。第一個(gè)碎片形成(第38 ms)以及碎片數(shù)量最多、速度最大(第82 ms)時(shí)球罐結(jié)構(gòu)的速度分布如圖11所示。

圖11 不同時(shí)刻球罐碎片速度分布

3.2.3 碎片的形狀、質(zhì)量分布

經(jīng)過(guò)統(tǒng)計(jì)，82 ms時(shí)共產(chǎn)生了6種不同形狀的碎片，總數(shù)量為14個(gè)，每種形狀碎片的速度分布如圖12所示，碎片的初始速度和質(zhì)量如表5所示。

圖12 第82 ms時(shí)球罐碎片速度分布

通過(guò)表5中的數(shù)據(jù)可以看出，質(zhì)量在3 000～5 000 kg范圍內(nèi)的碎片有8個(gè)，占碎片總數(shù)的57%;在10 000～20 000 kg的大質(zhì)量碎片共6個(gè)，占總數(shù)的43%?？傮w來(lái)看，碎片的速度隨著質(zhì)量的增大而減小，但是碎片c、碎片d分別位于球罐正上方、正下方，受重力的影響較大，導(dǎo)致了位于正下方的碎片d速度更大。

表5 球罐碎片速度-質(zhì)量統(tǒng)計(jì)結(jié)果

4 基于不同計(jì)算方法的內(nèi)爆碎片形成與拋射軌跡對(duì)比分析

4.1 基于多相耦合的碎片拋射軌跡數(shù)值模擬分析

在82 ms時(shí)刻，二維平面坐標(biāo)系中碎片b的初始狀態(tài)如圖13所示，碎片b的拋射軌跡計(jì)算過(guò)程如下。

圖13 第82 ms時(shí)碎片b初始狀態(tài)

(1)通過(guò)數(shù)值模擬結(jié)果得到了碎片的初始位置坐標(biāo)為(6.36，15.45)，質(zhì)量為4 563 kg，初始速度v0=72.4 m/s，拋射角θ=47°。

(2)選取與初始速度方向垂直的投影面，將碎片b向投影面投影后，得到碎片的迎風(fēng)面積S=29.25 m2?？諝饷芏圈補(bǔ)取1.25 kg/m3，空氣阻力系數(shù)C取0.3[22]，通過(guò)式(8)得到比例系數(shù)k=5.48 kg/m。

(3)在碎片軌跡方程式(19)中，令Y=0，可求得碎片b的拋射距離為548 m，通過(guò)公式(16)(20)得到球罐整體坐標(biāo)系下碎片b的拋射高度H=157 m。

經(jīng)計(jì)算，文中6種形狀碎片的最終拋射高度、拋射距離如表6所示，各碎片可能的拋射軌跡如圖14所示。

表6 碎片信息統(tǒng)計(jì)

圖14 碎片可能拋射軌跡

通過(guò)表6和圖12可以看出，在球罐爆炸中，a,d,e,f四類碎片向下拋射(θ為負(fù)值)共有9個(gè)，占總數(shù)的64%；b,c兩類碎片向上拋射(θ為正值)共有5個(gè)，占總數(shù)的36%。在向下拋射的碎片中，拋射高度為碎片的初始高度，拋射距離最大為58 m。對(duì)于向上拋射的碎片，c類碎片豎直向上拋射，所以拋射高度最高為231 m，拋射距離為0；b類碎片拋射角度為47°，拋射距離最大為548 m，拋射高度較高,為157 m。對(duì)比各類碎片的信息可以發(fā)現(xiàn)，向下拋射的碎片拋射高度和拋射距離都比較?。幌蛏蠏伾涞乃槠瑨伾浣嵌圈仍酱?、拋射高度越高；向上拋射且拋射角度θ在45°附近的碎片,拋射距離最遠(yuǎn)。

4.2 基于蒙特卡洛理論求解的碎片軌跡計(jì)算

蒙特卡洛計(jì)算方法以概率和統(tǒng)計(jì)理論為基礎(chǔ)，將所求解的問(wèn)題同一定的概率模型相聯(lián)系，通過(guò)大量的模擬統(tǒng)計(jì)以獲得問(wèn)題的近似解。文獻(xiàn)[24]采用蒙特卡洛方法對(duì)1 600 m3球罐爆炸碎片的拋射軌跡進(jìn)行了計(jì)算，并通過(guò)與實(shí)際事故數(shù)據(jù)相比較，驗(yàn)證了蒙特卡洛方法的有效性。

4.2.1 碎片的能量與速度

在蒙特卡洛方法中，儲(chǔ)罐爆炸總能量通過(guò)式(25)計(jì)算，然后通過(guò)式(26)(27)得到碎片獲得的動(dòng)能Ek和碎片的拋射速度vp。

(25)

Ek=Efk

(26)

(27)

式中，p0為環(huán)境壓力，MPa；p1為爆炸壓力，MPa；γ為氣體絕熱指數(shù)；VS為儲(chǔ)罐體積，m3；fk為能量轉(zhuǎn)化效率；m為碎片質(zhì)量，kg。

根據(jù)文獻(xiàn)[22]中沸騰液體擴(kuò)展為蒸氣云爆炸(BLEVE)事故球罐的參數(shù)，其爆炸壓力為1.34 MPa，絕熱指數(shù)為1.4，爆炸壓力p1在區(qū)間[0.9p1，1.1p1]內(nèi)均勻分布時(shí)，通過(guò)式(25)計(jì)算得出爆炸總能量為2.61×109～3.32×109J。能量轉(zhuǎn)化效率fk在區(qū)間[0.2，0.5]內(nèi)右三角分布，按照f(shuō)k=0.2計(jì)算，碎片總動(dòng)能為5.22×108～6.64×108J，碎片的初始速度分布在84～95 m/s之間。

4.2.2 碎片的拋射軌跡

當(dāng)碎片數(shù)量、拋射角度等隨機(jī)參數(shù)已知后，采用蒙特卡洛方法對(duì)不確定參數(shù)進(jìn)行大量的隨機(jī)抽樣，并代入碎片拋射軌跡方程(17)中。球罐爆炸100次的碎片軌跡空間分布情況如圖15[24]所示，可以看出，球罐爆炸碎片的最大拋射高度為263 m，最大射程為600 m。

圖15 100次模擬碎片軌跡

4.3 不同計(jì)算方法碎片結(jié)果對(duì)比分析

為了驗(yàn)證本文多相耦合計(jì)算方法的正確性，同時(shí)突出多相耦合計(jì)算方法的優(yōu)勢(shì)，將多相耦合計(jì)算方法的結(jié)果與文獻(xiàn)[23]中蒙特卡洛計(jì)算方法的結(jié)果進(jìn)行對(duì)比，兩種方法中碎片數(shù)量、質(zhì)量、能量等碎片初始信息對(duì)比如表7所示。

表7 不同計(jì)算方法碎片初始信息對(duì)比

由表7可以看出，在碎片的數(shù)量、形狀和質(zhì)量方面，蒙特卡洛方法都基于一些基本假設(shè)，不能具體描述碎片的差異性，而多相耦合方法得到了碎片的明確數(shù)量、質(zhì)量以及具體形狀。在實(shí)際中，碎片的數(shù)量、形狀和質(zhì)量具有隨機(jī)性和特異性，而相比之下，多相耦合方法更有優(yōu)勢(shì)。文獻(xiàn)[23]忽略了球罐的體積，假設(shè)所有碎片從同一位置出發(fā)，而多相耦合計(jì)算結(jié)果可以考慮球罐的體積，進(jìn)而得到碎片不同的初始位置。

對(duì)于爆炸總能量，蒙特卡洛方法中為2.61×109～3.32×109J，由于爆炸壓力p1在一定范圍內(nèi)分布，導(dǎo)致得到的結(jié)果范圍較大;多相耦合方法得到的總能量為2.85×109J，在蒙特卡洛方法得到的結(jié)果范圍之內(nèi)。按照能量轉(zhuǎn)化效率為0.2計(jì)算，蒙特卡洛方法中碎片總動(dòng)能為5.22×108～6.64×108J;多相耦合中碎片總動(dòng)能為3.68×108J，能量轉(zhuǎn)化效率為0.129。相比之下，多相耦合計(jì)算方法的能量轉(zhuǎn)化效率比蒙特卡洛方法中的最小值偏小35.5%。有研究表明[24],能量轉(zhuǎn)化效率越低其概率越高，能量轉(zhuǎn)化效率在0.14以下占89.6%，顯然文獻(xiàn)[23]中總能量轉(zhuǎn)化為碎片動(dòng)能的效率fk取值偏大，多相耦合方法中的能量轉(zhuǎn)化效率更符合實(shí)際。

由于文獻(xiàn)[23]蒙特卡洛方法中碎片動(dòng)能偏大，導(dǎo)致碎片的初始速度偏大，直接影響了碎片的拋射高度和拋射距離，結(jié)果偏保守。與蒙特卡洛的結(jié)果相比，多相耦合中碎片拋射高度偏小12.2%，拋射距離偏小8.7%。

綜上所述，蒙特卡洛方法選用的不確定參數(shù)都在一定范圍內(nèi)分布，也就是包含了碎片拋射過(guò)程中的所有可能性，因此在模擬爆炸碎片的分布概率上更有優(yōu)勢(shì)，但同時(shí)也導(dǎo)致了拋射距離、拋射高度等范圍偏大，結(jié)果偏保守。而采用多相耦合計(jì)算方法在得到球罐動(dòng)態(tài)破碎過(guò)程的同時(shí)，還能準(zhǔn)確得到碎片的數(shù)量、速度、拋射距離等信息，所以在一個(gè)具體的球罐內(nèi)爆碎片拋射問(wèn)題上，多相耦合計(jì)算方法更精準(zhǔn)、更有優(yōu)勢(shì)。

5 結(jié)論

(1)以球形儲(chǔ)罐為研究對(duì)象，考慮罐內(nèi)氣體-球罐(碎片)-罐外空氣的多相耦合以及爆炸載荷的傳遞方式，采用CEL流固耦合算法建立球罐內(nèi)爆多相耦合有限元模型，基于指數(shù)內(nèi)聚力模型建立損傷失效準(zhǔn)則，將罐體單元之間的直接聯(lián)系轉(zhuǎn)換為通過(guò)粘結(jié)單元聯(lián)系，并通過(guò)粘結(jié)單元的失效刪除間接描述球罐材料的破壞，保證了球罐材料破碎前后的質(zhì)量守恒。

(2)基于多相耦合的計(jì)算方法得到了球罐不同時(shí)刻的結(jié)構(gòu)響應(yīng)以及碎片的具體數(shù)量、形狀、質(zhì)量和速度等?？紤]碎片不同初始位置以及空氣阻力的影響，通過(guò)理論計(jì)算方法得到碎片的拋射軌跡和最大拋射距離。對(duì)比蒙特卡洛方法，多相耦合計(jì)算方法中能量轉(zhuǎn)化效率比蒙特卡洛方法中的最小值偏小35.5%，碎片的最大拋射高度偏小12.2%，最大拋射距離偏小8.7%。

(3)在空罐狀態(tài)下對(duì)1 600 m3球罐的多相耦合內(nèi)爆模型進(jìn)行計(jì)算，在82 ms時(shí)，已經(jīng)不再有新碎片產(chǎn)生且碎片速度達(dá)到最大值，經(jīng)統(tǒng)計(jì)此時(shí)碎片的數(shù)量為14個(gè)；質(zhì)量在3 000～5 000 kg范圍內(nèi)的碎片占碎片總數(shù)的57%，在10 000～20 000 kg的大質(zhì)量碎片占總數(shù)的43%；爆炸總能量到碎片動(dòng)能的能量轉(zhuǎn)化效率為0.129；碎片速度在61.2～75.6 m/s之間分布；碎片最大拋射高度為231 m，最大拋射距離為548 m。