分支型微通道水油-两相流数值模拟＊

王琳琳，胡洪萍，田 辉

（1．西安文理学院 数学与计算机工程学院，西安710065；2．西安交通大学 能源与动力工程学院，西安710049；3．承德石油高等专科学校 机械工程系，承德067000）

自 微 机 电 系 统 （Micro Electro Mechanical Systems，MEMS）概念提出以来，微流体混合器内气泡和液滴在连续相中形成已成为微流动动力学研究中的重要问题之一．和常规尺度通道相比，微流体混合器内两相流体流动具有大的面积／体积比，有利于提高流体间的传质、换热和反应率，安全性能优越［1-2］，可控性良好［3］，目前被应用于纳米粒子合成［4］、DNA 分析［5］、微流体乳化、神经造影剂和分子筛检［6］等领域．

常见的微流体混合器分为：同向流动微通道、汇聚型微通道、T型微通道［7］．同向流动微通道中离散相和连续相流动方向相同；汇聚型微通道的连续相通道和主通道同轴，离散相由主通道两侧进入；T型微通道有两种类型：连续相通道和离散相通道垂直，与主通道同轴，称为错流接触T型微通道［8］；连续相和离散相入口方向相反，两通道和主通道垂直，称为对流接触T型微通道［9］．毛细数Ca表示黏性力和表面张力的比值，当毛细数小于10-2，微流体混合器内离散相的形成受到通道壁面约束，错流接触T型微通道内离散相大小受到通道宽度和高度、入口流量比、毛细数、壁面条件等影响，但与两相黏度比无关［9-15］．而人们对于其它结构的微流体混合器研究成果较少．

微流体混合器通道宽度微小，其特征尺度通常小于1mm，实验中对微流体混合器两相流体进口流量的控制，观察气泡和液滴脱离和形成过程，测量脱离气泡和液滴长度、体积，及通道的压强变化和流场分布等，这些都对设备和器材的精度要求极高，微小条件的变化都可能对实验结果产生明显的变化，研究成本也较高．随着计算机技术的飞速发展，数值模拟技术已成为研究流体动力学问题的重要手段之一．常见模拟两相流体界面运动的数值方法有标记粒子单元方法（Marker and Cell，MAC）、流体体积函数方法（Volume of Fluid，VOF）、相场方法（Phase Field Method，PFM）和水平集方法（Level Set Method，LSM）等．水平集方法［16］通过引入水平集函数，用该函数是某一固定值对应的集合来表示移动的两相界面，这种方法具有对界面捕获精度高、不需重构界面、对界面复杂拓扑变化描述准确等特点，被广泛应用于图像处理、界面演化、流体力学、燃烧和材料力学等领域．

文中通过水平集方法，建立适合的控制方程，对分支型结构微通道内水-油两相流进行数值研究，考察入口角、通道宽度和滑移系数对离散相长度的影响，得到液滴长度随这些参数的变化情况，为精确控制微通道内液滴尺寸和微流体系统优化设计提供一定参考．

1 几何模型和数值方程

1．1 几何模型和物理性质

文中以分支型微通道作为研究的几何模型，如图1所示．离散相和连续相通道之间有一定夹角，称为入口角［17］，用α表示．分支型微通道的连续相通道和主通道水平，水和油分别从左侧离散相通道和连续相通道入口进入，在微通道内混合，最后由通道右侧出口处流出．离散相和连续相通道长度为3D1，其中D1＝111μm，主通道长度是28D1，离散相通道宽度D和主通道宽度W可取0．5D1、D1、1．5D1，连续相通道宽度和主通道宽度相同．设置水和油的密度分别是998kg·m-3和750kg·m-3、黏性系数分别是0．001Pas和0．001 34Pas、表面张力系数是0．072 8N·m-1．

图1 分支型微通道的几何模型Fig．1 The geometric model of the bifurcation micro-channel

1．2 控制方程

由于微通道内水和油互溶性差，压降小，可视为不可压缩流体；两相流体入口速度较小，使得雷诺数小［18］，设定通道内是层流流动；邦德数小于1，故模拟中可忽略重力影响［10，18-20］．因此，对微通道内两相流动进行模拟的连续性方程和动量方程为

式中：Fσ为表面张力；p为压强；U为速度向量，两相流体的密度和黏性系数由以下等式得出

