细微颗粒在方腔自然对流中运动特性的格子Boltzmann模拟

李勇

（华中科技大学煤燃烧国家重点实验室 湖北武汉 430074）

细微颗粒在方腔自然对流中运动特性的格子Boltzmann模拟

李勇

（华中科技大学煤燃烧国家重点实验室 湖北武汉 430074）

为了了解颗粒物在气体中的流动规律，为大气的颗粒物污染提供理论上的指导，本文将不可压双分布热格子模型与基于点源颗粒的拉格朗日跟踪法相结合，数值研究了封闭方腔自然对流中的细微颗粒物的运动特性。分析了瑞利数和粒径大小对颗粒群流动特性的影响。结果表明，瑞利数较小时，颗粒在流场中会形成一个准平衡的循环区域，此时瑞利数和粒径大小都会对颗粒的运动特性产生影响；而当瑞利数增加到一定临界值时，不会形成类似的准平衡的循环区域，此时粒径大小对颗粒的运动影响几乎可以忽略。

格子Boltzmann方法；气固两相流；方腔自然对流

1 引言

鉴于近年来频繁出现的雾霾天气，对于细微颗粒物对大气的污染治理及控制越来越紧迫，因而研究颗粒物在流体中的运动规律具有重大意义。颗粒的热泳沉积被视作一种有效的除尘手段，一些学者对这类问题进行过研究。1998年，Thakurta等人应用DNS对细微颗粒在管道流中的热泳沉积进行了数值研究[1]。Tsai等人2001年发展了一个简单的关系式来评估热泳作用对气溶胶颗粒沉积作用的影响[2]。2009年刘若雷等人研究了可吸入颗粒物在温度梯度场内的运动特性，得到了PM2.5热泳沉积的半经验公式[3]。同年，Akbar等人对细微颗粒在方腔自然对流中的运动进行了模拟[4]，结果表明热泳力和布朗扩散效应对颗粒的运动起主要作用。然而，他们的结果并没有指出当瑞利数和粒径大小变化时这两种作用力对颗粒运动特性的影响。因此，本文将对颗粒群在方腔自然对流中的运动进行数值模拟，并研究瑞利数和粒径大小对颗粒运动特性的影响。

随着大型计算机性能的不断提高，数值模拟方法被迅速应用到气固两相流领域的研究中[5-8]。相对于传统的数值计算方法而言，格子Boltzmann方法（Lattice Boltzmann Method，LBM）具有演化过程清晰、编程容易和适合并行计算等优势。因此LBM被广泛应用气固两相流的数值模拟[9]。郭照立等人2002年提出的不可压双分布热格子模型（DDF-LBE）具有计算简单、数值稳定性良好等优点[10]。因此本文使用上述热格子LBM模型，并与拉格朗日跟踪法结合，模拟颗粒群在方腔自然对流中的流动，研究瑞利数与粒径大小对颗粒群运动特性的影响，希望能为细微颗粒的污染提供一些理论上的指导。

2 数值方法

2.1 气相流动的格子Boltzmann方法

本文考虑二维问题，采用郭照立等人2002年提出的不可压双分布（DDF-LBE）热格子模型模拟流场。该模型用两组演化方程分别模拟速度场和温度场，其演化方程如下：

其中fi，Ti分别为速度场和温度场的离散分布函数，对应的平衡态分布函数定义为：

f(ieq)的参数设定为

模型中的作用力项定义如下：

流体的宏观速度、压力、温度分别计算如下：

各类参数的选取参照文献[10]中的设置。

2.2 颗粒运动的拉格朗日跟踪方法

对于固体颗粒的运动模拟，采用基于点源颗粒的拉格朗日跟踪方法，主要考虑重力，曳力，布朗力和热泳力对颗粒的作用，并且忽略颗粒对气相的作用以及颗粒间的相互作用。在上述假设下，颗粒的运动方程可以表示为：

其中，坎宁安校正系数C与颗粒的布朗力FB根据文献[11]确定，热泳力根据文献[4]中的半经验公式计算。

对方程(6)和(7)在一个时间步长内积分，可以得到颗粒速度和位移的计算格式：

3 模拟结果及讨论

运用上述LBE-拉格朗日跟踪法，我们对细微颗粒在二维封

闭方腔自然对流中的流动进行了模拟。模拟时所采用网格大小是256×256，封闭方腔的上下平板均是绝热壁面，左平板为高温壁面，右平板为低温壁面。方腔自然对流的瑞利数其中L是特征长度，v和D分别为流体的运动粘度和热扩散率，β为热体积膨胀系数。起始状态颗粒以初速度为零随机分散在方腔内，颗粒数目为2000。计算过程中，假设颗粒和壁面发生接触即被壁面吸附。

首先，我们模拟了Ra=1000、颗粒直径为1μm时的工况。图1给出了不同时刻颗粒在方腔中的运动状况。从图中可以发现，一部分颗粒在热泳力等的作用下很快与冷壁面接触而被吸附。部分颗粒运动到一定的区域内由于力的作用达到平衡而处于类稳定状态，这部分颗粒会一直在该区域内运动而不会被冷壁面吸附，形成一个准平衡的循环区域。另外，也有一小部分颗粒由于重力的作用会沉积在下平板上。

图1 Ra=1000时1um颗粒在不同时刻的位置分布: (a)0s,(b)5s,(c)25s,(d)50s

我们现在考虑不同瑞利数条件下颗粒群在封闭方腔内的流动状况。颗粒直径仍设定为1μm。图2给出了瑞利数从100增加到8×105颗粒群在t=50s时的位置分布。可以发现，当Ra较小时，颗粒最终均会形成一个准平衡的区域并在该区域内运动。并且随着Ra的增大，该准平衡区域越来越分散。当瑞利数为时，颗粒在方腔内随机分布而不会形成类似的准平衡区域。

