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

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

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

在數(shù)值反應(yīng)堆高保真建模模擬中，一般采用兩類高保真計(jì)算程序，分別為蒙特卡羅程序和確定論程序。蒙特卡羅程序主要以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)堆物理團(tuán)隊(duì)開發(fā)的確定論核反應(yīng)堆物理計(jì)算程序［7］。OpenMC是由美國麻省理工學(xué)院（Massachusetts Institute of Technology）的計(jì)算反應(yīng)堆物理小組開發(fā)的程序［8］，該程序可進(jìn)行臨界問題和多種固定源問題的計(jì)算，也包含了多群功能，主要應(yīng)用于反應(yīng)堆物理分析設(shè)計(jì)；可產(chǎn)生多群均勻化截面，也可以進(jìn)行緩發(fā)中子的產(chǎn)生計(jì)數(shù)。

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

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

1 程序介紹

1.1 NECP-X程序

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

1.2 OpenMC程序

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

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

2 基準(zhǔn)題描述

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

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

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

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

3 計(jì)算程序建模

3.1 NECP-X程序建模

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

圖1 模型1的堆芯布置（1/4對(duì)稱）(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對(duì)稱）(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 蒙特卡羅程序

在計(jì)算模型中，設(shè)置每一代中子數(shù)目為2×107個(gè)，總中子代數(shù)為800代，跳過非活躍代數(shù)為100代，同時(shí)測量統(tǒng)計(jì)每個(gè)組件的總裂變能沉積能量。OpenMC程序的計(jì)算采用64核并行，在上海交通大學(xué)JCLOUD（Jiaotong Cloud）云平臺(tái)上開展。以模型2為例，計(jì)算不同控制棒組插入情況下的模型2所需時(shí)間約為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在相同模型條件下對(duì)模型1進(jìn)行了建模計(jì)算，計(jì)算了熱態(tài)零功率（Hot Zero Power，HZP）狀態(tài)下的有效增殖因子keff。在HZP狀態(tài)下，反應(yīng)堆正常運(yùn)行時(shí)沒有控制棒插入堆芯，硼濃度為9.75×10-4。兩程序截面數(shù)據(jù)庫均采用ENDF/B-VII庫。有效增殖因子keff計(jì)算結(jié)果如表1所示，NECP-X與OpenMC計(jì)算結(jié)果符合較好，相差為1.37×10-3。

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

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

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

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

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

4.2 模型2

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)語

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

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