李鑄倫，謝金森,*，徐士坤，鄧年彪，苑旭東，于 濤

(1.南華大學(xué) 核科學(xué)技術(shù)學(xué)院, 湖南 衡陽 421001； 2.南華大學(xué) 湖南省數(shù)字化核反應(yīng)堆工程與技術(shù)研究中心, 湖南 衡陽 421001)

通過反應(yīng)堆物理實(shí)驗(yàn)進(jìn)行核設(shè)計(jì)驗(yàn)證是開發(fā)新型反應(yīng)堆的必要環(huán)節(jié)，其中與反應(yīng)性(或有效增殖因數(shù)keff)相關(guān)的物理量(如控制棒價(jià)值、停堆深度等)的實(shí)驗(yàn)測量及利用其對(duì)核設(shè)計(jì)程序進(jìn)行校核是一項(xiàng)重要的工作內(nèi)容。實(shí)驗(yàn)得到的反應(yīng)性通常是以緩發(fā)中子有效份額βeff為單位的相對(duì)反應(yīng)性(單位為$)，而理論計(jì)算得到的反應(yīng)性是以keff換算得到的絕對(duì)反應(yīng)性，兩者之間的相互驗(yàn)證必須依賴中子動(dòng)力學(xué)參數(shù)作為橋梁。隨著加速器驅(qū)動(dòng)次臨界反應(yīng)堆(ADS)向工程邁進(jìn)[1]，以及核動(dòng)力反應(yīng)堆在次臨界狀態(tài)下實(shí)現(xiàn)控制棒刻度等實(shí)際需求出現(xiàn)，次臨界(特別是較深次臨界)絕對(duì)反應(yīng)性的精確測量正逐步得到重視。

次臨界反應(yīng)性測量方法的一類經(jīng)典方法(脈沖中子源(PNS)法、Rossi-α法等)是基于瞬發(fā)中子基波衰減常數(shù)α的。這類方法可通過式(1)得到絕對(duì)反應(yīng)性ρa(bǔ)bs，其中Λ為中子代時(shí)間。謝金森等[2-3]從中子時(shí)空動(dòng)力學(xué)方程出發(fā)，指出用瞬發(fā)中子基波衰減常數(shù)α獲得絕對(duì)反應(yīng)性時(shí)，中子動(dòng)力學(xué)參數(shù)計(jì)算中應(yīng)采用瞬發(fā)α基波中子注量率。

(1)

蒙特卡羅方法具有良好的幾何適應(yīng)性，采用連續(xù)能量截面數(shù)據(jù)還可使其具備良好的能譜適應(yīng)性，被認(rèn)為是反應(yīng)堆中子輸運(yùn)模擬的一類終極方法。利用蒙特卡羅方法進(jìn)行瞬發(fā)α本征值及瞬發(fā)α基波中子注量率計(jì)算時(shí)，通常將α本征值問題轉(zhuǎn)換為keff本征值問題間接求解，即k-α迭代策略。該策略通過在keff本征值問題中引入時(shí)間吸收項(xiàng)α/v，通過迭代尋找keff=1時(shí)的α值求解瞬發(fā)α本征值問題[4]。Cullen等[5]采用k-α迭代思想，利用TART程序?qū)odiva及其衍生基準(zhǔn)題進(jìn)行了模擬計(jì)算。經(jīng)典的通用蒙特卡羅程序MCNP4C同樣基于k-α迭代求解瞬發(fā)α本征值問題，但它對(duì)α的初始猜測值較為敏感，同時(shí)選取中子平均裂變壽命的倒數(shù)作為迭代參數(shù)，導(dǎo)致在中子速度較低或次臨界度偏深時(shí)誤差較大，甚至計(jì)算無法進(jìn)行[6-7]。文獻(xiàn)[8-10]為更好地進(jìn)行反應(yīng)性引入事故下的多群時(shí)空動(dòng)力學(xué)計(jì)算，在OpenMC上開發(fā)了緩發(fā)α本征值計(jì)算模塊，對(duì)次臨界問題，該模塊強(qiáng)制α本征值逼近負(fù)的最小緩發(fā)中子先驅(qū)核衰減常數(shù)，無法用于次臨界瞬發(fā)α本征值問題計(jì)算。

