胡晟^1,2, 吴怀京¹, 吕晓永^1,2
1. 东北大学秦皇岛分校 控制工程学院, 河北 秦皇岛 066004;
2. 东北大学秦皇岛分校 河北省微纳精密光学传感与检测技术重点实验室, 河北 秦皇岛 066004
收稿日期：2020-03-26
基金项目：国家自然科学基金资助项目(61903069)；河北省自然科学基金资助项目(F2020501040)。
作者简介：胡晟(1984-)，男，云南景洪人，东北大学副教授。

摘要：采用有限元法进行粒子电极化强度的计算, 以此得到对应的介电泳力和次级介电泳力.粒子的极化强度依托软件COMSOL Multiphysics 5.3a的AC/DC模块计算得到, 粒子追踪模块用于实现介电泳力、次级介电泳力、流体力, 以及排斥力的粒子动力学仿真.AC/DC模块和粒子追踪模块的多物理场耦合新方法能够满足正、负介电泳效应的粒子运动分析和研究, 且粒子轨迹与常用的粒子仿真算法的运算结果有较高的相似性.同时仿真结果还表明, 为了更准确地预测实验观测结果, 粒子承受金属电极形成的正、负介电泳效应需要考虑次级介电泳力.
关键词：介电泳粒子链流体动力学COMSOL
A Novel Approach to Simulate Particle-Particle Interaction and Dielectrophoretic Dynamics
HU Sheng^1,2, WU Huai-jing¹, LYU Xiao-yong^1,2
1. School of Control Engineering, Northeastern University at Qinhuangdao, Qinhuangdao 066004, China;
2. Hebei Key Laboratory of Micro-Nano Precision Optical Sensing and Measurement Technology, Northeastern University at Qinhuangdao, Qinhuangdao 066004, China
Corresponding author: HU Sheng, E-mail: husheng@neuq.edu.cn.

Abstract: Electrical polarization strength of a particle is calculated by finite element method to obtain the dielectrophoretic and the secondary dielectrophoretic forces. Relying on the commercial software COMSOL Multiphysics 5.3a, the electrical polarization strength can be obtained in the AC/DC module from COMSOL, and then the module of particle tracing for fluid flow(PTF)is used to realize the dynamic simulation of particles on which the dielectrophoretic force, secondary dielectrophoretic force, fluid force and repulsive force are exerted. This novel method in which AC/DC module coupled PTF module in multi-physical fields provides an insight into the particle trajectories with positive or negative dielectrophoresis, which has a good agreement in other particle simulation methods. Whatever positive or negative dielectrophoresis is, these results imply that it is essential to consider secondary dielectrophoretic effect on metal electrodes to improve the prediction of experimental observation.
Key words: dielectrophoresisparticle chainhydrodynamicsCOMSOL
介电泳作为一种微纳粒子的操控技术已经广泛应用于微流控芯片的细胞分选、蛋白质富集、染色体输运和基因检测等领域^[1-2].近年, 介电泳力成功汇集大量胶体颗粒或生物个体, 使它们形成的微组织器官具有了特定生化和生理功能^[3-5].在微粒受介电泳力的驱动过程中, 由显微成像可以经常发现粒子链状结构的产生和演化^[6].当粒子远离较强场强区域的金属电极, 粒子会扭曲电场形成非均匀电场, 进而诱导周围粒子彼此吸引形成珍珠链结构.因此, 次级电泳力是电场作用下粒子成链的根本原因.
国内外已有大量****针对粒子成链和运动特征进行了理论建模和分析.Ai等^[7-8]提出了粒子-流场-电场的多物理场模型, 采用任意拉格朗日-欧拉算法对电场(介电泳力)和流场(微流动力)进行求解, 结果显示粒子相互作用力与电场强度和粒子大小成正比, 当粒子间的方向矢量与电矢量平行, 粒子的相互吸引力最大.Hossan等^[9]使用混合浸入边界法求解电解液中的电场和流场分布: 模型无需对粒子边界或球体进行闭合积分(麦克斯韦应变张量法), 也能实现相同或不同粒子的时域模拟.Kang^[10]采用有限体积法实现交流电场中麦克斯韦应力张量与平均时间介电泳力的计算.锐界面法涉及交流电场分布和强度的求解, 浸入边界法则可获得流场运动信息.仿真结果指出Clausius-Mossotti(CM)因子实部数值是粒子受介电泳力运动的重要影响因素, 不同粒子的CM因子实部符号相同, 则成链方向与电矢量方向相同; 反之, 它们的链型结构与电矢量方向垂直.Xie等^[11]采用等效电偶极矩法模拟粒子受介电泳力的运动轨迹; 为了更好地解决电偶极矩法对介电泳力计算精度低的问题, 他们提出一种迭代电偶极矩法, 该方法可以有效降低网格计算量和仿真时间, 粒子受力和运动轨迹都与麦克斯韦应力张量法计算结果保持了较高的近似程度.
上述建模与仿真研究为次级介电泳效应的粒子成链提供了坚实的理论基础; 然而, 上述有限元方法涉及网格形变量, 以及时间步长没有较为准确的说明, 且机器运算量较大, 计算时间较长.同时, 大多数仿真都是借助有限元或有限边界法进行电场的计算, 再根据电场进行介电泳力分析和粒子运动模拟.当前介电泳效应的电场计算主要依托COMSOL Multiphysics多物理场有限元软件进行偏微分方程求解^[12-13].有限元COMSOL软件内联的AC/DC模块可进行Laplace方程求解.另外COMSOL软件还包含了粒子追踪(particle tracing for fluid flow，PTF)模块, 可用于流体动力学中粒子的应力分析和运动模拟^[14], 便于直观和高效地学习和理解物理场的分布、强度和关联度.因此, 本文在总结前期介电泳力物理关系式的基础上, 结合AC/DC模块和PTF模块对介电泳的粒子运动进行仿真和研究.本文提出的方法不仅能够降低编程门槛, 更能高效地实现电极结构的前期研发设计, 为介电泳的粒子分选、排列、输运实验提供必要的理论指导.
1 模型和原理根据电偶极矩法描述, 复介电常数为=ε_p-jσ_p的球形粒子悬浮于复介电常数=ε_m-jσ_m的去离子水溶液中.其中ε_p, ε_m分别为粒子和溶液的介电常数; σ_p, σ_m为粒子和溶液的电导率.因此, 时间平均的交流介电泳力公式如下:

(1)

式中: a代表粒子半径; 是交流电场的均方根；K_CM即是CM因子的实部:

(2)

当K_CM>0时, 正介电泳力驱使粒子运动至电场强度最强位置; 反之K_CM＜0时, 负介电泳力排斥粒子至电场强度最弱区域.本文模拟粒子a为5 μm, 二维方形电场尺寸为L×L(100 μm×100 μm), 如图 1所示.激励电压为第一类狄利克雷边界条件, 满足; 上、下边界(n为边界的单位法矢)代表绝缘边界, 满足第二类纽曼边界条件.仿真粒子的电学参数ε_p为2.5ε0, σ_p分别为5.0×10^-3 S/m, 8.0×10^-4 S/m; 溶液电学参数ε_m为78ε₀, σ_m为2.0×10^-3 S/m.此外, 为了直观地研究介电泳效应, 本文假设电场变化的特征时间远远小于力学(包括流体力)弛豫时间, 忽略电解质溶液自由电荷对流所造成的电渗流效应^[15].由式(2)可得两种不同电导率相同尺寸微球的CM频谱实部曲线, 如图 2所示.频率范围为10³~10⁶ Hz, 电导率高的粒子受正介电泳力, 反之粒子受负介电泳力.
图 1(Fig. 1)

图 1 仿真粒子空间位置与电场边界条件Fig.1 Spatial position of simulated particles and boundary condition of electric field

图 2(Fig. 2)

图 2 两种电导率粒子的实部CM曲线Fig.2 Curves of the real of CM factor for two particles with different conductivities

然而在均匀电场作用下, 粒子成链主要由于粒子扭曲电场形成的偶极子静电互动力所致.次级介电泳力(即偶极子静电互动力)的表达式为^[12]

(3)

式中: p_i(或p_j)等于; R_ij是粒子i与粒子j的欧拉距离矢量, R_ij是R_ij的模.当粒子半径远大于特征长度(如电极长度)时, 该表达式可以忽略高阶极矩效应.同时, 介质内运动的粒子i还受流体动力的作用, 即

