水力学聚焦微通道中气液两相流动的数值模拟

吕建华，李品，高梦璠，殷明海

（河北工业大学化工学院，天津 300130）

引 言

自20 世纪90年代以来，微化工系统因其传质传递速率快、安全性高、集成度高、可控性强、放大效应小以及过程节能等优点[1]逐渐成为研究的热点。

学者们对微通道内气、液两相流动的流型和流型转变判据[2]、流动行为影响因素[3]、气泡生产机理[4-5]等进行了大量的研究，但这些研究都主要针对泡状流、子弹流，针对环状流（膜状流动）的研究较少。

随着高性能计算技术（HPC）和仿真技术的快速发展，数值模拟在微化工领域的应用越来越广泛。在建立合理的模型、提供正确的控制方程、边界条件和求解参数的基础上，数值模拟得到的数据可以很好地吻合实验结果[3,6-8]，这为微化工领域的发展提供了一种可靠而高效的研究方法。

然而，关于微设备内两相流的实验研究、数值模拟等都主要集中在了T 型[6,8-9]、Y 型[10-11]和十字形[12]微通道领域，有少数研究者[13-16]对聚焦型微通道内的流体流动进行了研究，但主要考察了气泡（液滴）的生产机制、粒径的影响因素等，而平行流（膜状流动）研究更为少见。Austin 等[17]通过实验发现，水力学聚焦微通道内可以将液体聚焦成薄膜，但研究者们仅研究了互不相容液液两相的流动状态，微通道内气、液两相的流动却报道较少。鉴于气、液两相的黏度、密度、表面张力等性质有很大差异，本文拟采用FLUENT 软件，以空气-水为介质，模拟气、液两相在水力学聚焦微通道内的流动状态，以确定气、液两相的流动规律。

1 建立模型

1.1 物理模型

锥形缩孔是Austin 等[17]开发的水力学聚焦微设备中的一个重要设计，也是流体可以被聚焦成薄膜的关键所在，但文献中并未给出微设备的全部尺寸，因此，本文根据文献[17]提供的数据设计了如图1 所示的水力学聚焦微通道模型。微通道有3 个入口和1 个出口，最左侧是气相（空气）入口，上下两端是液相（水）入口，右侧为混合相出口，气相通道终端是一个锥形缩孔。空气与水的物性参数列于表1 中。气相入口、液相入口与出口的微通道宽度均为10 μm，气相通道A-B（A′-B′）和液相通道E-F（E′-F′）的长度均为20 μm，两相混合区长200 μm，与气相通道连接的锥形缩孔B-C（B′-C′）长10 μm，缩孔最终的宽度（C-C′）是2 μm。气液 两相在正方形微通道内呈中心对称形状[18]，因此采用二维模型进行数值模拟。

图1 水力学聚焦微通道示意图Fig.1 Schematic diagram of hydrodynamically focused microchannels

表1 空气、水物性参数Table 1 Physical parameter of air and water

在上述模型的基础上，本文做出以下假设，以更好地对水力学聚焦微通道内气、液两相流进行数值模拟：① 本模型的Bond 数（Bo）小于1，可忽略重力影响[19-20]；②微通道尺寸很小，通道内的流体都处于层流状态[21]；③微通道内的压降值很小，气、液两相均按不可压缩流体处理；④本模型只研究气、液两相流动型态和水力学相关参数，不考虑热能变化。

1.2 控制方程

多相流VOF（volume of fluid）模型，可以在固定的欧拉网格下进行表面追踪，捕获到不相溶流体间的交界面。模型的控制方程如下：

连续性方程

动量方程

体积分数方程

式中，u 是速度矢量；Fσ是方程的源相，采用CSF 方法得到的表面张力的体积力形式[22]；αi是相的体积分数(i 为G 和L 时分别表示气相和液相)；ρ和μ 分别为混合相的密度与黏度，计算公式如下

1.3 边界条件与初始条件

1.3.1 边界条件 气相入口的边界条件类型为压强入口（pressure-inlet，PG），相含率αG=1；液相入口的边界条件类型为压强入口（pressure-inlet，PL），相含率αG=1；出口的压力条件类型为压强出口（pressure-out），压强值为环境压强（101325 Pa）。

1.3.2 壁面边界条件 无滑移壁面，壁面接触角为0°。