針對(duì)MCNP4C因k-α迭代參數(shù)及初始α猜想值導(dǎo)致的次臨界度較深時(shí)瞬發(fā)α本征值誤差較大甚至無法計(jì)算的問題，以及文獻(xiàn)[8-10]的工作無法實(shí)現(xiàn)次臨界狀態(tài)瞬發(fā)α本征值計(jì)算的問題，本文基于傳統(tǒng)的k-α迭代方法，首先對(duì)迭代參數(shù)進(jìn)行重新選擇，其次在粒子輸運(yùn)過程中對(duì)瞬發(fā)、緩發(fā)中子分別考慮，基于開源蒙特卡羅程序OpenMC的二次開發(fā)，保留其計(jì)算緩發(fā)α本征值問題的功能，同時(shí)增加單獨(dú)可選的瞬發(fā)α本征值問題計(jì)算功能。利用簡單幾何的Godiva衍生基準(zhǔn)題及復(fù)雜幾何和燃料的MUSE-4次臨界實(shí)驗(yàn)裝置對(duì)OpenMC-PA進(jìn)行驗(yàn)證。

1 OpenMC介紹

OpenMC是由麻省理工學(xué)院(MIT)計(jì)算反應(yīng)堆物理組于2011年開始研發(fā)、并在2012年12月發(fā)行第1版的并行蒙特卡羅計(jì)算程序，其主要目的是開發(fā)一種能根據(jù)研究需要易擴(kuò)展的高性能、免費(fèi)開源的模擬計(jì)算程序。OpenMC輸入文件是一系列的XML格式文件，以便于數(shù)據(jù)在不同程序和接口之間的交換。具體包括幾何結(jié)構(gòu)、材料定義、計(jì)算參數(shù)設(shè)置、計(jì)數(shù)類型、繪圖及加速收斂設(shè)置等輸入文件[11]。

2 瞬發(fā)α本征值問題理論計(jì)算方法

2.1 理論模型

無外源增殖系統(tǒng)的非穩(wěn)態(tài)中子輸運(yùn)問題可由式(2)描述：

(2)

式中：t為時(shí)間；v為中子速度(標(biāo)量)；f(r,E′,Ω′→E,Ω)為散射函數(shù)，表示碰撞前中子的能量為E′，運(yùn)動(dòng)方向?yàn)棣浮?，碰撞后中子能量變?yōu)镋而運(yùn)動(dòng)方向?yàn)棣傅母怕?；φ′與φ分別為碰撞前后的中子注量率，φ′=φ′(r,E′,Ω′,t)、φ=φ(r,E,Ω,t)；r為中子群體的空間坐標(biāo)；Σt為宏觀總截面；Σs為宏觀散射截面；Σf為宏觀裂變截面；χ為中子裂變能譜；ν為每次裂變釋放的平均中子數(shù)；λ為緩發(fā)中子先驅(qū)核有效衰減常數(shù)；Cj為緩發(fā)中子先驅(qū)核密度；下標(biāo)j表示緩發(fā)中子分組號(hào)，p和d分別表示瞬發(fā)中子和緩發(fā)中子。

將中子注量率和緩發(fā)中子先驅(qū)核密度隨時(shí)間的變化規(guī)律寫成如下指數(shù)形式：

Cj(r,t)=Cα,j(r)eαt

(3)

將式(3)代入式(2)可得到廣義α本征值方程[8]：

(4)

由于瞬發(fā)中子與緩發(fā)中子具有明顯不同的時(shí)間特性，因此式(4)中的α本征值可分為兩族：其中一族絕對(duì)值較小，即緩發(fā)α本征值；另一族絕對(duì)值較大，即瞬發(fā)α本征值。對(duì)于次臨界系統(tǒng)，緩發(fā)與瞬發(fā)α本征值均小于0；對(duì)于臨界與緩發(fā)超臨界系統(tǒng)，緩發(fā)α本征值≥0，而瞬發(fā)α本征值≤0；對(duì)于瞬發(fā)超臨界系統(tǒng)，緩發(fā)與瞬發(fā)α本征值均大于0。

