石墨慢化通道式熔鹽堆有效緩發(fā)中子份額計算方法研究

崔 勇，陳金根,*，何 帆，李曉曉，蔡翔舟,*

(1.中國科學院 上海應(yīng)用物理研究所，上海 201800； 2.中國科學院 先進核能創(chuàng)新研究院，上海 201800；3.中國科學院大學，北京 100049)

熔鹽堆是6種第4代核能系統(tǒng)候選堆型中唯一使用液態(tài)燃料的反應(yīng)堆，與固態(tài)燃料堆不同，熔融的氟化鹽混合物(LiF，BeF2，ZrF4，UF4或ThF4)兼作冷卻劑與燃料，充滿整個一回路(含堆芯)。燃料流動使得緩發(fā)中子先驅(qū)核(DNP)在堆內(nèi)的分布受到流場的影響。相應(yīng)地，緩發(fā)中子的價值也會發(fā)生改變。部分具有較長半衰期的DNP會隨燃料流出堆外并衰變[1]，從而導致有效緩發(fā)中子份額(βeff)降低。其次，裂變能直接釋放于冷卻劑中，因此傳統(tǒng)的流動-傳熱計算方法不再適用。最后，熔鹽堆適合采用釷鈾燃料循環(huán)[2]，因而易裂變材料中可能包含233U和錒系核素(MA)[3]，這些重核的緩發(fā)中子數(shù)據(jù)與采用鈾钚燃料循環(huán)的固態(tài)堆不同。這些特殊性使得熔鹽堆βeff的計算較為復雜。

βeff的精確計算對研究熔鹽堆的瞬態(tài)安全特性具有重要意義。如在失流事故中，緩發(fā)中子源項的再分布會向堆內(nèi)引入正反應(yīng)性，該反應(yīng)性的量化需對不同流動工況下的βeff進行評估[3]；在反應(yīng)性引入事故中，需確定所引入反應(yīng)性與βeff的關(guān)系以得到反應(yīng)堆的響應(yīng)特性[4-5]。近年來，國際國內(nèi)針對DNP的對流[6]、擴散[7]和湍流輸運效應(yīng)[8]對熔鹽堆穩(wěn)態(tài)與瞬態(tài)特性的影響開展了廣泛研究。對于罐式熔鹽快堆，由于其堆內(nèi)無慢化材料，可采用解析方法[3]、數(shù)值方法[3,8-9]和蒙特卡羅方法[3]計算其βeff與DNP分布。對于通道式熔鹽堆，由于存在石墨、氫化鋯等慢化材料，其流動-傳熱方式具有特殊性[5,10]，對快堆提出的βeff計算方法不再適用。

本文利用解析方法和基于四階多項式節(jié)塊展開法與并聯(lián)多通道模型的三維核熱耦合程序TMSR3D[11]，對熔鹽實驗堆(MSRE)進行穩(wěn)態(tài)計算，得到堆內(nèi)DNP分布情況，研究入口流量、燃料的堆外流動時間等參數(shù)對βeff的影響。

1 計算模型與方法

1.1 計算模型

MSRE堆芯與組件幾何結(jié)構(gòu)[12]如圖1所示。堆芯由正方形石墨組件構(gòu)成，組件四周開孔作為燃料鹽流道，所有孔道與上下腔室連通。堆芯四周為石墨反射層。表1、2分別列出MSRE的主要設(shè)計參數(shù)[4,12]和裝載235U燃料時的緩發(fā)中子數(shù)據(jù)[4]。詳盡的堆芯設(shè)計與熱工參數(shù)參見文獻[4]。

表1 MSRE設(shè)計參數(shù)Table 1 Design parameter of MSRE

a——堆芯縱截面；b——堆芯橫截面；c——組件幾何圖1 MSRE堆芯與組件幾何結(jié)構(gòu)Fig.1 Core and assembly configuration of MSRE

1.2 計算方法

熔鹽堆穩(wěn)態(tài)下中子學計算模型為：

(1)

(2)

1.2.1解析方法 利用圓柱裸堆幾何對MSRE的臨界通量進行近似，記MSRE的堆芯半徑和高度分別為R和H，可得到圓柱坐標下的通量[13]為：

(3)

其中：J0為零階貝塞爾函數(shù)；φ0為待定系數(shù)，由反應(yīng)堆功率決定。

僅考慮軸向流動，記流速為u，將式(3)代入(2)中得到：

(4)

對上式在(0，z)區(qū)間積分得到：

c0exp(-λiz/u)

(5)

其中，c0為堆芯入口z=0處的DNP濃度。

設(shè)堆芯中熔鹽占堆芯總體積的比值恒為f，且堆芯出口處先驅(qū)核濃度均勻混合，其出口濃度ci，out可表示為：

(6)

設(shè)熔鹽在堆芯外回路中時間為τ，先驅(qū)核在堆外回路僅發(fā)生衰變，因此入口處先驅(qū)核濃度ci，in可表示為：

ci,in=ci,outexp(-λiτ)

(7)

由式(3)可得堆芯內(nèi)中子產(chǎn)生率N為：

