池式钠冷快堆事故余热排出系统一回路仿真研究

姜 博，张智刚，于 洋，陈广亮，张志俭

（哈尔滨工程大学 核安全与仿真技术国防重点学科实验室，黑龙江 哈尔滨 150001）

池式钠冷快堆事故余热排出系统一回路仿真研究

姜 博，张智刚＊，于 洋，陈广亮，张志俭

（哈尔滨工程大学 核安全与仿真技术国防重点学科实验室，黑龙江 哈尔滨 150001）

池式钠冷快堆事故余热排出系统采用了非能动工作原理，依靠液态钠及空气的自然对流排出堆芯余热。为研究事故工况下余热排出系统一回路的换热能力，基于FORTRAN语言，建立堆芯单通道及盒间流模型，采用全隐二阶迎风差分格式及改进的欧拉法离散求解，对事故余热排出系统一回路系统进行数值模拟，并对全厂断电事故进行仿真计算验证。结果表明：该程序能较好地反映事故余热排出系统瞬态变化过程，并可达到超实时仿真。

余热排出系统；自然循环；盒间流模型；数值模拟

池式钠冷快堆（PSFR）事故余热排出系统是反应堆专设安全设施系统之一，主要由位于一回路的堆芯、独立热交换器及位于二回路的空气热交换器等组成。其中，独立热交换器（DHX）与一回路主要换热设备中间热交换器（IHX）均为逆流垂直式壳管型热交换器，因此二者可采用相同的仿真方法。丘林等［1］采用节点模型建立IHX的热工模型，若所划分的节点数大于50，则仿真结果较好；当节点数小于50时，其两侧流体流量失衡时会发生温度震荡。而王武军等［2］采用全隐二阶迎风差分格式，用45个节点解决流量失衡时的温度振荡问题，并通过对EBR-Ⅱ失热阱实验瞬态和失流实验的实际计算来验证程序效果。模拟结果表明该方法具有稳定性，可用较少的节点达到较高的精度。

由于事故停堆后较长时间内，盒间钠流动对堆芯冷却起到重要作用，但采用盒间流模型对事故余热排出系统进行模拟的研究却鲜见报道。本文采用盒间流模型，对池式钠冷快堆事故余热排出系统进行仿真研究。

1 系统描述

事故余热排出系统一回路主要包括两个封闭式流道，在全厂断电事故发生初期由主热传输系统依靠泵的惰转排出堆芯热量，而当风门打开后，由流经独立热交换器及燃料组件盒间隙的钠带走余热，如图1所示［3］。

图1 事故余热排出系统示意图Fig.1 Schematic of decay heat removal system

图2为事故余热排出系统一回路钠流动示意图［4］。堆容器上部热钠池内区中的钠通过溢流窗进入热钠池外区，然后由上向下流入独立热交换器壳程，与中间回路自下而上流入管程的钠换热。一回路钠流出独立热交换器后分为两路，一路沿堆芯内外屏蔽柱腔室向下流，经过围桶上的侧孔流入燃料组件之间，最后流出堆芯，进入热钠池。另一路沿出口方向流入热钠池，流量较小，本文不予考虑。

图2 事故余热排出系统一回路钠流动示意图Fig.2 Schematic of sodium flow for decay heat removal system in primary loop

堆芯反射层、乏燃料及热屏蔽组件位于堆芯边缘区域，其释热在主热传输系统有效工作时，依靠主热传输系统冷却剂带出，当主热传输系统冷却剂流量为零时，依靠流经独立热交换器及屏蔽的自然循环冷却剂带出。

2 数学仿真模型

本文所作的主要假设有：1）将池式钠冷快堆事故余热排出系统简化为由堆芯、热钠池、冷钠池、独立热交换器、堆芯内外屏蔽柱腔室等主要的5部分组成；2）堆芯内外屏蔽柱腔室等效为绝热流道；3）堆芯采用单通道模型进行仿真计算；4）独立热交换器管束简化为单管，并忽略一次钠、管壁及中间回路钠的轴向传热。程序中物性如焓、比热、密度等，使用文献［5］所提供的数学模型求解。

2.1 热工水力基本方程

利用动量守恒方程，求解回路冷却剂流量：

