VOF模型界面传质与体积传质的转换方法

谭思超，赵富龙，李少丹，高璞珍

（1.哈尔滨工程大学核安全与仿真技术国防重点学科实验室，黑龙江哈尔滨150001；2.清华大学核能与新能源技术研究院，北京100084）

VOF模型界面传质与体积传质的转换方法

谭思超1，赵富龙2，李少丹1，高璞珍1

（1.哈尔滨工程大学核安全与仿真技术国防重点学科实验室，黑龙江哈尔滨150001；2.清华大学核能与新能源技术研究院，北京100084）

针对流体体积模型（volume of fluid，VOF）在模拟相间传质过程时需要将界面质量流密度转换为单位体积传质速率的问题，对VOF模型中的界面传质与体积传质转换方法进行了改进。提出一种既解决了网格无关性问题又可以反映局部界面传质特性的转换方法，并给出了理论推导证明。推导出相应的网格无关性条件，通过划分3组不同尺寸的网格模拟了汽泡冷凝问题，并在汽泡生长方面对该转换方法加以简单应用，模拟结果与理论分析及实验结果符合良好，合理可行，可广泛应用到多相流中界面传热传质模拟的过程。

VOF模型；界面；局部传质；理论推导；网格无关性；多相流

过冷沸腾由于其换热效果强烈，能较好地满足核工业体积小、能量密度高的要求，受到越来越多研究者的关注；窄通道由于结构紧凑、换热性能好等特点，得到越来越多的重视。为了较好地研究强化换热特性，大量研究者对窄通道内的汽泡行为进行了研究，大部分是基于实验研究得出相关结论。汽泡行为的研究中，数值模拟是一种重要的手段，其中VOF模型凭借较好的界面追踪特性得到越来越多研究者的青睐，Jeon等［1⁃6］分别采用VOF模型对过冷沸腾下汽泡的生长及冷凝行为进行数值模拟。

用VOF模型进行相变过程计算时，采用体积传质速率进行计算，但一般已知多为传质速率或者质量流密度，因此必然涉及到界面传质速率与体积传质速率间的转换。目前，在使用VOF模型进行界面传热传质特性计算时，主要有蒸发冷凝模型和平均处理方法。蒸发冷凝模型，是Lee［7］提出的体积传质处理方法，主要思想是，温度高于液相的饱和温度时，液相蒸发，对应到相应的网格即为液相向汽相传质；反之，汽相冷凝。该模型形象直观，易于理解，但结果对网格尺寸有很强依赖性，假定界面层厚度不变，由于网格体积与尺寸之间为三次方关系，若划分的网格越细，则界面上的网格数量呈三次方变化。该模型单个网格的体传质速率保持不变，会使得相间界面总传质速率随网格数量呈现指数增加，在网格无关性方面存在一定的缺陷。平均处理方法是Jeon［1］在用VOF模型对过冷流动沸腾条件下汽泡冷凝行为模拟时提出的，进行汽泡冷凝传质的计算时，先计算汽液界面总传质量，随后在每个网格进行平均。尽管考虑了网格尺寸的作用，解决了网格无关性问题，但由于涉及平均处理过程，不能很好地反映汽液界面的局部传质过程。

综上所述，现有VOF模型传热传质计算方法均存在一定的不足。为此，在理论分析的基础上提出一种VOF模型的界面传质与体积传质转换的处理方法，用于计算涉及多相相间界面传热传质问题，以便更好的用于VOF模型相变过程的计算。

1 汽泡描述方法及传热传质模型

1.1 汽泡描述方法

VOF［8］通过各相的体积份额来描述不同的流体相，对于汽液两相流，单元中汽相体积份额αv为0代表液相，αv为1代表汽相，介于两者之间为汽液界面；在固定的欧拉网格下，求解各相体积份额的连续性方程来描述2种或者多种互不混合流体界面；各相共享单一动量、能量方程，各单元的密度、粘度等参数采用体积份额进行加权平均求得。

文中采用几何重构方法描述汽液界面，即采用分段线性方法描述汽液界面，假定两相间界面在每个单元内有一个线性斜面，使用此线性形状计算穿过单元面的对流及位置信息，示意图如图1，白色弧线代表汽液界面，弧线左上方黑色区域代表液相，弧线右下方区域代表汽相。

图1 几何重构方法描述的汽液界面形状Fig.1 Liquid⁃vapor interface shape respresented by the geometric reconstruction scheme

