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

喻文广，余钊圣，叶尚军，邵雪明，王博康

(浙江大学力学系，浙江 杭州310027)

0 引 言

颗粒分离技术是生物、医学检测等领域中的一项关键技术。确定性侧向位移分离技术(DLD)是一种基于粒子尺寸实现高效、连续分离颗粒的技术。

DLD分离技术自2004年被 Huang等人［1］提出后，由于其具有高效、高分辨率、连续分离等优点，已受到广泛的关注，目前文献中已有不少的实验和应用研究报道［2-9］。关于DLD分离技术的数值研究工作还相对缺乏。Loutherback等人［10］采用COMSOL软件求解了无颗粒存在时的单相流场，然后根据速度剖面和前述近似理论计算了临界半径，并采用该方法研究了三角形微柱的形状和尺寸对分离临界半径的影响。AIFandi等人［11］采用COMSOL软件求解了菱形和翼形微柱阵列的流场，然后研究了球形颗粒在其中的分离特性(不考虑颗粒对流场的影响)。Quek等人［12］采用浸没边界法数值模拟了弹性颗粒在圆柱阵列中的运动。

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

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

1 数值模型

1.1 虚拟区域方法

针对流固两相流动的典型数值模拟方法［14-15］有双流体模型、拉格朗日点粒子模型和直接数值模拟。在双流体模型中，固体颗粒被当成连续相处理，而后两个模型中颗粒是离散的，以拉格朗日方法来追踪轨迹。在点粒子模型中，颗粒和流体的相互作用力是按经验公式给的，对于颗粒很小并且浓度很稀的情况，这样的模型是合适的。在针对流固两相流的直接数值模拟方法中，考虑了颗粒的体积以及颗粒和流体界面上的运动学和动力学边界条件，是一种没有引入近似模型的方法，所以这种方法被称为“直接数值模拟”。本研究采用的直接力/虚拟区域(Direct Forcing/Fictitious Domain，简写为DF/FD)方法［16］属于直接数值模拟方法。DF/FD方法是通过改进最早由Glowinski等［17］提出的分布式拉格朗日乘子/虚拟区域(DLM/FD)方法得到的。该方法的中心思想是假设固体颗粒的内部也充满流体，然后通过虚拟体积力(即拉格朗日乘子)来使得这部分流体的运动符合刚体颗粒的运动。这样，流场区域就变为简单区域，可以采用固定的结构化网格来求解。

为了描述上的简要，以下公式中只考虑一个颗粒情况。假设颗粒的密度、体积、转动惯量、速度和角速度分别为 ρs，Vp，J，U 和 ωp。用 P(t)和 Ω 分别表示颗粒所占的区域和整个计算区域。采用下述特征量来无量纲化方程:特征长度Lc，特征速度Uc，特征时间Lc/Uc，特征压力，特征拟体力。无量纲化后的虚拟区域法控制方程为:

式中:u—流体的速度，p—流体的压力，λ—定义在粒子区域的拟体力(拉格朗日乘子)，r—颗粒中心为原点的位置矢径，ρr—颗粒密度 ρs与流体密度 ρf的比值。

ρr定义为:ρr= ρs/ρf，雷诺数定义为:Re= ρfUcLc/μ，弗劳德数定义为，无量纲的体积定义为:，无量纲的转动惯量定义为:J*=J/。

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

1.2 控制参数

DLD装置的示意图如图1所示。该装置由周期性布置的微柱列构成，假定圆柱间距为S，沿着流动方向每排竖直柱列与前一排柱列侧向错开一定位移量ΔS，则柱列的周期为N=S/ΔS，倾斜柱列的斜率为ε=1/N。DLD分离技术的原理(近似理论)如下:当流体流过周期为N的柱列时，竖直方向两相邻圆柱间的流体可被分成N条具有相同流量的窄流体束如图1所示。流体束的相对位置(或序号)在流经每列圆柱时会发生改变，但在经过N列圆柱之后又回其原来位置。假定与圆柱相邻的流束的宽度为Rc，当颗粒的半径小于Rc时，可以近似认为颗粒(中心)跟随流体在流体束内部运动，呈现“之”字形(zigzag)运动轨迹;当颗粒半径大于Rc时，如颗粒中心位于第2条流束，则可以近似认为颗粒在流经每列圆柱时都会处于两圆柱间的第2条流束位置，因此颗粒会沿着倾斜通道一直运动，导致侧向位移。根据不同尺寸的粒子运动轨迹不同的特性，就可以实现分离。在实际应用过程中，可以将不同临界半径的柱列组合起来，实现更大尺寸范围的粒子的分离。

