微粒受介電泳效應的動力學建模和仿真

胡 晟， 王 克， 蔡 露

(東北大學秦皇島分校 控制工程學院, 河北 秦皇島 066004)

介電泳(dielectrophoresis, DEP)作為一種微/納粒子操控、分離、富集和輸運技術，已經(jīng)受到國內、外學者的廣泛關注和學習[1-4].在非均勻電場作用下，粒子電極化的應力不平衡，因此會發(fā)生沿電場強度梯度方向的定向運動.微粒向電場強度較強區(qū)域運動，該介電泳力稱為正介電泳力(positive DEP, pDEP).反之，負介電泳力(negative DEP, nDEP)則驅使粒子運動到電場強度較弱的空間位置.上述正、負介電泳力的選取主要與粒子或溶液的介電常數(shù)、電導率和激勵頻率有關.伴隨數(shù)學物理方法和計算方法學的飛速發(fā)展，目前介電泳的理論研究也取得了長足的進步.根據(jù)經(jīng)典電偶極矩法[5-6]和Maxwell應力張量法[7-8]關于介電泳力的物理描述，粒子運動軌跡和拓撲形態(tài)的仿真結果，都有效地證實了實驗觀測和電極設計.

然而，從當前的數(shù)學建模和仿真分析來看，粒子承受介電泳力作用的動力學研究仍然存在一些問題需要進一步研究和解決，主要由于編程的復雜以及后期代碼的維護工作.特別對于粒子承受排斥力與介電泳吸附效應等問題時，粒子相互作用力的受力矩陣也導致了求解精度和迭代次數(shù)的成倍增加，而關于代碼的優(yōu)化和電場偏微分方程的求解仍需要一個更好的方法和工具進行優(yōu)化求解.基于COMSOL Multiphysics 5.3a有限元分析軟件通過各種內置的方法已經(jīng)極大解決了物理域的網(wǎng)格離散和偏微分方程求解.通過Laplace微分方程和邊界條件的設定，可以快速得到電極芯片內部的電壓分布[9-10]，從而掌握電場變化的方向和強度.COMSOL Multiphysics 5.3a軟件提供了粒子跟蹤模塊(particle tracing module, PTM)，該模塊極大擴大了場與粒子相互作用的研究內涵.可以根據(jù)求解得到的物理場進行粒子運動特征的分析和學習.Zhao等[11]最早采用PTM和AC/DC Module結合的方法進行了粒子受介電泳效應的動力學研究.但是他們研究所提供的相關信息相對較少，并且介電泳粒子成鏈效應仍存在研究不足等缺點.本文重新對粒子位于介電泳芯片的粒子受力運動問題進行了仿真和研究，豐富了COMSOL Multiphysics 5.3a軟件在介電泳粒子操控領域的理論研究，使該軟件提供的AC/DC和PTM模塊可以較好解決現(xiàn)有Velocity-Verlet或ODE算法編程困難、執(zhí)行效率低的問題.

1 原理與物理場模型

微粒的介電泳力產(chǎn)生機理歸因于電場分布的非均勻特性，從空間角度觀察粒子運動可分為全局和局部效應.全局介電泳效應主要是電極產(chǎn)生的電場在較大尺度空間誘導粒子輸運和富集.但是局部介電泳效應，則是微粒迫使均勻電場在局部或粒子附近發(fā)生扭曲，進一步誘導周圍粒子移動成鏈的過程.在無損耗介質溶液中，電極產(chǎn)生的介電泳力表達式[1]為

FDEP=(p·)E.

(1)

式中：E為外界激勵電場強度；為哈密頓算子；p為粒子的電偶極矩，其表達式[5,12]為

p=4πε0εma3βE.

(2)

式中：ε0和εm分別為真空介電常數(shù)和溶液的相對介電常數(shù)；a為粒子的半徑；β=(εp-εm)/(εp+2εm)，代表粒子與溶液的電極化率，εp是粒子的相對介電常數(shù).

然而，粒子承受的局部介電泳力與它們之間的電偶極矩動量密切相關，表達式[5，13]為

(3)

式中：pi和pj分別為i粒子和j粒子各自的電偶極矩，可由式(2)計算得出；rij是i粒子和j粒子的單位方向矢量；dij是標量，代表了i粒子和j粒子之間的歐拉距離.需要注意，COMSOL Multiphysics 5.3a的PTM中，關于Particle-Particle Interaction應力項提供了dest子函數(shù)用于求解剩余粒子的空間位置.通過Variables對話框編寫上述函數(shù)解析式作為內聯(lián)函數(shù)，以便于應力求解所需.

