基于相对运动模型的多颗粒与颗粒链对比研究

Comparative Study of Multi-Particle and Particle Chains Based on Relative Motion Model

摘要: 本文基于单颗粒相对运动模型，扩展构建出多颗粒相对运动模型和颗粒链相对运动模型。针对计算结果、计算效率和计算范围进行了对比分析，得出结论：多颗粒相对运动模型和颗粒链相对运动模型计算结果均与实验结果吻合，具有较高的计算精度；多颗粒相对运动模型的网格数量远高于颗粒链相对运动模型，计算效率较低；对于颗粒数量，多颗粒相对运动模型具有一定的局限性，颗粒链相对运动模型的使用范围相对更宽泛。本文旨在为研究者在不同工况下选择合适的模型提供参考。

Abstract: Based on the single-particle relative motion model, this study develops extended multi-particle and particle chain relative motion models. A comparative analysis was conducted on computational results, efficiency, and scope. The results indicate that both the multi-particle and particle chain relative motion models exhibit high computational accuracy, with predictions in good agreement with experimental data. However, the multi-particle relative motion model requires a significantly larger number of meshes compared to the particle chain relative motion model, resulting in lower computational efficiency. Furthermore, the multi-particle relative motion model demonstrates limitations in handling large particle populations, whereas the particle chain relative motion model offers broader applicability. This study aims to provide references for researchers in selecting appropriate models under varying operational conditions.

1. 引言

微通道多颗粒直线聚集是微流控技术中的核心研究之一，其关注微米尺度通道内颗粒在流体惯性效应作用下的定向迁移行为。这一现象在生物医学检测、细胞分选、药物输送及微纳颗粒操控等领域展现出重要的应用价值。微通道内颗粒的行为受通道几何形状、流场特性、颗粒间相互作用等多重因素影响。

早期研究主要集中于单颗粒的惯性聚焦现象。Segré和Silberberg [1]于20世纪60年代在实验中发现，在层流条件下，通道中刚性颗粒会迁移至距通道中心约60%半径的平衡位置(即Segré-Silberberg效应)，这一发现颠覆了传统斯托克斯流理论中颗粒仅随流线运动的认知，揭示了惯性力在微尺度流动中的不可忽视性。Di Carlo等[2]系统研究了通道截面形状(矩形、圆形)对颗粒平衡位置的影响，提出了“惯性聚集”的概念，并构建出“相对运动模型”，将原本非定常的颗粒聚集问题转化成为准定常问题，大大提高了计算效率。王企鲲课题组[3]-[5]继承并进一步发展了“单颗粒相对运动模型”，通过“试凑法”确定了通道截面不同位置上颗粒处于稳定状态下的平移速度Up以及旋转角速度ωp，进而推导出颗粒所受到的惯性升力的横向分布特征，并深入分析了构成惯性升力的各种影响因素。

微通道颗粒惯性迁移的研究已从单颗粒体系拓展至多颗粒复杂系统。Humphry等[6]观察到沿流动方向排列的稳定直线颗粒链，并将这种现象归因于与流体惯性相关的流体动力相互作用。随着雷诺数的增加，颗粒链中自组织颗粒的比例先增加后下降，其中最大值取决于颗粒大小和管道尺寸。Di Carlo等[7]通过实验发现了类似自排序现象的颗粒链，并验证了稀释后的全血细胞在直线微通道中具有与刚性颗粒相似的行为，其结论表明刚性颗粒直线排列的行为可以为流式细胞术等实际生物医学应用提供参考。Kahkeshani等[8]通过实验发现，颗粒迁移到平衡位置时会形成单线颗粒链。随着颗粒雷诺数从2.8增加到8.3，相邻球形颗粒之间的稳定间距逐渐减小。Gupta等[9]通过实验分析了单线颗粒链在牛顿流体中迁移时的稳定性，结论表明：稳定的单线颗粒链是由有限数量的颗粒形成的，最下游的颗粒会逐渐脱离颗粒链，同时雷诺数和堵塞比会影响单线颗粒链的形成。Gao等[10]通过实验研究发现颗粒只有达到平衡位置时，颗粒链才会形成，并分析了颗粒间的稳定间距，结果表明颗粒稳定间距随着颗粒雷诺数的增加而减小，并且与颗粒大小和浓度无关。Liu等[11]使用浸没边界格子玻尔兹曼方法研究了方管内直线颗粒链受颗粒雷诺数、长度分数和颗粒大小对颗粒运动的影响，其中颗粒雷诺数的增大会导致颗粒聚集位置更靠近墙壁，且引起颗粒稳定间距的更大波动。

