环形燃料栅元有效温度计算方法研究

陈 钊肖英杰赵鹏程彭梁兴

1（中广核研究院有限公司 深圳 518000）

2（南华大学核科学技术学院 衡阳 421001）

环形燃料作为一种新型的核燃料概念设计，其冷却剂同时流经内通道表面和外通道表面以进行双面冷却，并使功率密度增加原来值的30%以上，提供了更高的安全裕度［1］。为了进一步了解环形燃料的特性，国内外通过理论分析和实验对其中子特性、热工水力性能和力学性能开展了大量的研究工作［2−5］。根据燃料棒共振有效温度的物理意义，燃料共振有效温度的计算过程分为两个过程：一是通过中子学与燃料热力学行为的迭代，得到燃料棒中的真实的功率分布和温度分布；二是利用等效原理计算燃料共振有效温度。

针对压水堆环形燃料子通道的热工水力性能分析，韩国原子能研究院开发了程序MATRATHAF［6］、美国麻省理工学院开发了程序TAFIX［1］、并与国际知名多功能热工水力分析程序VIPRE-01进行了对比验证［7］；在国内，环形燃料热工水力分析的相关程序开发工作也正在逐步展开，如中国原子能科学研究院的刁均辉等［8］开发了程序SAAF，西安交通大学的吴攀等［9］也开发了子通道分析程序NACAF。为准确分析环形燃料栅元的热工水力性能，本研究基于两相均质流模型，开发了适用于压水堆的环形燃料单通道热工水力分析程序THCAFS（Thermal-Hydraulic Code of Annular Fuel with Single channel）。

此外，针对有效温度和径向功率分布等相关方面的研究，尹强等［10］基于燃料棒共振有效温度的机理性计算模型开发了FRET程序；Pirouzmand和Roosta［11］研究了VVER-1000燃料中燃耗和原子密度的径向分布；Lemes［12］分析了高燃耗水平下燃料芯块的径向结构；Yuan等［13］发现了UO2燃料中燃耗和反应速率径向分布的简单公式；Chen等［14］根据不同慢化剂燃料比和不同内外冷却剂比，研究了环形燃料的功率和燃料温度的径向分布。本研究以西屋公司四环路压水堆所设计的13×13排列的环形燃料元件［1］为基准，基于蒙特卡罗程序SERPENT模拟燃料棒中功率分布和燃耗过程，通过自主开发的环形燃料两相程序THCAFS模拟燃料棒中的热力学行为，研究不同燃耗下的径向功率分布和温度场分布，采用总反应率等效代替有效共振积分来建立环形燃料栅元的有效温度计算模型。

1 栅元温度场计算程序开发

对于环形燃料，裂变能量只有一部分热量被传递到内通道冷却剂，内通道与外界无质量、动量上的交换，另一部分能量则传递到外通道冷却剂中，而外通道可与外界进行质量、动量间的交换。在开发环形燃料单通道热工水力分析程序THCAFS时，由于主要研究对象为环形燃料单通道模型，故外通道不考虑质量、动量的交换，但需考虑到内、外通道的冷却剂流量分配和热量分配。

1.1 数学物理模型

1.1.1几何结构

针对高功率密度的压水堆堆芯的热工水力分析，THCAFS程序以典型的西屋公司四环路压水堆13×13排列的环形燃料元件作参考对象，功率为参考功率的150%，栅元的详细几何结构如图1所示。

图1 环形燃料结构示意图（r1=4.316 5 mm，r2=4.888 0 mm，r3=4.950 0 mm，r4=7.050 0 mm，r5=7.112 0 mm，r6=7.683 5 mm，S=16.51 mm）Fig.1 Schematic diagram of annular fuel rod structure(r1=4.316 5 mm,r2=4.888 0 mm,r3=4.950 0 mm,r4=7.050 0 mm,r5=7.112 0 mm,r6=7.683 5 mm,S=16.51 mm)

1.1.2两相计算模型

