爆炸煙云中粒子運(yùn)動(dòng)的數(shù)值模擬＊

袁 濤，羅永鋒，尚愛國，羅昆升，高洪泉

（1.第二炮兵工程大學(xué)，陜西 西安 710025；2.第二炮兵裝備研究院第二研究所，北京 100085）

爆炸裝置爆炸會(huì)產(chǎn)生煙云［1］，煙云的主要成分是被氣溶膠化的顆粒物。這些微小的顆粒物可以長(zhǎng)時(shí)間懸浮在空氣中，并隨著環(huán)境而擴(kuò)散或遷移。在這個(gè)過程中，它們極易被動(dòng)植物吸附在表面或吸入體內(nèi)，根據(jù)顆粒物的種類不同而造成不同程度的危害。研究這些顆粒物運(yùn)動(dòng)和沉積特性是評(píng)估爆炸產(chǎn)物對(duì)環(huán)境影響的理論基礎(chǔ)。實(shí)地實(shí)驗(yàn)從事顆粒沉積的研究非常困難，這主要是因?yàn)闊熢瞥跏紲囟雀摺⒀舆t時(shí)間短、變化性強(qiáng)，并且對(duì)于其中的物理過程和化學(xué)過程知道的非常少［2－3］。爆炸產(chǎn)物中顆粒污染物的排放依據(jù)顆粒尺寸可以分為2大類：第1類是粒徑小于50μm的顆粒，煙云在大氣運(yùn)動(dòng)時(shí)它們將隨氣相排放物而分散；第2類是粒徑大于50μm的較大顆粒，它們由于重力作用而沉降于地面上。本文中擬通過數(shù)值計(jì)算，分析爆炸煙云近場(chǎng)擴(kuò)散的粒子運(yùn)動(dòng)規(guī)律，選取煙云氣溶膠中具有代表性的粒徑大小1、10、50、100μm粒子來研究這2類顆粒的運(yùn)動(dòng)規(guī)律，以期得到的量化結(jié)果可用于實(shí)驗(yàn)設(shè)計(jì)，評(píng)估爆炸煙云對(duì)空氣、水體和土壤等的影響。

1 爆炸源項(xiàng)計(jì)算

圖1 t＝1ms時(shí)，源項(xiàng)的壓力分布Fig.1 Pressure distribution of source at t＝1ms

本文的研究對(duì)象是10kg TNT爆炸裝置在開闊無風(fēng)環(huán)境下爆炸的煙云運(yùn)動(dòng)問題。根據(jù)源項(xiàng)表征模型［4］，將高能炸藥的爆炸過程劃分為爆炸、浮力氣團(tuán)上升、大氣散布3個(gè)理化過程。第1階段爆炸過程是一個(gè)極其復(fù)雜的理化過程，整個(gè)特征時(shí)間尺度為幾皮秒到幾微秒［5］。就本文而言，第1階段的爆炸過程是第2階段浮力氣團(tuán)上升過程的源項(xiàng)，所以必須要得到源項(xiàng)的物理特征參數(shù)。

為了簡(jiǎn)化源項(xiàng)計(jì)算，本文中做了如下設(shè)定：

（1）環(huán)境溫度為300K，氣壓為標(biāo)準(zhǔn)大氣壓，空氣為理想氣體；

（2）炸藥為球狀裝藥，放置于地面中央，裝藥密度ρe＝1.63g／cm3，半徑r＝11.36cm；

（3）第1階段時(shí)間尺度為1ms［6］，即源項(xiàng)的初始時(shí)刻為1ms。

Fluent軟件不能直接計(jì)算爆炸裝置爆炸后的物理參數(shù)，利用Autodyn軟件計(jì)算爆炸第1階段后壓力、溫度等物理參數(shù)，為源項(xiàng)第2階段在Fluent軟件中的計(jì)算提供初始參數(shù)［7］。接下來在Autodyn軟件中建立幾何模型，輸入初始參數(shù)條件，計(jì)算得到1ms時(shí)源項(xiàng)壓力分布如圖1所示。

圖2 爆炸源項(xiàng)雙層模型Fig.2 Adouble－layer model for explosion source

