基于BEAVRS基準題高保真建模的OpenMC程序和NECP-X程序的對比驗證

沈芷睿孫啟政何東豪潘清泉張滕飛彭良輝楊偉焱

1（上海交通大學(xué)核能與核技術(shù)工程上海200240）

2（上海核工程研究設(shè)計院有限公司上海200233）

在數(shù)值反應(yīng)堆高保真建模模擬中，一般采用兩類高保真計算程序，分別為蒙特卡羅程序和確定論程序。蒙特卡羅程序主要以MCNP（Monte Carlo NParticle Transport Code）［1］、JMCT（J Monte Carlo Code）［2］、RMC（Reactor Monte Carlo code）和SuperMC（Super Monte Carlo）［3］等為代表；確定論程序主要以nTRACER［4］、Decart［5］和MPACT［6］等為代表。NECP-X是西安交通大學(xué)反應(yīng)堆物理團隊開發(fā)的確定論核反應(yīng)堆物理計算程序［7］。OpenMC是由美國麻省理工學(xué)院（Massachusetts Institute of Technology）的計算反應(yīng)堆物理小組開發(fā)的程序［8］，該程序可進行臨界問題和多種固定源問題的計算，也包含了多群功能，主要應(yīng)用于反應(yīng)堆物理分析設(shè)計；可產(chǎn)生多群均勻化截面，也可以進行緩發(fā)中子的產(chǎn)生計數(shù)。

反應(yīng)堆模擬評價和驗證基準題（BEAVRS benchmark）提供了典型的美國商用壓水堆堆芯的初始運行周期1和重新裝載的運行周期2的實測物理數(shù)據(jù)。該基準題不僅包含了詳細的堆芯參數(shù)、燃料組成和其他重要組成的參數(shù)，還提供了熱態(tài)零功率條件下和運行時堆芯的詳細結(jié)果和參數(shù)。

蒙特卡羅程序和MOC（Method of Characteristics）方法是高保真計算的代表，現(xiàn)階段就提高其效率問題展開了許多研究。例如將CAD（Computer-aided Design）轉(zhuǎn)換為蒙特卡羅幾何的轉(zhuǎn)換程序CMGC（The CAD to Monte Carlo geometry converter）可以準確地形成CSG（Constructive Solid Geometry），從而提高MCNP程序計算效率［9］，同時在SuperMC程序中采用線程級數(shù)據(jù)分解方法可滿足計算全堆芯高保真輸運-燃耗計算的存儲，又能保持較高的并行效率［10］。MOC方法常被用于確定論程序中進行輸運計算，采用低階角通量非線性差分方程應(yīng)用于MOC方法可提高計算速度［11］。提高高保真模擬計算的同時，也要保證高保真計算的精度，為驗證確定論一步法和概率論方法在高保真計算方面的精度水平。本研究利用確定論程序NECP-X和蒙特卡羅程序OpenMC，分別基于BEAVRS基準題進行高保真建模，通過對基準題模型的部分物理參數(shù)進行計算，初步驗證了兩個程序?qū)?fù)雜堆芯精細建模計算的可行性和準確性。同時，本研究還將比較兩個程序?qū)ν荒Ｐ偷南嗤锢韰?shù)的計算結(jié)果，分析偏差來源以及兩個程序的計算效率，為程序以后的應(yīng)用及完善提供參考。

1 程序介紹

1.1 NECP-X程序

NECP-X能夠?qū)θS全堆芯問題直接進行共振和輸運計算。該程序的數(shù)據(jù)庫被分為兩部分：多群數(shù)據(jù)庫NECL-MG和連續(xù)能量數(shù)據(jù)庫NECL-CE。共振計算中使用了偽共振核素方法，通過共振自屏蔽計算得到有效的自屏蔽截面。偽共振核素是每個燃料棒中所有共振核素的混合物，利用偽共振核素和中間共振近似得到子群概率，并使用丹可夫因子考慮空間自屏效應(yīng)得到每個燃料棒的等效一維圓柱形問題，將全堆芯問題分解成一組一維問題［12］。NECP-X的全堆芯非均勻輸運計算采用了二維/一維耦合輸運計算方法，計算二維問題時采用二維特征線法，計算一維問題采用SN方法。同時，還采用了多級并行策略，包括空間區(qū)域分解、角度區(qū)域分解和特征線并行提高計算效率。NECP-X還采用了MOC方法中常用的粗網(wǎng)有限差分法以加速輸運方程中外迭代的收斂速度［13］。

