基于CFD的竖直圆管内空气-水逆向流动界面行为研究

万　洁，孙　皖，李　想，朱隆祥，马在勇，张卢腾，潘良明

基于CFD的竖直圆管内空气-水逆向流动界面行为研究

万洁1,2，孙皖1,*，李想1,2，朱隆祥1,2，马在勇1,2，张卢腾1,2，潘良明1,2

（1. 重庆大学低品位能源利用技术及系统教育部重点实验室，重庆 400044； 2. 重庆大学核工程与核技术系，重庆 400044）

为研究空气-水逆向流动局部特性，本文基于VOF（volume of fluid）模型对=25 mm的圆管中气液逆向流动过程中界面波的演变开展了数值计算，分析了流动过程中界面波波形的变化，计算结果与实验中液膜厚度变化符合较好。对不同液相流速工况下的流场变化分析可得，流动过程中界面波前后产生的压差引入了削弱界面波稳定性的附加作用力，且界面波前后压力变化随液相流速的增大而增大；与此同时，气相流速越大，界面波波幅沿轴向发展越迅速，表明界面波的不稳定性随气相流速的增大而增大。由此可得，气相流速及液相流速的增大都会使界面波的稳定性降低，即更容易发生液泛。

CFD；液泛；VOF；界面波；液膜厚度

反应堆发生大破口事故时，应急堆芯冷却系统（ECCS）投入运行，应急冷却水与堆芯内产生的蒸汽在堆芯孔板、下降段或热管段等位置处发生液泛导致堆芯冷却效率降低甚至无法冷却从而影响反应堆安全。对于液泛的表面波动特性、气-液相界面等的局部特性研究，即使是空气-水逆向流动工况下都有很大不确定性，而数值模拟因其能较好反映局部特性而被广泛使用。Jayanti等[1]研究发现施加在波上的剪切力由波附近的压力变化引起且与流道尺寸有关，并基于管径提出了一种区分液泛机理的方法。Shearer和Davidson[2]则提出了一种可以预测竖直表面上波形及波幅的理论。Da Riva和Del Col[3]采用VOF[4]界面重构法对一个轴对称域进行了模拟，发现环状波的形成频率与Barbosa Jr.等[5]的实验结果符合较好。Murase等[6]，Minami等[7,8]和Utanohara等[9]在1/15 的PWR模型中对热段进行了空气-水三维逆向流动模拟，对界面摩擦、流型及逆向流动极限（CCFL）进行了研究并总结得出，采用具有合适的界面摩擦关系式的两流体模型比VOF模型结果更好。两流体模型虽在对逆向流动极限（CCFL）的模拟中表现较好，但本文关注界面波的界面行为，VOF模型对于相界面的识别更为准确，因此本文采用VOF模型展开模拟研究。

本文基于VOF模型采用ANANSYS FLUENT研究了逆向流动过程中的界面波动特性，并将模拟结果与Wan等[10]的实验结果进行了对比。本研究通过跟踪界面波行为以得到液膜的演变，同时讨论了气相流速对波形变化的影响并对液泛过程中以及液泛后的流场进行了分析，得出气液逆向流动过程中相间界面的相互作用，以更好地分析CCFL发生的机理。

1 数值模拟

1.1　几何模型及边界条件

本文中数值计算几何区域为一直径25 mm，管长156 mm的竖直圆管。考虑到流动过程呈对称分布所以本研究取管中心为对称轴，以管左侧区域的流动为研究对象如图1所示，注水口长度为7.9 mm，采用常规光管注水，网格数量为109 950。管的下部为气相速度入口，左侧注水口为液相质量入口，上部出口为两相压力出口，其余管道壁面为无滑移壁面条件。

图1　液泛数值计算几何模型

本文采用VOF模型对界面进行追踪，压力梯度耦合采用PISO算法，梯度项离散采用最小二乘法单元格方法。除空泡份额以外的变量都采用二阶迎风格式进行离散，空泡份额采用几何重构法离散，瞬态计算的时间离散采用一阶隐式法。此外，当所有变量残差小于10-4时认为计算达到收敛。本文共选取了包含液泛及液泛后过程的5个工况展开计算，其中液相流速根据Wan等[10]实验中两种液相流速下对应两种不同的液泛机理分别选取了0.028 m/s及0.121 m/s，气相流速分别设置为实验中液泛起始点及液泛后对应的气相流速工况，工况表如表1所示。所有物性参数均定义为常温（25 ℃）常压（0.1 MPa）下的数值。

表1 液泛数值计算工况

1.2　网格无关性验证

为验证计算结果的网格无关性，本文采用不同的网格数量在l=0.028 m/s，g=6.379 m/s工况下进行数值计算并对比其界面波生长特性。从图2中可看出，界面波波幅随时间的变化随着网格数量的增加变化不大，本文最终选取了网格数量109 950进行计算。

图2　不同网格数量下界面波变化特性（jl=0.028 m/s，jg=6.379 m/s）

图2　不同网格数量下界面波变化特性（jl=0.028 m/s，jg=6.379 m/s）（续）

2 结果与分析