按照几何重构的方法，对汽泡的描述也是通过体积份额来刻画，汽液界面内部即为典型的汽泡，如图2所示。

图2 VOF模型刻画的汽泡剖面图Fig.2 Bubble section map by VOF model

需要指出的是，典型汽泡的刻画是需要多个网格来共同组合描述的，Rabha等［9］认为网格数／汽泡直径（cells per buddle diameter，CPD）应至少为16，即当汽泡直径与网格尺寸比值大于16时才能较好的刻画汽泡，下文提出的体积转换方法的网格要求也在这一范围内，说明了改进方法的网格无关性条件符合VOF模型描述汽泡的基本要求。

1.2 汽液界面传热传质模型及UDF

过冷流动沸腾汽泡行为方面，目前公认的汽泡生长动力为微液层蒸发与汽液界面蒸发冷凝综合作用的结果［10］，本文采用经典的体现界面传质特性的Hertz⁃Knudsen公式来进行汽液界面传质过程计算［11］，其表达式如下：

式中：θ为界面蒸发／冷凝系数；R为气体常数，J／mol-1·K-1；M为摩尔质量，kg／mol；Tl为液相温度，K；ps为Tl对应的饱和压力，MPa；pv为汽泡内压力，MPa；Tv为汽相温度，K。

根据上述数学模型编写UDF程序，加入相应的质量源项和能量源项，能量源项是在质量源项的基础上考虑汽化潜热的作用得到的。

2 体积传质转换方法

2.1 基本假设

在理论推导之前作以下几点假设：1）汽泡为已从加热壁面脱离位于管道中央的孤立汽泡；2）所划分的网格为正方体结构性网格，各个网格的边长、表面积、体积均相等；3）传质过程为均匀传质过程；4）汽液界面厚度均匀；5）在计算过程中假设汽泡内为饱和状态。

汽泡内压力由Young⁃Laplace公式来确定：

式中：pv和pl分别为汽泡内外压力，σ为液体表面张力，汽泡的半径R是通过在计算过程中测定汽泡在x、y、z 3个方向的高度：

2.2 理论推导

VOF模型进行多相流计算时，汽液界面是通过多个网格间界面重构方法构建的。每个网格内都充满流体，若流场中只含有汽液两相，则

式中：αvi、αli分别为第i个网格中汽相和液相的体积份额。

从式（4）可知，用VOF模型对两相流进行计算，蒸发与冷凝相当于2个逆向过程，即在VOF模型中适合蒸发过程的算法也适用于冷凝过程。若过程为蒸发过程，则相间传质界面面积可近似为αviSi，其中Si为单个网格截面积。采用这种近似方法进行处理，是因为VOF模型各相靠体积份额描述的，汽液界面处的网格单元中既有汽相又有液相，相间有一定的接触面，传质过程就是在接触面上进行的。已知网格表面积，若已知传质界面面积与网格表面积的比例关系，便可得到传质界面面积，但在VOF模型中传质界面的面积不能通过上述方法得到，能得到的物理量只有各相的体积份额，为此采用体积份额对相间传质界面面积近似处理，即认为相间传质界面的面积约为αviSi。

按上述近似处理方法，则单个网格的传质速率近似为αviGiSi，汽相所占的体积为αviVi，其中Vi为网格的体积，根据VOF模型中的体传质速率的定义可得到单个网格的体传质速率mi为

式中：Gi为质量流密度，Li为网格尺寸。式（5）即为本文中提出的将质量流密度转换为体传质速率的处理方法。下面在理论推导基础上对提出的转换方法进行网格无关性的分析。

汽液界面上总的传质速率Mt表达式为

结合式（6）、（7）算汽液界面总传质速率Mt：

从式（8）可以看出，当汽泡尺寸一定时，要以一定的精确度ε确保汽液界面上总的传质速率Mt与网格尺寸无关，则需满足下面的条件：

即

式（10）即为改进方法的网格无关性条件。可以看出要保证计算精度为10%，则要求汽泡半径约为网格尺寸的10倍，下文用此精度进行方法验证。

3 方法验证

3.1 几何模型及方法验证

