基于直接力浸入边界法的线形排列多球颗粒沉降特性研究

丁 周，周 锟，孙 科，贺 铸

(武汉科技大学省部共建耐火材料与冶金国家重点实验室，湖北 武汉，430081)

颗粒在流体中的自由沉降属于固液两相流[1-2]，普遍存在于自然界中，如泥沙的沉积[3-4]、大气中水蒸气以及雾霾颗粒的运动等[5-6]，同时，颗粒沉降还具有广泛的工业生产背景，其相关研究在管道输运[7]、冶金化工流化床设备[8]等一些工程实际问题中也具有重要意义。颗粒沉降过程复杂，受到颗粒的粒径、形状、密度以及流体黏度等众多因素影响。

针对一般的颗粒流动问题，常用的数值模拟方法有欧拉-欧拉法(Euler-Euler method)[9]、欧拉-拉格朗日法(Euler-Lagrange method)[10]和颗粒解析法(particle-resolved method)[11]等三大类，这三类方法依次适用于颗粒非常小(如微纳尺寸的气溶胶颗粒)、颗粒比较小(如粉尘)以及颗粒相对比较大的情况。Peskin在研究心脏血管流动时提出的浸入边界法(immersed boundary method，IBM)[12-13]就是一种非常典型的颗粒解析法，该法的主要特点是利用统一、结构化的欧拉网格来离散整个流固区域以及采用拉格朗日网格跟踪流固界面，并借助delta函数实现二者的交互。传统浸入边界法基于弹性边界，不适合对刚性边界的颗粒流动现象进行数值模拟，而Uhlmann[14]则通过引入直接力所提出的直接力-浸入边界法(direct forcing-immersed boundary method, DF-IBM)很好地解决了这一问题，借助DF-IBM对二维条件下的双圆形颗粒沉降以及三维条件下单个和多个球形颗粒沉降进行研究时，模拟结果与相关的实验结果[15]基本吻合，Di等[16]采用DF-IBM对经典的二维双圆形颗粒的牵引、碰撞及翻滚(DKT)[17]进行模拟，也获得了较理想的研究结果，此外，还有宫兆新等[18]对浸入边界法进行了较全面的研究。一直以来，在DF-IBM中欧拉网格和拉格朗日网格的尺寸被认为应该相互匹配，并且对于每个拉格朗日节点都需要定义一个与欧拉网格匹配的权重积分，然而本课题组最新的研究[19]表明，拉格朗日体积权重无需特殊选择，可以在很大范围内自由取值。

当前，有关二维条件下双圆形颗粒和三维条件下双球颗粒的沉降实验与模拟方法已经比较完备，但针对多颗粒的沉降问题仍停留在宏观分析阶段，考虑到实际颗粒输运过程中往往是多个颗粒一起沉降或者提升，故本文分别以单个球形、线形排列的双球和多球颗粒为研究对象，借助直接力浸入边界法对相应的沉降过程进行数值模拟研究，并对刚体假设(rigid body assumption,RBA)[20]的结果准确性进行分析，重点探讨了线形排列多球颗粒的沉降行为。

1 模型与方法

1.1 物理模型

(a)模拟区域 (b) 球形颗粒表面离散点图1 物理模型Fig.1 Physical model

1.2 直接力浸入边界法

采用直接力浸入边界法对颗粒自由沉降进行数值模拟，其数理模型中不可压缩流体的控制方程为

(1)

(2)

式(1)～式(2)中:u、p分别为流体的速度和压力；为哈密顿算子；fIBM为数值模拟中欧拉网格点处的作用力。颗粒在沉降时的运动方程为

(3)

式中：mP为颗粒质量；t为时间；τ为应力张量，τ=-Ip+vf(u+uT)，其中I是二阶单位张量,T表示矢量矩阵的转置；n为颗粒表面法线方向；GP为颗粒所受重力，GP=VP(ρf-ρP)，其中VP为颗粒体积，ρP为颗粒密度。式(3)中的相互作用力项

(4)

式中，χ表示欧拉节点。采用RBA近似处理式(4)等号右边第2项，有

(5)

(6)

式中，fIBM可表示为

fIBM=BWFIBM

(7)

式中:B为插值函数，W为控制收敛的权重因子，FIBM为颗粒表面拉格朗日离散点上的作用力，也就是说fIBM可由颗粒表面拉格朗日离散点上的作用力FIBM经插值处理而获得，而FIBM又可由颗粒与流体在界面上的速度差与时间步长的比率来确定，即

BiI=ψ(χi-ΧI)

(8)

(9)

(10)

(11)

(12)

