基于改進缺陷因子方法的FCM燃料有效多群截面計算方法研究

易思宇，劉宙宇，尹 文

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

包覆燃料顆粒是一種廣泛應用于高溫氣冷堆的燃料類型，具有良好的安全性。近年來，為提高壓水堆中燃料的安全性，多個實驗室開展了將這種燃料顆粒彌散到碳化硅基體中形成燃料密實體應用于先進壓水堆的研究，這種燃料類型稱為全陶瓷微膠囊封裝(FCM)燃料。FCM燃料存在雙重非均勻性，包括燃料密實體中顆粒和基質(zhì)的微觀非均勻性和燃料棒之間的宏觀非均勻性，這給中子學計算帶來了挑戰(zhàn)。

對于確定論方法，雙重非均勻性給共振和輸運計算都帶來了困難。在確定論方法中，空間被離散成多個平源區(qū)，且每個平源區(qū)都被視為均勻材料，但FCM燃料中所使用的TRISO顆粒尺寸過小，難以進一步離散。為解決此問題，人們最初提出了基于體積守恒的均勻化方法，將顆粒與基質(zhì)打混處理，使雙重非均勻性系統(tǒng)轉(zhuǎn)化為單重非均勻性系統(tǒng)，但這樣會由于忽略顆粒的空間自屏效應而引入巨大誤差，為改善精度，一些基于確定論的雙重非均勻性處理方法被提出。

Luo等[1]提出了反應率等價物理轉(zhuǎn)換(RPT)方法，該方法將顆粒彌散在基質(zhì)內(nèi)的真實模型轉(zhuǎn)化為內(nèi)部為均勻燃料、外部包圍一層基質(zhì)的等效模型，并以此實現(xiàn)與真實問題相同的空間自屏效應。但這種方法需使用蒙特卡羅程序提前計算出真實問題的反應性，并根據(jù)反應性守恒來建立等效模型。

Sanchez等[2-5]和Pormraing[5]提出了基于統(tǒng)計學的方法，稱為Sanchez-Pormraning方法。該方法廣泛應用于中子學計算程序，如APPOLO、HELIOS、DeCART。其主要思想是使用概率論與數(shù)理統(tǒng)計中的數(shù)學期望的思想來建立通量的期望，然后引入弦長的概率分布來建立中子平衡方程，并將其稱為更新方程。然后使用更新方程計算獲得各平源區(qū)的均勻化截面和源項。該方法精度較高，但相比其他方法其計算效率較低。

其他研究者[6-8]提出了基于各種模型的等效均勻化方法，即缺陷因子方法。此方法易于實現(xiàn)且能有效處理雙重非均勻性，在VOSP[9]、SCALE[6]和NECP-X[10]等程序中得到了使用。基于等效一維球模型的缺陷因子方法是其中常見的方法，該方法首先在一維或二維等效模型中求解慢化方程或碰撞概率以獲得缺陷因子，然后通過缺陷因子將燃料棒中的顆粒與基質(zhì)等效均勻化，最終實現(xiàn)雙重非均勻性系統(tǒng)到單重非均勻性系統(tǒng)的轉(zhuǎn)化。盡管研究人員基于缺陷因子方法進行了很多研究，但該方法仍存在一些問題。首先，計算缺陷因子的等效模型是通過體積權(quán)重守恒來建立的，這引入了顆粒均勻分布在無限基質(zhì)中的假設，無法考慮宏觀的非均勻性；其次，非共振能區(qū)的缺陷因子是通過在三維穿透模型中保證首次飛行碰撞概率守恒獲得，但這種方法無法滿足FCM燃料中高精度計算的要求。為在計算缺陷因子時考慮宏觀非均勻性，并獲得更準確的多群缺陷因子，本文在先前研究[10]的基礎(chǔ)上提出一種改進的缺陷因子方法。

1 理論模型

1.1 概述

