上海理工大学，上海

收稿日期：2020年12月1日；录用日期：2021年1月30日；发布日期：2021年2月7日

摘要

基于“相对运动模型”，利用CFD技术数值计算在矩形截面通道中刚性球形颗粒运动的力学特性，获得球形颗粒所受惯性升力的空间分布特征。分析使其稳定聚集位置产生变化的因素，为矩形截面通道中颗粒的惯性聚集现象的研究及可控性商业应用提供指导。研究结果表明：在矩形截面通道中，颗粒受到惯性升力水平分量的作用，远离短边壁面的中心附近向长边壁面方向移动，使本应存在的两个稳定聚集位置消失；流体雷诺数的不同会使颗粒所受的惯性升力竖直分量产生差异性变化，雷诺数较小时会使矩形截面通道中四个直角处的惯性聚集位置消失；截面通道高宽比会影响惯性升力水平分量的分布特征，颗粒的自转效应是这一现象的主要影响因素。

关键词

惯性聚集，刚性颗粒，惯性升力，数值研究

The Mechanical Causes of Inertial Focus of Particles in Rectangular Section Channel

He Zou, Qikun Wang, Zhuangzhuang Xue

University of Shanghai for Science & Technology, Shanghai

Received: Dec. 1^st, 2020; accepted: Jan. 30^th, 2021; published: Feb. 7^th, 2021

ABSTRACT

Based on the “Galileo Principle of Relativity”, a numerical model is proposed to describe the steady and constant motion of spherical particles, and the spatial distribution of the inertial lift of spherical particles in rectangular cross-section channel is studied by combining with CFD technology. The factors that cause the change of the stable aggregation position are analyzed, which can provide guidance for the study of the inertial aggregation of particles in rectangular cross-section channel and the controlled commercial application. The result of research shows that in the rectangular cross-section channel, the particles are affected by the horizontal component of the inertial lift, moving away from the center of the short side wall and towards the long side wall, so as to make the two stable aggregation positions that exist should disappear. Different Reynolds number of fluid will cause different changes in the vertical component of the inertial lift of particles. The ratio of height to width of the cross-section channel affects the distribution characteristics of horizontal components of inertial lift, and the rotation effect of particles is the main influencing factor of this phenomenon.

Keywords:Inertial Focus of Particles, Rigid Particle, Inertial Lift, Numerical Investigation

Copyright © 2021 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

Segré & Silberberg于1961年通过实验 [1] [2] 研究了悬浮(颗粒与流体的密度相等)颗粒在圆管直通道内随流体运动的情况。研究结果表明：颗粒在聚集过程中表现出自发性，即无需额外施加电或磁场，以现象为基础制作的微流控芯片已成功实现了颗粒和细胞的分离。

已有研究表明：颗粒惯性聚集的形态完全决定于通道截面的几何形状 [3]：圆形截面通道内，颗粒会聚集成为一个同心圆环 [4] (如图1(a)所示)；方形截面通道内，颗粒聚集在通道壁面中心附近，呈对称分布的四个点 [5] (如图1(b)所示)；矩形截面通道内，颗粒的聚集区域仅分布于长壁面的中心附近，呈对称分布的两个点(如图1(c)所示)，而短边壁面中心附近并不存在聚集区域 [6]。若通道Re变大，圆形和方形截面通道中颗粒的聚集位置会向壁面方向靠近 [7]；而在矩形截面通道中，颗粒会聚集在平行于长边壁面的附近，呈现对称的两列聚集点分布(如图1(c)所示)。

图1. 颗粒惯性聚集位置变化示意图

事实上，已有研究表明，颗粒的聚集形态完全决定于其所受流体的惯性力空间分布特征 [8]。而这种空间分布特征会随着通道雷诺数、颗粒尺寸以及通道截面形状的改变而变化，这是导致在不同条件下颗粒聚集形态产生差异性的根本原因。这就能很好地解释图1中三种不同截面形状的通道内，颗粒聚集形式不同的力学成因。

