應(yīng)用DSMC方法研究煙幕顆粒的撞擊擴(kuò)散特性

劉國(guó)生，關(guān) 華，宋東明，張祺祺

（南京理工大學(xué) 化工學(xué)院，南京 210094）

0 引 言

煙幕撞擊噴撒是以氣體為動(dòng)力源，采用兩個(gè)及以上噴嘴相對(duì)撞擊的方式，進(jìn)行煙幕顆粒分散的過程［1］，其重要特征是氣固兩相及固相顆粒間存在大量碰撞現(xiàn)象。在煙幕成型階段，這種相互作用往往對(duì)擴(kuò)散起著重要作用［2－3］。目前，煙幕擴(kuò)散的數(shù)值模擬研究中，常用模型有高斯擴(kuò)散模型和隨機(jī)游走模型［3－5］。前者以K理論為基礎(chǔ)，假定湍流交換系數(shù)為常數(shù)，后者計(jì)算中未考慮顆粒間的相互作用。因此，兩模型在煙幕撞擊擴(kuò)散中的應(yīng)用受到一定限制。

在流化床及撞擊流等相關(guān)研究領(lǐng)域中，人們對(duì)顆粒間相互作用進(jìn)行了深入研究，并發(fā)展了顆粒碰撞動(dòng)力學(xué)理論。Savage［6］等將氣體分子動(dòng)理論進(jìn)行推廣，提出了顆粒相動(dòng)力學(xué)理論，引入顆粒碰撞恢復(fù)系數(shù)的概念，并運(yùn)用Boltzmann方程，對(duì)顆粒碰撞動(dòng)力學(xué)方程進(jìn)行了修正。Tsuji［7］等建立了硬球顆粒碰撞動(dòng)力學(xué)模型，在計(jì)算中忽略顆粒間的碰撞而僅考慮顆粒與壁面碰撞，該假設(shè)僅在流動(dòng)中固相體積比非常低時(shí)才成立。Lun［8］等考慮了顆粒碰撞的滑移與粘性，并將顆粒碰撞模型耦合到兩相流模型中。歐陽(yáng)潔［9］等針對(duì)氣固兩相流動(dòng)特征，將顆粒相的運(yùn)動(dòng)過程處理為：顆粒與顆粒、顆粒與流體間的相互作用，建立了離散的顆粒運(yùn)動(dòng)、碰撞分解軌道模型。杜敏［10］等和王淑彥［11］等分別建立了離散顆粒運(yùn)動(dòng)－碰撞解耦模型，在模型中運(yùn)用直接模擬Monte Carlo（DSMC）方法對(duì)顆粒間碰撞過程進(jìn)行了模擬。

本文在氣固兩相流的基礎(chǔ)上，建立煙幕的撞擊擴(kuò)散動(dòng)力學(xué)模型，并采用DSMC方法對(duì)煙幕噴撒過程進(jìn)行數(shù)值模擬。所得結(jié)果與實(shí)驗(yàn)吻合較好，驗(yàn)證了模擬方法的可靠性，為解決煙幕撞擊擴(kuò)散過程中的復(fù)雜流動(dòng)問題提供一種新的研究思路。

1 煙幕撞擊擴(kuò)散動(dòng)力學(xué)模型

1.1 氣相控制方程

氣相連續(xù)方程：

式中，vg、ρg分別為氣體的速度與密度。

氣相動(dòng)量方程：

式中，τg為氣相應(yīng)力張量，μg為氣相粘度，μlam，g為氣相層流動(dòng)力粘度，μt為氣相湍流動(dòng)力粘度，I為單位張量，Sp－g為兩相間的作用力，其表達(dá)式為：

式中，N為網(wǎng)格中煙幕顆粒數(shù)量，Num為單個(gè)仿真顆粒所代表的真實(shí)顆粒數(shù)，si為顆粒所在網(wǎng)格的球表面積。

1.2 煙幕顆粒的運(yùn)動(dòng)方程

煙幕顆粒在外力作用下的運(yùn)動(dòng)滿足牛頓第二定律。本文計(jì)算中考慮了重力、氣體曳力及Saffman升力的作用［12－13］。煙幕顆粒運(yùn)動(dòng)方程可表示為：

式中，F(xiàn)D為氣體曳力，F(xiàn)S為Saffman升力，二者均為合力的作用，計(jì)算顆粒分速度時(shí)需進(jìn)行分解，表達(dá)式分別為：

式中，μ為空氣粘度系數(shù)；dp為球形煙幕顆粒粒徑；vp為煙幕顆粒速度。

