加速器驱动次临界系统束流瞬态分析模型的开发

加速器驱动次临界系统束流瞬态分析模型的开发
陈钊;崔大伟;石秀安
【摘要】加速器驱动次临界系统利用散裂反应产生外源中子驱动次临界堆运行,具有次临界固有安全性,同时具备能谱硬、嬗变能力强等特点,被国际公认为核废料处理的最有效手段.ADS系统中外中子源由质子束流轰击散裂靶产生,束流的瞬态变化将直接引起次临界堆堆芯功率的波动,从而影响整个ADS系统的安全运行.本文在调研分析国际现有的ADS束流瞬态分析模型的基础上,提出一种新型的ADS束流瞬态分析模型.基于通用CFD程序FLUENT,通过用户自定义功能(UDF)将中子动力学模型(PKM)和燃料棒瞬态热分析模型(PTM)集成进入FLUENT软件中,完成FLUENT-ADS束流瞬态分析模型开发.采用OECD/NEA发布的ADS失束事故国际基准例题进行模型验证,关键校验参数与发布结果吻合较好,最大计算误差为
5.2％,与国际同类功能的计算程序相当,模型具有一定的可信度,可满足ADS束流瞬态特性初步分析研究要求.
【期刊名称】《核安全》
【年(卷),期】2018(017)004
【总页数】8页(P51-58)
【关键词】ADS;束流瞬态;模型开发
【作者】陈钊;崔大伟;石秀安
【作者单位】中广核研究院有限公司,深圳518026;中广核研究院有限公司,深圳518026;中广核研究院有限公司,深圳518026
【正文语种】中文
【中图分类】TL99
加速器驱动次临界系统(Accelerator Driven Sub-critical System, 简称ADS)利用散裂反应产生外源中子驱动次临界堆运行，具有次临界固有安全性，同时具备能谱硬、嬗变能力强等特点，是国际公认的核废料处理最有效手段之一[1,2]。

ADS系
统的外源中子由质子束流轰击散裂靶产生，束流的瞬态变化将直接引起次临界堆堆芯功率的波动，从而影响整个ADS系统的安全运行[3]。

因此，开展ADS束流瞬
态特性的模型分析研究对发展ADS技术，尤其是开展安全分析和安全审评工作具有重要意义。

国际上ADS束流瞬态分析研究主要集中于德国KIT、意大利ENEA和比利时SCK·CEN等单位[3-5]。

在我国，中科院合肥物质科学研究院、中国科学技术大学等单位在中科院ADS嬗变系统战略性先导科技专项的支持下也开展了相关研究[6-9]。

2015年，作为国家重大科技基础设施的“加速器驱动嬗变研究装置(简称CiADS)”已由国家发改委批准立项，致力于开展国际领先的ADS技术研究。

但是，由于国际上尚未有建成的ADS系统可提供实验数据验证，现阶段关于ADS束流
瞬态特性的分析模型和方法均处于理论研究阶段，急需开发多种不同原理的分析方法和模型进行相互校验，提高计算模型的可信度。

本文在调研分析国际上现有的分析模型和方法的基础上，提出一种新型的ADS束流瞬态分析模型。

所开发的模型基于通用计算流体动力学(Computational Fluid Dynamics，简称CFD)程序FLUENT，通过用户自定义功能(User Defined Functions, 简称UDF)将中子动力学模型(Point Kinetic Model，简称PKM)和燃
料棒瞬态热分析模型(Pin Thermal Model，简称PTM)集成进入FLUENT软件中，
从而实现ADS束流瞬态分析功能，并采用国际基准例题进行校验计算，开展模型适用性和准确性的评估。

与目前通用的基于集总参数法的系统程序相比，本文提供的基于通用CFD方法的ADS束流瞬态分析模型和方法具有更加灵活的一回路结构建模能力，能够更加精确地模拟一回路结构中复杂多维的热工水力现象，具有较好的创新性，可为ADS系统的设计和安全审评工作提供新型分析工具，具有重要的研究意义和应用价值。

