王 婧 ，王军武

（1.中国科学院过程工程研究所，多相复杂系统国家重点实验室，北京 100190；2.中国科学院大学化学工程学院，北京 100049）

气固两相流在工业中有广泛的应用，因此，研究气固两相流的流动机制对于反应器放大和设计具有重要的意义。近年来，计算机模拟已经成为研究气固两相流的主要手段之一，其中连续介质模型经常应用于工业级反应器的模拟[1-2]，它通常采用颗粒动理论封闭固相应力[2]。

在颗粒动理论中，颗粒速度分布函数是核心。而颗粒动理论关系到连续模型中本构关系的封闭，因此确定颗粒速度分布函数的形式是建立正确本构关系的关键。在气固两相流中，颗粒-流体间的黏性作用、颗粒-颗粒以及颗粒-固体壁面间的非弹性碰撞引起的能量耗散，使得系统处于非平衡状态，导致颗粒速度分布有别于常规气体中的麦克斯韦分布。研究表明，在气固两相流中，颗粒速度分布函数可以呈现麦克斯韦分布[3]、指数分布[4-5]、双峰分布[6-8]、t分布[9]等形式。在循环流态床的顶端等颗粒浓度较低的地方，颗粒速度分布函数基本满足麦克斯韦分布，而在环核结构界面等非均匀结构起主导的位置，颗粒速度函数一般为指数分布、t分布或双峰分布。不同的速度分布会影响到颗粒相本构关系的建立，但是从文献调研中我们发现，颗粒速度分布在不同工况下具有不同的函数形式，为了完善更加全面的颗粒动理论，我们将分析这几种速度分布函数在颗粒流体系统中的普适性。

利用基于连续介质的计算流体力学方法（CFD）和离散单元法（DEM）耦合在一起，通过CFD-DEM耦合法模拟研究气固流化床中粗颗粒的速度分布函数形式，并与以上各种分布函数进行比较，进而找到非平衡系统中相对合理的颗粒速度分布函数形式，为颗粒动理论的发展提供支撑。

1 数值实验

1.1 CFD-DEM方法

大量的研究表明，CFD-DEM方法可以模拟气固流化床内的流体力学行为，特别是密相粗颗粒流态化的行为[10-11]，因此，直接采用CFD-DEM方法，不对CFD-DEM方法本身的合理性进行论证。

CFD-DEM方法[12-13]采用平均化的Navier-Stokes方程描述气相的运动，采用牛顿第二定律描述颗粒相的运动。它在处理颗粒和颗粒之间的相互作用时，采用“软球模型”[14-15]，这种处理方式考虑到颗粒碰撞有一定的接触时间，通过弹性、阻尼和滑移的力学特征来描述颗粒碰撞过程的形变、缓冲和滑移的接触过程。另外，在计算相间耦合时，采用Beestra关联式[16]来计算气固相间曳力。CFD-DEM模型的控制方程为式（1）—式（10）。

颗粒运动受力方程：

式中，ma为颗粒质量，ra为颗粒位置，Va为颗粒体积，p为压力，β为相间曳力系数，εg为气相浓度，ug为气相宏观速度，vs为单个颗粒的速度，g为重力加速度，Fc，a为颗粒碰撞时受到的接触力。

两颗粒间碰撞受力：

式中，Fab，n为颗粒碰撞时法向接触力，kn为法向刚性系数，δn为颗粒碰撞时发生的形变量，nab为法向单位向量，ηn为法向阻尼系数，vab，n为法向相对速度，Ra和 Rb分别表示颗粒a和颗粒b的颗粒半径。

气相的质量守恒方程：

气相的动量守恒方程：

式中，ρg为气相密度，τg为气相应力。

气固相间曳力Sp：

Vc为计算网格的体积。

相间曳力系数β：

理想气体状态方程：

式中，ρg为气体密度，dp为颗粒直径，Mg为气体摩尔质量，R为摩尔气体常数，Tg为气相温度。

气相应力本构关系：

式中，μg为气相的剪切黏度，I表示单位矩阵。

1.2 模拟参数

采用CFD-DEM模型模拟D类颗粒的鼓泡流态化、湍动流态化、循环流态化这3种两相流中常见的流型，模拟采用拟三维流化床，即气相是二维流动，固相是三维流动，床层的厚度是颗粒直径的6.1倍（采用周期性边界）。我们将流化床的水平方向设置为无滑移的壁面，进口使用速度边界，出口使用压力边界。D类颗粒流态化的模拟参数见表1。

表1 D类颗粒流态化的模拟参数Tab.1 Simulation parameters for the fluidization of Geldart D particles

1.3 速度分布函数

如果颗粒流体系统是均匀流动或系统处在平衡状态，颗粒速度分布函数通常满足麦克斯韦分布[2]，

式中，c为颗粒微观速度，f（c）表示颗粒速度分布函数，θ为颗粒温度，u为颗粒的平均速度。

但是在气固两相流系统中，大多数情况下，体系处于非平衡态并且伴随着热传导、扩散、黏滞现象等传递过程。如果颗粒流体系统中，颗粒处于非平衡状态，颗粒的速度分布就会偏离麦克斯韦分布呈现出长拖尾的趋势，或者呈现多个麦克斯韦分布叠加的形状等。