窄通道一般是指间隙尺寸相对很小，进出口长度与常规通道相当的通道。Kandlikar［12］认为当量直径在200 μm～3 mm的通道称为窄通道或小通道。而根据工程实际，通常将当量直径在1～3 mm的通道称为窄缝通道。考虑到进行三维数值计算网格量及计算精度要求，本文选取窄矩形通道几何模型的尺寸为2 mm× 4 mm×2 mm，流体竖直向上流动，重力方向竖直向下。划分3组不同疏密程度的六面体结构性网格：1）50× 100×50，2）64×128×64，3）80×160×80，相应的六面体网格的边长分别为0.04、0.031 25、0.025 mm，并对3组不同网格尺寸的汽泡的冷凝速率及冷凝过程中的汽泡形状图像进行对比。通过对汽泡冷凝过程的三维数值模拟结果分析，对提出的转换方法进行验证。

边界条件的选择，进口边界条件用速度入口，为稳定工况，出口边界条采用压力出口。离散方法的选择，考虑到求解速度、稳定性、精度等因素，采用双精度进行计算，压力速度耦合采用PISO算法，使用PRESTO方法离散压力，动量和能量方程使用二阶迎风格式求解，瞬态方程的离散采用一阶迎风格式求解，容积比率方程用几何重构方法求解，非稳态时间步长为5 μs。

3.2 模拟结果及工况

为了证明模型的准确性，参考潘良明等［3］的实验工况进行了模拟，所选择的计算条件如下表1所示。选择该文献进行对比主要是因为其采用VOF模型对过冷流动沸腾条件下汽泡冷凝过程进行数值模拟，其模拟结果与相应条件下的实验结果吻合较好，误差在±20%之内，表明能够反映真实参数变化过程。

表1 方法验证模拟工况参数Table1 Simulation parameters for method verification

计算过程中，所有的液相参数根据过冷水的温度和相应的压力查得，汽相参数根据相应的饱和温度查得，随着冷凝过程的进行汽泡的直径会逐渐减少，为保证计算精度，同时验证网格无关性条件，需要对网格进行实时加密，保证网格与汽泡半径的比值与初始时刻的比值相同。所得到的冷凝过程中的3组不同疏密程度网格对应的汽泡冷凝曲线如下图3所示，其中，1、2、3分别代表上述3组网格对应的模拟结果，上标“'”表示文献［3］的实验结果。b工况冷凝过程中汽泡的形状变化图像如图4所示，在图4中（a）、（b）、（c）网格结构分别为：50×100×50、64×128×64、80×160×80。

图3 汽泡冷凝曲线Fig.3 Bubble condensation curves

图4 模拟得到的汽泡生长过程形状Fig.4 The simulated shapes during bubble growth

从图3中，可以看出本文的计算结果与PAN等［5］的实验结果趋势基本一致，冷凝时间差别不大，气泡直径最大误差为14.8%，表明本文所提出的界面传质与体积传质的转换方法准确。从图3可以看出，汽泡冷凝过程中，不同疏密网格计算的汽泡直径变化吻合良好，最大误差在±6.3%以内。对于b工况，通过对图3中“b⁃1，b⁃2，b⁃3”汽泡半径随着时间的变化曲线的对比，可以看出不同疏密的网格计算结果吻合良好，冷凝过程一致，说明了上述的VOF模型中汽液两相相间传质的转换方法中网格无关性条件正确。

图3中的模拟结果表现出线性特性，主要是因为数值模拟中的相间传热传质过程的实现是按照经验公式人为编写程序加入的，而现有的经验公式考虑的因素有限，最终会导致结果变现为一定的线性关系；另外，为了保证文中的网格无关性条件，在方法验证过程中，汽泡并没有完全冷凝，如图3中的曲线所示，要想进一步计算则需要划分尺寸更小的网格。

方法验证虽是在均匀传质假设下进行的，但同样适用于非均匀传质过程，如过冷流动沸腾，虽涉及微液层蒸发和经典界面传质公式2种模型作用，但可将整个汽泡划分为2个部分球缺形，分别进行单独计算，同样适用于上述的论证。对于非结构性网格，可根据划分非结构性网格的算法，通过计算给出相应的权重系数来进行分配。

4 方法应用

在转换方法验证的基础上，加以简单应用，将转换方法应用到竖直窄矩形通道中过冷沸腾条件下汽泡在热力生长特性的三维模拟仿真。

4.1 数学模型