本研究只考虑固体和流体密度相等的情况，即ρr=1。也不考虑重力的作用，即Fr=0。不考虑颗粒的布朗运动，在控制方程式(1～5)的推导中已经做了忽略。

如图1所示，(初始)颗粒直径为d(半径为a)，固定圆柱的直径为D，两个方向的圆柱中心距均为S，间隙为Lg，柱列的周期为N，斜率为ε=1/N。本研究只考虑D=Lg=S/2的情形。

图1 DLD装置示意图(周期N=3)

取均匀来流速度U0为特征速度，圆柱直径D为特征长度。因而:

对于实际DLD装置中的流动，雷诺数比较低，流体惯性可以忽略。本研究取Re=0.1，验证算例表明在Re＜1时，雷诺数取值不影响颗粒的运动。

2 边界条件和数值参数

关于边界条件，流向入口为均匀来流，出口为无反射出流条件，侧向取周期性边界条件，计算域以及边界条件如图2所示。计算区域大小为56D×28D，圆柱数量为24×14。流场网格数为2 048×1 024，一个圆柱直径覆盖约37个网格。时间步长为0.00 025D/U0。对于d/D=0.6的较大颗粒，验证算例表明:把上面给定的网格尺寸和时间步长加大一倍，对结果的影响很小。限于计算条件，本研究所研究的颗粒尺寸d＞0.2D，圆柱阵列周期N≤7。

图2 计算域以及边界条件

对于侧向边界条件，实际的DLD装置在上、下两侧为固壁，且侧向的圆柱数量较多(如数百个)，如果采用实际的边界条件和圆柱数量，计算量很大。周期性边界条件相当于侧向的圆柱数量无穷多，如果采用周期性边界计算量将大大减小。但是周期性边界与实际装置的固壁边界条件又会有所差异，下面就该问题进行探讨。

本研究分别模拟了上下边界采用固壁边界条件，圆柱数目分别为:24×14、24×28、24×57、24×115的柱列无颗粒流动，和上、下边界采用周期性边界条件时柱列数量为(24×28)柱列无颗粒流动。

通过分析发现，DLD中流场并非完全像上文提到的近似理论那样:假定DLD装置中的流线是周期性的，其周期等于圆柱阵列周期，因而平均意义上的流线是水平的，没有侧向偏移。事实上，由于倾斜通道的存在，流线有总体的倾斜，且倾斜角度随上、下壁面间距增大而增加(倾斜角度α定义为:与流线相切的直线与水平方向的夹角，倾斜流线图如图3所示)。靠近上壁面处的流体倾斜角度较小，离壁面越远倾斜角度增大。不同侧向位置的流线倾斜角度与位置关系如图4所示，其中横坐标为无量纲化侧向位置，纵坐标为不同位置流线的倾斜角度。从图4可以看出，在靠近壁面处由于壁面的存在使得流线倾斜较小，但随着离远离壁面增大，流线的倾斜角度迅速增大，当超过壁面距离一段距离后增长平缓，在中间区域内缓慢变化并达到最大值。因为在中间区域，壁面影响作用最小，倾斜度最大，本研究将选区中间区域的流场倾斜角度来分析边界的影响。

图3 周期N=3、上下边界为固壁边界条件、柱列数目为24×115的无颗粒流场中的倾斜流线

图4 不同侧向位置处流线的倾斜角度

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

模拟结果如图5所示，从图5可以看出，当侧向柱列数目为14时流线倾斜角度接近于0°，当随着侧向圆柱数量的增多，流线倾斜角度越来越大。侧向圆柱柱列数目为115时，其偏移角度α≈3.4°。而当采用周期性边界条件下(N=3)的流线倾斜角度α≈3.7°，两者非常接近。当圆柱数目继续增加时，倾斜角度会越来越接近周期性边界条件的计算结果。因此，当DLD装置中侧向圆柱数目大于100时，采用周期性边界条件是合理的。

图5 固壁边界下不同侧向柱列数目对倾斜角度的影响

3 刚性颗粒在DLD中的运动