与单颗粒相比，多颗粒能够更真实地反映颗粒之间的相互作用和集体行为。众多学者通过实验研究、理论推导和数值模拟等多种手段，对牛顿流体中多颗粒和颗粒链的形成机制、稳定性以及流动轨迹等方面进行了广泛探讨，但对于直线多颗粒和颗粒链在不同间距下的力学特性研究还较为缺乏。本文基于单颗粒相对运动模型，通过FLUENT商用软件扩展构建出多颗粒相对运动模型和颗粒链相对运动模型，准确描述颗粒在不同间距下的力学分布情况，并对两种模型的计算精度、计算效率和使用范围进行分析对比，可供读者选用合适的模型。

2. 控制方程与计算方法

2.1. 控制方程

通过FLUENT商用软件进行有限体积法三维数值模拟，主要涉及包含颗粒的管内流动，流动介质为水，控制方程是由不可压缩连续性方程和动量方程构成：

$\nabla \cdot W=0$ (1)

$\left(W\cdot \nabla \right)W=-\frac{1}{\rho }\nabla p+\upsilon {\nabla }^{2}W$ (2)

其中，W为平移坐标系中流体的相对速度，p为流体的压强，ρ为水的密度，取1000 kg/m³，υ为运动粘度。

2.2. 计算方法

在FLUENT中采用双精度求解器，压力和速度的耦合采用COUPLED算法，压力方程的离散采用标准格式，动量方程采用QUICK格式离散。两种模型均保持相同计算方法与格式。

3. 模型提出与设置

3.1. 多颗粒相对运动模型

在许多实际应用场景和自然现象中，颗粒并非孤立存在，而是以多颗粒的形式相互作用和影响，呈现出更为复杂和丰富的集体行为特征。因此，基于单颗粒相对运动模型，提出了多颗粒相对运动模型，可用于研究多个颗粒以链状直线排列形式相互作用。多颗粒相对运动模型是在单颗粒相对运动模型的基础上，增加多个颗粒，坐标轴放置于居中颗粒的中心，可将多颗粒相对运动模型的非定常问题转化为定常问题。将多个颗粒移动至颗粒稳定聚集位置，通过单颗粒相对运动模型的“试凑法”获得单颗粒升力曲线[3]，将稳定聚集位置处的平移速度U_ep与旋转角速度ω_ep，代入多颗粒相对运动模型中的每个颗粒，逐渐改变间距，计算获得不同间距下颗粒的力学特性，物理模型如图1所示。

Figure 1. Physical model of multi-particle relative model

图1. 多颗粒相对运动物理模型

边界条件设置与单颗粒相对运动模型相同：

1) 进、出口边界条件：在相对运动模型中，流体相对速度的方向与坐标系中轴向x正方向相反，流体由通道出口面(Outlet)流入，从进口面(Inlet)流出，以此形成相对运动。因此，将进口面设置为零背压出口(Pressure-Outlet)，出口面设置为速度入口(Velocity-Inlet)；

2) 壁面条件(Wall)：通道壁面设置为无滑移移动壁面，并以稳定位置处的速度U_ep沿与主流方向相反的方向移动；

3) 颗粒表面条件(Ball)：相对运动模型中，颗粒的平移速度转移给了壁面，只存在稳定位置处的旋转角速度ω_ep。

3.2. 颗粒链相对运动模型