下标w和o分别表示水和油，φ是水平集函数，文中使用守恒的水平集函数［21］，水相对应的水平集函数值为0，油相对应的函数值为1，两相界面的水平集函数值为0．5．水平集函数由以下方程计算

式中：γ为重新初始化参数；ε为界面厚度．式（2）中的表面张力为

1．3 初始条件和边界条件

初始时刻假定离散相通道内充满水，其它通道被油充满，两流体处于静止状态．

水、油和壁面三相接触，会出现浸润现象，在表面张力作用下，油相和壁面间形成接触角θ，接触角的大小反映壁面对油相吸附力的强弱，接触角和n之间满足

其中nw、tw分别为壁面单位法向量和单位切向量．

边界滑移现象对微／纳米尺度下的流动特性影响显著［22］．由于微通道内的流动主要受到压差驱动，壁面剪切率低，适应采用Navier滑移模型，即无滑移壁面距离实际壁面Ls处，流体在实际壁面上出现滑移行为，滑移速度为

其中：Ls为滑移长度；u为流体水平速度．定义壁面产生的摩擦力为

文中研究油相和通道壁面是完全浸润的，即壁面接触角是0°．水和油垂直于入口边界进入微通道，入口速度分别是0．027 3m·s-1，0．054 6m·s-1，出口压强为0Pa，壁面采用滑移边界，文中滑移长度取1～20μm．

文中主要研究微通道内两相流动在二维平面内的变化，因此在数值模拟中忽略沿微通道高度方向的流动．计算采用三角形网格进行计算．通过比较粗网格和细网格的计算结果，发现当网格数大于11 000，模拟结果和网格数无关，因此设置网格数是32 472个，取时间步长0．000 1s．

2 模拟结果与分析

2．1 液滴在通道内的形成

对主通道宽度和离散相通道宽度同是D1、入口角是60°、滑移长度为1μm的分支型微通道内两相流动进行数值模拟，得到通道内液滴的形成过程，如图2所示．

随着水的进口量不断增加，水前端逐渐进入通道混合区后，受到壁面约束向主通道下游移动，在通道拐角处液滴脱离，脱离后，除少量水弹回离散相通道出口外，主通道内液滴前后端收缩，形成形状稳定、长度大于通道宽度的Taylor液滴．

图2 液滴在分支型微通道内的形成过程Fig．2 The droplet formation in the bifurcation micro-channel

T型错流接触微通道Taylor液滴形成中通道内压强会产生明显波动［10］，液滴后端压强增大表明液滴受到破坏力增大，进而促进液滴的脱离．为研究分支型微通道内压强在液滴形成过程中的变化情况，取一段水平通道的轴线AB（如图1所示），在图2中四个时刻AB上的压强分布如图3所示．

图3 液滴形成过程中微通道内压强的变化Fig．3 The change of the pressure in the micro-channel during the droplet formation

由图3可知液滴内的压强大于其两侧油相的压强，随着水进入主通道，液滴压强由1 500Pa左右升高到1 680Pa（a～b）后降至约1 572Pa（b～d）．油相在流动中受到液滴阻挡使液滴上游压强迅速由600Pa升高至1 420Pa（a～b），压强增大的油相开始推动液滴后端向通道拐角处移动，消耗油相的压能，油相压强缓慢增至1 567Pa（b～c），造成液滴后端内外压降几乎相同，液滴脱离前，液滴前后端总压降达到最大值．脱离后（d），液滴内部压强保持脱离前大小，在表面张力作用下，液滴两端收缩呈圆形，其后端油相压强大幅下降至300Pa．

2．2 入口角对液滴长度的影响

入口角是两相流体入口速度方向的夹角，来流速度方向影响微通道内两相流动．为研究入口角对液滴尺寸的影响，取入口角为10°～170°，壁面滑移长度和通道宽度同文中2．1部分，计算结果如图4所示．相对于D1液滴长度变化范围约是2．4～5．5，入口角从10°增至100°时，液滴长度缩短；入口角由100°增加到170°时，液滴长度随之增加．对计算数据拟合出对应曲线，结果表明微通道内液滴长度随入口角的增大呈二次变化．

图4 入口角对液滴长度的影响Fig．4 The influence of the inlet angles on the droplet length

