基于非平衡热力学模型的管内沸腾过程模拟

冯留海， 聂永广， 王江云， 雷晓东， 王　娟， 毛　羽

(1.中国石油大学 重质油国家重点实验室， 北京 102249； 2.北京低碳清洁能源研究所， 北京 102209；3.新奥科技发展有限公司， 河北 廊坊 065001； 4.中机康元粮油装备(北京)有限公司， 北京 100083)

基于非平衡热力学模型的管内沸腾过程模拟

冯留海1,2， 聂永广1,3， 王江云1， 雷晓东4， 王娟1， 毛羽1

(1.中国石油大学 重质油国家重点实验室， 北京 102249； 2.北京低碳清洁能源研究所， 北京 102209；3.新奥科技发展有限公司， 河北 廊坊 065001； 4.中机康元粮油装备(北京)有限公司， 北京 100083)

摘要：为了更准确地模拟壁面沸腾过程，根据Lavievile非平衡热力学模型对沸腾流动的壁面热流率进行了划分，通过调节函数结合Sato模型修正了沸腾流动中的混合k-ε模型，建立了同时适用于管内过冷沸腾与整体沸腾流动阶段的曳力模型。使用上述模型对环状竖直管内沸腾流动过程进行数值模拟，计算结果与实验值吻合，可以用来模拟整个管内沸腾流动过程。数值模拟显示，沸腾流动可以明显增大管壁的传热系数，但当管壁处蒸气体积分数超过0.6时，会出现传热恶化现象。

关键词：沸腾流动； 非平衡热力学； 曳力模型； 传热系数； 数值模拟

管内沸腾流动是石油、化工、能源动力等工业领域的一种常见现象，例如石油化工行业常用的管式加热炉和管壳式换热器中的管内沸腾流动。沸腾流动可以显著提高管内传热效率，所以受到人们的广泛关注。早在20世纪60年代，Bartolomei等[1]就实验研究了空泡分率对垂直管内高压水沸腾过程的影响。随着测量技术的发展，Bartel[2]、Lee等[3]对各种工况下的管内沸腾过程进行了更加详细的测量。一方面希望实验获得管内沸腾过程中详细的两相流局部参数特性，另一方面希望通过实验测量结合理论方法归纳出气泡脱离直径、成核频率和气泡停留时间等管内沸腾流动中的关键参数，从而完善数值模型，使其可以准确模拟管内沸腾过程。Kurul等[4]最早基于双流体模型对壁面热流率进行了划分，成功建立了RPI模型。因其考虑了壁面热流率的分配情况，所以能更真实地反映管内高压水过冷沸腾过程，并在相关领域得到广泛采用。此后，大量学者不断完善RPI的子模型。Anglart等[5]采用RPI模型研究了核燃料分配器中的过冷沸腾流动过程，Kocar等[6-7]建立了低压工况下水的过冷沸腾模型，李祥东等[8-9]详细总结了面向液氮过冷沸腾流动的数值模型，窦从从等[10-11]通过数值模拟确定了垂直管内高压高过冷水过冷沸腾流动的子模型。研究表明，RPI模型能够较好地预测过冷沸腾流动特性，但是由于其假设蒸气均为饱和温度，从而影响了管内沸腾流动过程中整体沸腾阶段的模拟精度。笔者在上述基础上，并根据Lavievile的非平衡热力学方法重新划分了壁面热流率[12]，通过数值计算修正了相关湍流模型和相间动量传输模型，最后将相关模型通过UDF植入通用CFD计算软件Fluent，以期更好地模拟管内沸腾流动过程。

1管内沸腾过程的数学模型

1.1双流体模型

双流体模型假设气、液两相均为连续介质，根据牛顿力学法则建立各自的控制方程，两相之间的相互作用通过控制方程中的相间传输相描述。在管内沸腾模型中，将液相视为连续相，蒸气相视为拟连续相。双流体模型是整个管内沸腾数值模型的骨架，管内沸腾流动数值模拟中所有子模型旨在封闭双流体模型的控制方程组。适用于管内沸腾流动模拟的双流体模型的基本控制方程组为式(1)～(3)。

质量守恒方程：

(1)

动量守恒方程：

(2)

能量守恒方程：

(3)

式中所有与相变相关的项都需要采用沸腾模型进行计算。

1.2壁面沸腾模型

