燃料组件偏置及弯曲对沸腾临界影响研究

龚佳奇,丛腾龙,肖 瑶,顾汉洋,丁 铭

(1.哈尔滨工程大学 核科学与技术学院,黑龙江 哈尔滨 150001;2.上海交通大学 机械与动力工程学院,上海 200030)

临界热流密度限制了反应堆堆芯的最大功率密度,与反应堆的安全性、经济性密切相关[1-3],是热工设计中必须要确定的限值[4]。燃料组件在堆内由于制造的偏差、外力的作用,或在临界热流密度热工实验因为安装的偏差,导致燃料棒整体偏置或轴向弯曲,引起冷却剂通道流通截面的几何形状和大小偏离设计值,可能对沸腾临界造成影响。Nagino等[5]对4×4带定位格架棒束组件进行了临界实验,其中轴向和径向功率都不均匀,沿轴向位置布置2个搅混格架及1个位于中间的支撑格架,两搅混格架中间的1根高功率棒发生弯曲,在两搅混格架中点与附近的高功率棒及导向管发生接触(即弯曲程度达到最大)。结果表明棒弯曲在高压下对DNBR的影响很大、低压下没有影响。Tamai等[6]对37稠密栅带定位格架棒束组件进行了临界实验,测试段包含4个定位格架,在靠近出口两格架间发生燃料棒单棒弯曲,弯曲形状呈余弦形,即弯曲段中间达到弯曲最大值,并与相邻棒发生接触。发现弯棒使临界数值比未弯曲棒束降低10%,临界位置发生在棒接触点下游。Markowski等[7]分析了5×5棒束组件弯棒对于CHF的影响,结果表明,当弯曲使相邻加热棒间隙减至50%以上,才观察到CHF降低,未加热的弯曲棒(冷壁)甚至与加热的杆接触也没有不利影响。Ren等[8]对5×5带定位格架的弯棒进行了数值模拟,弯曲位置发生在两格架间。结果表明,当弯曲程度增加到燃料棒间隙变化的50%以上时会减小CHF,在0～50%时没有影响,预测的DNB位置发生在弯曲段。以上实验或数值模拟采用的燃料棒弯曲位置均为内部棒且弯曲仅发生在两格架位置间,弯曲并未对边角通道产生直接影响,并未关注棒束全高度整体弯曲、偏移对临界的影响。在早期的OECD开展的实验及近期上海交通大学开展的CHF实验中发现,燃料棒可能发生大跨度的弯曲或由于加工工艺出现小幅的偏移。

为研究棒弯曲和偏移对CHF的影响,本文基于欧拉两流体模型和壁面沸腾模型[9-12],研究发生弯曲或偏置的5×5燃料组件的CHF,并比较3种变形对于CHF数值和位置的影响,为燃料组件的安全分析提供依据。

1 数值模型

过冷沸腾和沸腾临界的液相与气相温度并不相同,即主流过冷,近壁面可能过热[13]。另外,由于密度差和相间作用,两相之间出现滑移速度,因此应考虑相间的不平衡[14]。两流体模型可考虑相界面质量、动量和能量的传输,基于欧拉两流体模型分别建立气、液两相守恒方程[15-17]。同时,引入壁面热流密度分配模型描述换热过程壁面向气、液两相的热量传递过程。此外,为使方程封闭,引入辅助模型,如气泡脱离直径、脱离频率等[18]。

1.1 欧拉两流体模型

质量守恒方程:

(1)

动量守恒方程:

(2)

能量守恒方程:

(3)

1.2 壁面沸腾模型

CFD方法可基于两种方法模拟两相流,即界面捕获法和相位平均法。第1种的计算网格尺寸需远小于气泡尺寸,同时为了保持数值稳定性的时间步长,该方法需要巨大的计算资源。第2种基于统计界面面积浓度和一系列辅助模型来模拟两相流[19]。壁面沸腾模型是受热壁面沸腾应用最广泛的经验模型,为了模拟临界,本研究采用拓展后的壁面沸腾模型,即CHF模型。该模型将壁面热流密度分为4部分:

qwall=(qc+qq+qe)f(αl)+(1-f(αl))qg

