胡晟, 吕江涛, 司光远
东北大学秦皇岛分校 控制工程学院, 河北 秦皇岛 066004
收稿日期：2018-01-24
基金项目：河北省自然科学基金资助项目(F2017501059)；辽宁省博士启动基金资助项目(20170520325)；中央高校基本科研业务费专项资金资助项目(N172304033)。
作者简介：胡晟(1984-)，男，云南景洪人，东北大学秦皇岛分校讲师，博士。

摘要：建立了双电层的离子分布模型, 基于经典Poisson-Boltzmann(PB)方程和改进型MPB(modified Poisson Boltzmann)方程对不同浓度和激励电压的离子分布进行了理论研究.结果发现在电压高于0.4 V, 且自由离子浓度小于10^-4 mol/L时, 双电层内部的扩散层厚度存在较大的误差.这直接导致了基于Debye长度模拟电渗流运动与实际观测不符, 主要因为Debye-Hückel公式具有线性关系不适用于仿真高电压条件下的电渗流运动.因此借助非线性MPB方程求解扩散层厚度, 更能精确得到正、负电极宽度为500 μm, 间距为25 μm, 在±1 V, 500 Hz电信号产生的最大电渗流速度为1 034.31 μm/s.
关键词：交流电渗流双电层PB方程MPB方程有限元法
Analysis and Modeling of Electro-osmosis Based on the Modified Poisson-Boltzmann Equation
HU Sheng, LYU Jiang-tao, SI Guang-yuan
School of Control Engineering, Northeastern University at Qinhuangdao, Qinhuangdao 066004, China
Corresponding author: HU Sheng, E-mail: husheng@neuq.edu.cn
Abstract: The ion distribution for different concentrations and applied voltage is studied by modeling the ion distribution in electric double layer(EDL)based on the Poisson-Boltzmann(PB)and modified Poisson-Boltzmann(MPB)equations. The results indicate that some errors exist in the thickness of diffusion in EDL if the applied voltage is more than 0.4 V and the ion concentration is less than 10^-4 mol/L. The simulated electro-osmotic flow by using Debye length is not in agreement with the practical observation result mainly because the Debye-Hückel equation is of linear relationship and not compatible with the calculating electro-osmotic flow under higher voltage conditions. By use of the MPB equation with non-linear characteristics, the maximal electro-osmotic flow generated by two co-planar metal electrodes with the width of 500 μm and the gap of 25 μm is 1 034.31 μm/s at ±1 V potential with the frequency of 500 Hz
Key words: AC electro-osmosiselectric double layer(EDL)PB equationMPB equationfinite element method
伴随微纳米加工技术的迅猛发展, 芯片实验室(lab-on-a-chip)已经广泛应用于体外细胞在线分析^[1]、抗癌药物基础研发^[2-3]、生态环境即时监控^[4]和临床医学诊断^[5]等生化领域.必须指出, 上述研究领域都离不开被控粒子的微操控与空间转移.其中, 在芯片实验室微尺度环境下实现微粒(如:细胞、DNA和蛋白质等)的输运, 主要依靠微流控技术, 如今该技术的实施方法相对较多, 如微流体作用力、光镊(optical tweezers)、光泳(photophoresis)、电泳(electrophoresis)、交流电渗流(AC electro-osmosis)和声泳(acoustophoresis)等^[6-7]物理方法.相较而言, 交流电渗流方法的使用领域较为广泛, 主要因为交流电渗流具有较好的电学可控性和较大空间尺度的输送特性.光镊能够精准操控微/纳粒子, 然而当被控对象的平移距离超过100 μm或更大时, 需要实时调整光源的工作距离, 由此提高了透镜对焦的研发成本.另外, 毛细电泳(>1 000 V/m)^[8]和光场(>100 kW/cm²)^[9]产生的热量都相对过高, 极容易造成生物样本的热创伤, 给后期细胞药物分析带来严重的干扰.
对比交流电渗流存在的优点后, 合理研究交流电渗流的流体运动规律, 以及微/纳粒子富集现象将是一件十分重要的基础性研究工作.然而, 当前关于交流电渗流的数学建模一般忽略纳米厚度的双电层(electric double layer, EDL)影响, 部分学者将芯片与溶液交界处边界条件近似性地等效为有效滑移速度(effective slip velocity).特别金属电极尺寸缩小到纳米量级, 或者激励电压过大都与实验观测相悖.主要因为双电层的阳离子与阴离子在Gouy-Chapman(GC)模型中被理想化地看作与几何形状无关的点电荷, 改进后的Poisson-Boltzman方程在GC模型中又增加了Stern层，更为客观地描述了阴离子和阳离子分布.根据溶液离子分布得到的Zeta电势能够更准确地描述电渗流速度场.
1 理论模型微槽与纳米管中的电渗流现象主要归因于外加激励电压驱动带电离子滑移运动加以实现, 相应的物理机理可采用库仑力来描述, 即F=E·q, F表示离子受力强度和方向, E为外界激励电场强度, q为离子的电量.一般情况, 传统的水溶液都保持电中性, 换句话说溶液内不存在一定密度ρ的自由电荷.然而从纳米尺度出发, 以金属电极为例, 施加一定电压信号, 在电极/溶液交界处会发生新的电化学平衡.电极排斥同性电荷, 使其远离电极表面；反之, 吸引异性电荷到其表面.Helmholtz最早提出了上述物理模型, 并且他将两边牢牢束缚的电荷形象地称为双电层(EDL), 如图 1a所示.Helmholtz模型描述了电极/电解液界面效应, 但是并没有考虑离子扩散与迁移等实际影响, 以及溶液偶极子与电极形成的相互作用.因此Gouy和Chapman改进了双电层模型, 使阴/阳离子具有空间可移动性, 对此保留了扩散和静电力分析特征.离子的分布称之为扩散层, 如图 1b所示.在GC模型中, 浓度为c_i的i离子被看作点电荷, 且满足Boltzmann分布, 如式(1)所示:

(1)

图 1(Fig. 1)

图 1 电极与电解液交界面阳离子与阴离子的分布示意图(虚线代表了阳极附近双电层的电势分布)Fig.1 Schematic diagram of the arrangement of anions and cations near the electrode/electrolyte interface (the dotted lines stand for the distribution of potential in the EDL near the positive electrode) (a)—Helmholtz模型；(b)—Gouy-Chapman模型；(c)—Gouy-Chapman-Stern模型.

式中:z_i, c_i∞分别代表i离子的化合价与本体浓度; T为溶液温度; e为单位电子电量(e=1.60×10^-19C); k_B为玻尔兹曼常数(k_B=1.381×10^-23 m²·kg·K^-1· s^-2).如图 1所示, ψ为电势, 随着逐渐远离电极, ψ逐渐趋于零.GC模型依托PB方程如式(2)所示:

(2)

其中:ε₀为真空介电常数(8.85×10^-12F/m); ε_r为相对介电常数.一般的对称性电解液(如, 氯化钾、氯化钠和硫酸镁), 它们的化合价满足z₁=-z₂=z, 又因为溶液呈现中性, 所以c₁=c₂=c_∞.N_A代表阿伏伽德罗常数(N_A=6.022× 10²³mol^-1).
当静电场能量远低于热能(, 一般室温T取值为298 K), 以单价对称离子溶液为例(z=1), ψ_T等于26 mV.那么Debye与Hückel将上述的GC表达式(2)近似为线性公式, 如下所示:

(3)

以一维模型描述溶液的电势变化, 那么求解常微分方程可得对应ψ的解析式如式(4)所示:

(4)

式中:ζ代表Zeta电势; k^-1常称作Debye厚度(λ_D), 用于估算双电层的厚度.从式(5)可以看出Debye长度的厚度仅与溶液的电化学特性有关:

(5)

后期, 学者Stern将Helmholtz与GC模型重新进行了整合, 修正了双电层的物理结构, 如图 1c所示.在GC模型基础上, 重新增加了Stern层(紧束层), 该物理层与Helmholtz层类似, 不具有自由可移动电荷.同时还保留了扩散层的离子游移现象, 那么相应的静电表达式应改写为

(6)

2 双电层仿真与分析通过上述双电层的物理模型描述, 于是对不同激励电压或浓度的溶液进行理论研究.求解式(2)~式(3)与式(6)需要有效的数值分析方法(如:改进欧拉法或龙格-库塔法), 但是必须指出当方程非一维模型时, 上述方程的求解将极其复杂.本文采用COMSOL Multiphysics软件, 基于有限元分析方法求解上述物理公式, 并且后续还用于求解Navior-Stokes流体方程, 仿真并得到相应的电渗流流场.双电层的求解域和边界条件可参考图 2, 求解域满足GC或GCS方程, 上、下边界条件满足绝缘边界条件(, n代表边界的单位法线矢量，s为边界); 左边界电压为ψ_E, 即激励电压; 右边界电压为零.温度T为298 K为恒定常数, 且相对介电常数ε_r=80.如图 3所示, 初始电压ψ_E(小于ψ_T≈26 mV)设为25 mV, 初始设定溶液为单价电解液(如KCl或NaCl), 那么化合价z⁺=|-z^-|=1.从曲线变化可发现, 激励电压相对较弱的条件下, Gouy-Chapman和Debye-Hückel(DH)变化趋势相同, 呈指数曲线衰减.但是正、负离子溶度较高的条件下(＜0.1 mol/L), 电势ψ的衰减速度更快, 因为大量的负离子迅速使正电势急剧减小.
图 2(Fig. 2)