從圖1中可以看出，1ms時(shí)，爆炸產(chǎn)物的外緣為擊波區(qū)域，壓力p高，形成了明顯的分層狀態(tài)，如果在源項(xiàng)設(shè)定中使用簡(jiǎn)單的球狀或橢球狀單層結(jié)構(gòu)是不合適的?；贏utodyn的計(jì)算結(jié)果，本文中建立了如圖2所示的爆炸源項(xiàng)雙層模型：內(nèi)層為短軸a1＝1.466m，長(zhǎng)軸b1＝1.806m，壓力 p1＝18kPa的負(fù)壓橢球；外層為短軸a2＝1.8 m，長(zhǎng)軸b2＝2.5m，壓力p2＝530kPa的超壓殼狀結(jié)構(gòu)。

在源項(xiàng)幾何模型中填充粒子要根據(jù)爆炸后的物質(zhì)密度ρ分布。圖3表明：密度ρ分布基本服從壓力的分布，其規(guī)律是內(nèi)部稀疏、外層緊實(shí)，并且外層粒子密度是內(nèi)層的4～6倍。結(jié)合源項(xiàng)雙層模型，本文中設(shè)置外層粒子密度為內(nèi)層的5倍。

圖4是1ms時(shí)粒子速度v與距離爆心距離s關(guān)系圖。爆炸后產(chǎn)生的能量將給予粒子初始運(yùn)動(dòng)速度，速度分布的規(guī)律即：擊波區(qū)域速度最高，擊波區(qū)域到爆心為線性衰減。

至此，本文源項(xiàng)所需的物理參數(shù)已全部得到。

圖3 t＝1ms時(shí)，源項(xiàng)密度分布Fig.3 Density distribution of source at t＝1ms

圖4 t＝1ms時(shí)，源項(xiàng)速度曲線Fig.4 Velocity curve of source at t＝1ms

2 幾何模型建立及控制方程設(shè)置

若要建立源項(xiàng)雙層模型，必須在Fluent前處理軟件GAMBIT中實(shí)現(xiàn)。首先在單位上設(shè)置以m為標(biāo)準(zhǔn)單位，按照計(jì)算好的源項(xiàng)尺寸繪制。這里有幾點(diǎn)需要注意，其一是在公共邊的設(shè)置上要注意屬性上設(shè)置為Connect部件，否則將會(huì)導(dǎo)致網(wǎng)格錯(cuò)誤；其二是外環(huán)境的總尺寸問題，原則上來說開闊無風(fēng)環(huán)境應(yīng)該是無邊界最合適。但是在有限元計(jì)算時(shí)，這個(gè)是不允許的，必須要設(shè)置一個(gè)邊界。所以尺寸的選定也很關(guān)鍵：太大會(huì)極大地增加網(wǎng)格數(shù)，增加計(jì)算負(fù)荷，降低計(jì)算精度；太小又會(huì)影響源項(xiàng)壓力的釋放，導(dǎo)致煙云計(jì)算不準(zhǔn)確。本文中根據(jù)壓力衰減的規(guī)律和煙云的尺寸預(yù)先估算出合適的模型大小，再經(jīng)過幾次粗網(wǎng)格快速計(jì)算幾個(gè)尺寸，比較壓力最遠(yuǎn)傳輸范圍，最后取其中的最優(yōu)值。

整個(gè)區(qū)域計(jì)算的空間設(shè)置在一個(gè)60m×60m的區(qū)域，源項(xiàng)在底部的正中央。按照1∶1的尺寸建立雙層模型，根據(jù)源項(xiàng)的物理參數(shù)控制源項(xiàng)的網(wǎng)格設(shè)置，輸出網(wǎng)格，導(dǎo)入到Fluent軟件中去。

以2ddp模式啟動(dòng)Fluent軟件，在File－Read－Case中選擇已經(jīng)創(chuàng)建好的網(wǎng)格文件（.msh）。使用Grid－Check檢查網(wǎng)格是否完整，在Grid－Scale中檢查是否以m為標(biāo)準(zhǔn)單位，檢查無誤后進(jìn)入基本控制方程的設(shè)置。

在Define－Models－Solver中選擇Time－Unsteady（非定常）和 Unsteady Formulation－2nd Order Implicit（2階方程），其他保持默認(rèn)。勾選 Define－Models－Energy中的 Energy Equation（能量守恒），在Define－Models－Viscous Model中選擇k－ε湍流方程。最后勾選Define－Models－Discrete Phase Model中的Interaction with Continuous Phase（粒子與空氣相互作用）和Brownian Motion（布朗運(yùn)動(dòng)）等選項(xiàng)。接下來進(jìn)行環(huán)境的參數(shù)設(shè)置。