因此，探索与研究通道层流流场中颗粒所受惯性升力的空间分布特征及其受通道截面几何形状的影响规律是解释颗粒惯性聚集现象一般规律的根本且有效的途径。然而，不能充分理解颗粒在微通道中惯性聚集机理的原因在于研究颗粒受力困难。甚至是在宏观尺度，颗粒受力也很难用实验测出，大多是数值结果。

Asmolov [9] 利用摄动法从理论上获得颗粒在微通道内的力学特性，然而使摄动解成立的前提是基于“点颗粒”假设，颗粒的存在对其周围流场无影响。但实际情况中并不能忽略颗粒的大小。Carlo et al. [10] 基于“相对运动模型”，数值研究了方形截面管道流动中颗粒的横向升力，通过数据拟合的方法获得了该升力各个部分的近似表达式；王企鲲等 [11] 数值研究了颗粒在微通道中的力学特性，但是在计算中忽略颗粒自身旋转。

因此本文采用“相对运动模型”并结合CFD技术，数值研究了在矩形截面通道内，颗粒所受惯性升力的空间分布特征及其受通道截面高宽比和通道雷诺数的影响规律，以此旨在揭示颗粒惯性聚集形态变化特征的力学成因。研究成果有助于完善低雷诺数固–液两相流动中的力学机理，为进一步深入研究惯性升力的力学特征及惯性微流控芯片的商业应用提供有效帮助和理论指导。

2. 计算模型与方法

2.1. 数值计算模型简介

该流动本质上为一个含稀相颗粒的液固两相流动，因此为了研究简便，本文研究对象取为悬浮于液相中的单个球状刚性颗粒在液相的层流通道内稳定运动(即颗粒浓度很低，颗粒与颗粒之间的相互作用几乎可被忽略 [12] )。但由于现有的两相流动数值方法中多包含了很多经验性的参数需人为设置，具有很大的近似性与不确定性，本文并不打算采用常规的两相流动的数值计算方法处理这个问题，而是直接将颗粒的表面视为流场的一个运动着的固体边界。如此，所研究的流场为内部挖去一个球的通道内流动，如图2(a)所示。这便转化为常规的单相流的问题。

颗粒周围流体对其作用力的合力沿垂直于主流方向的分量为颗粒的惯性升力。过CFD的方法计算出来，便不难获得颗粒表面的应力张量P，进而按式(1)即可获得颗粒所受的惯性升力F_L。

由于颗粒是运动的，在如图2(b)中所示的置于通道中心的O-XYZ坐标系，为典型的非定常、动边界流动，在商用软件中，往往会选动网格和时间推进等技术进行CFD求解，计算过程十分复杂且计算周期极长。因此本文利用运动相对性的基本思想建立一个相对运动计算模型，以此将原本复杂的非定常计算转化为定常的计算，使计算过程大为简化。具体实施过程如下文所述。

如图2所示，以球形颗粒的中心为坐标原点，并使之随颗粒等速平移。在此平移坐标系中：颗粒相对静止，壁面转化为速度为U的匀速直线运动，从而此非定常流场转变为相对定常的绕流场。以上处理方法为本文的数值研究提供了便利。对于颗粒的速度本文采用“试凑法”加以确定：假设一个颗粒速度U_P，进行试算后，获得颗粒在主流方向(即X轴方向)所受的合力F_x，不断调节U_P使F_x在计算精度范围内趋近于零(即F_x小于F_y两个量级以上)，此U_P即为所需的颗粒运动速度。利用相同的方法可以试凑出小球颗粒的旋转速度。

图2. 相对运动模型

$F_L={\displaystyle {∯}_{\Sigma }n}\cdot P\text{d}S$ (1)

$F_y=j\cdot F_L$ (2)

$F_z=k\cdot F_L$ (3)

式(1)中：P为应力张量，n为颗粒表面外法线单位向量，S为颗粒表面积。

式(2)、(3)中，j、k分别为O-XYZ坐标系中Y轴与Z轴的单位向量。

