脉冲反应堆脉冲运行瞬态的子通道分析

陈立新，朱 磊，马腾跃，江新标，陈 达

（1.南京航空航天大学 核科学与工程系，江苏 南京 210016；2.西北核技术研究所，陕西 西安 710024）

脉冲反应堆是一种以铀氢锆为燃料的池式研究堆，由于铀氢锆材料具有较大的瞬发负温度系数，因此该类型反应堆既可稳态运行，也可脉冲运行。对于脉冲运行工况，反应堆的功率变化快，脉冲功率峰的宽度达ms量级，功率变化的范围大，脉冲前的稳态功率与脉冲峰值功率相差6～7个量级，这种快速的功率变化对堆芯的传热安全提出了更高要求。准确计算脉冲运行状态下堆芯的热工参数是评价脉冲反应堆瞬态安全的重要内容。

子通道模型是目前核反应堆堆芯热工水力分析中较精细的计算模型。由于子通道模型能准确给出反应堆各组件或各燃料通道内的流量、温度和压力等热工参数的分布，故基于子通道模型的分析程序成为反应堆热工水力分析的重要工具。世界上主要核能国家分别推出了各系列的子通道分析程序，如THINC 系列［1］、COBRA 系列［2］等。这些较成熟的商用子通道程序大多针对动力反应堆设计，其适用范围也多集中在高温、高压、高质量流量的强迫循环冷却方式，对脉冲堆这类低温、常压、低质量流量的自然循环冷却反应堆缺乏适用性。由于以往的脉冲堆瞬态热工安全分析多采用较保守的单通道方法，无法准确给出堆芯内部详细的热工安全参数的分布，因此，本文基于子通道分析方法，开发脉冲堆瞬态特性分析的子通道程序，并以西安脉冲反应堆为研究对象，分析堆芯的瞬态热工安全特性。

1 子通道物理模型

1.1 子通道模型的数学方程

一般，子通道分析常采用均匀平衡模型，即假设气、液两相的相对速度为零，且具有相同的压力与温度。基于以上假设，可将冷却剂看作气、液两相的混合物而进行统一处理，其连续性方程及能量和动量方程［3］如下。

连续性方程为：

能量守恒方程为：

动量守恒方程为：式中：ρ为流体密度；τ 为时间；u 为速度矢量；ui、uk（i、k=1，2，3）分别为不同的速度分量；E 为流体具有的能量；T 为应力张量；f 为体积力；Q 为流体发热量。

以上方程组的求解非常困难，为简化方程，一般的子通道分析方法均需进行必要假设，认为通道中冷却剂的轴向流速度远大于横流速度，以至于横流在通过通道间的间隙后就失去了其原来具有的方向而随轴流运动。这样就可将动量方程分解为轴向动量守恒方程和横向动量守恒方程。以上子通道方程的推导参见文献［3］。

1.2 子通道方程的数值解法

对于子通道方程，可采用有限差分方法对其进行求解。将堆芯划分为N 个子通道，并将各通道在轴向上分为L 层，将每层看作1个控制体，如图1所示。在控制体内部，流体具有单一的温度、压力、密度、比容等参数。将子通道守恒方程应用于每个控制体，则可得到N×L组守恒方程。

对于离散化方程守恒方程，在每个时间步长Δτ内按轴向分层逐层对各通道进行迭代求解，得到本层的所有参数后，进入下一层的求解，其中，能量方程需同燃料元件导热方程联合求解。

图1 子通道网格划分Fig.1 Mesh generation of sub-channel

2 辅助模型

2.1 功率计算模型

瞬态分析时，由于功率是时间的函数，且堆芯温度的改变会影响功率的变化，因此需耦合热工分析模块与功率计算模块进行计算。为避免复杂的物理热工耦合，本文采用考虑温度反馈的点堆动态方程计算堆芯的功率变化，同时也可根据物理计算程序的功率模拟结果将各燃料元件或堆芯整体的功率-时间变化曲线作为外部参数输入。

2.2 燃料元件导热模型

对于燃料元件的导热，采用粗棒非稳态一维导热模型，该模型可用于模拟燃料元件稳态和瞬态导热问题，其通用控制方程为：

式中：r为与热量传递方向平行的坐标；F（r）为与导热面积有关的因子；Q 为内热源；λ为导热系数；T 为温度。对于式（4），可将其在时间和空间上进行离散，采用差分法进行求解。

2.3 热传输模型

