DLD装置中刚性颗粒运动的直接数值模拟研究

Direct numerical simulation of motion of rigid particlein a deterministic lateral displacement deviceYU Wen-guang，YU Zhao-sheng，YE Shang-jun，SHAO Xue-ming，WANG Bo-kang(Department of Mechanics，Zhejiang University，Hangzhou 310027，China)Abstract：In order to solve the problems of the accuracy of numerical simulation of deterministic lateral displacement(DLD)device andstudy of particle separation，a direct numerical simulation method(i.e.fictitious domain method)was used to simulate the motion of rigidparticle in DLD.After the analyzing the efects of lateral boundary condition，the relationship of internal flow field tilt angle and boundaryconditions were established.The periodic boundary conditions were presented to simulate DLD with a huge number of lateral arrays.The nu-merical method and boundary condition were evaluated on a characteristic of the internal flow field and the separation diameter.The criticaldiameter of a rigid particle in DLD was tested。The experimental results show that for the practical DLD devices with a huge number of lateralarrays，the periodic boundary condition with fewer arrays can be adopted to reduce the computational cost。

Key words：deterministic lateral displacement(DLD)；fictitious domain method；particle separation；critical particle diameter0 71 言 后， 冀 ，誊颗粒分离技术是生物、医学检测等领域中的-项 到广泛的关注，目前文献中已有不少的实验和应用研关键技术。确定性侧向位移分离技术(DLD)是-种 究报道 ]。关于DLD分离技术的数值研究工作还相基于粒子尺寸实现高效、连续分离颗粒的技术。 对缺乏。Loutherback等人 m 采用 COMSOL软件求解收稿日期 ：2013-01-O7基金项目：国家 自然科学基金资助项目(11072217，11202185)作者简介：喻文广(1988-)，男，安徽马鞍山人，主要从事微流控芯片内部流体力学相关问题方面的研究.E-mail：21024017###ziu.edu.cn通信联系人：余钊圣，男，教授，硕士生导师.E-mail：yuzhaosheng###zju.edu.cn第7期 喻文广，等：DLD装置中刚性颗粒运动的直接数值模拟研究 ·779·了无颗粒存在时的单相流场，然后根据速度剖面和前述近似理论计算了临界半径，并采用该方法研究了三角形微柱的形状和尺寸对分离临界半径的影响。AI-Fandi等人 采用 COMSOL软件求解了菱形和翼形微柱阵列的流场，然后研究了球形颗粒在其中的分离特性(不考虑颗粒对流场的影响)。Quek等人 采用浸没边界法数值模拟了弹性颗粒在 圆柱阵列 中的运动。

至今为止 ，关于 DLD装置中刚性颗粒的运动和临界半径的预测，都基于前述的近似理论。Inglis等人 通过实验研究了圆柱阵列中刚性颗粒的临界尺寸，发现前述近似理论的临界半径比实验值要校本研究首次采用高精度直接数值模拟方法对刚性颗粒的运动和临界半径进行研究。另外，由于实际 DLD装置中侧向圆柱数目较多(上百或者更多)，数值模拟如果采用与实际相同数量的柱列数 目则计算量太大，实际的计算只能取有限的侧向区域大小，侧向的边界为人工边界，对于人工边界如何给边界条件是非常重要的问题。本研究将首先研究该边界条件问题。

本研究首先通过计算采用侧向采用固壁边界条件时不同数量圆柱的流惩周期性边界下流场，并对流场倾斜角度进行分析，以论证侧向边界条件；然后，对刚性颗粒在 DLD中的运动情况进行模拟分析其临界分离直径，并与实验结果进行对比，验证了该方法的精确性。

1 数值模型1.1 虚拟区域方法针对流固两相流动的典型数值模拟方法 有双流体模型、拉格朗日点粒子模型和直接数值模拟。

在双流体模型中，固体颗粒被当成连续相处理，而后两个模型中颗粒是离散的，以拉格朗Et方法来追踪轨迹。

