压水堆燃料棒芯块事故后传热的简化分析

齐宇博，张　伟，于　江，南金秋，赵剑刚

压水堆燃料棒芯块事故后传热的简化分析

齐宇博，张伟，于江，南金秋，赵剑刚

（中广核研研究院有限公司，广东 深圳 518124）

压水堆核电厂在严重事故下，堆芯换热条件恶劣，此时包含衰变热的堆芯非稳态传热分析过程较为复杂。本文以大亚湾核电厂M310机组一回路热段双端断裂为假想工况，用简化分析方法研究事故后衰变热的传递情况，获得了含时间变量的温度场简化计算公式。采用ANSYS有限元分析软件，用非简化分析的方法计算分析，两种方法对比验证。研究结果显示：在一定条件下，简化计算的误差相对较小，可采取简化分析方法对事故后燃料棒芯块的传热进行计算。本文的简化方法不仅局限于压水堆，其他堆型的事故分析亦可借鉴。

瞬态导热；集总参数法；有限元法；严重事故

非稳态传热过程中，控制方程为几何空间的函数，是时间尺度一阶偏导，对较复杂工程问题，几乎不能获得解析解。然而，当求解域内热阻较小，可近似认为温度场在瞬态均匀一致，此时导热微分方程简洁且易于计算，该方法称之为集总参数法[1]。然而传统的集总参数法，多数情况下回避研究对象内部释热，原因是内热源随时间变化，难以用确定的函数关系式表征，即使某些情况下得到表达式，控制方程依然无法解析。许多工程问题又存在内热源，这使得集总参数法在实际应用过程中受到一定程度的限制。

压水堆核电厂一回路发生大破口失水事故（LOCA）时，在紧急停堆且未启动再循环阶段，燃料棒芯块的衰变热近似指数衰减。本文借助origin后处理软件，用非线性指数函数拟合该衰变关系，获得较理想的结果。将拟合关系式引入非稳态导热微分方程中，求解并获得分析事故后燃料芯块瞬态导热的简化公式。同时，用ANSYS有限元软件进行对比验证，在毕渥数（v）不超过0.113 75，最大误差为13.1%。

因此，该简化方法是可应用于事故后燃料元件传热的初步分析。

1　研究的对象

压水堆核电厂在LOCA事故后，两台低压安注泵和两台安全壳喷淋泵投运，从安注箱和换料水箱取水。大约持续1 200 s，上述水源枯竭，安全注入系统（RIS）和安全壳喷淋系统（EAS）切换到从地坑取水，进入长期再循环阶段。

在正常运行工况下，燃料芯块中心局部热点温度可达1 670 ℃[2]，为探索集总参数法在严重事故下应用可行性，假设UO2热物性参数不随富集度、辐照、温度等变化，忽略芯块肿胀和变形。

1.1　物性参数

UO2芯块物性参数[2]如表1所示。

表1　物性参数

1.2　热边界条件

LOCA事故后紧急停堆，假设燃料芯块初始温度分布均匀，取1 670 ℃；燃料棒环境温度为90 ℃[3]。

便于后续计算和分析，燃料棒芯块与环境间热量交换等效成表面复合换热系数，分别取20 W/（m2·℃）、200 W/（m2·℃）和500 W/（m2·℃）三种工况进行参数化分析。

事故后，燃料棒芯块衰变热如图1所示。

图1　LOCA事故下芯块体积释热率随时间曲线

2　公式的推导

非稳态、有内热源导热微分方程：

忽略温度场随几何空间的变化，导热微分方程简化为：

——换热面积；

——体积；

——边界等效换热系数；

导热微分方程可进一步写成：

为常数项，由初始条件进行决定。

LOCA事故过程中，前1 200 s堆芯的衰变热接近指数衰减，如图2所示。

图2　origin 拟合

非线性拟合结果：

将方程（7）带入（6），求解：

常数项C：

根据方程（8）和（9），代入物性参数，就可得到LOCA事故后1 200 s内，燃料芯块瞬态导热的近似解。

3　ANSYS验证

实际中，依托计算机技术进行数值分析，瞬态导热分析常采用有限差分法（FDM）和有限元法（FEM）。首先将求解域在空间上离散，对单元体进行高阶差值或线性插值，利用计算机获得较高计算精度和自适应性。

瞬态热力学分析一般方程：

式中：[]——传导矩阵，包括导热系数、对流系数及辐射系数和形状系数；

[]——比热矩阵；

使用ANSYS有限元软件时，基于燃料芯块对称性，取1/4圆柱体建模，网格划分如图3所示。

图3　网格划分

3.1　表面复合换热系数h=20 W/（m2·℃）

通常用毕渥数（v）来表征物体内部热阻与外部热阻间关系，其定义：

其中：——物体长度量纲；

——物体表面传热系数；

芯块的表面复合换热系数一般不会达到该量级，但当发生超设计基准事故，例如堆芯堵塞，到达堆芯的流量不足，燃料组件上部区域处核态沸腾，换热系数可低至该量级[3]。

简化公式（8）和FEM计算结果如表2所示，FEM稳态温度分布云图如图4所示。简化公式（8）和FEM计算结果体现在曲线中，如图5所示。

