ABSTRACT

Electrokinetic microfluidics is one of the main subcategories for sample transport in Micro-TotalAnalysis-System. Numerical simulation is an efficient method for calculating and analyzing transport properties of samples in electrokinetic microfluidics. However, since the strong discontinuity exists at the interface of sample zone, the numerical dissipation near the interface becomes a key problem that can seriously influence the computational accuracy. In this paper, the reasons of numerical dissipation are theoretically analyzed and the techniques based on free interface reconstruction are provided to weaken the numerical dissipation. Firstly, the electric field and the flow field are solved by using the finite volume method. Then, in order to constraint the numerical dissipation, several techniques based on the volume of fluid (VOF) are utilized to reconstruct the flow vectors and solve the sample transport. Lastly, with the developed algorithm, capillary electrophoresis and electrokinetically driven enhanced mixing in microchannels are simulated. The results show that both FCT-VOF and TVD-VOF can significantly constraint the numerical dissipation while keeping the interface smoothing and are suitable for the simulation of electrokinetic microfluidics and other applications.

Keywords:Free Interface Reconstruction, Electrokinetic Microfluidics, Finite Volume Method, Volume of Fluid

用于电动微流体计算模拟的自由界面
重构技术比较

徐  征^*，森本龙，陆佳庆，胡帅龙，王林刚

辽宁省微纳米技术及系统重点实验室，大连理工大学，大连

Email: ^*xuzheng@dlut.edu.cn

收稿日期：2014年4月21日；修回日期：2014年5月20日；录用日期：2014年5月27日

摘  要

电动微流体是微全分析系统中样品输运的主要形式，数值模拟是计算和分析电动微流体样品输运特性的重要手段。由于微通道中区带样品在流体内存在大浓度梯度的界面，在计算时很容易引起数值耗散，对计算精度影响严重。本文探讨了电动微流体样品输运计算中数值耗散的产生原因和消除措施，利用有限体积法求解电场分布和流场分布，在此基础上应用自由界面重构技术对数值流向量进行了处理并求解微通道中样品传输，有效消除了界面处数值耗散的影响，并给出微通道中电泳分离和电致混合增强的计算模拟结果，研究结果表明：采用FCT-VOF和基于积分平均型TVD格式的VOF方法对上述算例的样品界面进行重构，能够保持界面平滑的同时，有效抑制了数值耗散，获得较好的计算精度。

关键词

自由界面重构，电动微流体，有限体积法，VOF

1. 引言

微全分析系统(micro total analysis system, μTAS)是利用微加工技术，将混合、反应、预浓缩、分离等功能集成在芯片上，将以往在多台仪器上才能完成的分析工作，转移到平方厘米大小的芯片上集中实现。具有样品耗量少、功能集成等突出优势，在生化分析、临床诊断、环境监测等领域具有广泛应用前景。电动微流体是通过施加电场对微通道内的流体进行驱动和操纵，从而达到控制样品定向输运的目的，目前采用外加电场驱动流体传输是μTAS的主要驱动方式之一。但电动微流体是复杂的耦合流动现象，电场和流场的紧密耦合增加了实验测量的难度。而采用数值计算方法则可以获得丰富的电动微流体流动信息，为μTAS研制提供理论依据，从而有效降低μTAS研发成本。生化样品在μTAS中混合、分离、定位都涉及到样品传输问题，运动中的样品形态和分布是评价μTAS的重要指标，而计算过程中的离散化会舍入低阶项而带来数值耗散，由于μTAS中实际样品区带存在较大的浓度梯度，会放大数值耗散效应，导致计算误差增大[1] [2] 。

本文针对上述问题，理论探讨了电动微流体物质输运计算中数值耗散的产生原因，构建了电动微流体的理论模型和有限体积法离散化框架，编制了相应的计算代码，利用有限体积法求解电场分布和流场分布，应用VOF自由界面重构技术对数值流向量进行处理并求解了样品传输，抑制了样品传输方程中数值耗散的影响，并对微通道中电泳分离和电致混合增强两种电动微流体现象进行了计算，验证了算法效果。