需要指出的是，本研究模拟颗粒尺寸远大于流体分子的平均自由程，因此颗粒表面所受流体分子的随机碰撞处于平衡状态，也就是说，流体的分子热运动对颗粒的影响是各向均匀的，其宏观作用表现为流体的黏性。对于液体来说，其分子平均自由程为纳米量级，而对于常温下的空气，其分子平均自由程为亚微米量级，所以，本文模拟适用于颗粒尺寸为微米量级以上的情况，文中颗粒尺寸单位均为其尺寸量级所对应的国际单位制单位及其导出单位。

2 模拟方法验证

为进一步验证数值模拟对网格的依赖性，在3组不同网格密度条件下对类似的颗粒沉降问题进行了模拟，模拟结果见图3。其中NP表示颗粒拉格朗日离散点数目，不同的拉格朗日离散点数对应不同密度尺寸的欧拉网格，相应的误差分析如表2所示。结合图3及表2结果可知，当颗粒拉格朗日离散点为258，相应槽道欧拉网格密度为68×544×69即Dp/h=9.063时，模拟计算结果较其它网格密度条件下的模拟值更为准确，采取此网格进行模拟计算，模拟过程中Courant-Friedrichs-Lewy(CFL)数均小于1。

表1 算例的对比结果

(a)算例1

(b)算例2

(c)算例3

图3 不同网格密度条件下的颗粒沉降速度Fig.3 Particle sedimentation velocities under various grid densities

表2 误差分析

3 结果及讨论

3.1 刚体假设(RBA)的影响

基于文献[14]中的算例2，通过计算，分析了不同RBA条件下，式(5)中积分项(等号左端)及乘积项(等号右端)的数值随时间变化的规律，结果如图4所示。图4中，标记为RBA乘积值和非RBA乘积值的曲线分别反映了颗粒是否为刚体时颗粒体积与颗粒质心沉降速度的乘积随时间变化的趋势，标记为RBA积分值和非RBA积分值的曲线分别反映了颗粒是否为刚体时借助Monte Carlo随机积分法计算所得式(5)中积分项的值随时间变化的趋势，作为参照，标记为RBA实验乘积值的曲线则反映了颗粒为刚体时，由文献[15]实验值所得出的颗粒质心沉降速度与颗粒体积的乘积随时间变化的趋势。上述积分项、乘积项值随时间的变化趋势即为颗粒沉降速度随时间的变化趋势。由图4可见，将颗粒看作刚体时，相应的计算结果更接近于文献[15]所报道的实验值；当颗粒不是刚体时，计算值和实验值的误差较大。此外，将颗粒作为刚体处理时，相应的乘积值曲线相较积分值曲线更为平滑、波动更小。以上结果表明，在球状颗粒沉降过程中，当颗粒为刚体时，相应计算结果更加准确。

图4 不同RBA条件下方程(5)的计算值Fig.4 Calculated values of equation (5) under different RBA conditions

3.2 双球颗粒沉降

在双球颗粒沉降模拟中，设定2个球体颗粒间球心的初始间距分别为1DP、1.5DP和2DP。在颗粒沉降过程中，不同初始间距条件下，2颗粒相对间距Dr随时间t变化的模拟结果如图5所示。由图5可见，在球状颗粒沉降过程初始阶段，2球间发生相互牵引致使二者间距逐渐缩短，这与经典双球沉降模型[17]中的牵引过程一致，当2球相对距离减小到一定数值时则会发生2球相互触碰，初始间距不同的2球均能发生牵引和触碰，且发生触碰的时刻随2球初始间距的增大而推迟。随着颗粒沉降的继续进行，2球相互触碰后又发生翻滚运动导致二者分离，2球相对间距又开始逐渐增大(图5(a))。

(a)初始间距1DP

(b)初始间距1.5DP

(c)初始间距2DP

图6所示为双球颗粒沉降过程中的速度云图及压力云图。需要特别指明的是，在DF-IBM方法中，由于采用统一的欧拉网格来处理整个区域(包含颗粒所占据的空间)，因此也有流体填充在颗粒内部，而颗粒内部的速度场，并不完全满足颗粒的刚体运动(见图6(a))，此问题的根源还有待进一步深入研究。从双球的压力云图(图6(b))中可以清楚看出，单个球体的下方为正压区，上方为负压区。当2号球处于1号球的负压区时，2号球所受压差阻力减小并导致其加速沉降，进一步的分析表明，1号球的沉降速度较其单独沉降时要小，也就是说，2号球在1号球的牵引下加速运动，最终追上后者发生触碰。

(a)速度 (b)压力图6 双球颗粒沉降云图Fig.6 Sedimentation nephograms of double spherical particles

3.3 线形排列多球颗粒沉降

线形排列多球模型之于双球模型，可看作在双球上方等间距放入一定数量的小球，原有的双球模型在加入更多球形颗粒后成为具有多球系统的新模型。多球颗粒系统沉降过程中的速度云图及压力云图如图7所示。从图7中可以看出，当颗粒已经沉降一段距离后，在球形颗粒数分别为3、4、5、6的算例中，1号球和2球的间距均小于其它相邻两球之间的间距，这表明1号球和2号球率先发生牵引，并且在不同的球形颗粒数条件下，1号球下方的正压梯度均大于其它球的相应值。相比双球模型，多球系统中的2号球上方还存在由3号球造成的正压区，因此，在3号球下方正压区及1号球上方负压区的共同作用下，2号球会更快地追上1号球。

