微粒受介電泳效應(yīng)的動(dòng)力學(xué)建模和仿真

胡 晟， 王 克， 蔡 露

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

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

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

1 原理與物理場(chǎng)模型

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

FDEP=(p·)E.

(1)

式中：E為外界激勵(lì)電場(chǎng)強(qiáng)度；為哈密頓算子；p為粒子的電偶極矩，其表達(dá)式[5,12]為

p=4πε0εma3βE.

(2)

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

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

(3)

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

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

(4)

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

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

Fs=6πηau，

(5)

(6)

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

(7)

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

圖1 微粒受力示意圖

2 結(jié)果和討論

2.1 均勻電場(chǎng)

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

E=-V,

(8)

D=ε0εrE.

(9)

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

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

圖2 均勻電場(chǎng)下，35個(gè)粒子且半徑為5 μm的仿真結(jié)果

圖3 均勻電場(chǎng)下，15個(gè)半徑為5 μm和20個(gè)半徑為2.5 μm的粒子仿真結(jié)果

2.2 非均勻電場(chǎng)

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

另外采用電極陣列型結(jié)構(gòu)計(jì)算區(qū)域內(nèi)35個(gè)半徑5 μm的微粒運(yùn)動(dòng)結(jié)果，如圖5所示.負(fù)介電泳效應(yīng)時(shí)，靠近邊緣的粒子運(yùn)動(dòng)到電極與電極之間電場(chǎng)較弱區(qū)域.反之正介電泳力驅(qū)使粒子到電極邊緣位置，該模擬結(jié)果也與文獻(xiàn)[16-18]實(shí)驗(yàn)基本類似.

圖4 非均勻電場(chǎng)下，半徑為5 μm的35個(gè)粒子承受負(fù)、正介電泳力仿真結(jié)果

圖5 陣列型電極結(jié)構(gòu)中，粒子承受負(fù)、正介電泳力的仿真結(jié)果

3 結(jié) 論

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

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