王志强 刘志宏 张 莹 谢博伟

(南昌大学机电工程学院 江西 南昌 330031)

毛细管内气泡变形及脉动运动模拟

王志强 刘志宏 张 莹*谢博伟

(南昌大学机电工程学院 江西 南昌 330031)

运用VOF模型结合CSF模型对毛细管内气泡的运动变形进行研究分析。首先对毛细管内不同气泡直径的运动变形进行数值计算，考虑气泡受粘性力、惯性力和表面张力的作用，计算得到气泡经历圆形、弹状、箭头状，直至破裂成气泡流；然后对管内气泡运动时气泡破裂进行数值模拟，结果表明不同直径的气泡均存在破裂临界毛细值；最后对毛细管两边施加脉动载荷，分析脉动载荷下气泡运动情况，得出随着脉动压力波的增加，气泡在管内的脉动运动更加剧烈，气泡外形也是随压力的脉动变化发生周期性变化。

毛细管 流形变化 气泡破裂 脉动运动

0 引 言

毛细管具有换热面积大，体积小，传热效率高等优点，使其在民用，医药，建筑领域得到了广泛的应用。电脑中的散热器就是由毛细管组成，人体的血液循环系统就是一个非常庞大且复杂的毛细管网系统。毛细管内一般具有复杂的流体流动，并且随着热量的交换，毛细管内流动往往是伴有相变的两相流动或者是液-液流动。

区别于宏观的流体流动现象，毛细管的当量直径已经发展到微米级甚至纳米级，所以宏观的力学及传热定理不能够精确地描述其传热传质现象，加上毛细管内的流动伴随着相变的发生，在很大程度上影响着流体的传热传质。

Hetsroni[1]通过实验研究加热毛细管内气泡生长的规律及气泡形态的变化。测量了不同加热功率和流体流速情况下管内的压力、温度和速度，利用红外技术获取管内温度分布，获得了管内单相液体，气液两相流和单相气体的流动参数。之后Khandekar等[2]在可视化实验的基础上，利用理论分析不同管道弯头数下管内流型与OHP 倾斜角度等参数的关系，建立了针对闭式振荡热管的数理模型，能够较好地描述振荡热管内的流动状况。但是二者实验中，没有具体分析粘性力、惯性力和表面张力对于实验结果的影响。

基于前期的实验研究，任桂香等[3]对微米级T形管内气泡的生长规律进行了数值模拟研究。结合Matlab 图像处理软件和 Fluent 求解器，进行了大量的模拟。主要从表面张力、流体粘性、气液流速、气液固三相接触角、毛细管管径以及管壁材料物性等方面对毛细管内气泡的生长规律进行了分析，并获得了其所建模型下气泡流形成的临界值(毛细数Ca为0.16)，但他们只研究了气泡的形成和成长，没有涉及到气泡的最终破裂情况。

DehaoXu[4]和Nandan Saha[5]通过ANSYS软件模拟了脉动热管内气泡的生成状态，但是没有考虑添加于脉动有关的UDF程序。

本文运用VOF模型结合CSF模型对毛细管内气泡的运动变形进行了研究分析：首先考虑气泡受粘性力、惯性力和表面张力的作用对毛细管内不同气泡直径的运动变形进行数值计算；然后对管内气泡运动时气泡破裂进行数值模拟；最后对毛细管两边施加脉动载荷，分析脉动载荷下气泡运动情况。

1 两相流数值模型

1.1 模型的选取

在控制方程的基础上通过VOF方法来描述毛细管内两相界面的位置和运动变化，并添加CSF模型，以充分考虑表面张力的作用。模拟结果可较好地反应毛细管工作时出现的气泡聚并及破裂过程，以及各种复杂汽、液两相流型和流型间的相互转化过程。

CSF模型可表示为：

(1)

其中，σ为表面张力系数；K为界面曲率。K可表示为：

(2)

(3)

1.2 影响气泡流动的无量纲数

本文所讨论的毛细管内气泡及柱塞运动情况，受到各种物性参数的影响，需要考虑粘性力、惯性力和表面张力的影响，并且模拟环境以非重力场为主，本文采用雷诺数、Eo数、韦伯数和毛细数来分析气泡及液塞在毛细管内的运动变形。

雷诺数是表征流体粘性力与惯性力在流体流动中起到影响效果，表示为：

(4)

Eotvos数是表征流体惯性力与表面张力在流体流动中起到的影响效果，表示为：

(5)

韦伯数是表征流体惯性力与表面张力在流体流动中起到的影响效果，表示为：

(6)

毛细数作为气泡断裂界面断裂的一种判断参数，用于表征流体粘性力、表面张力以及惯性力三力共同作用下载流体流动中起到的影响效果。表示为：

