李雷鳴 江書益 陳珍平 苑旭東 楊 超 謝金森 于 濤
1(南華大學核科學技術學院 衡陽421001)
2(湖南省數字化反應堆工程技術研究中心 衡陽421001)
在反應堆運行過程中,堆內結構材料受到中子的輻照而發(fā)生活化反應,使其產生大量放射性核素,是反應堆屏蔽設計、檢修換料方案、退役策略中都需要考慮的重要內容[1]。由于反應堆技術的發(fā)展,各式的先進反應堆具有更加復雜的結構以及材料布置與中子能譜,導致中子活化分析也更加復雜,所以有必要開展核反應堆結構材料精細活化分析方法研究[2]。傳統(tǒng)的柵元活化方法將每個柵元的中子注量率取平均值進行活化,無法解決由于柵元內部結構材料活化的非均勻空間效應而引起的活化源項及活化γ源空間分布不夠精確的問題,進而無法對活化γ輻射場進行精確評價[3],對比見圖1。
圖1 核反應堆結構材料活化分析方法對比Fig.1 Comparison of structural material activation analysis methods
為了提高核反應堆結構材料活化分析精度,本文發(fā)展了基于網格活化的核反應堆結構材料精細活化分析方法。所開發(fā)程序在支持直角網格與柱坐標網格的基礎上支持多個非均勻網格拼接,使網格劃分更加靈活的同時也提升了網格的上限。本文基于嚴格兩步法(Rigorous 2-step method,R2S),基于蒙特卡羅核粒子輸運程序系統(tǒng)(Monte Carlo NParticle Transport Code System,MCNP)[4]與活化程序FISPACT[5]耦合,開發(fā)了由“中子輸運”到“材料活化”再到“活化劑量”進行耦合的反應堆結構材料精細化活化計算分析程序(MCNP-FISPACT coupled Mesh-based Activation code,MCFisMA)。
本文基于嚴格兩步法建立核反應堆結構材料網格精細活化分析,其主要包括三步:中子輸運網格計算、網格活化源項計算和衰變光子劑量計算,流程如圖2所示。
圖2 基于R2S的網格活化分析方法Fig.2 Mesh-based activation analysis method based on rigorous 2-step method
1)中子輸運網格計算:建立核反應堆結構材料蒙特卡羅中子輸運計算模型,對結構材料區(qū)域構建網格計算模型,確定中子能群結構,基于MCNP 程序開展蒙特卡羅中子輸運網格計算,得到結構材料區(qū)域的網格中子注量率和網格中子能譜。
2)網格活化源項計算:建立針對結構材料網格模型的網格活化計算模型,根據中子輸運計算得到的網格中子注量率和網格中子能譜,基于FISPACT開展網格活化計算,得到網格活化源項精細空間分布,主要包括網格主要活化核素及貢獻程度、活化光子源強、活化光子能譜等。
3)活化光子劑量計算:根據2)中得到的網格活化衰變光子源分布,基于MCNP 進行二次開發(fā)建立網格衰變光子源抽樣方法,構建衰變光子劑量計算蒙特卡羅光子輸運計算模型,調用MCNP 進行網格活化源項衰變光子活化劑量場計算。
在網格活化分析方法中,有兩大關鍵問題需解決:1)網格材料等效均勻化計算。針對結構材料網格模型中存在的非均勻材料網格,需要確定其準確材料成分才可實現后續(xù)的網格活化計算;2)網格活化光子源精確抽樣。對于結構材料的網格活化衰變γ源,需要建立網格衰變光子源抽樣方法以實現活化γ劑量場精確計算。
核反應堆結構材料模型的幾何結構通常較為復雜,幾何邊界多為異形結構。因此,基于網格模型開展結構材料活化分析時,在建立結構材料區(qū)域的網格模型時,網格模型邊界可能不一定與結構材料柵元的外邊界完美契合,存在一個網格包含多個非規(guī)則材料區(qū)域情況(圖3),這將導致一個網格由多種不同的材料組成,此時則無法直接得到其網格材料核素成分,進而無法直接開展網格活化計算。因此,需要對網格材料進行均勻化處理,以得到每個網格中精確的等效均勻材料,用于后續(xù)網格活化計算。具體均勻化如圖3所示。
圖3 網格材料等效均勻化計算方法Fig.3 Calculation method of equivalent homogenization for mesh material compositions
在計算網格等效均勻化材料時,以當前網格中各種材料的體積比例作為權重因子,通過體積加權,可以得到網格等效均勻化材料中各核素的質量密度或核子密度。網格等效均勻化材料計算公式如下:
式中:mmesh,j與ρmesh,j分別為網格內第j種核素的質量與密度;i為網格內第i種材料;n為當前網格內總材料數;mi,j與ρi,j為當前網格內第i種材料中第j種核素的質量與密度;Vi為當前網格內第i種材料的體積;Vmesh為當前網格的總體積。
網格等效均勻化材料計算流程如圖4所示。為了獲得每個網格中各核素成分的體積加權因子,本文采用源粒子軌跡統(tǒng)計法實現體積加權因子的計算。針對每個網格,在網格區(qū)域內進行一定數目源粒子均勻抽樣,并進行蒙特卡羅粒子輸運,統(tǒng)計粒子在網格內不同材料區(qū)域內的分布情況,進而得到網格內不同材料的體積比例,最后按照式(1)進行網格等效均勻材料的計算。
圖4 網格等效均勻材料的計算流程Fig.4 Calculation flow of equivalent homogenization material for mesh
MCFisMA程序進行網格活化計算之后,將產生精細分布的網格活化光子源,通常具有網格源數目較多(102~106個網格源區(qū)域)、幾何分布復雜(直角網格或扇形網格)等特點,此時MCNP 程序自帶的源定義功能將無法滿足網格源定義要求。因此,本文對MCNP 程序源粒子抽樣模塊進行二次開發(fā),以實現以大規(guī)模復雜網格γ 源抽樣,進而實現網格活化γ劑量場的精確計算。
網格γ源粒子抽樣時,主要考慮網格抽樣概率、空間位置抽樣、運動方向抽樣和能量抽樣。其抽樣流程如下:首先,根據網格的衰變光子源強與整個模型總光子源強比例,得到每個網格的歸一化抽樣概率,根據抽樣概率可抽樣并確定源粒子所在的網格;其次,根據圖5的流程進行源粒子坐標、運動方向以及能量的抽樣,進而實現網格光子源抽樣。
3.1.1 ITER基準模型介紹
為了驗證本文方法及MCFisMA 程序的正確性,將其應用到ITER 國際組織發(fā)布停堆劑量(Shutdown Dose Rate,SDR)基準題中。ITER 例題被廣泛應用于驗證活化程序的正確性,模型結構如圖6 所示。在基準題中,由于模型尺度較大且屏蔽材料以鐵水混合物為主,存在較為嚴重的深穿透效應。在蒙特卡羅輸運計算過程中采用了重要性減方差方法,首先將屏蔽體分為大小相等的圓環(huán),再將重要性從模型中子源端開始依次增加,首先嘗試以指數為倍數增加,再根據輸出的tracks entering 和population信息進行調整,使得進入每一個細分柵元的粒子數基本相同,即可以開始中子輸運計算。
圖6 ITER停堆劑量率基準題Fig.6 ITER shutdown dose rate benchmark
3.1.2 停堆劑量率結果分析
劑量率的計數區(qū)域為圖6中的4個計數柵元,在進行活化計算時,首先將活化區(qū)域劃分為圓柱網格,分別在r、θ、z坐標軸上劃分10、10、3 個網格,總計300 個圓柱網格進行活化計算。使用14 MeV 的聚變中子源對基準題模型進行輻照,采用SA2-ITER的輻照方案。計算輻照后冷卻106s時4個計數柵元的停堆劑量率,計算結果如圖7所示,其中核能與應用 實 驗 室(Nuclear Energy and Application Laboratory,NEAL)為本文結果。
圖7 不同國際機構ITER停堆劑量率對比Fig.7 Comparison of ITER SDR benchmark of different international organizations
圖7給出了由本文方法得到的停堆劑量率計算結果與國際上其他機構計算值[3]對比情況。計算結果顯示,本文結果趨勢吻合較好,數值上與英國卡拉姆聚變能研究中心(Culham Centre for Fusion Energy,CCFE)和日本原子能機構(Japan Atomic Energy Agency,JAEA)最為接近,其中CCFE 為ITER 總結報告的參考單位之一。其細微區(qū)別來自于截面庫以及活化計算方法的區(qū)別,以及不同網格劃分方法也會對結果產生一定偏差?;谝陨辖Y論,可以證明本文方法和程序的正確性。
3.2.1 NUREG/CR-6115壓水堆模型介紹
相比于傳統(tǒng)柵元活化程序,MCFisMA程序可以得到更為精細的結構材料活化光子源分布。為了證明本文方法及程序的精細程度和適用性,將其應用到NUREG/CR-6115壓水堆中,其模型見圖8。
圖8 NUREG/CR-6115壓水堆模型Fig.8 The NUREG/CR-6115 PWR reactor model
對NUREG/CR-6115 壓水堆主要結構部件如吊籃、圍板、熱屏蔽層、壓力容器以及壓力容器內襯進行了網格活化精細計算。其中,網格活化模型使用柱坐標網格,角向劃分30個網格,軸向劃分128個網格。如圖9 所示,MCFisMA 程序給出了堆芯圍板、堆芯吊籃、熱屏蔽層、壓力容器等結構材料區(qū)域的活化光子源的空間高分辨率精細分布。數據結果表明,活化光子源的源強隨其距堆芯距離的增加呈遞減趨勢,部件的活化程度也隨距堆芯的距離而降低。而壓力容器內襯的體積較小,故其衰變光子源的源強遠低于其他結構,四個部件的最大光子源強對比見圖10,其中橫坐標中的1、2、3、4 依次代表堆芯圍板、堆芯吊籃、熱屏蔽層、壓力容器。針對同一結構材料而言,在活性區(qū)高度的中心區(qū)域,由于該區(qū)域中子注量率較高,材料活化程度也較高,衰變光子源的源強也較高。
圖9 NUREG/CR-6115壓水堆結構部件活化光子源項精細分布 (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壓水堆結構部件最大光子源強對比Fig.10 Maximum photon source strength comparison of NUREG/CR-6115 PWR structural components
3.2.2 結構材料活化劑量場分布計算
基于MCFisMA 程序分析了不同網格精細程度對NUREG/CR-6115 壓水堆結構材料活化劑量分析的影響。主要對堆芯圍板、吊籃、熱屏和壓力容器等結構材料進行不同精細程度的軸向網格劃分,分析了不同軸向網格數目下活化劑量場分布變化情況(圖11)。隨著軸向網格數目(圖11 中z表示網格數目)的變化,劑量場分布也呈現變化明顯。軸向網格的數目較少時,整個劑量場分布相對均勻,同時也存在較大的誤差。隨著網格數目增加,劑量場分布曲線逐漸收斂,逼近精確值。其主要原因在于軸向方向的結構布置不均勻,盡管在堆芯上下方布置有反射層,但中子注量率分布仍不均勻。在軸向上,靠近堆芯活性區(qū)高度區(qū)域中子注量率較大,因此隨著網格數目的增加,靠近活性區(qū)高度區(qū)域內結構材料活化源項結果更加精確,使得活化γ 光子所致劑量分布將逐漸收斂于真實值。
圖11 不同軸向網格數目下的劑量場分布Fig.11 Dose distribution under different axial mesh numbers
同時,為了量化分析網格活化方法相比傳統(tǒng)柵元活化方法對計算結果的影響,選取了NUREG/CR-6115模型軸向方向的最低點、中間點和最高點進行了不同網格數目下劑量場對比分析,如表1 所示。在網格數為1 的情況下,可將其視為傳統(tǒng)柵元活化計算結果。與柵元活化結果相比,隨著網格數目增加,劑量相對誤差逐漸降低,在網格數達到8之后相對誤差將顯著降低,當軸向劃分為16 個網格時(網格邊長24 cm),得到較為理想的劑量率分布結果(綜合考慮結果精度與計算時長,因為隨著網格數提升計算速度降低)。
表1 不同網格數目(N)下軸向位置劑量率的相對誤差(%)Table 1 Relative error of dose rate under different mesh numbers(%)
本文針對核反應堆結構材料活化的問題,開展了基于兩步法的網格模型精細活化分析方法研究,考慮了中子注量率在空間中的非均勻分布,實現對結構材料活化源項和活化γ源的精細網格空間分布計算,進而實現活化源項精細分析和活化γ 輻射場的精確評價。同時,基于MCNP與FISPACT耦合開發(fā)了全自動耦合網格活化分析程序MCFisMA,通過ITER停堆劑量率基準題驗證程序的正確性,并基于NUREG/CR-6115 壓水堆模型實現了結構材料網格精細活化分析。通過網格活化分析數值結果可知,不同精細程度網格模型對活化光子主要貢獻核素、活化光子能譜、活化光子源強均會產生較大影響,尤其對活化劑量場分布影響顯著。相比于傳統(tǒng)柵元活化分析方法,本文網格活化分析方法可以獲得活化核素、光子源等源項的精細空間分布,可評估核反應堆復雜結構部件中任意空間位置處的源項分布情況,為反應堆屏蔽設計、檢修換料方案、退役策略提供理論與數據支撐。
作者貢獻聲明李雷鳴:主要負責網格活化計算理論研究、程序開發(fā)、起草文章等工作;江書益:主要負責程序數值基準驗證研究;陳珍平:負責程序在壓水堆中應用及網格活化分析研究,提供研究經費,對文章進行審閱;苑旭東:負責網格活化源方法及源抽樣程序開發(fā);楊超:負責對網格源程序驗證研究;謝金森:負責蒙卡粒子輸運計算理論理論指導;于濤:負責蒙卡粒子輸運減方差理論指導、提供研究經費及超算平臺,對文章內容進行審閱。