孫建文，張健宇，李博文，鄧永波*

（1.中國科學(xué)院 長春光學(xué)精密機械與物理研究所，吉林 長春 130033；2.中國科學(xué)院大學(xué)，北京 100049）

1 引言

微流控系統(tǒng)具有高通量、分析時間短、樣品消耗率低等優(yōu)勢，已廣泛應(yīng)用在生物、化學(xué)和醫(yī)藥等領(lǐng)域［1-3］。微混合器是微流控技術(shù)中主要的功能器件，根據(jù)驅(qū)動方式不同可分為被動混合器和主動混合器。與被動微混合器通過改變微通道幾何形狀或添加幾何障礙物來增強混合效果不同［4-5］，主動微混合器通過施加外部能量來擾動樣品，具有更短的混合時間，并且操作者可根據(jù)具體情況對混合器進行開關(guān)操作［6］。

電滲驅(qū)動是主動微混合器中一種常用的驅(qū)動混合方式，其工作原理是通過對微通道流體施加外部電場產(chǎn)生電滲流動現(xiàn)象，達到不同濃度溶液混合的目的。Matsubara 等通過交錯的電極排布來增強非定常微流體的混合效果［7］。Seo 等研究了環(huán)形微混合器的電極位置對混合性能的影響［8］。Xiong 等利用收斂和發(fā)散原理設(shè)計一種菱形電滲微混合器，并通過改變菱形角度與組合方案，研究了電極電壓、電極結(jié)構(gòu)和雷諾數(shù)對混合性能的影響［9］。Farahinia 團隊對微通道壁面非均勻分布電極電滲微混合進行研究［10］，隨后將電滲流與壓力流相結(jié)合開展進一步研究［11］。Lv 等針對非牛頓流體與牛頓流體的微混合展開研究，驗證兩種流體混合效果與電壓頻率成正比，與入口流速成反比［12］。Chen 等提出一種能夠快速響應(yīng)并且充分混合的交流電場微混合器結(jié)構(gòu)，并對混合單元、電極對數(shù)等參數(shù)進行討論分析［13］。上述工作通過改變微混合器結(jié)構(gòu)、電壓、雷諾數(shù)以及電極排布方式等影響因素提高了樣品的混合效率，然而缺乏對電滲流微混合器電極版圖的結(jié)構(gòu)設(shè)計研究。

電滲流是在雙電層、外部電場與流場相互作用下生成的。因此，施加外電場的電極版圖對控制電滲流起關(guān)鍵作用。為提升電滲流微混合器的混合性能，可對微混合器的電極版圖進行拓撲優(yōu)化設(shè)計。拓撲優(yōu)化能夠在缺乏設(shè)計經(jīng)驗的情況下，找到滿足設(shè)計目標與約束條件的材質(zhì)最佳分布，從而確定結(jié)構(gòu)的最優(yōu)拓撲。

針對上述問題，本文研究了電滲流電極的拓撲優(yōu)化方法，以實現(xiàn)兩種不同濃度溶液的完全混合為設(shè)計目標，建立電滲流微混合器模型，進而確定微混合器的電極版圖，并通過改變電壓、混合單元設(shè)計方案等影響因素驗證微混合器的混合性能，得到混合效果最佳的設(shè)計方案。

2 模型建立

在微混合器中，電滲是由與電解質(zhì)溶液接觸的微通道壁表面上凈電荷集聚所產(chǎn)生的。電荷聚集在壁面形成雙電子層［14］，其厚度遠小于微流體的特征尺寸（通常為10 nm），可將其近似為滑移壁面。因此，電滲微流體可以用具有滑移壁面邊界條件的Navier-Stokes 方程來描述，其中滑移速度滿足Helmholtz-Smoluchowski 公式：

式中：u與p分別是流體速度和壓力；V是電極受到的外部施加電勢；ρ和η分別是電解質(zhì)溶液的密度以及動力黏度；μeo=-εrε0(ξ0/η)為 電滲遷移率，εr，ε0和ξ0分別是相對介電常數(shù)、自由介電常數(shù)以及zeta 電位；n是?Ω上的向外單位法線。

如圖1 所示，電滲微混合器示意圖由入口通道、混合單元（混合設(shè)計區(qū)域）和出口通道組成。其中，Ω為微混合器設(shè)計域，Γia∪Γib=Γi為微混合器入口，Γo為出口，Γw和Γd是微混合器壁面。電極版圖在Γd所示混合單元上下壁面進行設(shè)計，上表面為施加電壓的陽極，下表面為接地陰極，V0為陽極所受電壓。微混合器模型滿足Γi∪Γw∪Γd∪Γo=?Ω。電滲流微混合器特征尺寸l∝100 μm，混合單元半徑R=n·l，其中n為混合單元的半徑系數(shù)。