sin(πz/H)2πrfdrdz

(8)

從而得到歸一化緩發(fā)中子源項Si，d為：

Si,d(r,z)=λici(r,z)/N

(9)

由于熔鹽流速u不隨半徑變化，則得到熔鹽在堆芯中的流動時間tc為：

tc=H/u

(10)

將式(5)、(7)和(8)代入式(9)中得到緩發(fā)中子源項分布：

Si,d(r,z)=α0,iexp(-λitcz/H)+

(α1,isin(πz/H)-α2,icos(πz/H)+

α2,iexp(-λitcz/H))J0(2.405r/R)

(11)

α0,i=βi·

(12)

(13)

(14)

其中，J1為一階貝塞爾函數(shù)。第i組緩發(fā)中子有效份額βi，eff為：

(15)

1.2.2數(shù)值方法 熔鹽堆穩(wěn)態(tài)G群中子擴散方程和I組DNP守恒方程[11]如下：

g=1,2,…,G;i=1,2,…,I

(16)

(17)

其中，Dg、Σt，g和Σg′g分別為g群擴散系數(shù)、總截面和散射矩陣。基于節(jié)塊展開法和有限差分方法[11]可求解得到穩(wěn)態(tài)下的中子通量分布和DNP分布。

熱工水力學計算采用并聯(lián)多通道模型[5,11,14]，包括燃料鹽的質(zhì)量、動量和能量守恒方程，即：

(18)

(19)

(20)

其中：i和j分別為通道編號和軸向節(jié)塊編號；M和N分別為軸向總節(jié)塊數(shù)和堆芯總通道數(shù)；G、A和W分別為質(zhì)量流速、等效通道流通面積和質(zhì)量流量；ρ、f和h為熔鹽的密度、摩阻系數(shù)和焓；g為重力加速度；D為通道水力直徑；Qf和Qg分別為熔鹽和石墨體積熱源；Δz和Δp分別為節(jié)塊高度和壓降。

對式(18)和(19)耦合求解[14]可得到各通道內(nèi)的流場和壓力場分布。對于單個燃料通道，石墨慢化劑可簡化為空心圓柱，只考慮徑向?qū)?，其能量守恒方程如下?/p>

(21)

其中，Tg和λg分別為石墨溫度和其導熱系數(shù)。石墨內(nèi)邊界采用熔鹽-石墨對流換熱條件，外邊界為絕熱邊界條件[5]。根據(jù)式(20)得到熔鹽溫度分布后，可解析求解式(21)得到石墨溫度分布[5]。

(22)

其中，V為積分域。

基于TMSR3D的有效緩發(fā)中子份額計算流程如圖2所示。具體為：1) 采用DRAGON5[15]計算各種組件的少群參數(shù)，并進行擬合；2) 物理熱工耦合計算；3) 共軛通量計算；4) 根據(jù)式(22)計算βeff。

圖2 基于數(shù)值方法的有效緩發(fā)中子份額計算流程Fig.2 Calculation scheme of effective delayed neutron fraction based on numerical method

2 計算結(jié)果與分析

2.1 緩發(fā)中子先驅(qū)核損失份額及其分布

額定工況下，MSRE裝載235U時DNP損失份額列于表3?？梢钥吹?，本文基于解析方法和數(shù)值方法計算所得結(jié)果與其他程序[16]均符合較好。對于第2組DNP，由于其具有較大的裂變份額和較小的衰變常量，故損失份額最大。

表3 熔鹽實驗堆額定工況下緩發(fā)中子先驅(qū)核損失份額Table 3 Loss of delayed neutron precursor fraction of molten salt reactor experiment under rated condition

基于數(shù)值方法計算得到額定工況下DNP濃度的分布，如圖3所示，由于熔鹽的流動效應(yīng)，每組DNP濃度的峰值均沿軸向有一定偏移，且流動對半衰期較長的DNP(如第2組)影響較大。這是因為半衰期越長，DNP在衰變前隨熔鹽行進的距離越長，甚至流出堆芯，造成反應(yīng)性的損失。

圖3 緩發(fā)中子先驅(qū)核濃度分布Fig.3 Concentration distribution of delayed neutron precursor

圖4示出額定工況下最熱通道內(nèi)DNP與緩發(fā)中子軸向的分布。其中，DNP濃度由數(shù)值程序TMSR3D計算得到，緩發(fā)中子分布由解析方法得到?？梢钥吹?，第2組DNP歸一化濃度高于其他各組DNP。這是由于DNP的穩(wěn)定分布不僅與裂變份額(β)有關(guān)，也與半衰期(1/λ)有關(guān)。裂變份額越大，半衰期越長，其穩(wěn)態(tài)濃度越高，即DNP濃度正比于β/λ。由于第2組DNP具有最大的β/λ值，因此其濃度高于其他各組。緩發(fā)中子由DNP經(jīng)β衰變得到，如無燃料流動效應(yīng)，穩(wěn)態(tài)下各組緩發(fā)中子的分布次序應(yīng)與DNP份額的大小順序一致。對于液態(tài)燃料熔鹽堆，還需考慮流動對緩發(fā)中子濃度的影響，如對于第5組，盡管其DNP濃度較小(圖4a)，然而其具有較短的半衰期，DNP在堆內(nèi)產(chǎn)生后短期內(nèi)即衰變?yōu)榫彴l(fā)中子，因流動導致流出堆芯損失的緩發(fā)中子較少，因而第5組緩發(fā)中子具有較大的濃度(圖4b)。