视液态钠流动为不可压流动，将式（1）沿一封闭路线积分，假定积分路径可分为m段，方程化简为：

式中：W为质量流量，kg／s；t为时间，s；L为沿流动方向长度，m；A为流道横截面积，m2；ρ为密度，kg／m3；f为分布阻力系数；De为当量直径，m；K为阻力系数；z为每段长度，m；pen为进口压力，Pa；pst为出口压力，Pa。

基于能量方程（式（3）），对堆芯、独立热交换器、冷钠池、热钠池等进行数学建模。

式中：cp为比定压热容，J／（kg·K）；T为温度，K；Q为换热功率，W；v为速度，m／s。

1）堆芯热工模型

余热排出系统运行时，堆芯衰变热依次经燃料芯块、包壳、盒内冷却剂、燃料盒、盒间冷却剂等导出，由盒内冷却剂和盒间冷却剂将热量带至热钠池。

图3为堆芯节点划分示意图。堆芯为圆柱形燃料棒，将燃料棒划分为燃料芯块中心处（节点1）、燃料芯块1／2半径处（节点2）、燃料包壳内表面（节点3）、燃料包壳外表面（节点4）等4个节点，燃料芯块内的热量经导热及对流换热传给燃料包壳、盒内冷却剂及盒间冷却剂。由于燃料棒长度远大于燃料棒直径，因此计算时仅考虑径向导热，用一维瞬态热传导方程来求解燃料温度。

（1）节点1温度。根据节点2的温度在节点1处的泰勒展开式、中心对称条件及洛必达法则可得：

（2）节点2温度。由节点1及节点3的温度泰勒展开公式得到节点2计算公式，化简后得到：

（3）节点3温度。芯块与包壳接触处的导热可用一等效换热系数来表示：

由式（6）、热流密度连续条件公式及节点2在节点3处的泰勒展开公式，化简后得到：

（7）乏燃料组件释热计算式。计算时，假设堆芯反射层、乏燃料、热屏蔽组件功率占堆芯总功率的1%，分配总流量的10%对其进行冷却，保证冷却性能。

式中：λ为导热系数，W／（m·K）；h为换热系数，W／（m2·K）；qv为体积释热率，W／m3；Δr为节点间距，m；LK1为堆芯流量分配系数；LK2为堆芯功率分配系数；HFO为乏燃料组件出口焓，J／kg；HFI为乏燃料组件入口焓，J／kg；i为节点；下标f1、f2分别为盒内、盒外。

2）冷、热钠池传热模型

热钠池位于堆芯组件及堆芯支承环腔上板的隔热层以上，并由生物屏蔽支撑筒将其分为内、外两区，电离室在生物屏蔽支撑筒内与堆芯平行，热钠池内区顶部为堆芯测量柱，一次钠循环泵部分以及两个独立热交换器位于热钠池外区，热钠池基本结构如图4所示［4］。冷钠池位于主容器下部、上支撑板以下，并被分为两个密闭的隔间，每个隔间内有2台中间热交换器及1台钠循环泵。

图4 热钠池剖面图Fig.4 Cross section of hot sodium pool

从堆芯及乏燃料组件流出的钠在热钠池内区混合，通过生物屏蔽支撑筒上的溢流窗流入热钠池外区，并与从独立热交换器流出的冷钠混合后，流入堆芯内外屏蔽柱腔室、中间热交换器等。

图5 独立热交换器节点图Fig.5 Node diagram of DHX

以独立热交换器二次侧入口为基准，向上划分节点，节点如图5所示。冷、热钠池传热模型与主热传输系统传热模型相同，模型如下［4］。

热钠池内区：

其中：δ为壁厚，m；下标hot、cold分别代表冷、热钠池。

3）独立热交换器传热模型

热钠池外区的钠流入独立热交换器，与中间回路钠进行换热。将独立热交换器划分为一次侧、二次侧及管壁等3个区，分别用p、s、m代表，其能量方程［2］分别为：

其中，S为换热湿周，m。

2.2 辅助模型

1）堆芯盒内流与盒间流对流换热方程

堆芯内盒内冷却剂与盒间冷却剂对流换热系数由式（19）确定：