裂变产生的热量经燃料芯块向内、外通道传递，依次通过气隙、包壳，最后传至冷却剂，而冷却剂流经加热通道时将可能发生流动沸腾。由于伴随相变所产生的气泡明显改变了冷却剂的传热和流动特性，大大增加了两相计算模型的复杂程度，目前通常采用假设两相流体的基本参数仅沿通道轴向变化的方法对模型进行简化［15］。如图2所示，根据空泡份额随轴向高度的变化关系，通常可分为4个传热区域：单相区、高过冷沸腾区、低过冷沸腾区及饱和沸腾区。对于单相区的对流传热，由于气泡没有脱离壁面进入主流，通常采用Dittus关系式，而内外通道两相参数的计算模型则包括沸腾换热模型、空泡份额模型、两相摩擦倍率模型等。沸腾换热模型采用Thom关系式，空泡份额模型采用Zuber-Findly关系式，两相摩擦倍率模型采用均匀流模型。

图2 加热通道传热分区示意图Fig.2 Heat transfer partition diagram of heated channel

1.1.3热量分配模型

由于环形燃料通过内、外通道流动的冷却剂同时冷却，在THCAFS程序中需要将热量同时分配给内通道和外通道。为求解环形燃料芯块区域外的温度场，假设在芯块内部的温度最高处为虚拟绝热面，将芯块划分为内、外两层，内层产生的能量用于内通道冷却剂加热，外层产生的能量用于外通道冷却剂加热。在计算环形燃料温度分布前，需假设绝热面半径rm，如果芯块裂变功率均匀分布，则内、外环芯块体积释热率之比为：

而对于环形燃料芯块区域，采用平均体积释热率qv，W·m−3，并通过导热微分方程求解芯块内部的温度场分布，得到燃料芯块内部径向温度解析式［16］：

式（2）边界条件为：r=r3，T=Ti和r=r4，T=To；此外，系数C1和C2的表达式如下：

式中：To和Ti分别为外通道和内通道的燃料芯块表面 温 度，K；r为 半 径，m；ku为 芯 块 热 导 率，W·(m·K)−1。

1.1.4流量分配模型

由于环形燃料存在内、外两个通道，内、外通道流量分配的准则是内、外通道压降基本一致，本研究根据流道形状和流动特点，考虑了摩擦压降∆Pf、重力压降∆Pg、加速压降∆Pa及形阻压降∆Pc等，并在有限次数的迭代中实现上述准则。因此，沿轴向的冷却剂通道的总压降∆Ptot为：

式（5）可表示为如下的具体形式：

式中：f是摩擦系数；ϕ2是两相摩擦倍增因子；k是形阻系数；V是冷却剂流速，m·s−1；ρ是冷却剂密度，

kg·m−3。

当将通道内的流动视为全液相流动时，两相摩擦压降可通过单相摩擦压降和两相摩擦倍增因子的乘积表示，采用的ϕ2的关系式如下：

式中：xf为流动含气率；νfg和μfg分别为饱和液体与饱和蒸汽的比体积之差，m3·kg−1，以及黏性系数之差，Pa·s；νfs和μgs分别是饱和液体的比体积，m3·kg−1，以及饱和蒸汽的黏性系数，Pa·s。

根据能量守恒定律，内、外通道冷却剂吸收的热量为：

式中：Cp是冷却剂定压比热容，J·(kg·K)−1；Mtot是总质量流量，kg·s−1；下标i和o分别表示内通道和外通道，in和out表示入口和出口。

1.2 程序开发流程

THCAFS程序开发流程如图3所示，首先输入环形燃料几何参数与运行条件，并假设绝热面rm的初始值，以获取内、外通道的热量分配比；然后，在假设热量分配的基础上开展内、外通道冷却剂的流量分配，判断通道内所存在的传热区域，同时计算不同通道的换热系数Hx、冷却剂温度Tf、DNBR（Department from Nucleate Boiling Ratio）和压降∆P等参数，再选择内、外通道的压降为收敛准则k1，若满足条件则进行下一步计算，否则重新分配流量；最后，在已知流量分配的基础上开展燃料栅元的温度场计算，寻找划分内、外环的绝热面位置rc，若满足收敛条件k2，则输出计算结果，否则需要修正rm，重新分配热量。

图3 THCAFS程序流程图Fig.3 Flow chart of THCAFS code