1.3 煙幕顆粒的碰撞方程

等粒徑的同種煙幕顆粒碰撞前后速度變化滿足以下方程：

式中，mi、mj分別為顆粒i與j的質(zhì)量，vi，0、vj，0與vi，1、vj，1分別為顆粒i與j碰撞前后的質(zhì)心速度；vij為顆粒i與j的相對(duì)速度。J為碰撞時(shí)顆粒j作用于顆粒i的沖量。

式中，n為煙幕顆粒i與j之間的質(zhì)心距離矢量；e為顆粒非彈性碰撞恢復(fù)系數(shù)。

2 模型的應(yīng)用

2.1 煙幕顆粒碰撞過程的實(shí)現(xiàn)

將煙幕顆粒運(yùn)動(dòng)分解為顆粒自由運(yùn)動(dòng)與顆粒間碰撞過程。1）在足夠小的時(shí)間步長(zhǎng)Δt內(nèi)求解顆粒運(yùn)動(dòng)方程，計(jì)算其速度并確定位移，此時(shí)不考慮顆粒碰撞；2）采用DSMC方法計(jì)算顆粒間的碰撞，若Δt內(nèi)顆粒發(fā)生碰撞，則確定碰撞后速度，此時(shí)，保持顆粒位置不變，將其速度更新為碰撞后速度。

對(duì)于取樣顆粒i，采用修正Nanbu算法［14］在所處網(wǎng)格內(nèi)搜索其碰撞對(duì)象。首先，生成0～1之間的均勻隨機(jī)數(shù)Rrand，得到顆粒j（j＝int（Rrand×N）＋1，其中，int（）表示取整數(shù)函數(shù)。則顆粒j即為顆粒i的碰撞對(duì)象，碰撞幾率Pij可表示為：

式中，n為當(dāng)?shù)卣鎸?shí)顆粒數(shù)；N為網(wǎng)格內(nèi)的仿真顆粒數(shù)。若滿足：Rrand＞（j／N－pij），則認(rèn)為發(fā)生碰撞，否則不發(fā)生。

2.2 計(jì)算網(wǎng)格及模擬參數(shù)

2.2.1 計(jì)算物理及網(wǎng)格模型

圖1為計(jì)算的物理模型及網(wǎng)格模型，為顯示網(wǎng)格內(nèi)部結(jié)構(gòu)，圖1（b）中采用球面對(duì)網(wǎng)格進(jìn)行了剖視。

圖1 算例的物理模型與網(wǎng)格模型Fig.1 Physical and mesh model of the example

如圖1（a）所示，坐標(biāo)系位于圓柱體底面中心，兩噴嘴成120°布置于圓柱底部，噴嘴軸線的延長(zhǎng)線交于z軸，圓柱體高0.5m，直徑0.4m，噴嘴間距l(xiāng)＝0.1m，噴嘴直徑D＝0.0011m，圓柱體各端面以敞開面作為出口。

根據(jù)文獻(xiàn)［15］，網(wǎng)格特征尺寸（Δs）應(yīng)滿足：λ≤Δs≤λ，λ為顆粒間的平均自由程，其值在計(jì)算前估算得到。所以，在網(wǎng)格剖分時(shí)，對(duì)顆粒濃度較大的噴嘴附近區(qū)域進(jìn)行加密處理，其余部分采用五棱柱網(wǎng)格離散，網(wǎng)格總數(shù)約為22.3萬(wàn)。

2.2.2 模擬參數(shù)及方法

模擬參數(shù)設(shè)置如下：氣體、煙幕顆粒入口速度為ug，0＝vg，0＝25m／s，顆粒直徑10μm、密度1.2g／cm3，氣相時(shí)間步長(zhǎng) Δtg＝1×10－4s，顆粒相時(shí)間步長(zhǎng)Δtp，max＝5×10－5s，實(shí)際顆粒數(shù)與取樣顆粒數(shù)之比Num≤3000，單噴嘴煙幕質(zhì)量流率Wp＝0.03kg／s。

2.3 計(jì)算方法

應(yīng)用Fluent軟件求解氣相流場(chǎng)，煙幕顆粒碰撞計(jì)算由編程實(shí)現(xiàn)，并采用自定義函數(shù)（UDF）形式嵌入到Fluent求解器中。計(jì)算時(shí)，通過交替循環(huán)的方式實(shí)現(xiàn)兩相間的耦合計(jì)算，并由自定義函數(shù)進(jìn)行相間動(dòng)量與能量的傳遞。計(jì)算過程中，顆粒經(jīng)出口離開計(jì)算域即停止對(duì)其跟蹤，當(dāng)計(jì)算域內(nèi)的跟蹤顆粒數(shù)基本不變時(shí)，認(rèn)為計(jì)算收斂。

