单个中性球形颗粒在三维方腔中的运动

李 洋，夏振华

(浙江大学 航空航天学院 流体工程研究所，杭州 310027)

0 引 言

含颗粒多相流现象广泛存在于自然界和工业生产中，例如流化床、煤燃烧、种子干燥技术、粉体悬浮液输运、颗粒分离[1-2]、岩土矿石开采中的多孔介质流动等[3]。理解这种类型的多相流对于基础研究和工程应用都是至关重要的。

很多学者对简单流动中的颗粒运动进行了大量的研究，诸如颗粒在管道中的沉降[4-5]、颗粒在库特流中的运动[6-8]和颗粒在管道中的惯性迁移[9-12]。后者称为惯性聚焦现象， 其广泛地应用于微通道设备中有限尺寸颗粒的控制和分离[13-14]。在与管道相连的方腔内，Khojah 等[15]发现大（小）颗粒在高雷诺数时向外（内）迁移，而低雷诺数时大（小）颗粒的轨迹靠近（远离）涡中心。Haddadi 等[16]和Jiang 等[17]也考虑了不同纵横比的方腔的工况，这些工况也可以用来实现颗粒的分离。

学者们采用实验方法[18-21]和数值方法[22-27]对顶盖驱动方腔流进行研究。方腔几何形状简单，但是里面的流动表现出丰富的物理特征，被认为是数值验证的基准算例。关于无悬浮颗粒的顶盖驱动方腔流的更多信息，可参阅Shankar 和Deshpande 的综述[28]。

所有边界上的壁面约束使得顶盖驱动方腔流动不同于简单的槽道流动、管道流动。关于含颗粒的槽道流、管道流的研究已经发表了很多，但是对含颗粒的顶盖驱动方腔流的研究却很少。Tsorng 等[29]采用立体成像方法实验研究了中性悬浮颗粒在顶盖驱动方腔中的长时间运动轨迹（塑料球形颗粒直径为3 mm，方腔高度为10 cm）：颗粒在方腔一侧做螺旋运动；此外，颗粒的长时间轨迹似乎沿着内部循环模式的优先路径聚集。后来，Tsorng 等[30]对这项工作进行扩展，研究了刚性颗粒与流体密度比、雷诺数之间的影响关系。随着雷诺数和斯托克斯数的增加，颗粒轨道向方腔外围靠近。他们将颗粒的行为解释为旋转诱导的力迫使颗粒向优先路径运动，这与泊肃叶流中的Segre-Silberberg 效应类似[9]。之后，Kosinski 等[1]给出了颗粒团在二维方腔中运动的数值结果，其中雷诺数为1 000，采用双向耦合方法处理固体点颗粒与流体的相互作用。结果表明颗粒倾向于向壁面迁移，且斯托克斯数越大，迁移行为越强烈。Sidik 和Attarzadeh[2]数值模拟了不同雷诺数下二维方腔内固体点颗粒的运动。他们将模拟结果与之前的实验结果、数值结果进行比较，以证明三次插值法的应用价值。Safdari 和Kim[31-32]采用格子Boltzmann 方法耦合拉格朗日点颗粒追踪法，数值研究了含颗粒的三维方腔流动。他们观察到颗粒的运动轨迹主要取决于颗粒大小、方腔内的涡旋行为和颗粒密度。最近，Hu[33]利用格子Boltzmann 方法详细研究了单个中性椭球颗粒在二维方腔中的运动。他指出，极限环和颗粒初始朝向与位置无关，而且极限环随着粒径的增大向方腔中心收缩，随着雷诺数的增加向右下角迁移。

综上所述，有限尺寸颗粒在三维顶盖驱动流中的动力学尚未得到充分理解。本文采用多松弛格子Boltzmann 方法耦合插值反弹格式模拟了单个中性球形颗粒在三维方腔中的运动。展向的边界是弱受限的对称边界或强受限的固体壁面。着重考虑了初始位置、颗粒大小以及雷诺数对球形颗粒运动的影响。颗粒方腔流的研究能够加深对颗粒在方管方腔中分离的理解。

1 数学模型

1.1 格子Boltzmann 方法

格子Boltzmann 方法是求解Navier-Stokes 方程的方法之一[34-37]。为了保证数值稳定和参数灵活[38-41]，本文采用了多松弛的格子Boltzmann 模型。不可压流体的D3Q27 模型对应的颗粒分布函数的演化方程为：

1.2 颗粒的运动

颗粒的平动和转动由牛顿第二定律和欧拉方程控制：

