壁面滑移条件下微尺度通道内两相流数值模拟

王琳琳， 李　泽

(1.西安文理学院 信息工程学院，西安 710065； 2.西安交通大学 能源与动力工程学院

热流科学与工程教育部重点实验室，西安 710049； 3.西安工程大学 理学院，西安 710048)



壁面滑移条件下微尺度通道内两相流数值模拟

王琳琳1,2， 李泽3

(1.西安文理学院 信息工程学院，西安 710065； 2.西安交通大学 能源与动力工程学院

热流科学与工程教育部重点实验室，西安 710049； 3.西安工程大学 理学院，西安 710048)

为研究Y型微通道内壁面滑移和连续相表观黏度对通道内液滴尺寸和流动特性的影响，采用水平集法进行数值模拟.结果表明：随着连续相通道滑移长度的增加，液滴长度增长，通道内流动更均匀.壁面滑移作用的范围占通道的46.8%.决定表观黏性系数的两参数增加，造成负滑移边界形成，促进液滴脱离，造成液滴长度缩短.表观黏度作用范围占通道的一半以上.对微流体混合系统的优化设计提供了参考.

微通道；数值模拟；水平集法；液滴长度；边界滑移；表观黏度

微型化是当今科学技术发展的一个重要方向.目前，气泡和液滴在微流体混合器内的形成已成为微流动动力学研究中的重要问题之一.与常规尺度通道相比，微流体混合器更利于提高流体间的传质、换热和反应率，具有良好的安全性和可控性[1]，已成功应用于微流体乳化[2]、化学和材料合成、分子进化的分割、药物输送、药品检测与开发、生物鉴定、食品工业和分子筛检[3]等众多领域.

常见的微流体混合器分为三类，分别是：T型微通道、汇聚型微通道、同向流动微通道[4].毛细数表示黏性力和表面张力的比值，毛细数小于0.01，微流体混合器的通道壁面约束气泡或液滴的形成，连续相的挤压力是造成其脱离的主要作用力，常见微通道内气泡或液滴的尺寸主要受通道宽度和高度、入口流量比、毛细数、壁面条件等因素影响[5-6].目前，对于其他结构微通道的研究成果相对较少.宏观尺度下，通常忽略壁面滑移对流动的影响.微通道内，流体在壁面的滑移行为对流动特性产生重要影响，然而对滑移行为的研究多集中在特定的物理模型，很少将研究滑移问题延伸到微通道内两相流中.

由于微通道尺度微小，测量微通道内气泡和液滴的形成过程中通道内压降和流场的变化、气泡和液滴长度，体积等数据对仪器的精度要求极高，外界极小地扰动都可能使实验结果发生显著变化，测量设备价格高昂.随着计算机运行速度和数值方法精度的提高，数值模拟已成为研究微流动问题的重要手段之一.模拟两相流界面运动数值方法有：标记粒子单元法、流体体积法、耗散粒子动力学法、相场法和水平集法等.水平集法是通过引入水平集函数，两相界面用该函数取固定值对应的点集来表示的一种方法，该方法具有捕获界面精度高、无需重新构造界面、能够准确描述拓扑结构发生复杂变化的界面等特点，广泛应用于图像处理、界面演化、燃烧、流体力学等领域.

文中基于水平集法，对Y型结构微通道内水-油两相流动进行数值模拟，分析壁面滑移长度、表观黏度对液滴长度和通道内流动的影响，为控制微通道内液滴尺寸和微流体设备的优化设计提供参考.

1　物理模型和控制方程

1.1几何模型和物性设定

图1　Y型微通道几何模型

文中以Y型微通道为研究对象，如图1所示.离散相和连续相通道之间的夹角，称为微通道入口角，用α表示，文中几何模型的微通道入口角是120°.Y型微通道的离散相通道和连续相通道关于主通道对称，水、油两相从左侧分别由离散相通道、连续相通道入口进入，在微通道内混合后，从右侧出口流出.由此看出，对流接触T型微通道是一种特殊的Y型微通道.各通道的宽度都是D1，D1=111μm，离散相和连续相通道长度均为3D1，主通道长度是35D1.水、油两相密度分别是998kg·m-3、750kg·m-3，黏性系数分别是0.001Pa·s、0.001 5Pa·s，界面张力系数是0.047 4N·m-1.

1.2控制方程

微通道内水和油互溶性差，压降小，流体视为不可压缩的；两相流体的入口速度较小造成雷诺数小[7]，通道内为层流流动；邦德数小于1，模拟中忽略重力的影响[7].因此，微通道内两相流动的连续性方程和动量方程分别是

·U=0

(1)

(2)