3 計(jì)算結(jié)果及分析

3.1 氣相流場(chǎng)分析

圖2為y＝0截面的氣相速度矢量分布。由圖2可見，兩股氣體射流1噴出后，運(yùn)動(dòng)至區(qū)域5（稱為碰撞區(qū)）時(shí)發(fā)生碰撞，由于兩者的擠壓、碰撞，致其速度發(fā)生改變，射流1轉(zhuǎn)向形成折射流2，沿z軸正方向擴(kuò)張，同時(shí)在碰撞區(qū)兩側(cè)形成了沿z軸反方向的低速回旋流4。隨高度增加，折射流2不斷擴(kuò)張，最終演變成上升流3和擴(kuò)散流6，流場(chǎng)合二為一，呈一定角度不斷向外擴(kuò)張。

圖2 y＝0截面的氣相速度矢量分布Fig.2 Velocity vector distribution of gas－phase in y＝0section

3.2 煙幕顆粒的運(yùn)動(dòng)行為

圖3為典型時(shí)刻的煙幕顆粒分布圖，圖中顆粒顏色代表其豎直方向的速度大小。文中所有煙幕運(yùn)動(dòng)狀態(tài)圖中帶顏色部分皆為煙幕顆粒，噴嘴未畫出。

圖3 煙幕顆粒的運(yùn)動(dòng)行為Fig.3 Motion behavior of smoke particle

由圖3可見，煙幕射流首先沿噴嘴軸線方向運(yùn)動(dòng)，當(dāng)射流碰撞后即發(fā)生轉(zhuǎn)向，改沿豎直方向向上運(yùn)動(dòng)，并不斷向外擴(kuò)散。在0.004s前（見圖3b），顆粒速度未發(fā)生明顯變化。這是因?yàn)轭w粒間的碰撞是由它們之間的速度差引起的，此時(shí)相鄰顆粒的速度相當(dāng)，因而發(fā)生碰撞的概率較小，即使發(fā)生碰撞，其速度改變也較小。在0.010s（見圖3c），射流轉(zhuǎn)向后速度有所減弱。分析認(rèn)為，類似于撞擊流，顆粒在慣性作用下滲入另一噴嘴形成的流場(chǎng)中，之后因流場(chǎng)的阻力作用減速，速度降為零后又被該流場(chǎng)加速，再次回到原噴嘴流場(chǎng)中，如此反復(fù)振蕩若干次，顆粒水平方向速度分量逐漸消失，最終隨氣流一道離開碰撞區(qū)向上運(yùn)動(dòng)；此外，若射流間的顆粒發(fā)生了碰撞，則更加快了這一進(jìn)程［2］。在0.030s時(shí)（見圖3d），煙幕顆粒逐漸上升，集中分布于計(jì)算域中部。

圖4為0.03s時(shí)左側(cè)單個(gè)噴嘴的煙幕顆粒分布圖，圖中顆粒顏色為其在計(jì)算域內(nèi)的停留時(shí)間。由圖可見，僅有少部分顆粒滲入到x＝0面右側(cè)的流場(chǎng)中。統(tǒng)計(jì)表明，碰撞區(qū)的濃度最大值可達(dá)25kg／m3，所以在碰撞區(qū)及其上方的顆粒必將發(fā)生劇烈碰撞，導(dǎo)致其向x＝0面另一側(cè)流場(chǎng)的擴(kuò)散運(yùn)動(dòng)受阻。

圖4 tp＝0.030s時(shí)單個(gè)噴嘴的煙幕顆粒分布Fig.4 Smoke particle distribution of single－nozzle at 0.030s

圖5為煙幕撞擊噴撒試驗(yàn)得到的顆粒運(yùn)動(dòng)狀態(tài)圖，相鄰圖幅的時(shí)間間隔為4ms。

對(duì)比圖3與圖5可知，煙幕顆粒的運(yùn)動(dòng)狀態(tài)很相似，計(jì)算結(jié)果與試驗(yàn)吻合較好。所以，本文建立的煙幕撞擊擴(kuò)散動(dòng)力學(xué)模型對(duì)煙幕顆粒運(yùn)動(dòng)與碰撞過程的計(jì)算結(jié)果是合理的。