(4)

式中:qwall为壁面热流密度;qc、qq、qe、qg分别为液体单相对流传热、淬灭对流传热、蒸发热流传热及气体单相对流传热部分的热流密度,表达式[19]如下:

qc=hl(Tw-Tl)(1-Ab)

(5)

(6)

qe=VdρgNwhfvf

(7)

qg=hg(Tw-Tg)

(8)

式中:hl为单相换热系数;Tw为壁面温度;T为液相温度;1-Ab为加热壁面液体所覆盖的面积;Vd为脱离壁面的气泡体积;Nw为成核密度;ρg为气相密度;hfv为蒸发潜热;f为气泡脱离频率;kl为液相导热系数;cpl为比定压热容。

f(αl)为加热壁面液相所占的份额,表达式为Tenter[20]关系式:

f(αv)=1-f(αl)=

(9)

其中,αv,1和αv,2为过渡点,分别为0.9和0.95。

1.3 界面传输模型

欧拉两流体模型守恒方程中,源项考虑了界面传输项,因此需使用界面传输模型来进行求解。对于相间传质,考虑了两种机制。第1个是加热表面的液体蒸发及气泡进入过冷主流的冷凝或主流中过热液的气化。对于相间动量传递,考虑了曳力和非曳力。曳力由流体的分子黏度和两相间的滑移速度所引起,作用于滑移速度的负方向。非曳力作用在与流动方向正交的方向上,包括速度梯度引起的升力、将气泡推向主流的壁面润滑力及展平径向空泡份额的湍流耗散力。考虑从界面到液相和从界面到气相的热阻,通过双热阻模型模拟相间传热。气泡对湍流强度的影响通过气泡尾流引起的附加湍流源项来模拟,相关模型列于表1,对于相间动量模型,Ishii模型适用于沸腾流动不同尺寸、形状的气泡的曳力计算。Moraga升力系数模型对液滴和气泡所受升力适用性良好。Antal对于近壁面气相速度分布差异大的区域能够适用。Burns基于相间曳力的计算,与实验值的依赖性相对Bertodano弱。针对相间传质,文中计算均为DNB型沸腾临界,主要的物理现象仅存在受热壁面蒸发及进入主流冷凝,根据壁面沸腾模型中的蒸发热流来推导相间传质。对于相间能量传递,通常假定两相界面温度为相应压力下饱和温度,分别计算界面向液相间的传热,和界面向气相间的传热。表1中:hfv为汽化潜热;Tsat为饱和温度;Tl为液相温度;Tv为气相温度;h为比焓;Ai为单位体积内气液两相间相界面面积;db为气泡直径;σ为表面张力;g为重力;CD为曳力系数;CL为升力系数;CWL为壁面润滑力系数;Ftd为湍流耗散力系数;yw为网格中心到最近壁面的距离;σvl为常数0.9;CTD为常数1;μl为液相动力黏性系数;αv为气相份额;αl为液相份额;Nub为努塞尔数;Pr为普朗特数;qvt为单位体积内气相和液相相界面间的换热量;δt为时间尺度默认值为0.05(FLUENT默认);νv为运动湍流黏度;k为湍动能;ε为湍流动能耗散率;Uv为气相平均速度;Ul为液相平均速度。

表1 界面传输相关模型Table 1 Correlation model for interfacial transfer

表2 各工况变形类型Table 2 Deformation types under various conditions

1.4 数值方法

通过将非线性偏微分方程离散化为线性代数方程来求解。本文求解稳态控制方程,因此不需要时间离散化,通过伪瞬态方法求解以提高收敛性。动量方程中的压力项采用二阶方法离散化,梯度项通过基于网格的最小二乘法得到。离散得到的线性代数方程用Guass-Seidel方法求解,质量和动量方程同时用COUPLED算法求解。

2 几何、网格划分及边界条件

考虑到本文研究重点为偏置和弯曲对CHF的影响,为排除定位格架影响,本文研究对象为5×5光棒的燃料组件,包壳外径为9.5 mm,棒束栅距为12.6 mm,横截面结构如图1所示,同时给出子通道划分以及燃料棒编号。组件加热高度设置为2 m,不加热的出口段长约130 mm,燃料组件三维几何结构如图2所示。

