張 斌，李云召，吳宏春，王冬勇，劉 勇，于穎銳

(1.中國核動力研究設計院 核反應堆系統(tǒng)設計技術重點實驗室，四川 成都 610213；2.西安交通大學 核科學與技術學院，陜西 西安 710049)

20世紀70年代以后，得益于粗網節(jié)塊法[1-2]的迅速發(fā)展，以組件均勻化理論[3-4]和粗網節(jié)塊方法為理論框架的兩步法計算方案逐漸成為壓水堆工程計算中普遍采用的燃料管理中子學數(shù)值計算方法。伴隨著粗網節(jié)塊法的大量工程應用，適用于組件均勻化的均勻化理論得到了長足的發(fā)展?；诰鶆蚧话阍韀5]，結合早期有限差分堆芯計算的傳統(tǒng)均勻化方法，Smith基于等效均勻化理論提出了能用于工程計算的廣義等效均勻化方法[6-7]，即在等效均勻化參數(shù)中引入不連續(xù)因子以達到各節(jié)塊的守恒條件，此方法很好地滿足了絕大部分商用壓水堆的工程應用需求。

隨著科學研究的不斷深入以及計算條件和對核設計計算精度要求的不斷提高，傳統(tǒng)的兩步法計算方案面臨著越來越嚴峻的挑戰(zhàn)。為提高堆芯設計計算精度且滿足對新型反應堆堆芯計算的精度要求，基于均勻化理論的改進型兩步法計算方案全堆芯Pin-by-pin計算[8]成為下一代堆芯數(shù)值計算方法的研究熱點。全堆芯Pin-by-pin計算減少了堆芯計算過程中的近似與假設，能更加精細地考慮堆芯布置的非均勻性，直接求出單棒功率分布，便于堆芯燃料管理計算和相應的堆芯安全分析。

有別于不連續(xù)因子在傳統(tǒng)組件均勻化計算中的廣泛應用，壓水堆堆芯Pin-by-pin均勻化計算中除廣義等效均勻化方法外，超級均勻化(SPH)方法亦是其主流的均勻化技術之一。SPH方法由Kavenoky[9]于20世紀80年代提出。Herbet基于均勻化原理中特征值、各能群反應率及中子泄漏率3個守恒關系不斷對SPH方法進行完善，最終使其成為Pin-by-pin計算中的主流均勻化方法[10-11]。SPH方法的核心思想是放寬中子泄漏率守恒這一約束條件，通過SPH因子直接調整柵元均勻化截面使得均勻化前后各能群反應率保持守恒。這種方法不需要均勻化少群常數(shù)中保存額外的均勻化常數(shù)(不連續(xù)因子)，在有很大內存需求量的全堆芯Pin-by-pin計算中具有明顯優(yōu)勢。近年來，國內外學者對SPH方法在程序實現(xiàn)、能譜干涉修正、中子泄漏效應處理等方面進行了大量研究[12-19]，日本AEGIS/SCOPE2、德國DYN3D等程序已采用SPH方法作為其Pin-by-pin計算的均勻化方法。SPH方法不受具體堆芯計算模型限制，既可保證均勻化前后輸運-擴散守恒，也可保證輸運-輸運守恒，對于不同的均勻化過程體現(xiàn)在SPH因子修正少群常數(shù)時的不同，是壓水堆堆芯Pin-by-pin計算中首選的均勻化技術之一。

本文采用SPH方法計算燃料組件的Pin-by-pin等效均勻化參數(shù)，針對存在中子泄漏現(xiàn)象的反射層參數(shù)計算研究與空間泄漏相關的SPH方法。在保證反應率守恒的基礎上，同時保證各柵元各能群的中子泄漏率守恒，解決存在中子泄漏時SPH因子迭代計算的不收斂問題，并分析壓水堆堆芯Pin-by-pin計算中應用SPH因子的堆芯計算精度。

1 方法

1.1 SPH方法

均勻化少群截面的計算通常采用通量體積權重方法，如式(1)所示。

x=a,f,s,…

(1)

式中：g、h為能群編號；i為柵元編號；V為柵元體積；φ為反應堆中子通量密度，cm-2·s-1；Σ為宏觀截面，cm-1；上標hom表示均勻化后，het表示均勻化前。

式(1)假設柵格計算得到的非均勻中子通量密度等于均勻中子通量密度。此假設在只有均勻化少群截面與擴散系數(shù)作為等效均勻化常數(shù)的情況下是不成立的。SPH方法通過調整均勻化少群截面使得式(1)成立。

(2)

在保證反應率守恒的條件下，有：

(3)

整理可得：

(4)

聯(lián)立式(2)、(4)可得SPH因子的計算公式：