1.2 OpenMC程序

OpenMC是一款三維的連續(xù)能量蒙特卡羅反應(yīng)堆物理開源程序［14］。該程序能夠在利用立體幾何構(gòu)造或由CAD構(gòu)造的三維模型中進行固定源、k-本征值問題等不同的類型的中子輸運計算，同時兼具燃耗計算的功能。

OpenMC內(nèi)置了基于構(gòu)造性實體幾何表示（Constructive Solid Geometry representation）的模型構(gòu)造方法，利用不同的曲面構(gòu)造出封閉的空間，即柵元（cell），能夠?qū)^大部分幾何體進行準確地構(gòu)建，避免了由于空間網(wǎng)格離散引入的誤差。OpenMC還采用了fission bank算法，減少了不必要的并行通信，能夠進行更高性能的并行計算［15-16］。

2 基準題描述

BEAVRS基準題來源于美國西屋公司的一個真實壓水堆，該基準題詳細描述了堆芯的幾何結(jié)構(gòu)以及堆芯的運行參數(shù)。

本文對基準題中描述的模型進行了適當(dāng)?shù)睾喕７磻?yīng)堆全高為435.444 cm，包括堆芯網(wǎng)格組件、堆芯圍板以及反射層等。每個組件均勻分成17×17個柵元，燃料棒根據(jù)235U富集度的不同分為1.6%、2.4%和3.1%三種類型。根據(jù)組件在堆芯中的具體位置，每個組件的可燃毒物的排列方式以及個數(shù)不同，根據(jù)數(shù)量可燃毒物可分為12個和16個兩種?？刂瓢舨牧蠟锳g-In-Cd［17］。

本文主要對兩種堆芯模型進行了高保真建模，模型1為5×5全堆芯模型，其中共含有21個燃料組件，在模型1建模時考慮了軸向支撐格架、定位格架和可燃毒物；模型2為15×15全堆芯模型，其中共含有193個燃料組件，在模型2建模時省略了軸向支撐格架、定位格架和可燃毒物。

利用較為簡單的模型1，先驗證NECP-X和OpenMC兩個程序的建模以及兩個程序計算設(shè)置的正確性，同時也驗證BEAVRS基準題模型進行簡化的可行性。此后，利用更接近于真實堆芯情況的模型2，驗證兩個程序?qū)φ鎸嵉膹?fù)雜堆芯精細建模計算的可行性和準確性。

3 計算程序建模

3.1 NECP-X程序建模

在模擬計算中，采用1/4堆芯對稱模型以降低計算量。三維算例計算模型：對于UO2非均勻柵格，對燃料區(qū)域徑向劃分5圈，幅角方向劃分8區(qū)；對于慢化劑區(qū)域徑向劃分3圈，幅角方向劃分8區(qū)。軸向采用SN求解器，徑向采用MOC求解器，MOC特征線線寬為0.025 cm，求積數(shù)組采用Chebyshev_Yamamoto求積組，每個卦限內(nèi)的幅角數(shù)與極角數(shù)分別為8和3。特征值與注量率的收斂準則分別為1.0×10-5和1.0×10-4。輸 運 修 正 采 用inflow方 法，利 用global_local共振計算方法進行共振自屏計算，共振干涉因子方法為pmm，采用能群SPH（Smoother Particle Hydrodynamics）修 正［18］。一 步 法 程 序NECP-X采用64能群結(jié)構(gòu)，軸向采用2個區(qū)并行運算，徑向采用257個區(qū)并行運算，程序在浪潮機架服務(wù)器NF5280M5運行，服務(wù)器共用39個運算節(jié)點，每個節(jié)點有40個CPU。以模型2為例，計算不同控制棒組插入情況下的模型2所需時間約為0.72 h。模型1的堆芯組件布置如圖1所示，模型2的堆芯組件布置如圖2所示。

