高慶瑜，宋英明，徐宇超，王 珂，楊永偉，張 璐

(1.南華大學，湖南衡陽421001；2.中國科學院近代物理研究所散裂靶室，甘肅蘭州730000)

ADS次臨界系統(tǒng)中子時空動力學計算與瞬態(tài)分析

高慶瑜1，宋英明1，徐宇超1，王 珂1，楊永偉2，張 璐2

(1.南華大學，湖南衡陽421001；2.中國科學院近代物理研究所散裂靶室，甘肅蘭州730000)

加速器驅(qū)動次臨界反應堆(ADS)中子時空動力學計算需要考慮外中子源和空間分布的影響，比臨界系統(tǒng)中子動力學計算要復雜得多。本文將改進準靜態(tài)(IQS)近似與蒙特卡羅(MC)方法相結(jié)合，對于帶外源的ADS次臨界系統(tǒng)中子時空動力學過程，形狀函數(shù)、動力學參數(shù)由MCNPX程序計算得到，幅度函數(shù)與集總參數(shù)熱工反饋模型進行耦合計算，并開發(fā)了IQS/MC計算程序可視化操作界面。針對CIADS靶堆耦合系統(tǒng)參考方案物理模型，對引入束流瞬變及無保護失流工況過程進行瞬態(tài)模擬計算分析，給出了堆芯相對功率、燃料溫度及冷卻劑出口溫度隨時間的變化曲線。同時，將中子注量率進行分群計算，得到了堆芯分能群的相對中子注量率網(wǎng)格分布隨時間的變化，模擬結(jié)果與理論分析一致。

ADS次臨界系統(tǒng)；中子時空動力學；改進準靜態(tài)近似；模擬計算；瞬態(tài)分析

加速器驅(qū)動次臨界系統(tǒng)(ADS)主要由散裂靶、加速器質(zhì)子束流及次臨界反應堆三者構(gòu)成，對比傳統(tǒng)的壓水堆，由于外中子源的影響，造成了ADS 次臨界反應堆芯部功率的不均勻分布[1-3]。帶外中子源的ADS系統(tǒng)的中子學瞬態(tài)安全分析是衡量ADS次臨界反應堆安全性的標志之一，故其占有十分重要的地位[4,5]。在加速器驅(qū)動次臨界系統(tǒng)中，失束是一種典型的ADS束流瞬變事故，當質(zhì)子束流中斷，ADS次臨界反應堆內(nèi)中子學過程將受到影響，同時，束流瞬變也將影響反應堆的壽期[6-8]。失流事故是較為嚴重的熱工瞬變工況，在這種工況下，反應堆在熱量傳輸上的安全性會受到考驗，進而影響反應堆系統(tǒng)的安全[9,10]。本文采用蒙特卡羅方法和改進準靜態(tài)近似相結(jié)合[11]的方法，針對中科院戰(zhàn)略先導專項CIADS靶堆耦合參考方案物理模型，對引入失束和無保護失流兩種工況進行模擬計算，分析這兩種工況下CIADS靶堆耦合參考方案物理模型的瞬態(tài)安全特性。

1 方法概述

外源驅(qū)動的次臨界系統(tǒng)中，描述中子時空動力學行為的中子輸運方程[12]為：

(1)

如果略去緩發(fā)中子先驅(qū)核在介質(zhì)中的運動，則相應的與時間相關的先驅(qū)核方程為：

(2)

式中：V——中子速度；

φ——中子注量率；

Σt——總截面；

Σs——散射截面；

Σf——裂變截面；

S——外中子源；

χp——瞬發(fā)中子能譜；

χi——第i組緩發(fā)中子能譜；

β——總的緩發(fā)中子份額；

υ——每次裂變的中子產(chǎn)額；

λi——第i組緩發(fā)中子先驅(qū)核衰變常數(shù)；

Ci——第i組緩發(fā)中子先驅(qū)核濃度；

NDG——緩發(fā)中子先驅(qū)核總的組數(shù)。

(3)

并引入幅度歸一化條件：

(4)

通過整理式(1)～式(4)可得幅度函數(shù)方程組為：

(5)

i=1,…,NGD

(6)

同時，得到的形狀函數(shù)方程表達式為：

(7)

由于形狀函數(shù)隨時間的變化比較慢，沒有幅函數(shù)隨時間變化重要，因此可以沿時間軸計算很多幅值，但只算幾次形狀函數(shù)。形狀函數(shù)的計算采用改進準靜態(tài)近似，設大時間步長為[tj,tj+1]，Δt=tj+1-tj，式(7)中把ψ對時間t的偏微商換成一階向后差分：

