袁顯寶 劉曾豪 張彬航 張永紅 唐海波 楊森權 肖云龍

1（三峽大學 機械與動力學院 宜昌 443002）

2（三峽大學 理學院 宜昌 443002）

3（中核武漢核電運行技術股份有限公司 武漢 430223）

反應堆數(shù)值計算是一種利用計算機模擬技術來預測反應堆在運行過程中關鍵安全參數(shù)的方法。它為核能系統(tǒng)的研究開發(fā)和安全分析提供了重要手段和依據(jù)。根據(jù)不同的需求可以選擇不同的反應堆數(shù)值計算方法，蒙特卡羅方法或基于確定論的全堆芯一步法可以提供準確的計算結果，因其較為高昂的計算成本，限制了其在工程中的廣泛應用。基于“組件-堆芯”的兩步法能夠同時兼顧精度與效率［1］，仍是工程計算中的主要方法之一。隨著堆芯物理計算向更高保真度的發(fā)展，兩步法受制于組件均勻化計算和堆芯擴散近似計算中引入的簡化假設，其對于精細化模型建立的正確性和計算的準確性需要進一步研究。BEAVRS是由麻省理工學院提出的三維全堆芯壓水堆基準題［2］，相比于其他基準題，BEAVRS包含了詳細的幾何結構和真實的實測數(shù)據(jù)，更具參考價值。

DRAGON/DONJON 是由加拿大蒙特利爾綜合理工大學開發(fā)的反應堆物理計算程序，程序采用兩步法進行計算，其中DRAGON 為組件計算程序，DONJON 為堆芯計算程序［3］。本文基于DRAGON/DONJON對三維全堆芯壓水堆基準題BEAVRS 2.02（Benchmark for Evaluation And Validation of Reactor Simulations Rev.2.0.2）［4］進行建模分析，同時對于非對稱燃料組件研究了三種不同的均勻化方案，計算熱態(tài)零功率下的各項參數(shù)，并與各項實驗值進行比較分析，包括：臨界硼濃度、控制棒價值、裂變率分布，用以研究兩步法在精細化模型計算中的適用性。

1 基準題概述

BEAVRS 2.02 由美國麻省理工學院發(fā)布，該基準題詳細描述了堆芯的幾何結構和運行參數(shù)。BEAVRS 反應堆總高度為460 cm，熱態(tài)零功率狀態(tài)下運行功率為25 MW。

該堆芯包含193 個燃料組件，組件寬度為21.417 28 cm，每個組件平均分為17×17 個柵元，柵元寬度為1.259 84 cm，每個組件包括264根燃料棒，1根是位于中心位置用于測量235U裂變率的儀表管，其余24 根可用于放置控制棒組件、中子源組件、可燃毒物組件或阻力塞組件。燃料棒根據(jù)235U富集度的不同分為1.6%、2.4%和3.1%三種不同類型［5-7］。

BEAVRS 堆芯包含5 種毒物布置的組件，根據(jù)組件中的燃料235U 富集度和組件在堆芯中的位置，每個組件的可燃毒物棒所占據(jù)的導管的數(shù)量和相對位置是不同的。按照數(shù)量的不同可燃毒物棒分別為6BA、12BA、15BA、16BA和20BA，如圖1所示，其中12BA、16BA 和20BA 是對稱組件，6BA、15BA 是非對稱組件。

圖1 可燃毒物棒布置Fig.1 Layout of burnable absorber assembly model

2 計算模型與方法

2.1 組件群常數(shù)計算

基于DRAGON 的組件群常數(shù)計算方法如圖2所示。中子截面庫采用基于ENDF/B 8.1制作的295群多群數(shù)據(jù)庫［8］，核素的有效自屏截面計算采用廣義斯丹姆勒方法（SHI 模塊），通過特征線方法（MCCGT 模塊）完成一步法二維輸運計算。使用FLU模塊中B1模型和SIGS模型分別完成基模修正以及擴散系數(shù)求解；并群后的組件群常數(shù)為2 群結構，以4.0 eV為界限。