2.1　波形分析

2.1.1小液相流速

因界面波生长具有周期性，本文对某一周期内变化进行了分析。如图3所示，由于气流的拖曳作用使界面产生的界面波波幅增大，在该工况下波沿径向的增大也更为显著。=0.031～0.056 s时，波峰沿着气流方向不断被拉长到一临界点后其克服自身表面张力后被撕碎，分离出液滴随气流向上流动，此后界面波在气流作用下不断被撕碎使波幅逐渐变小，直至下一个小波产生从而进入到下一个周期的界面波生长过程。图4中液膜随时间的变化曲线与模拟结果变化一致，液膜厚度随时间的推移先增大后逐渐减小，其达到最大值时对应模拟过程中波幅达到峰值。由=0.056 s及0.059 s可看出，波峰处液体被撕裂形成小液滴被气相夹带向上导致液膜变小，这是由于界面波在管内径向上不断增大使得气相在截面上所占流通面积减小而流速增大，气液间的相间剪切增大，对应图4中液膜减小阶段。

图3　jl=0.028 m/s，jg=6.379 m/s时界面波形变化

图4　jl=0.028 m/s，jg=6.379 m/s时液膜在0.60 s内的变化[10]

2.1.2大液相流速

如图5所示，大液相流速下界面波的发展十分迅速。=0.019～0.026 s时，界面波在径向上快速拉长与对称侧的界面波几乎在管内形成搭接，波峰向气流流动方向弯曲直至=0.033 s时被气相吹散，整个过程用时0.024 s，=0.057 s时波开始继续生长并重复以上过程。由图6可看出，大液相流速下管内液膜厚度显著增大，对比图5也可看出管内截面液相份额显著增大，进而使气相在管截面所占面积相应变小使得气相流速增大，在气流作用下易形成搅混流，对应图6中液膜厚度的随机增大或减小。搅混作用使液相入口处的波难以呈孤立波形式沿轴向及径向发展。

图5　jl=0.121 m/s，jg=4.076 m/s时界面波形变化

2.2　气相流速的影响

如图7（a）所示，/0代表最大波幅所在轴向位置，其中为波峰位置在圆管竖直方向上的位置，0为最大波幅初始所在位置，这里假设界面波初始时是在液相入口处呈一个对称分布的独立波，因此将波峰最初位置定于液相入口中点处，即初始时0=68.95 mm。g=6.379 m/s时，最大波幅值随时间变化有显著增大，而最大波幅所在位置沿轴向变化很小，这说明在该气相流速下界面波几乎是在其初始产生的位置不断增大且受气流的拖曳力较小，界面波主要是沿径向发展直至达到临界波幅值时，波峰受到气流作用沿气流方向拉伸后被撕碎形成液滴进而被气流带走或重新在管壁上沉积。如图7（b）及图7（c）所示，g=7.000 m/s及7.721 m/s时，最大波幅较g=6.379 m/s工况减小且最大波幅所在轴向位置随时间不断变化，这说明随着气相流速的增大，界面波沿轴向的发展更加显著且气液界面更加平坦，波峰处被气流撕碎形成液滴时对应的临界波幅更小。而气相流速越小界面波在径向上的发展更显著，即气液界面波峰更显著，波峰被撕碎时的临界波幅更大。

图6　jl=0.0.121 m/s，jg=4.076 m/s时液膜在0.60 s内的变化[10]

2.3　流场分析

2.3.1小液相流速

如图8所示，界面波的存在使得气相流通面积变小从而导致气相流速变大，因此界面波迎风侧气相流线较密，在波峰处速度分布变得平坦，而波峰后方为尾流低压区，这与Davide Del Col等[3]模拟结果一致，如图9所示，压力在界面波背风侧突然下降，而在气相流动的下游位置压力则逐渐升高直至到达一常数，这是因为气流在远离界面波的一段距离后扰动减小而不再受尾流区影响。界面波的背风侧与迎风侧之间的压力差使界面波受到了向上的附加力作用，加速了液泛的发生。除此之外，随着时间的推移，界面波在径向上不断生长使得压差不断增大，进而导致尾流区不断增大直至波峰处被气流撕碎，大尾流造成流动分离。

图7　jl=0.028 m/s时不同气相流速下最大波幅值及其轴向位置随时间变化图

图8　jl=0.028 m/s时界面波附近流场图：（a）jg=6.379 m/s（b）jg=7.000 m/s

图9　界面波达到最大波幅时的压力变化[3]

2.3.2大液相流速

如图10所示，界面波在初始形成时背风侧也会形成尾涡，而在大液相流速下界面波在气流作用下会在径向上快速发展，使得背风侧尾涡也快速增大进而导致波前后两侧压差迅速变化，因此在大液相流速下界面波极易在管内被气流吹散形成搅混流，尾涡耗散后界面波重新生长并重复以上过程。大液相流速下由于管内形成搅混流态，管内各截面上都有一定液相份额，随着液相的不断补充以及管内扰动的存在，注水口处的波难以稳定地呈孤立波生长。

3　结论