1.3.3 初始条件 全场初始化从气相入口开始，即模拟开始时，计算域内充满气相。

2 模型求解与验证

2.1 数值模拟方法

本文的数值模拟是在惠普Z620 工作站上进行的，该工作站包含2 个Intel Xeon E5-2670 型号的CPU，每个CPU 频率为2.6 GHz，8 核，双线程。首先采用Gambit 软件，选用结构化的四方形网格对计算域进行网格划分；然后用FLUENT 6.3.26 软件进行求解，求解器的设置如下：非稳态项采用一阶隐式时间步进行处理，时间步长设定为10-7s，压强与速度的耦合选用PRESTO！算法，压强、动量、相界面的离散分别选用PRESTO！、二阶迎风、几何重构方法；最后用Tecplot 软件对数值模拟得到的数据进行处理分析。

2.2 网格无关性检验

本模型采用二维四方形结构网格，网格划分类型为平铺，为查看网格尺寸对计算结果的影响，分别设定最小网格尺寸为通道宽度的1/10、1/20、1/40和1/80，网格个数依次为2768、22726、89507 和355404。设定气、液两相入口压强值分别为PG=10 Pa（表压）、PL=10 Pa（表压），分别进行模拟计算，取各个网格密度下模拟时间为0.35 s 时的数据分析，如图2～图4 所示。

图2 气、液两相相分布Fig.2 Phase distribution of gas-liquid two-phase

图3 通道中心轴压强分布Fig.3 Pressure distribution of inlet channel symmetry axis

图4 通道中心轴处X 向速度分布Fig.4 X-velocity distribution of channel symmetry axis

图2～图4 分别为4 种网格密度下的气液两相分布（红色代表空气，蓝色代表水）、通道中心轴（Y=0 μm）处的压强分布和气相的轴向（X 向）速度分布。由图中可以看出C、D 两种网格密度得到的两相流动型态、压强分布、速度分布趋于一致，由此可认为小于通道宽度1/40 的网格密度对计算结果影响不显著；随着网格的加密，达到稳定流态所需的计算时间急速增加，因此，本文采用网格数为89507 个的网格密度进行计算。

3 模拟结果与讨论分析

从图2 可以看出，空气和水进入微通道后在两相混合区都呈稳定的膜状流动，这与文献[17]得到的实验结果相符合。验证了模型的准确性。

3.1 气、液两相成膜过程分析

当PG=10 Pa（表压）、PL=10 Pa（表压）时，经过数值模拟得到的气、液两相在聚焦区的流动与发展过程如图5所示（坐标轴对应微通道X、Y值，μm）。图中红色代表空气，蓝色代表水，黑色细实线为等压线，不同的压强值用数字标示。通过分析，本文将两相成膜的过程分为以下3 个阶段。

（1）气相受阻阶段[图5(a)、(b)]，随着液相从入口不断注入，液相在微通道中不断发展，置换了通道中的气相最终进入聚焦区，气相流通受阻。

（2）液相偏移阶段[图5(c)、(d)]，液相继续沿纵向发展，同时液相前段向两相混合区偏移，随后进入两相混合区。为方便研究，称进入两相混合区的液相前段为头部，暴露在聚焦区的液相为颈部。

（3）两相成膜阶段[图5(e)、(f)]，液相在纵向的发展到达极限，而在两相混合区的发展有所增长，液相贴着微通道壁面不断置换气相，最终呈膜状流动，气相在通道中心处形成气膜。

图5 气、液两相在聚焦区的发展过程Fig.5 Propagation of gas-liquid phase in focused section[(a),(b) stage one; (c),(d) stage two; (e),(f)stage three]

3.1.1 压强变化过程分析 表2 列出了对应图5 中各个时刻液相表面的压强值，其中(a)～(c)对应的是液相前端左右两侧的压强，(d)～(f)对应的是液相头部与颈部表面的压强。

表2 液相表面压强值变化过程Table 2 Pressure change of liquid surface

由表2 可知，随着液相在微通道内流动，液相前段左侧（颈部）压强先增大后稍微减小，而右侧（头部）的压强不断减小，液相内部的压差推动着液相不断向出口流动。