在Define－Operating Conditions中定義壓力為標(biāo)準(zhǔn)大氣壓，重力加速度為9.8m／s2，空氣密度為1.225kg／m3。在Define－Boundary Conditions中分別定義外層邊界為標(biāo)準(zhǔn)大氣壓恒壓邊界，且為粒子的逃逸界面；底邊為地面，且為粒子反射界面。在Define－Materials中分別定義爆炸煙云中物質(zhì)的密度、傳熱系數(shù)等。

源項(xiàng)粒子的初始設(shè)定是通過Define－Injects來實(shí)現(xiàn)的，根據(jù)本文的要求，分別創(chuàng)建不同粒徑的源。在Injection Type中選擇Surface，在Release from Surface中選擇源項(xiàng)的2個(gè)面，在Diameter中輸入粒徑。如此反復(fù)，創(chuàng)建出1、10、50、100μm的粒子源。

最后對(duì)求解器進(jìn)行設(shè)置，在Solve－Controls－Solution Controls中選擇PISO求解方程，收斂因子保持默認(rèn)。使用Solve－Initialize中的Solution Initialize和Patch填充雙層壓力，在Solve－Monitor中設(shè)置誤差限值，并選定監(jiān)測(cè)誤差，使用Solve－Animate和Solve－Particle工具記錄粒子運(yùn)動(dòng)軌跡和圖像。在Solve－Iterate中使用0.01s為步長(zhǎng)進(jìn)行模擬。

3 粒子運(yùn)動(dòng)模擬結(jié)果

煙云中粒子的擴(kuò)散與沉積特性如表1所示，不同粒徑的粒子在運(yùn)動(dòng)時(shí)呈現(xiàn)出不同的規(guī)律。粒徑越小的粒子更容易受到浮力渦環(huán)的作用，在雙渦環(huán)的作用下翻滾著向上運(yùn)動(dòng)。構(gòu)成煙云上端的云團(tuán)主要是粒徑小于50μm的粒子，并且在向上翻滾拋灑的過程中，云團(tuán)的雙渦環(huán)由內(nèi)向外是粒徑不斷增大的。

從壓力演變的圖像中可以發(fā)現(xiàn)，爆炸后初期渦環(huán)上升的速度較快。隨著時(shí)間的推移，渦環(huán)上升的速度變慢，逐漸成為一個(gè)動(dòng)態(tài)平衡的狀態(tài)。渦環(huán)中的壓力釋放速度較快，且渦環(huán)會(huì)不斷膨脹，直至壓力與大氣平衡，失去束縛云團(tuán)中粒子的能力。渦環(huán)消失時(shí)，煙云即會(huì)進(jìn)入SCM模型的第3個(gè)階段，隨著大氣而進(jìn)行沉積和遷移［8］。

比較粒徑為50和100μm的粒子的運(yùn)動(dòng)規(guī)律可以發(fā)現(xiàn)，粒徑大的粒子構(gòu)成了煙云的莖部，并且粒徑越大，粒子的分布就越靠下。粒徑為50μm的粒子幾乎是均勻地構(gòu)成了煙云的莖部桿狀圖像，上部的煙云則是分布在外圍。而粒徑為100μm的粒子在重力的作用下，沉積得更快，絕大部分粒子已經(jīng)沉降到了底部的區(qū)域。與渦環(huán)的浮力作用相比，質(zhì)量大的粒子更易在重力的作用下沉積。

美國勞倫斯·利弗莫爾實(shí)驗(yàn)室利用實(shí)驗(yàn)數(shù)據(jù)和經(jīng)驗(yàn)公式開發(fā)了經(jīng)典的計(jì)算常規(guī)爆炸煙云上升的Hotspot程序［9］。這個(gè)程序系統(tǒng)地描述了煙云上升后的表征尺度，并得出煙云上升后直徑與高度的比值約為0.4。

提取本文計(jì)算的結(jié)果，煙云在5和10s的表征尺度分別為：Ht＝5＝41.4m，Dt＝5＝13.2m；Ht＝10＝49.8m，Dt＝10＝16.2m，即比值分別為0.32和0.33。這與 Hotspot軟件得到的結(jié)果相比偏小，主要原因是本文模擬的對(duì)象是被氣溶膠化的粒徑范圍，而實(shí)際中會(huì)產(chǎn)生粒徑大于100μm的粒子。粒子質(zhì)量越大，在上升的過程中就容易被拋灑得更遠(yuǎn)，使煙云的直徑增大。所以，本文的比值與Hotspot軟件基本吻合。