1 ADS系统概述
ADS的基本原理是利用加速器产生的高能质子束流轰击靶材引起散裂反应，以散裂中子作为外源中子驱动次临界堆芯中的核嬗变反应，如图1所示[1]。

从图1中可以看出，ADS系统由三大部分组成，分别是质子加速器、散裂靶和次临界堆。

在一般的设计方案中，次临界堆多采用铅或铅铋作为冷却剂。

本文以铅铋冷却ADS系统作为研究对象。

图1 ADS系统原理图Fig.1 Schematic diagram of ADS
2 模型开发
ADS束流瞬态过程中，堆芯区域发生极为复杂的物理-热工耦合现象，其理论分析模型需同时考虑三种物理模型的耦合作用，分别是：带有外源驱动的中子动力学模型(PKM)、带有内热源的燃料棒瞬态热传导模型(PTM)和一回路冷却剂热工水力模型(CFD)。

本文采用的多物理耦合分析方法是：以CFD程序为主程序计算液态铅铋冷却剂热工力水力学参数，通过软件用户自定义工具FLUENT-UDF将具有外源驱动的中子动力学模型和具有内热源的燃料棒瞬态热传导模型添加到FLUENT软件中，建立相应的数据传递通道，从而实现ADS束流瞬态分析计算。

本文采用的数据传递结构如图2所示。

图2 计算模型数据传递结构Fig.2 Data transferring structure of simulation
model
2.1 中子动力学模型
中子动力学模型(PKM)用于计算ADS堆芯功率，采用的是具有六组缓发中子的含
外中子源的点堆动力学模型以及考虑了燃料多普勒效应和冷却剂密度等反应性反馈模型。

为了简化计算，本文采用的点堆动力学模型尚未考虑外源项的空间分布效应，该效应将在后续研究中进行考虑[10,11]。

本文采用的具有六组缓发中子的点堆动
力学模型的数学描述如下：
至6)
(1)
式中：n——中子密度；
βi——缓发中子的份额(第i组)；
缓发中子的总份额；
Ci——缓发中子先驱核浓度(第i组)；
λi——缓发中子的衰变常数(第i组)；
ρ(t)——反应性；
Λ——中子代时间；
S(t)——外中子源。

反应性反馈模型的数学描述如下：
·[Tc(t)-Tc(0)]+Ka·[Tf(t)-Tf(0)]
(2)
式中：ρext ——外部反应性；
KD——燃料多普勒常数；
Kc——冷却剂反馈系数；
Kr——堆芯径向膨胀反馈系数；
Kf——堆芯轴向膨胀反馈系数；
Tf——燃料温度；
Tc——冷却剂温度。

外中子源与次临界度的关系表达式采用的是C. Rubbia等人给出的计算模型[12]：
(3)
式中，k0——稳态时的keff；
n0——稳态时的中子密度。

在数值求解方法方面，PKM采用采用后退欧拉有限差分法(Backward Euler finite difference scheme，简称BEFD)[13]求解点堆动力学方程。

2.2 燃料棒瞬态热分析模型
燃料棒瞬态热分析模型(PTM)用于计算燃料棒瞬态温度分布，采用了单通道具有内热源的圆柱体瞬态热传导模型，如图3所示。

图3 燃料棒热分析模型Fig.3 The thermal model of PTM
由于燃料棒的轴线导热影响较小，为简化计算，本文采用了忽略轴向热传导的具有内热源圆柱体瞬态热传导方程，其数学描述如下：
(4)
式中，ρm——密度；
cp,m——比热容；
km——热导率；
t——时间；
T——温度；
QV——内热源体功率密度。

下标m代表材料，f，g，和cl分别代表燃料材料，间隙材料和包壳材料。

方程(4)的边界条件为：
式中，Tcoolant——冷却剂温度；
Tcl——包壳外表面温度；
R——半径；
h——传热系数。