(7)

2 数值模拟与分析

2.1 基于VOF模型和CSF模型的两个模拟验证

几何模型为管径1mm的毛细圆管，由于圆管为轴对称结构，所以选用二维模型进行二维模拟即可说明问题。选用VOF模型，结合CSF连续表面张力模型。选用速度入口边界条件以及压力出口边界条件，且出口压力值设定为0Pa。壁面为无滑移边界条件。气液界面采用几何重构方案进行处理，选用PISO速度压力耦合求解法进行计算，压力项采用bodyforceweighted格式，动量方程采用二阶迎风格式。经过对比发现，网格尺寸为 0.05mm的正方形网格能够很好地反映界面形态，所以选择0.05mm的网格尺寸进行计算。

为验证模拟的准确性，首先，对于大空间液池中大气泡的浮升运动进的描述行了数值模拟(模拟环境：Eo=6.99)，并与文献[6]实验结果和文献[7]LBM模拟结果进行了对比，如图1所示为不同时刻的气泡运动形态图，图中，(a)为本文VOF模拟结果，(b)为实验结果，(c)为LBM模拟结果。如图所示，气泡在液体中由于浮力的影响上浮，由最开始的球状转变为冠状，之后转变成伞状。经过对比可知，模拟结果和实验结果相吻合，本文所建立模型可以精确大空间液池中大气泡上浮时的运动变形情况。

图1 不同时刻的气泡运动形态

区别于大空间液池，毛细管由于特征长度小，管壁对气泡运动的影响不可忽略。针对毛细管内不同气泡直径下管内的运动变形进行数值模拟,并与文献[8]实验结果和文献[9]LBM模拟结果进行了对比。如图2所示，为毛细管内气泡稳定的运动形态对比，图中，(a)为本文VOF模拟结果，(b)为实验结果，(c)为LBM模拟结果。在狭窄空间内，气泡会变成弹状或者柱塞状。经过对比可知，本文模型可以精确地描述毛细管内气泡的运动变形情况。

图2 毛细管内气泡稳定运动形态对比

2.2 毛细管内气泡破裂变形分析

2.2.1 毛细管内小气泡运动变形分析

如图3所示，当We=0.05，Ca=0.002，Re=100时，气泡直径为0.3mm时，各时刻毛细管内气泡外形示意图。初始时刻为圆形的气泡向出口端运动，气泡外形由最初的圆形逐渐变形为0.002 5s时的类似半圆形，再到0.00 5s时的子弹头形，到0.01s时变成了箭头形，之后，随着流动的进行，气泡逐渐拉长到0.04s时的长箭头形并保持稳定的外形，不再发生明显的变化。根据气泡在毛细管内的位置可见，气泡在管内的运动基本为匀速运动。

图3 各时刻气泡运动变形图

为研究直径为0.3mm的气泡在管内的运动变形，保持Re=100和气泡直径不变，对不同We数及Ca数的情况下气泡在管内的运动变形进行了模拟，并得出了各时刻气泡的外形，如图4所示。

图4 各时刻气泡外形图:气泡直径0.3 mm

分别列出了气泡在0.01、0.03和0.055s时的界面形状。随着We的增大，气泡外形变化越来越剧烈，由We=0.005时的椭圆变到We=0.01时子弹头形，再到We=0.05 时的箭头形，最后气泡在We=0.1时出现气泡分裂，此时气泡的界面张力已经小到不能克服粘性力和惯性力的影响，气泡开始破裂成尺寸更加微小的气泡。

根据图示，可以划分两个气泡变形的临界值：当Ca>0.004时，此时气泡的表面张力太小，气泡会破裂成尺寸更小的气泡；当0.000 4

2.2.2 毛细管内大气泡运动破裂分析

毛细管由于管径尺度较小，所以管内经常出现大气泡的运动，所以针对管内大气泡的运动变形非常具有现实意义。为研究管内大气泡的运动变形，取气泡直径为0.6mm，Re=100，对不同We数及Ca数的情况下气泡在管内的运动变形进行了模拟，并得出了各时刻气泡的外形，如图5所示。由于气泡尺度的增加，在We=0.01和We=0.02的时候，气泡逐渐变形为类似子弹头形状的块状，形成弹状流，当气泡管径继续增加，管内气相率增大，管内将出现柱塞流。

图5 各时刻气泡外形图：气泡直径0.6 mm