管内沸腾是由加热壁面加热液体产生的沸腾现象，包括RPI模型在内的适用于这一过程的模型总称为壁面沸腾模型。在这一模型中，管内的传质传热主要分成2部分，即壁面处的传质传热和主流区域的传质传热，其中计算壁面区域传质和传热的关键在于合理划分壁面热流率。

(1)壁面热流率划分

前人数值模拟研究采用的RPI模型大多将壁面热流率划分为3部分，即液相对流热流率Ql、激冷热流率Qq和蒸发热流率Qe。需要说明的是，Qq是由于加热产生的气泡脱离壁面时补充的冷液体被瞬间加热升温所产生的热流率。具体的表达式为式(4)～(6)。

Ql=hl(Tw-Tl)(1-Aq)

(4)

(5)

(6)

现有的RPI模型假设流场中的蒸气均处于饱和温度状态，不考虑蒸气在壁面处的对流换热。本研究采用的Lavievile壁面热流率划分方法是在RPI模型的基础上加入了蒸气对流热流率Qv，其表达式为式(7)。壁面热流率的划分和分配函数f(αl) 如式(8)和(9)所示。式(6)中的气泡脱离直径dbw使用Unal模型计算[13]，而主流区的dbw根据过冷度[8,10-11]确定。

Qv=hv(Tw-Tv)

(7)

Qw=f(αl)(Ql+Qq+Qe)+(1-f(αl))Qv

(8)

(9)

(2)相间传质和传热

相间质量传递分为2部分。一是壁面附近液体在加热作用下沸腾相变为蒸气，二是蒸气在过冷度较高的主流区域冷凝为液体。蒸气相的传质项表达式为式(10)，液相与蒸气相的相间传热项Qlv表达式为式(11)。

(10)

Qlv=hlv(Tl-Tv)Alv

(11)

式(10)中，hlv根据Ranz-Marshall模型计算。

与RPI模型不同，因为考虑蒸气处于非饱和温度状态，因此既存在液相向气相的热量传递，也存在气相向液相的热量传递。在壁面附近，壁面热流率中的对流热流率Ql和激冷热流率Qq最终加入到液相能量源项，蒸发热流率Qe和蒸气对流热流率Qv加入到气相能量源项。

(3)湍流模型修正

RPI模型模拟过冷沸腾过程时，通常假设液相为湍流流动，气相为层流流动。但随着沸腾流动过程的发展，当管内蒸气体积分数超过0.3时，这一假设明显欠合理。因此，笔者使用混合k-ε模型计算两相的混合湍流动能km和混合湍流耗散率εm，然后采用调节函数fα计算最终蒸气湍流黏度。调节函数fα的表达式为式(12)。

(12)

式(12)中，a决定了函数曲线的变化幅度，b决定了函数曲线的阈值。图1为调整函数fα曲线。从图1可以看出，取a=20、b=0.3比较合适。最终有效黏度表达式为式(13)。

(13)

考虑到气泡扰动对液相湍流的影响，使用Sato模型[14]修正液相有效湍流黏度，表达式为式(14)～(15)。

(14)

(15)

(4)相间曳力模型修正

在对管内过冷沸腾流动数值模拟中，Ishii模型[15]表现出了良好的适用性。Ishii模型原理上仅适用于气相体积分数小于0.5的沸腾流动过程，而实际过程中随着蒸气量的增加，管内气液两相流的流型则从泡状流逐渐过渡到雾状流。所以，有必要建立同时适用于泡状流动和雾状流动的曳力模型，以研究整个管内沸腾流动过程。采用式(12)中的调节函数，将系数a、b的值分别取为20和0.5，修正的曳力系数的表达式如式(16)所示。

图1　调整函数fα曲线

K=fα·Kbubble+(1-fα)·Kdroplet

(16)

式(16)中的Kbubble采用Ishii模型计算，而Kdroplet则采用适用性较好的Schiller模型计算。

2加热管几何模型及模拟工况

加热管数值计算的几何模型如图2所示。模拟对象为Bartel的管内沸腾过程[2]，其内层为加热棒，外层为玻璃管。加热棒直径为9.55 mm，加热部分的长度为1500 mm(浅灰色部分)，玻璃管内径为19.05 mm。操作工况为常压、入口温度367.9 K，工质为水，质量流率为700 kg/(m2·s)，壁面热流率为188 kW/m2，出口设置为压力出口。因为计算域为轴对称结构，计算在柱坐标系中进行，网格划分采用四边形网格。压力与速度耦合采用PCSIMPLE(Phase couple SIMPLE)算法，对流项采用QUICK差分格式。

