窄矩形通道内搅混流和环状流相界面参数的计算方法

程林海，谷海峰,*，刘安泰，虞 想，张君毅

(1.哈尔滨工程大学 核安全与仿真技术国防重点学科实验室，黑龙江 哈尔滨 150001;2.中国核动力研究设计院 中核核反应堆热工水力技术重点实验室，四川 成都 610213)

两流体模型是目前两相系统流体传热学中最详细、精确的宏观模型[1]。准确获取窄矩形通道内相界面参数是构建和验证窄矩形通道中两流体模型的必要条件[2-3]，对板状燃料的核反应堆堆芯传热分析具有重要意义。空泡份额与界面浓度是至关重要的相界面参数[4]，同时由于窄矩形通道横截面具有较大的长宽比，导致通道宽边方向上局部参数具有一定的分布规律[5-7]，局部参数反映局部位置的传热传质能力，局部位置发生传热恶化往往会威胁整个设备的安全运行。因此获取窄矩形通道内空泡份额和界面浓度，以及两参数在其宽边方向上的局部分布，对设备设计、安全评估、构建和验证两流体模型十分重要。

目前在窄矩形通道内相界面参数的实验研究中，因受通道尺寸限制，测量方法以非接触式的高速摄像法为主，且研究对象主要是空泡份额较低的泡状流和弹状流工况。张偲琪[7]对横截面为40 mm×4 mm的窄矩形通道进行拍摄，在研究脉动条件下泡状流的相界面参数过程中提出了相界面参数计算方法。Chalgeri等[8]利用高速摄像仪对横截面为66.6 mm×2.35 mm的窄矩形通道进行拍摄，利用Fij软件计算了泡状流下的相界面参数。Shen等[9]针对40 mm×0.993 mm的窄矩形通道内泡状流和帽状流工况，提出了空泡份额和界面浓度计算方法。然而在窄矩形通道内高含气率下的搅混流和环状流流型中，气液相界面较为复杂，且气相空间缺少封闭且规则的几何形状，难以通过可视化实验测量相界面参数。这导致该流型下所对应的相界面参数数据较少[10]，给两流体模型的验证带来困难。

本文通过在横截面为65 mm×2 mm的竖直窄矩形通道内进行空气-水向上两相流动实验，基于高速摄像技术拍摄搅混流和环状流条件下两相流动图像，提出一种获取该流型下通道内空泡份额和界面浓度以及两参数在其通道宽边方向上局部分布的计算方法。

1 实验装置及工况

实验回路示意图如图1所示，窄矩形通道试验段材质为高透光率的聚甲基丙烯酸甲酯(PMMA)，中部流道截面为65 mm×2 mm，可视化有效长度为2 m。

通过调节气、水回路的阀门，使气、水流量稳定在目标工况，通道内出现相应的流型。调节高速摄像仪，使其拍摄到清晰且完整的两相流视频图像。拍摄速率为6 000帧/s，分辨率为1 280×960 pixel。实验中对试验段上、中、下游分别进行拍摄，拍摄中央分别距混合腔出口600、1 200、1 800 mm。

图1 实验回路示意图Fig.1 Sketch of experimental loop

为使气相和液相在入口处充分混合，试验段入口设有多孔混合腔室，如图2所示。混合腔室内均匀设有直径为1 mm、间距为1 mm的孔洞。

图2 两相混合腔结构Fig.2 Two-phase mixing cavity structure

所开展两相流动实验中，气相折算速度jg=1～9 m/s，液相折算速度jf=0.1～1.5 m/s，涵盖搅混流、环状流流型区间。具体实验工况如图3所示。

2 相界面参数获取方法

2.1 相界面特征

实验拍摄图像如图4所示。从图4可看出，高空泡份额下气相空间内具有以下特征：1) 气相空间具有连续性，无论是搅混流(图4a)还是环状流(图4c)，图像中央位置均存在较长的气相空间；2) 气相空间边界较复杂，尤其是搅混流中某些局部位置(图4b)。

图3 实验工况流型图Fig.3 Flow pattern map of experimental case

图4 高速摄像仪拍摄的图像Fig.4 Image of high speed camera

气相空间的连续性导致有限的拍摄区域内存在不完整的气泡和气相空间，边界复杂性导致气相空间没有规则的几何形状，难以计算其表面积和体积，因此需探索合适的计算方法。同时复杂的相边界会对光产生反射与折射[11]，导致数据处理过程中无法准确获取相界面轮廓，这也给相界面边界识别带来了难度。

