液體燃料熔鹽堆三維穩(wěn)態(tài)分析程序開發(fā)

周建軍，張大林，秋穗正，蘇光輝，田文喜，巫英偉

(西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

1 物理模型

圖1 MOSART堆芯結(jié)構(gòu)

MOSART堆是一種先進的概念設(shè)計熔鹽堆，它在放射性廢物的有效焚化以及降解長壽命廢料毒性和燃料閉式循環(huán)等方面均有很大的優(yōu)勢，最早由俄羅斯RRC-KI研究中心的Ignatiev等[7]提出。圖1示出MOSART堆芯結(jié)構(gòu)的概念設(shè)計。為了計算模擬的方便，本文將堆芯結(jié)構(gòu)進行簡化，僅考慮堆芯分流板以上到堆芯出口的部分，簡化后的堆芯模型如圖2所示，堆芯直徑為1.9 m，在壁面附近有0.2 m厚的石墨反射層，采用錐形出口，燃料鹽從下部進入，從上部流出，MOSART堆的主要熔鹽物性和堆芯參數(shù)列于表1。

圖2 堆芯結(jié)構(gòu)

表1 MOSART堆芯初始設(shè)計參數(shù)

2 數(shù)學(xué)模型

2.1 中子物理模型

熔鹽堆是一種液體反應(yīng)堆，燃料鹽溶解于氟化鹽中在一回路循環(huán)，燃料處于流動狀態(tài)。這與固體反應(yīng)堆相比有很大差距，在進行物理計算時須考慮流動的影響。因此，為得到適用于液體堆的物理控制方程，本文從最原始的粒子守恒來進行推導(dǎo)，最終得到考慮了流動影響的中子擴散方程：

(1)

對于緩發(fā)中子先驅(qū)核同樣可用這種方法得到，緩發(fā)中子先驅(qū)核守恒方程可表示為：

(2)

由于在石墨反射層中無流動特性和增殖特性。因此，可得到石墨反射層內(nèi)的中子擴散方程和緩發(fā)中子方程：

(3)

Ci=0

(4)

式中：χp為瞬發(fā)中子能譜；χd為緩發(fā)中子能譜；φ為中子通量密度，cm-2·s-1；λ為緩發(fā)中子衰變常量，s-1；C為緩發(fā)中子先驅(qū)核濃度，cm-3；D為擴散系數(shù)，cm；Σf為宏觀裂變截面，cm-1；Σs為宏觀散射截面，cm-1；v為中子速度。其中，方程中的6群緩發(fā)中子份額和緩發(fā)中子衰變常量列于表2。

表2 緩發(fā)中子份額和衰變常量

2.2 流動和換熱模型

穩(wěn)態(tài)情況下熔鹽在堆芯內(nèi)的平均流速為0.5 m/s，文獻[8]分析表明，穩(wěn)態(tài)情況下堆芯內(nèi)的Re將超過100 000，堆芯內(nèi)的流動為湍流流動，本文選用標準k-ε湍流模型來模擬，控制方程如下。

質(zhì)量守恒方程：

(5)

動量守恒方程：

μ

(6)

式中：uj為速度張量形式；Sui為廣義源項的張量形式；μ為流體的動力黏性系數(shù)。

能量守恒方程：

(7)

式中，ST為源項，本文是堆芯的平均功率，包括快中子、慢中子及緩發(fā)中子產(chǎn)生的熱量。

湍流控制方程：

(8)

式中：φ為廣義變量，表示k或ε；Г為廣義擴散系數(shù)；S為源項。對于k方程和ε方程，廣義擴散系數(shù)可分別表示為：

Γk=η+ηt/σk

(9)

Γε=η+ηt/σε

(10)

式中：η為黏性系數(shù)；ηt為湍流黏性系數(shù)；湍流模型系數(shù)σk=1.0，σε=1.3。初始k可取初始動能的5%來計算，而ε可通過下式計算：

(11)

式中：Cμ為湍流模型中的經(jīng)驗系數(shù)，一般取0.09；l為湍流長度尺度。

3 數(shù)值計算方法

中子擴散方程、緩發(fā)中子先驅(qū)核方程和流動傳熱方程均可化為式(12)的模型方程，模型統(tǒng)一后各變量所對應(yīng)的參數(shù)列于表3。

div(ρUφ)=div(Γgradφ)+S

(12)

需編制程序?qū)Ψ匠踢M行數(shù)值求解。首先對方程進行離散，為保證方程的守恒性，采用有限容積積分法[9]對控制方程進行離散。圖3為三維坐標網(wǎng)格系統(tǒng)下的控制容積P，對其積分可得到控制方程的離散方程。