對(duì)于次臨界問題，若直接數(shù)值求解式(4)，由于占優(yōu)本征值是緩發(fā)α本征值，因此最終迭代結(jié)果將返回αd，隨著次臨界度的加深，緩發(fā)α本征值將回歸到負(fù)的最小緩發(fā)中子先驅(qū)核衰減常數(shù)，即-min{λj|j=1,2,…,J}，這符合物理上中子注量率的衰減速度不能超過壽命最長的緩發(fā)中子先驅(qū)核組的衰減常數(shù)[12]。

(5)

2.2 k-α迭代基本思想

瞬發(fā)α本征值的確定論方法迭代方案描述如下。

(6)

(7)

(8)

(9)

式中，l為迭代的代數(shù)索引。式(6)即k-α迭代的工作原理，裂變項(xiàng)滯后表示抽樣的裂變中子被存儲(chǔ)并用于下一次迭代時(shí)的固定源，歸一化后如式(9)所示。標(biāo)準(zhǔn)的k本征值模擬計(jì)算程序原理如式(7)所示。式(8)所示的瞬發(fā)α本征值更新迭代過程是通過使用新的特征函數(shù)估計(jì)在整個(gè)相空間上對(duì)式(5)進(jìn)行積分，使用歸一化條件和本征值的定義來替代吸收和泄漏的積分[8]。

在式(6)中，對(duì)于任意的α都有唯一的k與之對(duì)應(yīng)，即k=k(α)。將系統(tǒng)調(diào)整至“虛擬臨界”的過程，即是尋找一個(gè)α值，可以使得k=k(α)=1，也即求解非線性方程k(α)-1=0。求解該非線性方程的方法一般是用迭代法，比較常用的迭代格式為：

αl+1=αl+cl(kl-1)

(10)

式中，cl為迭代參數(shù)，單位為μs-1，會(huì)隨著迭代的進(jìn)行而改變。迭代參數(shù)的選取不僅會(huì)影響收斂速度，更會(huì)影響誤差估計(jì)[6]。在文獻(xiàn)[15]中，Brockway等建議cl=αl+gn，gn為中子速度與宏觀總截面乘積的期望值；而MCNP4C使用的迭代參數(shù)cl=αl+1/lf，lf為中子平均裂變壽命[16]。文獻(xiàn)[17]對(duì)3種迭代參數(shù)(cl=αl+1/lf、cl=1/lf、cl=400)進(jìn)行了比較，結(jié)果表明偏離臨界越遠(yuǎn)，MCNP4C給出的統(tǒng)計(jì)誤差越大，反映了其迭代參數(shù)的選取并不完全可靠。

在OpenMC-PA中選擇的迭代參數(shù)計(jì)算為：

(11)

式中：d為中子連續(xù)兩次反應(yīng)間穿行的距離；s為積分變量，ds表示穿行距離上的微分量；DTl為中子在平均徑跡長度上的時(shí)間吸收權(quán)重[18]，它反映了中子隨機(jī)游走過程中在徑跡長度上的衰減過程，更符合真實(shí)的物理過程。

3 基于OpenMC的程序?qū)崿F(xiàn)

根據(jù)確定論方法的迭代思想，在中子輸運(yùn)方程中引入k作為參數(shù)計(jì)算瞬發(fā)α本征值，當(dāng)k=1時(shí)，式(6)與式(5)完全一樣，則k-α迭代的思想即是加入時(shí)間吸收截面[α/v]，將α本征值問題轉(zhuǎn)化為k本征值問題的反問題，即通過調(diào)節(jié)α本征值對(duì)k本征值進(jìn)行迭代計(jì)算，直至k本征值收斂于1，此時(shí)求得的α本征值即是系統(tǒng)瞬發(fā)中子衰減常數(shù)。

根據(jù)以上計(jì)算思想，OpenMC-PA總體通過兩層迭代計(jì)算：內(nèi)迭代完成k本征值的迭代計(jì)算，外迭代完成瞬發(fā)α本征值的迭代計(jì)算。其中，瞬發(fā)α本征值計(jì)算模塊采用牛頓-拉夫遜迭代算法判定收斂條件。計(jì)算流程圖如圖1所示，其計(jì)算基本流程如下。

圖1 基于OpenMC的α本征值計(jì)算模塊流程Fig.1 Flow of α eigenvalue calculation module based on OpenMC