2 问题设置和数值验证

2.1 问题设置

球形颗粒在顶盖驱动方腔中的运动如图1 所示，一个半径为R的中性悬浮球形颗粒浸没在高度为S的立方体方腔中。这里，中性颗粒是指密度等于流体密度的颗粒。方腔的上壁面以恒定的速度U沿着y方向运动，其底壁面和沿y方向的两个壁面均是静止的。在z方向采用了两种不同的边界条件。对于第一种情况（P0，下称为展向弱受限情形），展向壁面的限制被弱化，在展向两侧我们使用对称边界条件；对于第二种情况（P1，下称为展向强受限情形），在展向固壁使用无滑移边界条件。两个相关的无量纲参数雷诺数Re和颗粒斯托克斯数St，分别定义为：

图1 球形颗粒在顶盖驱动方腔中运动的示意图Fig. 1 Sketch of dynamics of a spherical particle in a lid-driven cubic cavity flow

2.2 数值验证

为了验证目前程序的准确性，首先模拟了Re=100、400 和1 000 下的单相顶盖驱动流，网格为96×96×96。不同雷诺数下，z= 0.5 平面上沿着水平中心线的速度u、沿着竖直中心线的速度v和Ku 等[27]结果的比较分别展示在图2（a）和图2（b）中，可以发现模拟结果与Ku 等的结果均吻合良好。

图2 P1 情况下，z = 0.5 平面上，水平中心线上速度u以及竖直中心线上速度v 的分布Fig. 2 In P1 case, the distribution of u-velocity on the horizontal center line and v-velocity on the vertical center line on the plane z = 0.5

接下来我们对颗粒在方腔中的运动进行了网格无关性测试。在模拟过程中，首先得到稳定的单相流场；然后球形颗粒加入到方腔中特定的位置，并允许其自由运动。对于P0 情况，两套网格被用来模拟Re= 1 000 下 的 颗 粒 方 腔 流。第 一 组 的 网 格 为96×96×96，其中颗粒半径为6 个网格，初始网格位置为（77，50，48）；第二组的网格为144×144×144，相应的颗粒半径为9 个网格，初始位置为（115.5，75，72）。对于颗粒在方腔中的运动，颗粒最终限制在展向z=0.5 平面的极限环上运动。如图3 所示，两组模拟下的极限环轨迹很好地重合在一起。此外，我们也发现两套网格下颗粒受到的力和力矩也能够较好地吻合。这些结果表明，96×96×96 的网格足以满足Re= 1 000下的球形颗粒运动的计算精度。在接下来的模拟中，96×96×96 的网格被采用。

图3 P0 情况下，两套网格的颗粒最终运动轨迹比较Fig. 3 In P0 case, comparison between the final trajectories of particles between two sets of mesh grids

3 展向弱受限的情形

3.1 初始位置的影响

这一部分我们主要考虑了颗粒初始位置对最终运动的影响。对于P0 情况，雷诺数Re= 1 000 和半径R= 6 下，对应的St= 0.868。颗粒在展向z= 0.5 平面上的不同初始位置释放，我们模拟了41 组算例，将颗粒的初始位置展示在图4 中，虚线为无颗粒时的流函数等值线，1、2、3 为三种最终的稳定运动轨迹。初始位于红色圆形、蓝色菱形或黑色三角形位置的颗粒最终分别在极限环3、极限环2、稳定点1 上运动。根据目前的模拟结果，颗粒的初始位置影响着其最终的运动轨迹，这里存在三种不同的运动模态。当颗粒初始放置在涡中心区域时，颗粒最终迁移到稳定点1，即涡中心，此时几乎以恒定的角速度绕展向做旋转运动；当颗粒初始放置在外侧的蓝色菱形点处，其最终沿着极限环2 做周期性运动；而当颗粒初始位置位于红色圆点时，球形颗粒最后在极限环3 上做周期性运动。

图4 P0 情况下，z = 0.5 平面上，颗粒从不同位置释放的相图Fig. 4 In P0 case, plane z = 0.5, phase diagram of particles released from different positions