傳統(tǒng)的壓水堆燃料只需考慮燃料棒間的單重非均勻性，但FCM燃料還需考慮來自顆粒和基質(zhì)間的非均勻性。早期的研究中，通過簡單的體積權(quán)重守恒實現(xiàn)顆粒與基質(zhì)的均勻化，但顆粒的光學厚度大于中子的平均自由程，且顆粒內(nèi)部的通量會受空間自屏影響發(fā)生劇烈變化，簡單的體積均勻化忽略了顆粒的空間自屏效應，會帶來較大誤差。為實現(xiàn)更準確的均勻化，提出了可衡量顆?？臻g自屏效應的缺陷因子方法。缺陷因子方法的主要思想是通過缺陷因子來衡量顆粒與基質(zhì)，以及顆粒內(nèi)部的空間自屏效應，并實現(xiàn)等效均勻化。本文首先通過缺陷因子實現(xiàn)等效均勻化，然后使用Liu等[11]提出的全局-局部方法處理具有單重非均勻性的柵格系統(tǒng)中的共振自屏效應，最終獲得適用于FCM燃料的多群有效自屏截面。為保證等效均勻化時的準確性，本文提出一種改進的缺陷因子計算方法。等效均勻化是缺陷因子方法的基礎(chǔ)思想，其基于燃料密實體內(nèi)的反應率守恒，如式(1)所示。

(1)

通過對式(1)簡單轉(zhuǎn)化，可得到式(2)：

(2)

缺陷因子是子區(qū)域的通量和平均通量的比值，但真實問題的通量分布是無法提前獲得的，所以需構(gòu)建等效模型，在等效模型中獲得通量分布。

1.2 等效模型的建立

傳統(tǒng)的缺陷因子方法通過體積權(quán)重守恒來建立等效一維球模型，即顆粒尺寸保持不變，基質(zhì)外半徑通過體積比計算獲得，如式(3)所示。

(3)

其中：Rm為等效一維球模型中的基質(zhì)外半徑；Rt為TRISO顆粒的外半徑；F為TRISO顆粒的填充率。在白邊界條件下，體積權(quán)重等效模型可看作是由基質(zhì)中的顆粒組成的無限陣列，如圖1所示。該模型很明顯無法考慮真實柵格系統(tǒng)中燃料棒之間的包殼和慢化劑(圖2)。為考慮來自包殼和慢化劑的影響，本文提出一種基于顆粒丹可夫因子的等效一維球模型，以此來反映顆粒在真實柵格系統(tǒng)中的自屏效應。

圖1 燃料顆粒組成的無限陣列Fig.1 Infinite array of fuel particels

圖2 裝載FCM燃料的真實柵格系統(tǒng)(a)和等效一維球模型(b)Fig.2 Real lattice system loaded with FCM fuel (a) and equivalent 1D sphere model (b)

1) 雙重非均勻系統(tǒng)中顆粒丹可夫因子的計算

本文對顆粒丹可夫因子的計算參考Ji等[12]提出的方法，等效模型的推導都將引入顆粒包裹層和基質(zhì)的均勻化，并將均勻化后的非燃料材料統(tǒng)稱為基質(zhì)。顆粒的丹可夫因子包含兩部分：燃料棒間的丹可夫因子和燃料棒內(nèi)的丹可夫因子。燃料棒內(nèi)的丹可夫因子可通過式(4)計算獲得。

(4)

其他變量可通過以下公式獲得：

t=2(rp-rk)

(5)

(6)

(7)

(8)

(9)

其中：rp為顆粒的半徑；rk為燃料核的半徑；λ為基質(zhì)中的平均自由程；〈l〉為無限基質(zhì)中顆粒-燃料核之間的非燃料材料的平均弦長；frac′為燃料填充率；frac為顆粒填充率。

至此獲得了所有棒內(nèi)顆粒的丹可夫因子組成部分的表達式。燃料棒間顆粒的丹可夫因子可通過式(10)計算獲得。

(10)

其中，Cb為燃料棒的丹可夫因子，可將燃料棒當作黑體，然后使用中子流方法獲得。由此可通過式(4)計算獲得燃料棒內(nèi)顆粒的丹可夫因子，通過式(10)計算獲得燃料棒間顆粒的丹可夫因子，將兩者加和，即可獲得雙重非均勻性系統(tǒng)下顆粒的整體丹可夫因子，如式(11)所示。

CDH=Cintra+Cinter

(11)

2) 等效半徑的推導

假設中子是均勻且各向同性地從燃料核外表面射出，白邊界條件下等效一維球模型即可看作是由排布在基質(zhì)中的燃料核組成的無限陣列[13](圖1)，因此等效模型的丹可夫因子Ceq滿足式(12)，即：

