低流速下燃料包壳表面污垢沉积的传热计算分析*

张小英　尹盛晖　卢向晖　林维青　王婷　陈焕栋

(1.华南理工大学 电力学院， 广东 广州 510640； 2.中科华核电技术研究院， 广东 深圳 518026)

低流速下燃料包壳表面污垢沉积的传热计算分析*

张小英1尹盛晖1卢向晖2†林维青2王婷2陈焕栋1

(1.华南理工大学 电力学院， 广东 广州 510640； 2.中科华核电技术研究院， 广东 深圳 518026)

摘要：针对核电站发生破口失水事故后碎片在燃料包壳外表面形成污垢沉积的现象,采用一维模型建立了燃料包壳和周围冷却剂温度场的数值计算模型,与流体动力学分析软件CFX模拟结果进行对比,分析污垢厚度、污垢导热系数以及冷却剂流速对包壳温升的影响,并分析包壳温升瞬态变化.经计算发现，此计算模型计算结果与CFX模拟结果相近,最大相对偏差为2%.经分析发现污垢的产生导致包壳外表面温度突升,且温度升高的幅度随着污垢厚度增大而增大,随着污垢导热系数的减小而增大.存在一个临界冷却剂入口流速,当入口流速低于临界值时,堆芯内将产生沸腾危机.

关键词：燃料棒；污垢增长；沸腾传热流动；流量；临界热流密度

核电站发生破口失水事故后,为了保证堆芯连续的长期的冷却,就需要应急堆芯冷却系统(ECCS)能够长时间运行.当换料水箱到达低水位后,ECCS系统将由安全壳地坑取水.但是管道破口产生的碎片随着泄露的冷却剂或喷淋液迁移到安全壳地坑滤网，造成滤网堵塞.同时，碎片在较窄的燃料组件流道内会由于碎片堆积而堵塞流道,造成燃料包壳表面污垢沉积.燃料包壳表面污垢沉积会造成燃料棒传热恶化，可能威胁到燃料棒完整性,因此,深入研究棒状燃料组件在发生污垢沉积时的热工水力特性和安全性,对我国压水堆核电站的安全运行具有重大意义.

在国外,Schikorr等[1]运用修正的Rehme/Dalle Donne关系式计算了带定位格架的棒束间流场压降.Anglart 等[2]对沸水堆燃料组件进行了流动传热计算,提出了冷却剂回路的流动和传热的多维两相流物理模型.Holloway等[3]提出了计算带定位格架的燃料组件的单相对流换热关系式.Yang等[4]用CFX软件对超临界压水堆棒束通道的冷却剂的流动和换热特性进行了三维数值模拟.Sreenivas等[5]运用三区域一维模型提出了高质量的CHF预测关系式.Mohanty等[6]运用K-1模型对子通道内湍流流动与换热进行了数值模拟.在国内,罗磊等[7]对运行条件下燃料元件及其周围冷却剂温度场和流场进行了三维数值模拟.刘天才等[8]运用Relap5程序对CARR一盒燃料元件入口全部堵塞时堆功率、压力、流量和温度的变化进行了分析.张小英等[9]采用CFX软件对带定位格架的5×5燃料组件内流场进行了三维数值分析.

综上所述,对于燃料棒的热工问题,考虑低流速下且燃料棒包壳表面存在污垢沉积时的传热研究较少.文中将对低流速下燃料棒包壳表面存在污垢沉积时的传热问题进行研究,建立数学分析模型并计算,获得瞬态和稳态换热情况下燃料包壳外表面的温度分布,分析污垢厚度和污垢导热系数以及入口冷却剂流量对包壳温升的影响,并且采用CFX软件对燃料棒的流动与传热进行模拟,与数学分析结果相互验证.

1燃料包壳污垢沉积的几何模型