圖4 額定工況下最熱通道內(nèi)DNP與緩發(fā)中子軸向分布Fig.4 Axial distribution of DNP and delayed neutron in the hottest channel at rated condition

2.2 入口流量對βeff的影響

液態(tài)燃料熔鹽堆中燃料的流動導致引起DNP的空間再分布，部分DNP流出堆芯并衰變，從而對βeff產(chǎn)生影響。圖5示出采用解析方法和數(shù)值方法計算得到的不同額定流量百分比下堆芯βeff的變化?？梢钥吹?，兩種方法計算得到的βeff變化趨勢基本一致，即入口流量增加，βeff減小。這是由于流量增加，更多的DNP流出堆芯且在堆外衰變，導致了反應(yīng)性損失，從而引起了βeff的減小。除0%額定流量外，解析方法得到的計算結(jié)果始終高于數(shù)值程序計算結(jié)果，二者最大相對誤差σmax為6.26%，額定流量下相對誤差σnormal為4.18%。

圖5 不同流量下的有效緩發(fā)中子份額Fig.5 Effective delayed neutron fraction at different percentages of flow rate

2.3 堆外流動時間對βeff的影響

對于半衰期較長的DNP，其隨燃料流出堆芯后可能未發(fā)生衰變重新進入堆芯，從而使堆芯反應(yīng)性及βeff產(chǎn)生波動，這一波動主要與燃料在堆外回路的流動時間有關(guān)[17]。圖6示出采用解析方法和數(shù)值方法計算得到的堆芯βeff隨燃料在堆外流動時間的變化，此過程維持入口流量為額定值?？梢钥吹剑耬ff隨燃料在堆外的流動時間的增加而減小，80 s后趨于穩(wěn)定。這是因為隨著燃料在堆外流動時間的增加，更多DNP將在堆外衰變，返回堆芯的DNP減少，從而使βeff減小。由表2可知，第1組DNP的平均壽命最大，約為80.65 s，因此當堆外流動時間超過80.65 s時，所有DNP均將在堆外衰變，因而無任何DNP返回堆芯，βeff降低至最小值，趨于穩(wěn)定?？梢钥吹?，60 s前，兩種方法計算結(jié)果存在一定偏差，解析方法仍高于數(shù)值方法計算結(jié)果，二者最大相對誤差為10.61%，額定工況下相對誤差為4.33%。結(jié)合2.2節(jié)討論可發(fā)現(xiàn)，解析方法計算得到的βeff相對于數(shù)值方法偏大。這是由于，解析方法忽略了堆芯的非均勻性和反射層的作用，采用均勻的流速場，認為堆內(nèi)所有緩發(fā)中子具有相同的價值，即為1，而堆外所有緩發(fā)中子價值為0?？偟膩碚f，由于燃料的流動性，相對于先驅(qū)核衰變產(chǎn)生緩發(fā)中子的位置，先驅(qū)核產(chǎn)生的位置中子價值是比較高的[3]，這也就是說，解析方法總體上高估了堆內(nèi)緩發(fā)中子的價值，導致βeff計算結(jié)果高于數(shù)值方法的計算結(jié)果。

圖6 有效緩發(fā)中子份額隨燃料在堆外回路流動時間的變化Fig.6 Effective delayed neutron fraction change with fuel residence time out of core

表2 緩發(fā)中子先驅(qū)核份額及其衰變常量Table 2 Delayed neutron precursor fraction and precursor decay constant

3 結(jié)論

有效緩發(fā)中子份額βeff的精確計算可為反應(yīng)堆設(shè)計和安全分析提供合理準則和依據(jù)。燃料的流動性、特殊的傳熱方式和釷鈾燃料循環(huán)的使用導致液態(tài)燃料熔鹽堆βeff的計算方法與固態(tài)燃料反應(yīng)堆不同。本文針對石墨慢化通道式熔鹽堆，分別基于解析方法和數(shù)值方法提出了計算βeff的數(shù)學模型，計算了MSRE在額定工況下的DNP損失份額和堆內(nèi)DNP濃度分布，并分析了燃料在堆外流動時間和入口流量對βeff的影響，可得到以下結(jié)論。

1) 兩種方法均可對DNP行為提供合理描述。解析方法忽略了堆芯的非均勻性和反射層的作用，總體上高估了堆內(nèi)緩發(fā)中子的價值，導致βeff計算結(jié)果相對于數(shù)值方法偏大。

2) 固定燃料在堆外流動時間，βeff隨入口流量的增加而減小。

3) 固定入口流量，βeff隨燃料在堆外流動時間的增加而減小，80 s后趨于穩(wěn)定。