基于連續能量蒙特卡羅的快中子反應堆均勻化截面計算方法研究

郭 輝,沈宇陽,吳逸煒,陳萃嵐,宋去非,顧漢洋

(上海交通大學 核科學與工程學院,上海 200240)

快中子反應堆能有效提高鈾資源利用效率、減少高放廢物并有助于穩定核燃料循環,是第四代核能系統技術的重要組成部分,也是我國重點發展的先進堆型之一。截至目前,國內外共建造了近20臺各型快堆,積累了400余堆年的運行經驗。高效準確的反應堆物理分析方法對新型核能系統的研發至關重要。新設計的涌現與精度要求的提高,對快堆的物理計算方法提出了新的挑戰,如復雜結構的精確模擬、寬能譜和非均勻效應下的截面生成及高性能大規模并行計算等[1]。蒙特卡羅方法基于精細幾何和連續能量數據庫,具有較高的計算精度。但對于優化設計、燃料管理、瞬態計算及多物理耦合等任務,其對計算資源要求高的問題尤其顯著。快堆確定論兩步法通常由組件均勻化截面計算和堆芯擴散/輸運計算共同組成,具有較高的計算精度和效率,已廣泛應用于快堆工程設計與分析領域。但由于采用空間均勻化、能量分群和角度離散的數值方法,新型設計的復雜結構與能譜必須在均勻化截面生成中精確計算。

在中子能譜上,快堆中子能量主要位于keV～MeV區間,該能量區域共振現象復雜,非彈性散射與(n,xn)等閾能反應重要性增加。在空間上,快堆中子自由程長,堆芯空間耦合緊密。總體空間自屏效應弱,但吸收體等局部結構非均勻性強。快堆中子泄漏強且呈現明顯的各向異性。針對這些快堆物理計算的特點,國內外對快堆堆芯用均勻化截面開展了廣泛的研究。MC2-3、CONSYST、SARAX-TULIP、ERANOS-ECCO等零維、一維或準二維截面計算程序在傳統快堆中已得到廣泛驗證與應用。為提高燃料組件內非均勻結構與控制棒等復雜組件的計算精度,二維特征線(MOC)方法受到廣泛關注,如Guo等[2]在APOLLO3程序、Wei等[3]在SARAX程序中的研究。Faure等[4]基于APOLLO3進一步開發了三維特征線生成截面方法,試圖解決ASTRID-CFV等非均勻堆芯的均勻化截面生成問題,但其共振計算能力、計算效率和幾何處理能力有待提高。近年來,利用連續能量蒙特卡羅計算均勻化截面受到廣泛關注,Nikitin等[5]、Lin等[6]、Tran等[7]、Nguyen等[8]、Martin等[9]、Guo等[10]對蒙特卡羅生成快堆用截面開展了研究。蒙特卡羅方法可建立二維/三維組件與超組件乃至三維全堆模型生成均勻化截面,以充分考慮快堆長中子自由程引起的結構相互作用,精確處理局部復雜能譜變化。已有結果表明,蒙特卡羅產生均勻化截面與堆芯擴散計算結合具有良好的計算效率與精度,但需結合修正技術以提高控制棒等強吸收體的計算精度。Lin等[6]與Guo等[10]的研究表明,基于蒙特卡羅體積通量均勻化方法與堆芯輸運計算結合會高估堆芯反應性,需要發展適用于堆芯輸運計算的蒙特卡羅生成均勻化截面技術。

本文聚焦基于連續能量蒙特卡羅的快堆均勻化截面計算方法研究。首先,簡要梳理國內外蒙特卡羅生成快堆均勻化截面的研究進展;然后,介紹蒙特卡羅生成快堆均勻化截面的基礎理論模型,包括體積通量均勻化方法和針對堆芯擴散計算的超級均勻化等效修正方法(SPH),并提出針對堆芯輸運計算的蒙特卡羅通量矩均勻化方法(MHT);最后,基于1 000 MWth金屬燃料鈉冷快堆進行數值驗證與分析。

1 國內外研究現狀

