部分重力环境下水平方管内气液两相流动特性的数值研究

黄 娜，周云龙，高 聚

(1.华北电力大学能源动力与机械工程学院，北京102206;2.东北电力大学能源与动力工程学院，吉林吉林132012;3.中机国能电力工程有限公司，上海200061)

0 引言

气液两相流换热系统由于具备高效的传热能力，较小的温度差分布和对抽运功率要求低等优点，已经在航空航天和高海拔等实际应用中起到越来越重要的作用［1］。美国宇航局预计在未来启动两相流换热系统来担任主动热控制、生命支持、推进剂输送和原位资源利用等航空任务［2］。另外，新型两相流循环热控制系统将负责未来载人探月和探索火星工程中部分重力环境(partial gravity，10－2g0＜g ＜g0)下的能量输运工作［3］。同时，部分重力下的气液两相流动也是月球基地中许多问题都涉及到的关键技术问题［4］。由于部分重力环境下气液两相的流动特性对气液两相流设备的工作效率具有重要影响，因此，部分重力环境下的气液两相流动问题成为实际应用面临的一个重要挑战。

近年来，人们利用落管或落塔、驾驶飞机沿抛物线轨迹飞行、高空气球搭载、火箭搭载、空间轨道飞行器搭载以及地面模拟等方法制造失重环境，进行了一些失重环境下的气液两相流动研究。1999年8 月，中国科学院国家微重力实验室与俄罗斯Keldysh 研究中心展开合作，以“和平号”空间站为平台，首次在长期、稳定的微重力环境下进行不同重力环境下的气液两相流实验研究［5］。然而，由于受到空间实验时间的限制，在部分重力环境下仅进行了14 组实验，难以获取大量实验数据来建立确定性结论。

由于长时间、稳定的部分重力环境难以实现，且需要付出昂贵的代价，越来越多的学者们采用理论分析和数值模拟方法进行部分重力环境下的研究工作。赵建福等［6］以“和平号”空间站为背景研究了部分重力环境(0.1g0和0.014g0)下的气液两相流动特征，重点关注于流型和流型转变界限。Clarke 和Rezkallah 等［7］基于VOF 方法建立了二维非稳态模型来分析微重力环境下竖直管内泡状流的漂移速度，研究表明，液体Reynolds 数、管径及气泡的直径和径向位置对漂移速度均有重要影响。为了跟踪竖直管内的气弹运动，Gabriel［8］研发了一个二维轴对称弹状流模型，并对μg0和g0环境下弹状流的流场分布进行了数值研究。除此之外，Yang［9］开发了一个三维非稳态模型来研究微重力环境下球形玻璃容器内液体表面的运动。目前，关于水平方管的研究结果包含了微重力与常重力环境下，气液及液液的流动特性，本文通过数值模拟方法，建立水平方管内的气液两相流三维非稳态模型，对部分重力环境下空气－水和制冷剂R134a 蒸汽－液体两相的流动的压降特性进行数值研究。

1 物理模型

根据实验中应用成熟的典型设计建立了数值计算的物理模型，如图1。气体和液体由交叉部分流入并混合，随后流入水平主管道。设定液体沿z轴方向引入，气体沿y 轴方向引入(微重力环境下的两相流动不依赖于流动方向［6］)。由于当试验段长径比大于50 时，可以忽略进口效应对失重条件下气液两相流型的影响［10］，故计算过程重点关注于50D～60D 间管段的压降变化。

图1 物理模型Fig.1 Physical model

2 数值模拟

(1)控制方程

流体体积函数模型(volume of fluid，VOF)已经被广泛应用于非稳态气液两相流动计算［11］。本文基于VOF 模型，采用FLUENT 内置跟踪界面的界面重构法(piecewise linear interface construction，PLIC)计算气液交界面的平流流态［12］，通过连续表面张力模型(continuous surface force，CSF)计算表面张力对动量方程的影响［13］，应用RNG k－ε湍流模型结合标准壁面函数研究湍流的影响，k 和ε 由整个两相流域确定。研究中假定两相流体不易混溶，不可压缩，两相间没有质量输运，不考虑过程中的热交换。