微通道中的稳定直线颗粒链在宏观上具有周期性特征。通道结构和充分发展的流体流动在一定程度上是重复的，通过使用周期性边界条件(Periodic boundary condition)，便能模拟这种实际的周期性情况，使计算结果更符合真实物理系统的行为。因此，本节将“多颗粒相对运动模型”的基本思想与周期性边界条件相结合，构建了“颗粒链相对运动模型”，即相对运动模型坐标轴原点放在“颗粒链相对运动模型”中的单一周期颗粒上，在一个周期单元进出口实施平移周期性边界条件[4]，通过该单一周期颗粒可反映无限长颗粒链居中颗粒的力学特性，物理模型如图2所示。颗粒链相对运动模型的边界条件与多颗粒相对运动模型的不同之处在于实际计算域的进、出口处采用了平移周期性边界条件并且只考虑单一周期颗粒的力学特性，其余条件均保持相同，周期性边界条件具体设置如下：首先在FLUENT控制台输入TUI命令：mesh/modify-zones/make-periodic；依次输入出口、进口面对应的边界ID；随后输入“no”指令选择平移周期，否则将选择为旋转周期边界条件；输入两次“Enter”指令；最后在边界条件面板中输入流体相对速度对应的质量流量，方向为x负方向，即可创建完成平移周期性边界条件。

Figure 2. Physical model of particle chain relative model

图2. 颗粒链相对运动物理模型

使用周期边界条件后，唯一周期颗粒处于充分发展段流动稳定。颗粒链相对运动模型可以将颗粒间距∆x巧妙转换为实际计算域长度L，如图3所示：相邻两颗粒中心间距∆x，其长度为相邻两颗粒表面间距x和两个颗粒半径(0.5 a)的总和；将管道进出口移动至相邻颗粒表面间距的中心，实际计算域长度L宏观上为管道进出口之间的距离，其数值大小为颗粒直径a、管道进口到颗粒表面的距离0.5 x和出口到颗粒表面的距离0.5 x的总和，由此可得相邻两颗粒中心间距∆x和实际计算域长度L的数值大小均为：x + a，因此可以通过改变实际计算域长度L从而改变颗粒间距∆x。

Figure 3. Correspondence graph between particle spacing and computational domain length

图3. 颗粒间距与计算域长度对应图

多颗粒相对运动模型考虑到有限数量之间的力学表现，尤其居首、居尾颗粒。而颗粒链相对运动模型重点关注颗粒链中周期颗粒的力学特性，通过唯一周期颗粒即可反映出长颗粒链中整体颗粒的力学特性。在数值模拟计算后如何判断颗粒间距是否达到稳定状态是首要考虑的问题。在层流情况下，随着颗粒间距逐渐增大，颗粒间的相互作用力以及颗粒受到的流体粘性阻力、升力等流体动力作用力相互平衡抵消，达到了相对数值零点，系统总动能达到了最小状态[12]，此时可认为颗粒间距达到稳定。为了方便下文讨论和分析，本文定义以下无量纲参数：

1) 颗粒雷诺数：

$R{e}_p={\left(a/D\right)}^{2}Re=\frac{\rho V{a}^2}{\eta D}$ (3)

其中，V为管内流体平均流速，η为水的动力粘度，取0.001 Pa∙s，a为颗粒直径，D为管道水力直径。

2) 阻力系数：

$C_D=\frac{F_D}{\frac{1}{2}\rho V^{2}A}$ (4)

其中，F_D为颗粒所受到的阻力，A为颗粒迎风面积，$A=\pi \cdot {\left(a/2\right)}^2$ 。

3) 升力系数：

$C_{F_L}=\frac{F_L}{\rho V^2{a}^4/D^2}$ (5)

其中，F_L为颗粒所受到的升力。

4. 结果与讨论

4.1. 计算结果对比

为了验证和对比多颗粒相对运动模型和颗粒链相对运动模型的计算结果，本节分别采用两种模型对文献[8]的双工况(Re_p = 2.8、Re_p = 8.3)进行数值模拟，通过颗粒在不同间距下的力学特性(阻力、升力)表现，从而获得颗粒临界稳定间距，模拟结果如下图4、图5所示。通过模拟结果发现颗粒分别在无量纲间距∆x/a为5、2.5左右时：每个颗粒的阻力和升力同时达到相对数值零点。颗粒链相对运动模型的唯一周期颗粒计算结果与多颗粒相对运动模型的居中颗粒基本吻合。

