T型结构微通道气液两相流型的数值模拟

袁希钢，宋文琦

(天津大学化工学院化学工程联合国家重点实验室，天津 300072)

研究者通过大量的实验观测发现，随着气相和液相表观流速的不断变化，微通道内会出现泡状流、Taylor流、弹状-环状流、环状流、波状流和搅拌流等几种气液两相流型[1-6]．许多研究者[6-10]根据实验数据得出微通道中两相流的流型分布受到惯性力、通道几何形状、通道几何尺寸、通道壁面润湿性、流体的物性和通道入口处条件等因素的影响，已提出的有以气液相表观流速为坐标的流型图，以气液相韦伯数为坐标的流型图，以气含率为划分依据的流型图和流型转变的无因次经验公式[11-14]．Qian和 Gupta等[15-16]通过数值模拟对微通道气液两相流特性进行研究，但没有对气液两相流型分布及其影响因素进行分析和预测．黏度和表面张力的改变对气液两相流型变化的实验研究较少，仅得到一些规律[4,9,13-14]．笔者通过建立适当的计算模型，采用 FLUENT软件对T型微通道气液两相流进行数值模拟，对气液相表观流速和流体黏度、表面张力等物性对两相流型的影响进行了研究，给出气-液两相流型图，分析流体物性对流型图的影响，并与文献报道的实验结果进行对比．

1 T型结构微通道的模型化

1.1 几何模型和控制方程

图 1是 T型结构微通道数值模拟的几何模型．气相和液相分别从气相入口和液相入口进入系统，在两相接触器中间部位汇合后流入微通道中，最后由出口处流出．微通道的无量纲截面宽度为 d，气相和液相的入口处总长度为 6,d，微通道混合区域长度为30,d．

图1 T型结构微通道的几何模型Fig.1 Geometry model T-junction of microchannel

为了更好地对 T型结构微通道流型变化进行数值模拟研究，做出如下假设：

(1) 由于矩形截面的微通道尺寸非常小，使得雷诺数一般小于 200，因此对微通道中的流动按层流流动处理；

(2) 由于微通道内的压力降一般较小，故气相和液相一般视作不可压缩流体，即气相和液相的密度为恒定值；

(3) 气相和液相的物性组成在微通道中视为常数；

(4) Hazel和 Edvinsson等[17-18]通过研究指出，只有当Bo远远大于1时，重力对系统的影响将会很显著．本文模拟条件下，Bo＜1，所以模拟过程中不考虑重力的影响；

(5) 模拟过程中设定为常温条件．

本文采用了有限体积法对微通道中气液两相流进行模拟研究，应用多相流VOF(volume of fluid)模型[15-16,19]对微通道中气液两相间的相界面进行捕捉．那么应用 VOF模型对多相流进行模拟计算的控制方程如下所示.

VOF方法是通过追踪计算单元内各相体积分率的分布来得到相界面的位置，由于只考虑气液两相，因此气相体积分率方程如式(3)所示．在整个求解域内只求解一套控制方程，因求得的速度场和压力场是各相共享的．出现在控制方程中的物性参数也是由单元网格上的各相物性参数共同决定．因此式(1)和式(2)中的混合密度和混合黏度计算式为

对于气液两相流来说，气相的体积分率可以通过式(3)获得，由于在每个控制容积内，所有的相含率和为1，因此液相的体积分率为

气液两相界面处由于物性参数的突然改变，会受到表面张力的影响．表面张力是否对流动产生影响在数量上的判断主要取决于 2个无量纲数：雷诺数Re和毛细数Ca，或雷诺数Re和韦伯数We．

当 1Re≪ 时，如果毛细数 1Ca≫ ，表面张力效应可以忽略；当 1Re≫ 时，如果韦伯数 1We≫ ，表面张力效应可以不考虑．本文中气相和液相表观流速变化较大，这使得 Re数的波动较大，因此不能忽略表面张力的影响．笔者采用 Brackbill等[20]提出的连续表面张力(continuum surface force，CSF)模型，使用散度定理将表面张力表示为体积力，此体积力为动量守恒方程(2)中的源相．

1.2 初始条件和边界条件

对微通道两相流的模拟计算，就是对控制方程由偏微分方程组离散后成线性方程组，然后按照设定的迭代步长求解，因此计算时需要对整个流场进行初始化．