VOF 模型中，每个控制体内所有相的体积分数的和为1。将气相的体积分数记为αG，液相的体积分数记为αL，则:

质量守恒方程和动量守恒方程:

每个控制体内的密度和黏度分别定义为

式中:σ是表面张力;ρ 是体积平均密度;表面曲率κ 为

计算区域如图1所示，采用均匀六面体网格划分。作者对比了边长为0.05 mm，0.1 mm，0.2 mm 和0.5 mm 四种不同尺寸网格的计算精度和计算成本。图2 的对比结果显示，与0.2 mm 和0.5 mm 网格相比，0.1 mm 网格的计算精度显著提高;与0.1 mm 网格相比，0.05 mm 网格的计算精度提高不明显，且需要耗费更多的计算时间。故选择0.1 mm 网格，网格总数约为308 000 个，时间步长为10－5s。

(2)边界条件

入口边界条件选择方向垂直于入口边界的恒定速度入口，出口边界条件是压力出口。物面边界条件是无滑移固壁边界。初始时刻管道内充满处于静止状态的气体。

(3)控制方程的离散

控制方程的对流项采用二阶迎风格式离散，压力方程的差分格式选择二阶迎风格式来离散控制方程的对流项，选择交错压力格式(Pressure Staggering Option，PRESTO)作为压力方程的差分格式，选择压力的隐式算子分割法(Pressure Implicit with Splitting of Operators，PISO)进行压力－速度耦合计算。以残差值全部小于10－5作为收敛标准。

图2 网格尺寸分别为0.05 mm，0.1 mm，0.2 mm 和0.5 mm 时液体入口处的相分布图Fig.2 Contours of phase at the entrance of liquid with different mesh sizes: 0.05 mm，0.1 mm，0.2 mm and 0.5 mm

3 结果与讨论

(1)数值方法验证

表1 实验物性参数Tab.1 Physical parameters in experiments

表2 实验参数范围Tab.2 Parameter range in experiments

根据赵建福等人［6］实验中的操作工况(实验材料物性参数列于表1，实验参数范围列于表2)进行数值模拟。通过统计不同流型出现的工况得到了流型分布图，数值计算结果同实验结果的对比如图3所示，图中字母代表的流型列于表3。图中的点为数值模拟得到的结果，作者在0.1g0重力环境下计算了30 组数据用于数值方法验证。表2中，Bond 数、Weber 数和Reynold 数均为无量纲数，计算式依次为式中ρ 为密度;a 为剩余加速度;u 为表观流速;d为特征尺度;σ 为表面张力。除与赵建福等［6］的实验结果对比外，图中还标志了Weisman 模型［14］和Lee 模型［15］预测的微重力环境下气液流型的转变界限。由对比可知，模拟结果与实验结果和模型的预测结果都吻合较好。

图3 数值计算结果与实验结果对比Fig.3 Comparisons between numerical results and experimental results

表3 字母代表的流型Tab.3 Flow patterns represented by letters

(2)部分重力环境下的压降分布

模拟选用材料的物理参数列于表4，其中空气和水的温度是20℃，压力环境是0.10 MPa;R134a蒸汽和液体的温度是饱和温度30℃，压力环境是0.77 MPa;重力场分别选取g0(常重力，g0=9.8 m/s2)，0.38g0(马丁重力)和0.17g0(月球重力)进行计算。数值计算过程中气液两相工况参数范围列于表5，其中，U 为表观速度，I 为表观动量通量(I =ρU2)，We 为Weber 数，Re 为Reynolds 数，Fr 为Froude 数。

一般地，管道内的流动压降主要由沿程摩擦压降Δpf、局部压降Δpc、沿程加速压降Δpa和重力压降Δpg四部分组成［16］，即

表4 计算物性参数Tab.4 Physical parameters in numerical simulations

表5 数值计算选取的工况范围Tab.5 Parameter range in numerical simulations

由于部分重力环境下水平管内气液两相绝热流动可忽略Δpc，Δpa和Δpg，则Δp = Δpf。

