湍流模型对过冷沸腾计算影响分析

张 蕊，田文喜，丛腾龙，苏光辉，秋穗正（1．西安交通大学动力工程与多相流国家重点实验室，陕西西安 710049；2．西安交通大学核科学与技术学院，陕西西安 710049）

湍流模型对过冷沸腾计算影响分析

张 蕊1，2，田文喜1，2，丛腾龙1，2，苏光辉1，2，秋穗正1，2
（1．西安交通大学动力工程与多相流国家重点实验室，陕西西安 710049；2．西安交通大学核科学与技术学院，陕西西安 710049）

摘要：采用ANSYS FLUENT14．5分析湍流模型对过冷沸腾的影响，建立多套网格，每套网格采用不同的湍流模型、壁面函数及两相湍流处理方法，将计算结果与基准实验数据进行对比，分析网格、湍流模型、壁面函数及两相湍流处理方法对计算精度的影响。通过分析发现：kε－模型的计算精度高于k－ω模型的；Dispersed方法和Per Phase方法的计算精度相对于Mixture的无明显提升；增强壁面函数对Y＋≈1的网格无法获得满意的计算结果。

关键词：过冷沸腾；湍流模型；计算精度

压水堆堆芯及蒸汽发生器内均存在过冷沸腾，过冷沸腾可极大地提高换热表面的换热系数，但由于传热表面热流密度较高，很容易发生传热恶化，当过冷沸腾热发生偏离核态沸腾（DNB）现象时，由于加热表面被气膜覆盖，加热表面温度瞬间飞升并发生烧毁。因此，准确预测临界热流密度（CHF）、壁面DNB的发生对于燃料棒的完整性具有重要意义。

过冷沸腾和CHF研究通常基于实验进行，通过实验获得过冷沸腾换热系数和CHF数据、拟合关系式，为同工况下的设计分析提供参考。由于经验关系式存在误差，基于经验关系式的设计分析均需留有足够的裕量以保证设备的安全，且经验关系式仅适用于特定的实验区间及几何；另外，对于大型复杂工况，如燃料组件内的过冷沸腾现象，实验极其昂贵且难以测量。随着CFD技术的进步及对过冷沸腾机理的认识，研究者提出不同的计算过冷沸腾的CFD模型［1－2］，并对一些简单几何内的过冷沸腾现象进行数值模拟。Krepper等［3］对单管内的过冷沸腾进行计算，并尝试将PRI（Rensselaer Polytechnic Institute）过冷沸腾模型用于燃料组件的过冷沸腾分析；Yan等［4］对燃料组件内的过冷沸腾现象进行模拟，并得到DNB现象。以上分析将RPI过冷沸腾模型用于实际问题的计算，但这些文献中均未给出湍流模型的选取依据，本文基于RPI过冷沸腾模型和Bartolemei等［5－6］的实验数据，对过冷沸腾计算中的湍流模型进行敏感性分析。

1 数学物理模型

过冷沸腾的数值分析基于欧拉两流体模型进行，对气相和液相分别求解质量、动量、能量守恒方程，同时求解相间质量、动量、能量传递模型。

1．1 欧拉两流体模型

质量守恒方程：

动量守恒方程：

能量守恒方程：

其中：αi、ρi、vi、Si、p、τ＝i、g、vij、Fi、hi、qi和hij分别为i相的体积份额、密度、速度、源项、压力、应力张量、重力加速度、速度差、相间力、比焓、热流量和焓差；ij和Qij分别为i相至j相的质量和能量传递。应力张量表达式为：

其中：μi和λi分别为i相的剪切黏度和体积黏度；T表示湍流；I＝为单位对角矩阵。

1．2 相间动量传递模型

气液相间作用力包括相间拽力、升力、壁面润滑力、湍流耗散力，拽力FD为单位体积内气泡施加在液相上的力：

其中：CD为升力系数，由Ishii关系式［7］计算；μf为液相黏度；dg为气泡直径；vg和vf分别为气、液相速度；Ai，f为相界面密度，即：

Re为相间相对雷诺数：

其中：αg为气相体积份额；ρf为液相密度。

相间升力为由于液相速度梯度而施加在气相上的力：

其中，CL为升力系数，由Moraga关系式［8］计算。

壁面润滑力作用在气泡上，推动气泡向主流方向运动，由下式计算：