圖2 基于DRAGON的組件群常數(shù)計算流程Fig.2 Calculation flowchart of a few-energy-group cross-sections based on DRAGON

如圖3 所示，分別是含有12BA、16BA 和20BA的對稱燃料組件模型，綠色部分為燃料棒，藍色部分為可燃毒物棒，黃色部分表示空的導管。燃料組件采用全反射邊界的二維組件模型，棒束、包殼和水隙均顯示建模。對于非燃料組件包括反射層、冷卻劑以及某些組件軸向段，采用“超柵元”模型執(zhí)行均勻化計算，以獲得精確的均勻化常數(shù)。如圖4所示，將一個或多個燃料組件與反射層、冷卻劑等待測介質（最右邊純黃色區(qū)域）作為整個計算區(qū)域完成輸運計算，待測組件外側采用真空邊界條件，其余部分依然采用全反射邊界條件。

圖3 對稱燃料組件模型示意圖（彩圖見網(wǎng)絡版）Fig.3 Diagram of symmetrical fuel assembly model (color online)

圖4 “超柵元”模型示意圖Fig.4 Diagram of supercell model

2.2 全堆芯擴散計算

基于DONJON 的三維堆芯計算方法如圖5 所示。堆芯模型最外側采用真空邊界條件，堆芯模型在徑向上以組件為單位劃分節(jié)塊，采用節(jié)塊法對每個組件的軸向部分進行網(wǎng)格細分。模型的軸向結構如表1所示，包括帶格架和不帶格架兩種類型，在軸向兩種類型的模型交錯布置，控制棒下層材料為B4C，上層材料為Ag-In-Cd。每個組件在軸向上劃分了29 段。對不同燃料組件軸向網(wǎng)格數(shù)的堆芯模型進行計算后發(fā)現(xiàn)，當軸向節(jié)塊數(shù)大于等于29 時，進一步增加節(jié)塊數(shù)對有效增殖因數(shù)的影響可以忽略。因此，當燃料組件軸向被劃分為29 段時，堆芯計算模型可以滿足網(wǎng)格無關性的要求。徑向上，每個組件劃分為2×2 個節(jié)塊，接著將組件計算產(chǎn)生的每個節(jié)塊上的少群常數(shù)放入相應的節(jié)塊。

表1 組件軸向布置Table 1 Distribution of the assembly axis

圖5 基于DONJON的三維堆芯計算流程Fig.5 Calculation flowchart of the 3D full core based on DONJON

2.3 非對稱燃料組件的均勻化方案

在BEAVRS反應堆中包含5種可燃毒物布置的組件，分別為6BA、12BA、15BA、16BA和20BA。其中，6BA 和15BA 類型的燃料組件并不是呈中心對稱分布的，傳統(tǒng)的均勻化方法是將整個組件的材料打混均勻化，忽略了不對稱燃料組件在空間上的非均勻性，從而引入計算誤差［9-10］。

為了充分考慮6BA和15BA組件在空間上的非均勻性對計算結果的影響，本文采用了三種均勻化計算方案：

方案1：傳統(tǒng)均勻化方案，即建立完整的組件模型并對整個組件進行均勻化計算，如圖6（a）所示；

圖6 均勻化方案 (a) 方案1，(b) 方案2，(c) 方案3Fig.6 Homogenization scheme 1 (a) Scheme 1,(b) Scheme 2,(c) Scheme 3

方案2：一次多區(qū)均勻化方案，即建立完整的組件模型并將組件等分為4 個區(qū)域，分區(qū)域進行均勻化計算，如圖6（b）所示，獲得4 種均勻化少群常數(shù)，堆芯計算時將每個組件劃分為2×2 個節(jié)塊，將均勻化少群截面按照區(qū)域順序依次填入對應節(jié)塊中；