图4所示为水平方管内不同气液速度下，空气－水和R134a 蒸汽－液体两相流在部分重力和常重力环境下的压降分布。图中I 为表观动量通量，下标G 和L 依次代表气相和液相。由图可知，两种混合物的压降均随着气、液速度的增大而增大，这一现象与气液两相流压降的一般规律相符。相同工况时，部分重力下的压降大于常重力下的压降。其原因是与常重力环境相比，部分重力环境下的浮力减小，滑移比的降低导致液相流速的增大，气泡运动引起的紊动幅度减低。由于水的粘度是R134a 液体粘度的五倍，故与空气－水两相流相比，R134a 蒸汽－液体的压降较小，这种现象随着重力环境的减小而变得越明显。

(3)压降的对比及修正

本节将部分重力环境下空气－ 水两相流和R134a 蒸汽－液体两相流压降的数值结果与各模型的预测值进行对比分析。图5所示为数值结果和均相流模型［17］的对比。由图可知，模型的预测值整体偏大，在摩擦压降较小时其偏差更大。这是因为，均相流模型是在常重力环境下提出的压降预测模型，在部分重力环境下由于气液相对运动强度减弱，使得气液与壁面摩擦阻力显著减小，导致均相流模型的预测值出现了较大偏差。大约61.24 %的数据高于+30 %的偏差线，均相流模型不能准确地预测部分重力环境下气液两相压降的数值计算结果。

图4 不同重力环境下空气－水两相流、R134a 蒸汽－液体两相流的压降对比Fig.4 Comparisons of pressure drop for air-water two phase flow and R134a vapor-liquid two phase flow under different gravity

图5 模拟值与均相流模型预测值对比Fig.5 Comparison of numerical result with homogeneous flow model

数值结果和Friedel 模型［18］间的对比如图6所示。由图可知，Friedel 模型的预测值比数值结果偏大，约18.55 %的数据超出了+30 %的误差线。Friedel 模型也难以获得令人满意的预测结果。

数值结果与Chisholm 模型［19］预测值的比较结果如图7所示。相对均相流模型和Friedel 模型而言，Chisholm 模型的预测值较为理想，总体偏差分布在+30 %偏差线以内。为了得到较理想的模型预测结果，需要对现有的模型进行修正。

图6 模拟值与Friedel 模型预测值对比Fig.6 Comparison of numerical result with Friedel model

图7 模拟值与Chisholm 模型预测值对比Fig.7 Comparison of numerical result with Chisholm model

为方便应用，Chisholm 将Lockhart-Martinelli关联式简写为

式中:关联参数C 的取值取决于气、液流动状态，具体如表6所示。

表6 Tabular 常数的取值Tab.6 Values of tabular constant

图8 为数值计算结果与该关联式预测值的对比。由图可知，气、液两相流动均为层流流态时，约36.11 %的数值结果超出了关联式的预测值;气、液两相流动均为湍流流态时，约40.52 %的数值结果超出了关联式的预测值。因此，根据气液两相Reynolds 数计算C 值难以对部分重力环境下结果进行准确预测。大量文献研究中也显示，根据气液两相Reynolds 数对关联参数C 进行取值，难以全面准确地体现流动的实际情况和参数，故需要对现有模型进行修正。

图8 数值结果与Chisholm 的分区对比Fig.8 Comparison of numerical result with Chisholm correlation

按照分液相Reynolds 数对流动重新分区:层流区为ReL＜1 200;过渡区为1 200 ＜ReL＜2 100;湍流区为ReL＞2 100。由于Chisholm 关联式中两相相互作用项中X 指数并不为1［20］，故将该关联式修正为

对数值结果在层流区、过渡区和湍流区的分别进行最小二乘拟合，得到m =1.31，得到的关系式列于表7。

表7 Tabular 常数的取值Tab.7 Values of tabular constant

图9 根据液相Reynolds 数分区时ΦL2与X 的关系Fig.9 Relationship between ΦL2 and X divided according to liquid Reynolds number

图10 模拟值与Chisholm 模型预测值对比Fig.10 Comparison of numerical result with Chisholm model