表2　数据对比

注：① 指用简化公式（8）计算结果。

② 指用FEM方法计算的结果。

③ 指用FEM方法计算芯块表面温度和中心温度差值，与简化公式（8）计算的平均温度之比。

图4　1 200 s FEM计算结果

结果显示：考虑芯块导热热阻，1 200 s后中心最高温度为876.84 ℃，表面温度为872.48 ℃，温度差为3.56 ℃。相对于整体温度，温差非常小。

图5　结果比较

从变化曲线看，两种方法计算结果非常吻合，相对误差小于0.5%。

3.2　表面复合换热系数h=200 W/（m2·℃）

在该工况下，v=0.045 5。

简化公式（8）和FEM计算结果如表3所示，FEM稳态温度分布云图如图5所示。简化公式（8）和FEM计算结果体现在曲线中，如图6所示。

表3　数据对比

注：① 指用简化公式（8）计算结果。

② 指用FEM方法计算的结果。

③ 指用FEM方法计算芯块表面温度和中心温度差值，与简化公式（8）计算的平均温度之比。

图6　1 200 s FEM计算结果

结果显示：考虑芯块的导热热阻，在该工况下，1 200 s后中心最高温度为166.55 ℃，表面温度为163.22 ℃，温度差为3.33 ℃。相对于整体温度，该温差并不大。此时，简化分析方法与有限元计算结果的符合性较高。

从图7中可见，两种方法计算结果吻合，相对误差小于5%，从工程角度来讲是可接受的。

图7　结果比较

3.3　表面复合换热系数h=500 W/（m2·℃）

在该工况下，v=0.113 75。

简化公式（8）和FEM计算结果如表4所示，FEM稳态温度分布云图如图8所示。简化公式（8）和FEM计算结果体现在曲线中，如图9所示。

表4　数据对比

注：① 指用简化公式（8）计算结果。

② 指用FEM方法计算的结果。

③ 指用FEM方法计算芯块表面温度和中心温度差值，与简化公式（8）计算的平均温度之比。

图8　1 200 s FEM计算结果

结果显示：考虑芯块导热热阻，在该热边界条件下，1 200 s后中心最高温度为122.61 ℃，表面温度为119.28 ℃，温度差为3.33 ℃。

图9　结果比较

起始段相对误差较大，最大相对误差达到13.1%，60 s后两种方法计算结果相对吻合，误差小于5%，此时芯块已接近稳态。

4　结论

实际物体的几何形状复杂，通常控制方程难以获得解析解，而集总参数法有效地避开几何尺度影响，使分析过程大大简化，当毕渥数较小时，获得相对精确计算结果。有国外文献[3]报道，当堆芯面临堵塞，堆芯上部为核态沸腾区，对流换热系数处于102量级。此时采用集总参数法进行分析和计算是可行的，有限元计算也证明了本文所推导出公式有较高的计算精度。

本文将堆芯衰变功率拟合为指数形式，忽略芯块导热热阻，求解一阶非线性微分方程获得研究芯块导热简化公式，并采用ANSYS有限元数值解法对比验证，当毕渥数较小时，能够获得满意的计算结果。本文的研究思路并不仅限于压水堆事故分析，其他堆型亦可借鉴。

［1］ 杨世铭、陶文铨. 传热学［M］．北京：高等教育出版社，1998：66-69.

［2］ 广东核电培训中心. 900 MW压水堆核电厂系统与设备［M］．北京：原子能出版社，2004：50-53.

［3］ Westinghouse Electric Company. Evaluation of Long-Term Cooling Considering Particulate、Fibrous and Chemical Debris in the Recirculating Fluid，WCAP-16793-NP［R］．U.S，2011.

［4］ 黄厚诚，王秋良. 传导热问题的有限元分析［M］．北京：科学出版社，2009：1-3.

［5］ Westinghouse. Evaluation of Downstream Sump Debris Effects in Support of GSI-191，WCAP-16406-P-A［R］. U.S，2008.

Simplified Method on the Heat Conduction of PWR Fuel Rods During the Accident

QI Yubo，ZHANG Wei，YU Jiang，NAN Jinqiu，ZHAO Jiangang

（China Nuclear Power Technology Research Institute Co.，Ltd，Shenzhen of Guangdong Prov.518124，China）

During the serious accident of PWR nuclear power plants，such as a loss of coolant accident，the analysis of core decay heat transfer will be complicated。This paper is based on the hypothetical condition of a double-ended loop break of the Daya Bay M310 and uses simplified method to research the decay heat transfer，by which getting the simplified calculating formula of the temperature filed via the time variable.Meanwhile，the results of the calculation are compared with the results of ANSYS software to testify.The results have showed that the errors of calculations can be accepted when taking a simplified analytical method on a certain conditions.This simplified analytical method not only can be used for PWR，but is also available for the other types of reactors.

Transient conduction；Lumped parameter method；Finite element method；Severe accident

TLLT331

AA

0258-0918（2022）01-0053-06

2020-09-24

国家科技支撑计划（2011BAA06B01）资助

齐宇博（1986—），男，陕西定边人，工程师，学士，现主要从事核电、氢能方面研究