由图3知受到液滴阻碍，液滴后端的油相压强逐渐增大，由于油相的挤压作用对Taylor液滴形成至关重要，接下来研究在不同入口角下液滴后端的油相压强在液滴形成过程中的变化．由前可知，液滴的形成过程分为三个阶段：水的前端逐渐进入微通道，增大的液滴堵塞主通道，其后端油相压强迅速增加，这一阶段称为阻塞阶段；积压驱动液滴向主通道下游移动，液滴颈部形成并逐渐变细，油相压强增长程度减缓，这一阶段称为挤压阶段；液滴从颈部“断裂”后，油相压强迅速减小，通道内形成形状稳定的液滴，这一阶段称为脱离阶段．液滴大小受阻塞阶段和挤压阶段持续时间的影响．

连续相通道出口取点A1如图1用以监测液滴后端油相的压强，取入口角分别为30°、90°、150°，得到液滴形成中油相压强的变化，如图5所示．图中可见，在液滴形成的各阶段，30°入口角的通道内，油相压强增长相对缓慢，油相对液滴的破坏力增长平缓造成液滴脱离点向通道下游移动，延长液滴脱离时间如图5（b）；入口角增至90°，油相压强增长速度加快，造成油相的挤压力增加，使得液滴在通道拐角处“断裂”，加快液滴脱离如图5（c）；入口角增大到150°，阻塞阶段同时刻压强稍有增大，挤压阶段压强增长速度有所提高，但各阶段内油相压强的极大值减小，其原因是入口角增大，微通道拐角形状逐渐呈尖锐的楔形，在后端油相挤压下，更多油相挤入离散相通道内，液滴在油相压强还没有增至充分大就已在拐角处被“夹断”如图5（d），说明楔形通道拐角促进液滴的脱离．

图5 不同入口角下液滴形成中连续相通道压强和液滴脱离位置的变化Fig．5 The influence of the inlet angles on the pressure change in the continuous phase channel and droplet departure positions

2．3 通道尺寸对液滴长度的影响

2．3．1 离散相通道宽度对液滴长度的影响

由于微通道尺寸影响两相流动，为考察离散相通道宽度对液滴长度的影响，在主通道宽度W＝D1的条件下，设置离散相通道宽度分别为0．5D1、D1、1．5D1，得到液滴长度的变化，如图6所示．

由图6可知，同一入口角下，离散相通道越宽，液滴越长；不同入口角下，液滴长度变化幅度随离散相通道宽度变窄而减小；离散相通道宽度减小至0．5D1并且入口角大于100°，液滴长度几乎不随入口角变化．对计算数据进行拟合，得到三种离散相宽度的通道内，随入口角增加液滴长度呈二次变化；当离散相通道宽度大于0．5D1，长度最短的液滴在入口角取100°左右的通道内取得；离散相通道宽度越窄，计算结果和拟合曲线符合程度越好．

图6 离散相通道宽度对液滴长度的影响Fig．6 The influence of the width of the dispersed phase channel on the droplet’s length

2．3．2 主通道宽度对液滴长度的影响

主通道宽度变化也会使液滴长度发生变化，为考察主通道宽度对液滴长度的影响，固定离散相通道宽度D＝D1，设置主通道宽度分别为0．5D1、D1、1．5D1，模拟得到不同入口角下液滴长度的变化，如图7所示．

图7 主通道宽度对液滴长度的影响Fig．7 The influence of the width of the main channel on the droplet’s length

由图7中可看出，主通道越宽，入口角对液滴长度的影响越弱；入口角不大于150°，主通道宽度增加液滴长度随之增加；入口角越小，主通道宽度对液滴长度影响越大；随着入口角的增加，主通道宽度对液滴长度的影响逐渐减弱．将模拟数据进行拟合，得到在三种主通道宽度的通道内，液滴长度都随入口角呈二次变化；入口角大于150°，三条拟合曲线几乎重合；主通道宽度越宽，计算结果和拟合曲线符合程度越好．

2．4 滑移长度对离散相长度的影响

边界滑移是指流体和壁面之间存在相对滑动，由于微通道的特征尺度通常小于1mm，使得壁面对微流动有显著影响．Chang等［23］试验测得水的滑移系数为2～20μm，甘油的滑移系数为5～50 μm；刘赵淼等［22］研究了以油为介质、滑移系数为0．5～10μm的微流动．为了进一步研究壁面滑移作用对微通道内流动特性的影响，固定离散相通道壁面滑移长度为1μm，取主通道的滑移长度为1～20μm．入口角为60°的微通道内，随着主通道滑移长度的增加，液滴长度变化如图8所示．滑移长度大的主通道内形成长的液滴；随着壁面滑移长度的增加，液滴长度增长程度有所减弱．