图 2 求解域与边界条件物理方程示意图Fig.2 Schematic diagram of domain and boundary conditions in the physics equations

图 3(Fig. 3)

图 3 GC模型和DH模型的对比结果Fig.3 Comparison between Gouy-Chapman and Debye-Hückel physical models

根据式(5)可求解得到Debye厚度(λ_D)与离子浓度(c_0∞)的关系曲线, 如图 4所示.与上述ψ分布结论一致, 同样可证实溶液离子浓度过高, 双电层的厚度逐渐降低.同时根据式(1)Boltzmann的离子分布描述, 能够求解出GC和DH扩散层的厚度.必须指出, 针对图 4的激励电压增至0.4 V后, 此时的DH线性模型明显与GC得出的扩散层厚度不一致.无论考虑GC或DH模型, 它们求得的双电层厚度都比λ_D厚度偏大.
图 4(Fig. 4)

图 4 基于德拜解析式、GC和DH模型的双电层厚度求解Fig.4 EDL thickness solution by using Debye equation, GC and DH models

受限于激励电压大小的取值(ψ＜ψ_T), GC和DH模型都不能精确求解较大电压的带电离子分布.GC和DH模型主要以PB方程作为前提假设, 然而该方程理想性地把带电离子看作点电荷存在, 在现实物理空间中忽略了离子的几何特性.以Cl^-离子在金属阳极周围为例, 此时的Cl^-离子初始浓度为1 mol/L.在普通室温条件下, 当激励电压增大至10 V时, 采用式(1)得到的Cl^-离子的浓度应为1.34×10¹⁷² mol/L.一般NaCl溶液中的氯离子平均半径约为282.01 pm, 那么1 nm²面积内将会有数以万计的氯离子分布于阳极表面, 而此结果完全有悖于氯离子的真实几何形状.对于PB方程的不足, 需要对该表达式进行适当的物理修改, 在此重新引入Stern层.改进型PB方程如式(7)所示, 然而该数学公式具有非连续性，不符合物理求解需要.通过适当地修改后, 改进型PB方程如下所示:

(7)

对比式(1), 式(7)中引入了ν因子, 其等于(c_max为离子富集于电极表面的最大浓度, 以Cl^-为例, 其最大浓度为117.2 mol/L).由此GCS模型引入改进型PB方程的求解域表达式应满足式(8).采用式(1)和式(7), 图 5显示c_0∞初始为0.01 mol/L, 随着电压的不断增大也使阳极表面的负离子浓度不断增加.然而式(7)得出的结果显示表面Cl^-浓度维持在一个恒定值, 即便激励电压ψ_E怎么增大都不受影响.
图 5(Fig. 5)

图 5 PB与MPB方程的表面浓度对比图Fig.5 Comparison of surface concentration ratio between PB and MPB equations

(8)

3 电渗流仿真与分析如上所述, MPB方程能够更好地描述双电层内部的离子分布, 本文的电渗流建模首先通过式(8)求解得出双电层的厚度L_EDL, 即可得到对应的单位面积电容率C_EDL=ε₀ε_r/L_EDL.不同于纳米尺度下的电场求解, 电渗流一般在微米尺度下进行研究, 电解液保持电中性, 所以电势V^*(*代表复数)满足式(9):

(9)

式中:ω为激励电压角频率; σ为溶液电导率.金属电极/溶液界面交界处应满足式(10):

(10)

式中Δψ_EDL=ψ_E-V.对于微流控芯片而言, 其雷诺数Re相对较小(＜10^-3), 可以忽略惯性分量, 则Navior-Stokes方程可改写为

(11)

式中:p为溶液压强; u为流场速度矢量; η代表溶液黏滞系数(10^-3Pa·s).
在微米尺度下(远大于纳米尺度), 相应双电层间距近似为0, 那么双电层界面的电渗流速度满足式(12):

(12)