(8)

式(8)相當于具有非齊次方程形式的中子輸運方程[13,14]。

基于蒙特卡羅改進準靜態(tài)(IQS/MC)方法的中子時空動力學模擬流程如下：

第1步，利用MCNPX程序進行初始穩(wěn)態(tài)計算，獲得初始時刻堆芯能譜分布、歸一化形狀函數(shù)分布及計算所需的中子動力學參數(shù)。第2步，進入IQS/MC流程：將由MCNPX程序計算給出的動力學參數(shù)進行處理并輸入IQS/MC動力學程序(含集總參數(shù)熱工模型[15,16])，在引入瞬態(tài)工況過程或其他控制情況下進行小步長計算得到幅度函數(shù)，同時可得燃料溫度及冷卻劑出口溫度隨時間的變化。第3步，利用上一時間步長內(nèi)的形狀函數(shù)及其一階導數(shù)來預估下一時間步長內(nèi)的形狀函數(shù)。第4步，將形狀函數(shù)方程的非齊次項作為外源輸入MCNPX程序，計算得到本步長的形狀函數(shù)和中子動力學參數(shù)。第5步，將IQS/MC動力學程序計算得到的幅度函數(shù)與MCNPX程序計算得的相應時刻的形狀函數(shù)相乘，得到本步長內(nèi)堆芯相對中子注量率分布，同理可得到相對功率三維分布。第6步，重復1～5步，進行下一步長內(nèi)幅度函數(shù)與形狀函數(shù)的耦合計算至模擬結(jié)束。

2 物理模型及參數(shù)

給定的CIADS參考方案模型堆芯軸向結(jié)構(gòu)見圖1，該次臨界堆芯模型包含燃料組件、反射層組件、屏蔽層組件，中心區(qū)域為可發(fā)射散裂中子源的靶區(qū)。

圖1 MCNP堆芯模型Fig.1 Core model by MCNP

中子時空動力學瞬態(tài)計算過程，中子分4個能群描述，其量范圍和平均速率見表1。

表1 各能群能量范圍

進行瞬態(tài)分析模擬前，需要對CIADS靶堆耦合參考方案物理模型進行穩(wěn)態(tài)工況計算。表2給出反應堆在穩(wěn)態(tài)運行時由MCNPX程序計算[17,18]得出的中子動力學參數(shù)結(jié)果。

表2 中子動力學參數(shù)結(jié)果

幅度函數(shù)計算與含集總參數(shù)熱工反饋模型耦合，穩(wěn)態(tài)計算時，熱工參數(shù)選取見表3。

表3 堆芯熱工參數(shù)的選取

3 結(jié)果及分析

在反應堆穩(wěn)態(tài)運行的基礎上對失束和無保護失流兩種瞬態(tài)工況算例進行模擬計算。

3.1 失束工況計算結(jié)果與分析

假設t=0時刻，外源突然中斷，模擬0～10 s堆芯的瞬態(tài)變化。圖2～圖3給出了失束工況下反應堆相對功率、燃料溫度及冷卻劑出口溫度隨時間的變化曲線。從圖2可以看出，反應堆0s開始失束，0～4×10-5s內(nèi)堆芯相對功率驟降，4×10-5～10 s內(nèi)相對功率依然緩慢衰減。從圖3可以看出，由于失束引起功率的衰減，堆芯燃料溫度及冷卻劑出口溫度亦會逐漸降低，同時，冷卻劑出口溫度對時間的響應較燃料溫度對時間的響應慢，符合物理規(guī)律。

圖2 失束工況下相對功率隨時間的變化Fig.2 Relative power change over time under beam interruption condition

圖3 失束工況燃料溫度和冷卻劑出口溫度隨時間變化Fig.3 Fuel temperature and outlet temperature of coolant change over time under beam interruption condition(a) 燃料溫度隨時間的變化

圖3 失束工況燃料溫度和冷卻劑出口溫度隨時間變化(續(xù))Fig.3 Fuel temperature and outlet temperature of coolant change over time under beam interruption condition(b) 冷卻劑出口溫度隨時間的變化

圖4、圖5給出了失束工況下第一能群及第三能群相對中子注量率網(wǎng)格分布隨時間的變化，從圖中可以看出，相對中子注量率變化趨勢與圖2中相對功率變化趨勢一致，符合規(guī)律。