根据图4 的结果，可以发现极限环3 的轨迹更靠近外围，它的空间位置不同于极限环2。为了进一步说明颗粒在极限环2 和极限环3 上运动的区别，我们在图5 中展示了初始网格位置分别位于（77，50，48）和（15，50，48）的两个算例，在相应极限环3 和极限环2 上运动时的合速度大小以及展向的角速度分量大小随时间的变化规律（为了方便画图，时间轴进行了平移，t0不同）。其中，颗粒速度和角速度分别由上壁面速度U和U/S无量纲化。由图5 可以发现，在两个极限环上运动的颗粒，它们的速度和角速度在具体数值上体现出明显差异。另外，它们都呈现周期性变化，极限环3 和极限环2 对应的无量纲周期分别约为6.9 和7.2。这说明位于极限环3 上的颗粒虽然经历的周长更长，但是其运动速度更大（图5（a）），可以用更短的时间完成一个周期的运动。鉴于极限环3 和极限环2 的运动特征高度相似，接下来我们仅选取了其中一种进行详细分析。

图5 极限环2 和极限环3 上，颗粒运动的合速度和角速度在一段时间内的变化规律Fig. 5 Time evolutions of the particle’s velocity magnitude and the corresponding angular velocity for the limit cycle 2 and limit cycle 3

为了进一步探索颗粒的运动机理，我们以初始网格位置（77，50，48）为例进行分析。图6（a）展示了角速度随时间的演化，角速度由U/S标准化，此时 Ωx和Ωy等 于0， Ωz经过一段时间的不规则运动后周期性地变化。这意味着颗粒经过一段时间的不规则运动后到达极限环，最终在极限环3 上运动。此外颗粒动能E随时间的演化显示在图6（b），这里E=0.5mup2+0.5IΩ2，是平动能和转动能之和；动能E由 0.5mu2max无量纲化，其中umax代表颗粒最大速度。图6（b）中虚线框的局部放大图展示了一个周期上的变化，可以发现，动能E存在最小值，对应图7的d点。

图6 颗粒角速度和颗粒动能随时间的演化Fig. 6 Time evolutions of the angular velocity and kinetic energy for the particle

图7 为颗粒受力与颗粒速度的夹角，以及颗粒受力在其运动轨迹和其运动轨迹垂直方向上的分量随时间的演化曲线。图7（a）给出了极限环上的一个周期，这里a、b、c、d、e、f分别代表颗粒受力与速度垂直时的位置。这里指出当 θ大于90°时，颗粒做减速运动；反之颗粒做加速运动。紧接着我们把颗粒受到的总力F分解成沿着运动轨迹方向的分量Fup和垂直于该运动方向的分量Fuv，力Fup关系到颗粒速度大小的改变，而力Fuv控制着颗粒运动的方向，结果展示在图7（b）中。可以发现Fup明 显小于Fuv，而且颗粒在c-d和d-e阶段经历了更大的减速和加速运动。

图7 颗粒受力F 与速度之间的夹角以及相应颗粒在运动方向上的受力分量Fup 和在垂直于运动方向上的受力分量Fuv 随时间的变化Fig. 7 The angle between particle force F and velocity and the force component Fup in the direction of motion and the force component Fuv in the direction perpendicular to the direction of motion in function of time

为了更加清晰地观察颗粒在极限环上的运动规律，几个典型的位置画在图8 中，这里实线代表加速过程，虚线代表减速过程。从图中可以明显地看到颗粒交替的加速减速运动，这与图7 中的结果保持一致。球形颗粒上黑色实线用来表征其转动特征。由于上壁面的剪切作用，颗粒在方腔中平动的同时伴随着顺时针地转动。

图8 z = 0.5 平面上，颗粒在极限环3 上一个周期的运动Fig. 8 Orbiting and rotating of a particle at different positions on the limit cycle 3 in plane z = 0.5

此外我们还测试初始网格位置为（70，50，36）的情形，其他参数保持不变。图9 为Re= 1 000 时，颗粒位置在x-y平面和展向随时间的演化曲线。我们发现颗粒交替地向上壁面运动或远离上壁面。同时一个三维的颗粒轨迹图展示在图9（d），有趣的是，我们观察到球形颗粒只在方腔的一侧运动。颗粒首先在横向螺旋运动的伴随下向内部迁移，然后逐渐向外围移动。轨迹呈现两个明显的特点：x-y平面上的快速迁移和z方向上相对较慢的横向运动，这和Tsorng等[29-30]的结果相似。他们实验观察到球形颗粒在Re= 860 时被限制在方腔的一侧，然而颗粒在高雷诺数Re= 3 200 时在两侧往复运动。

图9 x-y 平面和展向，颗粒位置随时间的演化和三维视图Fig. 9 Evolution of x-y plane and spanwise particle positions over time and 3D views

3.2 雷诺数的影响