2.2 图像处理

根据已有图像处理的研究[12-15]，针对摄像法拍摄的两相流图像进行气相空间识别，识别过程如下：增强对比、阈值分割、孔洞填充和分水岭分割。

搅混流图像处理过程如图5所示，通过增强对比和阈值分割，可初步识别相边界(图5c)。阈值分割后的图像，由于边界复杂性导致某些位置相边界不连续，进而导致孔洞填充后存在不能识别的气相空间(图5d)。边界不连续的位置可称为“断点”，本文采用手动补全的方法进行优化调整，直到完全识别所有气相空间(图5e)。

手动补全过程如下：以“断点”前后相界面边界像素点宽度为补全线条的宽度，对照拍摄原图边界走向对“断点”进行补全，补全后与原图进行反复对比，以尽量减小误差。

采用分水岭分割的方法将相互粘连的气泡分隔开，由于相界面的复杂性，存在将不规则气相空间错误分割的现象(图5f)，本文通过手动调整错误分割的方法进行调整(图5g)。由于环状流工况图像同样存在“断点”的现象，使用手动补全的方法对阈值分割后的图像进行处理，图像处理过程如图6所示。

图5 搅混流图像处理过程Fig.5 Image processing of churn flow

图6 环状流图像处理过程Fig.6 Image processing of annular flow

2.3 空泡份额与界面浓度

通过图像处理可得到清晰、准确的相界面二值化图像。从图5、6可看出，拍摄区域内同时存在较小的气泡和较大的气相空间。通过程序识别方法可得到每个气泡、气相空间的投影面积和投影周长。根据投影面积可计算投影面积直径ds(mm)：

(1)

式中，S为识别出的每个气泡、气相空间的投影面积，mm2。

在窄矩形通道内的两相流动中，根据气泡的大小可将气泡分为两类[16]：1) 气泡投影面积直径ds小于或等于通道宽度s，即ds≤s(本文中s为2 mm)时，认为是球形气泡；2) 气泡或气相空间较大(ds>s)时，将会受到通道宽边的挤压，变为压扁的气泡或气相空间。由于二者形态不同，因此具有不同的体积和表面积计算方法。

1) 球形气泡的体积和表面积(ds≤s)

ds≤s时，球形气泡的体积Vb(mm3)和表面积Ab(mm2)由下式计算：

(2)

2) 压扁气泡和气相空间的体积和表面积(ds>s)

Shen等[16]在窄矩形通道内两相流动研究中，假设压扁气泡边缘为半圆形。同样地，本文假设受挤压的气泡和气相边缘为半圆形，如图7所示。

图7 压扁气泡(a)和气相空间(b)示意图Fig.7 Sketch map of compressed bubble (a) and gas-phase space (b)

图7中，长方形区域为高速摄像拍摄区域，通过程序识别出的投影面积S是图7a中蓝色实线围成的面积。图7b中的投影面积S是蓝色实线与上下拍摄边界围成的面积。程序识别出的投影周长L是图7a、b中为蓝色实线的周长。图中蓝色实线是压扁气相最外圈的边界，虚线是气相空间与窄矩形通道宽边内壁面直接接触的边界。从图7a、b可看出，假设压扁气相边缘为半圆形的区域为虚线与实线围成的部分。

从图5～7可看出，气相空间较复杂，且有不同的长径比。拥有相同投影面积而不同长径比的气相，其周长不同，直接影响气相体积和表面积。长径比对周长L的影响示于图8。由图8可计算得，在投影面积相同的情况下，长径比为4的规则椭圆形压扁气泡的周长，较投影面积为正圆形(长径比=1)压扁气泡的周长大46%。因此，在相界面复杂的工况中，若将拍摄图像内所有受挤压的气相都等效为投影面积为圆形的气泡进行计算，表面积和体积的计算结果将产生较大的误差。

为准确求得上述气泡、气相空间的体积和表面积，提出将其拆分求解的思想。将任意受挤压的气泡和气相拆分为中央部分和外圈部分，如图9所示，中央部分为直接与窄通道壁面接触的气相部分，外圈部分为受到挤压、横截面为半圆形的气相部分。

图8 相同投影面积下不同长径比气泡的周长对比Fig.8 Perimeter comparison of bubbles with different aspect ratio in same area