图2　加热管几何结构及尺寸

3结果与讨论

图3为Sato模型中不同Cμb取值时加热管z=1500 mm处液相湍流黏度沿径向的分布。由图3看到，Cμb决定了在湍流效应中是液体剪切造成的湍流还是气泡诱导产生的湍流起主导地位。当Cμb取值较小时，液体剪切造成的湍流占主导地位，这时液体湍流黏度在加热壁面与管壁的中间位置最高，而在管壁附近则较低；壁面沸腾时管壁附近加热产生的气泡会不断脱离壁面，破坏边界层，增加湍流扰动，促使壁面附近的气、液两相快速混合，因此液相湍流黏度应该在加热壁面附近较高。而Cμb最适宜的取值需要进一步分析其对管内温度场和蒸气体积分数分布造成的影响而定。

图3　不同Cμb取值时加热管内z=1500 mm处液

图4为不同Cμb时 加热管z=1500 mm处的液体温度的径向分布。从图4可以看出，Cμb取值越大，表明壁面附近液体的温度越低，说明气泡诱导的湍流作用增强了壁面附近的热传递，从而有更多的热量被带到主流区域。在主流区域，由于气相体积分数的减小，气泡产生的扰动作用减小，所以液体剪切作用产生的湍流重新占据主导地位。液体温度会直接影响加热壁面处气泡脱离直径的大小，进而影响到蒸发热流率的大小，因此液相湍流修正对准确预测加热管内蒸气体积分数起到了非常重要的作用。

图5为不同Cμb时加热管z=1500 mm处的蒸气体积分数沿径向的分布。从图5可以看出，Cμb为0.6时，计算得到的径向蒸气体积分数分布与实验值最为吻合，说明Cμb=0.6时计算得到的气泡诱导湍流作用最能反映实际情况。此外，加热管内蒸气体积分数最高处并非处于管壁上，而是在临近加热壁面的位置。这是因为在壁面润滑力的作用下，加热壁面附近生成的气泡会被推离壁面，最终在壁面附近形成蒸气富集区。

图4　不同Cμb取值时加热管内z=1500 mm

图5　不同Cμb时加热管内z=1500 mm处蒸气体积分数

图6为使用修正后的曳力模型和Ishii模型计算得到的加热管内蒸气体积分数轴向分布的对比。从图6可以看出，在过冷沸腾段，使用修正后的曳力模型计算得到的体积分数轴向分布与Ishii模型的计算结果基本吻合，而在蒸气体积分数超过0.4后(过冷沸腾向整体沸腾过渡)，修正后的曳力模型计算得到的体积分数高于Ishii模型的计算结果。修正后的曳力模型在管内沸腾的高体积分数区域与实验值吻合更好。

图6　不同曳力模型模拟的加热管内蒸气体积分数

图7是加热管壁面传热系数与壁面附近蒸气体积分数轴向分布。从图7可以看出，在加热段开始阶段，壁面传热系数随着蒸气体积分数的增加迅速增加；在蒸气体积分数略大于0.1时，壁面传热系数达到峰值，之后存在一个迅速下降阶段和一个稳定发展阶段，随后开始呈现线性递减趋势。分析认为，在管内沸腾的初始阶段，壁面产生的气泡破坏了边界层，增强了壁面附近的流体湍流作用，导致壁面热传递系数迅速升高；当壁面沸腾进一步发展时，气泡产生的湍流也增强了气泡和液体之间的热传递和质量传递，而这时壁面附近液体的过冷度依然较高，使得壁面上气泡的产生和壁面附近气泡的冷凝达到短暂的平衡状态，壁面附近蒸气体积分数增加趋势放缓，同时这种平衡状态也弱化了壁面附近的对流热传递，使得壁面传热系数减小；管内沸腾流动继续发展，管壁附近液体被冷凝的气泡充分加热使得液体过冷度明显降低，这时壁面附近的气泡将不会被液体完全冷凝，加热管壁附近的蒸气体积分数开始持续增加，而加热管壁附近不断增加的蒸气会使得壁面传热系数降低。从图7可以推断，当加热管壁面附近的蒸气体积分数超过0.6时，壁面传热系数甚至会低于单相液体加热时的壁面传热系数。

图7　加热管壁面传热系数和壁面附近

4结论

