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]给出了颗粒团在二维方腔中运动的数值结果，其中雷诺数为1000，采用双向耦合方法处理固体点颗粒与流体的相互作用。结果表明颗粒倾向于向壁面迁移，且斯托克斯数越大，迁移行为越强烈。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模型对应的颗粒分布函数的演化方程为：

$ \begin{split}& {\boldsymbol{f}}({\boldsymbol{x}} + {\boldsymbol{c_i}}{\delta _t},t + {\delta _t}) - {\boldsymbol{f}}({\boldsymbol{x}},t) =\\&\qquad - {{\boldsymbol{M}}^{ - 1}}{\boldsymbol{s}}\left[ {{\boldsymbol{m}}({\boldsymbol{x}},t) - {{\boldsymbol{m}}^{{\rm{eq}}}}({\boldsymbol{x}},t)} \right] + {{\boldsymbol{M}}^{ - 1}}{\boldsymbol{\psi}} \end{split}$

(1)

在这里， ${\boldsymbol{f}}({\boldsymbol{x}},t)$ 是在位置 ${\boldsymbol{x}}$ 和时间 $ t $ 上的分布函数； ${\boldsymbol{ M}}$ 是转换矩阵，其把 ${\boldsymbol{f}}$ 转换成矩向量 ${\boldsymbol{m}} = {\boldsymbol{Mf}}$ ； ${{\boldsymbol{m}}^{{\rm{eq}}}}$ 和 ${\boldsymbol{\psi}}$ 分别是平衡态矩向量和外力项向量。对角碰撞矩阵s为：

$ \begin{split} {\boldsymbol{s}} =& {{\rm{dig}}} (0,{s_2},{s_2},{s_2},{s_4},{s_5},{s_5},{s_7},{s_7},{s_9},{s_{10}},{s_{10}},{s_{10}}, \\& {s_{13}},{s_{13}},{s_{13}},{s_{16}},{s_{17}},{s_{18}},{s_{18}},{s_{20}},{s_{20}},{s_{20}},{s_{23}},{s_{23}},{s_{23}},{s_{26}}) \end{split} $

(2)

这里松弛因子分别设置为： ${s_2} = {s_5} = {s_7} = {s_9} = 1/\tau $ 、 ${s}_{4} = {s}_{17} = 1.4、{s}_{10} = 1.5、{s}_{13} = {s}_{16} = 1.74、{s}_{18} = {s}_{20} = 1.8$ 、 $ {s_{23}} = 1.83 $ 、 $ {s_{26}} = 1.61 $ ，松弛因子 $\tau $ 和流体运动学黏性有关： $\upsilon = c_{s}^2(\tau - 0.5){\delta _t}$ 。相应的权系数分别为 ${\omega }_{0} = 8/27、{\omega }_{1-6} = 2/27、{\omega }_{7-18} = 1/54、{\omega }_{19-26} =$ $1/216$ 。相应的离散速度 ${{\boldsymbol{c}}_i}$ 为：