(12)

由式(7)、(12)可獲得燃料的等效填充率frac′,eq，如式(13)所示。

(13)

而等效模型的外徑Req和燃料填充率有如下關(guān)系：

(14)

因此，Req是顆粒丹可夫因子的函數(shù)，寫作Req(Ceq)。一旦獲得雙重非均勻系統(tǒng)下顆粒的丹可夫因子，即可通過式(13)、(14)獲得Req(Ceq)。由此即可建立與真實雙重非均勻性系統(tǒng)具有相同丹可夫因子的等效一維球模型。

1.3 均勻化多群有效截面的計算

等效均勻化多群有效截面的計算分為共振能量段和非共振能量段兩部分。這兩個部分中，缺陷因子都是必需的參數(shù)，對于非共振能量段，本文提出了一種獲得缺陷因子的新方法。

1) 共振能量段的均勻化多群有效截面

共振能量段的多群有效截面計算由等效均勻化和全局-局部共振計算兩部分組成，具體流程如圖3所示。本文使用全局-局部方法進行共振計算，其中局部計算使用超細群方法，因此需計算超細群缺陷因子來均勻化共振段的超細群截面。

圖3 共振段等效多群有效截面的計算流程Fig.3 Calculation flowchart for equivalent multi-group effective cross sections in resonance energy range

(1) 等效均勻化

使用基于丹可夫因子的等效一維球模型，后續(xù)的計算與Yin等[10]提出的過程相同。通過碰撞概率法求解如式(15)所示的超細群慢化方程，以獲得超細群通量。

(15)

其中：下標g為超細群能群，i和j為子區(qū)域編號；Qjg為體積平均慢化源項；Pjig為j區(qū)產(chǎn)生的中子在i區(qū)發(fā)生第1次碰撞的概率。

超細群缺陷因子可通過式(16)獲得：

(16)

可通過式(17)獲得超細群均勻化截面：

(17)

(2) 全局-局部共振計算

將燃料棒等效均勻化后，帶有FCM燃料的柵格系統(tǒng)可被處理為典型的壓水堆柵格系統(tǒng)，即燃料棒中的燃料密實體為均勻材料。對于具有單重非均勻性的典型柵格系統(tǒng)，使用全局-局部共振方法進行自屏計算。

首先通過中子流方法獲得燃料棒的丹可夫因子Cb，在黑體假設下，有：

(18)

其中：φ0為孤立棒中燃料區(qū)的通量；φ為柵格系統(tǒng)中燃料區(qū)的通量。φ可通過MOC方法求解單群固定源問題獲得。在黑體假設下，φ0等于中子從燃料到慢化劑的逃脫概率，孤立棒問題的逃脫概率可通過Carlvik方法獲得。

然后根據(jù)丹可夫因子守恒搜索慢化劑的外徑，如式(19)所示。

(19)

其中：rM為慢化劑的外半徑；Pe為中子從燃料棒中逃脫的概率；Pf-m(rM)為在燃料中產(chǎn)生的中子在慢化劑中發(fā)生第1次碰撞的概率。

至此即可建立針對每根燃料棒的一維柵元模型，最后在等效一維圓柱柵元模型中通過超細群方法獲得共振能量段的多群有效截面。

2) 非共振能量段的均勻化多群有效截面

非共振能量段無需進行共振計算，所以該能量段的多群截面只需使用多群缺陷因子進行等效均勻化即能獲得均勻化多群有效截面。具體流程如圖4所示。

圖4 非共振段等效多群有效截面的計算流程Fig.4 Calculation flowchart for equivalent multi-group effective cross section in non-resonance energy range

在之前的研究中，非共振能量段缺陷因子的計算參考She等[8]的方法，其主要思想是在一個三維穿透模型中保證首次飛行碰撞概率守恒，所以只需計算中子在各子區(qū)域的首次飛行碰撞概率和穿透概率。She等的方法有較高的計算效率，但在填充率較低或較高的情況下計算精度不足。