通過蒙特卡羅方法生成的均勻化截面與多群堆芯計算結合的結果列于表1,其中主要包括Serpent/DYN3D[5]、Serpent/VARIANT[6]、MCS/RAST[7-8]、MCS/MCS(MG)[8]、Serpent/Griffin[9]、OpenMC/TRIVAC[11]、OpenMC/OpenMC(MG)[12-14]。驗證與應用算例的堆芯功率涵蓋65～3 600 MWth,冷卻劑涵蓋鈉與鉛鉍,燃料涵蓋UOX、MOX、金屬、碳化物燃料和氮化物燃料。根據表1結果,堆芯求解器類型可分為堆芯擴散求解器和堆芯輸運求解器。

表1 蒙特卡羅生成快堆均勻化截面對堆芯有效增殖因數的計算精度Table 1 Accuracy of Monte Carlo generated multigroup cross-section in predicting core reactivity

1.1 與堆芯擴散計算結合

蒙特卡羅體積通量均勻化方法與堆芯擴散計算結合的兩步法所得堆芯有效增殖因數與基準方法之間的偏差在-163～263 pcm之間。對于控制棒未插入堆芯,蒙特卡羅生成均勻化截面與堆芯擴散計算結合的兩步法表現出較好的計算精度與計算效率。然而,Nikitin等[5]、Nguyen等[8]、Guo等[11]的結果表明,直接使用該兩步法會普遍高估控制棒價值。使用SPH等效修正方法是減小控制棒價值高估的有效方法。

蒙特卡羅法生成均勻化截面在精細幾何構建上具有優勢。表1中均勻化截面生成模型主要為二維或三維的燃料組件模型、燃料組件與結構組件構成的超組件模型。Serpent/Griffin對MET-1000采用了三維全堆模型,并且對全堆芯截面進行了SPH修正,從而將與基準計算的反應性誤差降低到0 pcm。OpenMC/TRIVAC方法對CEFR堆芯采用了三維全堆模型和全堆芯SPH修正,從而將與實驗的反應性誤差降低到-105 pcm,該誤差主要為核數據庫引起的誤差[20]。全堆SPH等效修正方法在微型反應堆等強非均勻堆芯中亦可應用[21-22]。全堆SPH方法可為特定堆芯計算提供基準群參數,但其產生的截面在不同狀態下的魯棒性還有待進一步驗證。

1.2 與堆芯輸運計算結合

蒙特卡羅體積通量均勻化方法與堆芯輸運計算結合的兩步法會明顯高估堆芯有效增殖因數。Serpent/VARIANT兩步法高估二維OECD MET-1000堆芯反應性約643 pcm。OpenMC/OpenMC(MG)與MCS/MCS(MG)分別高估三維OECD MET-1000堆芯反應性約1 087 pcm和1 085 pcm。OpenMC使用體積通量均勻化方法與DRAGON5-TRIVAC中SP5簡化球諧函數輸運計算結合高估三維OECD MET-1000堆芯反應性約880 pcm。因此,需要發展針對堆芯輸運計算的均勻化方法提高其計算精度。蒙特卡羅通量矩均勻化方法生成截面可將OpenMC/OpenMC(MG)對三維OECD MET-1000堆芯反應性的高估從1 087 pcm降到374 pcm。

2 理論模型

2.1 蒙特卡羅體積通量均勻化方法

蒙特卡羅體積通量法均勻化截面計算方法已得到廣泛發展,如SERPENT[23]、MCNP6[24]、McCARD[25]、RMC[26]、JMCT[27]、OpenMC[28]、MCS[7-8]與MCX[29]等。本文以開源蒙特卡羅程序OpenMC[28]中的體積通量法均勻化截面計算方法為例,簡要介紹其均勻化截面計算方法[30],不同程序在一些細節上會略有區別。

1) 0階截面

蒙特卡羅體積通量法中0階宏觀截面定義如下:

(1)

〈Σx,i,φ〉k,g=

(2)

(3)

其中:x為反應類型,如總截面、吸收截面、裂變截面與有效裂變中子數等;i為核素;k為空間區域;g為能群;〈Σx,i,φ〉k,g和〈φ〉k,g為蒙特卡羅計數的反應率和標通量,OpenMC中采用徑跡長度估計。

2) 散射矩陣

OpenMC中可通過勒讓德多項式或直方分布的形式來表示散射矩陣。勒讓德多項式形式較為通用,其l階散射矩陣計算如下:

(4)

σs(r,E′→E,Ω′·Ω)ψ(r,E′,Ω′)

(5)

其中,〈Σsl,φ〉k,g′→g和〈φ〉k,g采用解析估計進行計數。總截面等采用徑跡長度估計,而高階散射反應率難以采用徑跡長度估計。為保持反應率的一致,OpenMC中可采用“一致性”散射矩陣:

Σsl,i,k,g′→g=Σs,i,k,gPsl,i,k,g′→g

(6)

其中:Σs,i,k,g為總散射截面,可由式(1)通過徑跡長度估計求得;Psl,i,k,g′→g為散射概率矩陣,可通過解析估計進行計算。

(7)

3) 輸運修正

輸運截面(Σtr)和輸運修正的散射矩陣(Σstr),簡稱TCP0截面,對減小擴散方程求解的誤差具有重要作用。OpenMC中采用限流近似方法進行輸運修正:

(8)

Σstr,i,k,g′→g=

(9)

其中,δg,g′為克羅內克函數。

擴散系數由下式計算:

(10)

嚴格意義上,擴散系數可用過累積徙動方法(CMM)[31-35]獲得。CMM可消除由于輸運修正近似引起的誤差,其對水堆的適應性已得到驗證[31-34]。

2.2 超級均勻化等效修正方法

超級均勻化等效修正方法(SPH)[36-37]已廣泛用于保持基準非均勻計算與均勻計算間的反應率守恒。在諸多文獻[7-8,11,38-39]中,SPH修正已應用于提高蒙特卡羅生成均勻化截面與堆芯擴散計算結合的計算精度。SPH等效修正法的主要流程如圖1所示。

圖1 超級均勻化等效修正因子計算流程[11]Fig.1 Iteration scheme of SPH[11]

首先,用基準非均勻模型計算均勻化截面,同時存儲非均勻體系下的中子通量。蒙特卡羅方法可基于連續能量和精細幾何進行組件、超組件及全堆芯模型的基準非均勻計算。

其次,基于以上所得均勻化截面,進行多群、均勻化堆芯計算,獲得均勻化通量。

再次,可獲得均勻化區域r、能群g的SPH等效修正因子:

(11)

(12)

然后,使用SPH等效修正因子修正截面,計算公式如下:

(13)

最后,將修正后的截面重新代入均勻化多群堆芯計算,獲得新的均勻化通量,產生新的SPH等效修正因子與修正截面,往復循環,直至收斂。

SPH等效修正方法屬于定點迭代方法,僅需進行一次計算資源要求較高的基準計算,在計算均勻化截面的同時存儲對應中子通量數據即可,從而有效降低超級均勻化等效修正所需的計算資源。

2.3 蒙特卡羅通量矩均勻化方法

通量矩均勻化方法在確定論程序APOLLO3的MOC快堆均勻化截面生成模塊中得到發展,并應用于快堆分析[41],用于處理總截面隨入射角的各向異性,普適于一般的堆芯輸運求解器。結合蒙特卡羅生成均勻化截面特點,本文提出了蒙特卡羅通量矩均勻化方法,以下是其簡要推導過程。

穩態波爾茨曼方程可寫為:

L(r,Ω)+T(r,Ω)=S(r,Ω)+F(r,Ω)

(14)

其中:L為泄漏;T為總反應;S為散射;F為裂變。

如假設截面不隨入射角的變化,即通量可分假設或P一致性假設,采用球諧函數表達為:

(15)

(16)

該假設普遍應用于確定論堆芯輸運求解器,有效減小了計算量。然而,在材料不連續而引起通量角度強相關的情況下將導致誤差。若考慮總截面的角度相關性,準確的方程應寫為:

(17)

入射角相關的總截面亦用球諧函數加以考慮:

(18)

為了保證所生成的截面形式在多數堆芯求解器中可用,保持總截面不隨入射角變化。

(19)

(20)

因此,等效的散射矩陣為:

(21)

以上散射矩陣通用性仍較低,因此Vidal等[41]提出了基于最小二乘法的歸并方法,最終等效的散射矩陣為:

(22)

(23)

由上述推導可知,MHT的基本原理是將總截面隨入射角的各項異性歸并入散射矩陣中。從而在考慮截面各項異性的同時,保持所生成總截面在堆芯輸運求解器中的泛用性。