h1与h2分别为盒内钠与燃料盒以及盒间钠与燃料盒之间的对流换热系数，由式（20）确定［6］：

式中：δm为燃料盒厚度；λm为燃料盒的传热系数；λ1、λ2分别为盒内、盒间冷却剂的导热系数；D1、D2分别为盒内、盒间当量直径。

2）独立热交换器热侧与冷侧对流换热方程

独立热交换器一、二次侧流体等效对流换热系数由式（21）确定：

h3与h4分别为一次钠与传热管以及二次钠与传热管之间的对流换热系数，由式（22）确定［6］：

式中：δg为传热管厚度；λg为传热管的热传导系数；λ3、λ4分别为管内、管外冷却剂的导热系数；D3、D4分别为管内、管外的当量直径。

燃料盒材质为铁素体-马氏体钢，独立热交换器材质为304不锈钢，物性参数取相应材料参数。

3）堆芯衰变功率模型

停堆后，堆芯衰变功率［6］为：

Q（t）＝Qn（t）＋Qd（t）＋Qc（t）（23）其中：Qn（t）＝Q0（0.1e－0.15t），为停堆后缓发中子引起的裂变功率，W；Qd（t）＝Q00.068［t－0.2－（t＋t0）－0.2］，为停堆后裂变产物的衰变功率，W；Qc（t）＝Q0（1.63×10－3e－4.92×10－4t－1.60× 10－3e－3.4×10－6t），为停堆后俘获产物衰变功率，W。

4）泵惰转模型

一、二回路系统两条支路各有1台离心式钠循环泵。根据文献［3］，一回路钠循环泵惰转时间为91s，泵从990r／min惰转到550r／min的时间为10～11s，设泵的半时间为11s；二回路钠循环泵惰转时间为131s，设泵的半时间为13s。由带有飞轮的离心泵的惰转模型［7］及泵实际惰转规律得出一、二回路泵及流量双段式（双曲线段及直线段）泵惰转模型：

其中：np为泵转速，r／min；k为系数；下标0为初始值。

3 离散方程及求解方法

针对式（2）～（25）中微分方程，所需求解量主要为冷却剂温度及流量，分别采用二阶迎风格式及改进的欧拉公式进行离散求解。

3.1 二阶迎风格式

在对流-扩散差分格式中，迎风格式是绝对稳定的。但一阶迎风格式当Peclet数较大时，假扩散较严重，需采用细密的网格进行计算［8］。虽然一阶迎风格式不会出现解的振荡，但其截差较低而限制了解的准确度，所以采用一阶导数具有二阶截差的偏差分格式，即二阶迎风差分格式。

利用二阶迎风差分格式，将独立热交换器划分为30个节点，采用内节点法进行计算，将方程（15）～（17）进行离散，离散方程如下。

节点1、2：

其中，k为时刻。对式（4）～（34）中微分方程采用Gear算法求解。

3.2 改进的欧拉公式

梯形公式具有较高的精度，但算法复杂，计算量大。为了简化算法，采用改进的欧拉公式［9］对动量方程（式（2））进行求解，其离散格式为：

4 全厂断电事故仿真验证

4.1 主要假设

基于全厂断电事故，对事故余热排出系统一回路进行仿真模拟。主要假设［3］为：

1）反应堆初始功率取102.5%额定功率；

2）堆芯入口温度为363℃（正常为360℃）；

3）当失去厂外电源后，柴油发电机组系统未能启动，一、二回路主泵直接惰转到0；

4）空气热交换器的应急排热系统600s后打开风门。

4.2 仿真结果

根据“电网电源丧失”信号，全厂断电事故发展序列列于表1［3］。

每台独立热交换器额定功率为0.525MW，两台独立热交换器总额定功率占反应堆装置功率的1.6%。备用工况（反应堆负荷为100%，风门开度为10%）下，仿真程序采用保守估计（即空气温度为50℃，流量为0.11kg／s）时独立热交换器二次侧的冷却能力（入口温度为485℃，流量为1.37kg／s）计算，求得独立热交换器一次侧入口温度为516.62℃，出口温度为491.10℃，与设计值（入口温度为516℃，出口温度为490℃）相比，相对误差在0.25%以内。当反应堆功率为两台独立热交换器功率之和，即1.05MW时，所得独立热交换器流量稳定在5.80kg／s左右，符合设计标准值。

