陳仁宗,周 琦,朱慶福,夏兆東,寧 通,馬驍?shù)?孫 旭

(1.清華大學(xué) 能源環(huán)境經(jīng)濟研究所,北京 100084;2.中國原子能科學(xué)研究院,北京 102413)

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

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

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

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

1.1 多群截面的蒙特卡羅統(tǒng)計方法

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

(1)

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

(2)

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

1) 輸運修正

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

(3)

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

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

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

(4)

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

(5)

式中,G為能群總數(shù)。

(6)

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

2) 散射矩陣

(7)

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

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

(8)

3) 需要統(tǒng)計的物理量

為了提供后續(xù)擴散計算所需要的擴散系數(shù)、移出截面、裂變截面等完整的輸入?yún)?shù),參考Boyd等[7]提出的計算思路實現(xiàn)。

1.2 利用OpenMC生成多群截面

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

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

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

3) 多群截面的應(yīng)用,主要目標(biāo)是對基礎(chǔ)多群數(shù)據(jù)進(jìn)行處理,具體過程將根據(jù)擴散計算采用的近似處理方法進(jìn)行處理,如進(jìn)行輸運修正,生成擴散系數(shù)、散射移出截面、散射移入截面等。具體的文件格式根據(jù)擴散求解器的構(gòu)造利用Python進(jìn)行適應(yīng)性的調(diào)整。

1.3 計算流程

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

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

2 多群中子擴散軟件開發(fā)

2.1 有限體積法與OpenFOAM

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

OpenFOAM[9]是一個完全面向?qū)ο蟮腃FD類庫,包含許多可執(zhí)行文件,也稱為應(yīng)用程序包,從文件組織結(jié)構(gòu)來說這些應(yīng)用程序大體可以分為兩大類:求解器和輔助工具。求解器用來求解連續(xù)介質(zhì)力學(xué)中的某個特定問題,而輔助工具主要用來進(jìn)行數(shù)據(jù)操作、輔助求解器完成計算任務(wù)。

2.2 DESOF開發(fā)

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

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

neutronGroups.Sp(neutron)

(9)

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

2.3 數(shù)值驗證

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

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

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

圖4 基準(zhǔn)題的組件歸一化功率分布

3 組件均勻化與驗證

3.1 MCDESOF的集成與功能性檢驗

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

3.2 小型LFR燃料組件模型

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

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

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

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

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

3.3 超級均勻化方法的應(yīng)用

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

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

(10)

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

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

(11)

(12)

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

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

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

4 臨界實驗驗證

4.1 臨界裝置與實驗方案

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

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

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

4.2 MCDESOF建模

1) 組件網(wǎng)格劃分

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

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

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

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

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

圖8 堆芯蒙特卡羅均勻化的區(qū)域設(shè)置

3) 堆芯擴散計算網(wǎng)格劃分

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

圖9 堆芯邊界處不規(guī)則網(wǎng)格的處理方式

4.3 臨界實驗結(jié)果與計算值對比

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

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

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

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

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

5 結(jié)論

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

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