(4)

式中: η为黏度; I为二阶单位矩阵; 为矢量R_ij的单位矢量; u_i和u_j分别是粒子i和j的速度.等号右侧第一项描述粒子承受流体阻碍并与运动速度方向相反; 等号右侧第二项则指出粒子j的运动会诱导粒子i运动, 其Stokes流体力的强度正比于a/R_ij.由此可知, R_ij距离越远, 该作用力对粒子j的影响越小.另外, 粒子半径越大, 周围粒子受到干扰和波动的程度就越大.大多数仿真仅考虑式(4)等式右侧第一项, 忽略了运动粒子使流体流动.将式(4)植入PTF模块, 结果如图 3所示.将粒子i初始速度设定为由左向右的直线运动, 在缺少等式右侧第二项的条件下, 粒子j位置和速度无任何变化, 如图 3a所示.考虑式(4)右侧第二项粒子的相互作用, 在相同初始条件下可得粒子i和j的运动轨迹和相对位置(见图 4), 发现粒子i的运动会使上方粒子j发生小范围移动, 两粒子有相互靠近的趋势, 这主要是由于两粒子内侧流速高于外侧流速, 导致内部流体压力低于外部, 使粒子相互靠拢.
图 3(Fig. 3)

图 3 忽略流体动力相互作用的粒子i和j运动仿真Fig.3 Motion simulation of particle i and j with the effect of hydrodynamic interaction neglected (a)—粒子i和j的轨迹; (b)—随时间变化的粒子i和j位置.

图 4(Fig. 4)

图 4 考虑流体相互作用的粒子i和j的运动仿真Fig.4 Motion simulation of particle i and j with the effect of hydrodynamic interaction considered (a)—粒子i和j的轨迹; (b)—随时间变化的粒子i和j位置

为保证粒子之间不会发生重叠等不合理现象, 应力仿真需要考虑短程排斥力F_rep, 排斥力的作用范围在埃米到纳米之间.因为粒子半径为微米量级, 所以排斥力选择硬球模型^[16], 短程相互作用势能公式如下:

(5)

由排斥能表达式(5)所知, 该数学式不具有物理连续性, 编程求解难以实现.因此排斥力的数学式采用指数衰减形式表达:

(6)

式中: F₀为排斥力常数, 数值大小与介电泳力相当; 1/κ为排斥力的作用范围; δ为两粒子之间的最小间距.忽略布朗运动力的影响, 粒子i(或粒子j)运动满足牛顿第二定律, 如式(7)所述:

(7)

根据式(7)等式右侧合力, 可对粒子受电场作用的次级介电泳力进行模拟.考虑仿真时间步长, 其与激励电压V₀成反比, 表示电场强度增大, 则为保证运算的收敛性, 使介电泳力和其他作用力在同一数量级, 时间离散度必须降低.
2 结果和讨论如图 1所示, 设定两粒子的初始距离R_ij为60 μm, 夹角θ为60°, 将其输入COMSOL软件中PTF模块, 并引入AC/DC电场频率求解器, 随后采用时域求解器可得粒子受次级介电泳力的运动轨迹, 如图 5所示.从图 5两个粒子的运动轨迹可发现, 其与文献[8, 10-11]具有很高的相似度.因为粒子i和粒子j具有较好的空间对称性, 且粒子半径和极化强度K_CM相同, 两者的运动轨迹也保持较好的旋转对称特性.本文的激励电场为10 V, 频率120 kHz, 根据图 2得到两种不同电导率的粒子K_CM分别为0.42和-0.45.
图 5(Fig. 5)

图 5 两粒子相对运动仿真结果Fig.5 Simulated result of two particles moving in opposite direction

不同激励频率下, 粒子i, j的位置固定不变, 它们的运动轨迹基本相同, 并且都与电场方向平行成链; 同时可得粒子间距离R_ij不同条件下的速度曲线, 如图 6所示.初始距离为60 μm时, 粒子i初始速度为0.随着粒子间距离的减小, 速度逐渐增大, 该仿真结果也与文献[8]一致.粒子分别受到正、负介电泳力时，其运动速度都表明, 当R_ij小于粒子直径2a后, 受排斥力影响速度开始急速降低.由于粒子间距离还未达到排斥力作用范围(纳米量级), 因此粒子速度的降低主要归因于粒子形成的流场阻力.
图 6(Fig. 6)

