基于OpenFOAM的过冷流动沸腾数值模拟

秦 浩，王明军，李林峰，田文喜，苏光辉，秋穗正

(西安交通大学 陕西省先进核能技术重点实验室，动力工程多相流国家重点实验室，陕西 西安 710049)

过冷流动沸腾指主流中液相温度低于饱和温度，而近壁面处发生局部沸腾的两相流动换热现象。相比于单相对流换热，过冷沸腾具有较高的换热效率，因此在能源与动力工程等行业中被广泛应用。国内外学者均在探索基于计算流体动力学(CFD)方法来研究两相流动换热问题[1-2]。Krepper等[3]采用Bartolomei等[4]开展的压力范围在1.5～15 MPa间的实验数据对商用软件CFX中的过冷沸腾模型进行了验证。Cheung等[5]研究了核化密度、气泡脱离直径和气泡脱离频率等不同模型组合对低压过冷流动沸腾实验的适用性，发现任一经验关系式组合的计算结果都不能与所有低压实验数据良好符合。Zhang等[6]研究了不同湍流模型对数值模拟结果的影响，认为k-ε模型较k-ω模型效果好。陈二峰等[7]和李松蔚等[8]修正了CFX中的气泡脱离直径模型，模拟了低压工况下的过冷沸腾。

目前主流的商用CFD软件中均包含两相流的求解模块，可基于欧拉两流体模型或VOF(体积分数)模型求解两相流场。但商用软件封装严格，开放程度弱，对用户的自主开发有明显的限制，在求解复杂的两相流问题时，很多模型难以修改或添加，从而不一定能满足用户的科研需求。而采用C++语言编写的开源CFD类库OpenFOAM以其面向对象编程、大规模并行能力及开放的编程环境等优势逐步在学术及工程领域得到广泛应用。目前，OpenFOAM在海洋与船舶工程、化学工程等方向的数值模拟已取得了良好的效果。基于该开源平台添加模型或开发新的求解器以模拟过冷流动沸腾现象将是两相流研究的一个重要方向。

本文以4.5 MPa下竖直圆管内过冷流动沸腾现象为研究对象，阐述数值模拟所采用的数学物理模型，包括欧拉两流体模型及相关辅助模型，基于OpenFOAM平台进行模拟。

1 数学物理模型

本文基于欧拉两流体六方程模型，考虑了气液两相间质量、动量和能量的传递，并引入壁面热流密度分配模型以描述壁面沸腾现象。

1.1 相间动量交换模型

由牛顿第三定律可知，气相对液相的作用力Fgl与液相对气相的作用力Flg大小相等、方向相反，即：

Flg=-Fgl

(1)

因此本节仅对液相对气相的作用力进行描述。弥散在液相中的气泡受到来自液相的作用力可分为曳力Md和非曳力，非曳力又包括升力Ml、壁面润滑力Mwl、湍流耗散力Mtd和虚拟质量力Mvm，则气泡总的受力情况可表示为：

M=Md+Ml+Mwl+Mtd+Mvm

(2)

1.2 壁面热流密度分配模型

广泛应用于流动沸腾模拟的壁面热流密度分区模型是RPI模型[9]，它由提出者Kurul和Podowski就职的伦斯勒理工学院(RPI)而得名。该模型将从壁面传递到流体中的热流密度分为：1) 单相过冷流体强迫对流换热Qc；2) 气泡脱离壁面时单相过冷液体重新覆盖壁面而引入的淬灭热流Qq；3) 液相蒸发引入的热流Qe。因此，壁面总热流Qtot表示为：

Qtot=Qc+Qq+Qe

(3)

但上述表达式未考虑气泡覆盖壁面时壁面与气相之间的单相对流换热Qv。为完善物理模型，提高模拟的准确性，改进的热流密度分配模型逐渐得到了广泛应用。

Qtot=f(αl)(Qc+Qq+Qe)+(1-f(αl))Qv

(4)

式中，f(αl)为关于液相份额αl的经验关系式，本文采用Lavieville等[10]提出的模型：

(5)

式中，αl,crit为液相份额的临界值，本文取为0.2。f(αl)函数的图像示于图1，当液相份额低于0.2即空泡份额大于0.8时，f(αl)的值迅速减小，壁面热流密度主要表现为气相单相对流换热。

单相液体的对流换热Qc表达式为：

Qc=(1-Aw)hc(Tw-Tl)

(6)

式中：Aw为气泡影响的区域面积；Tw为壁面温度；Tl为主流液体温度；hc为单相液体对流换热系数。

图1 f(αl)函数图像Fig.1 f(αl) function graph

壁面骤冷过程中的换热Qq可表示为：

Qq=Awhq(Tw-Tl)

(7)

式中，hq为骤冷换热系数。

(8)

式中：f为气泡脱离频率；twait为气泡等待时间；kl为液相导热系数；ρl为液相密度；cpl为液相比热容。气泡等待时间twait可近似表示为：

twait=0.8/f

(9)

蒸汽产生带走的热量Qe通过汽化过程中的传质速率表示：

Qe=mwhlg

(10)

式中：hlg为汽化潜热；mw为汽化的传质速率，其经验关系式为：

(11)

式中：ρg为气相密度；dw为气泡脱离直径；N为核化密度。