圖5 煙幕撞擊擴(kuò)散的試驗(yàn)圖片F(xiàn)ig.5 Experimental pictures of smoke particle distribution byimpact diffusion

3.3 煙幕速度場(chǎng)分析

圖6為煙幕顆粒的速度矢量分布。由圖6（a）可見，煙幕射流相遇后，由于氣動(dòng)流場(chǎng)的劇變及顆粒間的強(qiáng)烈碰撞作用，煙幕轉(zhuǎn)向后以6m／s～9m／s的速度發(fā)散，其速度約減小65%～75%，形成了上升顆粒流2和擴(kuò)張顆粒流1、3（見圖6c），同時(shí)在氣體渦流的作用下（見圖2），形成了回旋顆粒流4。此外，在煙幕射流的擴(kuò)散過程中，不同方向（如y＝0與x＝0截面）的顆粒速度矢量分布的差異較小，煙幕以近似球頂錐體的形態(tài)向外發(fā)散。由于射流碰撞后形成的新流場(chǎng)各向異性較小，以致顆粒在各個(gè)方向的受力較均衡，因而擴(kuò)散狀態(tài)也很相似。

圖6 y＝0與x＝0截面煙幕顆粒速度矢量分布Fig.6 Velocity vector distribution of smoke particle in y＝0and x＝0section

3.4 煙幕濃度分布

圖7為z軸不同高度處x方向的煙幕濃度分布。由圖7可見，不同高度處x方向的濃度梯度較大。當(dāng)z＝0.15m時(shí)，x方向的濃度最大值約達(dá)0.5kg／m3，煙幕集中分布于x軸兩側(cè)，由于此高度正位于噴嘴交匯處附近，顆粒大多在此處相遇，所以濃度較大；z＝0.25m時(shí)，濃度最大值約為0.1kg／m3；z＝0.35m、z＝0.45m時(shí)，濃度最大值約降至0.02kg／m3量級(jí)，隨高度增加煙幕濃度迅速降低。分析認(rèn)為，距離碰撞區(qū)越遠(yuǎn)時(shí)顆粒的碰撞效應(yīng)越弱，在重力、氣體曳力及隨機(jī)力的作用下，顆粒的隨機(jī)擴(kuò)散運(yùn)動(dòng)趨勢(shì)逐漸增強(qiáng)［16］。根據(jù)Fick擴(kuò)散定律，顆粒將向著濃度較低的地方不斷遷移，使煙幕濃度梯度逐步減?。?7］。所以，隨高度增加，x方向的濃度分布均勻性也逐漸提高。

圖7 z軸不同高度處x方向的煙幕濃度分布Fig.7 Smoke concentration distribution in x direction at different height

4 結(jié) 論

本文建立了煙幕撞擊擴(kuò)散動(dòng)力學(xué)模型，并采用DSMC方法對(duì)煙幕撞擊噴撒過程進(jìn)行了數(shù)值模擬，經(jīng)分析得到：

（1）煙幕撞擊擴(kuò)散的計(jì)算結(jié)果與實(shí)驗(yàn)吻合較好，DSMC方法能有效模擬煙幕撞擊噴撒過程。

（2）射流相遇后，由于氣流擠壓及顆粒間的碰撞作用，煙幕轉(zhuǎn)向向上運(yùn)動(dòng)，顆粒速度約減小65%～75%；在擴(kuò)張氣流的作用下，煙幕以近似球頂錐體的形態(tài)逐步向外擴(kuò)散。

（3）煙幕濃度極大值集中在射流及碰撞區(qū)，在此以外，由于碰撞效應(yīng)的減弱，煙幕顆粒隨機(jī)擴(kuò)散運(yùn)動(dòng)的趨勢(shì)逐步增強(qiáng)，使其濃度梯度隨高度增加迅速減小，空間分布也更趨均勻。

（4）本研究為煙幕撞擊擴(kuò)散運(yùn)動(dòng)的數(shù)值模擬研究提供了一種新的思路；但作為嘗試，計(jì)算中未考慮顆粒的凝并、碰并等效應(yīng)，有待于進(jìn)一步發(fā)展。

［1］PAN G P，YANG S.Principle of pyrotechnics［M］.Beijing：Beijing Institute of Technology Press，1997.（in Chinese）潘功配，楊碩.煙火學(xué)［M］.北京：北京理工大學(xué)出版社，1997.