為獲得更準確的缺陷因子，本文采用在等效一維球模型中進行特征值計算的方法獲得通量分布。首先需獲得進行特征值計算必需的多群有效截面，因為前文已獲得了共振區(qū)的超細群通量，所以只需歸并截面即可獲得有效共振截面，而非共振區(qū)可直接使用典型能譜下的有效截面。然后是求解特征值問題，本文使用碰撞概率法進行求解。

2 數(shù)值結(jié)果

從單柵元問題、二維單組件問題、三維全堆問題等多方面對本文提出的方法進行驗證。并在單柵元規(guī)模下，對不同顆粒尺寸、填充率、帶可燃毒物等工況進行驗證。參考解來自蒙特卡羅程序，為獲得準確結(jié)果，蒙特卡羅計算時采用顯示建模的方法。蒙特卡羅程序使用的連續(xù)能量庫和NECP-X使用的截面庫皆基于ENDF/Ⅶ.0。NECP-X的多群結(jié)構(gòu)為WIMSD69，其中共振能群為13～45群。各區(qū)域的溫度，如果沒有提及，均為600 K。

2.1 單柵元問題

單柵元和TRISO顆粒的幾何布置示于圖5，其尺寸和材料列于表1。材料的核子密度列于表2。在蒙特卡羅計算中，使用30萬個粒子，共500代進行模擬，其中前200代為非活躍代。由于TRISO顆粒在尺寸、填充率、富集度等方面都有多種選擇，為方便討論，設計了表3所列的一些常見工況的顆粒參數(shù)。

圖5 FCM燃料柵元和TRISO顆粒Fig.5 FCM pin cell and TRISO particle

表1 FCM燃料柵元的尺寸和材料Table 1 Dimension and material of FCM cell

表2 材料的核子密度Table 2 Atomic number density of each material

表3 一些常見工況的顆粒設計參數(shù)Table 3 Particle design parameters for some common cases

1) 顆粒尺寸和填充率

為更全面地反映方法的幾何適應性，在工況1的基礎(chǔ)上通過改變?nèi)剂习霃胶褪杷蔁峤馓紝雍穸仍O計了表4所列的系列不同顆粒尺寸工況，其中包含了工況2～5。驗證在15%、25%、35、45%的填充率下進行。

表4 用于幾何適應性測試的TRISO顆粒尺寸Table 4 TRISO particle sizefor geometric adaptivity test

為說明方法的改進，分別采用新舊方案對問題進行計算。其中新方案為使用本文提出的基于丹可夫因子的等效一維球模型與求解特征值問題獲得非共振區(qū)缺陷因子的計算過程，而舊方案為使用體積權(quán)重等效一維球模型和在三維穿透模型中計算非共振區(qū)缺陷因子的計算過程。兩方案計算結(jié)果與蒙特卡羅結(jié)果的偏差列于表5。

從表5可見，舊方案燃料核直徑為0.08 cm、填充率為0.45的工況中，kinf的絕對偏差超過300 pcm，即對于有較大燃料核直徑、較大填充率和較小疏松熱解碳層厚度的工況，舊方案無法獲得準確的計算結(jié)果。而新方案中這幾個工況都得到了極大的改善。對其中的工況1、3、4進行深入分析，具體計算結(jié)果列于表6。

表5 新舊方案中不同顆粒尺寸和填充率下kinf的絕對偏差Table 5 Absolute bias of kinf in different particle sizes and packing fractions with new and old schemes

表6 部分工況的kinf和等效半徑Table 6 kinf and equivalent radius in part of cases

在新方案中，基于丹可夫因子的等效一維球模型能考慮燃料棒外慢化劑的影響，因此相比舊方案有更大的基質(zhì)外徑來反映更加真實的屏蔽效應。將新舊方案中的238U和235U有效共振自屏截面與蒙特卡羅結(jié)果進行比較，結(jié)果如圖6所示。

圖6 238U和235U共振能量段有效吸收截面相對偏差比較Fig.6 Comparison of relative bias of effective absorption cross section of 238U and 235U