(a)速度

(b)压力

因为当线形排列的球状颗粒数超过3以后，伴随着1号、2号球发生DKT(牵引，触碰、翻滚)过程，1～3号球的运动状态基本一致，故而对1～3号球中相邻球间的相对间距变化进行分析，结果如图8所示，其中P2-P1表示1号球和2号球的相对间距，P3-P2表示2号球和3号球的相对间距。由图8可见，在颗粒沉降过程中，1号球和2号球先发生牵引和触碰，然后再分离，这与二者在双球模型中的表现一致，区别在于整个过程因3号球的加入而导致进度加快。同时，3号球和2号球仅仅发生牵引而不会产生触碰。此外，本研究模拟结果还表明，随着线形排列颗粒数目的不断增加，使用DF-IBM对相应沉降过程进行模拟计算的效率始终保持在较高水平。IBM模拟主要采用欧拉网格及拉格朗日网格，其中欧拉网格的位置和数量是固定不变的，而拉格朗日网格的变化则依赖于颗粒数目。在模拟计算中，计算量主要取决于欧拉网格的数量，这是因为在整个使用了压力投影法的模拟过程中，对压力场的泊松方程进行求解耗时最多，而泊松方程又是在欧拉网格上离散求解的。一般来说，颗粒越多，计算量越大，但是IBM法中欧拉网格数量固定，所以计算量的增加幅度非常有限，这也是IBM方法目前应用比较广泛的原因之一。

图8 多球颗粒模型的Dr-t曲线

3.4 线形排列双球与多球颗粒沉降对比

借助双球和三球颗粒模型对比分析了沉降过程中1号球和2号球的相对间距与碰撞时间的关系，结果如图9所示。由图9可见，当相邻颗粒之间的初始间距从1DP增大到4DP时，即随着Dr的增加，双球和三球模型中的1号球和2号球均能发生触碰，而且在三球模型中，因3号球的引入导致1号、2号球触碰时间提前。

当相邻两球颗粒初始间距为1DP时，在双球和三球模型中发生牵引的1号球和2号球的沉降速度随二者相对间距变化的曲线如图10所示，其中2P-P1、3P-P1分别表示双球和三球模型中的1号球，即牵引球；2P-P2和3P-P2分别表示双球和三球模型中的2号球，即受牵引球。从图10中可以看出，在沉降开始后的一段时间内，两种模型中1号球和2号球都在加速沉降且二者相对间距基本不变，紧接着双球模型中的两球间距稍有增大然后再减小，而与此同时，三球模型中发生牵引的1号球和2号球间的相对间距却一直在缩小，表明三球模型中的牵引现象较双球模型先发生，应归因于三球模型中3号球的引入加速了这一过程。

图9 双球和三球颗粒模型的碰撞时间Fig.9 Collision time of double and three spherical particles models

图10 双球和三球颗粒模型的沉降速度Fig.10 Sedimentation velocities of double and three spherical particles models

在球形颗粒沉降过程中，针对1号球和2号球在三球模型中较双球模型先发生牵引的问题，从球形颗粒所受阻力角度进行分析，两种模型中1号球及2号球所受阻力变化情况如图11所示。由图5与图8结果可知，沉降时间t超过0.5后，1号球和2号球之间的相对间距出现明显变化，而从图11(a)中又可以看出，受牵引的2号球在三球模型中所受阻力较其在双球模型中所受阻力要小，从而能在更短的时间内追上1号牵引球以实现牵引。在颗粒沉降初始阶段，1号牵引球在三球模型中所受阻力与其在双球模型中所受阻力几乎一样，当t超过0.7时，1号牵引球在三球模型及双球模型中所受阻力出现明显差别(图11(b))，这是因为相同时刻三球模型中1号球和2号球相对间距小于双球模型中的相应值，导致1号球所受斥力增大而阻力减小。

(a)受牵引球(2号球)

(b)牵引球(1号球)

4 结语

当球形颗粒初始间距介于一定范围内且视球形颗粒为刚体时，竖直线形排列的双球或多球颗粒(颗粒按位置从低到高依次排序)中的1、2号球均可进行DKT(牵引、触碰、翻滚)运动。在双球颗粒基础上引入数量不等的球形颗粒形成线形排列多球颗粒，能改变2号球的受力从而加速1号球和2号球DKT运动的发生，这一现象对引入颗粒的数量不敏感。本研究方法可直接扩展应用到任意非规则形状颗粒的沉降研究，能为颗粒沉降问题分析提供一定的参考。