式（12）表示的是颗粒速度分布呈指数分布的形式：

式中k、α、β均为修正系数。

式（13）是呈双峰分布的速度分布函数，它是基于能量最小多尺度模型（EMMS）思想得到的。

式中，f为密相体积分数，nc为密相数密度，（1-f）为稀相体积分数，nf为稀相数密度，且有（1-f）nf+fnc=n，n 是颗粒总的数密度，θf是稀相颗粒温度，usf为稀相颗粒的平均速度，θc为密相颗粒的温度，usc为密相颗粒的平均速度。

式（14）是呈t分布的速度分布函数，t分布与麦克斯韦分布相比具有长拖尾特征[17]。

这里 υ表示自由度，并且 υ＞2，（θ+Θ）是修正后的颗粒温度。式（14）中指的是函数，它的定义如下：

2 结果与讨论

2.1 水平方向的速度分布函数

通过CFD-DEM模型计算出颗粒x方向的速度，根据DEM数据统计得到x方向的颗粒速度分布函数，然后，将上述的几种颗粒速度分布函数和数值实验得到的速度分布进行对比。图1是鼓泡流态化下x方向颗粒速度分布函数。从图可以发现，t分布对颗粒速度函数的预测比较接近于模拟数据得到的速度分布。

图1 鼓泡流态化下x方向颗粒速度分布函数Fig.1 PDF of x velocity in bubbling fluidization

图2为湍动流态化下x方向的颗粒速度分布。而从图2可以看出，湍动流态化下x方向的速度分布函数更接近指数分布或双峰分布。与鼓泡流态化下的速度分布相比，随着气速的增加，可以发现速度分布的区间增大，在这种状态下，指数分布或双峰分布可以较为合理地回归出模拟得到的速度分布函数。

图2 湍动流态化下x方向的颗粒速度分布函数Fig.2 PDF of x velocity in turbulent fluidization

随着气速的继续增大，速度分布的区间明显增大，速度分布也出现了明显的长拖尾趋势。图3是循环流态化下x方向的颗粒速度分布函数。从图可以看出，t分布能够预测出这种长拖尾特征，但是与模拟结果仍存在一定的偏差。我们可以看到在颗粒速度达到3 m/s以上由模拟数据得到的速度分布更离散，这是因为在这个范围的颗粒速度出现比较少，统计样本数不足造成的。

图3 循环流态化下x方向颗粒速度分布函数Fig.3 PDF of x velocity in circulating fluidization

从图1—3可以发现，麦克斯韦分布整体上偏离模拟数据较多，这是因为麦克斯韦分布仅仅适用于描述平衡系统中固体颗粒的速度分布函数，但是气固流态化系统是一个典型的非平衡系统，因此需要用其他形式的分布函数来描述。

2.2 气流方向的速度分布函数

采用2.1节同样的方法，我们比较z方向的速度分布函数，见图4—6。从图4、5、6可知：CFD-DEM模拟数据统计出来的速度分布不再是对称分布，在不光滑的地方较多，这可能与实验本身的模拟条件有关。我们模拟的是一种拟三维的微型流化床，这有可能使得模拟出的速度分布呈现明显的边界效应。

图4 鼓泡流态化下z方向颗粒速度分布函数Fig.4 PDF of z velocity in bubbling fluidization

图5 湍动流态化下z方向颗粒速度分布函数Fig.5 PDF of z velocity in turbulent fluidization

图6 循环流态化下z方向颗粒速度分布函数Fig.6 PDF of z velocity in circulating fluidization

z方向的速度分布函数显示，在气体流动的方向也就是z方向上气相和固相的相互作用更加明显，颗粒和颗粒之间的相互作用也比x方向更加剧烈，因此系统内的流动非平衡特征较明显，这有可能导致速度分布不对称甚至出现一些波动。

根据图中显示的几种速度分布函数对比可以发现，双峰分布对于这种远离平衡态的速度分布有较好的预测，但是在分布的尾端，双峰分布也不能很好地拟合模拟得到的数据，表明基于EMMS模型的双峰分布仍需要进一步完善。

3 结论与展望

将DEM模型计算得到的颗粒速度分布和目前提出的各速度分布函数进行比较，发现麦克斯韦分布不能准确预测出系统内颗粒速度分布的形貌，颗粒流体系统中颗粒的速度分布有“长拖尾”的特征。基于EMMS模型的双峰分布，从物理意义上更贴近气固两相流的流动机制，适合于描述远离平衡态下的颗粒速度分布函数。

颗粒速度分布函数在不同方向上具有不同的分布特性，水平方向上气固间作用较弱时颗粒速度满足t分布，竖直方向上气固耦合作用明显，两相流动控制机制协调竞争，颗粒速度满足双峰分布。本文中的模拟研究对今后颗粒动理论的发展提出了新的思路，要求颗粒动理论中可以建立不同方向上满足不同速度分布特性的本构关系，从而得到更加合理的本构关系。