其中：Cwl为壁面润滑系数，由Antal模型计算［9］；nw为壁面法向矢量。

湍流耗散力用来描述两相间的湍流相互作用，由下式计算：

其中：CTD＝1；σfg＝0．9；μt，f为动力黏性系数。

1．3 相间能量传递

气相脱离壁面之后，进入过冷流体时向液体释热，单位体积内的传热量为：

其中：Nu采用Ranz－Marshall公式［10］计算；kf为液相的导热系数；Tg为气相温度；Ts为饱和温度。

1．4 相间质量传递－RPI过冷沸腾模型

采用RPI壁面沸腾模型计算壁面上的沸腾现象及气泡的冷凝等传质过程。RPI模型包括壁面热流密度分配模型和辅助模型。

1）壁面热流密度分配模型

RPI模型将壁面传递给液体的热流密度分为3部分，即单相液体对流带走的热流密度qc、淬灭热流密度qq和蒸发热流密度qe。

其中：hc为单相对流换热系数；Af为壁面上液相润湿面积；Tw为壁面温度；Tf为液相温度。壁面液相蒸发产生气泡时带走的qe可表示为：

其中：Vd为基于气泡脱离直径的气泡体积；Nw为气泡核化密度；hfg为汽化潜热；f为气泡脱离频率。qq为气泡从加热壁面脱离、液相重新覆盖壁面时所带走的热流密度，表示为：

2）辅助模型

RPI模型的辅助模型包括液相润湿面积Af、气泡脱离频率模型f、气泡核化密度Nw、气泡脱离直径Dw等模型，计算关系式分别为：

其中，Ja为过冷度的Jacob数。

1．5 湍流模型

本文研究目的为分析不同湍流模型对过冷沸腾计算的适用性，在本研究中，分别采用不同的网格数量、不同的湍流模型对基准实验进行计算，根据网格Y＋不同，所采用的湍流模型包括标准k－ε模型、RNGk－ε模型、可实现k－ε模型、标准k－ω模型、SST k－ω模型。对于k－ε模型，其近壁面网格需满足Y＋在对数率层，即Y＋＞11．225。根据ANSYS FLUENT手册，增强壁面函数可对Y＋＜11．225的工况进行计算，并获得合理的计算结果。考虑到随着流体被加热相变，流体流速逐渐增大，对于相同的壁面网格尺度，Y＋也逐渐增大，因此，采用壁面网格Y＋适应性很强的增强壁面函数。在计算中，还进行了近壁面函数的敏感性分析。另外，对于湍流的计算，可采用Mixture、Dispersed 和Per Phase方法。采用Mixture方法时，求解器将气液相视作混合物，求解1套相应的湍流方程；采用Dispersed方法时，求解器将气相视作离散相，采用上述指定的湍流模型求解液相湍流参数，采用Tchen理论关系式计算离散相湍流参数，同时考虑相间湍流相互作用；采用Per Phase方法时，对每一相分别求解湍流方程，同时考虑相间湍流相互作用。

不同湍流模型对网格的要求不同，k－ε模型［11－13］要求近壁面网格中心节点位于黏性支层外，需采用壁面函数对近壁面区处理，标准壁面函数、非平衡壁面函数要求壁面网格Y＋大于黏性支层的Y＋（约11．225）；ANSYS FLUENT14．5中引入增强壁面函数，使k－ε模型对Y＋无要求，原则上Y＋越小计算结果越精确［14］。k－ω模型［15－17］要求近壁面网格Y＋≈1，对黏性支层进行详细的网格划分；ANSYS FLUENT14．5中引入增强壁面函数，使k－ω模型对Y＋的要求降低，对于Y＋≈10～100的网格，也可计算并得到合理的结果［14］。

本文共采用4套网格，采用不同湍流模型进行计算，湍流模型计算矩阵列于表1。其中StdWF、EnhWF、NeqWF和ScaWF分别表示标准壁面函数、增强壁面函数、非平衡壁面函数和扩展壁面函数，M、D和P分别表示湍流采用Mixture、Dispersed和Per Phase方法计算。

2 实验工况及计算建模

