风洞管道中细微颗粒物浓度分布规律的数值模拟

Numerical Simulation of the Concentration of Fine Particles Distribution in Wind Tunnel Ducts

摘要: 本文利用计算流体动力学(CFD)研究风洞管道中气流驱动的细微颗粒物的扩散分布规律。基于欧拉–拉格朗日方法建立数值模型，使用随机轨道模型来考虑湍流扩散，求解颗粒物的扩散轨迹。分析气流风速和入口颗粒物质量流率变化对颗粒物扩散的影响，研究风洞管道内流场分布规律和颗粒物浓度场分布规律。结果表明，气流通过弯管之后的产生的二次流是影响颗粒物分布均匀性的主要原因。并且发现只改变管道内风速或者颗粒入口质量流率的情况下，颗粒物浓度随着风速增加而降低，随颗粒物质量流率增加而增加。研究结果有助于加深对颗粒物随气流运动规律及机理的研究，为颗粒物排放监控提供更合适的指导策略。

Abstract: In this study, computational fluid dynamics (CFD) was used to study the diffusion of particles driven by airflow in wind tunnel duct. A numerical model is established based on the Euler-Lagrange method, and discrete random walk model is used to consider the turbulence diffusion, and the trajectory of particles is solved. The influence of wind speed and mass flow rate of particles at the inlet on the diffusion was analyzed, and the distribution law of flow field and particles concentration field in wind tunnel ducts was studied. The results show that the secondary flow after the flow through the elbow is the main reason that affects the distribution uniformity of particles. It was also found that the concentration of particles decreased with the increase of wind speed and increased with the increase of mass flow rate of particles when only the wind speed in the ducts or the mass flow rate of particles at the inlet was changed. The research results are helpful to deepen the understanding on the movement law and mechanism of particles, and provide more appropriate guiding strategies for particles emission monitoring.

1. 引言

如今巨大的能源消耗造成污染物的大量排放，严重威胁着人们的身体健康。2013年，世界卫生组织(WHO)将细颗粒物(PM2.5)归为第一类致癌物质，人们对颗粒污染物PM2.5的关注不断增加 [1]。PM2.5超标是造成雾霾的主要原因，因此许多国家制定了相关标准来管制PM2.5的排放问题。随着环保标准的日益严格，环保局对各工厂烟气排放控制的要求也越来越高，因此监控颗粒物的排放有重要意义。

烟道中颗粒物的运动，从流动方式上来归类为气固两相流，目前国内外的学者做了许多与管道中气固两相流相关的研究。比如：Kliafas和Holt [2] 使用激光多普勒测速仪测量了90˚弯管方形截面内的气固两相流，并获得了平均流向速度和径向速度。Tsuji等人 [3] 和Huber、Sommerfeled [4] 等人发表了从相位多普勒风速测量(PDA)获得的关于粒子速度和粒子尺寸分布的空间发展的信息。Akilli [5] 等人使用光纤探头在水平管道的不同位置进行了颗粒浓度和颗粒速度的测量。宋勇 [6] 利用计算流体力学软件FLUENT对CLT/A型旋风分离器进行了数值模拟，验证了优化的可行性。杨瑞昌等人 [7] 提出了用于模拟多维湍流气固两相流的改进的扩散模型。李剑锋 [8] 利用Fluent对风机内部的气固两相流进行了数值模拟，与实验结果进行的比较表明模拟结果是可信的。Sa Z等人 [9] 利用Fluent软件模拟了硐室爆破后20分钟内通风采石场的粉尘浓度变化和不同风速下的粉尘轨迹，对减少采石场粉尘危害具有重要意义。Zhao M等人 [10] 通过CFD-DEM的方法模拟了流化床中湿颗粒的团聚和沉积现象。Yang D等人 [11] 针对大于5毫米的大尺寸颗粒，采用基于欧拉–拉格朗日方法和离散相模型的耦合方法进行了模拟研究，预测了大颗粒煤的水平气力输送特性。Bilal F S等人 [12] 研究了45˚和90˚弯管中颗粒物的侵蚀现象。