由圖6可見，相比舊方案，新方案在20～27群的238U吸收截面明顯減小且更接近參考值，整個共振能量段最大相對偏差小于2.5%。235U的吸收截面有類似情況，20～27群的截面得到了改善，雖然最后幾群的截面更加偏離參考值，但變化很小，且最大相對偏差都小于1%。新方案中的等效模型更大的基質(zhì)外徑使得其共振段吸收截面減小，增強了空間自屏效應，使燃料核中的缺陷因子減小，而缺陷因子的減小會反映到吸收截面上，最終吸收截面減小使得kinf增大。這也能解釋在舊方案中當燃料核直徑增大時會出現(xiàn)kinf往負方向偏離參考值的趨勢，且這種趨勢會隨著顆粒填充率的增大、疏松熱解碳層變薄更加明顯。這是因為增大燃料核直徑和填充率、減小疏松熱解碳層厚度都會使燃料占比增大，這時體積權(quán)重模型導致的忽略慢化劑所帶來的影響就會被放大。

除共振能量段的改進，新方案中熱能區(qū)的宏觀截面同樣也得到了改善，如圖7所示，燃料密實體中宏觀裂變截面和宏觀總截面相比舊方案都有所減小，且更接近參考值，這進一步說明了新方案中獲得多群缺陷因子方法的準確性。

圖7 燃料密實體宏觀裂變截面和總截面相對偏差比較 Fig.7 Comparison of relative bias of macro fission and total cross sections in fuel lump

2) 燃料溫度和富集度

工況1、6、7是關(guān)于溫度和富集度的敏感性測試，kinf的計算結(jié)果列于表7?？梢姡诓煌瑴囟群透患认耴inf與參考值符合得較好。不同工況下235U和238U的有效截面相對偏差比較示于圖8?？梢?，238U有效截面的最大相對偏差小于2%，235U有效截面的最大偏差小于1%，均與蒙特卡羅計算所得的參考解符合較好。

圖8 新方案中不同溫度及富集度下238U和235U共振能量段有效吸收截面相對偏差比較Fig.8 Comparison of relative bias of effective absorption cross section of 238U and 235U in different termperatures and enrichments in new scheme

表7 不同溫度和富集度下kinf和絕對偏差Table 7 kinf and absolute bias in different temperatures and enrichments

3) 燃料和毒物燃耗

為驗證本文提出的方法在燃耗計算時的準確性，基于工況4進行燃耗計算，并在工況4的基礎(chǔ)上設計了帶毒物的工況。在毒物工況中，使用QUADRISO顆粒替代TRISO顆粒。QUADRISO顆粒即在TRISO顆粒的燃料核外額外包裹1層毒物，本文使用10 μm厚的Er2O3。相關(guān)的顆粒尺寸列于表8。

表8 QUADRISO顆粒尺寸Table 8 Size of QUARISO particle

整個燃耗過程從0 GW·d/tU到72.5 GW·d/tU。整個燃耗過程中kinf及其絕對偏差的變化示于圖9。對于沒有毒物的工況，kinf的絕對偏差從最開始的-110 pcm到20 GW·d/tU的10 pcm，隨后緩慢增長到72.5 GW·d/tU的30 pcm。其中最大絕對偏差為-120 pcm。對于有毒物的工況，kinf的絕對偏差從最開始的-80 pcm增長到最后的-230 pcm。相比沒有毒物的工況，毒物的出現(xiàn)使得偏差略有增加，但仍在可接受范圍內(nèi)。

圖9 無毒物和有毒物燃料柵元在燃耗情況下的kinf及其絕對偏差Fig.9 kinf and its absolute bias in depletion for case without and with poison

2.2 二維組件問題

本文設計了2個二維的13×13組件，第1個組件由帶TRISO顆粒的燃料棒和導向管組成，第2個組件將第1個組件中部分燃料棒中的TRISO顆粒替換為QUADRISO顆粒，幾何布置如圖10所示。TRISO顆粒尺寸同工況4，QUATRISO顆粒尺寸如表8所列。材料組成和幾何參數(shù)列于表9，材料的核子密度如表2所列。

圖10 無毒物組件和有毒物組件的幾何布置Fig.10 Geometry configuration of fuel lattice without and with poison

表9 導向管尺寸和材料Table 9 Dimension and material of guide tube

將kinf和功率分布皆與蒙特卡羅結(jié)果進行比較，結(jié)果見表10和圖11?？煽吹?，兩個工況下kinf皆與參考值符合較好。對于不含毒物的組件，其柵元功率最大相對偏差為0.31%；在含毒物的組件中，柵元功率最大相對偏差為0.32%。