Bartolemei等［5－6］于1967年和1982年分别对不同直径的竖直圆管内向上流动的过冷沸腾进行研究，实验获得不同压力、热流密度、质量流密度、进口过冷度下的空泡份额沿轴向分布以及壁面温度沿轴向分布。实验工况列于表2。Bartolemei共给出20余个工况的实验数据，但其中仅有1个工况（d＝15．4mm，G＝900kg／（m2·s），q＝0．57MW，p＝4．5MPa）给出了沿轴向的内壁面温度分布及截面平均空泡份额分布，可用于模型的验证。

考虑到试验段的对称性，选取1／4区域作为计算流体域，建模不考虑固体域，几何模型及网格划分如图1所示。

表1 湍流模型计算矩阵Table 1 Computing matrix of turbulence model

表2 Bartolemei实验参数范围Table 2 Parameter range of Bartolemei’s experiment

图1 几何建模及网格划分Fig．1 Geometry and grid of calculation domain

3 结果分析

考虑到不同湍流模型及近壁面函数对网格要求不同，本文采用4套网格进行计算，网格数量、近壁面网格几何尺寸、计算得到的湍流Y＋范围列于表3。由表3可看出，ICM1和ICM2的网格满足传统k－ε模型的近壁面网格要求；ICM4满足传统k－ω模型的网格要求；ICM3部分区域近壁面网格中点位于黏性支层以内、部分位于黏性支层之外，属于传统湍流模型和壁面函数难以处理的范畴。

表3 计算选用网格的数量、近壁面网格尺寸及计算得到的Y＋范围Table 3 Grid number，dimension of near－wall mesh and range of Y＋

3．1 网格对计算结果的影响

由ANSYS FLUENT14．5理论手册［13］可知，当湍流模型和增强壁面函数同时使用时，无论是k－ε模型还是k－ω模型对近壁面网格Y＋的要求均降低，为评估该模型，对4套网格分别采用RNGk－ε方程和增强壁面函数进行计算，得到的截面平均空泡份额、流体截面平均温度、壁面温度沿轴向高度的分布如图2所示。由图2a可看出：对于截面平均空泡份额（VOF），ICM1～3的计算结果与实验值符合较好；而网格最密的ICM4的计算结果与实验值偏差较大，采用ICM4计算得到的沸腾起始点晚于实验值。由图2b可看出，流体截面平均温度的计算值均高于实验值，且随近壁面网格Y＋的减小，流体截面平均温度与实验值的偏差变大；总的来说，ICM1～3计算结果与实验值符合良好，ICM4计算值与实验值偏差较大。由图2c可看出，ICM1、2计算得到的壁面温度分布与实验值符合良好，ICM3计算得到的壁面温度分布在0．25～0．35m高度处呈轻微震荡，ICM4计算结果在0．1～1．2m高度处呈锯齿形震荡，且ICM3、4计算结果与实验值偏差较大。综合图2可得：虽然增强壁面函数可允许k－ε模型的近壁面网格Y＋＜11．225，但由分析结果可看出，对于过冷沸腾计算，当Y＋＜11．225时，与近壁面网格的相关的参数（壁面温度）出现不合理的锯齿形震荡。

图2 网格对计算结果的影响Fig．2 Effect of grid on calculation result

3．2 湍流模型对计算结果的影响

图3示出不同湍流模型对计算结果的影响。分析在ICM2的基础上进行，分别采用标准、可实现、RNGk－ε模型，标准和SSTk－ω模型和增强壁面函数法，采用Mixture、Dispersed、Per Phase方法计算湍流。由计算结果可看出，对于ICM2，所有湍流模型的计算结果与实验值均符合较好。在图3a空泡份额的分布中，采用Standard－k－w－EnhancedWF－D方法得到的计算结果与实验值偏差最大；采用不同模型计算得到的截面流体平均温度和壁面温度无明显区别；采用Dispersed和Per Phase方法时，计算得到的空泡份额较Mixture方法的大。根据文献［14］，采用Dispersed和Per Phase方法求解湍流模型相较Mixture方法需要更多的计算时间和内存，亦可获得更高的精度。通过计算分析发现，对于过冷沸腾工况，采用Dispersed 和Per Phase方法计算精度并不比Mixture方法的计算精度高。另外，采用Dispersed和Per Phase方法时计算收敛性降低且计算耗时成倍增加。

3．3 壁面函数对计算结果的影响