图 6 粒子间距离与速度的关系Fig.6 Relationship between the distance of two particles and velocity

为了进一步验证本文方法的有效性, 引入粒子k, 假设三种可能的粒子初始分布分别为三角形和团簇Ⅰ型和Ⅱ型, 如图 7、图 8所示.图 7中粒子初始位置以等边三角形分布, 仿真结果显示粒子j和粒子k最终相互吸引, 而粒子i被排斥到仿真区域正上方.该仿真结果与文献[14]结论一致.
图 7(Fig. 7)

图 7 初始正三角形分布的3个粒子轨迹Fig.7 Trajectories of three particles started from an equilateral triangle distribution

图 8(Fig. 8)

图 8 单粒子与粒子链的运动轨迹Fig.8 Trajectories of single particle and particle chain (a)—成链粒子与电场方向平行; (b)—成链粒子与电场方向垂直.

图 8a中粒子j和粒子k初始相连, 构成粒子链且平行于水平方向电场.粒子j和粒子k同时向粒子i靠拢, 最终3个粒子形成新的粒子链; 该结果主要是因为成链粒子形成的团簇结构对均匀电场的扭曲效应更强.图 8b中粒子j和粒子k初始相连, 构成粒子链且垂直于水平方向电场.单个粒子i向着成链粒子移动, 但是粒子j和粒子k难以同步向上与粒子i靠拢, 因为它们电偶极矩与电场方向垂直, 粒子k被排斥到求解域外, 而粒子j与粒子i继续相互吸引形成短链结构.此时粒子i和粒子j构成的短链型水平高度(y>0)高于图 8a长链(y＜0).
COMSOL的粒子追踪模块可以无编程、灵活地进行两个或多个粒子介电泳效应受力分析, 同时也能进行金属电极阵列的粒子运动模拟.本文对三角形金属电极进行了建模和仿真, 如图 9、图 10所示.由于电极提供的电场具有非均匀特性, 需要考虑电极形成的介电泳力, 因此将式(1)的介电泳力引入粒子追踪模块.激励电压和频率大小保持不变, 根据图 2两种电导率不同的粒子频谱, 分别进行受力仿真分析.同时, 针对模型是否考虑偶极子静电互动力进行对比研究.
图 9(Fig. 9)

图 9 电导率为5.0×10-3S/m粒子受正介电泳力影响的动力学仿真结果Fig.9 Simulated results of positive dielectrophoretic force effect on particles with conductivity of 5.0×10-3S/m (a)—t=0 s, 粒子的初始位置; (b)—忽略次级正介电泳效应的粒子分布; (c)—考虑次级正介电泳效应的粒子分布.

图 10(Fig. 10)

图 10 电导率为8.0×10-4S/m粒子受负介电泳力影响的动力学仿真结果Fig.10 Simulated results of negative dielectrophoretic force effect on particles with conductivity of 8.0×10-4S/m (a)—忽略次级负介电泳效应的粒子分布; (b)—考虑次级负介电泳效应的粒子分布.

正、负介电泳力分别将粒子驱使到电极附近或电极中心位置, 总体分布基本相似, 但是观察细节可发现, 次级介电泳效应能够使粒子排列得更加有序.正介电泳力吸引粒子到电极尖端位置, 介电泳力大于偶极子互动力, 粒子难以成链.远离电极的粒子受到的介电泳力较弱, 次级介电泳力使它们彼此成链, 但并非如图 5和图 8a完全水平排列, 而与弯曲电场方向一致.负介电泳力驱使粒子远离较强电场的电极区域, 它们的链簇数量比正介电泳更多.
3 结论1) 基于COMSOL Multiphysics 5.3a的AC/DC模块和PTF模块通过植入介电泳力表达式, 可以实现粒子的介电泳效应模拟, 并且粒子运动特征与先前的研究结果保持较大的相似性.
2) 为了更准确地预测实验观测结果, 粒子受正、负介电泳效应的数值模型应该考虑次级介电泳效应.
参考文献