2.4 方法討論

蒙特卡羅體積通量均勻化方法在許多蒙特卡羅程序中均有實現,超級均勻化等效修正方法在蒙特卡羅均勻化生成均勻化截面中已有應用。本文關鍵技術在于將通量矩均勻化方法與蒙特卡羅生成均勻化截面方法相結合,從而在通用截面形式下考慮截面隨入射角的各向異性。需要特別指出的是,本文工作較為局限,蒙特卡羅生成快堆均勻化截面還有許多問題需要深入研究,如不連續因子修正、基模修正、歷史效應處理方法等。不連續因子可保證凈流守恒,從而改善堆芯有效增殖因數和功率分布計算精度,蒙特卡羅生成不連續因子方法有待研究。基模修正理論可調整計算對象的泄漏率,從而使用考慮泄漏的燃料組件漸進能譜分布并修正群常數,在RMC[42]、Serpent[23]等蒙特卡羅程序中已有研究。快堆均勻化截面歷史效應處理方法還需進一步考慮,已有蒙特卡羅程序基于核素微觀截面,以燃耗深度和溫度為參數,可有效處理壓水堆的燃耗歷史效應問題[43]。

3 數值驗證與分析

3.1 MET-1000基準題

本文采用OECD/NEA/WPRS快堆數值對標[44]中的金屬燃料鈉冷快堆MET-1000作為對標算例驗證以上均勻化截面計算方法。MET-1000堆芯模型如圖2所示,該堆芯由72個內燃料組件、102個外燃料組件、114個反射組件、66個屏蔽組件和19個控制棒組件組成。在OECD/NEA/WPRS報告中,有7個組織使用20種不同的計算方法對MET-1000壽期初與壽期末的有效倍增因子、鈉空泡值、多普勒常數等關鍵堆芯參數進行了數值對標。近年來,該算例得到了更加廣泛的研究與應用。MET-1000算例基于先進金屬燃料快堆,有較詳實的參數描述與廣泛的研究,因此本文數值驗證與分析基于MET-1000開展。

a——徑向;b——軸向圖2 MET-1000堆芯示意圖Fig.2 Layout of MET-1000 core

本文利用OpenMC程序,基于三維全堆MET-1000模型計算其均勻化截面。堆芯計算分別采用DRAGON程序[45]中的TRIVAC擴散求解模塊[46]和OpenMC程序中的多群輸運計算模式[47]。TRIVAC是DRAGON程序包的堆芯求解模塊,包括基于擴散理論的Raviart-Thomas(RT)混合對偶有限元法和簡化球諧函數節塊法,本文采用其中的擴散求解器。同時,針對蒙特卡羅生成均勻化截面與堆芯擴散計算結合,生成SPH等效修正因子。多群蒙特卡羅方法是假設較少的堆芯輸運計算方法,因此本文采用OpenMC的多群模式作為多群輸運計算基準。針對蒙特卡羅生成均勻化截面與堆芯輸運計算結合,使用了MHT技術。

3.2 與堆芯擴散計算結合

在無控制棒插入時,連續能量精細幾何蒙特卡羅計算所得有效增殖因數為1.029 07±0.000 40,堆芯擴散計算所得有效增殖因數為1.031 63,偏差為+242 pcm。控制棒積分價值如圖3所示,若無SPH等效修正,堆芯擴散計算方法(TCP0/Diffusion)明顯高估了控制棒價值,高估隨控制棒插入逐步增加到完全插入時的+2 487 pcm,高估約為控制棒總價值的+13.5%。若對控制棒含吸收體部分的均勻化截面采用SPH等效修正方法(TCP0+SPH/Diffusion),隨著控制棒的插入其價值誤差在-141～74 pcm之間變化。在控制棒完全插入時,TCP0+SPH/Diffusion對價值預測較好,誤差比例僅為+0.35%。但在控制棒微插入時,控制棒價值較小,堆芯計算誤差可導致較大控制棒價值誤差比例,此時應考慮對相應位置燃料及其他結構物也采用SPH等效修正。

圖3 控制棒積分曲線Fig.3 S-curve of control rods