图8 主通道壁面的滑移长度对液滴长度的影响Fig．8 The influence of the slip length of the main channel on the droplet’s length

油相挤压力对液滴表面起到破坏作用．为了考察油相挤压作用受通道滑移参数的影响，取滑移长度分别是4、8、12、16、20μm，液滴形成中其后端油相压强的变化如图9所示．

由图9可知，滑移长度增加造成了油相压强减小；随着滑移长度增加，压强减小程度变弱．这些结果说明，滑移长度增加，液滴后端受到油相的挤压作用减弱，破坏力的减小延长了液滴脱离时间，使得液滴长度随壁面滑移程度的增加而增加；但随着滑移长度增加，油相挤压力减小程度变弱，造成液滴增长减弱，解释了图8中模拟结果形成的原因．

滑移长度影响通道内流体的速度分布．微通道内的流动是泊肃叶流动，主通道中单相流体的控制方程为

图9 不同滑移长度对油相压强的影响Fig．9 The influence of different slip length the oil pressure

对式（10）积分并结合壁面滑移速度是条件（8），得到微通道内压强梯度作用引起的单相流体速度分布为

图10中给出不同滑移参数下，主通道中油相水平速度的的模拟值和由式（11）得到的理论值．

图10 不同滑移参数下，油相水平速度的模拟值和理论值Fig．10 The simulation and theoretical results of the velocity of oil in the continuous phase channel with different slip lengths

由图10中可看出，滑移参数取4μm和20 μm，水平速度的模拟值和理论值几乎吻合；油相水平速度沿主通道轴线对称并且呈抛物线分布，随着滑移参数由4μm增加到20μm，水平速度最大值逐渐由0．116m／s减小至0．106m／s，壁面滑移速度的增加使通道内速度分布更加均匀，距离主通道上下壁面约24．25μm处，各速度曲线交于两个公共点，由此得到油相受到壁面滑移作用的厚度约为62．5μm，占主通道的56．3%．

3 结 论

文中借助水平集方法，采用Navier滑移模型，在油相和壁面完全浸润和条件下，对分支型微通道液滴在油内的形成过程进行模拟，得到液滴形成过程中两相流体压强的变化特点．针对通道几何形状、尺寸和壁面条件对两相流动的影响，通过改变入口角、离散相和主通道宽度、主通道滑移参数，得到液滴尺寸随这些参数的变化：入口角增大造成离散相长度呈二次变化，并促进液滴脱离；入口角小于160°，离散相或连续相通道宽度增加都能增加液滴长度；离散相通道变窄或主通道变宽能够提高模拟结果和拟合曲线的符合程度；固定入口角是60°，在两种滑移尺寸下，主通道内油相水平速度分布的模拟值和理论值都几乎一致，主通道滑移尺寸从1μm增至20μm，液滴长度随之增加，但增长程度逐渐渐弱，并发现主通道内流体受壁面作用的厚度约62．5μm．上述研究结果对进一步研究微通道几何结构、尺寸和壁面条件对液滴的影响、为有效控制微通道内液滴尺寸、优化设计微流体混合系统提供一定的参考．

［1］ 郑亚锋，赵阳，辛峰．微反应器研究及展望［J］．化工进展，2002，23（5）：461．ZHENG Ya-feng，ZHAO Yang，XIN Feng．Research on and Prospects of Microreactors［J］．Chemical Industry and Engineering Progress，2004，23（5）：461．（in Chinese）

［2］ 骆广生，王凯，吕阳成，等．微反应器研究最新进展［J］．现代化工，2009，29（5）：27．LUO Guang-sheng，WANG Kai，LV Yang-chen，et al．Advances in Research on Micro-reactors［J］．Modern Chemical Industry，2009，29（5）：27．（in Chinese）

［3］ WANG W，XIE R，JU X J，et al．Controllable Microfluidic Production of Multicomponent Multiple Emulsions［J］．Lab on a Chip2011，11（9）：1587．

［4］ GUNTHER A，KHANS A，THALMANN M．et al．Transport and Reaction in Microscale Segmented Gas-Liquid Flow［J］．Lab on a Chip，2004，4（4）：278．

［5］ BURNS A，JOHNSON N，BRAHMASANDRA N，et al．An Integrated Nanoliter DNA Analysis Device［J］．Science，1998，282（5388）：484．

