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

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

0 引 言

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

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

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

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

1.1 氣相控制方程

氣相連續(xù)方程：

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

氣相動量方程：

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

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

1.2 煙幕顆粒的運動方程

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

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

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

1.3 煙幕顆粒的碰撞方程

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

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

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

2 模型的應(yīng)用

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

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

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

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

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

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

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

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

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

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

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

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

2.3 計算方法

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

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

3.1 氣相流場分析

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

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

3.2 煙幕顆粒的運動行為

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

圖3 煙幕顆粒的運動行為Fig.3 Motion behavior of smoke particle

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

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

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

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

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

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

3.3 煙幕速度場分析

圖6為煙幕顆粒的速度矢量分布。由圖6（a）可見，煙幕射流相遇后，由于氣動流場的劇變及顆粒間的強烈碰撞作用，煙幕轉(zhuǎn)向后以6m／s～9m／s的速度發(fā)散，其速度約減小65%～75%，形成了上升顆粒流2和擴張顆粒流1、3（見圖6c），同時在氣體渦流的作用下（見圖2），形成了回旋顆粒流4。此外，在煙幕射流的擴散過程中，不同方向（如y＝0與x＝0截面）的顆粒速度矢量分布的差異較小，煙幕以近似球頂錐體的形態(tài)向外發(fā)散。由于射流碰撞后形成的新流場各向異性較小，以致顆粒在各個方向的受力較均衡，因而擴散狀態(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時，x方向的濃度最大值約達0.5kg／m3，煙幕集中分布于x軸兩側(cè)，由于此高度正位于噴嘴交匯處附近，顆粒大多在此處相遇，所以濃度較大；z＝0.25m時，濃度最大值約為0.1kg／m3；z＝0.35m、z＝0.45m時，濃度最大值約降至0.02kg／m3量級，隨高度增加煙幕濃度迅速降低。分析認為，距離碰撞區(qū)越遠時顆粒的碰撞效應(yīng)越弱，在重力、氣體曳力及隨機力的作用下，顆粒的隨機擴散運動趨勢逐漸增強［16］。根據(jù)Fick擴散定律，顆粒將向著濃度較低的地方不斷遷移，使煙幕濃度梯度逐步減?。?7］。所以，隨高度增加，x方向的濃度分布均勻性也逐漸提高。

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

4 結(jié) 論

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

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

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

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

（4）本研究為煙幕撞擊擴散運動的數(shù)值模擬研究提供了一種新的思路；但作為嘗試，計算中未考慮顆粒的凝并、碰并等效應(yīng)，有待于進一步發(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）姚祿玖，高鈞麟，肖凱濤，等.煙幕理論與測試技術(shù)［M］.北京：國防工業(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）金陶勝，楊錦，韓素芹，等.隨機游動模型在大氣污染擴散分析中的應(yīng)用探討［J］.安全與環(huán)境工程學(xué)報，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）歐陽潔，李靜海.模擬氣固兩相流動非均勻結(jié)構(gòu)的顆粒運動分解軌道模型［J］.中國科學(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é)報，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)對流化床內(nèi)氣固兩相流動特性的影響［J］.工程熱物理學(xué)報，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)，劉代志，王紅霞.重力作用下的顆粒擴散模型及數(shù)值模擬［J］.空氣動力學(xué)學(xué)報，2010，28（2）：209－212.

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