图4为采用不同壁面函数得到的计算结果。由图4可看出，当近壁面网格满足壁面函数的要求时，所有的壁面函数均可得到相当高精度的模拟结果。通过对比图4a中增强壁面函数Dispersed和Mixture方法得到的VOF亦可发现，Dispersed方法得到的VOF较Mixture方法得到的偏大。不同壁面函数得到的流体截面平均温度无区别。壁面温度分布稍有不同，采用非平衡壁面函数得到的壁面温度最低且与实验值偏差最大；增强壁面函数的计算结果略逊于标准壁面函数和扩展壁面函数的结果，但总体与实验值符合较好。

图3 湍流模型对计算结果的影响Fig．3 Effect of turbulence model on calculation result

3．4 k－ω模型网格依赖性分析

在k－ω模型［15－17］中，k－ω模型需近壁面网格Y＋≈1，而ANSYS FLUENT14．5［14］通过对k－ω模型引入增强壁面函数近壁面处理，使该模型对近壁面网格精度要求降低，允许Y＋≈10～100量级。本文分别在ICM2和ICM4的基础上采用标准和SSTk－ω模型进行计算，并将计算结果与实验值进行对比，结果如图5所示。由图5可看出，基于ICM2的计算结果与实验值符合较好，且优于基于ICM4的计算结果。在该计算中，两套网格的近壁面网格Y＋列于表4。ICM4满足常规k－ω模型的网格要求［15－17］，但其计算精度却低于ICM2，且壁面温度出现不合理的波动。这是由于当近壁面网格太小时，网格内气泡的瞬时影响无法被时均，从而导致时均计算方法出现波动。

3．5 模型总体性能评价

图4 壁面函数对计算结果的影响Fig．4 Effect of wall function on calculation result

图5 网格对k－ω模型计算结果的影响

Fig．5 Effect of grid on calculation result of k－ωmodel

根据文献［11－13，15－17］，对每组网格分别选取最佳湍流模型（ICM1～3，k－ε；ICM4，k－ω）进行计算，将计算结果与实验值进行对比，结果如图6所示。由图6可看出，对于ICM1～3，基于k－ε模型的计算结果与实验值符合较好；而对于ICM4，基于k－ω模型的计算结果与实验值偏差较大。又由图2可知，对于ICM4，基于k－ε模型的计算结果与实验值偏差也较大，即无论k－ω还是k－ε模型，ICM4的计算结果与实验值偏差均大于其他网格。由图6b可知，ICM1网格、标准k－ε模型计算得到的流体截面平均温度与实验值最为接近。由图6c可知，ICM2网格标准k－ε和可实现k－ε模型的计算壁面温度与实验值最为接近。即当近壁面网格Y＋满足传统k－ε网格要求（Y＋＞11．225）时，计算结果最优。文献［13］指出，当引入增强壁面函数之后，近壁面网格越密越能得到更为准确的结果；文献［18］中通过单相算例证实了这一结论，但在两相过冷沸腾计算中，过于精细的网格不能得到合理的计算结果。

表4 近壁面网格Y＋Table 4 Near wall mesh Y＋for grid

图6 各网格配以最佳湍流模型计算结果Fig．6 Result of each mesh with optimized turbulence model

4 结论

本文采用多组网格、多种湍流模型对过冷沸腾基准实验进行计算，分析湍流模型对过冷沸腾计算的影响，得出如下结论：

1）采用Dispersed方法和Per Phase方法计算两相的湍流参数，相较Mixture方法，计算收敛性和耗时增大，但计算精度并无明显提高；

2）当近壁面网格满足Y＋＞11．225时，采用各种k－ε模型计算结果与实验值均符合良好且差别不大；

3）采用k－ε模型时，非平衡壁面函数的计算精度较其他壁面函数的低；

4）当Y＋≈1时，增强壁面函数不能获得满意的预测计算结果；

5）采用k－ω模型时，近壁面网格Y＋≈1时预测精度较差；

6）对于过冷沸腾数值模拟，k－ε模型的预测精度整体高于k－ω模型。

参考文献：

［1］ KURUL N，PODOWSKI M Z．On the modeling of multidimensional effects in boiling channels［C］∥Proceedings of the 27th National Heat Transfer Conference．Minneapolis，USA：［s．n．］，1991．

［2］ OSE Y，KUNUGI T．Development of a boiling and condensation model on subcooled boiling phenomena［J］．Energy Procedia，2011，9：605－618．