液相从侧通道入口进入后沿纵向在液相微通道中发展，随着液相量不断增加，液相前端首先到达聚焦区[图5(a)]，此时，液相前端左右两侧的压强值相等，随后液相进入聚焦区内，气相流通受阻，致使气相积压，液相前端左侧压强略大于右侧[图5(b)]。不断注入的液相继续阻塞气相通道，加大了液相前端左右两侧压强的差值，在压差的驱动下，液相前端向两相混合区偏移[图5(c)]，随后进入两相混合区内[图5(d)]，形成了头部和颈部。如图5(e)所示，液相沿纵向的发展达到极限，此时颈部压强达到最大值，之后液相在两相混合区中呈膜状向前流动[图5(f)]，缓解了对气相的阻塞，颈部压强也略有降低。进入微通道内的液相量越多，推动液相发展所需的压强差值就越大，因此液相前端右侧和头部的压强呈现出一直减小的趋势。

在这个过程中，压强差一直推动着气、液两相界面不断发展并形成膜：图5(b)、(c)过程的压差阻碍了两个侧通道的液相在聚焦区汇合，保证了气相的连续性，因此形成气膜；图5(d)、(f)过程中，F（F′）点周围的液相等压线呈辐射发散状，液相在此压差下会不断由液相微通道流入两相混合区，而不会被气相在F（F′）点夹断，保证了液相的连续性，形成液膜；两相混合区的液相压强梯度平行于流动方向，液膜沿轴向流动，不会切断气膜。

3.1.2 液相表面黏性力变化过程分析 液相表面所受的黏性力与速度梯度呈正比，黏性力会影响液相变形[23]，从而对通道内气、液两相的流动状态产生影响。

图6 气、液两相轴向速度分布Fig.6 X-velocity distribution of gas-liquid two-phase

气相在微通道中流通受阻时流速会增加，速度梯度也会增大，为查看黏性力对成膜过程的影响，本文绘制了对应图5 中不同时刻气相流通空间最小 处两相的轴向（X 向）速度曲线，如图6 所示。由图中可看出，在整个过程中，随着|Y|减小，两相轴向速度逐渐增大，且在两相交界面处存在速度拐点，气相的速度大于液相，因此液相所受的黏性剪切力为X 正向，能推动液相向两相混合区发展。气相的轴向速度在a～d 阶段一直是增大的趋势，e 时刻减小，而后又略增大，d 时刻有最大值；液相的速度在a～e 阶段一直增大，f 时刻减小；两相界面处的速度梯度也在d 时刻达到最大值。

气、液两相轴向速度变化的原因是：液相刚进入聚焦区时，气相因流通空间被阻塞而速度增大，两相速度梯度也随之增大，在黏性力的作用下，液相的速度也增大，液相前端发生偏移；在d 时刻，液相前端刚好到达两相混合区的入口，此时液相在聚焦区积攒量最多，气相受阻也最严重，因此气相速度值最大，速度梯度最大，黏性剪切力也最大，在此作用下液相的轴向速度继续增大，液相迅速涌入两相混合区；在e 时刻，液相进入两相混合区，缓解了对气相的阻塞，气相速度略有下降，此时两相速度差值小，黏性作用力小；之后的f 时刻，气、液在两相混合区中以膜状流动并趋于稳定，速度也达到相应的稳定值。

3.1.3 表面张力变化过程分析 液相表面张力与液相表面曲率有关，曲率越大则表面张力越大。如图5 所示，(a)～(d)阶段，液相前端曲率逐渐增大，表面张力增大，阻碍了来自两个侧通道的液相在聚焦区汇合，保证了气相的连续性，促进气相成膜。而(e)～(f)阶段，液相头部与颈部曲率较小，所受表面张力很小。两相成膜后，膜界面与通道壁平行，曲率几乎为零，表面张力的作用可忽略。

3.2 压强对两相流量和膜厚度的影响

定义气、液两相压强比值为β

3.2.1 恒定β 值，改变两相压强 保持β 值恒为1，同时改变PG、PL，得到两相流量和膜厚度（δG、δL），相应数据列于表3。

表3 压强（β 恒为1）对两相流量和膜厚的影响Table 3 Effects of pressure （β was held constant at 1） on two-phase flux and film thickness

由表3 中数据可知，两相入口压强从9.000 Pa增大至11.000 Pa 时，气、液相流量分别增大了23.344%和23.664%。气膜厚度有小幅度（1.582%）的减小，而液膜厚度有小幅度（0.532%）的增加。两相膜厚度的微小变化是因气、液相流量的不同步增长导致的，由于通道总体积的守恒，流量涨幅大的液相的膜厚度必然会相应增大[17]。