圖11 有無毒物組件的功率分布及相對偏差Fig.11 Power distribution and relative bias for lattice with and without poison

表10 2D組件問題的kinf及其絕對偏差 Table 10 kinf and absolute bias of FCM lattice

2.3 三維全堆問題

根據(jù)Awan等[14]設計的裝載FCM燃料的小型模塊化堆，本文建立了1個有4類共37個組件的三維全堆問題，其中2類只裝載了TRISO顆粒，其燃料分別為14.3%和10%的UC。另外2類中分別將前兩者中部分燃料棒的TRISO顆粒替換為含有5 μm和10 μm厚毒物層的QUADRISO顆粒。具體組合列于表11。組件由13×13柵元構(gòu)成，具體幾何布置如圖10所示。顆粒和柵元的尺寸基于2.1節(jié)中的工況1，所用材料如表2所列。

表11 4類組件的配置Table 11 Configuration of 4 kinds of assemblies

為簡化計算，將巨大的環(huán)形反射層簡化為一圈20 cm厚的反射層，軸向反射層被簡化為上下各20 cm厚的反射層，如圖12所示。至此建立起了一個由160 cm高活性區(qū)、40 cm厚反射層組成的簡化全堆。在NECP-X計算中，將模型分為20層，每層10 cm。

圖12 全堆和簡化全堆的俯視圖和側(cè)視圖Fig.12 Plan view and side view of whole and simplified FCM fuel core

使用本文提出的改進缺陷因子方法進行計算，以蒙特卡羅結(jié)果為參考。在特征值上，采用改進缺陷因子方法所得的keff為1.227 81，蒙特卡羅結(jié)果為1.226 67，兩者間的絕對偏差為114 pcm，符合得很好。徑向和軸向功率分布與參考解的比較如圖13、14所示。軸向功率的最大相對偏差出現(xiàn)在最頂部，為2.28%，但兩端的功率較低，所以偏差是可接受的。組件徑向功率的相對偏差最大為-1.13%，出現(xiàn)在最中心的組件，整體符合較好。圖15為有最大棒功率相對偏差的組件中的功率分布，其中棒功率最大相對偏差為1.88%，整體與參考解符合得較好。

圖13 歸一化軸向功率分布和相對偏差Fig.13 Normalized axial power distribution and relative bias

圖14 全堆徑向組件功率分布和相對偏差Fig.14 Radial assembly power distribution and relative bias in core

圖15 組件內(nèi)徑向柵元功率分布和相對偏差Fig.15 Radial pin power distribution and relative bias in assembly

3 結(jié)論

本文提出了一種改進的缺陷因子方法用于獲得FCM燃料中的有效多群截面。首先提出了基于丹可夫因子的等效一維球模型來代替體積權(quán)重模型，傳統(tǒng)的體積權(quán)重模型忽略了燃料棒之間的非均勻性，因此當燃料密實體中燃料體積占比較大時會引入較大誤差，而較大的燃料體積占比在FCM燃料設計中很常見，新的等效模型能考慮燃料棒間的非均勻性，可獲得更精確的計算結(jié)果。其次在非共振段處理中，采用在等效模型中進行特征值計算的方法來獲得多群缺陷因子，相比穿透模型的碰撞概率計算，該方法可獲得更精確的結(jié)果。

該方法已在NECP-X中實現(xiàn)，并在不同工況和不同規(guī)模的問題中進行了驗證。在單柵元問題中，計算了不同顆粒尺寸、填充率等工況，kinf偏差均在100 pcm以內(nèi)，有效多群截面除少數(shù)能群相對偏差接近2.5%，大多數(shù)能群都在±1%以內(nèi)。在二維組件問題中，kinf偏差在103 pcm以內(nèi)，柵元功率最大相對偏差為0.32%。在三維全堆問題中，keff偏差為114 pcm，除兩端功率較低外，其余活性區(qū)層軸向功率相對偏差均在1.5%以內(nèi)，徑向組件功率最大相對偏差為1.14%，徑向積分柵元功率相對偏差為1.88%。所有計算結(jié)果與參考解符合較好，證明本方法能有效處理FCM燃料中的雙重非均勻性，并獲得了精確的有效多群截面。