1.3 程序验证

本研究所选用的栅元验证对象详见图1，其平均线功率密度为74.3 kW·m−1，轴向功率呈余弦函数（轴向峰因子1.55），冷却剂总质量流量为0.783 57 kg·s−1，入口压力为15.5 MPa，入口温度为567.87 K。程序中流量分配的迭代计算嵌套在热量分配迭代计算内，并有两个重要收敛条件，若在有限的迭代内同时满足流量分配和热量分配条件下，则输出内、外通道的计算结果。

内通道的对流换热系数如图4（a）所示。在环形燃料中下段区域为单相换热，4个不同程序计算的内、外通道换热系数变化平缓且吻合较好。在ONB（Onset of Nucleate Boiling）点之后，流体换热模式变为核态沸腾换热，且随着高度上升通道换热系数迅速上升，当换热系数达到最高值，由于内外通道热流密度降低，换热系数也随之迅速下降。而且内通道接近出口位置由于燃料呈余弦功率分布，其热流密度较低，流体换热重新进入单相换热区域。外通道的对流换热系数如图4（b）所示，在环形燃料上半段，外通道换热系数计算原理与内通道相似，主要区别在于外通道后半段存在饱和沸腾区域。

图4（c）和（d）分别计算了内、外通道的DNBR，呈先下降后上升的趋势，总体吻合度较高。对于内通道，无定位格架，故THCAFS程序采用W-3公式计算；对于外通道，则采用考虑定位格架修正的W-3L公式［17］计算。

图4 环形燃料热工程序结果对比(a)内通道换热系数，(b)外通道换热系数，(c)内通道DNBR，(d)外通道DNBRFig.4 Comparison of annular fuel thermal hydraulic code(a)Heat transfer coefficient of inner channel,(b)Heat transfer coefficient of outer channel,(c)DNBR of inner channel,(d)DNBR of outer channel

表1列出了各环形燃料程序之间的对比验证。结果表明，与其他环形燃料子通道程序相比，THCAFS程序开展的流量分配与压降计算结果在可接受范围之内。

表1 质量流量与压降计算验证Table 1 Verification of mass flow and pressure drop

程序THCAFS与NACAF在热点处的径向温度场如图5所示。自主开发的环形燃料栅元温度场计算程序THCAFS与VIPRE-01、TAFIX和NACAF等程序的计算值吻合度较高，而且与NACAF所对比的栅元温度场结果再次证明了THCAFS程序的可行性。

图5 热点处燃料栅元温度场结果验证Fig.5 Verification of the temperature field of fuel cell at the hot spot

2 栅元有效温度计算方法研究

本文从中子学模拟、燃耗过程模拟、热力学模拟等方面来进行模拟分析，研究不同燃耗下的径向功率分布和温度场分布，采用反应率等效代替有效共振积分来建立环形燃料的共振有效温度计算模型。栅元有效温度求解流程如图6所示，利用蒙特卡罗程序SERPENT构建计算对象的精细化模型（芯块区域划分12环），对于不同的材料区域进行初始温度赋值。然后，设置初始燃耗步长，分别开展中子学计算与燃耗过程计算，得到径向功率分布与主要核素密度变化。根据不同燃耗和不同核素密度下的径向功率分布，基于THCAFS程序进行热力学模拟，并将真实的栅元温度场计算结果重新赋值给精细化模型的各材料区域，以计算其真实总反应率。

图6 环形燃料栅元有效温度求解流程图Fig.6 Flow chart for solving effective temperature of annular fuel cell

2.1 中子学模拟

由于空间自屏效应，热中子通量在燃料芯块内部分布不均，易造成环形燃料径向功率的分布不均。本文基于三维蒙特卡罗中子输运的代码SERPENT和连续能量截面库ENDF/B-VⅡ.0进行中子学模拟，通过求解中子输运方程的方式模拟中子在燃料中的分布，进而得出环形燃料中的功率分布。为描述传统燃料棒中的径向功率分布，作为给定空间位置处燃耗的函数，Chen等［18］根据研究提出了多项式函数：

式中：s为燃耗深度，MWd·kgHM−1；系数取决于相对半径x。