(1) 初始时刻假定整个微通道被液相全部充满，即

(2) 初始时刻各个网格单元的速度为0，即

(3) 假定计算域中全局表压为0．

由于 T型结构微通道的物理模型存在气、液相入口和一个共同的出口，所以需要分别给出气相和液相在入口处的速度、相含率和出口处条件．

(1) 液相入口．

速度：uL＝uL，in，方向如图 1 所示．

相含率：αG＝0，αL＝1．

(2) 气相入口．

速度：uG＝uG，in，方向如图 1 所示．相含率：αG＝1，αL＝0．

(3) 出口及壁面边界条件．

出口处为环境压力，pout＝0(表压)；壁面处设置为无滑移，即紧贴壁面处流体的速度与壁面的速度相等；壁面接触角为 0，即认为壁面是被液相完全浸润的．

1.3 数值求算方法

本文的数值模拟研究是在曙光 TC4000工作站上进行的，该工作站包含 8个刀片，每个刀片装有 4个 AMD Opteron型号的 CPU，该 CPU的频率为2.2,GHz，双核，每个刀片的内存为 8,G．在模拟计算中选用的求解器为基于有限体积法的CFD商用软件FLUENT6.3.26，其运行环境为 Linux系统．采用FLUENT自带前处理软件 GAMBIT2.4对物理模型进行网格划分，然后将网格文件导入求解器FLUENT中，通过求解器得到最终结果，使用FLUENT自带后处理器进行分析和讨论．

为了便于与Kawaji等的实验结果[5]进行对比，取微通道的当量直径d为0.53,mm，采用均匀的网格划分方法，其最小单位的尺寸为当量直径d的1/25．因为整个流动为动态过程，因此采用非稳态求解器进行模拟计算．为了减小计算过程中的存储误差，采用双精度进行求解．在 FLUENT6.3.26的使用手册中[21]，对于气液两相流模拟在多相流模型中的体积力公式推荐使用隐式质量力，VOF模型中的体积分数守恒方程采用显式算法．表1给出了气液两相在模拟过程中所使用的物性参数值．控制方程的离散化选项中，压力项采用 PRESTO(pressure staggering option)算法，动量方程采用二阶迎风格，气液相界面的处理采用几何重构技术，压力-速度耦合方程采用 PISO(pressure-implicit with splitting of operators)算法．初始时刻各个变量的收敛残差值均设置为10-12，计算后适当进行调整，但始终小于 10-6．时间步长在 10-5～10-8,s之间选取[22]．

表1 在模拟过程中气液两相的物性参数(20,℃)[23]Tab.1 Physical parameters of gas and liquid used in simulation(20,℃)[23]

2 模拟结果与讨论

2.1 模拟结果合理性验证

为了验证不同网格密度对模拟结果的影响，采用均匀的网格划分方法，设定最小单元的网格尺寸分别为当量直径 d的 1/10、1/16、1/25和 1/32．图 2显示了在不同网格密度下同一时刻微通道中空气-水两相流流动情况．从图 2可知，选取的网格密度对流型发展没有影响，但在不同网格密度下，对气液相界面的模拟结果是有影响的，网格密度越小，气液相界面越模糊．因此考虑到模拟的精度和计算时间，采用四方形网格数为 22,525的网格密度对微通道中的气液两相流进行数值模拟．

根据本文中建立的计算模型，在不同的气-液相表观流速下进行了数值模拟，得到有不同流体特征的几种气液两相流型．如图 3(a)和(b)所示，在液相表观流速相对较大而气相表观流速相对较小的情况下，气体以小气泡的形式分布在液体中，气泡直径一般小于管径，该流型为泡状流；随着气相表观流速的增大，气泡在微通道内的含量在增加，在微通道内形成类似子弹形状的气泡(图 3(c)和(d))，因此该流型为弹状流或者 Taylor流；当气相表观流速继续增大，在微通道内会形成稳定的气柱，即为环状流，如图 3(e)和(f)所示；在气相和液相表观流速都很高的情况下，气液相的流动在微通道内会有比较大的波动，导致有微小液滴进入气泡中，同时也会有小碎气泡进入液相，这时搅拌流就会出现(图3(g)和(h))．