方案3：多次多區(qū)均勻化方案，將組件模型等分為1/4 組件計算進行多次均勻化計算，如圖6（c）所示，其余設置和方案2一致。

3 計算結果及分析

3.1 燃料組件計算

對于具有全反射邊界條件的燃料組件，設定溫度為566 K，硼濃度為9.75×10-4，分別計算在1.6%、2.4%和3.1%的富集度以及不同數(shù)量的可燃毒物棒對組件kinf的影響。計算結果如表2 所示，并與韓國首爾大學和韓國原子能研究院開發(fā)的確定論程序nTRACER和蒙特卡羅程序McCARD的計算值進行了比較［11-12］。結果表明，本文計算值與兩個參考值誤差均在5×10-3以內，滿足計算精度要求，驗證了組件模型計算的正確性。

表2 燃料組件kinf計算Table 2 Calculation of kinf for the fuel assembly

3.2 臨界硼濃度

采用三種均勻化方案分別計算控制棒組全提、控制棒D插入堆芯、控制棒C和控制棒D插入堆芯、控制棒A、B、C、D 插入堆芯以及控制棒A、B、C、D、SE、SD、SC 插入堆芯等5 種狀態(tài)下臨界硼濃度的變化情況，并與BEAVRS 提供的實驗值進行比較。如表3所示，與基準參考值相比，三種均勻化方案的計算誤差均在5×10-5以內。如表4 所示，硼濃度的keff計算值與keff的臨界值相比誤差在1×10-5以內，符合精度要求。進一步比較三種均勻化方案的計算結果，發(fā)現(xiàn)三種均勻化方案的計算結果精度相當，初步驗證了三維穩(wěn)態(tài)堆芯計算模型的正確性與可靠性。

表3 熱態(tài)零功率臨界硼濃度計算Table 3 Calculation of the critical boron concentration for HZP

表4 臨界硼濃度對應keff計算值Table 4 Calculation of the critical boron concentration for keff

3.3 控制棒價值

本文進一步對熱態(tài)零功率狀態(tài)下不同控制棒組的價值進行了計算。對于控制棒價值的計算，通常對每組棒進行兩次獨立的臨界本征值計算，第一次是完全抽出棒組，第二次是完全插入棒組，兩次計算應使用相同的硼濃度。在通過兩次獨立計算獲得臨界本征值之后，根據(jù)式（1）計算控制棒價值。

式中：kin是控制棒組插入時反應堆的臨界本征值，kout是控制棒組未插入時反應堆的臨界本征值。

控制棒價值的計算如表5 所示，與BEAVRS 提供的實驗值對比，誤差均在5×10-4以內；如表6所示，隨著更多的控制棒組的插入，反應堆越來越偏離臨界狀態(tài)。通過比較三種均勻化方案的計算結果，三種均勻化方案的計算對控制棒價值的計算精度影響不大，這是因為不同的均勻化方案主要體現(xiàn)在對含有非對稱布置可燃毒物的組件均勻化計算中，這些組件中并不含有控制棒。而含有控制棒的組件都采用方案一進行均勻化計算，因此對于控制棒價值以及臨界硼濃度的計算結果影響不大，但對于考慮空間分布的裂變率計算有顯著的影響。

表5 控制棒價值計算Table 5 Calculation of control rod bank worths

表6 不同控制棒組對應的keff計算值Table 6 Calculation of different control rods for keff

3.4 裂變率

BEAVRS 基準題中共有58 個包含探測器的組件，測量值是將235U裂變率的數(shù)據(jù)歸一化后的結果。圖7 展示了傳統(tǒng)均勻化方法（均勻化方案1）在堆芯徑向探測器的計算值和實驗值的比較，其中，在B13組件處具有16.67%的最大相對誤差。方案1 的計算結果與基準值整體趨勢符合較好，但是右下方誤差較大。在熱態(tài)零功率下，堆芯功率太低，測量值誤差較大，因此出現(xiàn)了徑向功率傾斜現(xiàn)象［13-16］。OpenMC 的計算結果與cosRMC 的計算結果具有相同的傾斜趨勢且在B13 組件同樣存在最大相對誤差［17-18］，由此可見，方案1 的計算結果在合理范圍之內。

