核反應(yīng)堆結(jié)構(gòu)材料網(wǎng)格精細(xì)活化分析研究

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

1（南華大學(xué)核科學(xué)技術(shù)學(xué)院 衡陽(yáng)421001）

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

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

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

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

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

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

2.1 嚴(yán)格兩步法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 數(shù)值驗(yàn)證

3.1 ITER停堆劑量率基準(zhǔn)驗(yàn)證

3.1.1 ITER基準(zhǔn)模型介紹

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

圖6 ITER停堆劑量率基準(zhǔn)題Fig.6 ITER shutdown dose rate benchmark

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

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

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

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

3.2 NUREG/CR-6115壓水堆精細(xì)活化分析

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

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

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

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

圖9 NUREG/CR-6115壓水堆結(jié)構(gòu)部件活化光子源項(xiàng)精細(xì)分布 (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)部件最大光子源強(qiáng)對(duì)比Fig.10 Maximum photon source strength comparison of NUREG/CR-6115 PWR structural components

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

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

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

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

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

4 結(jié)語(yǔ)

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

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