Re_p = 2.8：由图4可知，当∆x/a从1增加至5时，处于不稳定状态，五个颗粒的阻力大小逐渐减小，阻力方向始终为负，居中三个颗粒阻力基本相同且大于首尾两个颗粒。但升力变化与阻力并不相同，居首颗粒的升力方向为负，即指向管道中心方向；居尾颗粒的升力方向为正，即指向管道壁面方向。首尾两颗粒的升力大小逐渐单调减小。而居中的三个颗粒，升力方向先为正，升力大小逐渐减小为零值；在∆x/a = 2.5时，方向开始转为负，升力大小开始增大，∆x/a = 3时，升力大小开始减小，∆x/a = 5时，升力减小为零。当∆x/a从5增加至12，颗粒达到稳定状态，所有颗粒的阻力和升力方向不变，大小为相对零值且保持不变。

(a) (b)

Figure 4. Multi-particle mechanical characteristic distribution at Re_p = 2.8: (a) Drag; (b) Lift

图4. Re_p = 2.8时多颗粒力学特性分布：(a) 阻力；(b) 升力

(a) (b)

Figure 5. Multi-particle mechanical characteristic distribution at Re_p = 8.3: (a) Drag; (b) Lift

图5. Re_p = 8.3时多颗粒力学特性分布：(a) 阻力；(b) 升力

Re_p = 8.3：由图5可知，当∆x/a从1增加至2.5时，每个颗粒的阻力和升力变化与Re_p = 2.5的不稳定状态相同，不同之处为居首颗粒的升力方向在∆x/a = 1时为正，间距增大后转为负方向，升力大小先减小后增大。当∆x/a从2.5增加至6时，阻力和升力均保持不变，大小达到相对零值。

因此，可认为∆x/a = 5、2.5分别是Re_p = 2.8、8.3时多颗粒的临界稳定间距，这与Kahkeshani等的实验结果[8]基本一致。∆x/a小于临界稳定间距时，多颗粒间相互影响，处于不稳定状态。∆x/a大于临界稳定间距时，所有颗粒的阻力和升力保持不变，始终为数值零点，说明颗粒与颗粒逐渐远离直至独立，互不影响，达到了单颗粒稳定状态。综上所述，可认为通过多颗粒相对运动模型和颗粒链相对运动模型获得临界稳定间距以及各个颗粒在不同间距下的力学特性是可靠的。

4.2. 计算效率对比

网格数量直接影响着数值模拟计算效率。本文通过ICEM软件对模型进行结构化网格划分。鉴于在计算域中的数值计算主要是为了获得颗粒的受力情况，图6为颗粒网格划分示意图，为了提高计算的稳定性和精度，将颗粒及其周围的流体域单独切割出来，切成一个立方体包裹着颗粒，然后对该立方体域由内向外逐渐变疏进行网格加密处理。颗粒材料强度为刚性，不考虑变形，因此可去除颗粒内部网格。

(a) (b)

Figure 6. Schematic diagrams of the grids: (a) Grids around particle; (b) Girds around particle surface

图6. 网格示意图：(a) 颗粒附近网格；(b) 颗粒表面网格

为了保证计算的准确性，首先对网格无关性进行了验证分析，如图7所示。分别对颗粒链相对运动模型和多颗粒相对运动模型的三种颗粒间距(∆x/a = 2、4、6)进行网格划分验证。由图7可知，当颗粒链相对运动模型分别在30万、34万、40万左右时，多颗粒相对运动模型网格数量分别在140万、180万、220万左右时，颗粒的阻力都趋于稳定。为了在满足计算精度的基础上尽可能地节约计算时间，最终确定颗粒链间距∆x/a = 2、4、6，颗粒链相对运动模型网格数量控制在30万、34万、40万左右，多颗粒相对运动模型网格数量控制在140万、180万、220万左右。

图8为不同间距下两种模型网格数量对比图。由此图可知，不同间距下颗粒链相对运动模型的网格数量远少于多颗粒相对运动模型。随间距增大，多颗粒相对运动模型网格数量急剧增加，而颗粒链相对运动模型网格数增加相对平缓。其内在原因：多颗粒相对运动模型不仅要考虑间距增大使得计算域加长与不同工况下所需的颗粒数量，还要考虑到进出口对居首、居尾颗粒的影响，从而导致网格数量庞大。颗粒链相对运动模型只需要对单一周期颗粒周围进行加密处理即可，且周期性边界条件可忽略进出口段效应。与多颗粒相对运动模型相比，网格数量大大减少，有效提高了计算效率。

