李雷鳴 江書益 陳珍平 苑旭東 楊 超 謝金森 于 濤

1（南華大學核科學技術(shù)學院 衡陽421001）

2（湖南省數(shù)字化反應(yīng)堆工程技術(shù)研究中心 衡陽421001）

1 結(jié)構(gòu)材料活化分析方法

在反應(yīng)堆運行過程中，堆內(nèi)結(jié)構(gòu)材料受到中子的輻照而發(fā)生活化反應(yīng)，使其產(chǎn)生大量放射性核素，是反應(yīng)堆屏蔽設(shè)計、檢修換料方案、退役策略中都需要考慮的重要內(nèi)容［1］。由于反應(yīng)堆技術(shù)的發(fā)展，各式的先進反應(yīng)堆具有更加復雜的結(jié)構(gòu)以及材料布置與中子能譜，導致中子活化分析也更加復雜，所以有必要開展核反應(yīng)堆結(jié)構(gòu)材料精細活化分析方法研究［2］。傳統(tǒng)的柵元活化方法將每個柵元的中子注量率取平均值進行活化，無法解決由于柵元內(nèi)部結(jié)構(gòu)材料活化的非均勻空間效應(yīng)而引起的活化源項及活化γ源空間分布不夠精確的問題，進而無法對活化γ輻射場進行精確評價［3］，對比見圖1。

圖1 核反應(yīng)堆結(jié)構(gòu)材料活化分析方法對比Fig.1 Comparison of structural material activation analysis methods

為了提高核反應(yīng)堆結(jié)構(gòu)材料活化分析精度，本文發(fā)展了基于網(wǎng)格活化的核反應(yīng)堆結(jié)構(gòu)材料精細活化分析方法。所開發(fā)程序在支持直角網(wǎng)格與柱坐標網(wǎng)格的基礎(chǔ)上支持多個非均勻網(wǎng)格拼接，使網(wǎng)格劃分更加靈活的同時也提升了網(wǎng)格的上限。本文基于嚴格兩步法（Rigorous 2-step method，R2S），基于蒙特卡羅核粒子輸運程序系統(tǒng)（Monte Carlo NParticle Transport Code System，MCNP）［4］與活化程序FISPACT［5］耦合，開發(fā)了由“中子輸運”到“材料活化”再到“活化劑量”進行耦合的反應(yīng)堆結(jié)構(gòu)材料精細化活化計算分析程序（MCNP-FISPACT coupled Mesh-based Activation code，MCFisMA）。

2 基于R2S的網(wǎng)格活化分析方法

2.1 嚴格兩步法

本文基于嚴格兩步法建立核反應(yīng)堆結(jié)構(gòu)材料網(wǎng)格精細活化分析，其主要包括三步：中子輸運網(wǎng)格計算、網(wǎng)格活化源項計算和衰變光子劑量計算，流程如圖2所示。

圖2 基于R2S的網(wǎng)格活化分析方法Fig.2 Mesh-based activation analysis method based on rigorous 2-step method

1）中子輸運網(wǎng)格計算：建立核反應(yīng)堆結(jié)構(gòu)材料蒙特卡羅中子輸運計算模型，對結(jié)構(gòu)材料區(qū)域構(gòu)建網(wǎng)格計算模型，確定中子能群結(jié)構(gòu)，基于MCNP 程序開展蒙特卡羅中子輸運網(wǎng)格計算，得到結(jié)構(gòu)材料區(qū)域的網(wǎng)格中子注量率和網(wǎng)格中子能譜。

2）網(wǎng)格活化源項計算：建立針對結(jié)構(gòu)材料網(wǎng)格模型的網(wǎng)格活化計算模型，根據(jù)中子輸運計算得到的網(wǎng)格中子注量率和網(wǎng)格中子能譜，基于FISPACT開展網(wǎng)格活化計算，得到網(wǎng)格活化源項精細空間分布，主要包括網(wǎng)格主要活化核素及貢獻程度、活化光子源強、活化光子能譜等。

3）活化光子劑量計算：根據(jù)2）中得到的網(wǎng)格活化衰變光子源分布，基于MCNP 進行二次開發(fā)建立網(wǎng)格衰變光子源抽樣方法，構(gòu)建衰變光子劑量計算蒙特卡羅光子輸運計算模型，調(diào)用MCNP 進行網(wǎng)格活化源項衰變光子活化劑量場計算。

在網(wǎng)格活化分析方法中，有兩大關(guān)鍵問題需解決：1）網(wǎng)格材料等效均勻化計算。針對結(jié)構(gòu)材料網(wǎng)格模型中存在的非均勻材料網(wǎng)格，需要確定其準確材料成分才可實現(xiàn)后續(xù)的網(wǎng)格活化計算；2）網(wǎng)格活化光子源精確抽樣。對于結(jié)構(gòu)材料的網(wǎng)格活化衰變γ源，需要建立網(wǎng)格衰變光子源抽樣方法以實現(xiàn)活化γ劑量場精確計算。