圖1 電滲微混合器示意圖Fig.1 Schematic diagram of electroosmotic micromixer

在電滲流微混合器中，影響混合效果的兩個主要因素為擴散與對流。因此，兩種不同濃度流體的電滲流混合可以用對流擴散方程進行描述：

式中：c是濃度，D是質(zhì)量擴散系數(shù)。如圖1 所示，兩種不同濃度的溶液進入微混合器，其濃度設(shè)定為：

式中：ci是電滲微混合器入口處的已知濃度分布。微混合器壁面與出口的濃度場邊界條件為：

電滲流混合器模型所受的外部電勢可用Laplace 方程-?·(σ?V)=0 描述，σ是電導(dǎo)率。如圖1 所示，在電滲流微混合器模型中，壁面邊界被分為兩部分Γw與Γd。Γw為電絕緣邊界，Γd為混合單元電極設(shè)計區(qū)域。為了區(qū)分模型中的兩種邊界，將它表示為取值在［0，1］之間的物理密度的插值，其中0 和1 分別表示電勢與電絕緣邊界。因此，電極的設(shè)計區(qū)域的邊界條件可以表示為電絕緣與電勢的插值：

式中：V0施加電極電壓；α為該方程懲罰因子［15］。

式中：γfp，αmax與q分別是物理密度變量、懲罰參數(shù)，以及用于調(diào)節(jié)懲罰凸度的參數(shù)。當模型的物理密度為0 時，在式（6）中，αmax取值足夠大，以確保式（5）中的(V-V0)占據(jù)主導(dǎo)，使等式左側(cè)-(σ?V)·n項被忽略為0，即電極受到的電壓V=V0，此時壁面為電勢邊界；而當模型物理密度為1 時，懲罰因子為0，式（5）可以轉(zhuǎn)化為-(σ?V)·n=0，即電絕緣邊界。基于數(shù)值實驗，αmax與q分別選擇1×103與1×10-3。此外，在電滲流模型的入口與出口也施加為電絕緣邊界。因此，電滲微混合器中電極受到的外部電勢可以描述為：

物理密度變量是由定義在Γd上的設(shè)計變量通過Helmholtz 濾波器調(diào)整計算的［16］：

式中：γ是設(shè)計變量；γf是經(jīng)過過濾的設(shè)計變量；?s在Γd上局部坐標系定義的梯度；r是濾波半徑，用于控制電極版圖的特征尺寸；ns是在Γd向外的單位法線。在過濾計算后，使用閾值法去除過濾后的設(shè)計變量中0 與1 之間的中間值，并導(dǎo)出物理密度：

式中β與ξ是投影參數(shù)。

此時，電滲流模型的電極設(shè)計已轉(zhuǎn)變成求解0-1 的最優(yōu)材質(zhì)分布問題。本文的目標函數(shù)為模型出口處的實際平均濃度與預(yù)期濃度之間的歸一化最小二乘方差。電滲流微混合器的設(shè)計目標是通過找到合理的電極版圖設(shè)計，達到不同濃度流體的有效混合。因此，模型變分問題可構(gòu)造為：

約束條件下：式中Ψ0是出口處所測得濃度與理想目標濃度之間的最小二乘方差。

3 敏度分析

本文使用數(shù)值迭代的方法解決模型的變分問題。根據(jù)從等式的變分方程中提取的梯度信息，對設(shè)計變量進行最優(yōu)化求解?；趯α鲾U散方程、Navier-Stokes 方程、Laplace 方程和Helmholtz 濾波器的伴隨方程的弱形式，通過約束優(yōu)化問題的偏微分方程對變分問題進行伴隨分析：

式中：λf與λQ是拉格朗日乘子系數(shù)，用于在模型入口出口施加邊界條件；λc為入口濃度邊界條件的拉格朗日乘子系數(shù)；ca，ua，pa，Va，γfa，λfa，λQa和λca分別是其所對應(yīng)狀態(tài)變量的伴隨變量；是其對應(yīng)伴隨變量的測試函數(shù)；H(Ω)和L2(Ω)分別是定義在Ω的一階Hilbert 空間與二階 Lebesgue 可積空間 ；H-1/2(Γw∪Γd)是跡空間H1/2(Γw∪Γd)的對偶空間；H(Γd)是定義在Γd的一階Hilbert 空間；R 表示實數(shù)域。為了求解Navier-Stokes 方程和對流擴散方程，采用有限元法、廣義最小二乘法和流線逆風Petrov-Galerkin（SUPG）法。Ne是離散化Ω的有限單元數(shù)；τgls和τsupg為穩(wěn)定參數(shù)：

式中：hi表示的是在第i個計算域Ωi的大小。式（10）變分之后的一階伴隨導(dǎo)數(shù)為：

（1）用當前設(shè)計變量求解PDE 約束；

（2）基于PDE 約束的解求解伴隨方程；

