小型鉛冷快堆堆芯物理計算軟件的開發與臨界實驗驗證

陳仁宗,周 琦,朱慶福,夏兆東,寧 通,馬驍笛,孫 旭

(1.清華大學 能源環境經濟研究所,北京 100084;2.中國原子能科學研究院,北京 102413)

鉛冷快堆(LFR)是第4代核能系統論壇(GIF)選定的6種優選堆型之一[1],隨著材料和設備技術的進步,LFR在中子物理、傳熱、安全等方面具有的獨特優勢不斷凸顯[2],已成為世界各國新型核能系統研發的重點[3]。當前LFR堆芯物理計算主要以確保精度的蒙特卡羅方法為主,進行全堆芯的建模計算,隨著型號設計的不斷深入,效率較低的蒙特卡羅全堆計算設計方法已無法滿足LFR設計需要巨大計算需求。

如果參考壓水堆堆芯物理計算的兩步法[4],首先采用較為精確的中子輸運理論給出柵元和組件的等效均勻化參數,再將其應用于基于中子擴散理論的堆芯計算,相較于目前常用的全堆芯蒙特卡羅精確建模計算,能夠大幅提高計算效率。與壓水堆或鈉冷快堆(SFR)相比,小型LFR冷卻劑中的鉛、鉍等核素對中子的彈性、非彈性散射更強,輸運截面較大,但慢化能力更小,燃料元件柵距變化范圍大,空間異性強,燃料組件的等效均勻化遇到的問題更加特殊[5]。

本文參考國內外研究成果,針對小型LFR特有的材料與能譜環境,利用開源的蒙特卡羅計算軟件OpenMC計算多群截面,能夠更好地描述柵元與組件在空間與能量上的局部特性。調用開源的有限體積法(FVM)軟件OpenFOAM的各種求解器,開發中子擴散求解器DESOF,通過Python實現蒙特卡羅均勻化與擴散計算的耦合,形成完整的堆芯物理計算軟件MCDESOF,并開展臨界實驗,利用實驗數據驗證MCDESOF的適用性。

1 蒙特卡羅方法生成多群截面

1.1 多群截面的蒙特卡羅統計方法

蒙特卡羅均勻化的核心內容是將求解輸運問題得到的各種事件統計結果轉化為均勻化常數,這些事件主要包括散射、裂變、俘獲吸收、(n,xn)等反應,均勻化常數包括各種反應的群截面、輸運修正、散射矩陣、各階勒讓德項,次級中子產生矩陣、裂變產額、裂變能譜,以及動力學計算需要的中子速度倒數、緩發中子先驅核產額與衰變常量等。為了簡化表達式,引入內積算符〈·,·〉來代表能量、空間或角度的積分形式,這樣第x類核反應的反應率計數估計可表示為:

(1)

式中:Σx為第x類核反應的宏觀截面;φ為中子通量密度;r為中子的空間位置;E為中子能量;Ω為中子的運動方向;V為r的積分空間范圍;S為Ω的積分空間角范圍。

(2)

使用徑跡長度法,通過設置反應率和中子通量密度的計數箱可計算得到式(2)的分子與分母,分別是k空間(柵元)、第x類核反應、第g能群(能量區間)的核反應率R和k空間、第g能群的中子通量密度φ。

1) 輸運修正

式(2)考慮的散射是各向同性的,但實際散射存在各向異性,因此需要引入輸運修正處理各向異性散射。首先,定義以勒讓德多項式Pl(μ)展開的散射核為:

(3)

式中:Σsl為散射核;E′為散射前的能量;μ為散射角余弦。

對于k空間、Vk體積、g群能量范圍勒讓德多項式可定義為:

φ(r,E′,Ω)dE′dEdΩdr

(4)

式中,g′為散射前的能群。

(5)

式中,G為能群總數。

(6)

由于輸運修正項必須通過模擬估計法計算得到,因此總碰撞反應率與中子通量密度也需要通過模擬估計法得到,確保計算結果準確性。

2) 散射矩陣

(7)