这一部分研究雷诺数对固体颗粒运动的影响。初始网格位置选取（74，50，48），颗粒半径设置为6 个网格。和之前一致，在形成单相稳定流场后球形颗粒加入，并允许其自由运动。对于雷诺数Re= 100、400、1 000，相 应 的 斯 托 克 斯 数 分 别 为St= 0.087、0.347、0.868。对于展向弱受限的情形P0，图10（a～c）描述了不同雷诺数下的极限环轨迹。其中红色实心圆点代表颗粒初始位置，实线代表极限环轨迹。背景虚线表示z= 0.5 平面上单相时的流线。极限环上颗粒速度的大小用颜色表示，颗粒速度由相应的上壁面速度U无量纲化。如图10 所示，随着雷诺数的增大，极限环逐渐靠近外侧，逐渐变大。我们知道单相时，随着流体惯性的增加，主涡中心向方腔中心移动。极限环轨迹也有往方腔外围移动的趋势。在不同雷诺数下，方腔内流场结构存在较大差异，流固相互作用影响了极限环轨迹。

图10 不同雷诺数下的极限环轨迹Fig. 10 Limit cycles at different Re

此外，我们发现颗粒在极限环轨道上交替的加速和减速运动。当Re= 100 时，颗粒最大速度出现在极限环的顶部，靠近移动的上壁面，最小速度出现在左侧。与低雷诺数结果不同的是，Re= 400 和1 000 时的最小速度出现在区域对角线附近的极限环的右上角。如图7 所示，在d点之前，颗粒的减速周期相对较大且较长。

3.3 颗粒大小的影响

这一部分我们主要考虑了颗粒大小的影响。这里初始网格位置选取（74，50，48），当雷诺数取Re=1 000，颗粒半径R= 6、9、12，对应的斯托克斯数分别为St= 0.868、1.953、3.472，其他参数保持不变。对于展向弱受限的情形，雷诺数Re= 1 000 时，不同尺寸的颗粒对应的极限环轨迹展示在图11（a）。其中红色实心圆点表示初始位置，从内到外分别代表R=6、9、12 时的极限环。背景虚线表示z= 0.5 平面上单相时的流线。极限环上颗粒速度的大小用颜色表示，颗粒速度由相应的上壁面速度U无量纲化。由于初始位置位于平面z= 0.5，颗粒仅在该平面内运动，且在稳定轨道上交替地加速和减速运动。此外，随着St的增加，极限环向外围移动。在这里，St描述了方腔流动中悬浮固体颗粒的行为。如果St越大，则表示颗粒受惯性的影响较大，继续沿着原来的轨迹运动；反之则沿着相应的流体流线运动。如图11（a）所示，极限环轨迹随着颗粒大小的变化而变化，大致呈椭圆形。当颗粒运动到极限环的上部时，R= 12 时的St较大，导致与流体流线偏离。由于速度较低，它在底部几乎沿流线移动。

图11 不同雷诺数下，极限环大小随颗粒大小的变化规律Fig. 11 Limit cycles for particles of different sizes at different Re

此外，我们还模拟了Re相对较低的情形。这里采用Re= 100，R= 6、12，对应斯托克斯数分别为St=0.087、0.347。其 他 参 数 与Re= 1 000 的 情 况 相 同。图11（b）给出了不同颗粒半径下的极限环，其中红色实心点为初始位置，外侧实线和内侧实线分别为R=6 和12 的极限环。与高雷诺数结果不同的是，大的颗粒向靠近涡中心的轨道迁移，而小的颗粒向外极限环迁移，这与Khojah 等[15]的结论定性上一致。根据颗粒或细胞的尺寸特性，通过改变相应的流动惯性或腔体尺寸可以在微流控平台上实现惯性分离[15-17]。

4 展向强受限的情形

4.1 初始位置的影响

这一部分我们主要考虑了展向在强受限情况下初始位置对颗粒最终运动的影响。对于P1 情况，雷诺数Re= 1 000 和半径R= 6 的参数下，对应的斯托克斯数St= 0.868。颗粒在展向z= 0.5 平面上的不同初始位置释放，我们模拟了25 组算例，且颗粒的初始位置如图12 中空心菱形所示，黑色实线表示最终的稳定轨迹。可以发现，颗粒的初始位置对最终的极限环轨迹没有影响。这与Hu[33,48]的结果一致，他在文章中指出，二维方腔中颗粒的极限环轨迹与长椭球颗粒的初始朝向和初始位置无关。