图1 几何截面及子通道划分示意图Fig.1 Schematic diagram of geometric section and sub-channel division

图2 三维几何结构示意图Fig.2 Schematic diagram of three-dimensional geometric structure

在网格划分方面,首先利用ICEM对无偏置和弯曲的基准算例进行六面体网格划分,经网格独立性分析,最终网格总数409万。在基准算例的基础上通过动网格技术,使网格发生偏置与弯曲,共设计3种类型,即X型(横向偏置)、XY型(对角线偏置)、C型(轴向弯曲)。3种类型网格均为整体移动(即所有包壳都发生移动),以下在3种类型中各挑选1个算例网格进行展示,如图3所示,其中X、XY偏置示意图选择右上角4个子通道进行放大,图中展示结果均为变形1 mm工况的网格示意图。图4示出X、XY、C型变形具体示意图,3种变形均为1 mm,X、XY、C型变形方向为x轴右侧。图中对于中心坐标轴,X型偏置所有燃料棒向右侧位移,XY型偏置向对角位移。图中C型弯曲,几何进口位置对应弧度0,加热段末端对应弧度3.14,而弯曲最大值(1 mm)在加热高度中心位置(1 m),整体弯曲效果示于图3。另外动网格移动后的网格正交质量均在0.55以上。

图3 基准算例、X型、XY型及C型网格示意图Fig.3 Schematic diagrams of base case, X-type, XY-type and C-type meshes

图4 3种变形示意图Fig.4 Schematic diagrams of three kinds of deformation

本文在基准工况的基础上进行棒偏置和弯曲影响规律研究,基准工况中设计压力15.5 MPa、质量流量3 000 kg/(m2·s)、出口含气率0、进口温度534.66 K及CHF参考值2 424 kW/m2。根据压力、含气率及质量流密度查询CHF查询表(look-up table)[28],再对查询值通过水力直径和棒束修正因子进行修正,从而得到最终的CHF查询表数值,即参考值。在数值模拟过程中,从0.5倍参考值逐步提升功率,功率在50%～70%区间以10%提升功率,超过70%时以1%提升功率直至临界。

在基准工况基础上,通过动网格技术实现燃料棒的整体偏置或弯曲,获得新结构下的计算网格。设计以下工况,各工况类型列于表3。将需要计算的工况分两类,A类9个工况轴向功率均匀分布,又分3组,即X型、XY型、C型,每组3个变形量分别为0.5、0.8、1.0 mm;B类2个工况轴向功率呈余弦分布[29],均为C型。其中,1.0 mm的变形属于较大的变形,0.5 mm和0.8 mm分别为小变形和中等变形工况。目前公开的文献中,并没有对变形量取值的直接依据,实验中通常以两燃料棒弯曲至接触的局部现象出现,为了简化问题,文中采用沿1个方向变形的方式来处理。CHF数值对于边角通道间隙宽度很敏感,这样的结果使得变形1 mm对于CHF的数值及位置影响非常明显,同时0.5 mm的变形临界位置依然在热通道,综上选取变形量0.5、0.8、1 mm。

表3 边界条件Table 3 Boundary condition

由于燃料包壳外表面周向热流密度分布较为均匀,在建立燃料组件模型时可不考虑燃料棒内的导热过程,直接在包壳外表面上给定热流密度。轴向加热功率包括均匀和余弦分布两类,详细的边界条件列于表2、3。

3 模型验证

将上述工况预测的CHF与查询表数值进行对比,结果列于表4。查表法是将FLUENT计算临界高度位置发生临界的局部子通道水温、空泡份额、流量及压力取出,用水温及空泡份额转换为含气率,根据局部子通道压力、流量、含气率及当量直径,获得最终的CHF查询表数值。表中列出5个工况的CFD预测值、查询表数值及误差。其余工况,临界位置均由中心迁移到边缘,而在边、角通道位置处CHF查询表并不适用。这是因为CHF查询表基于8 mm的圆管制成,对于非8 mm的圆管或者棒束结构(如内部子通道)可通过水力直径修正因子进行修正,往往预测精度高。然而在棒束组件的边角通道位置,受热边界及湿周并不均匀,加上燃料棒偏置或弯曲对边角通道进行挤压,导致该位置受热边界及湿周分布更不匀衡,因此CHF查询表不再适用[30]。而对于临界位置发生在热通道的工况,CFD预测值与CHF查询表数值误差均在5%以内,预测精度良好。