（3）計算優(yōu)化目標函數(shù)的伴隨導(dǎo)數(shù)；

（4）使用移動漸近線法計算并更新設(shè)計變量；

（5）檢查迭代結(jié)果是否滿足收斂準則，滿足則變分計算結(jié)束；若不滿足，將返回（1），繼續(xù)迭代計算。

在電滲微混合器模型變分問題的收斂過程中，連續(xù)5 次迭代的目標函數(shù)變化滿足第k次迭代步驟中的，或者達到已設(shè)定的最大迭代次數(shù)300 時，迭代結(jié)束。其中Ψk是第k次迭代的目標函數(shù)。在迭代計算中，式（9）中的閾值參數(shù)ζ設(shè)為0.5；投影參數(shù)β初始值設(shè)置為1，每40次迭代后加倍，直到最大值128。

4 數(shù)值實驗

如圖1所示，不同濃度的溶液從微通道入口流入到混合單元腔體內(nèi)進行混合，并從微通道出口流出。在微混合器中實現(xiàn)混合的物理過程主要是擴散與對流，二者可通過改變流經(jīng)微混合器流體的佩克萊數(shù)（Pe）與雷諾數(shù)（Re）進行調(diào)節(jié)。其中，U為模型入口處的平均速度，D為質(zhì)量擴散系數(shù)。當Pe?1 時，影響混合的主要因素為對流，反之即為擴散。而流體的對流強度是由Re表征的，。當Re?1 時，對流在微混合器中起主導(dǎo)作用，反之黏性作用為主要作用。

模型數(shù)值實驗的具體參數(shù)如下：電解質(zhì)濃度ρ=1×，動力黏度η=1×10-3Pa·s，介電常數(shù)εr=80.2，電導(dǎo)率σ=0.12 (Ω·m)-1，Zeta 電勢ζ0=-0.1 V。電滲流微混合器入口和出口微通道的特征尺寸l=4×10-4m，混合單元半徑R=n·l。

微混合器選取V0=4 V，n=2，l=4×10-4m，Pe=1×103和Re=1，對式（10）進行迭代計算求解。

為直觀展示模型的混合性能，將目標函數(shù)Ψ作為混合評價指數(shù)，其數(shù)值與混合效果呈負相關(guān)關(guān)系。當混合評價指數(shù)小于0.05 時，兩種不同濃度的溶液已完全混合［17］。

圖2 是拓撲優(yōu)化模型設(shè)計變量的演化過程。當?shù)?00 次后，電滲流微混合器與混合單元表面電極版圖如圖2（f）所示，黑色區(qū)域表示優(yōu)化設(shè)計后的電極布局。圖3 是目標函數(shù)的收斂過程，其逐漸收斂平穩(wěn)的趨勢可表示該拓撲優(yōu)化方法的穩(wěn)定性。

圖2 電滲流微混合器物理密度演化過程Fig.2 Physical density evolution of electroosmotic flow micromixer

圖3 拓撲優(yōu)化目標函數(shù)的收斂過程Fig.3 Convergence history of optimization objective function values

為驗證電極拓撲優(yōu)化微混合器的混合性能，將未進行電極設(shè)計的微混合器與電極拓撲優(yōu)化微混合器的混合效果進行對比，如圖4 所示。結(jié)果表明，未進行電極設(shè)計的微混合器混合作用主要以擴散為主，在低雷諾數(shù)下只有較弱的平流，致使微混合器的混合性能大幅降低；而經(jīng)過電極拓撲優(yōu)化的電滲微混合器可對不同濃度流體進行有效混合。

圖4 微混合器濃度截面對比Fig.4 Concentration cross section comparison of micromixer

本文考慮影響微混合器混合性能的因素，通過改變混合單元設(shè)計方案、施加電壓、混合單元特征尺寸、雷諾數(shù)和佩克萊數(shù)等參數(shù)，研究它們對電滲流微混合器混合效果的影響，驗證電極拓撲優(yōu)化方法的有效性。

4.1 多混合單元設(shè)計

由于電滲流微混合器的混合單元體積較大，在混合溶液時，存在混合不充分的現(xiàn)象。因此采用多混合單元（N）的設(shè)計方案，通過增加混合驅(qū)動路徑的長度來提高溶液的混合性能［13，18］。

圖5 為N=9 時微混合器電極版圖和不同混合單元方案（N=1，…，9）的出口濃度分布。從圖5 中N=9 模型截面濃度分布可知，不同濃度溶液進入微混合器后，溶液流經(jīng)多個混合單元，混合效果隨混合驅(qū)動路徑的增加逐漸疊加，在出口處濃度分布呈現(xiàn)有效混合。從圖5 和圖6 可知，溶液混合效果隨混合單元數(shù)量的增加而增強，當N=9 時，混合評價指數(shù)達到0.05 以下，上述結(jié)果證明該方案可滿足不同濃度溶液的完全混合［17］。