对于燃料元件包壳与冷却剂之间的对流换热系数的计算，采用对流换热的分区处理方法，即传热曲线被划分为4个分区［4］。单相对流传热区采用Sieder-tate关系式，过冷沸腾区采用McAdams关系式，过渡沸腾传热采用对数坐标的沸腾曲线，取临界热流密度和最小热流密度线性内插法估算传热系数，膜态沸腾传热采用简化的处理方法。图2为本文计算的脉冲堆传热曲线。

图2 脉冲堆燃料包壳与冷却剂的传热曲线Fig.2 Heat transfer curve of pulsed reactor between clad and coolant

2.4 水和水蒸气物性的计算模型

选用水和水蒸气性质国际联合会（IAPWS）发布的IAPWS-IF97公式［5］计算水和水蒸气的物性。IAPWS-IF97公式将水和蒸气按性质的不同分成5个区：1区为常规水区，2区为过热蒸气区，3区为临界区，4区为湿蒸气区，5区为低压高温区。本文主要用到1、2 区的数据。IAPWS-IF97 适 用 范 围 为273.15 K ≤T ≤2 273.15K，压力p≤100MPa。每个分区的物性参数计算公式参见文献［5］。

3 计算程序的编制及验证

3.1 程序编制

在以上对子通道模型进行离散化的基础上，结合燃料元件导热模型、包壳冷却剂对流换热模型、瞬态功率计算模型、水及水蒸气物性计算模型等编制脉冲堆瞬态安全分析的子通道程序PRC-STAC。程序采用FORTRAN 语言编写，内部计算模块采用子函数形式，主程序通过输入输出控制对各子函数的调用。程序主要包括输入输出模块、功率计算与处理模块、燃料元件温度场分析模块、换热系数计算模块、水力计算模块、辅助参数计算模块等功能模块。

3.2 程序验证

本文通过对OSTR 堆［6］典型的瞬态计算结果开展比对计算，来验证PRC-STAC程序的可靠性。OSTR 反应堆属于美国GA 公司研制的TRIGA MARK Ⅱ型反应堆，与脉冲反应堆结构类似。该堆额定运行功率为1.1 MW，反应堆工作于常压环境，堆芯依靠水的自然循环冷却。本文选取了RELAP 5-3D 软件计算的OSTR 反应 堆 采 用LEU 燃 料（235U 富 集 度 为19.75%，与脉冲堆相同）的堆芯瞬态计算结果。为使二者更具可比性，选择文献［6］中8通道计算模型给出的计算结果进行比对。

图3 示出了两个程序计算的OSTR 堆堆芯燃料最高温度随时间的变化曲线计算比对结果。计算的燃料温度变化趋势一致，燃料到达最高温度的时间略有差别。图4为燃料温度达最大值时温度在燃料棒径向上的分布计算结果。由比对结果可看出，本文的计算结果与RELAP5-3D的符合较好。

图3 燃料最高温度随时间的变化Fig.3 Fuel maximum temperature vs time

4 脉冲堆瞬态特性分析

为准确模拟反应堆在脉冲工况下堆芯热工参数的时间-空间变化，采用三维时空动力学程序XAPR-HENKO 对脉冲堆在引入3.5＄脉冲时的堆芯三维功率分布进行计算［7］，脉冲工况的动力学参数列于表1。

图4 燃料元件温度沿径向的变化Fig.4 Fuel element temperature variation along radial direction

表1 脉冲工况的动力学参数Table 1 Kinetics parameters in pulse operation

分析时选取沿堆芯径向布置的1组燃料元件进行计算，其栅元编号为D2、E2、F2、G2、H2、I2，其在脉冲堆芯内的位置如图5 所示。图6为XAPR-HENKO 程序计算的以上几根燃料元件的功率随时间的变化曲线。

4.1 单棒计算

D2测温燃料元件位于脉冲堆芯脉冲运行时的热点位置，因此首先对D2元件及其周围的冷却剂进行分析。图7 为发射3.5＄脉冲后不同时刻D2 燃料元件内部燃料温度在不同时刻的径向分布。从图7可看出，在发射脉冲的最初阶段，燃料温度的分布呈边缘高、中间低的分布特征，这主要是由于粗棒元件自屏效应的影响。随时间的变化，燃料温度逐渐趋于一致。

图5 选取的燃料元件在脉冲堆芯的位置Fig.5 Position of fuel elements in core

图6 燃料元件功率随时间的变化Fig.6 Fuel element power vs time

图7 脉冲时D2元件燃料径向温度变化Fig.7 Fuel radial temperature variation of D2element in pulse operation