表4 查询表数值与预测值对比Table 4 Comparison of sub-channel look-up table value and prediction value

4 结果分析

临界热流密度主要受流量、压力、含气率参数影响。图5示出无变形的基准工况临界高度位置液相温度和空泡份额云图分布,图中可发现中心子通道温度最高,临界位置也发生在该处(黑色虚线框住包壳位置),边角通道温度次之。空泡份额峰值位于热棒最小间隙处,最大值约为0.94,触发沸腾临界。图6示出基准算例空泡份额及温度随流向分布的趋势,图中截面最大空泡份额(即包壳)及截面平均空泡份额均随流向而增加,在出口达到最大,趋势与主流温度曲线基本一致;主流温度出口位置仍低于饱和温度,为DNB型沸腾临界。包壳温度先升高,中间段达到624 K左右保持不变,在靠近出口位置发生飞升。

图5 基准算例临界高度处的液相温度及空泡份额云图分布Fig.5 Cloud maps of water temperature at critical height and vapor fractions distribution of benchmark

燃料棒的整体偏置或弯曲使不同子通道冷却剂流量重新分配,从而影响传热和沸腾临界。图7～10及表5比较了不同类型工况边、角及中心子通道总质量流量分布、沸腾临界位置和临界高度处温度云图。图中不同颜色代表不同变形程度,对应颜色标记了沸腾临界发生位置。由表5可知,对于同种工况,15号子通道质量流量最大,边通道(12号)次之,角通道(6号)最小。对于同类变形的不同工况,随变形程度的增大,中心子通道(15号)质量流量通道增加,边角子通道(12、6号)质量流量降低。这是由于燃料棒整体偏置或弯曲,对右侧边角通道造成挤压使得通道流通面积减小,冷却剂流量降低。对于X型和XY型相同偏置距离的工况,XY型12号边通道质量流量降幅更小,15号子通道上升幅度更小(除位移0.5 mm非常接近),X型偏置工况相比XY型有更多边通道的流量向中心通道流动;而在6号角通道质量流量降幅更大,在位移达到1 mm时,XY型6号角通道质量流量降幅比X型高6.5%,因为燃料棒沿对角线位移,若将其正交分解,燃料棒沿x和y方向同时挤压6号角通道,导致该处流通面积极窄,冷却剂流量降幅更大。即X型对于边通道的影响比XY型大,在角通道的影响弱于XY型偏置,对于整体流量的影响X型更强。对于C型弯曲工况,15号子通道除轴向功率均匀弯曲1 mm的工况比X型同变形程度稍小,其余C型工况中心子通道上升幅度比X和XY型都更大,而功率曲线为余弦分布时中心通道上升幅度最大。整体而言,结合边棒临界情况C型弯曲对于燃料组件子通道流场的影响最大。

图7 X型临界高度处径向临界位置及液相温度分布Fig.7 Radial critical position and water temperature distribution at X-type critical height

表5 不同工况子通道总质量流量分布Table 5 Total mass flux distribution of sub-channels under different conditions

这种流量分布的结果,导致温度云图的分布趋势发生改变,高温区集中分布在被挤压的边角通道,低温区集中远离燃料棒位移或弯曲方向的子通道。图7中可发现X偏置0.8、1 mm时边棒临界,偏置0.5 mm和基准工况(黑色)临界发生在热通道,温度云图高温区向位移方向移动,集中在右侧边棒位置。图8中最高温度区域随位移发生移动,集中在6号角通道,31号角通道冷却剂温度最低。XY偏置0.5、0.8 mm的临界位置发生在热通道,偏置1 mm边棒临界,仅在6号角通道位置触发临界,流量分析分布已解释该现象。