圖5 多混合單元（N=9）的電極版圖和不同混合單元模型方案的出口濃度分布Fig.5 Electrode layout of mixing units（N=9）and outlet concentration distribution of different mixing unit model schemes

圖6 不同混合單元模型方案的混合評價指數(shù)Fig.6 Mixed evaluation index of different mixed unit model schemes

本文繼續(xù)討論微混合器混合效果的影響因素，從而驗證微混合器的混合性能。為便于后續(xù)數(shù)值模擬，采用單個混合單元微混合器方案進行分析討論。

4.2 電極電壓的影響

電極電壓是影響電滲流微混合的主要因素之一。這里在保持其他參數(shù)不變的情況下，將施加電壓V0分別設(shè)置成1，2，3，4 和5 V，討論施加電壓對電滲流微混合器混合效果的影響。圖7 為不同施加電壓下得到的電極設(shè)計布局和出口濃度分布，圖8 是分別對應(yīng)的不同施加電壓下的混合評價圖。根據(jù)圖7 和圖8 可知，較高電勢會產(chǎn)生更強的電場，在電滲流微混合器中有效引起擾動以及對流，從而得到更好的混合效果，有效提升微混合器的混合性能。

圖7 不同施加電壓電極版圖和出口濃度分布Fig.7 Electrodes obtained by applying different voltages and export concentration distribution

圖8 不同施加電壓混合評價指數(shù)Fig.8 Mixed evaluation of different applied voltages

4.3 混合單元特征尺寸的影響

將微混合器混合單元半徑R作為特征尺寸，R滿足R=n·l，其中l(wèi)=4×10-4m。圖9 為改變控制半徑系數(shù)n后所得的電極設(shè)計布局與出口濃度分布。從圖10 可知，混合評價指數(shù)與模型混合單元特征尺寸的關(guān)系，當n從1 逐漸增到3 時，微混合器的混合驅(qū)動區(qū)域逐漸增大，混合效果也隨之增強。

圖9 不同半徑電極版圖和出口濃度分布Fig.9 Electrode and outlet concentration distribution obtained by different radii

圖10 不同半徑混合評價指數(shù)Fig.10 Mixed evaluation of different radii

4.4 佩克萊數(shù)與雷諾數(shù)的影響

佩克萊數(shù)的物理意義是對流速率與擴散速率的比值。從圖11 和圖12 可知，電滲流微混合器在保持Re=1 的情況下，佩克萊數(shù)越低，混合效果越好。雷諾數(shù)是流體力學(xué)中表征黏性影響的相似準則數(shù)，在微通道流體流動中，不考慮湍流的影響時，雷諾數(shù)主要表征對流的程度。在保持佩克萊數(shù)不變時，低雷諾數(shù)時混合過程主要受到擴散的影響，拓撲優(yōu)化設(shè)計的電極布局所產(chǎn)生的電滲流動可以加強流體的橫向流動，以此加強混合效果，如圖13 所示。

圖11 不同佩克萊數(shù)下的電極版圖和出口濃度分布Fig.11 Electrode and outlet concentration distribution obtained by different Pelect numbers

圖12 不同佩克萊數(shù)混合評價指數(shù)Fig.12 Mixed evaluation of different Pelect numbers

圖13 不同雷諾數(shù)電極版圖與出口濃度分布Fig.13 Electrode and outlet concentration distribution obtained by different Reynolds numbers

當微混合器內(nèi)流體雷諾數(shù)逐漸升高時，流速的增加導(dǎo)致兩種不同濃度溶液在混合器中的混合時間減短，致使微混合器的混合效果不佳，如圖14 所示。

圖14 不同雷諾數(shù)混合評價指數(shù)Fig.14 Mixed evaluation of different Reynolds numbers

5 結(jié)論

本文基于拓撲優(yōu)化方法，以加強微流體混合效果為設(shè)計目標，通過在微混合器壁面施加電絕緣與電勢的插值的方法，建立電滲流微混合器模型，設(shè)計優(yōu)化電極版圖。數(shù)值模擬結(jié)果表明，電滲流微混合器的混合評價指數(shù)在多混合單元（N=9）設(shè)計下可達0.047，滿足完全混合的預(yù)期效果，從而驗證本文設(shè)計方法的合理性。本文構(gòu)造了由Navier-Stokes 方程、對流擴散方程所約束的變分問題，引入Helmholtz-Smoluchowsk 將雙電子層進行近似約束處理，采用基于拉格朗日乘子的伴隨方法和基于梯度信息的迭代方法對變分問題進行分析和求解。該電滲電極設(shè)計方法擺脫了經(jīng)驗設(shè)計的限制，在化學(xué)與生物等領(lǐng)域具有潛在的應(yīng)用價值。