表1 全厂断电事故发展序列Table 1 Time series of station blackout accident

事故工况下，余热排出系统空气冷却器风门全部打开时间为63s。仿真计算时，以中间回路钠进口温度及流量为边界条件，计算独立热交换器二次侧出口温度、一次侧入出口温度、流量及堆芯出口温度等。其中，独立热交换器二次侧入口温度采用保守设计假设，在事故发生后，流量由1.37kg／s升至2.93kg／s，二次侧温度由485℃降至345℃，并在其后的时间内以每1 000s降低1℃的冷却能力运行，由于事故余热排出系统空冷器空气侧入口温度设计值为50℃，钠侧出口温度可降至250℃，因此在很长的一段时间内，独立热交换器二次侧入口温度是保守的。

图6 堆芯功率变化Fig.6 Power of core

图6为堆芯功率变化。当全厂断电事故发生时，反应堆紧急停堆，泵惰转，堆芯功率急剧下降。主冷却系统冷却剂流量及二次侧流量随泵的转速下降而不断降低，最终趋于零。

图7示出独立热交换器一次侧入、出口温度变化。事故发生600s后，事故余热排出系统风门在63s内开度由10%升至100%，独立热交换器换热功率迅速上升，独立热交换器出口温度急剧下降。由于前600s内独立热交换器二次侧入口处温度为定值，因而随着一次侧入口温度降低且逐渐接近于二次侧入口温度，独立热交换器功率下降，独立热交换器一次侧出入口温差趋于0。事故发生600s后，在63s内，空冷器风门打开，独立热交换器二次侧入口温度急剧下降，其一次侧出口温度随之剧烈降低（图7、8）。随着功率平稳下降及独立热交换器二次侧入口温度趋于稳定，独立热交换器出口温度下降逐渐平缓。仿真计算中，双环路对称，模拟结果趋于一致。

图8 独立热交换器二次侧入、出口温度变化Fig.8 Temperatures of DHX inlet and outlet for the second loop

图9为堆芯盒间冷却剂入、出口温度变化。由于主热传输系统无盒间流模型，因此堆芯盒间流在风门打开后引入。在风门打开后，由于堆芯衰变热在短时间内依然大于独立热交换器能带出的热量，故盒间流堆芯出口温度持续升高，并在896s时达到最高温度527.67℃。随反应堆功率的下降，堆芯出口温度开始下降，并在前4 000s下降速度较快，而后趋于稳定降低。堆芯盒间冷却剂出口温度与盒内冷却剂流量及温度、堆芯余热、盒间冷却剂流量、盒间冷却剂入口温度等相关，所受影响较为复杂。

图9 堆芯盒间冷却剂入、出口温度变化Fig.9 Coolant temperatures of core inlet and outlet between fuel assemblies

图10为独立热交换器流量变化。由于主泵转动停止后主冷却系统流量急剧下降，此时衰变热功率依然较高，仅依靠主冷却系统流量不能带出，因而反应堆堆芯温度不断上升，与独立热交换器内的钠形成密度差而建立自然循环，流量上升。当风门打开后，由于中间回路钠温下降，独立热交换器换热能力增强，其出口钠温降低，进入堆芯的冷却剂温度降低，导致堆芯平均温度降低，因而短时间内流量下降。而随着更多的热量需由盒间钠带出，自然循环流量再次上升，在686s时达到最大值6.23kg／s。此后，随热钠池温度的降低，二回路温度稳定，以及堆芯功率下降缓慢，流经独立热交换器的流量下降趋于平缓，7 500s时降至4.8kg／s，10 000s时降至4.6kg／s。

图10 独立热交换器流量变化Fig.10 Mass flow of DHX

5 结论

本文采用二阶迎风格式和改进的欧拉法，建立堆芯平均通道及盒间流模型，对池式钠冷快堆事故余热排出系统一回路进行了仿真计算，结论如下。