由于在脉冲发射时，燃料边缘的温度最高，图8示出了燃料边缘和包壳温度随时间的变化曲线。

脉冲进行时，燃料边缘的温度随反应堆功率的上升迅速升高，而包壳温度的上升则滞后于燃料温度的变化，显然，这是由于热量在向包壳传递过程中存在时间延迟。燃料内部的能量释放来自核裂变反应，可认为瞬时释放全部能量，热量在短时间内来不及向包壳传递；因此包壳的温度变化会滞后于燃料边缘，同时由于包壳与冷却剂直接接触，其温度的最大值也较燃料温度低很多，但最高温度也接近500 ℃。这一瞬间高温持续时间极短，不会造成燃料元件不锈钢包壳的烧毁，但包壳温度的快速变化与将炙热的金属快速浸入冷水的情况相似，会使不锈钢包壳表面产生类似烧灼淬火的痕迹，在对西安脉冲堆燃料元件进行表面观测时，也确实发现了包壳表面呈蓝黑色的现象。

图8 D2元件燃料边缘和包壳的温度变化Fig.8 Temperature variation of fuel and clad for D2element

图9示出了D2元件周围冷却剂通道轴向主流温度变化情况，在开始时刻，冷却剂轴向最高温度并不在通道的出口处，而是位于出口下方某一位置，随时间的增加，出口处冷却剂温度逐渐升高。总体分析，冷却剂主流温度在脉冲后上升并不显著，这主要是由于脉冲发射持续的时间很短，实际释放到冷却剂中的积分能量并不高。

图9 D2燃料元件周围冷却剂主流温度随时间的变化Fig.9 Main flow temperature variation of coolant around D2element

4.2 堆芯分析

选取D2、E2、F2、G2、H2、I2燃料元件分析反应堆堆芯。图10示出了6根燃料元件最高温度的对比。D2、E2、F2元件的最大瞬时温度要显著高于其他元件，这是因为这3根元件布置于脉冲棒周围，发射脉冲时，其功率峰因子最大。图11示出了冷却剂温度沿堆芯径向的分布，此处给出的是临近D2、E2、F2、G2、H2、I2燃料元件的冷却剂温度分布图。在堆芯内侧，由于有中央水腔，冷却剂温度较低；堆芯外侧，因为布置了大量的燃料元件，脉冲发射后，通道内的水温变化则较大。

图10 燃料元件最高温度分布Fig.10 Maximum temperature distribution of fuel element

图11 冷却剂温度沿堆芯径向的分布Fig.11 Coolant temperature distribution along radius of core

5 结论

通过本文建立的脉冲堆子通道分析方法，对西安脉冲堆的脉冲运行瞬态进行了理论分析，结果表明：在3.5＄脉冲工况下，堆芯热棒（D2）温度和包壳最高温度较大，但并未突破燃料堆芯温度的限值，在脉冲开始阶段燃料最高温度会形成1个小的温度峰值，但其持续时间很短，由于整个脉冲过程释放的能量有限，不会对燃料元件造成不可接受的安全影响。

［1］ CHU P T，HOCHREITER L E，CHELEMER H，et al.THINC-Ⅳ：An improved program for thermal hydraulic analysis of rod bundle cores［M］.US：Westinghouse Electric Corporation，1973.

［2］ STEWART C W，WHEELER C L，CENA R J，et al.COBRA Ⅳ：The model and the method，BNWL 2214［R］.US：Pacific Northwest Laboratory，1997.

［3］ MERROUNA O，ALMERSB A，BARDOUNIA T，et al.Analytical benchmarks for verification of thermal-hydraulic codes based on sub-channel approach［J］.Nuclear Engineering and Design，2009，239：735-748.

［4］ CHEN Lixin，TANG Xiaobin，JIANG Xinbiao，et al.Theoretical study on boiling heat transfer in the Xi’an Pulsed Reactor［J］.Science China：Technological Sciences，2013，56（1）：137-142.

［5］ IAPWS.Release on the IAPWS international formulation 1997for the thermodynamic properties of water and steam［R］.Erlangen，Germany：IAPWS，1997.

［6］ WADE R M.Thermal hydraulic analysis of the Oregon State TRIGA Reactor using RELAP5-3D［D］.US：Oregon State University，2008.

［7］ 赵柱民，陈立新，缪正强，等.脉冲堆脉冲工况下堆芯功率时空分布模拟［C］∥中国核学会2009年学术年会.北京：中国核学会，2009.