3.2 無保護失流工況計算結(jié)果與分析

假設廠房突然斷電，冷卻劑主泵全部停運，冷卻劑流量迅速下降。模擬0～100 s內(nèi)堆芯的瞬態(tài)變化。圖6、圖7給出了無保護失流工況下反應堆相對功率、燃料溫度及冷卻劑出口溫度隨時間的變化曲線。從圖6中可以看出，相對功率隨冷卻劑流量降低而減小，而且流量變化越大，相對功率變化越激烈。從圖7中可以看出，當冷卻劑流量降低時，燃料溫度和冷卻劑出口溫度呈上升趨勢，然后逐漸趨于穩(wěn)定。冷卻劑流量變化越大，燃料溫度及冷卻劑出口溫度隨時間的變化越劇烈，而且，冷卻劑溫度對時間的響應慢于燃料溫度對時間的響應，符合物理規(guī)律。

圖4 第一能群(0.400

圖5 第三能群(2.154×10-3

圖6 無保護失流工況下相對功率隨時間的變化Fig.6 Relative power change over time under loss-of-flow without scram condition

選取冷卻劑流量為穩(wěn)態(tài)時的30%這一工況模擬分能群相對中子注量率網(wǎng)格分布隨時間的變化。圖8、圖9給出了相應的第二能群及第四能群t=0時刻相對中子注量率三維網(wǎng)格分布情況。對于第二能群，t=0s，30s，60s，90s時，相對中子注量率峰值依次為0.0514，0.0512，0.0511，0.0510。對于第四能群，t=0s，30s，60s，90s時，相對中子注量率峰值依次為0.0388，0.0386，0.0386，0.0385。由上述結(jié)果可以看出，在無保護失流工況下相對中子注量率變化趨勢與圖6中相對功率變化趨勢一致，符合物理規(guī)律。

圖7 失束工況燃料溫度和冷卻劑出口溫度隨時間變化Fig.7 Fuel temperature and outlet temperature of coolant change over time under beam interruption condition(a) 燃料溫度隨時間的變化；(b) 冷卻劑出口溫度隨時間的變化

圖8 第二能群(4.642×10-2

圖9 第四能群(0

4 結(jié)語

本文將改進準靜態(tài)(IQS)近似與蒙特卡羅(MC)方法相結(jié)合，對引入失束和無保護失流兩種瞬態(tài)工況進行了模擬，分析這兩種工況下CIADS靶堆耦合參考方案物理模型的瞬態(tài)安全特性。文中給出了堆芯總的相對功率、燃料溫度及冷卻劑出口溫度隨時間的變化曲線，同時將中子注量率進行分能群描述，給出了分能群相對中子注量率網(wǎng)格分布隨時間的變化。由計算結(jié)果可知，在ADS中子時空動力學計算中，當外源變化時，堆芯相對功率能夠?qū)ζ渥鞒鏊查g的響應，而燃料溫度和冷卻劑出口溫度隨時間的變化趨勢基本相同，但冷卻劑出口溫度對時間的響應慢于燃料溫度對時間的響應。同時須指出，當堆芯冷卻劑流量變化微小時，反應堆仍能正常運行，但是當冷卻劑流量發(fā)生較大變化時，燃料溫度和冷卻劑出口溫度瞬間上升，所以為保護反應堆安全須采取一定措施。

[1] 詹文龍, 徐瑚珊. 未來先進核裂變能——ADS嬗變系統(tǒng)[J].中國科學院院刊,2012,27(3)：375-381.

[2] 陳森. 加速器驅(qū)動鉛鉍冷卻自然循環(huán)次臨界堆無保護瞬態(tài)分析研究[D].合肥： 中國科學技術大學,2014.

[3] Gandini A, Salvatores M. The Physics of Sub-critical Multiplying System [J]. Journal of Nuclear Science and Technology，2002，39(6)：673-686.

[4] 付龍舟.反應堆穩(wěn)定性分析[J].原子能科學技術,1964,02：144-151.

[5] 陳森, 金鳴, 陳志斌, 等.鉛鉍冷卻次臨界堆反應性引入下的中子學動態(tài)特性分析[J]. 核安全, 2014, 13(2)：45-49.

[6] Liu P, Chen X N, Rineiski A, et al. Transient analysis of the 400 MWth-class EFIT accelerator driven transmuter with the multi-physics code： SIMMER-III[J]. Nuclear Engineering and Design, 2010, 240(10)： 3481-3494.

[7] Sugawara T, Nishihara K, Tsujimoto K. Transient analyses for lead-bismuth cooled accelerator-driven system[J]. Annals of Nuclear Energy, 2013, 55： 238-247.

[8] 汪振,王剛,辜峙钘,等.加速器驅(qū)動鉛鉍冷卻自然循環(huán)次臨界堆束流瞬變事故研究[J].核技術, 2015, 38(1)：91-96.

[9] Shan J, Wang W, Wang X, et al. Loss of Flow without Scram and Safety Analysis of LMFBR[J]. Chinese Journal of Nuclear Science & Engineering, 2002, 22(2)：113-118.

[10] Sutanto, Oka Y. Analysis of anticipated transient without scram of a Super Fast Reactor with single flow pass core[J]. Annals of Nuclear Energy, 2015, 75：54-63.

[11] 徐琪, 王侃, 李世銳, 等. 反應堆時空動力學蒙特卡羅改進準靜態(tài)方法研究[J].原子能科學技術, 2013, 47(s1)： 275-279.

[12] 黃祖洽. 核反應堆動力學基礎[M]. 北京大學出版社, 2014.

[13] 宋英明,馬遠樂,周志偉. 高溫氣冷堆中子時空動力學實時仿真計算[J]. 原子能科學技術,2010, 44(02)：188-192.

[14] 宋英明,周志偉,馬遠樂. 模塊式高溫氣冷堆三維中子動力學計算分析[J]. 核科學與工程,2010, 30(02)：135-143.

[15] Cammi A, Luzzi L, Porta AA, et al. Modelling and control strategy of the Italian LBE-XADS[J]. Progress in Nuclear Energy, 2006,48(6)：578-589.

[16] Ponciroli R, Cammi A, Lorenzi S, et al. A preliminary approach to the ALFRED reactor control strategy[J]. Progress in Nuclear Energy, 2014, 73(5)：113-128.

[17] 張良, 陳偉, 趙柱民, 等. 瞬發(fā)中子密度衰減法計算中子代時間[J]. 強激光與粒子束, 2013, 25(1)： 237-240.

[18] Rubbia C. ACCELERATOR-DRIVEN NUCLEAR SYSTEM WITH CONTROL OF EFFECTIVE NEUTRON MULTIPLICATION COEFFICENT： EP, EP2561514[P]. 2014.

CalculationandTransientAnalysisofNeutronTime-spaceKineticsforAcceleratorDrivenSub-criticalSystem

GAOQing-yu1，SONGYing-ming1，XUYu-chao1，WANGKe1，YANGYong-wei2，ZHANGLu2

(1. Institute of Nuclear Science and Technology, University of South China, Hengyang of Hunan Prov. 421001, China；2. Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou of Gansu Prov. 730000, China)

Calculation of neutron space-time kinetic for Accelerator driven sub-critical reactor (ADS) need to consider the impact of external neutron source and spatial distribution, and it was much more complex than critical system in neutron kinetics calculation. In this paper, the improve quasi-static approximation (IQS) and Monte Carlo (MC) method were combining, the shape function and kinetic parameters were calculated by MCNPX program, and the amplitude function was coupled with the lumped parameter thermal feedback model to calculate for neutron space-time kinetic process of ADS critical system with external source, besides, a visual operating interface of IQS/MC calculation procedure was developed. For CIADS physical model of target and reactor coupled reference scheme, the transient simulation analysis was carried out on inserting beam interruption and loss-of-flow without scram, and the core relative power, fuel temperature and outlet temperature of coolant along time were given. At the same time, multigroup neutron flux were calculated, and the core gird distribution of multigroup relative neutron flux along time were obtained, and simulation results were consistent with theoretical analysis.

Accelerator-driven sub-critical system； Neutron time-space kinetics； Improved quasi static approximation；Simulation calculation；Transient analysis

2016-12-20

高等學校博士學科點專項科研基金(20134324120003); 湖南省教育廳重點項目基金(14A120);中國科學院戰(zhàn)略性先導專項“ADS嬗變系統(tǒng)”子課題(XDA03030102); 南華大學研究生科研創(chuàng)新項目(2016XCX21)

高慶瑜(1991—)，女，內(nèi)蒙古人，碩士研究生，現(xiàn)從事核能科學與工程方面研究

宋英明：songyingming@tsinghua.org.cn

TL327

：A

：0258-0918(2017)04-0532-08