三維功率誤差最大值與最小值隨控制棒插入深度的變化如圖4所示,可見,SPH等效修正提高了功率預測的準確性。在控制棒未插入時,無SPH等效修正的三維功率分布與連續能量蒙特卡羅對比,誤差在-4.4%～+3.2%之間,主要誤差原因為控制棒吸收能力被高估導致堆芯頂部功率被低估,而堆芯底部功率被高估。經過SPH修正后,三維功率分布的誤差減小到-2.3%～2.2%之間,且誤差的空間分布更加均勻。隨著控制棒的移動,無SPH等效修正的局部功率誤差可達-10.5%,而SPH等效修正的局部功率誤差始終在-4.2%～+3.0%以內。

圖4 三維功率誤差最大值與最小值隨控制棒插入深度的變化Fig.4 Variation of relative maximal and minimal error with insertion of control

3.3 與堆芯輸運計算結合

對于MET-1000堆芯,OpenMC/OpenMC(MG)計算方法的精度如圖5所示。使用體積通量均勻化及3階勒讓德散射矩陣時多群輸運計算比基準方法高估1 192 pcm。通過入射角直接離散方法證明了截面隨入射角的各向異性是計算誤差的主導因素,約772 pcm。入射角直接離散方法難以普適確定論輸運求解器。采用2.3節中的MHT,可降低高估約698 pcm。剩余74 pcm可能是MHT均勻化未考慮其他截面類型隨入射角的各項異性。

圖5 蒙特卡羅生成均勻化截面與堆芯輸運計算結合的誤差分解Fig.5 Decomposition of reactivity bias of Monte Carlo generated homogenized cross-sections with transport core calculation

在假定散射出射各向同性的系統中,OpenMC/OpenMC(MG)仍高估反應性約938 pcm,證明散射出射角各向異性影響約為254 pcm。將散射矩陣勒讓德階數從3階提高到7階,可減小誤差120 pcm。通過進一步在多群蒙特卡羅中采用完全非均勻幾何以及更加精細的能群結構,可進一步減小誤差172 pcm。剩余128 pcm可能是因為蒙特卡羅生成多群截面時在式(5)中使用的不是對應l階的通量。

功率誤差如圖6所示。體積通量均勻化方法誤差分布不均,趨向于低估內堆芯頂部功率,而高估外堆芯底部功率,誤差在-3.63%～+4.02%間。MHT減小功率預測誤差,誤差在-2.39%～+2.76%之間,且誤差分布更加均勻。

圖6 蒙特卡羅體積通量(a)和通量距(b)均勻化方法功率誤差Fig.6 Power errors without MHT (a) and with MHT (b)

4 結論

連續能量蒙特卡羅方法計算均勻化多群截面方法可精細模擬2D/3D幾何,同時可減少共振計算的假設,可適用于復雜幾何與能譜,對先進快堆非均勻堆芯的高精度計算具有重要作用,近年來在高性能計算機不斷發展的背景下得到廣泛研究。本文聚焦使用連續能量蒙特卡羅的快堆均勻化截面計算方法研究。

1) 蒙特卡羅體積通量均勻化方法與堆芯擴散計算結合時,對無強吸收體堆芯計算精度良好。需要通過SPH等效修正方法減小對控制棒價值的高估。SPH等效修正方法亦可實現三維全堆修正,以實現特定堆芯狀態點的高精度計算。但SPH修正因子在不同狀態點的魯棒性有待深入研究。

2) 蒙特卡羅體積通量均勻化方法與堆芯輸運計算結合時,明顯高估堆芯有效增殖因數。對于MET-1000對標算例,通過定量分解分析誤差分布,明確了其主要原因是未考慮總截面隨入射角的各項異性。MHT將總截面隨入射角的各向異性并入散射矩陣中,從而提高了計算精度并普適堆芯輸運求解器,為蒙特卡羅生成適用于堆芯輸運計算截面提供了一種新的方法。該方法殘余誤差仍有待進一步消除,在類似ASTRID-CFV等三維強非均勻堆芯中的使用仍有待進一步驗證。

蒙特卡羅均勻化截面計算方法在功率分布、反饋系數、擾動狀態和計算效率等方面的特性仍需深入研究。其在強非均勻先進快堆和小型快堆方面的誤差與改進方法有待系統性研究。