文中计算考虑全高度的单根燃料棒,沿燃料棒高度有6个宽格架和3个窄格架,格架形状为正方形,窄格架的高度为宽格架的一半.假设污垢均匀地覆盖在燃料包壳后面部分,约1/3高度.考虑几何的对称性,选取燃料棒及流体通道横截面的1/4为研究对象,假定格架为贴覆在包壳/污垢层外表的指环状结构.由此,污垢沉积部分的燃料棒几何模型简化如图1所示.图中,包壳内径为8.36 mm,包壳外径为9.5 mm,燃料棒中心距为12.6 mm,包壳外污垢层的厚度分别为0、δ、2δ、3δ、4δ和6δ共6组数值,其中δ=0.225 mm.

图1　燃料棒横截面示意图Fig.1　Cross section diagram of fuel rod

2流场与传热计算数学模型

针对文中研究的问题,建立流体的质量、能量和动量守恒方程,计算沿燃料棒高度流体的焓、速度、温度和压力.然后判断通道的流动形态,计算燃料棒与流体的对流换热系数以及燃料棒外表的温度.

2.1　流场的守恒方程

采用均相流模型计算沿通道高度流体参数随时间的变化,根据均相流模型可列出流体的能量、质量和动量守恒方程:

(1)

(2)

(3)

采用时间的向后差分和空间的迎风差分,得到流场守恒方程的差分格式:

(4)

(5)

(6)

式中，i表示时间节点,n表示位置节点.

求解上述差分方程时,为确定冷却剂物性参数,需首先假设一组流场压力,然后采用能量守恒方程计算沿燃料棒高度的冷却剂焓,再求解质量守恒方程计算冷却剂的流速,最后根据动量守恒方程计算冷却剂压力,完成一次迭代.采用新的流场压力重新确定冷却剂物性参数,并开始下一次迭代计算,直至计算得出的出口压力等于所给定的出口压力时结束迭代.

2.2　流型判断

在文中冷却剂以低速垂直向上的流动过程中,冷却剂将经历单相过冷区、过冷沸腾区、饱和沸腾区、过度沸腾区,甚至发展为单相蒸汽区.这些区域流体与壁面间的对流换热强度有很大不同,适用的对流换热系数计算关系式也不同.针对文中的问题,需要计算的流型转换点有：泡核沸腾起始点(ONB)、气泡脱离壁面点(FDB)、饱和沸腾点(SAT)、临界热流密度点(CHF)和蒸干点(DRY).

泡核沸腾起始点是流动自单相向两相的转变点,该点既具有单相对流换热的特征,又具有泡核沸腾换热的特征,因此可以采用Jens-Lottes公式[10]计算.对于气泡脱离壁面点的计算采用Saha-Zuber模型[11].对于热平衡饱和沸腾起始点,文中采用通道中流体比焓达到当地压力下的饱和液体比焓所在的位置.

临界热流密度点的计算采用Kandikar[12]提出的低压低流量工况下CHF关系式确定:

(7)

式中,qCHF为对应节点处的临界热流密度，hfg为冷却剂汽化潜热，θR为气泡与壁面接触角，ρV为气相密度，ρL为液相密度，ρm为混合相密度，Dh为燃料组件内通道当量直径，σ为表面张力，x为含气率，μL为液相动力黏度.

式(7)适用于G=23~53 000 kg/(m2·s),p=0.047~1.55 MPa.蒸干点视为在含气率xg>0.99时达到.

2.3　对流换热系数的计算

从通道入口至FDB点之间的传热为冷却剂单相对流换热,由于冷却剂流速很低,热流密度很高,若流动处于层流状态,需考虑自然对流的影响,因此文中采用米海耶夫[10]推荐的关系式计算该区域的对流换热系数:

(8)

式中,Cp为比定压热容,g是重力加速度,k为冷却剂导热系数,αv是液体的体积膨胀系数,ΔT=Tw-Tf,Prf是以流体平均温度Tf作为定性温度计算的普朗特数,Prw则是以壁面温度Tw作为定性温度计算的普朗特数,μf为冷却剂动力黏度.

若流动处于湍流状态,文中采用Dittus-Boelter公式计算FDB点之前的换热系数.