图12 P1 情况下，z = 0.5 平面上颗粒从不同位置释放下的相图Fig. 12 Phase diagram for a spherical particle inside lid-driven cavity flow in plane z = 0.5 for cases P1

4.2 雷诺数的影响

这一部分研究雷诺数对固体颗粒运动的影响。初始网格位置选取（74，50，48），颗粒半径设置为6 个网格。对于雷诺数Re= 100、400、1 000，相应的斯托克斯数St= 0.087、0.347、0.868。对于展向强受限的情形P1，图13 给出了展向z= 0.5 平面上不同雷诺数下的极限环轨迹，其中红色实心圆点代表颗粒初始位置。发现雷诺数对极限环的拓扑结构有明显的影响。随着雷诺数的增加，除了方腔的左上区域外，极限环轨道有向外迁移的趋势。Hu[48]在文章中指出，随着雷诺数的增加，左上角的涡旋发展，球形颗粒在二维方腔中的极限环被推向方腔右下角。然而由于展向强受限（P1），极限环的拓扑结构与P0 情形下的结果有很大不同。

图13 P1 情况下，Re = 100、400、1 000 的极限环轨迹Fig. 13 Limit cycles at Re = 100、400 and 1 000 for the cases P1

4.3 颗粒大小的影响

最后我们主要考虑了颗粒大小对最终运动轨迹的影响。这里初始网格位置选取（74，50，48），当雷诺数取Re= 1 000，颗粒半径R= 6、12 对应的斯托克斯数St分别为0.868、3.472，其他参数保持不变。对于展向强受限的情形P1，雷诺数Re= 1 000 时，不同尺寸的颗粒对应的极限环轨迹展示在图14（a）。其中红色实心圆点表示初始位置。可以发现小颗粒的极限环在外围，大颗粒的极限环在内侧。这种现象和展向不受限情形P0 下Re= 1 000 的结果恰好相反。

此外，我们还模拟了Re相对较低的情形。这里采 用Re= 100，R= 6、12，对 应St数 分 别 为0.087、0.347。其他参数与Re= 1 000 的情况相同。图14（b）给出了颗粒半径R= 6 和12 时的极限环，其中红色实心点为初始位置，外侧实线和内侧实线分别为R=6 和12 的极限环，虚线为颗粒到达极限环之前的运动轨迹。与高雷诺数Re= 1 000 的结果一致，大的颗粒向靠近涡中心的轨道迁移，而小的颗粒向外极限环迁移。这与展向弱受限情形P0 下的不同，可能是由于展向壁面的约束导致的。

图14 不同雷诺数下，极限环大小随颗粒大小的变化规律Fig. 14 The variation of limit cycle size with particle size at different Re

5 结 论

本文采用多松弛格子Boltzmann 方法模拟了单个中性球形颗粒在三维方腔中的运动。首先通过单相顶盖驱动流和网格无关性测试验证了程序的准确性。其次主要考虑了方腔展向弱受限（对称边界条件）和强受限（固壁无滑移边界条件）两种情况，着重研究了初始位置、颗粒大小以及雷诺数的影响，得到以下结论：

对于展向弱受限的情形P0，发现颗粒的初始位置强烈地影响着最终的运动轨迹。根据展向z=0.5 平面的相图，可以得到其被划分为三个区域：外层稳定区、内层稳定区以及涡中心区域。通过对颗粒受力的分解，解释了其在极限环上运动的机理。此外，还详细介绍了球形颗粒在极限环上的顺时针旋转运动，呈现交替性地加速和减速运动。也模拟了初始位置位于展向z/S= 0.375 平面的情形，观察到球形颗粒只在方腔的一侧运动，伴随着在展向z方向上相对较慢的运动。随着雷诺数的增加，颗粒逐渐向外围靠近，不断旋转达到最优路径，这里最优路径是指相应参数下对应的极限环轨迹，运动路径几乎遵循相应的流线。对于选定的初始位置，我们观察到大颗粒的极限环在高雷诺数时向外迁移，而在低雷诺数时更靠近涡中心。

对于展向强受限的情形P1，极限环轨迹与颗粒的初始位置无关。随着雷诺数的增加，除方腔左上角外，极限环轨迹有向外迁移的趋势。在Re= 1 000 和100 两种情况下，随着颗粒尺寸的增大，极限环都会相应缩小。

通过与参考文献的对比分析，在一定程度上验证了本文结果，希望目前的结果对管道和空腔流中颗粒和细胞的过滤和分离提供一些参考。