应用DSMC方法研究烟幕颗粒的撞击扩散特性

刘国生，关 华，宋东明，张祺祺

（南京理工大学 化工学院，南京 210094）

0 引 言

烟幕撞击喷撒是以气体为动力源，采用两个及以上喷嘴相对撞击的方式，进行烟幕颗粒分散的过程［1］，其重要特征是气固两相及固相颗粒间存在大量碰撞现象。在烟幕成型阶段，这种相互作用往往对扩散起着重要作用［2－3］。目前，烟幕扩散的数值模拟研究中，常用模型有高斯扩散模型和随机游走模型［3－5］。前者以K理论为基础，假定湍流交换系数为常数，后者计算中未考虑颗粒间的相互作用。因此，两模型在烟幕撞击扩散中的应用受到一定限制。

在流化床及撞击流等相关研究领域中，人们对颗粒间相互作用进行了深入研究，并发展了颗粒碰撞动力学理论。Savage［6］等将气体分子动理论进行推广，提出了颗粒相动力学理论，引入颗粒碰撞恢复系数的概念，并运用Boltzmann方程，对颗粒碰撞动力学方程进行了修正。Tsuji［7］等建立了硬球颗粒碰撞动力学模型，在计算中忽略颗粒间的碰撞而仅考虑颗粒与壁面碰撞，该假设仅在流动中固相体积比非常低时才成立。Lun［8］等考虑了颗粒碰撞的滑移与粘性，并将颗粒碰撞模型耦合到两相流模型中。欧阳洁［9］等针对气固两相流动特征，将颗粒相的运动过程处理为：颗粒与颗粒、颗粒与流体间的相互作用，建立了离散的颗粒运动、碰撞分解轨道模型。杜敏［10］等和王淑彦［11］等分别建立了离散颗粒运动－碰撞解耦模型，在模型中运用直接模拟Monte Carlo（DSMC）方法对颗粒间碰撞过程进行了模拟。

本文在气固两相流的基础上，建立烟幕的撞击扩散动力学模型，并采用DSMC方法对烟幕喷撒过程进行数值模拟。所得结果与实验吻合较好，验证了模拟方法的可靠性，为解决烟幕撞击扩散过程中的复杂流动问题提供一种新的研究思路。

1 烟幕撞击扩散动力学模型

1.1 气相控制方程

气相连续方程：

式中，vg、ρg分别为气体的速度与密度。

气相动量方程：

式中，τg为气相应力张量，μg为气相粘度，μlam，g为气相层流动力粘度，μt为气相湍流动力粘度，I为单位张量，Sp－g为两相间的作用力，其表达式为：

式中，N为网格中烟幕颗粒数量，Num为单个仿真颗粒所代表的真实颗粒数，si为颗粒所在网格的球表面积。

1.2 烟幕颗粒的运动方程

烟幕颗粒在外力作用下的运动满足牛顿第二定律。本文计算中考虑了重力、气体曳力及Saffman升力的作用［12－13］。烟幕颗粒运动方程可表示为：

式中，FD为气体曳力，FS为Saffman升力，二者均为合力的作用，计算颗粒分速度时需进行分解，表达式分别为：

式中，μ为空气粘度系数；dp为球形烟幕颗粒粒径；vp为烟幕颗粒速度。

1.3 烟幕颗粒的碰撞方程

等粒径的同种烟幕颗粒碰撞前后速度变化满足以下方程：

式中，mi、mj分别为颗粒i与j的质量，vi，0、vj，0与vi，1、vj，1分别为颗粒i与j碰撞前后的质心速度；vij为颗粒i与j的相对速度。J为碰撞时颗粒j作用于颗粒i的冲量。

式中，n为烟幕颗粒i与j之间的质心距离矢量；e为颗粒非弹性碰撞恢复系数。

2 模型的应用

2.1 烟幕颗粒碰撞过程的实现

将烟幕颗粒运动分解为颗粒自由运动与颗粒间碰撞过程。1）在足够小的时间步长Δt内求解颗粒运动方程，计算其速度并确定位移，此时不考虑颗粒碰撞；2）采用DSMC方法计算颗粒间的碰撞，若Δt内颗粒发生碰撞，则确定碰撞后速度，此时，保持颗粒位置不变，将其速度更新为碰撞后速度。

对于取样颗粒i，采用修正Nanbu算法［14］在所处网格内搜索其碰撞对象。首先，生成0～1之间的均匀随机数Rrand，得到颗粒j（j＝int（Rrand×N）＋1，其中，int（）表示取整数函数。则颗粒j即为颗粒i的碰撞对象，碰撞几率Pij可表示为：