下标f，g，cl分别代表燃料，间隙和包壳。

在数值求解方法方面，PTM采用Crank-Nicholson差分格式[14]建立代数方程组，采用三对角矩阵算法(The tridiagonal matrix algorithm，简称TDMA)[15]求解
代数方程组。

2.3 冷却剂热工水力分析模型
冷却剂热工水力分析模型采用CFD方法进行模拟。

CFD计算区域可根据具体的研究对象而定，对于堆内结构呈对称性布置的反应堆，可采用二维轴对称方法进行简化，而对于堆芯结构布置不具有对称性的反应堆，则需要采用三维模拟。

对于采用池式布置的典型ADS次临界堆，其二维轴对称计算模型如图4所示。

图4 采用池式布置的典型ADS次临界堆二维轴对称CFD计算模型Fig.4 CFD simulation model of a typical pool-type ADS subcritical reactor
由于反应堆堆芯区域几何结构极为复杂，采用CFD方法进行精细模拟需要花费巨
大的计算资源，因此需要采用合理的方法进行模型简化。

本文参考SIMMER等程序，采用宏观等效模拟方法进行模型简化[16,17]，具体如下。

采用等效流通面积方法将具有六边形结构或四边形结构的燃料组件等效为圆环型结
构，如图5所示；采用多孔介质方法将由于冷却剂流动造成的压降等效为动量守恒方程中的外动量源，如式(5)所示。

通过多孔介质方法，在CFD软件中建立起与雷诺数相关的堆芯压降计算方法。

图5 CFD堆型简化模型Fig.5 Simplified CFD core model
(5)
式中，Si——外动量源；
νi——流速；
μ(T)——冷却剂动力黏度；
ρ(T)——冷却剂密度；
Di——黏性系数；
Ci——惯性系数。

2.4 数据传递结构
基于通用CFD软件FLUENT，通过用户自定义功能(UDF)将PKM和PTM集成进入FLUENT软件中，形成FLUENT-ADS耦合分析模型。

FLUENT-ADS模型计算流程图如下：
图6 FLUENT-ADS程序计算流程图Fig.6 The calculation flowchart of FLUENT-ADS
3 模型验证
为了验证FLUENT-ADS模型的适用性和准确性，本文采用OECD/NEA发布的ADS失束事故国际基准例题进行模型验证，校验参数包括功率、燃料芯块中心温度、燃料芯块外表面温度和冷却剂活性区出口温度四项重要参数。

在OECD/NEA ADS失束事故基准例题研究项目中，共有9家研究机构参加并提交了计算结果，具体信息见表1[4]。

在OECD基准例题中，大部分学者采用的方法是在基于集总
参数法的系统程序上进行ADS模型开发，如SIM-ADS、LOOP2、SAS4ADS、TRAC-M等，而本文采用的是基于通用CFD方法进行ADS模型开发，具有更加灵活的一回路结构建模能力，能够更加精确地模拟一回路结构中复杂多维的热工水力现象。

表1 OECD/NEA ADS失束事故基准例题项目参与单位Table 1 Participant organisations of OECD/NEA ADS beam trip benchmark参与方国家程序
A.DAngeloENEA (意大利)TIESTE-MINOSSEG. Van den Ey-nde
B. ArienSCK·CEN (比利时)SITHER-PKSK. TsujimotoJAERI (日本)EXURS-MM. ErikssonKTH(瑞典)SASSYS/SAS4AM. SchikorrKIT(原FZK,德国)SIM-ADSP.
C. Konstantin PSI (瑞士)TRAC-MM. P. Codding-ton and K. Miki-tyukPSI (瑞士)LOOP2P. Wakker and J. KuijperNRG (荷兰)TRAC MODR. Dagan and C. BroedersKIT(原FZK, 德国)SAS4ADSY. KimKAERI ( 韩国)DESINUR
3.1 基准例题
基准例题的分析对象选用的是欧盟第五框架协议项目下开发的80 MWth铅铋冷却ADS概念设计方案(80 MWth-XADS)，主要设计参数见表2 [18]。