另外粒子相互碰撞和遷移還應包含相互排斥力，防止粒子發(fā)生重合現(xiàn)象，表達式[13-14]為

(4)

式中：δmin為兩個粒子之間的最短距離；1/κ為排斥力作用范圍；Frep0為標量值，通常等于靜電力大小.

粒子在運動過程中還受Stokes拖曳力，并且粒子所受的重力分別滿足式(5)和式(6).

Fs=6πηau，

(5)

(6)

式中：u為粒子移動速度；η為溶液黏度；ρp和ρm分別為粒子和溶液的密度；g則是重力加速度.對于i和j粒子的應力模型，如圖1所示.此外，粒子下方的芯片底面也需要設置排斥力，防止粒子運動至求解域外.

(7)

式中：Fw0是墻壁給予粒子的最大排斥力；ri是i粒子的空間位置(xi,yi,zi).式(4)和式(7)的初始值Frep0，F(xiàn)w0和1/κ分別取值2.3×10-12N，1×10-10N和0.01.同時，求解域中剩余5個墻壁面，設定Freeze邊界條件，有利于粒子的觀測和學習.

圖1 微粒受力示意圖

2 結果和討論

2.1 均勻電場

首先進行平行極板產(chǎn)生均勻電場的求解，圖1中，左、右邊界分別設定+3.5 V和接地邊界條件，其余邊界為絕緣邊界條件.通過求解Laplace方程，即式(8)和式(9)，可得到施加電壓后芯片內部的電場分布[15].

E=-V,

(8)

D=ε0εrE.

(9)

式中：V為電勢的標量場；εr為溶液與粒子的相對介電常數(shù).相關的仿真參數(shù)如表1所示.根據(jù)上述應力表達式的簡要說明，采用PTM獲取AC/DC的電場數(shù)值進行粒子應力模擬.因此得到相關的粒子運動結果，如圖2所示.在均勻電場中，粒子受局部介電泳效應將發(fā)生定向遷移，最終形成一定長度的珍珠鏈結構.同時，PTM模塊的靈活可編程特點，也能進行不同尺寸粒子在均勻電場應力下的運動特征分析.圖3為粒子半徑為 5 μm 和2.5 μm時的仿真結果.從圖中發(fā)現(xiàn)，混合了不同尺寸粒子很難實現(xiàn)珍珠鏈結構，主要歸因于偶極子之間空間對稱性不匹配.

表1 仿真模型關鍵參數(shù)

圖2 均勻電場下，35個粒子且半徑為5 μm的仿真結果

圖3 均勻電場下，15個半徑為5 μm和20個半徑為2.5 μm的粒子仿真結果

2.2 非均勻電場

設置金屬電極產(chǎn)生非均勻電場，誘導粒子向電場較強或較弱方向運動.激勵電極如圖4a所示，右側負電極接地.電極電壓仍為+3.5 V.聚苯乙烯粒子εp為2.25,受負介電泳力運動到電場強度較弱位置，并且局部形成短鏈結構，如圖4c所示.然后對粒子表面進行化學修飾,假設微粒εp為120，高于水溶液80，電極針尖附近的粒子將運動到其附近，如圖4d所示.但是遠離針尖的粒子所受正介電泳力相對較弱，則逐漸沉淀至芯片底部.

另外采用電極陣列型結構計算區(qū)域內35個半徑5 μm的微粒運動結果，如圖5所示.負介電泳效應時，靠近邊緣的粒子運動到電極與電極之間電場較弱區(qū)域.反之正介電泳力驅使粒子到電極邊緣位置，該模擬結果也與文獻[16-18]實驗基本類似.

圖4 非均勻電場下，半徑為5 μm的35個粒子承受負、正介電泳力仿真結果

圖5 陣列型電極結構中，粒子承受負、正介電泳力的仿真結果

3 結 論

1) 本文采用了COMSOL Multiphysics 5.3a軟件，通過求解Laplace方程獲得求解域內的電場分布；后期耦合了粒子追蹤模塊，模擬芯片內微粒受介電泳效應的運動行為.

2) 仿真結果表明，本文方法能有效降低計算程序的繁瑣程度，提高動態(tài)模擬的人機可視化效果.較好解決現(xiàn)有Velocity-Verlet或ODE算法編程困難、執(zhí)行效率低的問題.