3.2.2 恒定气相压强，改变β 值 保持PG为10.00 Pa 恒定，改变β，得到相关数据如表4 所示。

表4 β 对两相流量和膜厚的影响Table 4 Effects of β on two-phase flux and film thickness

由表中数据可知，随β 的增大，气相流量逐渐减小，液相流量逐渐增大，且气相流量的变化幅度远大于液相，相应地，气膜厚度也大幅度减小。

对比表3 和表4 的数据可知，保持β 值恒定不变，同时改变气、液相的入口压强对两相的流量和膜厚度的影响不明显，而保持PG不变，改变β 值，两相流量和膜厚度则变化显著。

3.3 β 值对流动型态的影响

Austin 等[17]指出，对特定的微通道结构，只有当0.48<β<1.28 时，两相可以形成膜状流动，β 大于上限值，从液相通道进入的液体通过缩孔进入气相通道中，β 小于上限值时，从气相通道进入的气体会进入液相通道并流出。通过数值模拟发现，本模型β 的范围是0.6<β<1.05，由于本模型的流体介质和微通道尺寸与文献[17]中不同，因此得到的膜状流动的β 范围也不相同。

图7(a)～(d)分别为β 为0.6、0.5、1.05、1.1 时对应的聚焦区两相流动情况。β=0.6 时，液相已经不能形成稳定的膜，液相在聚焦区间歇性地涌入两相混合区内；β=0.5 时气相从液相通道入口流出；而β 增大到1.05 时，气相不再形成连续的膜，而是生成气泡被液相夹带至出口流出；β=1.1 时，气相不能进入聚焦区，液相则通过缩孔进入气相微通道中，最终从气相通道入口流出。图7(a)、(d)所展示的流动状况与文献 [17]的实验结果相吻合，这也进一步说明了本研究模型的可靠性。

图7 非膜状流动的两相分布Fig.7 Phase distribution of non-film two-phase flow

4 结 论

采用FLUENT 软件对水力学聚焦微通道内气、液两相流动进行了数值模拟计算，分析了气、液相接触成膜的过程中液相表面压强、黏性力和表面张力的变化过程，并考察了两相入口压强变化对两相流动的影响，结论如下。

（1）水力学聚焦微通道中，随液相的注入，液相表面左右两侧（头部与颈部）的压强差增大，促进液相由入口通道进入两相混合区；液相颈部压强呈辐射状分布，保证了液相的连续性；液膜内压强沿轴向减小，保证了气相的连续性。

（2）黏性力在成膜过程中先增大后减小，带动液相沿轴向发展，对成膜有重要影响。

（3）表面张力在气相受阻、液相偏移阶段作用显著，在两相成膜阶段中影响较小，可忽略。

（4）保持β 值恒定不变，同时改变PG、PL对两相的流量和膜厚度的影响不明显，而保持PG不变，改变β 值，两相流量和膜厚度则变化显著。

（5）当0.6<β<1.05 时，本模型设计的微通道中两相呈膜状流动，超出此范围则不能形成膜。

[1]Chen Guangwen(陈光文), Zhao Yuchao(赵玉潮), Dong Zhengya(董正亚), Cao Haishan(曹海山), Yuan Quan(袁权).Transport phenomena in micro-chemical engineering [J].CIESC Journal(化工学报), 2013, 64(1): 63-75.

[2]Triplett K A, Ghiaasiaan S M, Abdel-Khalik S I, Sadowski D L.Gas-liquid two-phase flow in microchannels(Ⅰ): Two-phase flow patterns [J].Int.J.Multiphase Flow, 1999, 25(3): 377-394.

[3]Hou Jingxin(侯璟鑫), Qian Gang(钱刚), Zhou Xinggui(周兴贵).Effects of gas inlet angle and cross-section aspect ratio on Taylor bubble behavior in microchannels [J].CIESC Journal(化工学报), 2013, 64(6): 1976-1982.

[4]Garstecki P, Fuerstman M J, Stone H A, Whitesides G M.Formation of droplets and bubbles in microfluidic T-junction scaling and mechanism of bresak-up [J].Lab Chip, 2006, 6(3): 437-446.