圖1 模型1的堆芯布置（1/4對稱）(a)堆芯燃料裝載方式，(b)毒藥裝載方式Fig.1 Core assembly layout of model 1(quarter symmetry)(a)Core fuel loading pattern,(b)Poison loading pattern

圖2 模型2的堆芯布置（1/4對稱）(a)堆芯燃料裝載方式，(b)堆芯控制棒布局Fig.2 Core assembly layout of model 2(quarter symmetry)(a)Core fuel loading pattern,(b)Layout of core control rod

3.2 蒙特卡羅程序

在計算模型中，設(shè)置每一代中子數(shù)目為2×107個，總中子代數(shù)為800代，跳過非活躍代數(shù)為100代，同時測量統(tǒng)計每個組件的總裂變能沉積能量。OpenMC程序的計算采用64核并行，在上海交通大學(xué)JCLOUD（Jiaotong Cloud）云平臺上開展。以模型2為例，計算不同控制棒組插入情況下的模型2所需時間約為21.1 h。模型1的堆芯組件布置如圖3所示，模型2的堆芯組件布置如圖4所示。

圖3 模型1的堆芯布置(a)軸向幾何，(b)徑向幾何Fig.3 Core assembly layout of model 1(a)Axial geometry,(b)Radial geometry

圖4 模型2的堆芯布置(a)軸向幾何，(b)徑向幾何Fig.4 Core assembly layout of model 2(a)Axial geometry,(b)Radial geometry

4 模擬結(jié)果分析

4.1 模型1

利用確定論高保真程序NECP-X和蒙特卡羅程序OpenMC在相同模型條件下對模型1進行了建模計算，計算了熱態(tài)零功率（Hot Zero Power，HZP）狀態(tài)下的有效增殖因子keff。在HZP狀態(tài)下，反應(yīng)堆正常運行時沒有控制棒插入堆芯，硼濃度為9.75×10-4。兩程序截面數(shù)據(jù)庫均采用ENDF/B-VII庫。有效增殖因子keff計算結(jié)果如表1所示，NECP-X與OpenMC計算結(jié)果符合較好，相差為1.37×10-3。

表1 熱態(tài)零功率下NECP-X與OpenMC有效增殖因子keff結(jié)果對比Table 1 Comparison of keff calculated by NECP-X and OpenMC in HZP condition

同時，NECP-X和OpenMC分別計算了HZP狀態(tài)下的模型1的組件功率分布，HZP狀態(tài)下的硼濃度相同，均為9.75×10-4。OpenMC預(yù)估中子通量的碰撞來計算用戶自定義的空間和能量的總反應(yīng)率，通過總裂變率乘以每次裂變釋放的能量得到能量沉積。NECP-X根據(jù)不同材料進行分區(qū)，計算得到不同材料區(qū)的每次裂變生成的能量，最后乘以通量與裂變截面得到功率密度。選取模型1中燃料富集度為1.6%的有17×17個柵元的組件，對1/4組件模型的棒相對功率進行了計算，1/4組件模型的棒相對功率分布如圖5所示，相對功率為零的地方插入的是導(dǎo)向管。

圖5 1/4組件棒相對功率分布Fig.5 Relative power distribution of rod(1/4 assembly)