式中，n为当地真实颗粒数；N为网格内的仿真颗粒数。若满足：Rrand＞（j／N－pij），则认为发生碰撞，否则不发生。

2.2 计算网格及模拟参数

2.2.1 计算物理及网格模型

图1为计算的物理模型及网格模型，为显示网格内部结构，图1（b）中采用球面对网格进行了剖视。

图1 算例的物理模型与网格模型Fig.1 Physical and mesh model of the example

如图1（a）所示，坐标系位于圆柱体底面中心，两喷嘴成120°布置于圆柱底部，喷嘴轴线的延长线交于z轴，圆柱体高0.5m，直径0.4m，喷嘴间距l＝0.1m，喷嘴直径D＝0.0011m，圆柱体各端面以敞开面作为出口。

根据文献［15］，网格特征尺寸（Δs）应满足：λ≤Δs≤λ，λ为颗粒间的平均自由程，其值在计算前估算得到。所以，在网格剖分时，对颗粒浓度较大的喷嘴附近区域进行加密处理，其余部分采用五棱柱网格离散，网格总数约为22.3万。

2.2.2 模拟参数及方法

模拟参数设置如下：气体、烟幕颗粒入口速度为ug，0＝vg，0＝25m／s，颗粒直径10μm、密度1.2g／cm3，气相时间步长 Δtg＝1×10－4s，颗粒相时间步长Δtp，max＝5×10－5s，实际颗粒数与取样颗粒数之比Num≤3000，单喷嘴烟幕质量流率Wp＝0.03kg／s。

2.3 计算方法

应用Fluent软件求解气相流场，烟幕颗粒碰撞计算由编程实现，并采用自定义函数（UDF）形式嵌入到Fluent求解器中。计算时，通过交替循环的方式实现两相间的耦合计算，并由自定义函数进行相间动量与能量的传递。计算过程中，颗粒经出口离开计算域即停止对其跟踪，当计算域内的跟踪颗粒数基本不变时，认为计算收敛。

3 计算结果及分析

3.1 气相流场分析

图2为y＝0截面的气相速度矢量分布。由图2可见，两股气体射流1喷出后，运动至区域5（称为碰撞区）时发生碰撞，由于两者的挤压、碰撞，致其速度发生改变，射流1转向形成折射流2，沿z轴正方向扩张，同时在碰撞区两侧形成了沿z轴反方向的低速回旋流4。随高度增加，折射流2不断扩张，最终演变成上升流3和扩散流6，流场合二为一，呈一定角度不断向外扩张。

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

3.2 烟幕颗粒的运动行为

图3为典型时刻的烟幕颗粒分布图，图中颗粒颜色代表其竖直方向的速度大小。文中所有烟幕运动状态图中带颜色部分皆为烟幕颗粒，喷嘴未画出。

图3 烟幕颗粒的运动行为Fig.3 Motion behavior of smoke particle

由图3可见，烟幕射流首先沿喷嘴轴线方向运动，当射流碰撞后即发生转向，改沿竖直方向向上运动，并不断向外扩散。在0.004s前（见图3b），颗粒速度未发生明显变化。这是因为颗粒间的碰撞是由它们之间的速度差引起的，此时相邻颗粒的速度相当，因而发生碰撞的概率较小，即使发生碰撞，其速度改变也较小。在0.010s（见图3c），射流转向后速度有所减弱。分析认为，类似于撞击流，颗粒在惯性作用下渗入另一喷嘴形成的流场中，之后因流场的阻力作用减速，速度降为零后又被该流场加速，再次回到原喷嘴流场中，如此反复振荡若干次，颗粒水平方向速度分量逐渐消失，最终随气流一道离开碰撞区向上运动；此外，若射流间的颗粒发生了碰撞，则更加快了这一进程［2］。在0.030s时（见图3d），烟幕颗粒逐渐上升，集中分布于计算域中部。

图4为0.03s时左侧单个喷嘴的烟幕颗粒分布图，图中颗粒颜色为其在计算域内的停留时间。由图可见，仅有少部分颗粒渗入到x＝0面右侧的流场中。统计表明，碰撞区的浓度最大值可达25kg／m3，所以在碰撞区及其上方的颗粒必将发生剧烈碰撞，导致其向x＝0面另一侧流场的扩散运动受阻。

图4 tp＝0.030s时单个喷嘴的烟幕颗粒分布Fig.4 Smoke particle distribution of single－nozzle at 0.030s