跟直径为0.3mm的小气泡类似，毛细管内大气泡的运动变形也具有气泡破裂临界值，但是由于气相率的增加，管内呈现出了不一样的流型。当Ca>0.000 66时，此时气泡的表面张力太小，管内形成的气弹会破裂成尺寸更小的气泡，形成泡状流；当Ca<0.000 66时，气泡会随着We的增加，变形越来越剧烈，但不会出现气泡破裂现象。可见， 在Re=100，气泡直径为0.6mm时，Ca=6.6×10-4是气泡破裂临界值。

2.2.3 气泡直径对毛细管内气泡破裂临界值的影响

为研究气泡直径对毛细管内气泡运动变形的影响，在毛细管内，保持Ca=2×10-3不变，分别对气泡直径为0.3、0.4、0.5和0.6mm的情况进行了模型，得出了各时刻气泡的外形，如图6所示。在Ca数不变的情况下，随着气泡直径的增加，气泡由最开始不破裂，到破裂程度逐渐递增。初始直径为0.3mm直径的气泡，运动到稳定状态不会发生破裂，外形类似为箭头形；而初始直径0.6mm的大气泡，运动初期就会发生破裂，随着流动的进行，破裂成多个小气泡，形成了泡状流。

图6 不同直径气泡各时刻外形图：Ca=2×10-3

对不同的气泡直径在管内的运动变形进行了模拟计算，得出不同气泡直径时，气泡的破裂临界值Ca，如图7所示。随着气泡直径的增大，气泡的破裂临界值Ca逐渐减小，可见，在毛细管中，大气泡更难保持稳定的形态，容易破裂成小气泡。

图7 不同直径气泡的破裂临界值曲线图圆管内层流流动的速度分布规律

(8)

在管中央，层流的切应力最小，切应力的大小与偏离管子中心的距离是线性上升的关系。可见，气泡越大，其所在范围的切应力也就越大，促使气泡更容易分裂成小气泡。所以，大气泡更容易分裂的原因是由于管壁处的大粘性力导致的。

2.3 毛细管内气泡脉动运动分析

2.3.1 脉动压力载荷及UDF编写

在毛细管两端施加脉动压力载荷，创造出热管中的脉动压力驱动流的环境，分析毛细管中的气泡在脉动压力载荷下的运动变形情况。为简化模型，对毛细管两端施加正弦脉动压力载荷，两端脉动压力波的相位角为π，具体如公式所示：

(9)

(10)

其中，A为振幅，T为脉动周期，t为流动时间。所以施加脉动压力载荷后，毛细管两端压差为：

(11)

在FLUENT求解器边界条件设置中，只能设置边界压力载荷为定常值，而不能设定随时间变化的脉动值。所以，为施加脉动压力载荷，需借助用户自定义函数模块UDF。编辑的UDF代码，可以用来作为解释函数，也可以用来作为编译函数。

使用UDF模块，定制边界条件，能够达到在毛细管两端施加脉动压力载荷的要求。运用C语言编写程序包，编写毛细管入口处脉动压力载荷UDF函数，其中设置脉动频率为31.4，即脉动周期为0.2s，压力值的单位选用Pa。

2.3.2 模拟结果分析

如图8所示，气泡在毛细管内的运动到0.25T(0.25个脉动周期)时，气泡变形程度达到最大值，大致为指向出口端的箭头方向；此后，管子两端压差开始下降，管内流体速度开始减小，气泡外形又逐渐恢复到圆形，在0.5T时，气泡速度减小为0，受压力差的影响开始往反方向运动；外形变化也开始呈现倒置现象，变为指向入口端口端的箭头方向，在到达0.75T时，气泡外形与0.25T时的变形程度一样；到单个周期结束，气泡回到初始位置，速度变为零，外形基本为圆形。如此，在脉动压力驱动下做往复运动。

图8 一个脉动周期内气泡在毛细管内的运动变形图

图9展示了在不同脉动压力振幅下，气泡在管内的瞬时位置及外形图。由图可知，随着脉动压力波强度的加大，气泡在管内的变形程度加大，逐渐向弹状转变。

图9 不同脉动压力振幅下气泡在管内的瞬时位置及外形图

图10为不同脉动压力振幅下，气泡中心的速度变化示意图。由图可见，气泡在管内的速度也呈现为正弦变化，且随着脉动压力振幅的加大，气泡速度的脉动也更加剧烈。

图10 不同脉动压力振幅下气泡中心的速度变化曲线图

图11为气泡脉动运动的振幅随脉动压力振幅的变化曲线，由图可知，脉动压力强度的增加直接导致气泡脉动运动更加剧烈。

图11 气泡运动振幅随脉动压力振幅的变化曲线