(5)

式(5)中非均勻柵元平均中子通量密度由柵格高階輸運計算直接獲得；均勻柵元平均中子通量密度的求解須與下游堆芯計算中所使用的中子學求解方法保持一致。SPH因子通過均勻柵元平均中子通量密度計算得到，而獲得中子通量密度計算所需的均勻化少群常數(shù)與擴散系數(shù)又需通過SPH因子的修正，因此SPH因子的計算是一反復迭代的計算過程，具體計算步驟如下。

1) 在單組件全反射邊界條件下進行高階輸運計算，產生非均勻柵元的平均中子通量密度、各柵元的少群均勻化截面和擴散系數(shù)。

2) 初始化SPH因子。

3) 利用SPH因子修正步驟1中產生的少群均勻化截面與擴散系數(shù)，并用下游堆芯計算中所使用的中子學求解方法對全反射邊界條件下柵元均勻化的組件問題進行計算，產生均勻化柵元的平均中子通量密度。

4) 根據(jù)式(5)，利用步驟1、3中的中子通量密度計算結果更新SPH因子。在更新SPH因子前，需對步驟3中產生的平均中子通量密度進行歸一化處理，使得整個組件每一能群下的平均中子通量密度與步驟1中的中子通量密度保持一致，歸一化公式如下：

(6)

5) 判斷SPH因子是否收斂，收斂條件如式(7)所示，若滿足上述收斂條件，則計算停止；否則，進入下一次迭代計算，返回步驟3。

(7)

式中，n為迭代次數(shù)。

1.2 與空間泄漏相關的SPH方法

在柵格計算中，對單燃料組件采用全反射邊界條件進行計算，整個問題在均勻化前后不存在中子泄漏，即利用SPH因子保證反應率守恒的同時還保證了中子泄漏率的守恒，燃料組件特征值也因此而守恒。當SPH方法應用于真空邊界，即反射層參數(shù)計算時，SPH因子迭代計算的收斂性無法保證。這是因為在真空邊界下入射中子流為0的邊界條件無法保證出射中子流密度的守恒，迭代過程中改變均勻化少群截面的同時，影響著整個問題的中子泄漏率及其特征值。因此，此類邊界條件存在時，SPH因子迭代計算的收斂性無法保證，且在收斂的情況下依舊無法保證整個問題中子泄漏率和特征值的守恒。

針對此類邊界下SPH因子的迭代計算，與空間泄漏相關的SPH方法能使迭代計算很好地收斂并保證各柵元中子泄漏率的守恒，其核心思想是以柵元級少群反照率邊界條件保證均勻化前后各柵元的中子泄漏率守恒。與空間泄漏相關的SPH方法將柵元級少群反照率的計算更新加入到SPH因子的迭代求解過程中，具體計算流程如下。

1) 真空邊界條件下柵格高階輸運計算，產生非均勻柵元的平均中子通量密度、各柵元的少群均勻化截面和擴散系數(shù)、組件在真空邊界上各柵元的入射和出射中子流密度和少群反照率。反照率計算公式如下：

(8)

式中：βi,g為第i個柵元第g群反照率；h為非均勻輸運計算的細群編號；J為中子流密度；角標in、out分別表示入射、出射。

2) 初始化SPH因子并以式(8)計算的柵元級少群反照率作為Pin-by-pin均勻化組件的首次邊界條件。

3) 利用SPH因子修正少群均勻化截面與擴散系數(shù)，并用下游堆芯計算中所使用的中子學求解方法計算柵元級少群反照率邊界下的柵元均勻化組件問題，產生均勻化柵元的平均中子通量密度及邊界上各柵元的入射和出射中子流密度。

4) 按照式(6)歸一化處理平均中子通量密度，根據(jù)式(5)更新SPH因子，迭代過程中真空邊界上的柵元級少群反照率的更新如式(9)、(10)所示。

(9)

(10)

5) 判斷SPH因子是否收斂，收斂條件為式(7)。若滿足上述收斂條件，則計算停止；否則，進入下一次迭代計算，返回步驟3。

與空間泄漏相關的SPH方法通過SPH因子保證各能群各柵元的反應率守恒的同時，通過柵元級少群反照率邊界保證各群中子泄漏率的守恒，燃料組件特征值也因此而守恒，計算流程如圖1所示。

圖1 與空間泄漏相關的SPH方法計算流程圖Fig.1 Calculation flow of spatial neutron-leakage dependent SPH method

2 數(shù)值驗證

為驗證SPH方法在壓水堆堆芯Pin-by-pin計算中的計算精度，以KAIST基準題2A問題下的二維全堆芯問題進行驗證與分析，堆芯布置和幾何結構如圖2所示。