2. 数值耗散的原因分析

样品在微通道中传输取决于缓冲溶液的扩散、流体对流、样品电泳迁移作用，假设所研究样品为小c分子稀溶液，电场作用下的样品传输方程可用式(1)描述。

(1)

式中，为样品浓度，为时间，为方向坐标，为方向速度分量，为扩散系数，对上式任一方向的点对流项微分部分采用一阶迎风格式离散化，则有

(2)

不考虑二阶以上舍入余项的影响，则数值耗散主要是二阶项离散引起的。以微通道中的毛细管电泳为例，由于实际样品在微通道中传输时，样品扩散系数小，样品与缓冲溶液界面之间呈不连续或间断状态。如按真实初始状态设置计算初值，数值耗散使原本不连续的界面趋于平滑，出现样品区带的虚假增宽，如图1所示。对方程中对流项采用高阶格式离散化能减小数值耗散，但高阶格式都是条件稳定和非有界的，其应用范围受到一定限制。以QUICK高阶格式为例，仅当满足网格Peclet数条件时，一维对流-扩散方程才是稳定的，、、和分别为流体密度和速度、网格宽度和扩散系数。当数超过这一范围，容易出现负浓度等计算越界问题，如图2所示。

当与对流作用相比，样品扩散至可以忽略时，可以认为样品和缓冲溶液是检测周期内无相互扩散和渗透的两种流体，流体交界面称为自由界面。对类似流体体系，在界面附近的数值耗散最为突出，可以采用波前追踪方法、VOF方法、LevelSet方法等自由面计算方法，通过在样品界面附近的特殊处理，抑制数值耗散，保持自由面的锐利和准确，提高样品传输的计算精度[3] 。

3. 电动微流体输运理论模型与离散化方法

电动微流体理论由外加电场模型、双电层模型、流体流动模型和样品输运模型组成。为简化计算，这里忽略了流体流动对电场分布的影响和样品输运对流体粘度等的影响。

图1. 电泳计算中的虚假样品区带增宽现象

图2. 样品电泳迁移运动数值模拟中的负浓度现象

外加电场模型可用泊松方程描述，如式(3)；双电层模型可用Stern模型和泊松方程计算(式(4)和式(5))；电场作用下的流体流动用式(6)描述，其中包含电场力项，电场力可通过式7求得；带电荷物质在电场作用下的电泳速度通过式(8)求得；式(9)是电动流体中样品传输方程。首先通过求解式(3)、式(4)、式(5)求得电场分布，然后求解流场分布和电泳速度，并将二者合成，带入式(9)，求解样品分布[4] [5] 。

(3)

(4)

(5)

(6)

(7)

(8)

(9)

式中，为外加电场产生的电势，为电荷密度，为流体介电常数，为双电层zeta电势引起的电势，为方向双电层厚度，为电渗流速度矢量，为压强，为流体黏度，为电场力矢量，为电场强度，为样品电泳速度，为样品浓度，为样品扩散系数。考虑到微通道存在弯曲和突起等不规则形状，选择了贴体坐标系建立上述控制方程，并基于有限体积法离散化。

4. VOF自由面重构技术

对微通道中区带样品而言，在两相界面处存在大浓度梯度。可表示为：，；。忽略分子运动产生的弱自由扩散，样品输运方程可简化为

(10)

令，在贴体坐标系下上式守恒形式为

(11)

式中，、和为贴体坐标系坐标轴方向，、和为贴体坐标系速度分量。

采用有限体积法进行离散化的形式为

(12)

式中，为当前时间步，为通过控制体边界的各方向数值流向量，如图3，各方向数值流向量的积分形式为

，，

，，

图3. 计算平面内贴体网格控制体的界面与节点分布