多数学者使用实验或数值模拟的方法研究了气固两相流动的运动规律，在颗粒物运动方面取得了大量进展，但很少学者研究管道中颗粒物浓度分布均匀性和稳定性，因此本文研究了颗粒物在风洞管道中的运动规律，分析内部流场和浓度场分布，为颗粒物排放监测提供指导。

2. 数学模型

2.1. 气相模型

本次模拟采用欧拉-拉格朗日方法 [13]。气体被视为连续相，用纳维尔-斯托克斯方程求解，而颗粒物被视为离散相，用牛顿第二定律求解。包括两个步骤：连续相流场模拟和颗粒跟踪。

由于风洞管道中流动雷诺数很高，处于紊流状态，因此使用基于雷诺平均方程的RNG k-ε [14] 湍流方程求解流场。

当气流处于稳定状态时，空气被认为是不可压缩的，连续性方程和动量方程表示为：

$\frac{\partial \left(\rho u\right)}{\partial x}+\frac{\partial \left(\rho v\right)}{\partial y}+\frac{\partial \left(\rho w\right)}{\partial z}=0$ (1)

$\frac{\partial }{\partial t}\left(\rho u_i\right)+\frac{\partial }{\partial x_j}\left(\rho u_i{u}_j\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\left(\mu +{\mu }_t\right)\left[\frac{\partial u_j}{\partial x_i}+\frac{\partial u_i}{\partial x_j}\right]\right]$ (2)

对于RNG k-ε方程，k方程为：

$\frac{\partial \left(\rho k{u}_i\right)}{\partial x_i}=\frac{\partial }{\partial x_j}\left({\alpha }_k{\mu }_{eff}\frac{\partial k}{\partial x_j}\right)+G_k+G_b-\rho \epsilon$ (3)

其中，k代表湍动能，m²/s²， ${\alpha }_k$ 代表湍动能的有效普朗特数的倒数，取1； ${\mu }_{eff}$ 为黏性系数；G_k，G_b分别为由平均速度梯度和浮力引起的湍动能；ε为湍动能耗散率，m²/s³。

方程为：

$\frac{\partial }{\partial x_i}\left(\rho \epsilon u_i\right)=\frac{\partial }{\partial x_j}\left({\alpha }_{\epsilon }{\mu }_{eff}\frac{\partial \epsilon }{\partial x_j}\right)+\frac{C_{1\epsilon }\epsilon }{k}\left(G_k+C_{3\epsilon }{G}_b\right)-C_{2\epsilon }\rho \frac{{\epsilon }^2}{k}-R_{\epsilon }$ (4)

${\mu }_{eff}={\mu }_t+\mu$ (5)

${\mu }_t=\rho C_{\mu }\frac{k^2}{\epsilon }$ (6)

其中C_1ε，C_2ε，C_3ε为经验常数，C_1ε = 1.44，C_2ε = 1.92，C_3ε = 0.09；μ_t，μ为湍流和层流的粘性系数；α_ε为耗散率的有效普朗特数的倒数，取0.769。

2.2. 颗粒相模型

颗粒物被认为是离散相，颗粒物的求解由牛顿第二定律得到 [15]。在笛卡尔坐标系下表示为：

$u_p=\frac{\text{d}x}{\text{d}t}$ (7)

$m\frac{\text{d}{u}_p}{\text{d}t}=F_D+\frac{mg\left({\rho }_p-{\rho }_f\right)}{{\rho }_p}+\sum F$ (8)

式中，u_p为颗粒的速度，m/s；u_f为气流速度，F_D为颗粒所受的曳力，N；g为重力加速度，m/s²；ρ_p为颗粒密度，kg/m³；ρ_f为流体密度，kg/m³； $\sum F$ 为颗粒受到的其他作用力的合力，因为本文使用的颗粒密度远大于空气，所以考虑曳力和重力，而忽略其他次要力 [16]。

$F_D=\frac{18\mu }{{\rho }_P{d}_p^2}\frac{C_{D}R{e}_P}{24}$ (9)