图9 压扁气泡/气相空间分解图Fig.9 Decomposition diagram of compressed bubble/gas-phase space

对拆分后的中央部分采用面积等效法(即等效前后投影面积相等)，等效为圆柱体进行参数计算。对拆分后的外圈部分采用周长等效法(即等效前后周长相等)，等效为圆环进行参数计算。对等效后的圆柱体和圆环分别计算体积和表面积，再将结果相加，即可得到拆分前气泡或气相空间的体积和表面积。

这样既保留了原气相中央部分的实际体积与表面积，又保留了因气相复杂形态对气泡体积和表面积产生的影响。

经过拆分等效后的圆柱体和圆环的体积和表面积计算方法如下。

1) 圆柱体的体积和表面积

等效后圆柱体的体积和表面积与拆分前原气相中央部分的体积和表面积相等，计算公式如下：

(3)

(4)

式中：V柱体为圆柱体的体积，mm3；A柱体为圆柱体的前后表面积之和，mm2；r为等效后圆柱体的半径，mm。本实验中窄矩形通道宽度s为2 mm。

2) 圆环的体积和表面积

等效后圆环的周长与原气相周长相等，因此求得的圆环体积和表面积即为等效前气相外圈部分的体积和表面积。圆环的体积和表面积通过旋转体积分公式求解：

V圆环=π2RLs2/4+πs2/6

(5)

A圆环=π2RLs+πs2

(6)

式中：V圆环为圆环的体积，mm3；A圆环为圆环的外侧表面积，mm2；RL为圆环的外圈半径，mm。

(7)

式中，L为气泡投影周长，mm。

因此，结合式(3)～(6)，则压扁气泡的总体积与表面积为：

(8)

通过图像识别程序，可得到每幅图像内每个气泡和气相的周长和投影面积。以投影面积直径2 mm为分界线，将气泡进行分类，分别代入式(2)和(8)，可求得每个气泡的体积和表面积。

对某张图像内所有气泡体积求和，然后除以通道的体积即为此张图像的空泡份额，对所有气泡表面积求和后除以通道体积即为此张图像的界面浓度。选取某工况下某位置处某段拍摄时间内的多张照片，对每张图片的空泡份额和界面浓度求平均值，即为此位置在此工况下的时均空泡份额和界面浓度。

(9)

(10)

式中：n为每张图像内气泡的数目；N为某段时间内图像的数目；Vb,i为第i个气泡的体积，mm3；Ab,i为第i个气泡的表面积，mm2；Vchannel为拍摄区域内通道的体积，mm3。

2.4 局部参数模型

由于窄矩形通道横截面具有较大的长宽比，因此在通道宽边方向上，空泡份额和界面浓度具有一定分布特性。为研究其分布特性，需将拍摄图像沿通道宽边方向划分为多个区域，计算每个区域内的局部参数，如图10所示，图中a、b区域分别为搅混流和环状流工况下的一处划分区域。

图10 划分区域示意图Fig.10 Schematic diagram of area division

通过图像识别程序，可获得局部区域内每个气相空间的面积Sloc和周长Lloc。图10中黑色区域为液相，白色区域为气相。局部区域内每个白色区域面积为Sloc，每个白色区域与黑色区域接触界面的周长为Lloc。

从图10可看出，划分区域内的气泡往往不是完整的，当投影面积直径小于等于窄矩形通道宽度(ds≤s)时，假设气泡为球形，此时气泡体积和面积按照球形公式(式(2))计算。此假设可能将大气泡的一角当作是球形气泡进行计算，会导致结果存在一定的偏差。

当局部气相的投影面积直径大于窄矩形窄边宽度(ds>s)时，其体积和表面积计算方法如下：

Vloc=sSloc-(s(s/2)-π(s/2)2/2)Lloc

(11)

Aloc=2Sloc+(πs/2-s)Lloc

(12)

式中：Vloc为局部气相的体积，mm3；Aloc为局部气相的表面积，mm2。

由于假设中气泡边缘为半圆形，局部气相的体积若按照柱体体积来计算，会忽略受挤压的气相边缘为半圆形，将部分液相区域算作气相，导致计算结果偏大。因此需减去多余的部分体积，如式(11)所示。局部气相的表面积若按照2倍投影面积来计算，同样忽略了气相边缘为半圆形的部分，导致计算结果偏小，因此需加上部分表面积，如式(12)所示。