对式12求解，关键在于如何避免或消除自由面附近的数值振荡和数值耗散。处理自由面的数值方法有多种多样，VOF法是主要方法之一，其特点是在固定网格内将自由面定义为流体体积函数，并构造流体函数的演化方法，从而确定自由面的位置和形状。VOF一般由自由面捕捉和自由面重构两部分组成，本文选择了Hirt & Nichols的施主–受主方法、FCT方法和基于积分平均型TVD格式的方法进行研究。

•     Hirt & Nichols的施主–受主方法

Hir和Nichols提出的施主–受主(Donor-Acceptor)格式的数值流向量构造方法是最早出现的VOF方法。该方法采用交错网格，体积函数和压力定义在控制体单元中心处，而方向速度定义在控制体左右边界中点，方向速度定义在上下边界的中点。贴体坐标系下，对空间网格坐标(i,j,k)点的式12积分可得

即：

、和分别对应、和方向的数值流向量差。基于施主-受主思想，可构造如下数值流向量形式，，下标和分别表示受主和施主控制体，下标表示可能取受主控制体值，也可能取施主控制体值，具体要根据界面速度矢量方向和自由面法向方向确定，取最小值是为了限制通过界面的值不超过施主控制体内流体体积(对应电泳中样品溶液)，取最大值是为约束通过界面的空体积流量(对应电泳中缓冲溶液)不超过施主控制体内空体积。和构造形式类似。实质上，这种构造方式是通过重构界面流向量替代了原来通过线性插值计算界面流向量的方法，从而消除了对流项离散化引入的数值耗散。

•     通量校正传输(Flux Corrected Transport, FCT)VOF方法

FCT-VOF方法是基于通量修正输运算法的思想来求解流体体积函数。低阶格式的优点是能保持单调性和无振荡，但耗散性强。而高阶格式分辨率高，但容易引起色散振荡。通量修正输运的思路是先利用低阶通量计算中间值，再利用高阶通量值进行修正。在构造数值流向量时，先取迎风格式构造低阶数值流向量，再采用逆风格式构造高阶数值流向量，以e方向边界数值流向量构造为例：

,  

最终求得的数值流向量为，FCT-VOF的具体实施方法有多种形式，本文主要采用Zalesak[6] 和Gopala[7] 等建立的数值流向量构造形式。

•     基于积分平均型TVD格式的VOF方法

基于积分平均型的VOF方法是源于变差不增(Total variation diminishing, TVD)的概念。对于一维传输

方程，其一般离散格式形式为。总变差定义为，总变差

不增即：。TVD性质能保证格式单调，能抑制跳跃和间断处的振荡，并保持数值解的锐利形态。而对于不具备TVD性质的格式，可结合Superbee等函数修正，使其具有TVD性质，TVD可用库朗数(Courant Number)作为判据：。利用TVD构造数值流的方法，可以构造TVD格式的高精度数值流向量，进行自由面重构。对一维传输方程，有：。

将方程转换为分片形式，利用特征线法推进时间步求解上述方程，并带入控制体进行积分，可得

当速度，

当速度，

取，就获得Lax-Wendroff格式，该格式容易振荡，可用限制器对斜率值加以约束。这种方法不需要对自由面形态进行重构过程，而是结合时空平均的TVD格式和限制器，来调整数值流向量的数值，实施相对简单。

5. 计算与讨论

本文对微通道中样品电泳迁移和电致混合增强效应进行了模拟，以验证应用VOF方法计算电动微流体的效果。

5.1. 微通道中样品电泳迁移的计算与讨论

微通道取为，网格数取为，速度场，扩散系数设为零，时间步，样品区带初值取为立方体(图4(a))，图中颜色表示样品的浓度。由计算结果可知：Hirt-VOF方法(图4(b))、FCT-VOF方法(图4(c))、基于积分平均型TVD方法(图4(d))都可以给出区带电泳运动趋势和位置，FCT-VOF方法的效果更好一些，界面清晰，无明显虚假区带增宽；而积分平均TVD方法在界面处有所光滑，区带宽度也略有增加；采用Hirt-VOF方法效果最差(区带增宽约8个网格)，这是因为Hirt-VOF方法仅在低阶的迎风格式和逆风格式间进行切换，数值耗散严重，而基于积分平均型TVD-VOF和FCT-VOF方法均采用了高阶抑制数值耗散的措施(表1)。