［2］TAMIR A.Impinging－stream reactors：principle and application［M］.WU Y，edit.Beijing：Chemical Industry Press，1996：9－10.（in Chinese）TAMIR A.撞擊流反應(yīng)器－原理和應(yīng)用［M］.伍沅，譯.北京：化學(xué)工業(yè)出版社，1996：9－10.

［3］YAO J L，GAO J L，XIAO K T，et al.Smoke theory and testing technology［M］.Beijing：National Defense Industry Press，1997.（in Chinese）姚祿玖，高鈞麟，肖凱濤，等.煙幕理論與測(cè)試技術(shù)［M］.北京：國(guó)防工業(yè)出版社，1997.

［4］ZHU C G，PAN G P，WU X Y，et al.Measurement and account method for the available standing［J］.Propellants，Ex－plosives，Pyrotechnics，2006，31（3）：226－229.

［5］JIN T S，YANG J，HAN S Q，et al.Practical discussion on the random floating model with the atmospheric pollutant dispersion［J］.JournalofSafetyandEnvironment，2008，8（1）：73－75.（in Chinese）金陶勝，楊錦，韓素芹，等.隨機(jī)游動(dòng)模型在大氣污染擴(kuò)散分析中的應(yīng)用探討［J］.安全與環(huán)境工程學(xué)報(bào)，2008，8（1）：73－75.

［6］SAVAGE S B，JEFFERY D J.The stress tensor in a granular flow at high shear rates［J］.JournalofMathematicalFluid Mechanics，1981，110：255－272.

［7］TSUJI Y，MORIKAWA Y，TANAKA T.Numerical simulation of gas－solid two－phase flow in a two－dimensional horizontal channel［J］.InternationalJournalofMultiphaseFlow，1987，13：671－702.

［8］LUN C K K.Numerical simulation of dilute turbulent gas－solid flows［J］.InternationalJournalofMultiphaseFlow，2000，26：1707－1736.

［9］OUYANG J，LI J H.The simulation of decomposition orbit and motion model of particle of gas－solid two－phase flow in Non uniform structure［J］.ChineseScience（BSeries），1999，29（1）：29－39.（in Chinese）歐陽(yáng)潔，李靜海.模擬氣固兩相流動(dòng)非均勻結(jié)構(gòu)的顆粒運(yùn)動(dòng)分解軌道模型［J］.中國(guó)科學(xué)（B輯），1999，29（1）：29－39.

［10］DU M，HAO Y L，LIU X D.Application of DSMC method in large－scale gas－solid two－phase impinging streams［J］.Journal ofchemicalindustryandEngineering，2009，60（8）：1950－1958.（in Chinese）杜敏，郝英力，劉向東.DSMC方法在大規(guī)模氣固兩相撞擊流中的應(yīng)用［J］.化工學(xué)報(bào)，2009，60（8）：1950－1958.

［11］WANG S Y，ZHAO Y H，JAING J，et al.Numerical simulation of effect of particle rotation on gas and solid flow behavior in fluidized beds［J］.JournalofEngineeringThermoPhysics，2007，02：84－86.（in Chinese）王淑彥，趙云華，姜健，等.數(shù)值模擬顆粒旋轉(zhuǎn)對(duì)流化床內(nèi)氣固兩相流動(dòng)特性的影響［J］.工程熱物理學(xué)報(bào)，2007，02：84－86.

［12］KONG L.Two phase fluid mechanics［M］.Beijing：Higher Education Press，2003.（in Chinese）孔瓏.兩相流體力學(xué)［M］.北京：高等教育出版社，2003.

［13］SAFFMAN P G.The lift on a small sphere in a slow shear flow［J］.JournalofMathematicalFluidMechanics，1965，22：385－400.

［14］ILLNER R，NEUNZERT H.On simulation methods for the Boltzmann equation［J］.TransportTheoryandStatistical Physics，1987，16（2／3）：141－154.

［15］BIRD G A.Molecular gas dynamics and the direct simulation of gas flows［M］.Oxford：Clarendon Press，1994.

［16］SONG Z B，LIU D Z，WANG H X.Mathematical model and numerical simulation of particles dispersion with gravity［J］.ACTAAerodynamicaSinica，2010，28（2）：209－212.（in Chinese）宋仔標(biāo)，劉代志，王紅霞.重力作用下的顆粒擴(kuò)散模型及數(shù)值模擬［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2010，28（2）：209－212.

［17］LU Z Y.Aerosol science and technology［M］.Beijing：Atomic Energy Press，1999.（in Chinese）盧正永.氣溶膠科學(xué)引論［M］.北京：原子能出版社，1999.