[5]Abate A R, Mary P, van Steijn V, Weitz D A.Experimental validation of plugging during drop formation in a T-junction [J].Lab Chip, 2012, 12(8): 1516-1521.

[6]Dang Minhui(党敏辉), Ren Mingyue(任明月), Chen Guangwen(陈光文).Effect of microchannel inlet configuration on Taylor bubble formation in microreactors [J].CIESC Journal(化工学报), 2014, 65(3): 806-812.

[7]Tatineni M, Zhong X.Numerical study of two-phase flows in microchannels using the level set method, AIAA-2004-0929//42nd AIAA Aerospace Sciences Meeting and Exhibit[C].Reno, Nevada, 2004.

[8]Dong Hefei(董贺飞), Zhang Deliang(张德良), Zhao Yuchao(赵玉潮), Chen Guangwen(陈光文), Yuan Quan(袁权).Numerical simulation of immiscible two-phase flow in T-shaped microchannel [J].Journal of Chemical Industry and Engineering(China)(化工学报), 2008, 59(8): 1950-1957.

[9]Deng Jin(邓金), Zheng Xiaobin(郑晓斌), Wu Jizhou(吴纪周), Chen Wen(陈文), Liu Xiong(刘雄), Dong Lichun(董立春), Li Junhong(李俊宏), Xu Yuting(徐玉婷).CFD simulation of droplet formation in microfluidic T-junction [J].Chemical Engineering, 2012, 40(11): 39-43.

[10]Cong Zhenxia(丛振霞), Zhu Chunying(朱春英), Fu Taotao(付涛涛), Ma Youguang(马友光).Bubble breakup and distribution in asymmetric Y-bifurcating microchannels [J].CIESC Journal(化工学报), 2014, 65(1): 93-99.

[11]Qian D Y, Lawal A.Numerical study on gas and liquid slugs for Taylor flow in a T-junction microchannel [J].Chemical Engineering Science, 2006, 61(23): 7609-7625.

[12]Zhao Y, Hemminger O, Fan L S.Experimental and lattice Boltzmann simulation of two-phase gas-liquid flows in microchannels [J].Chemical Engineering Science, 2007, 62(24): 7172-7183.

[13]Gañán-Calvo A M, Gordillo J M.Perfectly monodisperse microbubbling by capillary flow focusing [J].Physical Review Letters, 2001, 87: 274501.

[14]Garstecki P, Gitlin I, DiLuzio W, Whitesides G M, Kumacheva E, Stone H A.Formation of monodisperse bubbles in a microfluidic flow-focusing device [J].Applied Physics Letters, 2004, 85: 2649-2651.

[15]Garstecki P, Stone H A, Whitesides G M.Mechanism for flow-rate controlled breakup in confined geometries: a route to monodisperse emulsions [J].Physical Review Letters, 2005, 94: 164501.

[16]Anna S L, Bontoux N, Stone H A.Formation of dispersions using “flow focusing” in microchannels [J].Applied Physics Letters, 2003, 82: 364-366.

[17]Knight J B, Vishwanath A, Brody J P, Austin R H.Hydrodynamic focusing on a silicon chip: mixing nanoliters in microseconds [J].Physical Review Letters, 1998, 80(17): 3863-3866.

[18]Taha T, Cui Z F.CFD modelling of slug flow inside square capillaries [J].Chemical Engineering Science, 2006, 61 (2): 665-675.

[19]Hazel A L, Hell M.The steady propagation of a semi-infinite bubble into a tube of elliptical or rectangular cross-section [J].J.Fluid Mech., 2002, 470: 91-114.

[20]Edvinsson R K, Irandoust S.Finite-element analysis of taylor flow [J].AIChE J., 1996, 42(7): 1816-1823.

[21]Ehrfeld W, Hessel V, Löwe H.Microreactors: New Technology for Modern Chemistry[M].Beijing: Chemical Industry Press, 2004: 70.

[22]Brackbill J U, Kothe D B, Zemach C.A continuum method for modeling surface tension [J].Journal of Computational Physics, 2000, 162(2): 301-337.

[23]Wang Linlin, Li Guojun, Tian Hui, Ye Yanghui.Numerical simulation of gas-liquid two-phase flow in a T-junction micro-channel [J].Journal of Xi’an Jiaotong University, 2011, 45(9): 65-69.