图2 不同网格密度下空气-水两相流流型的模拟结果(d＝0.53,mm，uG＝0.1,m/s，uL＝0.15,m/s，t＝0.053,81,s)Fig.2 Simulation results of air-water two phase flow at different mesh resolutions(d＝0.53,mm，uG＝0.1,m/s，uL＝0.15,m/s，t＝0.053,81,s)

Triplett、Waelchli和 Kawaji等[1，4，6]对 T 型结构的微通道进行过气液两相流的实验研究，观察到了泡状流、Taylor流、环状流、弹状-环状流和搅拌流等不同流型．本文对 T型结构微通道内空气-水体系两相流型的模拟结果与文献[1]在水平放置的1.097,7,mm三角形截面微通道内的空气-水两相流的实验观测到的流型结果非常吻合，模拟结果能够准确反映出泡状流、Taylor流、环状流、搅拌流的基本流体特征．

图3 不同气液相表观流速下流型图的模拟结果Fig.3 Representative pictures of flow patterns in different gas-liquid superficial velocities

将T型结构微通道内空气-水体系的两相流流型的模拟结果，即不同气液相表观流速下微通道内数值模拟的流型结果，绘制成横纵坐标分别为气相表观流速(uG)和液相表观流速(uL)的 uG-uL流型图．如图 4所示，根据不同气液相表观流速的数值模拟结果，可以给出不同流型间的过渡线，与 Kawaji等[6]在0.53,mm 圆柱形微通道中氮气-水两相流实验得出的流型间的过渡线进行对比，可以看出 Taylor流向搅拌流和环状流的过渡线与Kawaji实验中得出的过渡线比较吻合．在模拟中所建立的几何模型与对比实验中所使用的入口处结构不同，再加上模拟过程中不能对通道壁面材料的物性进行模拟，这使得模拟结果得到的泡状流与 Taylor之间的过渡线与对比实验的过渡线有偏差．

图 4 文献[6]的实验流型图过渡线与本文模拟的流型图过渡线对比Fig.4 Comparison of transitions of flow regimes between simulated results and experimental data[6]

通过改变气相和液相表观流速，模拟出空气-水两相流各流型的流动特征，进而与 Triplett等[1]和

Kawaji等[6]的两相流的实验数据进行对比，从图3的两相流模拟结果和图 4的流型图实验对比结果可以看出，模拟结果与对比实验结果吻合得很好，说明本文所建立的计算模型和选取的求解算法是合理的，能够准确地反映微通道内两相流的实际流动情况．

2.2 黏度对两相流流型影响

选取了与已有实验不同的气液相流体，从模拟的角度考察黏度对气液两相流流型产生的影响．依据本文所建立的几何模型和选取的算法，对液相质量分数分别为 20%、60%和 80%的丙三醇溶液，气相为空气的3组气液相流体进行模拟．从表1可以看出，随着丙三醇溶液质量分数的不断增加，表面张力的变化非常小，可以近似认为是恒定的，但运动黏度却在不断增大，因此可以通过这3组不同的气液相流体衡量黏度变化对两相流流型产生的影响．将 3组不同组成的气液相流体的 uG-uL流型图绘制在一起，并分别画出 3组气液相流体各个流型的过渡线，从图 5可得，随着模拟过程中黏度的不断增加，Taylor流与泡状流之间的过渡线基本没有发生变化，而 Taylor流向搅拌流和环状流的过渡线在模拟结果中也是随着黏度的增加向更高的气相表观流速(uG)方向移动，但移动得非常小，可见黏度变化对空气-丙三醇溶液的气液两相流型变化影响很小．这与 Zhang等[14]用水和不同浓度的三甲基纤维素钠分别为液相，氮气为气相流体在Y型结构的微通道(当量直径496,µm)内进行实验得出的流型变化的结论是一致的．但从图 5可以看出，不同质量分数丙三醇溶液的黏度变化对Taylor流向搅拌流和环状流的过渡线影响是非常小的，而 Zhang等[14]实验结果中过渡线的移动距离要比本文模拟结果的距离大，这种差距可能是模拟和实验中使用了不同的气液相流体，物性黏度的差异造成的．

图5 黏度的改变对T型结构微通道两相流流型图的影响Fig.5 Influence of viscosity on the gas-liquid two phase regimes in T-junction microchannel

2.3 表面张力对两相流流型影响