［3］ KREPPER E，KONCˇAR B，EGOROV Y．CFD modelling of subcooled boiling：Concept，validation and application to fuel assembly design［J］．Nuclear Engineering and Design，2007，237（7）：716－731．

［4］ YAN J，SMITH L D，KAROUTAS Z．Departure from nucleate boiling modeling development for PWR fuel［C］∥2013 21st International Conference on Nuclear Engineering．Chengdu：［s．n．］，2013．

［5］ BARTOLEMEI G G，CHANTURIYA V M．Experimental study of true void fraction when boiling subcooled water in vertical tubes［J］．Thermal Engineering，1967，14（2）：123－128．

［6］ BARTOLEMEI G G，BRANTOV V G，MOLOCHNIKOV Y S，et al．An experimental investigation of true volumetric vapor content with subcooled boiling in tubes［J］．Thermal Engineering，1982，29（3）：132－135．

［7］ ISHII M，ZUBER N．Drag coefficient and relative velocity in bubbly，droplet or particulate flows［J］．AIChE Journal，1979，25（5）：843－855．

［8］ MORAGA F，BONETTO F，LAHEY R．Lateral forces on spheres in turbulent uniform shear flow［J］．International Journal of Multiphase Flow，1999，25（6）：1 321－1 372．

［9］ ANTAL S，LAHEY R，Jr，FLAHERTY J．Analysis of phase distribution in fully developed laminar bubbly two－phase flow［J］．International Journal of Multiphase Flow，1991，17（5）：635－652．

［10］RANZ W，MARSHALL W．Evaporation from drops［J］．Chem Eng Prog，1952，48（3）：141－146．

［11］LAUNDER B E，SPALDING D B．Lectures in mathematical models of turbulence［M］．France：Academic Press，1972．

［12］SHIH T H，LIOU W W，SHABBIR A，et al．A newk－εeddy viscosity model for high reynolds number turbulent flows［J］．Computers ＆Fluids，1995，24（3）：227－238．

［13］ORSZAG S A，YAKHOT V，FLANNERY W S，et al．Renormalization group modeling and turbulence simulations［J］．Near－wall Turbulent Flows，1993：1 031－1 046．

［14］ANSYS FLUENT theory guide［M］．Canonsburg，USA：ANSYS Inc．，2011．

［15］MENTER F R．Two－equation eddy－viscosity turbulence models for engineering applications［J］．AIAA Journal，1994，32（8）：1 598－1 605．

［16］MENTER F，KUNTZ M，LANGTRY R．Ten years of industrial experience with the SST turbulence model［J］．Turbulence，Heat and Mass Transfer，2003，4：625－632．

［17］WILCOX D C．Turbulence modeling for CFD ［M］．CA：DCW Industries La Canada，1998．

［18］ANSYS FLUENT tutorial guide［M］．Canonsburg，USA：ANSYS Inc．，2011．

Effect of Turbulence Model on Prediction of Subcooled Boiling

ZHANG Rui1，2，TIAN Wen－xi1，2，CONG Teng－long1，2，SU Guang－hui1，2，QIU Sui－zheng1，2
（1．State Key Laboratory of Multiphase Flow in Power Engineering，Xi’an Jiaotong University，Xi’an710049，China；2．School of Nuclear Science and Technology，Xi’an Jiaotong University，Xi’an710049，China）

Abstract：The ANSYS FLUENT14．5was employed to simulate the two－phase forced convection subcooled boiling in a vertical heated pipe．Different turbulence models，wall functions and two－phase turbulence treatments were utilized to each of four sets of grids．Results were compared with experiment data to analyze the effects of grids，turbulence model，wall functions and two－phase turbulence treatments．The results of k－εmodel show better agreement with experiment data than those of k－ωmodel．Dispersed and Per Phase treatments for two－phase turbulence parameters perform no better than Mixture treatment．The enhanced wall function can’t achieve satisfactory results for the grid with near wall function Y＋≈1．

Key words：subcooled boiling；turbulence model；calculation accuracy

作者简介：张 蕊（1990—），女，陕西咸阳人，博士研究生，核能科学与工程专业

基金项目：国家杰出青年科学基金资助项目（11125522）

收稿日期：2014－04－21；修回日期：2014－06－07

doi：10．7538／yzk．2015．49．08．1366

文章编号：1000－6931（2015）08－1366－08

文献标志码：A

中图分类号：TL334