［6］ BERNARD T．KellyJean-Christophe Baret Valerie Taly-Andrew D．Griffiths．Miniaturizing Chemistry and Biology in Microdroplets［J］．Chem．Commun，2007，14（18）：1773．

［7］ ZHAO C X，ANTON P M．Two-Phase Microfluidic Flows［J］．Chemical Engineering Science，2011（66）：1394．

［8］ WANG L D，WANG F，XIN C F．Numerical Study on Bubble Formation of a Gas-Liquid Flow in a T-Junction Microchannel［J］．Chemical Engineering ＆Technology，2009，32（12）：1984．

［9］ VAN S，DER GRAAF T，NISISAKO C G，et al．Lattice Boltzmann Simulations of Droplet Formation in a T-Shaped Microchannel［J］．Langmuir，2006，22（9）：4144．

［10］ MENECH M D，MENCH P，GARSTECKI F，et al．Transition from Squeezing to Dripping in a Microfluidic T-Shaped Junction［J］．Journal of Fluid Mechanics，2008（595）：141．

［11］ GARSTECKI P，FUERSTMAN M J，STONE H A，et al．Formation of Droplets and Bubbles in a Microfluidic T-Junction-Scaling and Mechanism of Break-up［J］．Lab on a Chip，2006，6（3）：437．

［12］ XU J H ，LI S W，TAN J，et al．Correlations of Droplet Formation in T-Junction Microfluidic Devices：From Squeezing to Dripping［J］．Microfluidics and Nanofluidics，2008（5）：711．

［13］ STEEGMANS M L J，SCHROEN C G P H，BOOM R M．Generalised Insights in Droplet Formation at T-Junctions Through Statistical Analysis［J］．Chemical Engineering Science，2009，64（13）：3042．

［14］ VAN V，STEIJIN C R，KLEIJIN M，et al．Flow A-round Confined Bubbles and Their Importance in Triggering Pinch-off［J］．Phys．Rev．Lett，2009，103（214501）：1．

［15］ VAN V，STEIJIN C R，KLEIJIN M，et al．Predictive Model for the Size of Bubbles and Droplets Created in Microfluidic T-junctions［J］．Lab Chip，2010，10：2513．

［16］ OSHER S，SETHIAN J A．Fronts Propagating with Curvature-Dependent Speed：Algorithms Based on Hamilton-Jacobi Formulations［J］．Journal of Computational Physics，1988，79：12．

［17］ 朱春英，付涛涛，高兴群，等．微通道内液液两相流流型化工进展［J］，2011（30）：65．ZHU Chun-ying，FU Tao-tao，GAO Xing-qun，et al．Study on the Liquid-Liquid Two-Phase Flow Patterns in Micro-Channels［J］．Chemical Industry and Engineering Progress，2011，30：65．（in Chinese）

［18］ 袁希钢，宋文琦．T型结构微通道气液两相流型的数值模拟［J］．天津大学学报，201245（9）：763．YUAN Xi-gang，SONG Wen-qi．Numerical Simulation of Gas-liquid Two Phase Flow Pattern in T-junction Microchannel［J］．Journal of Tianjin University，2012，45（9）：763．（in Chinese）

［19］ EDVINSSON R，KLRANDOUST S．Finite-Element Analysis of Taylor Flow［J］．AIChE Journal，1996，42（7）：1815．

［20］ HAZEL A L，HEIL M．The Steady Propagation of a Semiinfinite Bubble into a Tube of Elliptical or Rectangular Cross-Section［J］．Journal of Fluid Mechanics，2002，470（1）：91．

［21］ OLSSON EKREISS G．A Conservative Level Set Method for Two Phase flow［J］．Journal of Computational Physics，2005（210）：225．

［22］ 刘赵淼，王国斌，申峰．基于Navier滑移的油膜缝隙微流动特性数值分析［J］．机械工程学报，2011，47（21）：104．LIU Zhao-miao，WANG Guo-bin，SHEN Feng．Numerical Analysis of Oil Film Flow in Micro Gap with NAvier Slip Boundary Conditions［J］．Journal of Mechanical Engineering，2011，47（21）：104．（in Chinese）

［23］ CHANG H C，CHANG J K．Large Sip of Aqueous Liquid Flow Over a Nanoengineered Superhydrophobic Surface［J］．Physical Review Letters，2006，96（6）：1．