可见，随着脉动压力波的增加，气泡在管内的脉动运动更加剧烈。虽然在一定压力值范围内，管内流体依然保持为层流流动，但是流动的紊乱程度更大，流体流动的稳定性减小，更易导致相变及传热的加剧；其次,两相流动与单向流动相比，由于界面的存在，其流动特性表现出极大的区别，两相之间的相间力会导致强烈的界面扰动,甚至造成复杂多变的流型变化，由此而形成两相流动本身固有的波动性,形成界面波。

毛细管正是在这种热传导产生的压力脉冲作用下，自激式地产生流体流动脉动和界面波的加剧，导致其表现出相变的加剧和传热的优良性。

3 结 语

运用VOF模型对毛细管内气泡的运动变形进行了研究分析，结果表明：

(1) 在毛细管内，气泡受粘性力、惯性力和表面张力作用，表现出复杂多变的运动形态，主要经历圆形、弹状、箭头状，直至破裂成气泡流的一个过程。

(2) 在毛细管内，气泡均存在一个破裂临界毛细值，在Re=100,1mm管内气泡直径为0.6mm的气泡破裂临界毛细值Ca=6.6×10-4，而气泡直径为0.3mm的气泡破裂临界毛细值Ca=4×10-3。在毛细管内，随着气泡直径的增加，其破裂临界毛细值增加，这是由于在毛细管内流体切应力梯度比较大造成。

(3) 在施加脉动压力载荷的毛细管内，随着脉动压力波的增加，气泡在管内的脉动运动更加剧烈。这种压力脉动产生的流体流动脉动和界面波的加剧增加了流动的紊乱程度，更易导致相变的产生及传热的加剧。气泡的外形也是随压力的脉动变化发生周期性变形。

[1]HetsroniG,GurevichM,MosyakA,etal.Boilingincapillarytubes[J].InternationalJournalofMultiphaseFlow,2003,29:1551-1563.

[2]KhandekarS,CharoensawanP,GrollM,etal.ClosedlooppulsatingheatpipesPartB:visualizationandsemi-empiricalmodeling[J].AppliedThermalEngineering,2003,23(16):2021-2033.

[3] 任桂香.T型微米级通道内微液滴形成过程的数值模拟[D].重庆:重庆大学,2009.

[4]DehaoXu.Influenceofprocessvariablesonthehydrodynamicsandperformanceofasinglelooppulsatingheatpipe[J].InternationalCommunicationsinHeatandMassTransfer,2012,39:504-508.

[5]SahaN,DasPK,SharmaPK.Influenceofprocessvariablesonthehydrodynamicsandperformanceofasinglelooppulsatingheatpipe[J].InternationalJournalofHeatandMassTransfer,2014,74(7):238-250.

[6] 武博,郝宗瑞,陈涛,等.水下气泡运动的数值模拟[J].中国科技论文在线,2010,5(8):647-650.

[7] 马强,张莹,曾建邦,等.毛细管内多相脉动流动的格子-Boltzmann模型[J].化工学报,2011,11:3101-3018.

[8]LiJ,BulusuV,GuptaNR.Bouyancy-drivenmotionofbubblesinsquarechannels[J].ChemicalEngineeringScience,2008,63(14):3766-3774.

[9]AmayaBowerL,LeeT.NumericalsimulationofsinglebubblerisinginverticalandinclinedsquarechannelusinglatticeBoltzmannmethod[J].ChemicalEngineeringScience,2011,66(5):935-952.

SIMULATION OF BUBBLE DEFORMATION AND PULSATINGMOVEMENT IN CAPILLARY

Wang Zhiqiang Liu Zhihong Zhang Ying*Xie Bowei

(SchoolofMechanicalandElectricalEngineering,NanchangUniversity,Nanchang330031,Jiangxi,China)

The movement and deformation of the bubbles in capillary are analyzed by VOF model combined with CSF model. Taking viscous force, inertial force and surface tension effect into consideration, we proceed numerical calculations of the deformation of different bubbles’ diameters in capillary and achieve the result that bubbles break into a bubble flow through the form of circle, bullet shape and arrow shape. Then, we proceed a numerical simulation when the bubbles in capillary burst while moving. The result shows that there all exist a critical value of capillary rupture. Finally, we apply a pulsating load on both sides of the capillary and analyze the situation of bubble movement under pulsating load. We achieve that pulsating movement of bubbles in capillary comes to be more vigorous with the increasing pulsating pressure wave, and the bubble’s shape also changes periodically with the pulsating of pressure.

Capillary Manifold evolution Bubble breakup Pulsating movement

2016-03-21。国家自然科学

11562011)；江西省自然科学

基金项目(20151BAB202002)。王志强，硕士生，主研领域：动力工程及工程热物理。刘志宏，本科。张莹，教授。谢博伟，硕士生。

TP31

A

10.3969/j.issn.1000-386x.2017.03.009