一般而言，周围流体对颗粒表面的作用力存在沿Y轴和Z轴两个分量，本文记为F_y和F_z，它们可以式(2)、式(3)被确定。

2.2. 数值计算方法

经网格独立性检查后，模型网格数量确定为70万。计算模型采取三维、相对定常、不可压缩N-S方程组，如式(4)所示：

$\{\begin{array}{l}\nabla \cdot W=0\\ \left(W\cdot \nabla \right)W=-\frac{1}{\rho }\nabla p+\upsilon {\nabla }^{2}W\end{array}$ (4)

式(4)中，W为平移坐标系中的速度(即相对速度)，p为流体压强；ρ、υ为流体密度及运动粘度，均为常数。

若记静止坐标系O-XYZ中，流体速度为V、颗粒平移速度为U_p，则根据运动相对性原理，有：

$V=W+U_p$ (5)

当数值求解式(4)后，即可获得相对速度场W及压力场p，于是便可根据式(5)计算流场的绝对速度，进而根据式(6)、式(7)，便可获得流场的变形速率张量S及应力张量P。

$S=\frac{1}{2}\left[\left(\nabla W\right)+{\left(\nabla W\right)}^{\text{T}}\right]$ (6)

$P=-pI+2\mu S$ (7)

上式中，I为二阶单位张量、μ为流体的动力粘度(为常数)。

在求解过程中，压力、速度的耦合采用精度更高的“Coupled”算法进行计算；边界条件：将边界条件设置为速度入口(Velocity-inlet)，速度的具体数值取为流体在不同雷诺数下的运动速度；将右侧端面的边界条件设为压力出口(Pressure-outlet)，并将参考压力设定在计算域出口的中心位置处。

2.3. 相关物理参数

为便于下文表述，本文定义无量纲参数如下：

升力系数

$C_{Fy}=\frac{F_y}{\rho U^2{a}^4/D^2}$ (8)

$C_{Fz}=\frac{F_z}{\rho U^2{a}^4/D^2}$ (9)

式中，U为通道中流体的平均流速；a为球形颗粒的直径；D为通道截面的直径或当量直径；

无量纲颗粒直径

$a^{+}=\frac{a}{D}$ (10)

无量纲横向位置

$y^{+}=\frac{Y}{H/2}$ (11)

$z^{+}=\frac{Z}{W/2}$ (12)

式中，Y、Z为O-XYZ坐标系中坐标轴Y、Z的坐标值，具体坐标系见图2(b)。

通道Re数

$Re=\frac{\rho UD}{\mu }$ (13)

矩形截面通道的高宽比RA

$RA=\frac{H}{W}$ (14)

式中，H、W分别为通道横截面的高度与宽度，如图2(b)所示。

3. 结果与分析

3.1. 短边中线上惯性聚集点消失的力学分析

在不同高宽比RA的矩形截面通道中，图3给出了粒径a⁺ = 0.22的颗粒在短边中线上所受惯性升力竖直分量F_y的分布曲线。

由图3可以看出，在不同高宽比的通道中，颗粒的惯性升力分布特征存在明显差异，当RA = 1时，方形截面中线上颗粒所受的惯性升力竖直分量F_y随y⁺的改变呈现正负波动，具有单一零点。但是只有当升力曲线的空间分布上表现为升力曲线在该点为零点且在该点沿径向的斜率(一阶导数)为负值这才是颗粒稳定的聚集点。当RA > 1时，矩形截面短边中线上颗粒所受的惯性升力竖直分量F_y随y⁺位置的变化呈现不规则的分布情况，且具有多个零点位置。根据这些零点位置相邻区域内F_y的方向判断，这些零点位置中存在F_y的稳定零点位置，同时，靠近短边壁面的零点位置必然为稳定零点位置。然而现有的实验表明矩形截面通道中，在短边中线上是不存在颗粒的稳定聚集点，因此极有可能同时存在沿Z轴方向的惯性升力水平分量F_z作用于颗粒上，迫使其离开短边中线位置。

图3. 短边中线上颗粒所受F_y的分布情况

为应证上述猜想，图4给出了惯性升力水平分量F_z的空间分布曲线。

图4. 短边中线上颗粒所受F_z的对比情况

图4说明，方形截面通道(RA = 1)中，通道截面中心线上惯性升力水平分量几乎为零，远远小于对应的惯性升力竖直分量，使颗粒在Y轴方向进行聚集运动时，惯性升力水平分量几乎不起作用，因而颗粒被稳定地聚集在方形通道截面四边的中点附近；在矩形截面通道(RA > 1)中，截面短边中线上颗粒所受的惯性升力存在明显的水平分量。处于短边中线上的颗粒在惯性升力竖直分量F_y的作用下进行沿Y轴方向的惯性聚集运动时，会在水平分量F_z的作用下偏离短边中线位置，向长边壁面方向移动。

综上所述，矩形截面通道的短边中线上的惯性升力除了存在竖直分量F_y外，还存在水平分量F_z，这是与方形截面通道上惯性升力空间分布特征最显著的差异。正是在这个F_z的作用下，颗粒会偏离短边中线位置，向长壁方向移动。这是导致矩形截面短边中线上颗粒不存在稳定聚集点的力学原因。

3.2. 颗粒惯性聚集形态的力学分析

为进一步研究在不同高宽比RA的矩形截面通道中，颗粒惯性聚集点不同的分布情况的力学成因。本文选取在矩形截面通道短边中线上F_y稳定零点位置处，数值计算颗粒在此y⁺位置的右侧不同z⁺位置上所受的惯性升力水平分量F_z，如图5所示。

图5. 颗粒在不同y⁺位置的侧向位置所受C_FZ

图5表明，矩形截面通道中(RA > 1)，当颗粒偏离短边中线位置时，受到的惯性升力水平分量先指向长边壁面方向，使其继续远离短边中线向长边壁面移动；F_z的数值大小随短边中线距离的增加呈现出先上升后下降的趋势，到零值后改变受力方向，即指向短边中线方向；F_z在Z轴方向存在单一零点，且根据相邻区域的受力方向判断，此零点位置为F_z的稳定零点位置；随着高宽比的增大，稳定零点位置会向长壁方向移动。

有上述分析易知：在矩形截面通道内壁，颗粒惯性聚集运动是在F_y和F_z的共同作用下完成的；矩形截面上颗粒聚集点的形成，必须同时存在F_y和F_z的稳定零点位置。

为研究矩形通道截面直角区域中颗粒所受惯性升力竖直分量F_y的分布情况，本文选取高宽比RA = 4的矩形截面通道为例，在图6中给出惯性升力竖直分量沿截面Z轴正向聚集点上侧的不同y⁺位置的分布曲线。

图6中可以看出，在y⁺较小的区域，颗粒所受的惯性升力竖直分量F_y的方向始终指向下方，即通道截面长边中线方向；在y⁺较大的区域，雷诺数Re较大的工况下，惯性升力竖直分量F_y存在稳定零点位置；随着雷诺数Re的减小，惯性升力竖直分量F_y的稳定零点位置也在逐渐向长边壁面的中线方向靠近。在雷诺数Re较小(Re = 50)的工况下，惯性升力竖直分量不存在零点位置，其方向始终指向长边壁面的中线方向。

图6. 颗粒在长边壁面附近不同位置所受C_Fy

结合之前提到的惯性升力水平分量F_y的分布规律，可以推断：当雷诺数Re较大时，矩形截面通道中的颗粒在惯性升力竖直分量F_y和惯性升力水平分量F_z的共同作用下，于通道截面四个直角区域附近存在稳定的惯性聚集区域。当雷诺数Re较小时，尽管惯性升力水平分量F_z在z⁺较大区域仍然存在稳定零点位置，但是此区域附近的颗粒在惯性升力竖直分量F_y的作用下，会向截面长边壁面的中线方向移动。最终，此工况下的矩形截面通道中，颗粒只会在长边壁面中线上的特定位置，形成两个稳定的惯性聚集点，而在通道截面四个直角区域附近，已不存在稳定的颗粒惯性聚集区域。这个结论与Zhou等人 [6] 的实验结果相一致。

3.3. 惯性升力水平分量差异性分布成因分析

在前文分析中可以看出：不同高宽比的矩形截面通道中，颗粒所受的惯性升力水平分量F_z是导致其形成不同聚集形态的决定性因素；图5中还可得出结论：当RA = 1时，即在方形截面通道中，F_z始终指向通道壁面中线方向；当通道高宽比RA > 1时，截面短边中线附近(z⁺较小时)，颗粒所受F_z的方向会发生突变，由指向短边中线方向转为指向长边壁面方向；在长边壁面附近(z⁺较大时)，F_z的分布情况则较为相似。

为研究在z⁺较小的情况下，通道中颗粒所受惯性升力水平分量F_z发生突变的原因，本文将颗粒所受惯性升力水平分量F_z分解为两部分：一部分为颗粒自转产生的旋转诱导升力F_ZW，另一部分为剪切诱导升力与壁面诱导升力共同组成的非旋转诱导惯性升力F_ZS。在相对运动模型中，非旋转诱导惯性升力F_ZS可以通过忽略颗粒的自转效应获得；旋转诱导升力F_ZW可以通过前文得到的颗粒所受惯性升力水平分量F_z与非旋转诱导惯性升力F_ZS的差值获得。

将图5中不同高宽比RA的矩形截面通道模型中，按照上述方法得到颗粒所受的惯性升力水平分量F_z不同部分的分布特征后，分别进行对比，其结果如图7所示。

研究表明：在不同高宽比RA的矩形截面通道中，非旋转诱导惯性升力F_ZS的方向始终指向通道短边壁面中线方向，随着z⁺的增大，其数值大小也大致相同；旋转诱导升力F_ZW的方向始终指向通道长边壁面方向，随着z⁺的增加，其数值大小先逐渐增大，至最大值后开始减小；在高宽比RA > 1的矩形截面通道中，短边壁面的中线附近(即z⁺较小的区域)，旋转诱导升力F_ZW的差异性分布与惯性升力水平分量F_z的分布情况具有一致性。因而可以推测：随着矩形截面通道高宽比RA的增大，由颗粒自转效应产生的旋转诱导升力F_ZW的分布特征发生变化，进而影响了颗粒所受F_z的分布情况。

图7. 颗粒所受F_z拆分后不同部分的分布情况

4. 结论

1) 矩形截面通道中，在不可忽视的惯性升力水平分量Fz的作用下，颗粒仅于长壁附近存在稳定聚集点，在短壁中心附近已不存在稳定聚集位置。

2) 矩形截面通道中，流体雷诺数会影响惯性升力竖直分量F_y的分布特征。当Re较小时，颗粒只会在长边壁面中线上的特定位置，形成两个稳定的惯性聚集点；当Re增大时，原本仅存在与长边壁面中心处的两个惯性聚集点增加为六个，通道截面的四个直角区域也出现了颗粒的惯性聚集点。

3) 对比方形和矩形截面通道，高宽比RA会影响惯性升力水平分量F_z的分布特征。当通道高宽比RA > 1时，惯性升力水平分量F_z的分布特征会发生突变，之后随着高宽比RA的增大，颗粒的惯性聚集位置向长边壁面方向移动；高宽比RA对于惯性升力水平分量F_z的影响，主要是通过颗粒的自转效应实现的。

基金项目

国家自然科学基金(51776128)和国家教育部博士点基金(20113120120003，20130073110059)。

文章引用

邹 赫,王企鲲,薛壮壮. 矩形截面通道中颗粒惯性聚集现象的力学成因
The Mechanical Causes of Inertial Focus of Particles in Rectangular Section Channel[J]. 建模与仿真, 2021, 10(01): 83-92. https://doi.org/10.12677/MOS.2021.101009

参考文献