由于换热机理不同，单位加热面积可划分为核态沸腾气泡影响区Ab和非影响区1-Ab[11]。Valle等[12]研究认为，气泡影响面积份额可写成单位加热壁面上核化数目和气泡投影面积的表达式：

(12)

式中，K为沸腾模型常系数。Valle等[12]推荐的K计算关系式为：

K=4.8e-Ja/80

(13)

Ja=ρlcpl(Tsat-Tl)/ρghlg

(14)

式中，Tsat为液相饱和温度。

1.3 辅助模型

1) 气泡脱离直径模型

气泡从加热壁面上脱离时的直径称为气泡脱离直径，其大小取决于气泡长大过程中的受力情况。本文采用Tolubinsky等[13]根据实验观测得出的经验模型计算气泡脱离直径dw。

dw=max(drefexp(-(Tsat-Tl)/Tref)，dmax)

(15)

式中：参考直径dref和气泡最大直径dmax分别取0.6 mm和1.4 mm；参考温度Tref通常取45 K。

2) 壁面核化密度模型

单位加热壁面面积上气化核心的数目称为核化密度。本文采用Lemmert-Chawla模型[14]计算，该模型认为核化密度N与壁面过热度ΔTsup有关：

N=(CNΔTsup)n

(16)

式中：参数CN一般取210；n=1.805。

3) 气泡脱离频率模型

气泡脱离频率是加热壁面上气泡脱离速度的量度。Cole等[15]通过高速摄像技术研究了临界热流密度附近的沸腾现象，并推荐气泡脱离频率f的关系式为：

(17)

2 算例与结果分析

2.1 算例描述

Bartolomei等[16]于1967年设计搭建了竖直加热圆管内流动沸腾实验回路，进行了不同压力、热流密度、质量流量以及入口过冷度下的过冷沸腾实验，测量获得了壁面温度、流体截面平均温度、截面平均空泡份额等沿圆管高度的分布。文献[16]中公开的实验段及实验条件示于图2。本文将以此实验为研究对象，进行数值模拟。

图2 实验段示意图Fig.2 Schematic diagram of experimental section

2.2 计算方法

为减少计算量，本文取计算区域为二维模型，即取圆管通道的1个竖直矩形截面，上下面分别为速度入口和压力出口边界条件，左侧为对称边界条件，右侧为加热面。液相考虑其物性随温度的变化，气相采用饱和温度下的物性值。曳力采用Schiller-Naumann模型，壁面润滑力采用Antal模型，湍流耗散力采用Burns模型。由于升力和虚拟质量力对该问题的计算结果影响不大，且添加升力后会使程序计算的收敛性变差，因此本文在计算中不添加升力模型。液相湍流模型采用标准k-ε模型。梯度离散采用Gauss linear格式，其他采用一阶迎风差分格式。计算残差要求小于10-4。

2.3 网格无关性分析

过冷沸腾由于其局部加热的特性，近壁面参数与网格相关，因此不能采用近壁面参数作为网格无关性判断的依据。本文选取沿轴向长度变化的截面平均空泡份额作为网格无关性的考核依据，计算结果如图3所示。由图3可知，轴向控制体网格数对计算结果的影响不明显，而当径向控制体网格数为30与35时，计算结果相差不大。所以，经过网格无关性计算，在径向及轴向分别取30和2 000个节点，总网格数为60 000。此时，网格尺度y+在40～90之间，满足两相流计算要求(20

图3 网格无关性分析Fig.3 Grid independence analysis

2.4 计算结果分析

截面平均空泡份额、液相平均温度和壁面温度的数值模拟结果与实验值的对比示于图4。空泡份额与反应堆的稳定性及堆芯中子动力学密切相关，因此空泡份额的准确预测在核反应堆的设计及运行中非常重要。由图4可知，数值模拟结果与实验值符合良好，说明OpenFOAM中嵌入的欧拉两流体模型及辅助模型可很好地预测该实验条件下的流动沸腾现象。

图4 计算值与实验值的对比Fig.4 Comparison between calculation results and experimental data

图5 计算结果云图Fig.5 Cloud map of simulation result

计算的液相流速、液相温度及空泡份额示于图5。可看出，在流道出口处气相占比较大，挤占了大部分流通面积，致使液相流速提高；加热管道中，在0.9 m处开始产生气泡，1.3～2.0 m处横截面上空泡份额的峰值并不是出现在壁面上，而是在距离加热壁面约R/5位置处，这是气液两相间作用力导致的，其中，壁面润滑力推动气泡向流道中心运动，湍流耗散力使气泡沿空泡份额梯度方向运动，故二者均有展平空泡份额在径向上分布的效果。

3 结论

本文基于开源CFD平台OpenFOAM，采用欧拉两流体模型对4.5 MPa下竖直圆管内的过冷流动沸腾现象进行了数值模拟。计算中考虑了气相与液相间的曳力、壁面润滑力和湍流耗散力，热流密度分配模型采用了考虑气相传热的四分模型，气泡脱离直径、核化密度和气泡脱离频率模型分别采用Tolubinsky模型、Lemmert-Chawla模型和Cole模型。计算得到了空泡份额、液相平均温度和壁面温度沿圆管轴向的分布，计算值与实验值符合良好，证明上述模型对该工况下过冷流动沸腾现象有较好的模拟能力。