在环形燃料中，相对半径x可用式（10）表示：

而系数a(x)、b(x)、c(x)可表示为：

其中：i(x)=a(x)，b(x)，c(x)；常数A、B、C、D、E由其模拟计算的径向分布可确定。系数A、C和E通过模拟结果的最小二乘拟合获得，B和D的优化值由式（12）给出：

2.2 燃耗过程模拟

对于环形燃料来说，其核素密度随着燃耗的加深而发生变化，芯块内部不同位置的功率密度并不相同，导致径向不同位置的温度也会不同。选取的计算对象详见图1，以富集度为4.95%的UO2作芯块材料，包壳材料为锆-4合金，冷却剂为水，芯块与包壳的间隙填充氦气，且平均线功率密度为74.3 kW·m−1，比功率为103.3 MW·tU−1。计算时采用的燃耗范围为0~60 MWd·kgHM−1，初始温度设置为：燃料900 K、包壳600 K、冷却剂600 K。根据不同的燃耗深度条件下，获得燃料中不同位置处不同核素的核子密度的变化过程。如图7所示，在燃料芯块区域沿径向均匀划分12环，考虑主要核素的密度变化，进行精细化建模和计算。

图7 环形燃料栅元精细化建模Fig.7 Refined modeling of annular fuel cell

2.3 热力学模拟

热力学模拟是在考虑上述的功率分布、燃耗深度以及核素密度变化等因素的基础上，基于THCAFS程序求解环形燃料的真实温度场分布。在计算温度场时，首先进行轴向控制体划分，而考虑燃料芯块内部功率分布不均等影响，需进一步对燃料进行径向的控制体划分，如图8所示。

图8 环形燃料径向节点划分Fig.8 Radial node division of annular fuel

1）第j个冷却剂控制体：

2）第j个包壳控制体外表面：

3）第j个包壳控制体内表面：

4）第j个芯块控制体外表面：

式中：j为轴向控制体标号；ql(z)和qs(z)分别是在轴向位置z处的线功率密度，W·m−1，以及面积释热率，W·m−2；Cp为比热容，J·(kg·K)−1；G为质量流密度，kg·(m2·s)−1；De为冷却剂通道的当量直径，m；hx（z）为轴向控制体的对流换热系数，W·(m2·K)−1；kclad为包壳热导率；kgas为间隙气体热导率，W·(m·K)−1；dcs为包壳外表面外径；dci为包壳内表面直径；du为芯块直径，m；Tf(z)为轴向冷却剂控制体温度；Tcs(z)为包壳控制体外表面温度；Tci(z)为包壳控制体内表面温度；Tu(z)为芯块控制体表面温度，K。

3 结果与分析

图9为不同燃耗下的环形燃料径向相对功率分布。除中间区域布置12个探测器外，两侧边缘处分别增加了1个探测器。由SERPENT所模拟计算的径向功率呈U型趋势，两侧高中间低，而且随着燃耗的加深不均匀程度更大，两侧越高，中间越低。这是由空间自屏效应效应引起的，热中子通量在燃料芯块的中间较为平坦，而在芯块边缘迅速增加。在燃耗过程中，燃料元件的空间自屏效应造成功率分布不均匀的同时，由于边缘处的235U和238U更易吸收中子产生239Pu和241Pu及其他裂变产物导致燃料中核素成分的分布也不一致，从而对燃料芯块中径向功率的分布也造成一定的影响。

图9 不同燃耗下的径向功率分布Fig.9 Radial power distribution under different burnups

图10为燃料芯块中235U、238U、239Pu和241Pu的平均原子密度随燃耗的变化曲线。在燃耗过程中235U的平均原子密度大体呈下降趋势，这种变化随燃耗加深而逐渐变缓；238U主要通过俘获中子形成239Pu逐渐消耗，下降趋势基本一致；239Pu和241Pu随着燃耗的增加逐渐增大，并且在平均燃耗为60 MWd·kgHM−1时，燃料芯块中239Pu的原子密度高于235U。在燃耗初期产生的热能主要来自235U吸收中子发生的裂变反应，燃耗后期的热能则来自235U、239Pu和241Pu三种核素的裂变。