基准例题所分析瞬态条件为在满功率稳态运行的情况下，将外源质子束流断开，断束时间分别为1s，3s，6s，12s，随后恢复束流，观察不同断束时间情况下的功率、燃料温度的变化。

表2 80 MWth-XADS主要设计参数Table 1 Main design parameters of 80 MWth-XADS设计参数数值热功率80MWth冷却剂液态铅铋堆芯温度(进口/出口)300℃/400℃燃料MOXkeff0.97多普勒常数-0.00580缓发中子份额0.00314中子代时间1.43μs径向膨胀反馈系数-0.00000499轴向膨胀反馈系数-
0.00000157
3.2 验证结果
图7给出的是FLUENT-ADS计算结果中反应堆功率变化曲线与OECD/NEA发布计算结果的对比情况，图8为局部放大后的曲线。

从图7和图8中可以看出，堆
芯功率在断束后迅速下降至10%，随后开始缓慢下降，当束流恢复后，堆芯功率
重新恢复。

值得注意的是，由于负反馈的作用，此时的堆芯功率会略超过初始功率，然后缓慢下降，最终恢复为初始功率。

FLUENT-ADS计算结果整体趋势与基准例
题发布的其他计算结果吻合较好。

图7 不同失束事故瞬态下的功率变化曲线Fig.7 Power changing curve under various beam interruption accidents
图8 不同失束事故瞬态下的功率变化曲线(局部放大)Fig.8 Power changing curve under various beam interruption accidents (partly enlargement)
图9和图10给出的是FLUENT-ADS计算结果中燃料芯块中心温度和燃料表面温
度与OECD/NEA发布计算结果的对比情况。

从图9和图10中可以看出，当断开外源中子后，燃料中心温度和燃料表面温度均迅速下降，在束流恢复后又重新返回初始值。

FLUENT-ADS结果中，失束1s，3s，6s，12s时燃料芯块中心温度与OECD/NEA发布计算结果的最大偏差分别为-3.61%，0.29%，5.241%，0.9%，见表3。

图9 不同失束事故瞬态下燃料芯块中心温度变化曲线Fig.9 Temperature changing curve of fuel pellet centre under various beam interruption accidents
图10 不同失束瞬态事故下燃料芯块外表面温度变化曲线Fig.10 Temperature changing curve of pellet surface under various beam interruption accidents 表3 燃料芯块中心温度计算误差Table 3 Calculation error of fuel pellet centre temperature
失束时间/s燃料芯块中心温度/KOECD/NEA发布计算结果(各程序平均
值)FLUENT-ADS结果相对误差/%17067.5-
3.63180180.50.3630028
4.7
5.212370373.30.9
图11给出的是FLUENT-ADS计算结果中堆芯活性区冷却剂出口温度与
OECD/NEA发布计算结果的对比情况。

从图11中可以看出，当断开外源中子后，堆芯活性区冷却剂出口温度均迅速下降，在束流恢复后又重新返回初始值。

FLUENT-ADS结果中，失束1s，3s，6s，12s时堆芯活性区冷却剂出口温度与OECD/NEA发布计算结果的最大偏差分别为1.68%，4.56%，3.06%，2%，见表4。

图11 不同失束瞬态事故下堆芯活性区冷却剂出口温度变化曲线Fig.11 Outlet coolant temperature under various beam interruption accidents表4 堆芯活性区冷却剂出口温度计算误差Table 4 Calculation error of outlet coolant temperature
失束时间/s堆芯活性区冷却剂出口温度/KOECD/NEA发布计算结果(各程序平均值)FLUENT-ADS结果相对误差
/%11313.2181.733536.5964.666061.8363.1128586.7042
4 结论
本文完成基于通用CFD方法的ADS束流瞬态分析模拟初步开发研究，通过用户
自定义功能(UDF)将中子动力学模型(PKM)和燃料棒瞬态热分析模型(PTM)集成进
入通用CFD程序FLUENT中形成FLUENT-ADS分析模型，并采用OECD/NEA
发布的国际基准例题进行模型验证。