表1 不同時(shí)間的粒子分布與壓力分布Table1 Distribution of particles and pressure at different times

4 結(jié) 論

（1）爆炸裝置爆炸后會(huì)迅速產(chǎn)生雙渦環(huán)結(jié)構(gòu)，在負(fù)壓渦環(huán)的作用下粒徑小的顆粒（50μm以下）極易受渦環(huán)的影響而翻滾上升運(yùn)動(dòng)。（2）在煙云運(yùn)動(dòng)的拋灑過程中，由于慣性的作用，質(zhì)量大的顆粒分布在煙云的邊緣，質(zhì)量小的顆粒則分布在頂部煙團(tuán)中。（3）煙云上升的過程中，構(gòu)成頂層云團(tuán)的粒子是粒徑在50μm以下的粒子，粒徑在50μm以上的大顆粒在重力的作用下沉積到地面或形成桿狀莖部。

本文研究結(jié)果能夠揭示爆炸后煙云中粒子運(yùn)動(dòng)的規(guī)律，可為爆炸后產(chǎn)物對(duì)環(huán)境的污染和人類的危害的研究提供了理論基礎(chǔ)，同時(shí)，該結(jié)果也能為實(shí)驗(yàn)的設(shè)計(jì)提供參考。

［1］鄭波，陳力，丁雁生，等.溫壓炸藥爆炸拋撒的運(yùn)動(dòng)規(guī)律［J］.爆炸與沖擊，2008，28（5）：433－437.Zheng Bo，Chen Li，Ding Yan－sheng，et al.Dispersal process of explosion production of thermobaric explosive［J］.Explosion and Shock Waves，2008，28（5）：433－437.

［2］Morton B R，Taylor G，Turner J S.Turbulent gravitational convection from maintained and instantaneous source［J］.Proceedings of the Royal Society A：Mathematical，Physical and Engingeering Sciences，1956，234（1196）：1－23.

［3］Makhviladze G M，Roberts J P，Yakush S E.Modelling of atmospheric pollution by explosions［J］.Environmental Software，1995，10（2）：117－127.

［4］Brown R C，Kolb C E，Conant J A，et al.A predictive capability for the source terms of residual energetic materials from burning and／or detonation activities［R］.ARI－RR－1384，2004.

［5］Zhang X J，Ghoniem A F.A computational model for the rise and dispersion of wind－blown，buoyancy－driven plumes：Ⅰ：Neutrally stratified atmosphere［J］.Atmospheric Environment：Part A：General Topics，1993，27（15）：2295－2311.

［6］劉庚冉，曹保榆，周凱元，等.FAE野外空中爆炸效應(yīng)實(shí)驗(yàn)研究及其一次起爆機(jī)理探討［J］.防化研究，2007（1）：7－11.Liu Geng－ran，Cao Bao－yun，Zhou Kai－yuan，et al.Experimental investigation of FAE blast field effect in the free air and the discussion of the single ignition mechanism［J］.Anti－chemical Warfare Research，2007（1）：7－11.

［7］周保順，張立恒，王少龍，等.TNT炸藥爆炸沖擊波的數(shù)值模擬與實(shí)驗(yàn)研究［J］.彈箭與制導(dǎo)學(xué)報(bào)，2010，30（3）：88－90.Zhou Bao－shun，Zhang Li－h(huán)eng，Wang Shao－long，et al.Numerical simulation and experimental research on TNT explosion shock wave［J］.Journal of Projectiles，Rockets，Missiles and Guidance，2010，30（3）：88－90.

［8］鄭毅，應(yīng)純同，李磊.穩(wěn)定大氣條件下爆炸煙云的數(shù)值模擬研究［J］.力學(xué)與實(shí)踐，2009，31（3）：31－34.Zheng Yi，Ying Chun－tong，Li Lei.Numberical study on explosion clouds in stably stratified atmosphere［J］.Mechanics in Engineering，2009，31（3）：31－34.

［9］HOTSPOT 2.06help manual［M］.Lawence Livermore National Laboratory，2005.