有限容積法的三維計算控制容積如圖3所示，對控制方程可采用控制容積積分法進行離散，三維坐標下的控制方程離散格式可寫為：

表3 物理模型中的廣義變量

圖3 控制容積

aPφP=aEφE+aWφW+aNφN+

aSφS+aBφB+aTφT+b

(13)

其中：

aP=aE+aW+aN+aS+aB+aT-SPΔV,

b=ScΔV

式中：下標W、S、E、N、B、T為節(jié)點；w、s、e、n、b、t為界面。

在離散方程中，界面上的當(dāng)量擴散系數(shù)De、Dw、Dn、Ds、Db、Dt按調(diào)和平均計算。得到的方程組為5對角方程組，采用TDMA算法并輔以塊修正方法求解。耦合計算流程示于圖4。

圖4 耦合計算流程

4 程序驗證及計算結(jié)果分析

4.1 程序驗證

由于本文開發(fā)的是一三維物理熱工耦合程序，而這種應(yīng)用于熔鹽堆的成熟的三維耦合程序和基準題尚很少見，但本文的程序既可考慮緩發(fā)和流動也可不考慮緩發(fā)和流動，因此對于程序的驗證可從中子物理模型驗證和流動傳熱模型驗證分別進行。中子物理模型可用常見的三維基準題——Small LWR Core來驗證[10]，計算結(jié)果列于表4，而對于流動傳熱模型可利用研究較多的三維自然對流模型來驗證，計算值與實驗值[11]的對比結(jié)果如圖5所示。

表4 中子物理基準題驗證結(jié)果

圖5 溫度隨無量綱長度的分布

通過與三維基準題對比可發(fā)現(xiàn)，中子物理模型的計算結(jié)果與基準題的計算結(jié)果符合較好，證明了本文中采用的中子物理型的準確性。另外，通過計算三維自然對流并與文獻中的實驗結(jié)果進行對比發(fā)現(xiàn)，除個別點有一定的誤差外，本文中的計算結(jié)果與實驗值符合較好，這也證明了本文中采用的流動傳熱模型的準確性，從而也驗證了開發(fā)的三維耦合程序的準確性。

4.2 計算結(jié)果分析

以快譜熔鹽堆MOSART作為研究對象，分析了在穩(wěn)態(tài)情況下入口速度效應(yīng)對堆芯物理熱工特性的影響，得到入口速度分別為0.1、0.4、0.5 m/s時的快中子分布、熱中子分布、緩發(fā)中子分布及堆芯內(nèi)的溫度在高度方向的中截面沿徑向分布。

1) 堆芯快中子

不同入口流速下快中子在堆芯高度方向中截面處沿徑向分布曲線如圖6所示。

圖6 快中子分布

從圖6可看出，作為一個均勻液體堆，MOSART堆的快中子通量密度在堆芯內(nèi)的分布比較對稱，最大值處于堆芯中心位置，流速對于快中子的分布影響不大。

2) 堆芯熱中子

在不同入口速度下沿堆芯高度方向中截面處的熱中子徑向分布曲線如圖7所示。由圖7可見，熱中子的分布在靠近堆芯壁面處有一增大的區(qū)域，這主要是因石墨反射層會使熱中子裂變截面升高而造成的1個突變。通過對比快中子和慢中子的通量密度分布可看到，在本文研究的對象MOSART堆中快中子通量密度要較熱中子通量密度高約1個數(shù)量級，因此，說明MOSART是1個快譜反應(yīng)堆。且入口流速對于熱中子在堆芯內(nèi)部的分布影響也很小。

圖7 熱中子分布

3) 堆芯緩發(fā)中子先驅(qū)核

對于液體燃料熔鹽反應(yīng)堆，流動對緩發(fā)中子先驅(qū)核的分布有很大影響，因此研究不同流速下緩發(fā)中子先驅(qū)核的堆芯分布對分析堆芯物理特性是非常有意義的，在兩種不同流速下6組緩發(fā)中子先驅(qū)核沿徑向分布曲線如圖8所示。