(1)根据Lavievile的非平衡热力学方法对加热管壁面沸腾中的壁面热流率进行划分，使其可以适用于整个管内沸腾过程的数值计算。同时，对壁面沸腾中的两相湍流模型修正，建立了同时适用于过冷沸腾和整体沸腾的双流体相间曳力模型。数值模拟显示，在使用混合k-ε模型时，Sato模型的修正系数Cμb取0.6时，计算结果与实验值吻合较好。

(2)在竖直管内沸腾中，壁面润滑力的作用使得蒸气体积分数最大值出现在临近加热壁面的位置。沸腾现象对管壁传热系数有明显的影响。在壁面处的蒸气体积分数介于0.1～0.3之间时，壁面的传热系数达到最大；如果壁面处蒸气体积分数继续增大，壁面传热系数会随之减小；当壁面处蒸气体积分数超过0.6时，壁面传热系数有可能会低于单相对流传热的传热系数。

符号说明：

Aq——气泡占据的壁面面积，m2；

Alv——相间传热面积，m2；

a,b——系数；

cp——液相比容，W/(m·K)；

Cμb——常系数；

db——气泡直径，m；

dbw——气泡脱离直径，m；

f——气泡脱离频率，Hz；

fα——调节函数；

Flv——两相间相互作用力，包括曳力、升力、湍流耗散力和壁面润滑力，N；

h——焓，J；

hfv——蒸发潜热，J/kg；

hl——液相对流传热系数，W/(m2·K)；

hlv——相间传热系数，W/(m2·K)；

hv——液相传热系数，W/(m2·K)；

hw——壁面传热系数，W/(m2·K)；

k——导热系数，W/(m·K)；

km——混合湍流强度，m2/s2；

K——曳力系数；

Kbubble——气泡流曳力交换系数；

Kdroplet——雾状流曳力交换系数；

mlv——液相向气相传质，代表液相的沸腾和蒸气相的冷凝，kg；

Na——气泡成核密度，m-2；

p——压力，Pa；

Qe——蒸发热流率，W/m2；

Ql——对流传热热流率，W/m2；

Qlv——相间传热速率，是非相变造成的相间能量传递，W/m3；

Qq——激冷传热热流率，W/m2；

Qw——壁面总热流率，W/m2；

Qv——蒸气单相对流传热热流率，W/m2；

r——径向位置，m；

R0——加热管半径，m；

R1——玻璃管内径，m；

t——时间，s；

T——温度，K；

Tsat——饱和温度，K;

u——速度，m/s；

z——轴向，m；

α——体积分数，%；

εm——混合湍流耗散率，m2/s2；

μ——黏性系数，Pa·s；

μeff——有效湍流黏性系数，Pa·s；

μv——气相黏性系数，Pa·s；

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

下标:

f，w——界面，壁面；

l，v——液相，蒸气相。

参考文献

[1] BARTOLOMEI G, CHANTURIYA V. Experimental study of true void fraction when boiling subcooled water in vertical tubes[J].Thermal Engineering, 1967, 14(2):123-128.

[2] BARTEL M D. Experimental investigation of subcooled boiling[D].West Lafayette: Purdue University, 1999.

[3] LEE T H, PARK G C, LEE D J. Local flow characteristics of subcooled boiling flow of water in a vertical concentric annulus[J].International Journal of Multiphase Flow, 2002, 28(8):1351-1368.

[4] KURUL N, PODOWSKI M Z. Multidimensional effects in forced convection subcooled boiling[C]//Jerusalem： Proceedings of the Ninth International Heat Transfer Conference, 1990:19-24.

[5] ANGLART H, NYLUND O, KURUL N, et al. CFD prediction of flow and phase distribution in fuel assemblies with spacers[J].Nuclear Engineering and Design, 1997, 177(1):215-228.

[6] KONCAR B, MAVKO B. CFD simulation of subcooled flow boiling at low pressure[C]//Ljubljana：Proceedings of International Conference Nuclear Energy in Central Europe, 2001:10-13.

[7] KONCAR B, KLJENAK I, MAVKO B. Modelling of local two-phase flow parameters in upward subcooled flow boiling at low pressure[J].International Journal of Heat and Mass Transfer, 2004, 47(6):1499-1513.