圖2 KAIST基準題二維全堆芯布置和幾何結構Fig.2 Layout and geometry of KAIST benchmark 2D whole core problem

采用堆芯Pin-by-pin計算和堆芯組件均勻化粗網計算兩種方案進行堆芯計算并加以分析和比較。

1) 依據(jù)SPH方法及改進型SPH方法的理論模型，針對擴散方法(P1)和簡化球諧函數(shù)方法(SP3)利用SPH因子分別產生燃料組件及反射層組件的兩群和七群的等效均勻化少群常數(shù)，并在指數(shù)函數(shù)展開節(jié)塊法中進行全堆芯Pin-by-pin計算。

2) 根據(jù)傳統(tǒng)組件均勻化計算方法，產生兩群帶組件不連續(xù)因子(ADF)的組件等效均勻化常數(shù)并采用先進變分節(jié)塊方法進行堆芯計算。變分節(jié)塊方法在空間處理上直接針對中子通量密度分布進行離散，因此不需要功率重構過程而能直接得到中子通量密度分布，通過形狀因子的求解即可得到堆芯棒功率分布。

二維堆芯參考解由特征線一步法計算獲得，特征值參考解為0.979 48，棒功率分布如圖3所示。堆芯兩步法計算的特征值誤差和棒功率誤差列于表1。由表1可知：堆芯Pin-by-pin簡化球諧函數(shù)的計算精度普遍高于組件均勻化的堆芯粗網節(jié)塊的計算精度；除因誤差抵消導致計算精度提高的兩群擴散計算外，堆芯Pin-by-pin擴散計算的計算精度均低于組件均勻化的堆芯粗網節(jié)塊計算。因此，簡化球諧函數(shù)方法是堆芯Pin-by-pin計算的必要條件；隨著能群數(shù)的增加，堆芯Pin-by-pin簡化球諧函數(shù)的計算精度顯著提高。七群簡化球諧函數(shù)計算中特征值誤差小于50 pcm，棒功率相對誤差的最大值(MAX)低于4%，棒功率相對誤差的平均值(AVG)、均方根值(RMS)、平均相對值(MRE)都在1%以下；而在組件均勻化的粗網節(jié)塊計算中，特征值誤差為219 pcm，最大棒功率相對誤差接近10%。因此，七群是適用于堆芯Pin-by-pin計算精度的能群數(shù)目。

組件均勻化粗網節(jié)塊計算的棒功率誤差分布如圖4所示，堆芯Pin-by-pin擴散計算和簡化球諧函數(shù)計算的棒功率相對誤差分布分別如圖5、6所示。由圖5、6可發(fā)現(xiàn)，最大棒功率相對誤差出現(xiàn)在組件與組件交界面以及燃料柵元與反射層交界面附近。這是堆芯內組件真實環(huán)境與單組件計算中全反射邊界條件之間的不同所致。

圖3 KAIST基準題二維全堆芯棒功率分布Fig.3 Pin-power distribution of KAIST benchmark 2D whole core problem

表1 堆芯兩步法計算的特征值誤差和棒功率相對誤差Table 1 Two-step calculation results of kinf and pin-power relative error

圖4 組件均勻化粗網節(jié)塊計算的棒功率相對誤差分布Fig.4 Pin-power relative error distribution of assembly-homogenized nodal calculation

3 結論

本文在壓水堆堆芯Pin-by-pin計算中采用SPH方法作為均勻化技術，對燃料組件計算了SPH因子，產生了Pin-by-pin等效均勻化參數(shù)；針對存在中子泄漏現(xiàn)象的反射層組件，采用與空間泄漏相關的SPH方法產生了反射層組件的等效均勻化參數(shù)。基于KAIST基準題，分析了壓水堆堆芯Pin-by-pin計算中應用SPH因子的堆芯計算精度，與傳統(tǒng)組件均勻化計算方法相比，應用SPH方法的壓水堆堆芯Pin-by-pin的計算精度更高?；谄呷旱暮喕蛑C函數(shù)方法是高精度堆芯Pin-by-pin計算方案的實施策略，與KAIST基準題一步法計算結果相比，其特征值誤差小于50 pcm，最大棒功率相對誤差小于4%，棒功率相對誤差的AVG、MRS、MRE都在1%以下。

圖5 堆芯Pin-by-pin擴散計算的棒功率相對誤差分布Fig.5 Pin-power relative error distribution of Pin-by-pin diffusion calculation

圖6 堆芯Pin-by-pin SP3計算的棒功率相對誤差分布Fig.6 Pin-power relative error distribution of Pin-by-pin SP3 calculation