對于群間散射,不僅要考慮群間散射的概率,還要考慮中子散射前后的能量和角度分布,蒙特卡羅方法可通過設置能量與角度計數箱,統計群間散射率與散射后的角分布。參考Redmond[6]提出的計數箱設置方法,首先模擬用于統計群間散射的散射事件,根據入射粒子的能量通過散射事件計算出射粒子的能量、散射角余弦,其次入射及出射能量決定了哪個群間散射率計數箱儲存信息,散射角余弦則決定了哪個群間散射角余弦接受信息。

類似輸運修正(式(6)),散射矩陣的修正可表示為:

(8)

3) 需要統計的物理量

為了提供后續擴散計算所需要的擴散系數、移出截面、裂變截面等完整的輸入參數,參考Boyd等[7]提出的計算思路實現。

1.2 利用OpenMC生成多群截面

利用OpenMC生成多群截面計算步驟如下。

1) 問題的預處理,主要目標是建立蒙特卡羅計算模型,根據問題所包含的幾何、材料與均勻化計算的要求,根據所選蒙特卡羅軟件的要求,利用Python按照特定的語法建立計算模型,準確定義幾何、材料、截面、計算參數與計數箱的設置。

2) 蒙特卡羅計算,主要目標是生成多群截面所需的原始數據。調用蒙特卡羅軟件進行計算后,將生成可直接讀取的ASCII格式與便于處理的二進制格式的包含計算結果的文件,繼續利用Python將其數據提取、處理并計算得到吸收、裂變、散射截面與散射矩陣、不同階的勒讓德散射項、裂變能譜與緩發中子參數等基礎數據。

3) 多群截面的應用,主要目標是對基礎多群數據進行處理,具體過程將根據擴散計算采用的近似處理方法進行處理,如進行輸運修正,生成擴散系數、散射移出截面、散射移入截面等。具體的文件格式根據擴散求解器的構造利用Python進行適應性的調整。

1.3 計算流程

以OpenMC為例,本文利用蒙特卡羅方法計算多群截面并為多群中子擴散計算提供均勻化參數的計算流程,如圖1所示。

圖1 蒙特卡羅方法為中子擴散計算提供均勻化參數的流程

2 多群中子擴散軟件開發

2.1 有限體積法與OpenFOAM

有限體積法又稱為控制體積法,其基本思路與有限元方法類似,也是節點劃分、方程離散、方程求解散步,不同在于將計算區域劃分為一系列不重復的控制體積,并使每個網格點周圍有1個控制體積,將待解的微分方程對每一個控制體積積分,得出1組離散方程,其物理意義就是因變量在有限大小的控制體積中的守恒原理[8]。

OpenFOAM[9]是一個完全面向對象的CFD類庫,包含許多可執行文件,也稱為應用程序包,從文件組織結構來說這些應用程序大體可以分為兩大類:求解器和輔助工具。求解器用來求解連續介質力學中的某個特定問題,而輔助工具主要用來進行數據操作、輔助求解器完成計算任務。

2.2 DESOF開發

中子擴散方程可看作一個基礎的偏微分方程,OpenFOAM針對中子擴散方程的求解通過式(9)的結構進行,主要包含拉普拉斯項和源項兩類離散項,左邊第1項泄漏項與第2項吸收項(實際上包含吸收與群件散射移出兩部分)采用自帶的隱性finiteVolumeMethod離散方式實現右邊的源項neutronGroups包含群間散射移入項與裂變產生項。

-fvm∷Laplacian(D,φ)+fvm∷Sp(ΣA,φ)=

neutronGroups.Sp(neutron)

(9)

基于C++語言,根據OpenFOAM求解穩態中子擴散方程采用的代碼形式,本文開發了基于OpenFOAM的擴散方程求解器DESOF。DESOF的基本計算流程如圖2所示,主要分為3部分:OpenFOAM輸入、網格生成及擴散方程的求解。

2.3 數值驗證