1) 開始內(nèi)迭代，給出初始估計(jì)值kTl=1、初始值αeff=0、收斂判定誤差ε=1×10-8、迭代的相對(duì)誤差與標(biāo)準(zhǔn)殘差均為1，確定粒子的初始能量、初始位置與初始運(yùn)動(dòng)方向，l=0時(shí)進(jìn)行第1次迭代，總共計(jì)算L代。

2) 采用系統(tǒng)徑跡長度估計(jì)模式統(tǒng)計(jì)粒子移動(dòng)軌跡預(yù)估子代粒子總數(shù)Nl。

3) 確定下一個(gè)碰撞點(diǎn)位置并判斷粒子是否穿出材料，穿出材料則粒子死亡停止追蹤該粒子。

4) 根據(jù)截面確定碰撞對(duì)象。

5) 選取[0,1)范圍內(nèi)的隨機(jī)數(shù)ξ對(duì)粒子抽樣判斷是否發(fā)生[α/v]反應(yīng)。發(fā)生[α/v]反應(yīng)且α>0則判定時(shí)間吸收，粒子死亡；發(fā)生[α/v]反應(yīng)且α<0則判定δ散射，粒子能量與方向不變，權(quán)重變?yōu)樵瓉淼?倍；不發(fā)生[α/v]反應(yīng)，則按照α=0進(jìn)行碰撞抽樣。

7) 子代粒子數(shù)與父代粒子數(shù)相除得到徑跡長度估計(jì)有效增殖因數(shù)kTl。

8) 根據(jù)瞬發(fā)α的初始估計(jì)值，通過式(6)進(jìn)行外迭代計(jì)算新的瞬發(fā)α本征值。

9) 通過牛頓-拉夫遜迭代算法收斂計(jì)算瞬發(fā)α本征值，更新迭代的相對(duì)誤差與標(biāo)準(zhǔn)殘差，并作為下一代的α初始估計(jì)值代入下一輪計(jì)算。對(duì)中子注量率進(jìn)行更新。

10) 將初始αeff的值替換更新為迭代計(jì)算的α本征值，輸出當(dāng)前第l代的α本征值(αeff)。

11) 若當(dāng)前已完成預(yù)設(shè)最大迭代代數(shù)L，則程序結(jié)束。

4 基準(zhǔn)題問題檢驗(yàn)

選取Godiva衍生基準(zhǔn)題[5]以及MUSE-4基準(zhǔn)題[19]分別對(duì)OpenMC-PA進(jìn)行超臨界系統(tǒng)以及次臨界系統(tǒng)的模擬計(jì)算。MUSE-4基準(zhǔn)題模型的材料及參數(shù)參見文獻(xiàn)[19]。

4.1 Godiva衍生基準(zhǔn)題

Godiva1模型是一個(gè)高富集度的純鈾裸球系統(tǒng)，Godiva2與Godiva4模型是在Godiva1的基礎(chǔ)上更改的均勻超臨界裸球系統(tǒng)。Godiva2是密度翻倍的快中子超臨界裸球系統(tǒng)，Godiva4是原子密度之比為100∶1的水與鈾直接混合的熱中子超臨界裸球系統(tǒng)，模型及參數(shù)列于表1。Godiva3模型是在Godiva1的基礎(chǔ)上更改的非均勻快中子超臨界系統(tǒng)，在燃料的外側(cè)增加了30 cm的水作為反射層。Godiva3模型及參數(shù)列于表2。

表1 Godiva1、Godiva2、Godiva4參數(shù)Table 1 Parameter of Godiva1, Godiva2 and Godiva4

表2 Godiva3參數(shù)Table 2 Parameter of Godiva3

分別以50代非活躍中子、300代活躍中子、每代106個(gè)中子進(jìn)行計(jì)算，結(jié)果列于表3。由表3可見，OpenMC計(jì)算結(jié)果與TART的計(jì)算結(jié)果吻合較好，最大相對(duì)誤差小于0.5%，以簡單模型驗(yàn)證了超臨界問題中本文程序計(jì)算結(jié)果的準(zhǔn)確性。

以Godiva1基準(zhǔn)題模型為參考，通過持續(xù)降低質(zhì)量密度使其從臨界狀態(tài)轉(zhuǎn)為次臨界狀態(tài)，隨著次臨界度的加深，分別以50代非活躍中子、300代活躍中子、每代104個(gè)中子，在OpenMC與MCNP4C中計(jì)算該次臨界問題。計(jì)算結(jié)果列于表4。