[8] 李祥东， 汪荣顺， 黄荣国， 等. 垂直圆管内液氮流动沸腾的理论模型及数值模拟[J].化工学报， 2006, 57(3):491-497.(LI Xiangdong, WANG Rongshun, HUANG Rongguo, et al. Modelling and numerical simulation of boiling flow of liquid nitrogen in vertical tube[J].Journal of Chemical Industry and Engineering (China), 2006, 57(3):491-497.)

[9] LI Xiangdong, WEI Wei, WANG Rongshun, et al. Numerical and experimental investigation of heat transfer on heating surface during subcooled boiling flow of liquid nitrogen[J].International Journal of Heat and Mass Transfer, 2009, 52(5):1510-1516.

[10] 窦从从， 毛羽， 王娟， 等. 高压高过冷度下过冷流动沸腾数值模拟[J].化工学报， 2010, 61(3):580-586.(DOU Congcong, MAO Yu, WANG Juan, et al. Numerical simulation of subcooled boiling flow under high pressure and high subcooling condition[J].Journal of Chemical Industry and Engineering (China), 2010, 61(3):580-586.)

[11] 窦从从， 毛羽， 王娟， 等. 高过冷度下直管内气液两相传热特性[J].化工学报， 2010, 61(9):2314-2319.(DOU Congcong, MAO Yu, WANG Juan, et al. Gas-liquid two phase flow heat transfer at high subcooling in vertical tube[J].Journal of Chemical Industry and Engineering (China), 2010, 61(9):2314-2319.)

[12] MIMOUNI S, BOUCKER M, LAVIEVILLE J, et al. Modelling and computation of cavitation and boiling bubbly flows with the NEPTUNE_CFD code[J].Nuclear Engineering and Design, 2008, 238(3):680-692.

[13] UNAL H C. Maximum bubble diameter, maximum bubble growth time and bubble growth rate[J].International Journal of Heat Mass Transfer, 1976, 19(6):643-649.

[14] SATO Y, SADATOMI M, SEKOGUCHI K. Momentum and heat transfer in two-phase bubble flow—I Theory[J].International Journal of Multiphase Flow, 1981, 7(2):167-177.

[15] ISHII M, ZUBER N. Drag coefficient and relative velocity in bubbly, droplet or particulate flows[J].AIChE Journal, 1979, 25(5):843-855.

Simulation of Boiling Flow in Vertical Tube Based on Non-Equilibrium Thermal Model

FENG Liuhai1,2, NIE Yongguang1,3, WANG Jiangyun1, LEI Xiaodong4, WANG Juan1, MAO Yu1

(1.StateKeyLaboratoryofHeavyOilProcessing,ChinaUniversityofPetroleum,Beijing102249,China;2.NationalInstituteofClean-and-LowCarbonEnergy,Beijing102209,China;3.XinAoScienceandTechnologyDevelopmentCo.Ltd.,Langfang065001,China;4.ChinaMachineryKangyuanCerealsandOilEquipment(Beijing)Co.,Ltd.,Beijing100083,China)

Abstract：To simulate the boiling process, the wall heat flux was divided based on Lavievile’s thermal non-equilibrium heat flux partition algorithm. The mixture k-ε turbulent model was improved by combination of Sato model and a modifier function, by which a new drag model was presented to simulate the full boiling flow in tube. The simulated result of boiling flow in vertical circular tube agreed with experiment data well, indicating that the model was valid in simulation of boiling flow in whole tube. It is found from the simulation that the wall heat transfer coefficient was evidently enhanced in boiling, however, heat transfer deterioration appeared when the volume fraction of vapor near the wall reached about 0.6.

Key words：boiling flow; non-equilibrium thermal dynamics; drag force model; heat transfer coefficient; numerical simulation

收稿日期：2015-04-20

基金项目：国家重点基础研究发展计划“973”项目(2010CB226902)和中国石油大学(北京)科研基金项目(2462015YQ0303)资助

文章编号：1001-8719(2016)03-0591-06

中图分类号：TK124

文献标识码：A

doi:10.3969/j.issn.1001-8719.2016.03.021

第一作者： 冯留海，男，博士研究生，从事多相流动的数值模拟与实验方面的研究

通讯联系人： 王娟，女，副教授，博士，从事多相流动与燃烧过程的数值模拟与实验研究；Tel：010-89733293；E-mail：wangjuan@cup.edu.cn；毛羽，男，教授，博士，从事多相流动及燃烧、气固分离及液体雾化技术、化工过程装备优化等方面的研究；Tel：010-89733293；E-mail：maoyu@cup.edu.cn