選取二維IAEA基準題對DESOF進行驗證?；鶞暑}是一簡化的兩群PWR基準題[10],它由241個組件、4種不同的材料構成。堆芯按雙區布料方案布置,組件幾何尺寸為20 cm×20 cm,堆芯外邊界條件為真空邊界?；鶞暑}1/4堆芯幾何布置、DESOF構建的幾何布置如圖3所示。圖3中,Jin為入射界面流。

圖3 基準題(a)與DESOF(b)建立的幾何模型

基準題keff的參考值為1.029 585,DESOF的計算值為1.029 400,二者偏差為18.5 pcm。將組件歸一化功率分布計算值與參考值進行比較,如圖4所示。由圖4可見,最大相對偏差出現在邊緣處的組件,為1.00%。從對比結果看,只要輸入準確的均勻化參數,DESOF的多群擴散計算準確性能夠得到保證。

圖4 基準題的組件歸一化功率分布

3 組件均勻化與驗證

3.1 MCDESOF的集成與功能性檢驗

為了提高計算效率,本文基于Python開發了鏈接蒙特卡羅均勻化計算與多群擴散方程求解器的Pylink程序,實現便捷的蒙特卡羅均勻化建模、計算與數據處理,傳遞并進行擴散求解,形成了完整的MCDESOF軟件。其中,蒙特卡羅均勻化建模模塊將自動生成OpenMC運行所需的幾何、材料、設置、圖形和計數等部分的XML格式的輸入文件,調用OpenMC完成蒙特卡羅輸運計算;多群截面處理模塊處理蒙特卡羅輸運計算結果,將其處理運算得到多群擴散需要的均勻化參數;多群擴散計算模塊對指定網格的幾何排列方式逐一進行賦值、設定網格的邊界條件以及圖形化顯示,并調用DESOF完成擴散求解;結果處理模塊將擴算計算結果轉化為keff和堆芯功率分布或組件分布,便于后續處理。

3.2 小型LFR燃料組件模型

參考各種LFR燃料組件的特點,本文構建了4種類型的組件,如圖5所示。選取金屬鈾芯塊進行計算驗證,燃料富集度為15%,燃料包殼為不銹鋼,冷卻劑材質選擇鉛鉍合金。燃料棒芯塊外徑為0.9 cm,包殼內徑為0.92 cm,包殼外徑為1 cm,燃料棒對邊距為1.40 cm,燃料棒慢化劑外邊界條件為全反射邊界。

使用蒙特卡羅方法對4個計算案例進行輸運計算,給出k∞與各群相對中子通量密度計算結果,作為本文提出的等效均勻化方法的對比參考值。蒙特卡羅計算考慮keff統計誤差1σ小于10 pcm。采用3群能群結構(分界能分別為5.53×10-3MeV和0.8 MeV)與網距為1/2柵距的網格劃分方式,利用MCDESOF進行了計算。蒙特卡羅均勻化計算考慮的k∞統計誤差小于10 pcm,擴散計算以特征值偏差|ε|<1 pcm作為收斂判據。將MCDESOF的計算值與蒙特卡羅參考值對比分析,結果列于表1、2。

表1 MCDESOF計算值與參考值的k∞對比

表2 案例1 MCDESOF計算值與參考值的3群中子通量密度對比

由表1、2可見,k∞最大偏差超過400 pcm,最小偏差也超過了200 pcm,各群的歸一化中子通量密度最大相對偏差超過1%,符合得不夠理想,因此有必要對蒙特卡羅多群計算得到的均勻化參數的等效性進行核實,并進行適當修正。

3.3 超級均勻化方法的應用

為保證精度,均勻化過程必須保證一些重要的、反映堆芯物理特性的參量在均勻化過程中守恒,這些參量一般包括反應率、界面流和特征值。在實際情況中,全部滿足3個參量守恒是不現實的,一般通過放寬約束條件或增加均勻化參數的自由度來達到或大致滿足等效均勻化的要求,當前應用廣泛的是廣義均勻化理論(GET)和超級均勻化(SPH)方法。廣義均勻化理論通過組件邊界處的不連續因子保證界面流守恒,由Henry等[11]提出。Gunow等[12]提出不連續因子的具體實現方法,該方法目前在輕水堆工程計算中應用廣泛。Koebke[13]提出超級均勻化方法直接調整截面保持反應率守恒。Hébert等[14]先后發展了三代超級均勻化方法。

