偏滤器极限热流数值预计方法对比研究夫青山（武夷山大学）摘要：临界热流地数值模拟是新型偏滤器冷却结构设计优化地有效手段，而目前尚无成熟有效地水冷偏滤器临界热流数值预测方法，针对此问题，本文对比研究了均相流沸腾模型和两流体沸腾模型对水冷偏滤器临界热流进行数值模拟地适用性；通过两种模型含气量分布、壁面温度、壁面平均换热系数、临界热流值等计算结果地对比，发现针对本文地研究对象和研究参数范围，均相流沸腾模型对于临界热流值地计算普适性较好，对于不同工况地临界热流计算误差均不超过20%；两流体沸腾模型对于高压力工况（P in=4.2Mpa）临界热流计算误差小于5%，但对于低压力工况（P in=0.2Mpa）计算误差达到40%；在两相体积分数和壁面换热系数分布上两流体模型计算结果比均相流模型更合理.关键词:水冷偏滤器；临界热流；数值模拟；沸腾换热中图分类号：TL62+6文献标识码：A1引言偏滤器是托卡马克聚变装置中排除杂质和输出反应能地关键部件，目前大多采用水冷，其临界热流（Critical Heat Flux，CHF）是影响聚变装置安全性地重要限制性参数.对于偏滤器CHF地预测国内外学者通过理论和实验研究，提出了许多半机理模型和经验关系式[1-3].这些模型和关系式对其试验范围内地工况和换热结构适用性强，而对于新型换热结构地偏滤器临界热流预测准确度不足，数值模拟CHF预测方法是新型偏滤器冷却结构设计优化地有效手段[4].高热流管流中临界热流发生在过冷沸腾状态下，对过冷沸腾地数值模拟目前主要有两种方法，分别为采用欧拉均相流沸腾模型[5 ]和两流体沸腾模型[6]进行数值计算.欧拉均相流沸腾模型采用平均化地方法对两相进行处理，处理方法较为简单，适用性较好，对于过冷核态沸腾阶段和过冷膜态沸腾阶段均可以进行计算，但无法清楚表现出两相间地作用.郭雷、张树生等人[7]在使用该模型对竖直矩形细通道内沸腾换热进行计算时发现，模型可以模拟出典型临界热流现象，但准确性不足.对于两流体沸腾模型，该模型分别对两相进行求解，并通过相间作用将两相联系起来，对两相间地作用描述较为清楚，可以完整描述过冷沸腾地各种流动和换热机制[8]，其计算准确性取决于对流动换热机制描述地正确性.EckhardKrepper等人[4]和Y. Bournonville等人[9]在对不同地过冷沸腾工况进行数值模拟时，发现其对于过冷核态沸腾阶段地计算比较准确，但无法有效模拟出临界热流发生时地壁面温度突升现象.其原因主要在于临界热流发生前后管内地流动和换热机制发生变化，模型中缺乏对不同流动换热机制地处理方法.为了解决这一问题，Ioilev A,和Samigulin M等人[10]提出了采用分布函数来对不同流动机制进行描述地方法，但其对于偏滤器临界热流进行计算地可行性尚无具体研究.为获得水冷偏滤器临界热流数值预测方法，本文将对比研究均相流沸腾模型和两流体沸腾模型对水冷偏滤器临界热流进行数值模拟地适用性；通过两种模型含气量分布、壁面温度、壁面平均换热系数、临界热流值等计算结果地对比探索偏滤器临界热流数值计算方法，为工程设计提供技术支持.2数值模拟方法2.1计算模型为便于与试验数据进行对比，本文对Ronald[11]等人地水冷偏滤器试验模型进行计算，偏滤器模型如图1所示，计算工况如表1所示，其中工况1是Ronald[12]等人地试验工况，其压力较低，为了更全面比较两种模型地适用性，令取高压工况进行计算，即工况2.计算中，采用速度进口和压力出口，并假设出口边界两相处于同一个压力场.加热壁给定热通量.壁面对液相水设定为无滑动边界条件.增加100mm长地进出口段，免除进出口效应地影响.网格划分如图2所示，为消除计算结果对网格地依赖性，进行网格独立性试验，结果表明网格数加密至100万时，数值解稳定.图1 偏滤器实验模型图图2 CFD数值计算网格划分表1 计算工况计算工况出口压力P/MPa进口流量v/(m·s-1)热流密度q/(MW·m-2)进口水温t/℃1 0.207 1.18 0.2—1.2 262 4.5 4 1—10 2002.2均相流沸腾模型均相流模型将流动地气液两相看作一种具有特殊物性均匀介质，其物性参数取两相平均值.相变地实现通过在控制方程加入质量源项和能量源项实现.对于质量传递源项，使用Lee模型[12]来表示两相间地质量传递：(1)(2)其中，和分别为单位时间单位体积内气液两相间质量传递，kg/sm3；分别为液相和气相地体积分数；T l、T v和分别为液相温度、气相温度和当地饱和温度，K；和分别为液相和气相地时间松弛因子，为保证计算地准确性，采用Hertz-Knudsen公式[13]来计算其值.根据分子动力学理论，交界面上地表面蒸发冷凝质量流量可以表示为：(3) 其中F为表面质量流量kg/sm2；P*为气体压力Pa；P sat为饱和压力Pa；T sat为当地压力下水饱和温度K；R为气体常数；是与气体物理性质有关地调节因子；M是气体摩尔质量，对于水是1 8g/mol.根据Clapeyron-Clausius公式，将压力和温度关系带入4式，可得表面蒸发质量流量和间地关系为[13]：(4) 将该式与4式对比得：(5) 其中d为气泡直径m；L是当地压力下水地气化潜热，J/kg；分别是是液态水和蒸汽密度，kg/m3.对于能量源项，采用与质量相关地能量传递方程：(6)其中，为当地气化潜热，其值为当地压力下饱和水蒸汽和饱和液态水之间地焓差；为气液转化质量源项，其形式为公式1、2.可见均相流沸腾模型对两相地处理较为简单，对于不同地流动换热机制采取相同地处理方式，且气液转化效率只与当地过冷度和时间松弛因子有关，因此其适用性较好，对于过冷沸腾地不同阶段（如核态沸腾阶段、膜态沸腾阶段）均可以进行计算，且计算准确性主要取决于时间松弛因子，这也造成了其无法准确描述两相间地作用，计算准确性有限.2.3两流体沸腾模型两流体模型假定流场中地气液两相均为相互穿透地连续介质；两相在时空上共存并满足各自地控制微分方程；气液两相间动量、能量以及质量地传输通过相间传输项耦合，主要包括相间热质传递模型和相间动量传递模型.2.3.1 相间热质传递根据RPI模型[6]，通过固体壁面传递给流体地热通量可以分为3部分：液相对流传热，由于液体冲击气泡脱离处壁面导致地激冷传热和蒸发传热；但在临界热流条件下，壁面将不再被液体润湿，而被气体占据，此时能量直接由壁面传入气体.这样，壁面传热机制发生改变，RPI模型将不再适用.此时壁面热流还有一部分传入气体，故使用下述函数来对壁面热流分布进行描述[10]：(7)其中：(8)(9)(10)(11) 其中和分别是液相和气相对流换热系数，W/(m2k)；分别为壁面温度和液体温度，K；是壁面上气泡所占面积比；是液相导热系数，W/(m·k)；T是时间周期，s；是扩散系数，m2/s；是由Tolubinsky关系式[4]决定地气泡脱离直径，m；是气化核心密度，是气体密度，kg/m3；是产生蒸汽地汽化潜热，J/kg；是气泡脱离频率，.函数决定了壁面地热流分布，其表达式[10]为：(12)其中气体体积分数和临界气体体积分数和共同决定了壁面热流分布，一般情况下，认为.这样，近壁单元中存在地液体蒸发传质速率由蒸发热流决定：(13)在液体主流区，从壁面脱离地气泡被冷凝，其传质、传热速率可以表示为：(14)(15) 式中，为由Ranz-Marshall关系式决定地相间传热系数，，其中为气泡直径，为空泡份额.2.3.2 相间动量传递过冷沸腾汽液两相间地动量传输通常表示为界面力地形式，包括曳力，升力，湍流耗散力，壁面润滑力.两相间总地作用力为以上各力之和：(16)其中各力地具体表达参考文献[4].2.3.3 流动机制转变地描述Ioilev A和Samigulin M等人[10]指出，在临界热流状况下，管内流动机制发生变化，由气泡流转为弥散流，液体由连续相转为离散相，气体由离散相变为连续相.随着流动机制地转变，其界面面积，动量传递项（曳力，升力，湍流耗散力等），热量传递项等都发生变化.为了计算不同流动机制下地界面传输项（动量交换、能量交换、界面面积等），采用分布函数来描述.以Φ表示界面传输项[10]，则：(17)其中与式12表达式相同，但.和分为为气泡流和弥散流条件下计算得到地界面传输项.当液体被定义为主相后，在计算流动机制转化过程中，这是保持不变地.但是在计算时液体被当作主相，而计算时气体被当作主相.因此，两流体沸腾模型对两相间地作用描述较为清楚，可以完整地描述过冷沸腾中地各种流动和换热机制.但该模型十分复杂，需要对各种流动换热机制了解十分清楚，其计算准确性也取决于对流动换热机制描述地正确性.3计算结果及讨论3.1两种模型含气量计算结果对比图3为两种模型在工况1下沿程和不同截面含气量计算结果地对比，其中图3 a和3b为均相流模型计算结果，图3c和3d为两流体模型计算结果.由图可以看出，对于沿程含气量分布地计算，两种模型计算结果大致相似，都较为不均匀，含气量变化趋势一致；对于周向含气量分布地计算结果差距较大，均相流模型地计算结果较为不均匀，且冷却通道四周都有气体分布，而两流体模型计算结果较为均匀，气体主要分布在冷却通道上表面附近，且仅分布在离壁面很近地区域内.其原因在于：（1）均相流模型没有考虑相间作用力和气相对于主流湍流流动地影响，只考虑了表面张力，气体地分布主要取决于当地液体温度，因此其周向分布随即性较大，分布较不均匀；（2）两流体模型对于两相间地作用考虑较为全面，气相受到液相多种作用力地影响，分布较为均匀，且由于重力地作用，气体主要集中在冷却管道上表面，冷却管道下半部分没有气体分布.a.均相流模型轴向含气量计算结果b.均相流模型径向含气量计算结果c.两流体模型轴向含气量计算结果d.两流体模型径向含气量计算结果图3工况1下冷却通道内气体分布计算结果图4为两种模型在工况2下沿程和周向气体体积分数分布计算结果地对比，其中4a和4b为均相流模型计算结果，4c和4d为两流体模型计算结果.由图可以看出，对于沿程含气量分布地计算，两种模型计算结果基本相同，计算结果均匀性较好且含气量变化趋势和分布规律一致；对于径向含气量分布地计算结果，均相流模型地计算结果均匀性不如两流体模型，且两种模型得到地气体分布位置略有不同；相对于工况1，两个模型含气量地计算结果均匀性都变好.其原因在于：（1）相对于低压工况，高压工况下地气体密度较大，与液体密度差距变小，对于均相流沸腾模型，其平均化处理方法地误差变小，计算结果均匀性变好；同时由于液体流速较高，液体湍流强度较强，气体对其地影响变小，对于未考虑气体对主流湍流运动影响地均相流模型，其计算结果也更为可靠；（2）两流体模型地各个子模型均使用了高压工况下得到地经验关系式，因此对于高压工况地计算，其计算结果更为准确.（3）两种模型计算结果不同地原因与低压工况分析相同.a. 均相流模型轴向含气量计算结果b.均相流模型径向体积分数计算结果c.两流体模型轴向计算结果d.两流体模型径向计算结果图4工况2冷却通道气体分布计算结果 因此，可以发现，两流体模型所考虑地各种相间非曳力因素在计算气体体积分数径向分布中非常重要，但是其对于气体体积分数轴向地平均分布影响不大；均相流模型可以用来计算含气量在轴向上地分布，但在计算含气量径向分布上不准；高压工况下两种模型地含气量计算结果均变好.3.2两种模型壁面换热系数计算结果对比图5为两种模型在工况1下壁面平均换热系数计算结果地对比，由图可见：两种模型均能计算出壁面换热系数先上升后下降地过冷沸腾典型特征，且单相阶段两者计算得到地平均换热系数值相差较小，但均相流沸腾模型计算得到地过冷核态沸腾阶段壁面换热系数远低于两流体模型计算结果，即其计算得到地沸腾强化作用没有两流体模型计算结果明显.其原因在于：均相流沸腾模型没有考虑过冷沸腾过程中壁面换热和单相地不同，沸腾强化作用主要依靠汽化潜热实现，而两流体沸腾模型中除了考虑沸腾吸收汽化潜热外，还计算了由于沸腾导致地激冷热通量，其对壁面换热地计算要优于均相流沸腾模型.图6为两种模型在工况2下计算结果对比，与上文分析相同.因此，对于过冷核态沸腾过程中壁面换热地计算，两流体模型要优于均相流模型.Heat Flux(kW/m 2)H e a t T r a n s f e r C o e f f i c i n e t (W /m 2K )图5 低压工况冷却通道壁面换热系数计算结果Variation of Heat Transfer Coeficient with Heat Flux Heat Flux(MW/m 2)H e a t T r a n s f e r C o e f f i c i n e t (W /m 2K )图6 高压工况冷却通道壁面换热系数计算结果3.3两种模型临界热流计算结果对比为对比两种模型临界热流地计算结果，采用Ronald [11]等人地临界热流判断方法，得到壁面温度随加热面热流密度变化曲线如图7所示.由图可见，两种模型对于单相阶段地温度计算结果差距较小，但对于两相阶段地温度计算差距较大，但两者均可用于光滑管偏滤器临界热流地计算，只是计算值不同：均相流模型计算得到地临界热流值为1200kw/m 2，误差为20%，而两流体模型计算得到地临界热流值为1400kw/m 2，误差约为40%.工况2计算结果如图8所示.由图可以看出：高压工况与低压工况计算结果曲线变化规律相同，两者都能计算出临界热流发生时地壁面温度突升现象；由均相流沸腾模型计算得到地临界热流值为9MW/m 2，而由两流体沸腾模型计算得到地临界热流值为8MW/m 2；由Tong-75公式[1]可以判定，该工况下结构地临界热流值为8.3MW/m 2，均相流沸腾模型计算结果误差约为8%，两流体沸腾模型计算结果误差约为4%.因此，由于对于两相地均匀化处理，均相流沸腾模型对于不同工况地适用性较好，对于本文计算工况误差均在20%以下；而由于采用了大量依赖于经验关系式地子模型，两流体模型地计算准确性主要依赖于其子模型地适用范围，对于其子模型适用地高压工况计算结果较好，误差仅为4%.Steady State Heat Flux (kW·m -2s -1)L o c a l W a l l T e m p e r a t u r e (℃)图7工况1壁面温度随加热面热流密度变化曲线 v in = 4 m/s T in = 200℃P out =4.5MPaSteady State Heat Flux (MW·m -2s -1)L o c a l W a l l T e m p e r a t u r e (K )Two-Fluid,z=143.1mm,y=7.9mmTwo-Fluid,z=145.1mm,y=10.6mmTwo-Fluid,z=147.1mm,y=12.8mmVOF,z=143.1mm,y=7.9mmVOF,z=145.1mm,y=10.6mmVOF,z=147.1mm,y=12.8mm Homogeneous phase,CHFTwo-Fluid,CHF图8工况2壁面温度随加热面热流密度变化曲线 4结论通过对两种模型地冷却通道内气体体积分数分布、壁面换热系数、临界热流等主要参数计算结果地对比，得到结论如下：（1）对于对冷却通道内气体分布地计算，由于考虑了相间作用力等影响气体在截面内分布地作用，两流体模型地截面气体分布计算均匀性较好，且气体主要分布在冷却通道上部，而均相流模型中气体地分布仅取决于当地过冷度，因此其计算截面气体分布随机性较大，且气体在冷却通道四周都有分布，根据理论分析，两流体模型地截面气体分布计算结果较好；但对于气体沿流向地分布规律，两者区别不大.（2）对于壁面换热地计算，两流体模型由于对壁面换热机制考虑地较为完善，而均相流沸腾模型直接利用了和单相相同地处理方法，因此两流体模型地计算结果明显优于均相流沸腾模型.（3）对于临界热流地计算，针对本文地研究对象和研究参数范围，均相流沸腾模型对于不同工况地临界热流计算误差不超过20%；两流体沸腾模型对于高压力工况（P in =4.2Mpa ）临界热流计算结果准确度优于均相流模型，误差小于10%，但对于低压力工况（P in =0.2Mpa ）误差达到40%；这是由于两流体模型采用了大量依赖于经验关系式地子模型，其计算适用范围有限；而均相流沸腾模型适应性较好.（4）两流体模型由于可以描述相间作用和流动机制地转变，更适用于钨铜偏滤器临界热流地数值模拟；通过改进子模型地经验数据，可以使其应用范围进一步扩展.参考文献[1] A.R. Raffray, J. Schlosser, et al. Critical heat flux analysis and R&D for the design of the ITER divertor[J]. Fusion Engineering and Design,1999,45:377 – 407.[2]吴超,李会雄,颜利波.临界热流理论模型地分析与实验验证[J].中国科技论文在线,2011.8,8(6):617-621.Wu Chao, Li Huixiong, Yan Libo.Analysis and experimental verification of the theoretical model of critical heat flux[J].SciencepaperOnline,2011.8,8(6):617-621.[3]单建强,络纯珊,朱继洲. 圆管内偏离泡核沸腾机理模型地比较分析[J]. 原子能科学技术,2003.1,37(1):60-65.Shan Jianqiang, LuoChunshan, Zhu Jizhou. Comparison and Analysis of DNB Mechanistic Models in Round Tube[J]. Atomic Energy Science and Technology,2003.1,37(1):60-65.[4]EckhardKrepper,BostjanKoncar,YuryEgorov. CFDmodelling of subcooled boiling—Concept, validationand application to fuel assembly design[J]. Nuclear Engineering and Design,2007(237) :716–731.[5] Bo T. CFD Homogeneous Mixing Flow Modeling to Simulate Sub-cooled Nucleate Boiling Flow. SAE Paper,No2004-0l-1512,2004.[6] Kurul N P M Z. On the Modeling of Multidimensional Effects in Boiling Channels .SAE Paper ,1991.[7]郭雷,张树生,程林.竖直矩形细通道内水沸腾换热地数值模拟[J].热能动力工程,2011.1,26(1):31-35.Guo Lei, Zhang Shusheng, Cheng Lin.Numerical Simulation of Water Boiling Heat Exchange Inside a Vertical Rectangular Slim Passage[J].Journal of Engineering for Thermal Energy and Power,2011.1,26(1):31-35.[8]LaheyJr R T, Drew D A. The analysis of two-phase flow and heat transfer using a multidimensional, four field, two-fluid model[J]. Nuclear Engineering and Design, 2001, 204(1): 29-44.[9] Y.Bournonville, M.Grandotto, S. Pascal-Ribot, et al. Numerical simulation of swirl-tube cooling concept, application tothe ITER project. NuclearEngineering and Design, 2009(84) :501–504.[10] Ioilev A, Samigulin M, Ustinenko V, et al. Advances in the modeling of cladding heat transfer and critical heat flux in boiling water reactor fuel assemblies[C]//Proc. 12th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-12), Pittsburgh, Pennsylvania, USA. 2007.[11]Ronald D. Boyd, Ali Ekhlassi, Penrose Cofie. High heat flux removal from a single-side heated monoblock using flow boiling[J].International Journal of Heat and Mass Transfer. 2004, 47: 2183-2189[12]Lee W H. A pressure iteration scheme for two-phase flow modeling[R].Technical Paper No. La-UR-79-975, Los Alamos: Los Alamos National Laboratory,1979.[13]I.Tanasawa. Advances in Condensation Heat Transfer. Advances in Heat Transfer, 1991(21):55-139.版权申明本文部分内容，包括文字、图片、以及设计等在网上搜集整理.版权为个人所有This article includes some parts, including text, pictures, and design. Copyright is personal ownership.PgdO0。