式中：U是速度向量，p是压强，Fσ是界面张力，两相流体密度和黏性系数

ρ=ρw+(ρo-ρw)φ

(3)

μ=μw+(μo-μw)φ

(4)

下标w和o分别代表水和油，两相流体对应的水平集函数φ不同，文中水平集函数采用守恒型，水平集函数为0和1的区域分别表示水和油，两相界面由水平集函数值为0.5的点组成.水平集函数满足方程

(5)

式中：γ是重新初始化参数，ε是界面厚度.式(2)中的界面张力

Fσ=σκδn

(6)

1.3初始条件和边界条件

初始时刻设定离散相通道充满水，其他通道充满油，通道内的流体静止.

微/纳米尺度下滑移边界对流动的影响显著.微通道内的流体主要受压差的驱动，壁面剪切率低，适合采用Navier滑移模型[8]，即距实际壁面外侧Ls处为无滑移壁面，Ls称为滑移长度.

水、油、壁面接触时，分子间引力的作用使相界面上出现浸润现象，油相和壁面间形成接触角θ，接触角反映油相对壁面的吸附能力.接触角和向量n满足关系

n=nwcosθ+twsinθ

(7)

式中：nw、tw分别是壁面单位法向量和切向量.文中接触角取0°，即油相和通道壁面完全浸润.

水油两相流体分别以ud=0.017 8m·s-1、uc=0.035 6m·s-1的速度垂直于通道入口流入，通道出口压强是0Pa，设置通道壁面是滑移边界.

文中主要研究Y型微通道内两相流动在二维平面的变化特性.采用三角形网格进行计算.通过对比不同粗细网格下的计算结果，发现网格数大于11 000，模拟结果与网格数无关，设置文中模型网格数是32 472，时间步长取0.000 02s.

2　模拟结果与分析

2.1滑移长度对液滴长度和连续相流动的影响

速度滑移是指流体在固体表面运动时，流体和壁面之间在界面处的切向速度差，此时壁面为滑移边界，滑移边界有助于减小壁面对流体的拖拽作用.微通道宽度通常小于1mm，壁面对微流动的影响尤其显著.Choi等通过流变系统测得水的滑移长度为2～20μm，体积分数30%甘油的滑移长度为5～50μm.刘赵淼等数值模拟研究了油膜厚度为50μm、滑移长度为0.5～10μm的微流动.为进一步研究Y型微通道内滑移边界对液滴长度和流动的影响规律，取离散相通道壁面滑移长度为1μm，连续相和主通道的滑移长度在1～20μm之间变化.随着滑移长度的增加，液滴长度L逐渐增加，但滑移长度由14μm继续增加，液滴长度增长减缓，如图2所示.

液滴后端油相的挤压力对液滴表面起破坏作用，挤压力增大促进液滴脱离.为研究滑移长度对油相挤压作用的影响，在连续相通道内取点A1，见图1，取滑移长度分别为4、8、12、16和20μm，液滴形成过程中点A1压强的变化如图3所示.图中可见，滑移长度增加，油相压强减小；滑移长度由14μm继续增大，压强减小程度变弱.这些结果表明，滑移长度增加造成液滴后端受到油相的挤压作用减弱，破坏力的减小使液滴脱离时间延长，从而造成液滴长度随壁面滑移程度的增加而增加；然而，随着滑移长度增加，油相对液滴后端挤压力减小程度变弱，造成液滴长度的增长程度减弱，由此解释了图2中模拟结果的原因.

图2　液滴长度随壁面滑移长度的变化

图3　连续相通道内压强随滑移长度的变化

为研究滑移长度对通道内速度分布的影响，在连续相通道内建立坐标系，如图4(a)，x1轴正方向为油相入口方向，y1轴垂直于该方向.连续相通道内油相的流动是泊肃叶流动，沿x1轴正方向的流速uo满足

(8)

结合滑移壁面条件对式(8)积分，得出油相的流速

(9)

图4(b)显示了不同滑移长度下，通过数值模拟和理论值式(9)得到的连续相通道内油相沿径向的速度分布.图中可见数值模拟和理论公式得到的油相速度几乎相同，速度沿连续相通道轴线呈抛物线对称分布，各速度曲线交于公共点，交点分别在y1=0.026mm和y1=0.083mm，说明两交点距壁面占通道的46.8%左右，得到连续相通道内壁面滑移作用的范围.随滑移长度增加，油相速度在通道中心线上取得的最大值相对于油相入口速度由1.348减至1.24，而壁面上的滑移速度相对于油相入口速度由0.301增至0.519，这些结果表明，滑移长度的增加促进了通道内连续相均匀流动.