式中, Λ为双电层内部电容比值(, C_stern和C_diffuse分别为Stern和扩散层电容), 取值约等于0.25^[10].电场与流场的求解域和边界条件如图 6所示, 其余电场边界条件满足第二类纽曼边界条件(n·Δψ=0).流体周围边界为非滑移边界条件(u=0).AC电渗流的两侧激励金属电极宽度为500 μm, 间距25 μm, 激励电压频率为500 Hz.本文仿真的溶液电导率σ为2.1 mS/m, 根据科尔劳施无限稀释摩尔电导率可推得初始KCl的离子浓度c₀^∞应为1.4×10^-4 mol/L^[11].式(8)可得L_EDL为89.82 nm, 该值比线性Debye长度25.95 nm更宽.通过以上的偏微分方程求解得出芯片中的电渗流运动, 如图 7所示.
图 6(Fig. 6)

图 6 电场与流场求解域与边界条件示意图Fig.6 Solved domain and boundary conditions for electric and fluidic fields

图 7(Fig. 7)

图 7 德拜与MPB方程求解的电渗流运动Fig.7 Electro-osmotic flows by using Debye and MPB equations (a)—基于Debye公式模拟电渗流；(b)—基于MPB公式模拟电渗流.

4 结论基于MPB方程描述了电渗流内部双电层的正/负离子分布模型, 该模型不受溶液离子浓度和外界激励电压大小的影响, 可以更加准确地得到扩散层厚度.当稀释的KCl溶液电导率为2.1 mS/m, 且左、右激励电压为±1 V时, 求解得到的电渗流最大速度为1 034.31 μm/s, 液体对流范围较宽, 相比传统线性Debye求解更接近实验观测.
参考文献

[1]	Chen Y C, Li P, Huang P H, et al. Rare cell isolation and analysis in microfluidics[J].Lab on a Chip, 2014, 14(4): 626–645.DOI:10.1039/c3lc90136j
[2]	徐章润, 唐小燕, 王建华, 等. 重力驱动小型滴汞泵的研究[J].东北大学学报(自然科学版), 2008, 29(3): 449–452. ( Xu Zhang-run, Tang Xiao-yan, Wang Jian-hua, et al. Development of a gravity-driven mercury-dropping pump[J].Journal of Northeastern University(Natural Science), 2008, 29(3): 449–452.DOI:10.3321/j.issn:1005-3026.2008.03.037)
[3]	Dittrich P S, Manz A. Lab-on-a-chip:microfluidics in drug discovery[J].Nature Reviews Drug Discovery, 2006, 5(3): 210–218.DOI:10.1038/nrd1985
[4]	Zou Z W, Jang A, Macknight E, et al. Environmentally friendly disposable sensors with microfabricated on-chip planar bismuth electrode for in situ heavy metal ions measurement[J].Sensors and Actuators B:Chemical, 2008, 134(1): 18–24.DOI:10.1016/j.snb.2008.04.005
[5]	Kolluri N, Klapperich C M, Cabodi M. Towards lab-on-a-chip diagnostics for malaria elimination[J].Lab on a Chip, 2018, 18(1): 75–94.DOI:10.1039/C7LC00758B
[6]	Zhao B S, Koo Y M, Chung D S. Separations based on the mechanical forces of light[J].Analytica Chimica Acta, 2006, 556(1): 97–103.DOI:10.1016/j.aca.2005.06.065
[7]	Shields C W, Reyes C D, Lopez G P. Microfluidic cell sorting:a review of the advances in the separation of cells from debulking to rare cell isolation[J].Lab on a Chip, 2015, 15(5): 1230–1249.DOI:10.1039/C4LC01246A
[8]	Pething R. Dielectrophoresis theory, methodology and biological applications[M]. New York: Wiley, 2017: 1-448.
[9]	Hwang H, Park J K. Optoelectrofluidic platforms for chemistry and biology[J].Lab on a Chip, 2011, 11(1): 33–47.DOI:10.1039/C0LC00117A
[10]	Green N G, Ramos A, Gonzalez A. Fluid flow induced by nonuniform AC electric fields in electrolytes on microelectrodes.III.observation of streamlines and numerical simulation[J].Physical Review E, 2002, 66(3): 026305(1-11).
[11]	Pham P, Howorth M, Anne P C, et al.Numerical simulation of the electrical double layer based on the Poisson-Boltzmann models for AC electroosmosis flows[C]//Proceedings of the COMSOL Users Conference 2007.Grenoble: COMSOL, 2007: 1-10.