选取了水、质量分数为 10%的乙醇溶液和纯乙醇液体分别与空气组成气液两相流体，从模拟的角度考察表面张力对气液两相流型产生的影响．依据所建立的模型对这 3组气液相流体进行数值模拟．从表 1可以看出，随着水、质量分数 10%的乙醇溶液和纯乙醇液体表面张力的逐渐减小，而液相流体运动黏度的变化非常小，可以近似认为是恒定的，因此可以通过这 3组不同的气液相流体衡量表面张力变化对两相流流型产生的影响．如图6所示，将3组模拟结果的uG-uL流型图绘制在一起，并分别画出 3组气液相流体各个流型的过渡线，随着水、质量分数 10%的乙醇溶液和纯乙醇液体表面张力的减小，泡状流的区域逐渐向Taylor流扩张，搅拌流的区域逐渐向Taylor流和环状流扩张，环状流的区域也逐渐向 Taylor流扩张，而 Taylor流的区域则随着表面张力的减小而缩小．

Pohorecki等[9]用水和乙醇作为液相流体进行实验来考察表面张力对流型图中过渡线的影响，得知随着表面张力的增加，弹状-环状流向环状流的过渡发生在更小的气相表观流速(uG)下．Waelchli等[4]通过实验得出与 Pohorecki等[9]相反的结论，还发现随着表面张力的增加，泡状流向 Taylor流的转变发生在更小的液相表观流速(uL)下．本文的模拟结果与Zhang等[14]用水、0.060,8%(质量分数)十二烷基硫酸钠(SDS)、0.261,0%(质量分数)十二烷基硫酸钠(SDS)和纯乙醇为液相、氮气为气相流体、在 Y型结构的微通道(当量直径 496,µm)内进行实验得出的流型变化结论是一致的．从模拟结果可以说明，不同质量分数乙醇溶液表面张力的变化对气液两相流型变化的影响很大；随着表面张力的减小，Taylor流的控制区域在减小，这说明表面张力对 Taylor流的形成有较大的影响．

图 6 表面张力的改变对 T型结构微通道两相流流型图的影响Fig.6 Influence of surface tension on the gas-liquid two phase flow regimes in T-junction microchannel

3 结 论

(1) 建立了适用于微通道气液两相流的计算模型，利用 FLUENT求解器模拟出不同气相和液相表观流速下的泡状流、Taylor流、环状流和搅拌流，在今后的科学研究中可以用于微通道气液两相流型的分析和预测．

(2) 对质量分数分别为 20%、60%和 80%的丙三醇溶液进行模拟计算，发现随着丙三醇溶液质量分数的不断增加，黏度也不断的增加，Taylor流与泡状流之间的过渡线基本没有变化，Taylor流向搅拌流和环状流的过渡线随黏度的增加向更高的 uG方向移动．但丙三醇溶液黏度的改变对流型图中过渡线的影响很小．

(3) 通过对水、质量分数为10%乙醇溶液和纯乙醇液体进行模拟计算，得出表面张力的变化对微通道气液两相流型图中的过渡线会产生比较大的影响．随着表面张力的不断减小，Taylor流向泡状流、搅拌流和环状流的区域在逐渐缩小，即表面张力对Taylor流的形成有比较大的影响．

符号说明：

D—微通道宽度，m；

ρ—流体的密度，kg/m3；

t—时间，s；

g—重力，kg/m2；

F—体积力，N；

v —速度矢量，m/s；

α—体积分率；

u—表观速率，m/s；

μ—流体的黏度，Pa·s；

σ—表面张力，N/m；

Bo—邦德数，Bo＝(ρL-ρG)gd2/σ；

Re—雷诺数，Re＝ρL(uG+uL)d/μL；

Ca—毛细数，Ca＝μL(uG+uL)/σ；

We—韦伯数，We＝ρLduB2/σ；

下标：

L—液相；

G—气相．

［1］ Triplett K A，Ghiaasiaana S M，Abdel-Khalik S I，et al. Gas-liquid two-phase flow in microchannels (Part I)：Two-phase flow patterns [J].International Journal of Multiphase Flow，1999，25(3)：377-394.

［2］ Liu H，Vandu C O，Krishna R. Hydrodynamics of Taylor flow in vertical capillaries：Flow regimes，bubble rise velocity，liquid slug length，and pressure drop [J].Industrial & Engineering Chemistry Research，2005，44(14)：4884-4897.