图10 不同燃耗下的核素密度变化Fig.10 Changes of nuclide density under different burnups

本文以40 MWd·kgHM−1的燃耗深度为例，通过燃耗过程模拟得到计算对象中的主要核素密度变化，同时由中子学模拟计算其径向功率分布，并对模拟结果进行二阶多项式函数的拟合。图11为多项式中系数a（x）、b（x）和c（x）关于相对位置x的拟合函数，而且拟合结果较好。由以上系数可确定与环形燃料空间位置和燃耗相关的径向相对功率解析表达式，而不同燃耗下功率的拟合函数值与模拟计算值的比值见图12，最大偏差不超过2%，表明该函数可用于计算径向功率分布。在考虑上述径向功率分布、燃耗深度及核素密度变化等因素基础上，通过THCAFS程序求解的燃料栅元温度场分布如图13所示。

图11 多项式系数a(x)、b(x)和c(x)的拟合Fig.11 Fitting of polynomial coefficients a(x),b(x)and c(x)

图12 不同燃耗下的拟合函数与模拟结果的比值Fig.12 Ratio of fitting function value to simulation value under different burnup

图13 热点处燃料栅元温度场（40 MWd·kgHM−1）Fig.13 Temperature field of fuel cell at the hot spot(40 MWd·kgHM−1)

考虑到环形燃料栅元中真实温度分布对核素截面的影响，基于SERPENT程序主要针对可分辨共振能区进行了在线截面生成处理。为求解燃料共振的有效温度，首先假设燃料棒中燃料温度为某一平均温度，使采用平均温度后燃料棒的有效共振积分与采用真实的温度分布时的有效共振积分在数值上是相等的。

本文采用总反应率等效代替有效共振积分来计算环形燃料栅元的有效温度。具体求解过程为：在得到燃料棒中的真实燃料温度后，假设燃料棒中燃料温度为某一平均温度，通过使用假设的平均温度进行中子学计算，如果得到的反应率与采用真实燃料温度的反应率不一致，则调整该平均温度，直至反应率等效，最终求得的平均温度即为燃料共振有效温度。如图14所示，在精细化建模（芯块均匀划分12环）的基础上，将图13中热点处的实际温度分布赋值给燃料芯块后得到真实总反应率5.366 81×1017n·cm−3·s−1，再不断调整原始模型（芯块未分环）的材料温度计算反应率。若计算反应率与真实反应率的相对误差不超过0.01%，输出在条件范围内的平均等效温度1 008 K。

图14 热点处燃料栅元有效温度计算（40 MWd·kgHM−1）Fig.14 Calculation of effective temperature for fuel cell at the hot spot(40 MWd·kgHM−1)

4 结语

本研究以西屋公司设计的环形燃料栅元为计算对象，开发了可预测内、外通道压降、换热系数、偏离泡核沸腾比及栅元温度场等参数的环形燃料单通道热工水力分析程序THCAFS。利用蒙特卡罗程序SERPENT开展中子学模拟和燃耗过程模拟，分析其径向功率分布不均和主要核素的密度变化，再通过自主开发程序THCAFS模拟环形燃料中的热力学行为，研究不同燃耗下的温度场分布，采用反应率等效代替有效共振积分来建立环形燃料栅元的有效温度计算模型。得出主要结论如下：

1）通过THCAFS与其他环形燃料热工水力程序 如VIPRE-01、TAFIX和NACAF进 行Code-to-Code对比，其预测结果基本一致，可初步应用于环形燃料设计以及热工水力分析；

2）从中子学模拟、燃耗过程模拟、热力学模拟方面，较真实地模拟环形燃料在燃耗过程中发生的物理现象，并采用反应率等效代替有效共振积分建立了其有效温度计算的方法模型，可为相关环形燃料共振有效温度的机理性研究提供重要的参考性价值。

作者贡献声明陈钊：起草文章，分析/解释数据；肖英杰：采集数据；赵鹏程：获取研究经费，对文章的知识性内容作批评性审阅，行政、技术或材料支持；彭梁兴：分析/解释数据。