1）采用盒间流模型建模，符合物理过程，有效地模拟出事故余热排出系统工作方式。仿真结果表明，全厂断电事故过程中，堆芯出口温度达到峰值527℃（0.2h）后下降，峰值温度与钠沸腾温度相比有很大的裕度；3h后流量下降，趋于4.5kg／s；独立热交换器出口温度随堆芯功率的下降而下降，即事故余热排出系统可完成事故工况下堆芯余热排出，并能保证堆芯安全。

2）采用二阶迎风格式对独立热交换器进行离散，可用较少的节点达到仿真效果而不产生温度震荡，提高了计算速度，实现超实时仿真；该程序已经移植到SimExec实时仿真平台中，与快堆其他系统相连接，运行情况良好，可对快堆运行进行实时仿真研究。

［1］ 邱林，桑维良.池式钠冷快堆堆芯系统动态特性和稳定性的仿真研究［J］.核科学与工程，1998，18（4）：311-318.QIU Lin，SANG Weiliang.Research on stability and dynamic character for primary system of FFR［J］.Chinese Journal of Nuclear Science and Engineering，1998，18（4）：311-318（in Chinese）.

［2］ 王武军，单建强，朱继洲，等.全隐二阶迎风差分格式在快堆模拟中的应用［J］.核科学与工程，1999，19（4）：309-314.

WANG Wujun，SHAN Jianqiang，ZHU Jizhou，et al.Mehtod of fully implicit second-order upwind difference in fast reactor simulation［J］.Chinese Journal of Nuclear Science and Engineering，1999，19（4）：309-314（in Chinese）.

［3］ 徐銤.快堆主热传输系统及辅助系统［M］.北京：中国原子能出版社，2011：56，76-77，171-186.

［4］ 于洋.池式钠冷快堆一回路系统实时仿真模型与程序研究［D］.哈尔滨：哈尔滨工程大学，2013.

［5］ 饶彧先，崔满满，郭赟.中国实验快堆一回路热工水力稳态计算程序开发［J］.原子能科学技术，2012，46（9）：1 067-1 073.

RAO Yuxian，CUI Manman，GUO Yun.Development of thermal-hydraulic steady-state analysis program for primary loop of China Experimental Fast Reactor［J］.Atomic Energy Science and Technology，2012，46（9）：1 067-1 073（in Chinese）.

［6］ 徐銤.快堆热工流体力学［M］.北京：中国原子能出版社，2011：72-75，93.

［7］ 朱继洲.核反应堆安全分析［M］.西安：西安交通大学出版社，2011：70-73.

［8］ 陶文铨.数值传热学［M］.西安：西安交通大学出版社，2001：162-176.

［9］ 李庆扬，王能超，易大义.数值分析［M］.北京：清华大学出版社，2008：283-284.

Simulation Research on Decay Heat Removal System in Primary Loop of Pool-type Sodium-cooled Fast Reactor

JIANG Bo，ZHANG Zhi-gang＊，YU Yang，CHEN Guang-liang，ZHANG Zhi-jian
（Fundamental Science on Nuclear Safety and Simulation Technology Laboratory，Harbin Engineering University，Harbin150001，China）

The decay heat removal system in pool-type sodium-cooled fast reactor（PSFR）is the passive safety system，which depends on the natural circulation of sodium and air to keep the reactor coolant cooled.In order to verify the characteristics of the heat transfer of decay heat removal system in primary loop for accident condition，the core single-channel model and the flow between fuel assemblies model were established to simulate the decay heat removal system of primary loop and testify the program on station blackout accident，by using fully-implicit second-order upwind scheme and ameliorative Eular method to solve the equations based on FORTRAN.The calculation results show that the program could reflect the transient characteristics of the decay heat removal system，and it could reach excess real-time simulation.

decay heat removal system；natural circulation；flow between fuel assemblies model；numerical simulation

TL43

：A

：1000-6931（2015）05-0863-08

10.7538／yzk.2015.49.05.0863

2014-01-17；

2014-05-05

国家自然科学基金资助项目（51476040，11105034）；核能开发项目资助；第46批教育部留学回国人员科研启动基金资助项目

姜 博（1991—），女，辽宁本溪人，硕士研究生，核能科学与工程专业

＊通信作者：张智刚，E-mail：zzgbest＠hotmail.com