FLUENT-ADS模型验证结果显示，关键校验
参数发布结果吻合较好，最大计算误差为5.2%，与国际同类功能的计算程序相当，模型具有一定的可信度，可满足ADS束流瞬态的初步分析研究。

与国际通用的基于集总参数法的系统程序分析方法相比，本文提供的基于通用CFD方法的ADS
束流瞬态分析模型和方法具有更加灵活的一回路结构建模能力，能够更加精确地模
拟一回路结构中复杂多维的热工水力现象，对于深入研究ADS束流瞬态机理现象具有重要意义，同时可为我国ADS系统的设计和安全审评工作提供新型分析研究工具，具有工程应用价值。

参考文献
【相关文献】
[1]詹文龙,徐瑚珊.未来先进核裂变能——ADS嬗变系统[J].中国科学院院刊, 2012, 27(3):375-381.
[2]赵志祥. 加速器驱动的洁净核能系统国际研究进展 [J]. 原子核物理评论, 1997, 14(2):121-124.
[3]Konings R J , Marucci G, Abram T, et al.A european roadmap for developing accelerator driven systems (ADS) for nuclear waste incineration [M]. Roma:ENEA, 2011.
[4]D’Angelo A, Gabrielli F. Finalreport on phase I calculations - benchmark on beam interruptions in an accelerator driven System [R].Paris:OECD Publications,2003.
[5]Arien B, Heusdains S, Abderrahim A H, et al. Safety analysis of the MYRRHA facility with different core configurations [C]. Canada:PHYSOR, 2006.
[6]冯竟超.先进反应堆热工水力与多物理场耦合程序开发及应用研究 [D].合肥:中国科学技术大学, 2017.
[7]Wang Zhen,Wang Gang,Gu Zhixing,et al. Benchmark of neutronics and thermal-hydraulics coupled simulation program ntc on beam interruptions in XADS[J]. Annals of Nuclear Energy, 2015,77:172-175.
[8]王强龙.加速器驱动次临界堆始发事件及非能动安全特性初步研究[D].合肥:中国科学技术大学, 2013.
[9]Rineiski A,Maschek W. Kinetics models for safety studies of accelerator driven systems[J]. Annals of Nuclear Energy,2005,32(12):1348-1365.
[10]沈峰,王苏. ADS次临界反应堆的点堆中子动力学方程[J]. 原子能科学技术, 2011, 45(11):1300-1304.
[11]王苏,沈峰.ADS次临界反应堆的中子共轭方程[J].原子能科学技术, 2011, 45(7):775-779.
[12]Rubbia C, Rubio J A, Buono S, et al. Conceptual design of a fast neutron operated high power energy amplifier, Tech. Rep. CERN/AT/95-44 (ET) [R].Meyrin:European Organization for Nuclear Research , 1995.
[13]Nahla A A. An efficient technique for the point reactor kinetics equations with newtonian temperature feedback effects[J]. Annals of Nuclear Energy,2011, 38(12):2810-
2817.
[14]顾樵. 数学物理方法[M].北京：科学出版社, 2012.
[15]陶文铨. 数值传热学(第2版)[M]. 西安：西安交通大学出版社, 2001.
[16]Suzuki T,Chen Xuenong,Rineiski A,et al. Transient analyses for accelerator driven system pds-xads using the extended simmer-iii code[J]. Nuclear Engineering and Design,2005,235(24):2594-2611.
[17]Chen Zhao,Chen Xuenong,Rineiski A, et al. Coupling a CFD code with neutron kinetics and pin thermal models for nuclear reactor safety analyses[J]. Annals of Nuclear Energy, 2015,83: 41-49.
[18]Luciano C,Giuseppe G. The pb-bi cooled xads status of development[J]. Journal of Nuclear Materials, 2002,301(1):8-14.。