从气泡脱离点FDB至热平衡饱和沸腾起始点SAT之间的传热为欠热泡核沸腾,对该区域的对流换热系数,文中采用适用于低流速下欠热泡核沸腾换热的扩展Chen公式[13]计算:

(9)

式中,h为对流换热系数，hl为单相液体强迫对流换热系数,hnb为自然对流换热系数.

式(9)考虑单相液体强迫对流与自然对流两种作用.

流体通道的热平衡饱和沸腾起始点SAT至CHF点之间为饱和泡核沸腾传热区,该区域的对流换热系数采用饱和沸腾区的Chen公式[10]计算,考虑在饱和泡核沸腾区发生着泡核沸腾和强制对流两种传热机理,并且这两种机理对传热的作用可以叠加,因此

h=hlF(Xtt)+hnbSRetp

(10)

式中,Xtt为马蒂内里参数，Retp为雷诺数.

在CHF点之后至液体蒸干,文中用李虹波等[14]推导出的公式计算两相区传热:

(11)

在冷却剂中的液滴完全蒸干后,通道中流体将变成单相蒸汽的强制对流.该区域的对流换热系数采用Dittus-Boelter公式计算.

2.4　燃料包壳外表面温度计算

根据前面计算的流场温度和对流换热系数,可以采用牛顿对流换热公式计算燃料棒外壁温度:

(12)

燃料包壳外表面到燃料棒的热量是通过污垢/格架、气隙以及棒芯导热进行的,采用圆柱坐标下的一维导热方程模拟,由此得到

(13)

3流场与传热的CFD计算模型

在CFD分析中,先用UG软件对1/4燃料棒及冷却剂通道进行建模,然后将模型导入ICEM中进行网格划分.所划分的网格为结构化网格,得到的节点数为2 207 520.计算时入口边界条件设置为入口流体的速度为2u0,温度为120 ℃.包壳内表面热流密度取30 000 W/m2.出口类型设置为opening,压力为4×105Pa.两相流模型选取相变模型,湍流模型为k-e模型.求解采用Simple算法,收敛残差为0.000 1.

用最好的大地测量滑动模型作为参考，由KiK－Net－BH数据集导出的模型3（见图4）非常令人满意，并可认为优于由KiKNet－SF数据集导出的模型（图4模型2），后者大大低估了地震的震级。值得注意的是，50m的滑动峰值只有在使用海底GPS数据时才能获得。没有这种数据，由GPS导出的滑动模型（模型4）与KiK－Net－BH数据集得到的模型（模型3）近乎相同，二者波形相似，滑动峰值均小于25m。模型3与岸上GPS数据的拟合与模型1的拟合同等好。

4计算结果与分析

4.1　计算结果与CFD分析结果对比

图2　冷却剂温度与包壳外表面温度结果验证Fig.2　Results verification of coolant temperature and cladding temperature

4.2　低流量燃料包壳传热计算的计算结果

4.2.1最不利工况下瞬态和稳态计算结果分析

图3(a)中示出了冷却剂温度沿燃料棒高度分布的稳态结果,可以看出在0.25 m之前,冷却剂温度随着高度的增加而升高.当冷却剂温度到达饱和温度以后,基本维持不变.图3(b)中示出了冷却剂温度随时间的变化情况.可以看出,在25 s之前冷却剂温度随着时间的推移而增加.在25 s之后,当温度达到饱和液体温度时,冷却剂温度曲线基本不随时间变化,维持在饱和温度.

图4(a)中示出了冷却剂压力的稳态结果.可以看出,在初始阶段,冷却剂密度较大,重位压降较大,冷却剂压力下降速度很快.随着轴向高度的增大,冷却剂密度减小,压力的下降速度也变慢.当到达污垢区域时,流速增大,压力下降速度随之加快.图4(b)中示出了冷却剂压力瞬态计算结果,可以看出,初始时,冷却剂为液体,密度大,重位压降大,所以初始时压降最大.随着时间的增加,密度减小,压降随着时间的增加而减小,并在100 s趋于稳定,最终压降为3 000 Pa左右.