相比廣義均勻化理論,超級均勻化方法的優勢在于無需額外的均勻化參數,這使得其在精細柵元計算中成為首選。本文使用超級均勻化方法,將中子擴散方程的右邊簡化為固定源項,在空間k、能群g中引入SPH因子μk,g,擴散方程通過修正移出項實現平衡,如式(10)所示。

(10)

式中:φg為g群中子通量密度;Σt,k,g為空間k處的g群總截面;Qk,g為空間k處的g群源項。

本文利用迭代算法得到SPH因子,首先將式(10)轉化為迭代的形式:

(11)

(12)

針對案例1,采用超級均勻化方法對各群截面進行修正,在式(16)中設定ε1為50 pcm,ε2為0.5%。

用超級均勻化方法對4個案例的各能群均勻化參數進行修正,得到修正后各組件的k∞計算結果,并與蒙特卡羅參考值進行對比,結果列于表3。由表3可見,經過修正后k∞的偏差基本下降了1個數量級,最大偏差不超過50 pcm,充分顯現出該方法對于確保均勻化等效性的作用。

表3 修正后MCDESOF計算值與參考值的k∞對比

4 臨界實驗驗證

4.1 臨界裝置與實驗方案

為了驗證MCDESOF對于小型LFR堆芯物理計算的適用性,在中國原子能科學研究院鈾棒柵臨界實驗裝置的鉛堆上開展了臨界實驗。該裝置是為開展先進反應堆堆芯物理研究而構建的鉛鉍快堆零功率裝置,陸續開展了鉛基ADS次臨界反應堆物理特性實驗研究(啟明星Ⅱ號)[15]和鉛鉍快堆基準性零功率實驗研究[16],已成為獲取小型LFR核心堆芯物理實驗參數的重要實驗平臺[17]。

本次臨界實驗方案的堆芯中,以金屬鈾為核燃料,主要以金屬鈹為內部反射層,以石墨為外部反射層,以鉛鉍合金作為堆芯介質體,燃料棒、鈹棒與控制棒等布置于打孔的鉛鉍介質體中。燃料元件按六角形柵格排布,4圈燃料構成1個組件,單個組件包含37根燃料棒,整個堆芯包含19個燃料組件。安全棒采用富集度為92%的碳化硼作為吸收體材料,布置在外圍反射層組件中。臨界實驗堆芯布置如圖6所示。

圖6 鉛堆臨界實驗堆芯布置

4.2 MCDESOF建模

1) 組件網格劃分

在徑向上,采用菱形網格劃分方式,針對每個不同的組件,通過預先制定好的蒙特卡羅均勻化截面,將其填充到六角形的組件中。對于六角形組件的網格劃分,由于基于OpenFOAM的MCDESOF本身具有完整的六面體網格劃分方式,燃料組件和鈹反射層組件的劃分方式如圖7所示。

圖7 燃料組件(a)與粗鈹棒組件(b)的網格劃分

在軸向上,活性段劃分100層,上下反射層各劃分了20層水平方向的平行網格。

2) 蒙特卡羅均勻化與修正方案

根據堆芯中的各種組件的排列方式以及堆芯的對稱性,組件可分為6種,如圖8所示。其中:1號和2號組件,它們擁有較多數量的燃料元件,采用反射邊界的組件計算模型,考慮SPH因子,為擴散計算提供修正后的多群截面;3號、4號組件,它們不含燃料元件,計算時將2號組件的1/2與3號、4號組件以及外圍的石墨反射層組成“超組件”計算模型,內部采用反射邊界,外部考慮泄漏邊界,考慮SPH因子,為擴散計算提供修正后的多群截面;5、6號組件,控制棒不插入時,這兩種組件等價于4號組件,不需要額外計算,一旦控制棒插入,采用的相同的方法,引入2號組件的1/2組成“超組件”計算模型,考慮SPH因子,為擴散計算提供修正后的多群截面。