图9 为根据分液相Reynolds 数对流动进行重新分区的结果。根据修正后的Chisholm 关联式，有89.7 %的数值结果分布在模型预测范围内。数值计算结果与拟合公式预测值的比较如图10所示，有93.1 %的点落在±25 %的误差带以内，与数值计算结果有较高的吻合度。

4 结论

(1)与常重力环境相比，部分重力环境下的气液两相流压降较大;与常重力相同，部分重力环境下两相流的压降随气、液入口速度的提高而增大;相同工况下，由于R134a 的粘度仅为水的1/5，故R134a 蒸汽－液体的压降明显低于空气－水两相流的压降。

(2)对于部分重力环境下的气液两相流动，均相流模型和Friedel 模型的预测值与数值计算结果存在较大偏差，Chisholm 模型的预测值偏差较小，但仍然难以得到满意的预测结果。综合分析数值计算结果，以分液相Reynolds 数为基础重新对流动分区，同时对Chisholm 关联式修正。对比显示，修正后的关联式能较好地预测部分重力环境下气液两相流的压降。

［1］Rezkallah K S.Recent progress in the studies of twophase flow at microgravity conditions［J］.Advances in Space Research，1995，16:123–132.

［2］Hurlbert K M，White L C，Best F R，et al.Scaling two-phase flows to Mars and Moon gravity conditions［J］.International Journal of Multiphase Flow，2004，30:351–368.

［3］Park C，Sunada E.Vapor compression hybrid twophase loop technology for Lunar surface applications［J］.Journal of AIP Conference Proceedings，2008，969，:37–43.

［4］赵建福.微重力条件下气/液两相流流型的研究进展［J］.力学进展，1999，29 (3):369－382.

［5］Zhao J F，Xie J C，Lin H，et al.Experimental studies on two-phase flow patterns aboard the Mir space station［J］.International Journal of Multiphase Flow，2001，27 (11):1931－1944.

［6］赵建福，解京昌，林海等.部分重力条件下气液两相流型研究［J］.工程热物理学报，2004，25(1):85－87.

［7］Clarke N N，Rezkallah K S.A study of drift velocity in bubbly two-phase flow under microgravity conditions［J］.International Journal of Multiphase Flow，2001，27:1533–1554.

［8］Gabriel K S.Microgravity two-phase flow and heat transfer［M］.Springer，Berlin，2007:.

［9］Yang A.Investigation of liquid–gas interfacial shapes in reduced gravitational environments［J］.International Journal of Mechanical Science，2008，50:1304－1315.

［10］赵建福，林海，解京昌，等.失重飞机搭载气液两相流实验研究［J］.空间科学学报，2000，20(4):340－347.

［11］Hirt C W，Nichols B D.Volume of fluid (VOF)method for the dynamics of free boundaries［J］.Computer Physics，1981，46:201－225.

［12］Morton K W，Baibnes M J.Numerical methods for fluid dynamics［M］.New York:Academic Press，1982:273－274.

［13］Brackbill J U，Kothe D B，Zemach C.A continuum method for modeling Surface tension［J］.Computer Physics，1992，100:335–354.

［14］Weisman J，Kang Y.Flow pattern transitions in vertical and upwards inclined lines［J］.International Journal of Multiphase Flow，1981，7 (3):271－291.

［15］Lee J.Scaling analysis of gas-liquid two-phase flow pattern in microgravity.In:31stAerospace Sciences Meeting ＆ Exhibit，Reno，NV，1993.

［16］周云龙，孙斌，陈飞.气液两相流型智能识别理论及方法［M］.北京:科学出版社，2007.112－117.

［17］周云龙，郭婷婷.高等流体力学［M］.北京:中国电力出版社，2008.171－173.

［18］Friedel L.Improved friction pressure drop correlations for horizontal and vertical two-phase pipe flow［C］.Ispra，Italy:European Two-phase Group Meeting，1979.

［19］Chisholm D.Two-Phase Flow in Pipelines and Heat Exchanger［M］.New York:Longman，1983.

［20］徐建鸿，骆广生，陈桂光，等.液－液微尺度混合体系的传质模型［J］.化工学报，2005，56 (3):435–440.