图4　连续相通道内建立的x1-y1坐标系和油相速度分布随滑移长度的变化

由滑移壁面条件知，连续相通道壁面对油相施加的黏性切应力

(10)

图5　壁面切应力随滑移长度的变化

式中：us是壁面的滑移速度.不同滑移长度下，由式(10)得到壁面相对切应力如图5所示.图中可见，滑移长度增加，油相受到壁面黏性切应力减小，减小程度随滑移长度的增加而减弱，滑移长度由14 μm继续增加，该切应力随滑移长度由二次减小变为线性减小.由此表明，滑移长度增加，造成油相受到壁面的黏性阻力减弱，即滑移壁面具有减阻的效果，但滑移壁面减阻作用会随着滑移长度的增加而减弱.

2.2固壁分子作用对液滴长度的影响

宏观尺度下流体黏性称为常规黏性，由分子间的吸引力决定，其大小与分子间吸引力成正比.微尺度流动中，液体黏性除常规黏性外，还有一部分是由壁面和近壁面处液体分子间相互吸引而产生的附加黏性，受壁面影响的液体黏性系数为表观黏性系数μa

μa=μ+Δμ

(11)

式中：μ为宏观尺度下流体的黏性系数，Δμ为附加黏性系数.

附加黏性系数随流体和壁面间距减小而增加.对于特征尺度较大的流动通常忽略壁面对流体分子吸引产生的黏性.微尺度的流动中不能忽略这种壁面的吸引作用，在二维微流动中，考虑到通道壁面对流体的共同作用，将附加黏性系数定义为

(12)

式中无量纲数ξ和n由流体和固壁确定，d1是流体质点到通道轴线的距离，为了防止壁面上流体黏度无限大，取a的值为1.08.可以看出，附加黏性系数受n和ξ影响.分别改变两参数值，以详细研究表观黏性对微通道内液滴尺寸和流动的影响规律.

2.2.1n对液滴长度和连续相流动的影响

为研究参数n对液滴长度的影响，连续相通道滑移长度是8 μm的通道内，固定无量纲数ξ为0.1，n分别取1、2、3，得到长度分别为3.02D1、3.01D1、2.94D1的液滴.由此可见，参数n增加造成液滴长度缩短.

图6　连续相通道内速度分布随参数n的变化

对于无附加黏度和不同n对应的附加黏度，连续相通道y1轴上油相的速度分布变化如图6所示.图中可见，无附加黏度，速度分布呈抛物线状；随着n增加，壁面上的速度由0.218uc减小至0.008uc，连续相通道内油相速度最大值明显增加，由1.449uc增至2.037uc，速度分布不均匀程度增加，油相流动越来越集中在连续相通道的中间区域.图中的速度曲线基本在相同位置相交，交点分别位于y1=0.03 mm和y1=0.081 mm，交点至壁面的距离占通道的54.1%左右，得到了连续相通道内受参数n影响的油相流体表观黏性作用范围.

图中显示n大于1，油相速度分布曲线上出现拐点，两拐点间的区域为油相主流区，油相在主流区速度分布呈抛物线状，在拐点处做主流区速度分布曲线的切线，该切线和纵轴负半轴相交，表明此时连续相通道壁面上的滑移速度为负值，定义此边界为负滑移边界，它在模拟宏观和纳米尺度的流动中有较多应用.图中可见，n增大，拐点逐渐远离通道壁面，造成油相在壁面上负滑移的程度增大，油相的主流区变窄.速度分布曲线上拐点的出现，导致油相的流动容易失稳.n越大，越容易使流动失稳，造成n取3时，形成的液滴长度最短.

2.2.2ξ对液滴长度和连续相流动的影响

为考察ξ对液滴长度的影响，固定n为2，连续相通道壁面滑移长度为8 μm，参数ξ取不同值得到不同长度的液滴，如图7所示.图中可见，随着ξ由0增至0.1，液滴长度由3.09D1减小至3.01D1；ξ由0.1增加到0.5，液滴长度仅由3.01D1减小至2.99D1；ξ继续增加，液滴长度几乎不变.

图7　液滴长度随参数ξ的变化

图8　连续相通道内速度分布随参数ξ的变化

改变ξ造成连续相通道内y1轴上油相速度分布发生变化，如图8所示.图中可见，速度分布曲线基本在y1=0.032 mm和y1=0.079 mm处相交，表明连续相通道内受ξ影响的油相流体表观黏性作用范围是57.7%，比参数n的影响范围略大.随着ξ增加，油相在壁面上的速度由0.301uc减小至0.028 4uc，通道中心线上的最大速度由1.35uc增至1.972uc，但速度分布曲线的变化逐渐减弱，当ξ由0.5增至0.6，速度分布几乎不变.ξ大于零，速度分布曲线出现拐点，随着ξ增大，拐点远离通道壁面，造成油相在壁面上负滑移程度增长，油相的主流区变窄，但ξ由0.5增至0.6，油相的负滑移程度几乎不变，由此造成ξ由0.5继续增加，液滴长度基本相同.

