吳明昌，吳亮亮，郝麗娟，鄭華慶，宋 婧，F(xiàn)DS團(tuán)隊(duì)

（1.中國(guó)科學(xué)技術(shù)大學(xué)，安徽 合肥 230027；2.中國(guó)科學(xué)院核能安全技術(shù)研究所，安徽 合肥230031）

托卡馬克核聚變裝置在進(jìn)行等離子體放電實(shí)驗(yàn)時(shí)，將產(chǎn)生聚變中子。裝置停止運(yùn)行后，結(jié)構(gòu)材料受到這些聚變中子照射產(chǎn)生衰變光子，從而引起核聚變裝置的停機(jī)輻射劑量。近年來隨著核聚變裝置的告訴發(fā)展、實(shí)驗(yàn)功率水平的不斷提高，D－T或D－D等離子體實(shí)驗(yàn)產(chǎn)生大量聚變中子，裝置結(jié)構(gòu)材料的活化情況嚴(yán)重，導(dǎo)致裝置周圍的停機(jī)劑量難以在較短的時(shí)間內(nèi)衰減到可以忽略的程度。裝置運(yùn)行停止后，實(shí)驗(yàn)人員有時(shí)需要進(jìn)入裝置內(nèi)部進(jìn)行實(shí)驗(yàn)測(cè)量工作，較高的輻射劑量會(huì)對(duì)實(shí)驗(yàn)人員的人身安全產(chǎn)生影響。因此，需要準(zhǔn)確地評(píng)估實(shí)驗(yàn)停止以后裝置周圍的停機(jī)劑量率水平，為裝置輻射屏蔽系統(tǒng)的設(shè)計(jì)、實(shí)驗(yàn)方案的制訂以及實(shí)驗(yàn)人員允許進(jìn)入裝置內(nèi)部的時(shí)間等提供重要的參考[1］。

隨著對(duì)停機(jī)劑量率分析研究工作的深入，文獻(xiàn)[2］提出了一種基于網(wǎng)格計(jì)數(shù)的停機(jī)劑量率計(jì)算方法——嚴(yán)格二步法，用于精確計(jì)算幾何結(jié)構(gòu)復(fù)雜的托卡馬克裝置的停堆劑量率。根據(jù)該方法，F(xiàn)DS團(tuán)隊(duì)開發(fā)了一套精確的停機(jī)劑量率計(jì)算程序，并將其集成進(jìn)FDS團(tuán)隊(duì)開發(fā)的大型集成中子學(xué)計(jì)算分析系統(tǒng) VisualBUS[3－4］中。該程序通過耦合三維蒙特卡羅程序MCNP[5］和歐洲活化計(jì)算程序 FISPACT[6］，經(jīng)過一系列中間計(jì)算，最終獲得裝置停機(jī)后周圍空間的劑量率分布。

由于聚變裝置幾何和材料的高度復(fù)雜性，其基于網(wǎng)格計(jì)數(shù)的停機(jī)劑量率計(jì)算中涉及的光子源的柵元數(shù)目超過了1 000個(gè)，難以使用MCNP的標(biāo)準(zhǔn)源卡描述。為了解決SDEF卡對(duì)源柵元數(shù)目的限制，本文使用了MCNP自定義源子程序方法對(duì)衰變光子源進(jìn)行直接抽樣。

1 停機(jī)劑量率的計(jì)算

聚變堆的停機(jī)劑量率是由聚變中子活化的結(jié)構(gòu)材料所發(fā)射出的光子造成的。根據(jù)直接二步法，停機(jī)劑量率的計(jì)算過程如下：

（1）中子輸運(yùn)計(jì)算，計(jì)算運(yùn)行時(shí)中子注量率的空間分布；

（2）活化計(jì)算，使用第一步得到的中子注量率分布，求出衰變光子在時(shí)間、空間及能量上的分布；

（3）光子輸運(yùn)計(jì)算，根據(jù)第二步計(jì)算得到的衰變光子分布，求出裝置停止運(yùn)行后，其系統(tǒng)在不同時(shí)間段的空間劑量率分布。

2 停機(jī)劑量率源子程序的設(shè)計(jì)實(shí)現(xiàn)

2.1 總體設(shè)計(jì)

中子注量率分布與光子分布的結(jié)果的精確性直接影響著劑量率計(jì)算結(jié)果的精確性。為了精確計(jì)算劑量率，停機(jī)劑量率計(jì)算利用MCNP網(wǎng)格計(jì)算能力，將停機(jī)劑量計(jì)算所涉及的幾何劃分成若干網(wǎng)格，光子計(jì)算時(shí)網(wǎng)格被視為光子源柵元。當(dāng)網(wǎng)格劃分的越多越精細(xì)時(shí)，中子注量率分布以及光子分布也就計(jì)算得越精確，從而劑量率結(jié)果便就越精確。當(dāng)網(wǎng)格數(shù)目（也即光子計(jì)算時(shí)源柵元數(shù)目）過多時(shí)，MCNP的SDEF卡無法描述，系統(tǒng)需要通過MNCP源子程序?qū)庾拥姆植贾苯舆M(jìn)行抽樣。

每個(gè)網(wǎng)格中的材料因受到的中子輻射程度不同，被活化的程度也不同，光子源的分布特性也不相同。因此，在使用源子程序?qū)庾釉催M(jìn)行抽樣時(shí)，應(yīng)先根據(jù)每個(gè)網(wǎng)格的光子強(qiáng)度抽樣決定光子所在的網(wǎng)格（每個(gè)網(wǎng)格的光子強(qiáng)度由第二步計(jì)算得到，存儲(chǔ)在sourceShape與source文件中，其計(jì)算方法與格式見文獻(xiàn)[7］）。每個(gè)網(wǎng)格內(nèi)的材料是均勻的，這些材料被活化后出射的衰變光子也是在網(wǎng)格內(nèi)均勻分布，因此光子的位置在網(wǎng)格內(nèi)均勻抽樣。大部分衰變核素都具有各向同性特點(diǎn)，光子的方向根據(jù)該特點(diǎn)按照各向同性處理。然后根據(jù)第二步計(jì)算得到的光子能譜抽樣光子的能量，最后根據(jù)光子的位置計(jì)算光子所在的柵元。源子程序的流程如圖1所示。

2.2 隨機(jī)抽樣的基本方法

源子程序方法是根據(jù)源的分布規(guī)律抽樣而得到源粒子的位置、方向、能量、權(quán)重、壽命的方法。

由已知分布的隨機(jī)抽樣，是借助于隨機(jī)數(shù)進(jìn)行的，MCNP提供了RANG（）函數(shù)產(chǎn)生隨機(jī)數(shù)。該函數(shù)可以產(chǎn)生[0，1）均勻分布的偽隨機(jī)數(shù)。用戶在編寫源子程序依據(jù)源的分布規(guī)律對(duì)源的狀態(tài)變量進(jìn)行隨機(jī)抽樣時(shí)，可以調(diào)用該函數(shù)。

圖1 源子程序流程Fig.1 Procedure of Source subroutine

對(duì)于任意離散型分布

其中x1，x2，…為離散型分布函數(shù)的跳躍點(diǎn)，P1，P2，…為相應(yīng)的概率。F（x）的直接抽樣方法為對(duì)于連續(xù)型分布，如果分布函數(shù)F（x）的反函數(shù)F－1（x）存在。則抽樣方法是

其中ξ是[0，1）均勻分布的隨機(jī)數(shù)，用戶可以在源子程序中調(diào)用RANG（）產(chǎn)生該隨機(jī)數(shù)。

2.3 源子程序變量的抽樣

停機(jī)劑量率計(jì)算中，衰變光子的狀態(tài)取決于光子出射的網(wǎng)格，所以源子程序所涉及第一個(gè)隨機(jī)抽樣過程是確定源粒子所處的網(wǎng)格。這是一個(gè)離散型隨機(jī)分布，參考抽樣公式（2），其抽樣方法為

當(dāng)確定了源粒子所在的網(wǎng)格后，就根據(jù)該網(wǎng)格的位置和尺寸在網(wǎng)格均勻抽樣以確定源粒子的具體位置（xxx，yyy，zzz）。

下一步即是確定光子的方向，由于放射性核素的衰變各向同性特點(diǎn)，故而光子的方向按照各項(xiàng)同性處理，即方位角φ根據(jù)U[0，2π］分布，極角的余弦值μ遵循U[－1，1］分布，根據(jù)公式（3）可以得到方向矢量（uuu，vvv，www）為

然后利用網(wǎng)格的光子譜和光子強(qiáng)度信息，計(jì)算出光子的能量erg。最后對(duì)粒子賦初始權(quán)重值（通常為1.0）和壽命（通常為0.0，表示粒子剛剛產(chǎn)生）。至此粒子的位置（xxx，yyy，zzz）、方向（uuu，vvv，www）、能量（erg）、權(quán)重（wgt）以及壽命（tme）都已經(jīng)抽樣得出。

MCNP要求在任何時(shí)刻粒子所在的柵元都必須是明確的，并由icl保存該柵元的編號(hào)。因此，用戶編寫源子程序時(shí)，還需要確定粒子所在的幾何體的編號(hào)。為確定粒子所在的柵元，源子程序會(huì)根據(jù)粒子的位置（xxx，yyy，zzz）逐個(gè)檢查所有柵元，最終確定粒子所在的柵元。MCNP提供了CHECEL子程序逐個(gè)檢查柵元，用戶可以調(diào)用該子程序?qū)崿F(xiàn)幾何檢測(cè)。該子程序的用法參考文獻(xiàn)[5］。

3 基準(zhǔn)測(cè)試

本文采用ITER（國(guó)際熱核聚變實(shí)驗(yàn)堆，International Thermonuclear Experimental Reactor）組織發(fā)布的停機(jī)劑量率基準(zhǔn)例題[9］以及ITER－T426停機(jī)劑量率實(shí)驗(yàn)例題[10］對(duì)源子程序進(jìn)行校驗(yàn)。

本文使用核與輻射輸運(yùn)計(jì)算自動(dòng)建模軟件MCAM[11－14］對(duì)ITER 停機(jī)劑量率基準(zhǔn)例題的模型的幾何進(jìn)行劃分，共劃分為1 210個(gè)網(wǎng)格，采用源子程序?qū)υ撃Ｐ偷耐６褎┝柯蔬M(jìn)行計(jì)算，結(jié)果如圖2所示。本程序的計(jì)算結(jié)果（FDS）與英國(guó)卡拉姆聚變中心（CCFE）以及荷蘭核研究與咨詢研究集團(tuán)（NRG）的結(jié)果相近。

圖2 ITER基準(zhǔn)測(cè)試?yán)}模型的網(wǎng)格劃分圖Fig.2 Mesh of ITER benchmark

ITER－T426例題的幾何被劃分為5 000多個(gè)網(wǎng)格。計(jì)算程序?qū)υ摾}的計(jì)算結(jié)果如圖3所示，程序計(jì)算結(jié)果與實(shí)驗(yàn)測(cè)量結(jié)果趨勢(shì)相吻合。

圖3 不同程序?qū)TER基準(zhǔn)例題停機(jī)劑量計(jì)算結(jié)果比較[7］Fig.3 Shutdown dose rate calculated by different codes of ITER benchmark[7］

從該圖中可以看出，當(dāng)冷卻時(shí)間比較短時(shí)，程序計(jì)算停堆劑量率結(jié)果比實(shí)驗(yàn)結(jié)果最大處約低估9.4%，隨著冷卻時(shí)間的增加，誤差穩(wěn)定在3%左右。根據(jù)前文可知，計(jì)算的精確性主要依賴于中子通量和光子分布的計(jì)算結(jié)果，而其計(jì)算又依賴于網(wǎng)格的劃分精細(xì)程度，所以源子程序計(jì)算結(jié)果的誤差與網(wǎng)格劃分的精細(xì)度程序有關(guān)。

在冷卻階段的早期，停機(jī)劑量主要由放射性核素55Mn的衰變光子產(chǎn)生，其半衰期為2.58h。隨后58Co成為光子的主要來源，其半衰期為70.86d。理論上，劑量率的在早期階段衰減得比較迅速，隨后漸漸變緩。實(shí)驗(yàn)結(jié)果與理論計(jì)算都符合預(yù)期。

圖4 ITER－T426網(wǎng)格劃分Fig.4 Mesh of ITER－T426

圖5 ITER－T426劑量率隨時(shí)間變化圖[7］Fig.5 Curves of shutdown dos rate of ITER－T426in time[7］

4 結(jié)論

針對(duì)基于網(wǎng)格計(jì)數(shù)的停機(jī)劑量率計(jì)算時(shí)MCNP通用源卡柵元個(gè)數(shù)受限制問題，本文采用MCNP自定義源子程序技術(shù)克服這一問題。文章系統(tǒng)地介紹了源子程序的編寫方法，并將該方法應(yīng)用于網(wǎng)格計(jì)數(shù)的停機(jī)劑量率計(jì)算程序中。通過國(guó)際基準(zhǔn)的停堆劑量率實(shí)驗(yàn)例題對(duì)源子程序的可用性與正確性進(jìn)行了驗(yàn)證。

[1]吳亮亮，應(yīng)棟川，邱岳峰，等.三維停機(jī)劑量率計(jì)算程序研發(fā)及其在EAST上的應(yīng)用[J］.核科學(xué)與工程，2011，21（1）：80－82.

[2]陳義學(xué)，吳宜燦，U.Fischer.聚變裝置停機(jī)劑量率計(jì)算的嚴(yán)格兩步（R2S）法[J］.核技術(shù)，2003，26（10），763－764.

[3]吳宜燦，李靜驚，李瑩，等.大型集成多功能中子學(xué)計(jì)算與分析系統(tǒng)VisualBUS的研究與發(fā)展[J］.核科學(xué)與工程，2007，27（4）：365－368.

[4]曾勤，邱岳峰，王國(guó)忠，等.大型集成中子學(xué)程序系統(tǒng)VisualBUS研發(fā)進(jìn)展[A］.第四屆全國(guó)反應(yīng)堆物理與核材料學(xué)術(shù)研討會(huì)論文集（C），2009.

[5]X－5Monte Carlo Team. MCNP5－A General Monte Carlo N－Particle Transport Code，Volume II：Users Guide[R］.New Mexico：Los Alamos National Library，2003.

[6]Forrest R A.FISPACT 2007：User Manual [R］.UKAEA Fusion，Report UKAEA FUS 534，March 2007.

[7]吳亮亮.基于網(wǎng)格計(jì)數(shù)的停機(jī)劑量率計(jì)算程序設(shè)計(jì)研究[D］.中國(guó)科學(xué)技術(shù)大學(xué)碩士論文，2012.

[8]許淑艷.蒙特卡羅方法在實(shí)驗(yàn)核物理中的應(yīng)用[M］.北京：原子能出版社，2006.

[9]應(yīng)棟川.ITER偏慮器中子學(xué)性能及其更換維修輻射劑量場(chǎng)分析研究[D］.中國(guó)科學(xué)院等離子體物理研究所，2011.

[10]P.Batistoni. Experimental validation of shutdown dose rates calculations inside ITER cryostat[J］.Fusion Engineering and Design，2011，58－59：613－616.

[11]吳宜燦，李瑩，盧磊，等.蒙特卡羅粒子輸運(yùn)計(jì)算自動(dòng)建模程序系統(tǒng)的研究與發(fā)展[J］.核科學(xué)與工程，2006，26（1）：20－27.

[12]曾勤，盧磊，李瑩，等.蒙特卡羅粒子輸運(yùn)計(jì)算自動(dòng)建模程序MCAM在ITER核分析建模中的應(yīng)用[J］.原子核物理評(píng)論，2006，23（2）：138－141.

[13]王國(guó)忠，黨同強(qiáng)，熊健，等.MCAM4.8在ITER建筑大廳中子學(xué)建模中的應(yīng)用[J］.核科學(xué)與工程，2011，31（4）：352－353.

[14]吳宜燦，胡麗琴，龍鵬程，等.先進(jìn)核能系統(tǒng)設(shè)計(jì)分析軟件與數(shù)據(jù)庫(kù)研發(fā)進(jìn)展[J］.核科學(xué)與工程，2010，30（1）：55－64.