体重增加的速度和多少是影响肚中黑线产生的重要因素，所以预防肚中黑线的有效方法就是不要快速发胖，孕期体重应该渐进式地增加。建议整个怀孕过程中的体重增加控制在11～14千克。

当分区内无气泡边界，即Lloc=0时，如图10中环状流条件下的b区域，式(11)、(12)变为：

Vloc=sSloc

(13)

Aloc=2Sloc

(14)

因此窄矩形通道宽边方向上局部空泡份额αloc与局部界面浓度ai,loc(m-1)的计算式为：

(15)

(16)

式中：Vloc,i为分区内第i个气相空间的体积，mm3；Aloc,i为分区内第i个气相空间的表面积，mm2；Vchannel,loc为局部分区内的通道体积，mm3。

3 结果验证

对每个实验工况下的两相流动视频，以0.1 s为时间间隔进行图像提取和数据处理，每个工况处理数据时长为60 s。对得到的空泡份额、界面浓度和局部参数进行验证。

1) 空泡份额

Ishii[17]提出的适用于窄距形通道的空泡份额经验关系式如下：

jg/α=C0j+Vgj

(17)

(18)

式中，Vgj为漂移速度，m/s。

Ishii给出的适用于圆管的漂移速度经验关系式如下。

搅混流：

(19)

环状流：

Vgj=(1-α)·

(20)

Woldesemayat等[18]证明了此经验关系式适用于窄矩形通道。将实验得到的上、中、下游空泡份额与经验关系式进行对比，如图11所示。可看出，实验得到的空泡份额与经验公式的相对偏差大部分在10%以内。

图11 空泡份额实验值与经验公式计算值的对比Fig.11 Experiment value vs calculated valueof empirical formular for void fraction

2) 界面浓度

由于暂无窄通道内搅混流和环状流界面浓度经验关系式，自定义系数K来验证界面浓度计算方法的合理性。K表示界面浓度与空泡份额的相关性，定义如下：

(21)

从图5、6可看出，在搅混流和环状流区间内，空泡份额越大，相界面边界越简单，中央气空间的几何形状越接近于圆柱体。当拍摄空间内充满气体时，界面浓度和空泡份额可近似按照柱体的计算方法计算。则K的表达式如下：

(22)

当空泡份额较小时，搅混流相界面较复杂，此时相界面面积中侧面积部分占比增加，因此K应大于1，且搅混程度越深，K值偏离1越大。

将不同工况下的K值作图，如图12所示。

图12 空泡份额与界面浓度的相关性Fig.12 Correlation between void fraction and interfacial area concentration

从图12可看出，随着空泡份额的增加，K>1逐渐向K=1靠近，且总是大于1，符合上述理论分析，因此认为界面浓度计算方法是合理的。

3) 局部参数

为证明分区求解局部参数计算方法的合理性，将拍摄区域沿宽边方向等分为10份。对每份区域内计算得到的空泡份额和界面浓度求平均值，与整体图像所得结果进行对比，如图13所示。

空泡份额和界面浓度分区计算结果与未分区计算结果整体相对偏差在5%以内。存在偏差的主要原因在于分区参数计算过程中，将所有投影面积直径小于窄矩形窄边宽度的气相(ds≤s)，均假设为球形来计算。此假设下，可能将大气泡的一角当作球形气泡，偏离了实际情况。

图13 整体图像计算结果与分区计算结果的对比Fig.13 Whole image calculation compared with partition calculation results

4 结论

本文通过搭建窄矩形通道内空气-水向上流动实验回路，获取了搅混流、环状流流型下相界面图像，提出了适用于复杂相界面工况下的空泡份额和界面浓度的计算方法，得到如下主要结论。

1) 基于气相拆分的思想，实现了窄矩形通道内搅混流和环状流下空泡份额的计算。得到的空泡份额与经验公式计算值的相对偏差在10%以内，证明了空泡份额计算方法的合理性。

2) 基于气相拆分的思想，实现了窄矩形通道内搅混流和环状流下界面浓度的计算。自定义界面浓度与空泡份额的相关系数K，实验中K值随着空泡份额的增加逐渐趋近于1，符合理论分析情况。

3) 基于窄矩形通道的特点和分区计算的思想，提出了适用于计算分区通道内，局部位置空泡份额和界面浓度的计算方法。将拍摄区域沿通道宽边等分为10份，求得分区内空泡份额与界面浓度的平均值，与未分区计算得到的参数进行了对比，相对偏差整体在5%以内。