圖8 堆芯蒙特卡羅均勻化的區域設置

3) 堆芯擴散計算網格劃分

對于堆芯擴散計算,考慮到對稱性,采用菱形網格讀取蒙特卡羅均勻化計算得到的堆芯各組件多群截面,進行全堆擴散計算。對于堆芯外圍與石墨反射層接觸的不規則結構的網格,如圖9所示,以任意一處邊界為例,首先計算出組件與邊界處的交點1、3、5的坐標,然后使用三棱柱prism關鍵字對(1 2 3)以及(3 4 5)兩個面所在的棱柱進行定義。由于1-3以及3-5的邊分別為圓弧形,因此在進行網格定義的輸入時,除點、面的定義外,還需增加對邊的定義。

圖9 堆芯邊界處不規則網格的處理方式

4.3 臨界實驗結果與計算值對比

在啟明星Ⅳ號上開展了臨界裝載實驗和安全棒價值測量實驗。臨界裝載實驗堆芯共裝載198根燃料元件,1個空孔道,其余位置被金屬鈹反射元件填充,此時反應堆倍增周期為82.9 s,反應性為47.1 pcm。安全棒價值測量實驗采用了基于逆動態方法的落棒法[18]。落棒后使用逆動態方法得到1號安全棒價值為1 383.8 pcm。

利用MCDESOF對臨界附近的反應性與安全棒價值進行了計算,同時利用OpenMC采用蒙特卡羅方法直接對全堆進行建模計算,核截面數據庫選擇ENDF/B-Ⅶ.0。計算的硬件平臺為Dell T7920服務器,64核心并行計算,蒙特卡羅均勻化統計誤差1σ為10 pcm。實驗測量值、蒙特卡羅全堆建模計算值、MCDESOF計算值三者的對比列于表4。

表4 實驗測量值、蒙特卡羅全堆建模計算值與MCDESOF計算值的對比

由表4可見,MCDESOF計算值、蒙特卡羅全堆建模計算值與實驗測量值間的偏差相當,最大偏差不超過200 pcm。但是兩種理論方法在安全棒、調節棒插入前的計算值大于實驗值,而在安全棒或調節棒插入后計算值小于實驗值,這是因為安全棒及其配套的結構比較復雜,需要提高建模的精度,特別是吸收體密度、填充范圍與尺寸需要進一步核實,減少建模誤差。

從計算時間上看,OpenMC采用蒙特卡羅方法直接對全堆進行建模計算需要約2 h完成,MCDESOF內部實施的組件蒙特卡羅均勻化計算,截面處理與擴散計算總共不超過25 min,計算時間小于蒙特卡羅方法直接計算的25%,計算效率有顯著的提升。

5 結論

本文開展了小型LFR堆芯物理分析軟件的開發與臨界實驗驗證。利用開源蒙特卡羅軟件OpenMC,實現柵元與組件的蒙特卡羅均勻化?；陂_源有限體積法軟件OpenFOAM開發了中子擴散求解器DESOF,通過Python實現蒙特卡羅均勻化與擴散計算的耦合,形成完整的堆芯物理計算軟件MCDESOF,利用超級均勻化方法實現了組件的等效均勻化,通過4種LFR典型的燃料組件模型對MCDESOF進行了數值驗證。在鈾棒柵臨界實驗裝置鉛堆堆芯上開展了臨界裝載實驗與控制棒價值測量實驗,利用MCDESOF進行了針對性的建模,進行了全堆擴散計算,計算結果與測量結果進行對比,臨界附近的反應性偏差小于100 pcm,安全棒價值相對偏差小于200 pcm,計算準確度達到蒙特卡羅全堆建模計算的水平,所需的計算時間則遠小于蒙特卡羅全堆建模所需時間。

下一步將完善控制棒均勻化計算方法,引入誤差傳遞公式,定量分析蒙特卡羅均勻化截面的不確定度對擴散計算的影響,以及開展三維時空相關的瞬態計算功能開發,支持小型LFR瞬態安全分析。