圖7 方案1計算熱態(tài)零功率下235U裂變率與實驗值比較Fig.7 Calculation of the fission rate of 235U using scheme 1 compared with the BEAVRS benchmark

圖8展示了方案2和方案3分別與方案1計算結果的對比，誤差一為方案2 與方案1 的相對偏差對比，誤差二為方案3與方案1的相對偏差對比。與采用傳統(tǒng)均勻化方法的方案1相比，方案2和方案3在位于非對稱燃料組件附近的組件上數(shù)值有了明顯的變化，方案1 在非對稱燃料組件附近的組件的平均誤差為5.62%，方案2和方案3的平均誤差為3.345%和3.498%，與方案1 相比，方案2 和方案3 的誤差均有降低。圖9 展現(xiàn)了方案2 和方案3 與基準值的對比，誤差一為方案2與基準值的相對偏差對比，誤差二為方案3 與基準值的相對偏差對比。方案2 和方案3的最大誤差為13.38%和13.76%，仍在B13組件中。相較于方案1，方案2和方案3在堆芯邊緣的計算精度有了很大的提升，這是由于方案2和方案3在均勻化的時候考慮了空間上的非對稱性，因此能夠更加精確的提供裂變率分布；可以看出，考慮組件的非對稱性可以有效的提高計算精度。

圖8 方案2、方案3與方案1的計算結果對比示意圖Fig.8 Comparison diagram of calculation results between Scheme 2,Scheme 3 and Scheme 1

圖9 方案2、方案3與BEAVRS實驗值的對比示意圖Fig.9 Comparison diagram between calculation results of Schemes 2,3 and the BEAVRS benchmark

圖10顯示了方案2與方案3的計算結果之間的差異，不難發(fā)現(xiàn)，兩種方案的最大誤差為1.21%，其余區(qū)域的誤差均為1%以內，說明兩種方案的計算結果差異不大；但相比于方案2，方案3 需要單獨建立4 個1/4 模型，計算成本相較于前者更高，因此方案2 更適合作為非對稱燃料組件均勻化計算的方案。

圖10 方案2與方案3的計算結果對比示意圖Fig.10 Comparison diagram of calculation results between Scheme 2 and Scheme 3

4 結語

本文基于DRAGON/DONJON 對精細化壓水堆基準題BEAVRS 2.02 進行了三維建模計算與分析。同時進一步對于非對稱燃料組件的均勻化方案進行了研究，包括傳統(tǒng)均勻化方案、一次多區(qū)均勻化方案和多次多區(qū)均勻化方案。計算結果表明：對于臨界硼濃度和控制棒價值的計算，臨界硼濃度的誤差均在5×10-5內，控制棒價值的誤差均在5×10-4以內；對于全堆芯的裂變率計算，一次多區(qū)均勻化方案和多次多區(qū)均勻化方案相比于傳統(tǒng)均勻化方案的計算精度顯著提升，對于非對稱組件及相鄰組件的平均誤差能從5.62%分別降低至3.345%和3.498%。通過本文的計算結果表明“組件-堆芯”兩步法對于精細化模型建立的有效性和計算的正確性。后續(xù)將考慮對BEAVRS 2.02基準題中的多循環(huán)算例進行建模計算與對比驗證。

作者貢獻聲明袁顯寶負責建立模型；劉曾豪研究并分析數(shù)據(jù)以及起草論文；張彬航負責設計具體研究內容和方向，提出合理化研究方案以及提供理論支撐，同時指導程序的使用以及數(shù)據(jù)的采集，對論文的知識性內容做審閱并且負責對論文的修改；張永紅負責論文審閱；唐海波負責論文審閱；楊森權提供技術支持；肖云龍?zhí)峁┘夹g支持。