过冷流动沸腾下汽泡生长动力一方面来自汽液界面的传热传质，另一方面汽泡生长过程中微液层不断蒸发也会成为汽泡生长动力。为此，在计算中加入界面传热传质模型及微液层蒸发模型，界面传热传质模型采用式（1）来计算。微液层蒸发模型，采用Addlesee等［13］对滑移气泡的流场分析得到的滑移汽泡微液层厚度表达式：

式中：vl、H、u分别为液相运动粘度、汽泡高度和滑移速度。则转化的质量流密度为

式中：λ为液相导热系数，hγ为汽化潜热。

传质量以及传热量通过提出的转化方法转化为UDF程序语言，加入到计算过程中。

4.2 结果分析

由于不能实现自主产生汽泡，需要初始给定汽泡，初始直径为0.1 mm。考虑计算速度及精度要求选择管道尺寸为1 mm×1 mm×1 mm，不同区域划分不同粗细网格，既保证满足网格无关性条件，又能适当减少计算量，经过测试网格数量为173 781。

工况参数为工作压力0.101 325 MPa、液体过冷度5 K、壁面过热度10 K、流速0.1 m／s。边界条件及离散方法等与“方法验证”一节中相同。计算得到的汽泡生长过程如图4所示，重力方向竖直向下。

对汽泡生长曲线进行拟合，得到的拟合曲线见图5，拟合曲线的相关性系数为0.994，表明模拟的汽泡生长遵循指数规律，与大量实验得出的汽泡后期生长阶段规律一致。Thorncroft等［14］对FC⁃87工质进行的大量实验发现过冷沸腾下汽泡生长的指数介于1／3～1／2之间，本文数值模拟结果中指数为0.416 39，吻合较好，该指数n小于0.5，是因为初始给定汽泡尺寸，忽略了汽泡的惯性生长过程。

从图4、5可以看出，本文汽泡传质模型能较好的模拟汽泡热力生长过程；但前期的惯性控制阶段差别较大，是因为在计算时初始给定汽泡尺寸，不是经过自然核化点核化来产生汽泡，无法准确模拟汽泡惯性生长阶段的行为。另外，由于文中更加关注两相间的传质过程，而对汽泡的受力考虑不是很充分，不能直观看出汽泡的滑移运动。

上述简单的应用，表明提出的VOF模型界面传质与体积传质转换方法合理可行，可较好应用到沸腾条件下汽泡演化特性的模拟仿真过程中。

图5 汽泡生长曲线Fig.5 The bubble growth curve

5 结束语

本文在理论分析基础上，提出了一种适用于VOF模型的界面传质与体积传质的转换方法，数值模拟结果与理论推导得到的结果吻合较好，误差在±6.3%以内。该转换方法可以通过理论推导出相应的网格无关性条件，解决了网格无关性问题又，并可以描述局部界面传质特性，是较为适用于VOF模型的界面传质与体积传质转换的处理方法。

通过对竖直窄矩形通道中过冷沸腾条件下汽泡在热力生长阶段行为的三维数值模拟，对方法进行了简单的应用，与现有的方法相比精度较高，为今后多相相间界面相变过程提供了一种准确的数值仿真方法。

［1］JEON S，KIM S，PARK G.Numerical study of condensing bubble in subcooled boiling flow using volume of fluid model［J］.Chemical Engineering Science，2011，66（23）：5899⁃5909.

［2］魏敬华，潘良明，袁德文，等.过冷流动沸腾相变过程汽泡特性的VOF方法模拟［J］.核动力工程，2012，33（6）：65⁃71.WEI Jinghua，PAN Liangming，YUAN Dewen，et al.VOF simulation of bubble characteristics of subcooled flow boiling［J］.Nuclear Power Engineering，2012，33（6）：65⁃71.

［3］潘良明，谭智威，闫晓，等.窄流道内过冷流动沸腾汽泡凝结过程的数值模拟研究［J］.核动力工程，2011，32（6）：86⁃90. PAN Liangming，TAN Zhiwei，YAN Xiao，et al.Numerical investigation of bubble condensation of subcooled boiling in narrow rectangular channels［J］.Nuclear Power Engineer⁃ing，2011，32（6）：86⁃90.

［4］潘良明，谭智威.过冷流动沸腾汽泡凝结变形及流场特性的数值模拟［J］.重庆大学学报，2012，35（6）：53⁃57.PAN Liangming，TAN Zhiwei.Numerical investigation of bubble deformation and flow field characteristics of subcooled boiling during condensation［J］.Journal of Chongqing Uni⁃versity，2012，35（6）：53⁃57.