本文基于VOF模型，对=25 mm的圆管中气液逆向流动过程中界面波的演变开展了CFD计算。得到如下结论：

（1）小液相流速下界面波达到临界波幅后，波峰在气流作用下克服自身表面张力被撕碎；大液相流速下界面波发展迅速，难以孤立波稳定生长且波幅较小，易形成搅混流。

（2）小液相流速下随着气相流速的增大，界面波发展过程中界面的最大波幅不断减小且最大波幅所在轴向位置沿气流方向向下游发展。

（3）小液相流速下界面波在受到气流作用后波幅逐渐增大，波前后产生的压差对界面波有附加力的作用，界面波在气流及附加力双重作用下变得不稳定进而促使液泛发生。大液相流速下界面波沿径向发展十分迅速，其背风侧的尾涡快速增大从而导致界面波两侧的压差也迅速增大，加大了波的不稳定性。

致谢

本工作得到了国家重点研发计划（2018YFB1900400）的支持。

［1］ Jayanti S，Tokarz A，Hewitt GF. Theoretical investigation of the diameter effect on flooding in countercurrent flow［J］. International Journal of Multiphase Flow，1996，22（2）：307-324.

［2］ Shearer C J，Davidson. The investigation of a standing wave due to gas blowing upwards over a liquid film；its relation to flooding in wetted-wall columns［J］. Journal of fluid mechanics. 1965，22（2）：321-335.

［3］ Enrico D R，Davide D C. Numerical simulation of churn flow in a vertical pipe［J］. Chemical engineering science，2009，64（17）：3753-3765.

［4］ Hirt C W，Nichols B D. Volume of fluid（VOF）method for the dynamics of free boundaries［J］. Journal of computational physics，1981，39（1）：201-225.

［5］ Barbosa J F，Lima E C，Zanão R A，et al. The use of a W-Rh permanent modifier for direct determination of bismuth in urine and whole blood by electrothermal atomic absorption spectrometry［J］. Journal of Analytical Atomic Spectrometry，2001，16（8）：842-846.

［6］ Murase M，Kinoshita I，Utanohara Y，et al. Numerical calculations for air-water tests on CCFL in different-scale models of a PWR hot leg；International Conference on Nuclear Engineering，2010.

［7］ Minami N，Murase M，Tomiyama A. Countercurrent Gas-Liquid Flow in a PWR Hot Leg under Reflux Cooling（Ⅱ）Numerical Simulation of 1/15 Scale Air-Water Tests［J］. Journal of Nuclear Science Technology，2010，47（2）：149-155.

［8］ Minami N，Nishiwaki D，Nariai T，et al. Countercurrent gas-liquid flow in a PWR hot leg under reflux cooling（Ⅰ）air-water tests for 1/15-scale model of a PWR hot leg［J］. Journal of Nuclear Science Technology，2010，47（2）：142-148.

［9］ Utanohara Y.，Kinoshita I.，Murase M.，et al. Numerical simulation using CFD software of countercurrent gas–liquid flow in a PWR hot leg under reflux condensation［J］. Nuclear Engineering Design，2011，241（5）：1643-1655.

［10］ Wan J，Sun W，Deng J，et al. Experimental study on air-water countercurrent flow limitation in a vertical tube based on measurement of film thickness behavior［J］. Nuclear Engineering Technology，2021，53（6）：1821-1833.

Study on the Interfacial Behavior of the Air/Water Counter-current Flow in the Vertical Tube Based on the CFD

WAN Jie1,2，SUN Wan1,*，LI Xiang1,2，ZHU Longxiang1,2，MA Zaiyong1,2，ZHANG Luteng1,2， PAN Liangming1,2

（1. Key Laboratory of Low-Grade Energy Utilization Technologies and Systems，Ministry of Education，Chongqing University，Chongqing 400044，China；2. Department of Nuclear Engineering and Technology，Chongqing University，Chongqing 400044，China）

In order to study the local characteristics of the air-water counter-current flow, a numerical simulation of interfacial wave in a vertical tube with 25 mm diameter has been carried out based on the VOF (volume of fluid) model. The change of interfacial wave was analyzed and showed a good agreement with the experimental film thickness. Through the analysis of the variation of the flow field under different liquid flow rates, it can be concluded that the pressure difference generated between the windward and leeward of the interfacial wave introduces an additional force, which weakens the stability of the interfacial wave. And the pressure difference increases with the increase of the liquid flow rate. Moreover, the greater the gas flow rate, the faster the amplitude of the interfacial wave develops along the axial direction, indicating that the instability of the interfacial wave increases with the increase of the gas flow rate. Consequently, the increase of the gas and liquid flow rate reduces the stability of the interfacial wave, namely, the flooding is more likely to occur.

CFD; Flooding; VOF; Interfacial wave; Liquid film thickness

TL334

A

0258-0918（2023）05-1105-08

2022-08-04

国家重点研发计划资助项目（2018YFB1900400）

万 洁（1996—），女，贵州三穗人，博士研究生，现主要从事两相逆向流动极限方面研究

孙 皖，E-mail：sunwan@cqu.edu.cn