图8 XY型临界高度处径向临界位置及液相温度分布Fig.8 Radial critical position and water temperature distribution at XY-type critical height

对于轴向均匀功率及非均匀功率的C型弯曲工况结果如图9、10所示,图中可发现温度云图分布趋势与X型工况类似。对于轴向功率均匀的弯曲工况,0.8、1 mm时边棒临界,弯曲0.5 mm时临界位置发生在热通道;而轴向功率呈余弦形,弯曲0.5、0.8 mm时边棒临界。即弯曲0.5 mm不同功率曲线,径向临界位置发生变化。

图9 C型临界高度处径向临界位置及液相温度分布Fig.9 Radial critical position and water temperature distribution at C-type critical height

图10 C型(cosine)临界高度处径向临界位置及液相温度分布Fig.10 Radial critical position and water temperature distribution at C-type (cosine) critical height

CHF限制反应堆功率输出水平,它关系到反应堆堆芯的安全性和经济性。将上述12个工况预测的CHF和临界高度汇总成数据图,如图11所示。从图11可看出,无论是均匀加热或非均匀加热功率,对于整体趋势,X、XY、C型变形程度越大,CHF越小。轴向非均匀功率与弯曲共同作用下的工况CHF最低。对于均匀加热功率,相同变形程度,XY型对CHF影响最小;以基准工况为参考,位移1 mm 的X型工况预测的CHF比弯曲1 mm偏差更大;对于0.5、0.8 mm变形,X型预测的CHF比C型更接近基准工况CHF。综上,在某个变形量阈值前,C型对CHF影响更大,在阈值后X型对CHF影响更大。

图11 不同工况CFD预测临界值对比Fig.11 Comparison of CFD predicted critical values under different conditions

图12示出不同工况临界高度对比。对于轴向均匀功率,X、XY型和弯曲0.5 mm工况不影响临界高度,仅0.8、1 mm的C型弯曲改变临界高度,随弯曲程度的增大,临界位置向上游小幅迁移。而轴向非均匀功率的两个C型工况中,临界高度大幅前移,位于中部偏下游位置。

图12 不同工况轴向临界位置对比Fig.12 Comparison of axial critical positions under different conditions

5 结论

本文基于欧拉两流体模型对发生偏置或弯曲的5×5棒束燃料组件CHF进行了数值模拟研究,获得了横向偏置(X型)、对角偏置(XY型)和弯曲(C型)对CHF的影响,通过分析可得到如下结论。

1) 对于轴向均匀功率的工况,X型偏置及C型变形0.8、1 mm时,出现边棒临界;XY型偏置1 mm出现边棒临界。对轴向余弦功率的工况,弯曲0.5、0.8 mm均发生边棒临界;弯曲0.5 mm余弦与均匀功率的工况对比,带余弦功率工况触发边棒临界,即余弦功率影响边棒临界。当变形程度增加时,边棒临界的可能性增加。

2) 对轴向功率均匀的工况,X、XY型和0.5 mm的C型弯曲对于CHF位置没有影响。0.8、1 mm的C型弯曲工况,随弯曲程度的增大,临界位置向上游移动;而对于轴向余弦功率的C型弯曲,趋势也一致。

3) 3种变形程度越大,预测的CHF越小;XY型变形对CHF敏感性最弱,轴向功率余弦的C型工况对CHF影响最大,相比同等弯曲0.5 mm均匀功率工况,临界数值降低9%。对于均匀功率的X、C型变形工况,0.5、0.8 mm C型弯曲的临界数值偏差更大(与基准工况CHF对比),而变形1 mm的X型临界数值偏差更大,即两类变形对临界数值的影响在0.8～1 mm区间存在交点,在交点前C型对临界数值影响更大,交点后X型对临界数值影响更大。

4) 燃料棒的整体偏置或弯曲使子通道流量重新分配,在边角通道受燃料棒偏置弯曲影响大,XY型对6号角通道影响最大,而对燃料组件整体流量分配影响最小;C型弯曲对燃料组件整体流量分配影响最强。当边角通道缝隙小到冷却剂流量不能有效带走燃料棒表面气泡,造成气泡局部拥塞,将触发边棒临界。