兩個程序的組件功率計算值進行了相同的歸一化處理，兩個程序計算所得的組件功率分布趨勢相同，均呈現(xiàn)出內(nèi)部組件功率較高外部組件功率較低的趨勢。同時，還以NECP-X程序計算值為參考，計算并統(tǒng)計了兩個程序的組件功率計算值的相對偏差和平均絕對偏差。組件功率的相對偏差分布情況如圖6所示，相對偏差的最正偏差為0.79%，最負偏差為-0.80%，平均絕對偏差為0.55%。NECP-X與OpenMC的組件功率計算結(jié)果符合較好，組件功率相對偏差較大的組件多為堆芯的外圍組件，其功率水平較低。在模型1中，兩個程序關(guān)于HZP狀態(tài)下的有效增殖因子keff和組件功率的計算結(jié)果符合較好，誤差在可接受范圍內(nèi)，驗證了模型和兩個程序的計算參數(shù)設(shè)置的正確性。

圖6 組件功率的相對偏差分布Fig.6 Relative deviation distribution of assembly power

4.2 模型2

模型1的驗證結(jié)果，為模型2建立模型和設(shè)置計算參數(shù)奠定了基礎(chǔ)。由于模型2更為復(fù)雜，因此適當(dāng)增加了每代粒子數(shù)目以及總代數(shù)。

利用確定論程序NECP-X和蒙特卡羅程序OpenMC在相同模型條件下對模型2進行了建模計算，計算了HZP狀態(tài)下的有效增殖因子keff。在HZP狀態(tài)下，反應(yīng)堆正常運行時僅有控制棒組件D在213步，硼濃度為9.75×10-4。兩程序截面數(shù)據(jù)庫均采用ENDF/B-VII庫。有效增殖因子keff計算結(jié)果如表2所示，NECP-X與OpenMC的keff計算結(jié)果符合較好，相差為4.5×10-4。

表2 熱態(tài)零功率下NECP-X與OpenMC有效增殖因子keff結(jié)果對比Table 2 Comparison of keff calculated by NECP-X and OpenMC in HZP condition

NECP-X和OpenMC還計算了7組控制棒組在分別完全插入狀態(tài)下的有效增殖因子keff，結(jié)果如表3所示。

表3 插入不同控制棒組件NECP-X和OpenMC keff計算結(jié)果Table 3 Results of control rod bank worths

硼濃度與HZP狀態(tài)下的硼濃度相同，均為9.75×10-4，唯一變量為不同控制棒組插入［19］。有效增殖因子keff隨控制棒組插入而減小，反應(yīng)性也隨之降低。不同控制棒插入時keff減小程度不同，不同的控制棒組對反應(yīng)堆反應(yīng)性影響的大小不同。對兩個程序的keff計算結(jié)果進行比較，NECP-X的計算結(jié)果相對于OpenMC的計算結(jié)果均偏小，其中計算結(jié)果最大偏差為5.9×10-4和6.0×10-4，說明兩個程序的keff計算結(jié)果符合較好。表4給出了對應(yīng)的7組控制棒價值的計算結(jié)果，對反應(yīng)堆反應(yīng)性影響最大的控制棒組均為控制棒組D，比較兩個程序的控制棒價值計算結(jié)果，OpenMC對于控制棒組C的棒價值計算結(jié)果相對于NECP-X的計算結(jié)果較大，且該計算結(jié)果的偏差遠遠大于其他幾組控制棒組棒價值的偏差。兩個程序的計算結(jié)果中控制棒價值最大偏差為4.9×10-4，不超過5.0×10-4。NECP-X與OpenMC的控制棒價值計算結(jié)果符合較好。

表4 控制棒價值結(jié)果Table 4 Results of control rod bank worths

同時，NECP-X和OpenMC還計算了該7組控制棒組分別在完全插入下的組件功率分布。對兩個程序的功率計算值進行了相同的歸一化處理。本文以NECP-X計算值為參考，計算統(tǒng)計了兩個程序的組件功率的相對偏差和平均絕對偏差。圖7列出了7組控制棒插入情況下兩個程序組件功率計算值的相對偏差的分布。

圖7 組件功率的相對偏差分布Fig.7 Relative deviation distribution of fuel assembly power