2.2 網(wǎng)格材料等效均勻化計算

核反應(yīng)堆結(jié)構(gòu)材料模型的幾何結(jié)構(gòu)通常較為復雜，幾何邊界多為異形結(jié)構(gòu)。因此，基于網(wǎng)格模型開展結(jié)構(gòu)材料活化分析時，在建立結(jié)構(gòu)材料區(qū)域的網(wǎng)格模型時，網(wǎng)格模型邊界可能不一定與結(jié)構(gòu)材料柵元的外邊界完美契合，存在一個網(wǎng)格包含多個非規(guī)則材料區(qū)域情況（圖3），這將導致一個網(wǎng)格由多種不同的材料組成，此時則無法直接得到其網(wǎng)格材料核素成分，進而無法直接開展網(wǎng)格活化計算。因此，需要對網(wǎng)格材料進行均勻化處理，以得到每個網(wǎng)格中精確的等效均勻材料，用于后續(xù)網(wǎng)格活化計算。具體均勻化如圖3所示。

圖3 網(wǎng)格材料等效均勻化計算方法Fig.3 Calculation method of equivalent homogenization for mesh material compositions

在計算網(wǎng)格等效均勻化材料時，以當前網(wǎng)格中各種材料的體積比例作為權(quán)重因子，通過體積加權(quán)，可以得到網(wǎng)格等效均勻化材料中各核素的質(zhì)量密度或核子密度。網(wǎng)格等效均勻化材料計算公式如下：

式中：mmesh，j與ρmesh，j分別為網(wǎng)格內(nèi)第j種核素的質(zhì)量與密度；i為網(wǎng)格內(nèi)第i種材料；n為當前網(wǎng)格內(nèi)總材料數(shù)；mi，j與ρi，j為當前網(wǎng)格內(nèi)第i種材料中第j種核素的質(zhì)量與密度；Vi為當前網(wǎng)格內(nèi)第i種材料的體積；Vmesh為當前網(wǎng)格的總體積。

網(wǎng)格等效均勻化材料計算流程如圖4所示。為了獲得每個網(wǎng)格中各核素成分的體積加權(quán)因子，本文采用源粒子軌跡統(tǒng)計法實現(xiàn)體積加權(quán)因子的計算。針對每個網(wǎng)格，在網(wǎng)格區(qū)域內(nèi)進行一定數(shù)目源粒子均勻抽樣，并進行蒙特卡羅粒子輸運，統(tǒng)計粒子在網(wǎng)格內(nèi)不同材料區(qū)域內(nèi)的分布情況，進而得到網(wǎng)格內(nèi)不同材料的體積比例，最后按照式（1）進行網(wǎng)格等效均勻材料的計算。

圖4 網(wǎng)格等效均勻材料的計算流程Fig.4 Calculation flow of equivalent homogenization material for mesh

2.3 網(wǎng)格活化γ源粒子抽樣

MCFisMA程序進行網(wǎng)格活化計算之后，將產(chǎn)生精細分布的網(wǎng)格活化光子源，通常具有網(wǎng)格源數(shù)目較多（102~106個網(wǎng)格源區(qū)域）、幾何分布復雜（直角網(wǎng)格或扇形網(wǎng)格）等特點，此時MCNP 程序自帶的源定義功能將無法滿足網(wǎng)格源定義要求。因此，本文對MCNP 程序源粒子抽樣模塊進行二次開發(fā)，以實現(xiàn)以大規(guī)模復雜網(wǎng)格γ 源抽樣，進而實現(xiàn)網(wǎng)格活化γ劑量場的精確計算。

網(wǎng)格γ源粒子抽樣時，主要考慮網(wǎng)格抽樣概率、空間位置抽樣、運動方向抽樣和能量抽樣。其抽樣流程如下：首先，根據(jù)網(wǎng)格的衰變光子源強與整個模型總光子源強比例，得到每個網(wǎng)格的歸一化抽樣概率，根據(jù)抽樣概率可抽樣并確定源粒子所在的網(wǎng)格；其次，根據(jù)圖5的流程進行源粒子坐標、運動方向以及能量的抽樣，進而實現(xiàn)網(wǎng)格光子源抽樣。

3 數(shù)值驗證

3.1 ITER停堆劑量率基準驗證

3.1.1 ITER基準模型介紹