$ {{\boldsymbol{c}}_i} = \left\{{\begin{array}{*{20}{l}} {\left( {0,0,0} \right)c,}&{i = 0} \\ {\left( { \pm 1,0,0} \right)c,\left( {0, \pm 1,0} \right)c,\left( {0,0, \pm 1} \right)c,}&{i = 1 ～ 6} \\ {\left( { \pm 1, \pm 1,0} \right)c,\left( { \pm 1,0, \pm 1} \right)c,\left( {0, \pm 1, \pm 1} \right)c,}&{i = 7 ～ 18} \\ {\left( { \pm 1, \pm 1, \pm 1} \right)c,}&{i = 19 ～ 26} \end{array}} \right. $

(3)

在这里，i表示离散速度的方向，格子速度 $ c = {{{\delta _x}} \mathord{\left/ {\vphantom {{{\delta _x}} {{\delta _t}}}} \right. } {{\delta _t}}} $ ， $ {\delta _x} $ 和 $ {\delta _t} $ 分别为格子单位和时间单位，格子声速为 $ {c_s} =c/ {\sqrt 3 } $ 。在目前的模拟过程中 $ {\delta _x} = {\delta _t} = 1 $ 。局部的密度波动 $\delta \rho $ 和流体速度 ${\boldsymbol{u}} $ 通过以下公式计算：

$ \delta \rho = \sum\limits_{{i} = 0}^{26} {{f_{i} }} ,\;\;\;{\rho _0}{\boldsymbol{u}} = \sum\limits_{{\alpha} = 0}^{26} {{c_{i} }} {f_{i} } + \frac{{{\boldsymbol{F}}{\delta _t}}}{2} $

(4)

其中 $\boldsymbol{F}$ 为作用在颗粒上的合力。通过Chapman-Enskog分析，目前格子Boltzmann模型可以准确地恢复到不可压的Navier-Stokes方程^[40]。

1.2 颗粒的运动

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

$ {m_{\rm{p}}}\frac{{{\rm{d}}{{\boldsymbol{u}}_{\rm{p}}}(t)}}{{{\rm{d}}t}} = {\boldsymbol{F}}(t) $

(5)

$ {\boldsymbol{I}}\frac{{{\rm{d}}{\boldsymbol{\varOmega}} (t)}}{{{\rm{d}}t}} = {\boldsymbol{T}}(t) $

(6)

这里，I是惯性张量；颗粒质量为 ${m_{\rm{p}}} = 4{\rho _{\rm{p}}}\text{π} {{R}^3}/3$ ，其中 $ \;{\rho _{{\rm{p}}}} $ 和R分别为球形颗粒密度和半径；T、 $ {{\boldsymbol{u}}_{\rm{p}}} $ 和 $ {\boldsymbol{\varOmega}} $ 分别表示作用在颗粒上的力矩、速度和角速度。

Bouzidi等^[42]提出的插值反弹格式被用来处理颗粒表面，详细的形式如文章所述。本文仅考虑了中性悬浮的球形颗粒，颗粒所受到的力F包括两部分：水动力 ${\boldsymbol{F}_{\rm{H}}}$ 和润滑力 ${\boldsymbol{F}_{{\rm{lub}}}}$ 。动量交换法^[43-46]能够方便地求出水动力和力矩，Wen等^[46]提出的满足伽利略不变性的格式被采用。具体的形式如下所示：

$ {{\boldsymbol{F}}_{\rm{H}}} = \sum\limits_{{x_{\rm{f}}}} {\sum\limits_{{i_{{\rm{bl}}}}} \Big[ } {\tilde f_i}({{\boldsymbol{x}}_{\rm{f}}},t)({{\boldsymbol{c}}_i} - {{\boldsymbol{u}}_{\rm{w}}}) - {f_{\bar i}}({{\boldsymbol{x}}_{\rm{f}}},t + {\delta _t})({{\boldsymbol{c}}_{\bar i}} - {{\boldsymbol{u}}_{\rm{w}}})\Big] $

(7)

$ {\boldsymbol{T}} = ({{\boldsymbol{x}}_{\rm{w}}} - {{\boldsymbol{Y}}_{\rm{c}}})\sum\limits_{{{x}_{\rm{f}}}} {\sum\limits_{{i_{{\rm{bl}}}}} \Big[ } {\tilde f_i}({\boldsymbol{x}_{\rm{f}}},t)({{\boldsymbol{c}}_i} - {{\boldsymbol{u}}_{\rm{w}}}) - {f_{\bar i}}({{\boldsymbol{x}}_{\rm{f}}},t + {\delta _t})({{\boldsymbol{c}}_{\bar i}} - {{\boldsymbol{u}}_{\rm{w}}})\Big] $

(8)

这里， ${{\boldsymbol{Y}}_{\rm{c}}}$ 是球形颗粒的质心位置；下脚标f代表流体，w代表壁面，bl代表边界， ${\bar i} $ 代表i的反方向；边界点上的速度 ${{\boldsymbol{u}}_{\rm{w}}}$ 由颗粒的平动速度和角速度控制：

$ {{\boldsymbol{u}}_{\rm{w}}} = {{\boldsymbol{u}}_{\rm{p}}} + {\boldsymbol{\varOmega}} \times ({{\boldsymbol{x}}_{\rm{w}}} - {{\boldsymbol{Y}}_{\rm{c}}}) $

(9)

当球形颗粒非常靠近壁面的时候，润滑力 ${{\boldsymbol{F}}_{{\rm{lub}}}}$ 被引入，Feng和Michaelides^[47]提出的斥力模型被使用，形式如下：

$ {{\boldsymbol{F}}}_{{\rm{lub}}} = \left\{\begin{array}{l}0,\;\;\;\;\;\;\;\;\,\quad\qquad\qquad\qquad \Vert {{\boldsymbol{x}}}_{ij}\Vert > \zeta \\ \dfrac{{c}_{{\rm{f}}}}{{\varepsilon }_{{\rm{w}}}}{\Bigg(\dfrac{\Vert {{\boldsymbol{x}}}_{ij}\Vert -\zeta }{\zeta }\Bigg)}^{2}\dfrac{{{\boldsymbol{x}}}_{ij}}{\Vert {{\boldsymbol{x}}}_{ij}\Vert },\;\;\;\;\;\Vert {{\boldsymbol{x}}}_{ij}\Vert \leqslant \zeta \end{array}\right. $

(10)

其中： ${c_{\rm{f}}}$ 表示流体中特征力的大小； ${\varepsilon _{\rm{w}}}$ 是一个可调参数，参考Huang等^[5]取0.0025； $\Vert {{\boldsymbol{x}}}_{ij}\Vert$ 表示颗粒与壁面之间的距离；阈值 $ \zeta $ 小于3个网格时引入润滑力。

2 问题设置和数值验证 2.1 问题设置

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

$ Re = \frac{{US}}{\nu } $

(11)

$ S t = \frac{{{\tau _{\rm{p}}}}}{{{\tau _{\rm{f}}}}} = \frac{{\dfrac{{\rho _{\rm{p}}}d_{\rm{p}}^2}{18{\rho _{\rm{f}}}\nu }}}{\dfrac{S}{U}} = \frac{1}{{18}}\frac{{{\rho _{\rm{p}}}}}{{{\rho _{\rm{f}}}}}{\left( {\frac{{{{d}_{\rm{p}}}}}{S}} \right)^2}Re $

(12)

其中 ${\tau _{\rm{p}}} = {\rho _{\rm{p}}}d_{\rm{p}}^2/(18{\rho _{\rm{f}}}\nu )$ 和 ${\tau _{\rm{f}}} = S/U$ 分别是直径为 ${d_{\rm{p}}}$ 、密度为 $\;{\rho _{\rm{p}}}$ 的颗粒的弛豫时间和特征时间尺度。

2.2 数值验证

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

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

3 展向弱受限的情形 3.1 初始位置的影响

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

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

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

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

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

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

3.2 雷诺数的影响

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

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

3.3 颗粒大小的影响

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

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

4 展向强受限的情形 4.1 初始位置的影响

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

4.2 雷诺数的影响

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

4.3 颗粒大小的影响

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

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

5 结　论

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

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

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

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