在点粒子模型中，颗粒和流体的相互作用力是按经验公式给的，对于颗粒很小并且浓度很稀的情况，这样的模型是合适的。在针对流固两相流的直接数值模拟方法中，考虑了颗粒的体积以及颗粒和流体界面上的运动学和动力学边界条件，是-种没有引入近似模型的方法，所以这种方法被称为直接数值模拟”。本研究采用的直接力/虚拟区域(Direct Forcing/Fictitious Do。

DF/FD方法是通过改进最早由 Glowinski等 提出的分布式拉格朗日乘子/虚拟区域(DLM/FD)方法得到的。该方法的中心思想是假设固体颗粒的内部也充满流体，然后通过虚拟体积力(即拉格朗日乘子)来使得这部分流体的运动符合刚体颗粒的运动。这样，流场区域就变为简单区域，可以采用固定的结构化网格来求解。

为了描述上的简要，以下公式中只考虑-个颗粒情况。假设颗粒的密度、体积、转动惯量、速度和角速度分别为P ， ，J，U和 。用 P(t)和 分别表示颗粒所占的区域和整个计算区域。采用下述特征量来无量纲化方程：特征长度 ，特征速度 ，特征时间 t/，特征压力 p， ，特征拟体力 JDr /。。无量纲化后的虚拟区域法控制方程为：Ⅷ - p n (1) a 如 。。- 、lf：U0s×rinP(t) (2)·zfOin/2 (3)( 1)V(d～U- Fr詈)-dx (4)(p -1) - ×Adx (5)式中：zf-流体的速度，p-流体的压力，A～定义在粒子区域的拟体力(拉格朗日乘子)， 颗粒中心为原点的位置矢径，p -颗粒密度 p 与流体密度 p 的比值。

p 定义为：p p ，雷诺数定义为：RePfUcL /，弗劳德数定义为：F，gL。/ ，无量纲的体积定义为： ：，无量纲的转动惯量定义为：J ∥p 。

数值离散格式的主要特点是：采用算子分裂格式将耦合方程组式(1～5)分解为两个子问题：流体子问题和拉格朗日乘子子问题。流体子问题采用基于投影格式的二阶有限差分方法求解；拉格朗日乘子问题采用直接力格式求解。

1.2 控制参数DLD装置的示意图如图 1所示。该装置由周期性布置的微柱列构成，假定圆柱间距为 ，沿着流动方向每排竖直柱列与前-排柱列侧向错开-定位移量AS，则柱列的周期为 NS/AS，倾斜柱列的斜率为 1/N。DLD分离技术的原理(近似理论)如下：当流体流过周期为Ⅳ的柱列时，竖直方向两相邻圆柱间的流体可被分成 Ⅳ条具有相同流量的窄流体束如图 1所示。流体束的相对位置(或序号)在流经每列圆柱时会发生改变，但在经过Ⅳ列圆柱之后又回其原来位置。假定与圆柱相邻的流束的宽度为R ，当颗粒的半径小于R。时，可以近似认为颗粒(中心)跟随流体在第7期 喻文广，等：DLD装置中刚性颗粒运动的直接数值模拟研究离-段距离后增长平缓，在中间区域内缓慢变化并达到最大值。因为在中间区域，壁面影响作用最小，倾斜度最大，本研究将选区中间区域的流场倾斜角度来分析边界的影响。

图4 不同侧向位置处流线的倾斜角度下面将分析不同侧向柱列数量的分离装置中间区域流场的倾斜程度，并以此来表征侧向柱列数 目以及边界条件对流场的影响。

模拟结果如图5所示，从图5可以看出，当侧向柱列数目为 14时流线倾斜角度接近于 0。，当随着侧向圆柱数量的增多，流线倾斜角度越来越大。

周期性边界条件下(N3)的流线倾斜角度两者非秤近。当圆柱数目继续增加时，倾越来越接近周期性边界条件 的计算结果。

斜角度因此 ，DLD装置中侧向圆柱数 目大于 100时，采用周期性边界条件是合理的。

图5 固壁边界下不同侧向柱列数目对倾斜角度的影响3 刚性颗粒在 DLD中的运动Inglis等人测量了刚性颗粒在 DLD分离装置中临界直径，其结果表明近似理论给出的临界直径显著低于实验值，而至今尚未有采用直接模拟方法计算刚性颗粒临界直径的报道。本研究将计算刚性颗粒的临界直径，以验证虚拟区域方法的可靠性。

· 781-在模拟中，较小的颗粒会沿着流线运动，当颗粒直径接近临界直径时，颗粒相邻两次下拐的间距会随着颗粒直径增加而增加，直至没有下拐。本研究定义颗粒在计算区域内没有出现下拐为侧向迁移模态，否则皆为之”字形模态。两种运动模态如图6所示，由图 6可知，小于临界尺寸的颗粒沿之”字形路线运动，大于临界尺寸的颗粒沿倾斜通道做侧 向迁移运动。

本研究把计算得到的临界半径和Inglis等人的实验结果做了对比，其结果如图7所示。图7中横坐标为柱列的斜率 ，即柱列周期的倒数( 1/N)；纵坐标为颗粒直径与圆柱间距的比值。从图 7可看出，数1.00.90.80 70 60.50.40 30.20 10 0O。 ， o . 磁露 oe .S。毒 面 - 饿 a1.，Lateraldisp。l ac酬em哪 O0 0.05 0lO 0.15 O.2O 0 25 0.30 0 35 O.40图7 DLD装置中刚性颗粒临界直径数值结果和实验结果的对比0 0 O 0 0 0 O 0 0 0 0 nO O o 0 0 O 0 0 0 O 00 o o o o 0 0 0 0 0 ov O 0 O O 0 O 0 0 O 0 0u 0 0 0 0 0 O 0 0 0 0 0O 0 0 O 0 0 O 0 0 0 0 n0 0 0 0 0 0 0 0o0 0 O O0 0 0 0 O 0 O 0 0 0 Ov O 0 0 O 0 0 O 0 0 0 0u 0 0 0 0 0 O 0 Q 0 0 0 0嘲 8～ ～～ ～～～～- -会 当· 782- 机 电 工 程 第 30卷值结果和实验结果吻合较好。需要说明的是：图中给出的实验拟合曲线主要根据 <0.2(N>5)的数据，对于 >0.2(N<5)，侧向迁移模态的数据很少 ，因而无法确定较精确的临界半径。对于周期 3(占1/3)阵列，模拟中没有观察到侧向迁移模态(d/D≤0.98)，而 -O.32的实验数据也没有显示存在侧向迁移模态。

4 结束语本研究用虚拟区域法对刚性颗粒在 DLD装置中的运动进行模拟，并分析了临界半径，得到了与 Inglis的实验结果非常贴近的临界直径分布，验证了该方法在计算该问题时的准确性。此外，还对侧向边界条件问题进行了研究，其结果显示：当上、下边界采用周期性边界条件时的模拟结果，与上、下边界采用固壁边界的侧向数目较多(N>100)柱列的计算结果非秤近。

所以对于模拟数量较多的实际 DLD装置，可以采用周期性边界条件简化计算量。

在验证了该方法的准确性之后，笔者将在以后的工作中模拟弹性颗粒在 DLD中的运动，研究分析颗粒的变形特性、初始形状以及初始取向等因素对临界半径的影响。