Figure 7. Grid independence verification

图7. 网格无关性验证

Figure 8. Comparison of number of grids

图8. 网格数量对比

4.3. 计算范围对比

颗粒数量是影响颗粒聚集行为和精度的关键因素。由于多颗粒主要是集中在几条相对狭窄的流线上，除了颗粒间的相互作用外，还存在多颗粒与管道的空间拥挤效应，因此Kahkeshani等[8]和Di Carlo [13]定义了无量纲多颗粒长度分数，如下式(6)所示：

$L_f=\frac{\text{6 WH}\varphi }{\pi a^2}=\frac{Na}{L}$ (6)

其中，$\varphi$ 是实验溶液中的颗粒浓度体积分数，W是通道宽度，H是通道高度，N为数值模拟中的颗粒数量，L为数值模拟中的管道长度。

该经验公式可以将实验中的颗粒浓度体积分数$\varphi$ 计算获得数值模拟中的无量纲多颗粒长度分数$L_f$ 。根据式(6)可解释4.1节中模拟结果与文献的实验结果吻合的原因：文献[8]的实验颗粒浓度体积分数$\varphi =0.004$ ，可获得颗粒链长度分数$L_f=0.111$ ，实验中监测的管道长度约为L = 40 a，因此计算可得到模拟中颗粒数量应为N ≈ 5，与实验平均每帧观测到的颗粒数量相同，因此采用多颗粒相对运动模型进行模拟，其结果与实验结果可以吻合。而文献[6]颗粒体积浓度较大$\varphi =0.1$ ，根据式6计算得到颗粒链长度分数$L_f=0.47$ 。若采用数值模拟中管长L = 50 a，再根据式6计算得到对应的颗粒数量N ≈ 24，采用多颗粒相对运动模型进行数值模拟复现网格数量过于庞大，并不现实。图9为在不同颗粒浓度体积分数($\varphi$ < 0.01)时，管道长度L与颗粒数量N对应图。颗粒数量N随着管道长度L的加长而逐渐增加，且在较高的体积分数($0.004\le \varphi \le 0.01$ )时，随管长增加，颗粒数量急剧爬升；较低体积分数($0.001\le \varphi \le 0.004$ )时，颗粒数量随管长增加较为平缓。L = 40 a、N = 5、∆x/a = 6时，网格数量已达到220万，考虑到稳态层流数值模拟中网格数量对计算效率的影响，网格数量不宜过多。而颗粒数量直接影响到网格数量，所以采用多颗粒数值模拟时，应控制颗粒数量($N\le 5$ )。

综上所述，对于颗粒数量，多颗粒相对运动模型存在一定的局限性；而颗粒链相对运动模型采用了周期性边界条件，可实现不同浓度体积分数的数值模拟[11]。

Figure 9. Correspondence diagram of multi-particle number and pipeline length

图9. 多颗粒数量和管道长度对应图

5. 结论

本文基于单颗粒相对运动模型，构建出多颗粒相对运动模型与颗粒链相对运动模型，主要对比分析了两个模型的计算结果、计算效率和计算范围，可供读者针对不同工况选用模型提供参考，结论如下：

1) 使用多颗粒相对运动模型可获得各个颗粒(尤其是首尾颗粒)在不同间距下的力学特性分布，颗粒链相对运动模型可获得周期颗粒的力学特性分布，且与多颗粒相对运动模型居中颗粒相同，两种模型的模拟结果与实验结果吻合较好，均具有较高的计算精度；

2) 不同颗粒间距下，多颗粒相对运动模型的平均网格数量是颗粒链相对运动模型的约5.3倍。随间距增大，多颗粒相对运动模型网格数量急剧增加。通过网格数量分析，多颗粒相对运动模型的计算效率远远低于颗粒链相对运动模型；

3) 根据颗粒体积长度转换式，可获得不同体积分数下，不同长度管道对应的颗粒数量。对于多颗粒相对运动模型，颗粒数量存在一定的局限性，并不适用颗粒数量较多的工况；而颗粒链相对运动模型采用了周期性边界条件，可实现高浓度体积分数的数值模拟，范围更加宽泛。

NOTES

^*通讯作者。