基于拓撲優化的電滲流微混合器電極

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

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

1 引言

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

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

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

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

2 模型建立

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

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

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

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

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

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

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

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

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

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

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

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

式中β與ξ是投影參數。

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

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

3 敏度分析

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

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

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

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

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

（3）計算優化目標函數的伴隨導數；

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

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

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

4 數值實驗

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

模型數值實驗的具體參數如下：電解質濃度ρ=1×，動力黏度η=1×10-3Pa·s，介電常數εr=80.2，電導率σ=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）進行迭代計算求解。

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

圖2 是拓撲優化模型設計變量的演化過程。當迭代300 次后，電滲流微混合器與混合單元表面電極版圖如圖2（f）所示，黑色區域表示優化設計后的電極布局。圖3 是目標函數的收斂過程，其逐漸收斂平穩的趨勢可表示該拓撲優化方法的穩定性。

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

圖3 拓撲優化目標函數的收斂過程Fig.3 Convergence history of optimization objective function values

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

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

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

4.1 多混合單元設計

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

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

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

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

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

4.2 電極電壓的影響

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

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

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

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

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

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

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

4.4 佩克萊數與雷諾數的影響

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

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

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

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

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

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

5 結論

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