6組緩發(fā)中子先驅(qū)核分布表明，燃料鹽的流動作用，使緩發(fā)中子先驅(qū)核沿流動方向整體向下游移動，峰值不像固體燃料反應(yīng)堆中的會出現(xiàn)在堆芯的中心。緩發(fā)中子先驅(qū)核向下游移動的程度反映了流動對其影響程度，比較6組緩發(fā)中子先驅(qū)核濃度的分布可看出，緩發(fā)中子先驅(qū)核的衰變常量越小，流動對其影響越大。

圖8 緩發(fā)中子先驅(qū)核分布

4) 堆芯溫度

當(dāng)入口速度為0.1、0.4和0.5 m/s時，不同流速下堆芯內(nèi)高度方向中截面沿軸向溫度分布曲線如圖9所示。

圖9 溫度分布

從圖9可看出，流速對于堆芯內(nèi)的溫度分布影響較大，特別是在堆芯中間位置這種影響更明顯，當(dāng)流速較小時，堆芯中心線處的溫度較其他區(qū)域的溫度要高，但當(dāng)流速較大時整個堆芯除了壁面附近的區(qū)域之外的其他區(qū)域溫度分布相對比較均勻。

5 結(jié)論

本文以MOSART堆型作為研究對象，建立了適用于液體燃料熔鹽反應(yīng)堆的中子擴散方程，并結(jié)合流動傳熱方程，采用有限容積法對方程進行離散，編制了符合液體反應(yīng)堆的三維耦合穩(wěn)態(tài)分析程序，通過與中子物理和流動換熱基準題的對比，驗證了本程序的正確性。最后，利用本程序分析了MOSART堆芯的穩(wěn)態(tài)特性，通過分析入口流速對堆芯物理熱工特性的影響，可看出入口流速的變化對熱中子和快中子的分布影響較小，而對于緩發(fā)中子的分布和堆芯內(nèi)溫度的分布影響較大。

參考文獻：

[2] 張大林,秋穗正,劉長亮,等. 新概念熔鹽堆堆芯穩(wěn)態(tài)熱工水力計算[J]. 工程熱物理學(xué)報,2008,29(6):979-982.

ZHANG Dalin, QIU Suizheng, LIU Changliang, et al. Steady thermal-hydranlic analysis for the core of a new concept molten salt reactor[J]. Journal of Engineering Thermophysics, 2008, 29(6): 979-982(in Chinese).

[3] 張大林,秋穗正,劉長亮,等. 新概念熔鹽堆物理計算方法研究及程序設(shè)計[J]. 原子能科學(xué)技術(shù),2008,42(12):1 104-1 108.

ZHANG Dalin, QIU Suizheng, LIU Changliang, et al. Nuclear calculation and program development for moltensalt reactor[J]. At Energy Sci Technol, 2008, 42(12): 1 104-1 108(in Chinese).

[4] ZHANG D L, QIU S Z, SU G H, et al. Development of a steady state analysis code for a molten salt reactor[J]. Ann Nucl Energy, 2009, 36(5): 590-603.

[5] ZHANG D L, QIU S Z, SU G H, et al. Analysis on the neutron kinetics for a molten salt reactor[J]. Prog Nucl Energy, 2009, 51(4-5): 624-636.

[6] YAMAJI B, CSOM G, ASZDI A. Three-dimensional numerical investigation of a molten salt reactor concept with the code CFX-5.5[C]∥Nuclear Energy for New Europe 2002. Kranjska Gora: [s.n.], 2002.

[7] IGNATIEV V, FEYNBERG O, GNIDOI I, et al. Progress in development of Li, Be, Na/F molten salt actinide recycler & transmuter concept[C]∥Proceedings of ICAPP 2007. Nice, France: [s.n.], 2007.

[8] YAMAMOTO T, MITACHI K, SUZUKI T. Steady state analysis of small molten salt reactor: Effect of fuel salt flow on reactor characteristics[J]. JSME International Journal B: Fluids and Thermal Engineering, 2005, 48(3): 610-617.

[9] 陶文銓. 數(shù)值傳熱學(xué)[M]. 西安:西安交通大學(xué)出版社,2001.

[10] TAKEDA T, IKEDA H. 3D neutron transport benchmarks[J]. J Nucl Sci Technol, 1991, 28(7): 656-669.