$R{e}_p=\frac{{\rho }_f{d}_p\left|u_p-u_f\right|}{\mu }$ (10)

$C_D=\{\begin{array}{l}\frac{\text{24}}{R{e}_p}\left(R{e}_p<1\right)\\ \frac{24}{R{e}_p}\left(1+\frac{1}{6}R{e}_p^{2/3}\right)\left(1\le R{e}_p<{10}^3\right)\\ 0.44\left(R{e}_p\ge {10}^3\right)\end{array}$ (11)

式中，μ为流体的动力粘度，Pa·s；Re_p为颗粒的雷诺数；C_D为阻力系数。

3. 数值模拟方法

3.1. 物理模型及网格划分

本次模拟基于上海北分风洞实验室的实际布置方式，管道分为矩形入口扩散段、过渡段、弯管段1、直圆管测量段、弯管段2和排出段，其中矩形管道尺寸为0.4 m × 0.56 m，弯管角度为90˚，曲率为R/D = 2，圆形管道尺寸为D = 0.56 m，将第二个弯管对测量段的影响考虑在内，所以在solidworks中建立何模型，在ICEM软件 [17] 中采用结构网格方式划分，对矩形管道向圆形管道过渡的地方以及两个弯管和壁面附近进行加密处理。图1(a)为风洞管道的整体网格示意图，图1(b)为矩形管段截面和圆管截面网格示意图。具体尺寸见图2。测量段法兰位置具体安排见图1(a)和表1。其中 $X^{*}=\left(X-0.84\right)/0.56$ 代表法兰距离弯管段1出口的距离。

(a)(b)

Figure 1. The whole grid diagram of the wind tunnel duct

图1. 风洞管道的整体网格示意图

Figure 2. Duct pressure drop and flange ⑨ average pressure change with the mesh number

图2. 管道压降和法兰⑨平均压力随网格数量的变化

Table 1. The specific location of the flange

表1. 法兰的具体位置

3.2. 网格无关性验证

为了在确保仿真模拟的准确性的前提下，提高计算效率。在网格质量符合要求的情况下划分了232,224、309,504、864,020、1,382,776和2,292,569个网格进行网格无关性验证。根据图2，发现在网格数量大于1382776时，管道整体压降和法兰⑨监测的平均压力变化均小于1%，表示计算结果趋于稳定，所以选取数量为1,382,776的网格进行计算。

3.3. 边界条件和算法

使用速度入口边界条件，壁面采用无滑移，出口定义为自由流出。颗粒在入口和出口设置为逃逸，壁面条件设置为反弹。压力–速度耦合求解采用Simple算法，压力，动量，湍动能，湍动能耗散率均采用二阶迎风差分格式。先计算稳态流场至收敛，再加入颗粒物继续计算。

3.4. 模型验证

为了验证模拟模型的准确性，将入口风速为8.76 m/s，颗粒质量流率为13 mg/s的条件下的模拟结果和已有的实验数据进行对比。图3为直管测量段管道中心轴线上颗粒物浓度值的数值模拟结果和实验结果的对比。可以发现，在直管测量段的中心轴线上，无论是数值模拟还是实验结果，颗粒物浓度随着管长都是先增加后下降，再趋于平稳，发展趋势高度吻合。数值模拟结果的颗粒物浓度曲线平稳时候是处于震荡的，这主要是因为拉格朗日方法的限制，颗粒物被认为是离散相，颗粒物与管壁的每一次碰撞都会造成计算结果的颗粒物浓度变化。图4为法兰⑥处沿管道Z轴方向(Z*为在管道径向测点距离管道内壁的距离)颗粒物浓度和带误差棒(误差范围为30%)的实验结果对比图。发现模拟结果在误差范围内，说明模拟方法是可行的。

Figure 3. Model validation: Particle concentration on the central axis of the straight pipe measurement section

图3. 模型验证：直管测量段中心轴线上的颗粒物浓度

Figure 4. Model validation: Particle concentration in the Z axis direction on the flange ⑥ section

图4. 模拟验证：法兰⑥截面上Z轴方向的颗粒物浓度分布

4. 仿真结果及分析

4.1. 流场分布

直管测量段是此风洞用于标定的工作段，因此本文主要关注此段的速度和浓度分布。图5为入口风速为8.76 m/s时的速度云图和速度矢量图。气流从风洞方形入口进入，形成壁面速度低，管道中心速度大的湍流充分发展状态，随后进入对气流影响最大的弯管部分。由于正压力梯度，管道内壁附近的气流被加速，而管道外壁附近的气流由于负压力梯度而减速，在弯曲的中点45˚附近后，这种行为开始逆转。并且在弯管出口处气流由于惯性力和离心力的共同作用向管道外壁面运动，靠近内壁附近形成气流低速区。在弯管出口形成了由两个反向旋转的漩涡，这两个漩涡被称为迪恩涡 [18]。随着管道的延长，漩涡面积才能减小并逐渐消失。

Figure 5. ZX section and YZ section velocity cloud diagram and velocity vector diagram of pipeline at different positions

图5. ZX截面和不同位置管道YZ截面速度云图和速度矢量图

选取截面€ (弯管段1的出口)和法兰①、④、⑦、⑨为目标，图6和图7显示了沿管道径向的速度和压力分布图，可以看到截面处沿管道径向气流速度和压力波动十分明显，这主要是因为此处二次流强度大，两个反向旋转的漩涡扰乱了流场，法兰①处速度和压力波动虽然减小，但还未稳定下来，法兰④正好处于向稳定流动状态发展的过渡阶段，法兰⑦和法兰⑨速度，压力变化基本消失，流场已经稳定下来。

Figure 6. Velocity distribution along the radial direction of the pipe at different positions

图6. 不同位置沿管道径向的速度分布

Figure 7. Pressure distribution along the pipe radial direction at different positions

图7. 不同位置沿管道径向的压力分布

4.2. 浓度场分布

图8为颗粒物浓度分布云图，颗粒物进入弯管后，在离心力，惯性力和气流曳力的共同作用下，颗粒物聚集在弯管出口外壁附近，而在内壁的低速区，几乎看不到颗粒物的存在。颗粒物的聚集现象在离开弯管在之后就逐渐消失，这主要是因为二次流漩涡的存在扰动流场，颗粒物跟随漩涡运动，沿管道截面半径呈对称分布，漩涡核心区域颗粒物难以进入，所以形成两个低浓度区。结合图5的速度分布云图和图8的浓度分布云图看出，在法兰①~④，漩涡逐渐消失，漩涡核心处的低浓度区域面积逐渐减小，颗粒物浓度分布均匀性也显现出来。从法兰⑤~⑨可以看到，漩涡消失后，颗粒物均匀的分布在管道中。因此提出了一个评估因子Φ来判断整体颗粒物浓度均匀性的好坏，定义如下：

$\Phi =\frac{\left|C_1-C_2\right|}{C_1}\times 100\%$ (12)

Figure 8. Concentration distribution cloud map of ZX section and different sections of duct

图8. 管道ZX截面和不同处管道截面浓度分布云图

其中C₁代表截面的颗粒物面积加权平均浓度，kg/m³；C₂代表截面上line-z和line-y上得到的颗粒物平均浓度，kg/m³。图9为管道截面上line-z和line-y的位置分布关系，图10为评估因子Φ的变化曲线，法兰①~⑨所处截面Φ值是逐渐下降，然后趋于平稳，距离弯管出口越远，Φ值越小，表明浓度分布逐渐均匀，法兰④之后，Φ都低于5%，说明整体分布已经达到初步均匀，从图8的浓度分布云图的也能看出这一点。

图11为法兰①管道截面上管道截面上line-z和line-y的浓度分布图，结合图8发现法兰①的line-z正好处于两个漩涡中间，而颗粒主要随着漩涡运动，所以浓度并不高，管道壁面附近浓度大于管道中心处。且line-y正好穿过两个漩涡，可以发现两个浓度波谷正好处于漩涡中心处，颗粒很难进入漩涡核心，所以浓度接近于0 mg/m³。图12为法兰⑨截面上line-z和line-y的浓度分布图，此时流场已经发展稳定下来，外壁面附近整体速度相比于内壁偏高一点，因此从line-z上可以看到内壁面浓度更高，从line-y上可以看到因为重力的存在，管道下部颗粒聚集得更多，浓度更高。并且法兰⑨截面上漩涡已经完全消失，z方向和y方向上在管道中心附近浓度差别不大，整体均匀性很好。

Figure 9. Position of line-y and line-z

图9. line-y、line-z的位置分布

Figure 10. Along the axial direction of duct on the duct section Φ value size

图10. 沿管道轴向方向上各管道截面Φ值大小

Figure 11. Particle concentrations of line-y and line-z on the flange ① section

图11. 法兰①截面上line-y和line-z浓度分布曲线

Figure 12. Particle concentrations of line-y and line-z on the flange ⑨ section

图12. 法兰⑨截面上line-y和line-z浓度分布曲线

4.3. 气流风速和颗粒质量流率的影响

图13为颗粒质量流率固定为13 mg/s时，不同风速下，不同法兰截面上line-z的平均浓度规律基本一致，先上升再下降后趋于平稳。先上升是因为弯管出口二次流很强，大量颗粒物随着漩涡移动道管道中心部分。随着漩涡强度减小，颗粒物开始扩散到整个管道，当漩涡完全消失，颗粒物充满整个管道，所以line-z上的平均浓度会再下降然后趋于平稳。

在不同风速下可以发现，随着风速的提高，颗粒物平均浓度下降，这是湍流增强的原因，由于流体速度的增加，湍流强度增加，颗粒受到的曳力增加，高速区的颗粒物更容易被流体带走。并且随着流速的增加，平均浓度下降趋势也逐渐减弱，这主要是因为气流速度的增加，湍流强度增加，颗粒受到的的曳力增加，高速区的颗粒物更容易被流体带走。并且随着流速的增加，平均浓度下降趋势也逐渐减弱，这主要是因为气流速度的增加，湍流强度增加，漩涡寿命也越长，完全消失的距离也越远，所以line-z上的平均浓度下降得更慢。

图14为风速固定为8.76 m/s时，随着颗粒物质量流率的增加，颗粒物平均浓度上升，这也与质量守恒相符。并且在越小的颗粒质量流率情况下发现平均浓度上升和下降趋势越缓慢，特别是在3.25 mg/s时，平均浓度的下降趋势特别小。原因主要是颗粒质量流率减少，相同时间进入管道内的颗粒量也降低，减少了颗粒与气流的动量交换，更有利于颗粒跟随气流运动，从而使得颗粒浓度场紧随流场发展趋势，更快的发展开来，从而看到line-z上的平均浓度下降趋势减弱。

Figure 13. Variation of mean concentration at different flange sections line-z at different wind speeds

图13. 不同风速下不同法兰截面line-z上平均浓度变化曲线

Figure 14. Variation of mean concentration of line-z of different flange sections at different particle mass flow rates

图14. 不同颗粒质量流率下不同法兰截面line-z的平均浓度变化曲线

5. 结论

1) 使用RNG k-ε湍流模型和欧拉–拉格朗日模型模拟了细小颗粒在风洞管道中的运动扩散情况，将模拟结果与实验数据比较，误差在可接受范围内，说明本次的模拟方法具有可靠性。

2) 基于气流流场、速度云图，速度矢量图和颗粒物浓度图说明了细微颗粒物的扩散规律。影响颗粒物分布均匀性的最大因素是经过弯管之后形成二次流的漩涡，漩涡起到了混动颗粒的作用，但也是造成颗粒分布不均匀的直接原因，颗粒分布的均匀程度与漩涡的强度呈反比关系。

3) 颗粒物浓度随着气流速度的增加而下降；随着入口颗粒物质量流率的增加而增加。为了保证颗粒浓度分布的均匀性，在风速8.76 m/s和颗粒质量流率3.5 mg/s时，平均颗粒浓度曲线变化更加平缓。

NOTES

^*通讯作者。