图2 不同瑞利数Ra下颗粒群的位置分布: (a)100,(b)103,(c)104,(d)8×105

图3 dp=1um时不同瑞利数的颗粒保留率

为了更好地理解颗粒群的运动特性，我们统计了方腔中未被壁面吸附的颗粒数占颗粒总数的比率，即颗粒保留率RF(Remaining Fraction)。不同瑞利数Ra下的颗粒保留率RF随时间的变化曲线如图3所示。初始阶段所有颗粒随机分布在方腔内，在热泳力等的作用下，处于前述准平衡区域之外的外部颗粒迅速向冷壁面运动并被吸附，颗粒保留率RF随之减小；而运动到准平衡区域之内的内部颗粒将保持在该准平衡区域内运动而不会被吸附，RF最终达到一个稳定值。并且Ra越小，RF减小得越快，最终达到的稳定值也越小。这是因为Ra越小，颗粒达到受力平衡的准平衡区域越接近方腔中心，颗粒更加集中，因此更多的颗粒处于准平衡区域外部，它们在热泳力等的

作用下向冷壁面运动并被吸附。随着Ra增大，该准平衡区域逐渐分散和模糊，最终消失。对于图中Ra最大的情况，由于方腔内没有明显的准平衡区域存在，颗粒将逐渐运动到壁面上并被吸附，因而颗粒保留率呈单调下降趋势。

同时，考虑粒径对流动过程中颗粒保留率RF变化曲线的影响。我们统计了粒径分别为1μm,200nm,50nm时不同Ra下的颗粒保留率RF随时间变化，如图4所示。由图4中可知，Ra较小时，粒径大小对颗粒保留率的影响比较明显。粒径越小，流动过程中颗粒保留率RF减小得越快，并且最终达到的稳定值也越小。这是因为对于细微颗粒，布朗运动的剧烈程度较大，布朗扩散的加剧使颗粒的沉积速度增加。而随着瑞利数增加，热泳力的作用越来越强烈，粒径大小的影响则逐渐减弱；当瑞利数增加到足够大时，热泳力对颗粒的影响起主导作用，粒径大小的影响几乎可以忽略，如图4(d)。

图4 颗粒直径对颗粒保留率的影响: Ra=(a)100,(b)103,(c)104,(d)8×105

4 结语

为了研究颗粒物在气体中的运动规律，为颗粒污染物的治理提供一些理论指导。本文将不可压双分布（DDF-LBE）热格子模型与拉格朗日跟踪法相结合，对颗粒群在方腔自然对流里的流动进行了数值研究。结果表明，当瑞利数较小时，颗粒在流场中会形成一个准平衡的循环区域，在这个区域内的颗粒保持准平衡状态，不会被壁面吸附。随着瑞利数的增大，该准平衡区域逐渐分散、模糊，最终消失，此时颗粒逐渐向壁面运动并被吸附。当瑞利数较小时，粒径大小对颗粒群在流体中的流动特性影响较为明显。粒径越小，在整个流动时间段里的颗粒保留率减小得越快，最后达到的稳定值也越小；随着瑞利数增大，粒径大小的影响逐渐减弱；当瑞利数逐渐增大到某一个临界值时，粒径大小对颗粒运动特性的影响基本可以忽略。

[1] Thakurta D.G., Chen M., Mclaughlin J.B. etc. Thermophoretic depositionof small particles in a direct numerical simulation of turbulentchannel flow [J]. International Journal of Heat and Mass Transfer,1998,41,4167- 4182.

[2] Tsai R., Liang L.J. Correlation for thermophoretic deposition ofaerosol particles onto cold plates [J].Journal of Aerosol Science,2001,32,473- 487.

[3]刘若雷，杨瑞昌，由长福等.温度梯度场内可吸入颗粒物运动特性及热泳沉积[J].化工学报,2009,60,1623-1628.

[4] M.K.Akbar, M.Rahman,S.M. Ghiaasiaan. Particle transport in asmall square enclosure in laminar natural convection [J]. Aerosol Science,2009,40,747- 761.

[5] U. Lantermann, D. Haanel. Particle Monte Carlo and lattice- Boltzmannmethods for simulations of gas- particle flows [J]. Comput. Fluids,2007, 36∶ 407- 422.

[6] I. Alexander. Effects of the variation of mass loading and particledensity in gas- solid particle flow in pipes [J].Powder Technol,2009, 193∶176-181.

[7] Balachandar S, Eaton JK. Turbulent dispersed multiphase flow[J].Annual Reviewof Fluid Mechanics,2010,42∶ 111- 133.

[8] Durst F, Milojevic D, Schonung B. Eulerian and Lagrangian predictionsof particulate two- phase flows∶ A numerical study [J].AppliedMathematical Modelling, 1984,8(2)∶101- 115.

[9]郭照立,郑楚光.格子Boltzmann方法的原理及应用 [M].北京∶科学出版社,2009.

[10] Z.L.Guo, B.C.Shi, and C.G.Zheng. A coupled lattice BGK modelfor the Boussinesq equations[J]. Int.T.Numer.Fluids,2002,39, 325- 342.

[11] Friedlander, S. K. Smoke dust and haze. 2nd ed [M]. Oxford UniversityPress,2000.

[12] M.Honrrmann and M.Peric. Finite volume multigrid prediction oflaminar natural convection∶Benchmark solutions [J]. Int. J. Numer.Method.Fluids,1990,11,189- 207.

李勇（1989—），男，汉族，湖北省武汉市人，华中科技大学煤燃烧国家重点实验室，硕士研究生，主要研究领域为气固两相流动的计算数值模拟。