表3 Godiva衍生基準(zhǔn)題計(jì)算結(jié)果對(duì)比Table 3 Comparison of calculation result of Godiva derivative benchmark

表4 Godiva基準(zhǔn)題次臨界計(jì)算結(jié)果對(duì)比Table 4 Comparison of subcritical calculation result of Godiva benchmark

4.2 MUSE-4次臨界實(shí)驗(yàn)裝置

MASURCA反應(yīng)堆是CEA在法國卡達(dá)拉赫研究中心負(fù)責(zé)運(yùn)營的零功率裝置。MUSE-4是基于MASURCA改造的次臨界實(shí)驗(yàn)裝置，堆芯模型如圖2所示，不同類型MUSE-4組件的軸向視圖如圖3所示。該裝置在976組燃料組件裝載時(shí)的有效增殖因數(shù)keff=0.97[19]，緩發(fā)中子有效份額βeff=0.003 4，中子代時(shí)間Λ=0.74×10-6s。MUSE-4實(shí)驗(yàn)裝置參數(shù)列于表5。

圖2 MUSE-4的柵元裝載形式Fig.2 Cell loading configuration of MUSE-4

圖3 MUSE-4組件1～8和11的軸向視圖Fig.3 Axial view of MUSE-4 tube 1-8 and 11

對(duì)MUSE-4實(shí)驗(yàn)裝置模型在每代10 000粒子、100代非活躍代、1 000代活躍代下進(jìn)行模擬計(jì)算，結(jié)果列于表6。

表5 MUSE-4實(shí)驗(yàn)裝置參數(shù)Table 5 Parameter of MUSE-4 experimental device

表6 MUSE-4的模擬結(jié)果與基準(zhǔn)題對(duì)比驗(yàn)證Table 6 Comparison and verification of MUSE-4 simulation result and benchmark result

由表6可看出，OpenMC得到的計(jì)算結(jié)果與基準(zhǔn)題相比相對(duì)誤差均在0.5%以內(nèi)，驗(yàn)證了OpenMC在該模型下的計(jì)算精度滿足工程設(shè)計(jì)需求。同時(shí)，由OpenMC-PA得到的α，利用式(1)回歸計(jì)算得到的keff為0.970 1±0.000 1，與實(shí)驗(yàn)值吻合良好，相對(duì)誤差小于0.02%，從另一個(gè)角度驗(yàn)證了OpenMC-PA計(jì)算瞬發(fā)α本征值的正確性。

5 結(jié)論

本文采用k-α迭代方法，選取徑跡長度上裂變中子權(quán)重的平均修正參數(shù)DTl作為迭代參數(shù)，在粒子輸運(yùn)過程中對(duì)瞬發(fā)、緩發(fā)中子分別考慮，在開源蒙特卡羅程序OpenMC上開發(fā)了OpenMC-PA模塊，增加了單獨(dú)可選的瞬發(fā)α本征值問題計(jì)算功能?；贕odiva衍生基準(zhǔn)題以及MUSE-4次臨界實(shí)驗(yàn)裝置，分別對(duì)超臨界及次臨界狀態(tài)下的模擬計(jì)算結(jié)果進(jìn)行了驗(yàn)證。結(jié)果表明，在OpenMC上模擬計(jì)算得到的瞬發(fā)α本征值結(jié)果與基準(zhǔn)題結(jié)果吻合較好，在一定次臨界度下可保持穩(wěn)定運(yùn)行，瞬發(fā)α本征值計(jì)算結(jié)果的相對(duì)誤差穩(wěn)定在0.5%以內(nèi)，計(jì)算范圍與計(jì)算時(shí)間均優(yōu)于MCNP4C，驗(yàn)證了開發(fā)的OpenMC-PA模塊的可靠性和準(zhǔn)確性。后續(xù)可基于OpenMC繼續(xù)開展改進(jìn)的k-α迭代等方法優(yōu)化瞬發(fā)α本征值的計(jì)算。

感謝密歇根大學(xué)Ilham Variansyah博士對(duì)本論文研究提供的幫助。