图 3　冷却剂温度计算结果Fig.3　Calculation results of coolant temperature

图 4　冷却剂压力计算结果Fig.4　Calculation results of coolant pressure

图5(a)中示出了包壳外表面温度的稳态结果.可以看出,在初始阶段,包壳外表面温度沿着燃料棒高度方向升高,当冷却剂温度到达对应压力下的饱和温度后,包壳外表面温度也基本不再变化,到达污垢区域时,在包壳和冷却剂之间存在很大的导热热阻,包壳外表面温度升高将近200 K.在格架存在区域,由于格架的导热系数大于污垢的导热系数,故导热热阻小于不存在格架区域的热阻,包壳外表面温度也较无格架区域低.图5(b)中示出了包壳外表面温度的瞬态计算结果,可以看出,包壳外表面温度很快达到稳定状态,并且25 s时刻的温度还高于稳定后的温度,这是因为,在25 s时冷却剂温度虽然低,但对应的对流换热系数也小,根据对流换热公式,燃料棒外表面温度可能反而更高.

图5　包壳外表面温度计算结果Fig.5　Calculation results of cladding temperature

4.2.2污垢导热系数不变时污垢厚度的影响分析

图6(a)中示出了冷却剂温度随污垢厚度的变化情况.从图中可以看出,当污垢厚度增加时,通道内的阻力增大,流体流速减小,导致在入口段的温度上升速率变快.在饱和区,由于冷却剂压力相差不大,各污垢厚度下的冷却剂温度基本一致.图6(b)中示出了传热系数随污垢厚度的变化情况.从图中可以看出,在污垢存在区域,传热系数突然下降,且污垢越厚,传热系数越小.由于污垢的存在,在冷却剂与包壳外表面间存在污垢导热热阻,厚度越大,热阻越大,传热系数越小.图6(c)中示出了冷却剂压力随污垢厚度的变化情况.厚度减小时,阻力减小,流体流速增大,到达饱和区的距离增大,在单相液体区的压力下降速度比饱和区的压力下降速度大,这是因为单相液体区流体的密度大,重位压降大.到达污垢区域后,通道面积减小,流体流速增大,导致压力下降速度增大.图6(d)中示出了燃料包壳外表面温度随污垢厚度的变化情况.从图中可以看出,在污垢存在区域,包壳外表面温度突然升高,污垢越厚,包壳外表面温度越高.污垢越厚,综合换热系数越小,导致包壳外表面温升越大,当污垢厚度为δ时,包壳温升到达200 K左右.当污垢厚度为δ时,由于在定位格架处的污垢厚度比无定位格架处薄,所以包壳温升比无格架处低.

图6　污垢厚度的影响Fig.6　Influence diagram of dirt thickness

4.2.3污垢厚度不变时污垢导热系数的影响分析

在入口流速为2u0、污垢厚6δ时,改变污垢导热系数,分析导热系数对包壳温升的影响.

图7污垢导热系数影响图Fig.7　Influence diagram of dirt thermal conductivity

4.2.4冷却剂入口流速的影响分析

分别计算入口流速为u0、2u0、3u0、4u0,5u0和6u0时最不利工况下各参数的变化情况.