[1]	Viefhues M, Eichhorn R. DNA dielectrophoresis: theory and applications a review[J]. Electrophoresis, 2017, 38(11): 1483-506. DOI:10.1002/elps.201600482
[2]	Zhang H Q, Chang H L, Neuzil P. DEP-on-a-chip: dielectrophoresis applied to microfluidic platforms[J]. Micromachines, 2019, 10(6): 423-444. DOI:10.3390/mi10060423
[3]	Albrecht D R, Tsang V L, Sah R L, et al. Photo-and electropatterning of hydrogel-encapsulated living cell arrays[J]. Lab on a Chip, 2005, 5(1): 111-118. DOI:10.1039/b406953f
[4]	Ho C T, Lin R Z, Chang W Y, et al. Rapid heterogeneous liver-cell on-chip patterning via the enhanced field-induced dielectrophoresis trap[J]. Lab on a Chip, 2006, 6(6): 724-734. DOI:10.1039/b602036d
[5]	Brimmo A T, Menachery A, Qasaimeh M A. Microelectrofluidic probe for sequential cell separation and patterning[J]. Lab on a Chip, 2019, 19(24): 4052-4063. DOI:10.1039/C9LC00748B
[6]	Velev O D, Gangwal S, Petsev D N. Particle-localized AC and DC manipulation and electrokinetics[J]. Annual Reports on the Progress of Chemistry C: Physical Chemistry, 2009, 105(1): 213-246.
[7]	Ai Y, Qian S Z. DC dielectrophoretic particle-particle interactions and their relative motions[J]. Journal of Colloid and Interface Science, 2010, 346(2): 448-454. DOI:10.1016/j.jcis.2010.03.003
[8]	Ai Y, Zeng Z P, Qian S Z. Direct numerical simulation of AC dielectrophoretic particle-particle interactive motions[J]. Journal of Colloid and Interface Science, 2014, 417(1): 72-79.
[9]	Hossan M R, Dillon R, Roy A K, et al. Modeling and simulation of dielectrophoretic particle-particle interactions and assembly[J]. Journal of Colloid and Interface Science, 2013, 394(15): 619-629.
[10]	Kang S. Dielectrophoretic motions of multiple particles under an alternating-current electric field[J]. European Journal of Mechanics B: Fluids, 2015, 54: 53-68. DOI:10.1016/j.euromechflu.2015.06.014
[11]	Xie C C, Chen B, Liu L, et al. Iterative dipole moment method for the interaction of multiple dielectrophoretic particles in an AC electrical field[J]. European Journal of Mechanics B: Fluids, 2016, 58: 50-58. DOI:10.1016/j.euromechflu.2016.03.003
[12]	Zhao Y, Brcka J, Faguet J, et al. Elucidating the DEP phenomena using a volumetric polarization approach with consideration of the electric double layer[J/OL]. Biomicrofluidics, 2017, 11(2)[2020-12-18]. https://doi.org/10.1063/1.4979014.
[13]	Lo Y J, Lei U. Measurement of the real part of the Clausius-Mossotti factor of dielectrophoresis for Brownian particles[J]. Electrophoresis, 2020, 41(1/2): 137-147.
[14]	Hu S, Fu R R. Expanding the flexibility of dynamics simulation on different size particle-particle interactions by dielectrophoresis[J]. Journal of Biological Physics, 2019, 45(1): 45-62. DOI:10.1007/s10867-018-9514-7
[15]	Liu W Y, Ren Y K, Xue R, et al. On ion transport regulation with field-effect nonlinear electroosmosis control in microfluidics embedding an ion-selective medium[J]. Electrophoresis, 2020, 41(10/11): 778-792.
[16]	Lin Y, Amberg G, Aldaeus F, et al. Simulation of dielectrophoretic motion of microparticles using a molecular dynamics approach[C/OL]//Fourth International Conference on Nanochannels, Microchannels and Minichannels. Limerick, Ireland, 2006[2020-12-18]. https://www.researchgate.net/deref/http%3A%2F%2Fdx.doi.org%2F10.1115%2FICNMM2006-96095.