Inglis等人测量了刚性颗粒在DLD分离装置中临界直径，其结果表明近似理论给出的临界直径显著低于实验值，而至今尚未有采用直接模拟方法计算刚性颗粒临界直径的报道。本研究将计算刚性颗粒的临界直径，以验证虚拟区域方法的可靠性。

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

图6 不同尺寸颗粒的运动模态

图7 DLD装置中刚性颗粒临界直径数值结果和实验结果的对比

本研究把计算得到的临界半径和Inglis等人的实验结果做了对比，其结果如图7所示。图7中横坐标为柱列的斜率ε，即柱列周期的倒数(ε=1/N);纵坐标为颗粒直径与圆柱间距的比值。从图7可看出，数值结果和实验结果吻合较好。需要说明的是:图中给出的实验拟合曲线主要根据ε＜0.2(N＞5)的数据，对于ε＞0.2(N＜5)，侧向迁移模态的数据很少，因而无法确定较精确的临界半径。对于周期3(ε=1/3)阵列，模拟中没有观察到侧向迁移模态(d/D≤0.98)，而ε≈0.32的实验数据也没有显示存在侧向迁移模态。

4 结束语

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

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

［1］HUANG L R，COX E C，AUSTIN R H，et al.Continuous particle separation through deterministic lateral displacement［J］.Science，2004，304(5673):987-990.

［2］LOUTHERBACK K，PUCHALLA J，AUSTIN R H，et al.Deterministic microfluidic ratchet［J］.Physical Review Letter，2009，102(4):45301-45305.

［3］BEECH J P，HOLM S H，ADOLFSSON K，et al.Sorting cells by size，shape and deformability［J］.Lab on a Chip，2012，12(6):1048-1051.

［4］DAVIS J A，INGLIS D W，MORTON K J，et al.Deterministic hydrodynamics:taking blood apart［J］.Proceedings of the National Academy of Sciences of the United States of America，2006，103(40):14779-14784.

［5］NAGRATH S，SEQUIST L V，MAHESWARAN S，et al，I-solation of rare circulating tumour cells in cancer patients by microchip technology［J］.Nature，2007，450(7173):1235-1239.

［6］INGLIS D W，MORTON K J，DAVIS J A，et al.Microfluidic device for label-free measurement of platelet activation［J］.Lab on a Chip，2008，8(6):925-931.

［7］HOLM S H，BEECH J P，BARRETT M P，et al.Separation of parasites from human blood using deterministic lateral displacement［J］.Lab on a Chip，2011，11(7):1326-1332.

［8］FRECHETTE J，DRAZER G.Directional locking and deterministic separation in periodic arrays［J］.Journal of Fluid Mechanics，2009(627):379-401.

［9］BALVIN M，SOHN E，IRACK T，et al.Directional locking and the role of irreversible interactions in deterministic hydrodynamics separations in microfluidic devices［J］.Physical Review Letters，2009，103(7):78301-78304.

［10］LOUTHERBACK K，CHOU K，NEWMAN J，et al.Improved performance of deterministic lateral displacement arrays with triangular posts［J］.Microfluidics and Nanofluidics，2010，9(6):1143-1149.

［11］Al-FANDI M，Al-ROUSAN M，JARADAT M，et al.New design for the separation of microorganisms using microfluidic deterministic lateral displacement［J］.Robotics and Computer-Integrated Manufacturing，2011，27(2):237-244.

［12］QUEK R，LE D V，CHIAM K H.Separation of deformable particles in deterministic lateral displacement devices［J］.Physical Review E，2011，83(5):56301-56307.

［13］INGLIS D W，DAVIS J A，AUSTIN R H，et al.Critical particle size for fractionation by deterministic lateral displacement［J］.Lab on a Chip，2006，6(5):655-658.

［14］解婧陶.板料成形数值模拟中回弹问题的处理［J］.机电技术，2011(6):91-93.

［15］杨 阳，赵建平.EPS聚合反应釜内混合过程的数值模拟［J］.轻工机械，2012，30(6):1-4.

［16］YU Z，SHAO X.A direct-forcing fictitious domain method for particulate flows［J］.Journal of Computational Physics，2007，227(1):292-314.

［17］GLOWINSKI R，PAN T W，HESLA T I，et al.A distributed Lagrange multiplier/fictitious domain method for particulate flows［J］.International Journal of Multiphase Flow，1999，25(5):755-794.