3　结语

基于水平集法，对Y型微通道内水-油两相流动进行数值模拟，通过改变连续相通道的滑移长度和油相表观黏度，得到下面的结果：

1)滑移长度增加，促进了液滴长度和壁面滑移速度的增加，使连续相通道内的速度分布更加均匀.连续相通道内壁面滑移作用范围占通道的46.8%.

2)增加影响油相表观黏度的参数n或ξ，液滴长度均缩短.参数n大于1或ξ大于零，连续相通道出现负滑移边界，通道内流动失稳，促进液滴脱离.

3)连续相通道内受参数n影响的油相流体表观黏性作用范围占连续相通道的54.1%，受参数ξ影响的作用范围略大，占57.7%.

[1]WANGWei,XIERui,JUXiao-jie,etal.Controllablemicrofluidicproductionofmulticomponentmultipleemulsions[J].LabonaChip,2011,11(9):1587-1592.

[2]JANKOWSKIP,OGOCZYKD,DERZSIL,etal.Hydrophilicpolycarbonatechipsforgenerationofoil-in-water(O/W)andwater-in-oil-in-water(W/O/W)emulsions[J].MicrofluidicsandNanofluidics,2013,14(5): 767-774.

[3]PAWELJ,DOMINIKAO,LADISLAVD,etal.Hydrophilicpolycarbonatechipsforgenerationofoil-in-water(O/W)andwater-in-oil-in-water(W/O/W)emulsions[J].MicrofluidNanofluid,2013,14(1):597-604.

[4]ZHAOChun-xia,ANTONPJ.Two-phasemicrofluidicflows[J].ChemicalEngineeringScience,2011,66:1394-1411.

[5]VANSTEIJNV,KLEIJNCR,KREUTZERMT.PredictivemodelforthesizeofbubblesanddropletscreatedinmicrofluidicT-junctions[J].LabChip,2010,10:2513-2518.

[6]WANGXiao-da,ZHUChun-ying,WUYi-ning,etal.DynamicsofbubblebreakupwithpartlyobstructioninamicrofluidicT-junction[J].ChemicalEngineeringScience,2015,132:128-138.

[7]袁希钢,宋文琦.T型结构微通道气液两相流型的数值模拟[J].天津大学学报,2012,45(9):763-769.

[8]刘赵淼,王国斌,申峰.基于Navier滑移的油膜缝隙微流动特性数值分析[J].机械工程学报，2011,47(21):104-110.

[责任编辑马云彤]

Numerical Simulation of the Two-phase Flow in a Micro-scale Channelunder the Boundary Slip Condition

WANG Lin-lin1,2, LI Ze3

(1.SchoolofInformationEngineering,Xi’anUniversity,Xi’an710065,China; 2.KeyLaboratoryofThermalFluidScienceandEngineeringofMOE,SchoolofEnergyandPowerEngineering,Xi’anJiaotongUniversity,Xi’an710049,China; 3.SchoolofScience,Xi’anPolytechnicUniversity,Xi’an710048,China)

InordertoinvestigatetheinfluenceofthewallslipandapparentviscosityofthecontinuousphaseonthedropletsizeandflowcharacteristicsinaY-junctionmicro-channel,thelevelsetmethodisusedfornumericalsimulation.Theresultsshowthatwhenthesliplengthincreasesinthecontinuousphasechannel,thedropletlengthincreases,andthevelocitydistributionbecomesevenly.Therangeoftheboundaryslipeffectaccountsfor46.8%ofthechannel.Theapparentviscosityisdeterminedbythetwoparameters.Thenegativeboundaryslipisformed,thedropletdetachmentisaccelerated,anddropletlengthdecreasesforincreasingthetwoparameters.Therangeoftheapparentviscosityeffectaccountsformorethanhalfofthechannel.Itprovidesareferencefortheoptimaldesignofthemicrofluidmixingsystem.

micro-channel;numericalsimulation;levelsetmethod;dropletlength;boundaryslip;apparentviscosity

1008-5564(2016)03-0028-06

2016-01-04

国家自然科学基金资助项目(51076126);西安市科技计划项目(CXY1443WL19)

王琳琳(1981—)，女，河南长垣人，西安文理学院信息工程学院讲师，博士，主要从事微通道内两相流动研究.

TQ021.1

A