為了驗證本文方法及MCFisMA 程序的正確性，將其應(yīng)用到ITER 國際組織發(fā)布停堆劑量（Shutdown Dose Rate，SDR）基準題中。ITER 例題被廣泛應(yīng)用于驗證活化程序的正確性，模型結(jié)構(gòu)如圖6 所示。在基準題中，由于模型尺度較大且屏蔽材料以鐵水混合物為主，存在較為嚴重的深穿透效應(yīng)。在蒙特卡羅輸運計算過程中采用了重要性減方差方法，首先將屏蔽體分為大小相等的圓環(huán)，再將重要性從模型中子源端開始依次增加，首先嘗試以指數(shù)為倍數(shù)增加，再根據(jù)輸出的tracks entering 和population信息進行調(diào)整，使得進入每一個細分柵元的粒子數(shù)基本相同，即可以開始中子輸運計算。

圖6 ITER停堆劑量率基準題Fig.6 ITER shutdown dose rate benchmark

3.1.2 停堆劑量率結(jié)果分析

劑量率的計數(shù)區(qū)域為圖6中的4個計數(shù)柵元，在進行活化計算時，首先將活化區(qū)域劃分為圓柱網(wǎng)格，分別在r、θ、z坐標軸上劃分10、10、3 個網(wǎng)格，總計300 個圓柱網(wǎng)格進行活化計算。使用14 MeV 的聚變中子源對基準題模型進行輻照，采用SA2-ITER的輻照方案。計算輻照后冷卻106s時4個計數(shù)柵元的停堆劑量率，計算結(jié)果如圖7所示，其中核能與應(yīng)用 實 驗 室（Nuclear Energy and Application Laboratory，NEAL）為本文結(jié)果。

圖7 不同國際機構(gòu)ITER停堆劑量率對比Fig.7 Comparison of ITER SDR benchmark of different international organizations

圖7給出了由本文方法得到的停堆劑量率計算結(jié)果與國際上其他機構(gòu)計算值［3］對比情況。計算結(jié)果顯示，本文結(jié)果趨勢吻合較好，數(shù)值上與英國卡拉姆聚變能研究中心（Culham Centre for Fusion Energy，CCFE）和日本原子能機構(gòu)（Japan Atomic Energy Agency，JAEA）最為接近，其中CCFE 為ITER 總結(jié)報告的參考單位之一。其細微區(qū)別來自于截面庫以及活化計算方法的區(qū)別，以及不同網(wǎng)格劃分方法也會對結(jié)果產(chǎn)生一定偏差?；谝陨辖Y(jié)論，可以證明本文方法和程序的正確性。

3.2 NUREG/CR-6115壓水堆精細活化分析

3.2.1 NUREG/CR-6115壓水堆模型介紹

相比于傳統(tǒng)柵元活化程序，MCFisMA程序可以得到更為精細的結(jié)構(gòu)材料活化光子源分布。為了證明本文方法及程序的精細程度和適用性，將其應(yīng)用到NUREG/CR-6115壓水堆中，其模型見圖8。

圖8 NUREG/CR-6115壓水堆模型Fig.8 The NUREG/CR-6115 PWR reactor model

對NUREG/CR-6115 壓水堆主要結(jié)構(gòu)部件如吊籃、圍板、熱屏蔽層、壓力容器以及壓力容器內(nèi)襯進行了網(wǎng)格活化精細計算。其中，網(wǎng)格活化模型使用柱坐標網(wǎng)格，角向劃分30個網(wǎng)格，軸向劃分128個網(wǎng)格。如圖9 所示，MCFisMA 程序給出了堆芯圍板、堆芯吊籃、熱屏蔽層、壓力容器等結(jié)構(gòu)材料區(qū)域的活化光子源的空間高分辨率精細分布。數(shù)據(jù)結(jié)果表明，活化光子源的源強隨其距堆芯距離的增加呈遞減趨勢，部件的活化程度也隨距堆芯的距離而降低。而壓力容器內(nèi)襯的體積較小，故其衰變光子源的源強遠低于其他結(jié)構(gòu)，四個部件的最大光子源強對比見圖10，其中橫坐標中的1、2、3、4 依次代表堆芯圍板、堆芯吊籃、熱屏蔽層、壓力容器。針對同一結(jié)構(gòu)材料而言，在活性區(qū)高度的中心區(qū)域，由于該區(qū)域中子注量率較高，材料活化程度也較高，衰變光子源的源強也較高。

圖9 NUREG/CR-6115壓水堆結(jié)構(gòu)部件活化光子源項精細分布 (a)堆芯圍板，(b)堆芯吊籃，(c)熱屏蔽層，(d)壓力容器Fig.9 High-resolution distribution of activation photon sources for NUREG/CR-6115 structural components(a)Baffle,(b)Core barrel,(c)Thermal shield,(d)Pressure vessel

圖10 NUREG/CR-6115壓水堆結(jié)構(gòu)部件最大光子源強對比Fig.10 Maximum photon source strength comparison of NUREG/CR-6115 PWR structural components