在控制棒全提工況下，兩組程序計算得到的組件功率分布均呈現(xiàn)內(nèi)部和外部組件功率較低，與BEAVRS基準題所給出的組件功率分布趨勢符合較好。與控制棒全提工況下的組件功率分布相比，插入控制棒組的組件相對功率明顯減少且為最低，同時其周圍組件的相對功率也明顯減少，而其他組件的相對功率顯著增加。NECP-X程序計算得到的組件功率分布嚴格對稱，而OpenMC程序由于計算結(jié)果本身存在標(biāo)準差，組件功率分布不完全對稱。

不同控制棒組插入工況下的組件功率相對偏差分布中，相對偏差最大為1.62%，出現(xiàn)在控制棒組SE組全插入的工況下；相對偏差最小為-1.65%，出現(xiàn)在控制棒組SB組全插入的工況下。

組件功率分布相對偏差并不嚴格對稱，OpenMC計算得到的組件功率結(jié)果存在標(biāo)準差，自身的計算過程中會引入一定的偏差。同時，蒙特卡羅計算過程中所設(shè)置的每代粒子數(shù)和總代數(shù)也會影響計算結(jié)果，增加每代粒子數(shù)和總代數(shù)可減少計算結(jié)果的標(biāo)準差但同時也會增加計算所需要的時間。同時，當(dāng)堆芯功率很低時也會發(fā)生功率傾斜現(xiàn)象。

除此現(xiàn)象之外，組件功率相對偏差較大的組件多為堆芯的外圍組件，即富集度為3.1%的燃料組件，這與模型1中組件功率相對偏差的分布情況相似。外圍組件靠近反射層，NECP-X模擬計算時，靠近反射層部分的組件相對于內(nèi)部的組件計算和模擬更為復(fù)雜，但因為外圍組件的絕對組件功率值較低，所以對keff的影響較小。

為了更好地分析組件功率分布相對偏差分布現(xiàn)象的原因，本研究同時對組件功率相對偏差的最正偏差、最負偏差以及平均絕對偏差進行計算統(tǒng)計，結(jié)果如表5所示。

表5 組件相對功率計算值偏差(%)Table 5 Results of full assembly power distribution deviation(%)

組件功率平均絕對偏差在0.60%以內(nèi)，插入控制棒組D組時平均偏差最大為0.60%，插入控制棒組B組和SA組時平均絕對偏差最小為0.44%。雖然部分工況下部分組件的相對功率偏差超過1%，但各個工況下的組件相對功率平均絕對偏差均小于0.60%，說明NECP-X與OpenMC的組件功率計算結(jié)果符合較好。

5 結(jié)語

本文使用了NECP-X和OpenMC程序?qū)Υ笮蛪核鸦鶞暑}BEAVRS進行了高保真建模模擬計算。建模對象包括HZP狀態(tài)下的有效增殖因子、7組控制棒的控制棒價值以及不同控制棒組件分別插入情況下的組件功率分布。計算結(jié)果表明：兩個程序的計算結(jié)果偏差在可接受范圍內(nèi)，有效增值因子最大偏差為1.37×10-3，控制棒價值最大偏差為4.9×10-4，組件功率評價偏差在0.60%以內(nèi)。通過對BEAVRS基準題的建模計算，初步驗證了NECP-X和OpenMC在高保真建模、臨界計算、功率分布預(yù)測等方面的準確性與可靠性，同時也為兩個高保真程序的計算效率和實際反應(yīng)堆堆芯模擬方面提供了借鑒和參考。后續(xù)將對BEAVRS基準題的多循環(huán)算例進行高保真模擬計算和對比驗證。

作者貢獻聲明沈芷睿直接參與論文研究，實施研究，采集、分析并解釋計算數(shù)據(jù)，并負責(zé)撰寫論文；孫啟政協(xié)助使用程序；何東豪指導(dǎo)程序的使用以及數(shù)據(jù)的采集，對論文的知識性內(nèi)容做審閱；潘清泉負責(zé)論文的修改；張滕飛負責(zé)論文的修改，對論文的知識性內(nèi)容做審閱；彭良輝提供技術(shù)支持；楊偉焱提供技術(shù)支持。