［3］ 马友光，付涛涛，朱春英. 微通道内气液两相流行为进展研究 [J]. 化工进展，2007，26(8)：1068-1074.

Ma Youguang，Fu Taotao，Zhu Chunying. Progress of studies on gas-liquid two-phase flow in microchannels[J].Chemical Industry and Engineering Progress，2007，26(8)：1068-1074(in Chinese).

［4］ Waelchli S，von Rohr P R. Two-phase flow characteristics in gas-liquid microreactors [J].International Journal of Multiphase Flow，2006，32(7)：791-806.

［5］ Gupta R，Fletcher D，Haynes B S，et al. Taylor flow in microchannels：A review of experimental and computational work [J].Journal of Computational Multiphase Flows，2010，2(1)：1-31.

［6］ Chung P M-Y，Kawaji M. The effect of channel diameter on adiabatic two-phase flow characteristics in microchannels [J].International Journal of Multiphase Flow，2004，30(7/8)：735-761.

［7］ Shao N，Gavriilidis A，Angeli P. Flow regimes for adiabatic gas-liquid flow in microchannels [J].Chemical Engineering Science，2009，64(11)：2749-2761.

［8］ Yue J，Luo Lingai，Gonthieret Y，et al. An experimental investigation of gas-liquid two-phase flow in single microchannel contactors [J].Chemical Engineering Science，2008，63(16)：4189-4202.

［9］ Pohorecki R，Sobieszuk P，Kula K，et al. Hydrodynamic regimes of gas-liquid flow in a microreactor channel [J].Chemical Engineering Journal，2008，135：185-190.

［10］ Rebrov E V. Two-phase flow regimes in microchannels[J].Theoretical Foundations of Chemical Engineering，2010，44(4)：355-367.

［11］ Saisorn S，Wongwises S. The effects of channel diameter on flow pattern，void fraction and pressure drop of two-phase air-water flow in circular micro-channels [J].Experimental Thermal and Fluid Science，2010，34(4)：454-462.

［12］ Sur A，Liu D. Adiabatic air-water two-phase flow in circular microchannels[J].International Journal of Thermal Sciences，2012，53：18-34.

［13］ Yang C-Y，Shieh C-C. Flow pattern of air-water and two-phase R-134a in small circular tubes [J].Int J Multiphase Flow，2001，27(7)：1163-1177.

［14］ Zhang Tong，Cao Bin，Fan Yilin，et al. Gas-liquid flow in circular microchannel(Part I)：Influence of liquid physical properties and channel diameter on flow patterns [J].Chemical Engineering Science，2011，66(23)：5791-5803.

［15］ Qian D，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.

［16］ Gupta R，Fletcher D F，Haynes B S. On the CFD modelling of Taylor flow in microchannels [J].Chemical Engineering Science，2009，64(12)：2941-2950.

［17］ 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-114.

［18］ Edvinsson P K，Irandoust S. Finite-element analysis of Taylor flow [J].AIChE Journal，1996，42(7)：1815-1823.

［19］ Kumar V，Vashisth S，Hoarau Y，et al. Slug flow in curved microreactors：Hydrodynamic study [J].Chemical Engineering Science，2007，62(24)：7494-7504.

［20］ Brackbill J U，Kothe D B，Zemach C. A continuum method for modeling surface tension [J].Journal of Computational Physics，1992，100(2)：335-354.

［21］ FLUENT 6. 3 23. 3. 1 Overview and Limitations of the VOF Model [Z].

［22］ 黄 悦. 带有化学反应的微通道气液两相流体力学特性模拟研究 [D]. 天津：天津大学化工学院，2010.

Huang Yue．Numerical Study on Fluid Dynamic Characteristics of Gas-Liquid Two-Phase Flow with Chemical Reaction in Microchannel[D]. Tianjin：School of Chemical Engineering and Technology，Tianjin University，2010(in Chinese)．

［23］ 刘光启，马连湘，刘 杰，等. 化学化工物性数据手册：有机卷 [M]. 北京：化学工业出版社，2002.

Liu Guangqi，Ma Lianxiang，Liu Jie，et al．Handbook of Physical Property Data for Chemistry and Chemical Engineering(Vol Organic)[M]. Beijing：Chemical Industry Press，2002(in Chinese)．