3.2.2 結(jié)構(gòu)材料活化劑量場分布計算

基于MCFisMA 程序分析了不同網(wǎng)格精細程度對NUREG/CR-6115 壓水堆結(jié)構(gòu)材料活化劑量分析的影響。主要對堆芯圍板、吊籃、熱屏和壓力容器等結(jié)構(gòu)材料進行不同精細程度的軸向網(wǎng)格劃分，分析了不同軸向網(wǎng)格數(shù)目下活化劑量場分布變化情況（圖11）。隨著軸向網(wǎng)格數(shù)目（圖11 中z表示網(wǎng)格數(shù)目）的變化，劑量場分布也呈現(xiàn)變化明顯。軸向網(wǎng)格的數(shù)目較少時，整個劑量場分布相對均勻，同時也存在較大的誤差。隨著網(wǎng)格數(shù)目增加，劑量場分布曲線逐漸收斂，逼近精確值。其主要原因在于軸向方向的結(jié)構(gòu)布置不均勻，盡管在堆芯上下方布置有反射層，但中子注量率分布仍不均勻。在軸向上，靠近堆芯活性區(qū)高度區(qū)域中子注量率較大，因此隨著網(wǎng)格數(shù)目的增加，靠近活性區(qū)高度區(qū)域內(nèi)結(jié)構(gòu)材料活化源項結(jié)果更加精確，使得活化γ 光子所致劑量分布將逐漸收斂于真實值。

圖11 不同軸向網(wǎng)格數(shù)目下的劑量場分布Fig.11 Dose distribution under different axial mesh numbers

同時，為了量化分析網(wǎng)格活化方法相比傳統(tǒng)柵元活化方法對計算結(jié)果的影響，選取了NUREG/CR-6115模型軸向方向的最低點、中間點和最高點進行了不同網(wǎng)格數(shù)目下劑量場對比分析，如表1 所示。在網(wǎng)格數(shù)為1 的情況下，可將其視為傳統(tǒng)柵元活化計算結(jié)果。與柵元活化結(jié)果相比，隨著網(wǎng)格數(shù)目增加，劑量相對誤差逐漸降低，在網(wǎng)格數(shù)達到8之后相對誤差將顯著降低，當軸向劃分為16 個網(wǎng)格時（網(wǎng)格邊長24 cm），得到較為理想的劑量率分布結(jié)果（綜合考慮結(jié)果精度與計算時長，因為隨著網(wǎng)格數(shù)提升計算速度降低）。

表1 不同網(wǎng)格數(shù)目(N)下軸向位置劑量率的相對誤差(%)Table 1 Relative error of dose rate under different mesh numbers(%)

4 結(jié)語

本文針對核反應(yīng)堆結(jié)構(gòu)材料活化的問題，開展了基于兩步法的網(wǎng)格模型精細活化分析方法研究，考慮了中子注量率在空間中的非均勻分布，實現(xiàn)對結(jié)構(gòu)材料活化源項和活化γ源的精細網(wǎng)格空間分布計算，進而實現(xiàn)活化源項精細分析和活化γ 輻射場的精確評價。同時，基于MCNP與FISPACT耦合開發(fā)了全自動耦合網(wǎng)格活化分析程序MCFisMA，通過ITER停堆劑量率基準題驗證程序的正確性，并基于NUREG/CR-6115 壓水堆模型實現(xiàn)了結(jié)構(gòu)材料網(wǎng)格精細活化分析。通過網(wǎng)格活化分析數(shù)值結(jié)果可知，不同精細程度網(wǎng)格模型對活化光子主要貢獻核素、活化光子能譜、活化光子源強均會產(chǎn)生較大影響，尤其對活化劑量場分布影響顯著。相比于傳統(tǒng)柵元活化分析方法，本文網(wǎng)格活化分析方法可以獲得活化核素、光子源等源項的精細空間分布，可評估核反應(yīng)堆復雜結(jié)構(gòu)部件中任意空間位置處的源項分布情況，為反應(yīng)堆屏蔽設(shè)計、檢修換料方案、退役策略提供理論與數(shù)據(jù)支撐。

作者貢獻聲明李雷鳴：主要負責網(wǎng)格活化計算理論研究、程序開發(fā)、起草文章等工作；江書益：主要負責程序數(shù)值基準驗證研究；陳珍平：負責程序在壓水堆中應(yīng)用及網(wǎng)格活化分析研究，提供研究經(jīng)費，對文章進行審閱；苑旭東：負責網(wǎng)格活化源方法及源抽樣程序開發(fā)；楊超：負責對網(wǎng)格源程序驗證研究；謝金森：負責蒙卡粒子輸運計算理論理論指導；于濤：負責蒙卡粒子輸運減方差理論指導、提供研究經(jīng)費及超算平臺，對文章內(nèi)容進行審閱。