图5为烟幕撞击喷撒试验得到的颗粒运动状态图，相邻图幅的时间间隔为4ms。

对比图3与图5可知，烟幕颗粒的运动状态很相似，计算结果与试验吻合较好。所以，本文建立的烟幕撞击扩散动力学模型对烟幕颗粒运动与碰撞过程的计算结果是合理的。

图5 烟幕撞击扩散的试验图片Fig.5 Experimental pictures of smoke particle distribution byimpact diffusion

3.3 烟幕速度场分析

图6为烟幕颗粒的速度矢量分布。由图6（a）可见，烟幕射流相遇后，由于气动流场的剧变及颗粒间的强烈碰撞作用，烟幕转向后以6m／s～9m／s的速度发散，其速度约减小65%～75%，形成了上升颗粒流2和扩张颗粒流1、3（见图6c），同时在气体涡流的作用下（见图2），形成了回旋颗粒流4。此外，在烟幕射流的扩散过程中，不同方向（如y＝0与x＝0截面）的颗粒速度矢量分布的差异较小，烟幕以近似球顶锥体的形态向外发散。由于射流碰撞后形成的新流场各向异性较小，以致颗粒在各个方向的受力较均衡，因而扩散状态也很相似。

图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方向的浓度梯度较大。当z＝0.15m时，x方向的浓度最大值约达0.5kg／m3，烟幕集中分布于x轴两侧，由于此高度正位于喷嘴交汇处附近，颗粒大多在此处相遇，所以浓度较大；z＝0.25m时，浓度最大值约为0.1kg／m3；z＝0.35m、z＝0.45m时，浓度最大值约降至0.02kg／m3量级，随高度增加烟幕浓度迅速降低。分析认为，距离碰撞区越远时颗粒的碰撞效应越弱，在重力、气体曳力及随机力的作用下，颗粒的随机扩散运动趋势逐渐增强［16］。根据Fick扩散定律，颗粒将向着浓度较低的地方不断迁移，使烟幕浓度梯度逐步减小［17］。所以，随高度增加，x方向的浓度分布均匀性也逐渐提高。

图7 z轴不同高度处x方向的烟幕浓度分布Fig.7 Smoke concentration distribution in x direction at different height

4 结 论

本文建立了烟幕撞击扩散动力学模型，并采用DSMC方法对烟幕撞击喷撒过程进行了数值模拟，经分析得到：

（1）烟幕撞击扩散的计算结果与实验吻合较好，DSMC方法能有效模拟烟幕撞击喷撒过程。

（2）射流相遇后，由于气流挤压及颗粒间的碰撞作用，烟幕转向向上运动，颗粒速度约减小65%～75%；在扩张气流的作用下，烟幕以近似球顶锥体的形态逐步向外扩散。

（3）烟幕浓度极大值集中在射流及碰撞区，在此以外，由于碰撞效应的减弱，烟幕颗粒随机扩散运动的趋势逐步增强，使其浓度梯度随高度增加迅速减小，空间分布也更趋均匀。

（4）本研究为烟幕撞击扩散运动的数值模拟研究提供了一种新的思路；但作为尝试，计算中未考虑颗粒的凝并、碰并等效应，有待于进一步发展。

［1］PAN G P，YANG S.Principle of pyrotechnics［M］.Beijing：Beijing Institute of Technology Press，1997.（in Chinese）潘功配，杨硕.烟火学［M］.北京：北京理工大学出版社，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.撞击流反应器－原理和应用［M］.伍沅，译.北京：化学工业出版社，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）姚禄玖，高钧麟，肖凯涛，等.烟幕理论与测试技术［M］.北京：国防工业出版社，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］.安全与环境工程学报，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）欧阳洁，李静海.模拟气固两相流动非均匀结构的颗粒运动分解轨道模型［J］.中国科学（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方法在大规模气固两相撞击流中的应用［J］.化工学报，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）王淑彦，赵云华，姜健，等.数值模拟颗粒旋转对流化床内气固两相流动特性的影响［J］.工程热物理学报，2007，02：84－86.

［12］KONG L.Two phase fluid mechanics［M］.Beijing：Higher Education Press，2003.（in Chinese）孔珑.两相流体力学［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）宋仔标，刘代志，王红霞.重力作用下的颗粒扩散模型及数值模拟［J］.空气动力学学报，2010，28（2）：209－212.

［17］LU Z Y.Aerosol science and technology［M］.Beijing：Atomic Energy Press，1999.（in Chinese）卢正永.气溶胶科学引论［M］.北京：原子能出版社，1999.