图8(a)中示出了冷却剂温度随入口流速的变化情况.从图中可以看出,在单相液体区,入口流速越大,冷却剂温度上升速度越慢,冷却剂到达饱和区所经过的距离越长.流速为u0时,由于冷却剂流量过小,导致冷却剂的焓升很快,在轴向位置为2.5 m左右处变成单相蒸汽,在单相蒸汽区,温度急剧上升,最终到达960 K.其余5种流速下冷却剂不会出现蒸干的现象,最终温度均为对应压力下的饱和温度.图8(b)中示出了对流换热系数随入口流速的变化情况.从图中可以看出,流速越小,壁面产生气泡的速率越快,流动越快到达泡核沸腾换热区,对流换热系数突然升高的位置越靠前.在污垢产生处之前的一小段距离处,流速为u0的工况下冷却剂会变成单相蒸汽,在此点之前会出现沸腾危机,传热恶化,对流换热系数突然降低.其它5种流速下,不会出现沸腾危机,对流换热系数一直处于很高的状态.在污垢产生处,由于流通截面减小,流速增大,对流换热系数会发生突增.图8(c)中示出了包壳外表面温度随入口流速的变化情况.从图中可以看出,流速为u0的曲线与其它5条曲线差别巨大.流速为u0的工况下,由于液体冷却剂全部蒸干,最终温度很高,导致冷却能力下降,壁面温度急剧升高,在沸腾危机处的温升达到300 K,在出口处的温度超过1 200 K,这在核反应堆的运行当中均是不允许出现的.其他5种流速下,包壳温度在污垢处升高200 K左右,出口处温度在600 K左右.图8(d)中示出了冷却剂压力随入口流速的变化情况.从图中可以看出.入口流速越大冷却剂压降越大.

图8　冷却剂入口流速影响图Fig.8　Influence diagram of inlet coolant velocity

5结论

文中针对燃料包壳外表面存在污垢沉积的现象,对低流速下包壳温度进行了稳态和瞬态计算,分析污垢厚度和导热系数以及冷却剂流量对温升的影响.通过对结果分析得出如下结论:

(1) 在相同的边界条件下,燃料组件通道内冷却剂的温度计算结果与CFX模拟结果最大相对偏差为2%,包壳外表面温度几乎相同，证明了程序模拟方法的可行性.当燃料包壳表面存在污垢时,燃料包壳外表面产生很高的温升,对燃料包壳完整性构成威胁.

(2) 冷却剂温度随着污垢厚度的增加而降低,不随污垢导热系数变化而变化.包壳外表面温度随着污垢厚度的增加而升高,随着污垢导热系数的减小而升高.保持其他边界条件不变,减小冷却剂入口流速时,冷却剂温度上升速度将加快,冷却剂与壁面间的扰动将减小,对流换热系数减小,包壳外表面温度将有少许上升.另外,存在一个临界入口流速,当冷却剂入口流速低于此临界流速时,在燃料组件通道出口之前,冷却剂将出现沸腾危机现象,传热系数突降,传热恶化,包壳外表面温度将快速升高.

参考文献：

[1]Schikorr M,Bubelis E,Mansani L,et al.Proposal for pressure drop prediction for a fuel bundle with grid spacers using Rehme pressure drop correlations [J].Nuclear Engineering and Design,2010,240(7)：1830-1842.

[2]Anglart H,Nylund O,Kurul K,et al.CFD prediction of flow and Phase distribution in fuel assemblies with spacers [J].Nuclear Engineering and Design,1997,177(1)：215-228.

[3]Holloway Mary V,Beasley Donald E,Michael E Conner.Single-phase convective heat transfer in rod bundles [J].Nuclear Engineering and Design,2008,238(4):848-858.

[4]Yang Xing-bo,Su G H,Tian Wen-xi,et al.Numerical study on flow and heat transfer characteristics in the rod bundle channels under super critical pressure condition [J].Annals of Nuclear Energy,2010,37(12):1723-1734.

[5]Sreenivas Jayanti,Michel Valette.Calculation of dry out and post dry out heat transfer in rod bundles using a three field model [J].International Journal of Heat and Mass Transfer,2005,48(9):1825-1839.

[6]Mohanty A K,Sahoo K M.Turbulent flow and heat transfer in rod-bundle subchannels [J].Nuclear Engineering and Design,1988,106(3):327-344.

[7]罗磊,陈文阵,陈志云,等.单个燃料元件热工水力三维数值模拟 [J].海军工程大学学报,2012,23(1):63-72.