［5］PAN L，TAN Z，CHEN D，et al.Numerical investigation of vapor bubble condensation characteristics of subcooled flow boiling in vertical rectangular channel［J］.Nuclear Engi⁃neering and Design，2012，248：126⁃136.

［6］袁德文，潘良明，陈德奇.竖直窄流道内过冷流动沸腾的单汽泡生长模型［J］.化工学报，2009，60（11）：2723⁃2728.YUAN Dewen，PAN Liangming，CHEN Deqi.Model for sin⁃gle bubble growth of subcooled flow boiling in vertical narrow rectangular channel［J］.Journal of Chemical Industry and Engineering，2009，60（11）：2723⁃2728.

［7］LEE W H.Pressure iteration scheme for two⁃phase flow mod⁃eling［J］.Multiphase Transport：Fundamentals，Reactor Safety，Applications，1980，1：407⁃432.

［8］ANSYS Inc.FLUENT theory guide［Z］.Pittsburg：ANSYS Inc，2012.

［9］RABHA S S，BUWA V V.Volume⁃of⁃fluid（VOF）simula⁃tions of rise of single／multiple bubbles in sheared liquids［J］.Chemical Engineering Science，2010，65（1）：527⁃537.

［10］郭烈锦.两相与多相流动力学［M］.西安：西安交通大学出版社，2002：323⁃410.

［11］王遵敬，陈民，过增元.蒸发与凝结现象的分子动力学研究［J］.西安交通大学学报，2001，35（11）：1126⁃1130.WANG Zunjing，CHEN Min，GUO Zengyuan.Molecular dynamics study on evaporation and condensation［J］.Jour⁃nal of Xi'an Jiaotong University，2001，35（11）：1126⁃1130.

［12］KANDLIKAR S G.Fundamental issues related to flow boil⁃ing in minichannels and microchannels［J］.Experimental Thermal and Fluid Science，2002，26（2⁃4）：389⁃407.

［13］ADDLESEE A J，CORNWELL K.Liquid film thickness a⁃bove a bubble rising under an inclined plate［J］.Chemical Engineering Research and Design，1997，75（7）：663⁃667.

［14］THORNCROFT G E，KLAUSNER J F，MEI R.An experi⁃mental investigation of bubble growth and detachment in vertical upflow and downflow boiling［J］.International Jour⁃nal of Heat and Mass Transfer，1998，41（23）：3857⁃3871.

The transformation method of mass flux and mass transfer rate per volume at the interface in VOF model

TAN Sichao1，ZHAO Fulong2，LI Shaodan1，GAO Puzhen1
（1.National Defense Key Subject Laboratory for Nuclear Safety and Simulation Technology，Harbin Engineering University，Harbin 150001，China；2.Institute of Nuclear and New Energy Technology，Tsinghua University，Beijing 100084，China）

The interface mass flux needs to be transferred to mass transfer rate per volume，when simulating mass transfer process between the phases by volume of fluid（VOF）model.In order to solve this problem，the method for transformation from interface mass flux to mass transfer rate per volume in the VOF model is improved.A new trans⁃formation method was proposed in this paper，which could solve the network independence.However，it also re⁃flects the mass transfer characteristics at a local interface.The transformation method was proved by theoretical in⁃ference.The corresponding condition that it is independent of mesh size was deduced and subsequently the simula⁃tion of bubble condensation was conducted by dividing three groups of meshes with different sizes.The transforma⁃tion method was simply applied in the aspect of bubble growth.The simulation results coincided with the theoretical analyses and experimental studies very well，proving that this method is feasible and can be widely applied in the simulation of the mass and heat transfer processes at the interface between the different phases.

VOF model；interface；local mass transfer；theoretical derivation；mesh independence；multiphase flow

10.3969／j.issn.1006⁃7043.201312053

http：／／www.cnki.net／kcms／detail／23.1390.U.20150109.1525.013.html

TL331

A

1006⁃7043（2015）03⁃0317⁃05

2013⁃12⁃17.网络出版时间：2015⁃01⁃09.

核反应堆系统设计技术重点实验室基金资助项目（KZA⁃KA⁃1101）；教育部留学归国基金资助项目（2012⁃1707）；黑龙江省青年学术骨干支持计划资助项目（1254G017）.

谭思超（1979⁃），男，教授，博士，博士生导师.

谭思超，E⁃mail：tansichao＠hrbeu.edu.cn.