5.2. 电致混合增强效应的计算与讨论

生物样品和试剂的混合是μTAS重要环节，混合效率是微混合器的关键，许多生物操作都需要快速和均匀的样品混合。而利用电场良好可控性，能够实现了电致对流混合增强的效果，控制壁面zeta电势是其中一种实施方法[8] [9] 。通过设计壁面zeta电势分布，可以使流体受力不均而诱导流线折叠和拉伸，增加实际混合面积。本文设计了等直微通道(：，网格数取为)，对zeta电势均匀分布、对称间隔分布和交错间隔分布的对混合影响进行数值模拟。均匀分布是指壁面zeta电势处处等于−0.1 V；对于对称间隔分布，侧壁zeta电势均以3为间隔长度，分别取为−0.1 V和+0.1 V，且左右侧壁分布对称；对交错间隔分布，侧壁zeta电势也以3为间隔长度，分别取−0.1 V和+0.1 V，但左、右侧壁zeta电势分布呈相差一个间隔长度的交错形式。其余壁面处zeta电势均取−0.1 V。采用物性参数见表2。

计算得到的电场和流场分布见图5，图中颜色及等值线图为zeta电势分布，图中矢量箭头为流速，

图4. 利用VOF方法求解样品电泳迁移问题。(a) 初始样品区带；(b) Hirt-VOF方法；(c) FCT-VOF方法；(d) 基于积分平均型TVD格式的VOF方法

图5. 微混合通道内的电渗流速度形态。(a) 壁面zeta电势一致分布；(b) 壁面zeta电势对称间隔分布；(c) 壁面zeta电势交错间隔分布

图6. 不同zeta电势分布对混合效果的影响。(a) zeta电势一致分布；(b) zeta电势对称间隔分布；(c) zeta电势交错间隔分布

可见：对zeta电势一致分布的情况(图5(a))，流场形态为塞状流，无侧向拉伸；而对zeta电势对称间隔分布情况(图5(b))，流场在电场力方向交替作用下，呈周期性聚合与扩展变化，但扭曲和折叠不明显；对zeta电势交错间隔分布的情况(图5(c))，不对称电场力对流体有剪切作用，流动方向呈明显周期性扭曲变化。

样品在某截面的均匀分布程度体现了混合效果。设在入口处中间截面的样品浓度设为1.0，入口的其它节点设为0.0，观察内部某一截面的样品侧向分散，如图6所示(图中颜色为样品浓度)：壁面zeta电势分布一致的情况下，样品在横截面方向分散是非常有限的，混合效率很低(图6(a))；壁面zeta电势对称间隔分布可以促进样品侧向分散，在距入口34处，达到完全均匀分布(图6(b))；壁面zeta电势交错间隔分布侧向分散增强作用最为明显，在距入口20处，就达到完全均匀分布(图6(c))。可见，采用交错间隔的zeta电势分布能有效增强侧向分散，有利于提高混合效率。

6. 结论

本文在用有限体积法计算电动微流体电场和流场分布基础上，研究了利用Hirt-VOF方法、FCT-VOF方法和基于积分平均型TVD格式VOF方法求解电动流体中样品输运问题，并给出对微通道中常规毛细管电泳分离和电致对流混合增强的计算结果。研究结果表明：采用FCT-VOF方法和基于积分平均型TVD格式的方法能抑制了迎风格式的耗散效应，在样品区带界面附近能抑制振荡和保持数值解的锐利形态，并易于与电动微流体其他计算模块组合使用。

致  谢

本文研究工作得到中央高校基本科研业务费(No. DUT14LAB07)资金资助，特此感谢！

参考文献 (References)

NOTES

^*通讯作者。