Luo Lei,Chen Wen-zhen,Chen Zhi-yun,et al.Numerical simulation of thermal hydrodynamic of single reactor fuel rod [J].Journal of Naval University of Engineering,2012,23(1):63-72.

[8]刘天才,金华晋,袁履正.中国先进研究堆堵流事故分析 [J].核动力工程,2006,27(5):32-44.

Liu Tian-cai,Jin Hua-jin,Yuan Lü-zheng.Flow blocage accident analysis for china advanced research reactor [J].Nuclear Power Engineering,2006,27(5):32-44.

[9]张小英,孙庆友,乔磊,等.全结构的5×5定位格架及棒束通道的三维流场分析 [J].华南理工大学学报:自然科学版,2014,42(12):104-110.

Zhang Xiao-ying,Sun Qing-you,Qiao Lei,et al.Analysis of 3B flow in unabridged 5×5 spacer grid and rod bundle channel [J].Journal of South China University of Technology:Natural Science Edition,2014,42(12):104-110.

[10]郝老迷.核反应堆热工水力学 [M].1版.北京：原子能出版社,2010：36-58.

[11]于平安,朱瑞安,喻真烷,等.核反应堆热工分析 [M].上海：上海交通大学出版社,1979：36-43.

[12]Kandlikar Satish G.A scale analysis based theoretical force balance model or critical heat flux (CHF) during saturated flow boiling in icrochannels and minichannels [J].Journal of Heat Transfer,2010,132(8):1-13.

[13]任泽霈.对流换热 [M].1版.北京：高等教育出版社,1998：243-246.

[14]李虹波,陈炳德,赵华,等.DNB后膜态沸腾传热计算模型 [J].核科学与工程,2011,31(2):146-153.

Li Hong-bo,Chen Bing-de,Zhao Hua,et al.Numerical model of post-DNB film boiling heat transfer [J].Chinese Journal of Nuclear Science and Engineering,2011,31(2):243-246.

Calculation and Analysis of Heat Transfer Characteristics of Fuel Cladding with Dirt at a Low Flow Rate of Coolant

ZhangXiao-ying1YinSheng-hui1LuXiang-hui2LinWei-qing2WangTing2ChenHuan-dong1

(1.School of Electric Power, South China University of Technology, Guangzhou 510640, Guangdong, China;

2.China Nuclear Power Technology Research Institute, Shenzhen 518026, Guangdong, China)

Abstract:Aiming at the fragments of deposit on the fuel cladding caused by the loss of coolant accident of nuclear power stations, a numerical calculation model of the temperature field on the fuel cladding and the coolant is constructed by using the one-dimensional model. Then, the numerical calculation results are compared with the CFX simulation ones, so as to analyze the influences of the thickness of the dirt, the thermal conductivity of the dirt and the flow rate of the coolant on the temperature rise of the fuel cladding and discuss the transient change of the temperature rise. Calculation results show that the numerical results are close to the CFX simulation ones, with a maximum relative deviation of 2%.Through an analysis, it is also found that the dirt deposition leads to a temperature rise on the fuel cladding, and the extent of temperature rise increases with the increase of the dirt thickness as well as with the decrease of the thermal conductivity of the dirt. In addition, it is found that there exists a critical inlet flow rate of coolant, and a boiling crisis occurs in the reactor core when the inlet flow rate of coolant is below the critical value.

Key words:fuel rods;fouling growth;heat transfer boiling fluid;flow rates;critical heat flux

中图分类号：TL 352

doi:10.3969/j.issn.1000-565X.2015.08.001

文章编号：1000-565X(2015)08-0001-08

作者简介：张小英(1973-)，女，教授，博士生导师，主要从事反应堆热工水力研究.E-mail: zxy1119@scut.edu.cn†通信作者： 卢向晖(1978-)，男，高级工程师，主要从事核安全分析研究.E-mail: luxianghui@cgnpc.com.cn

*基金项目：国家自然科学基金资助项目(51176052，51376065)

收稿日期：2014-10-28

Foundation items: Supported by the National Natural Science Foundation of China(51176052，51376065)