鄭 劍，李 斌，鄒 俊，楊 琪

(中國科學院核能安全技術研究所，中國科學院中子輸運理論與輻射安全重點實驗室，合肥230031)

基于網(wǎng)格計數(shù)的停機劑量率計算方法研究

鄭 劍，李 斌，鄒 俊，楊 琪

(中國科學院核能安全技術研究所，中國科學院中子輸運理論與輻射安全重點實驗室，合肥230031)

本文根據(jù)嚴格二步法的計算理論，研究了基于網(wǎng)格計數(shù)的停機劑量率計算方法，設計并實現(xiàn)了基于網(wǎng)格計數(shù)的停機劑量率計算程序。該程序能夠支持圓柱坐標下的網(wǎng)格計算。本文使用源子程序進行復雜源描述。為了加快計算速度，本文采用了多節(jié)點和多線程等技術。本文利用國際熱核聚變實驗堆(ITER)停堆劑量基準實驗ITER-T426進行測試，計算結果與實驗值吻合良好，證明了該方法和程序的正確性和可用性。

網(wǎng)格計數(shù)；停機劑量率；嚴格二步法；自動耦合

核裝置運行停止或運行間歇期間，材料活化所釋放的衰變光子對需要靠近或進入裝置內(nèi)部進行實驗測量或維修檢測工作的人員造成嚴重的輻射傷害。因此準確地評估核裝置的停堆劑量率水平，對于核裝置屏蔽設計、維修計劃的制訂以及核裝置的退役有著重要的參考意義。

本文在FDS團隊自主研發(fā)的大型集成多功能中子學計算分析系統(tǒng)VisualBUS[3]及已經(jīng)開發(fā)的停機劑量率計算程序的基礎上[4-5]，研究了基于網(wǎng)格計數(shù)的嚴格二步法，設計并實現(xiàn)了基于網(wǎng)格計數(shù)的停堆劑量率計算程序。該程序支持圓柱坐標下的網(wǎng)格計算，使用源子程序及并行技術，解決了傳統(tǒng)停機劑量率計算過程中的精度不夠、模型劃分困難、源柵元數(shù)受限等問題，實現(xiàn)了核裝置停機劑量率的精確計算，計算結果可以為裝置的屏蔽設計、實驗方案的制訂以及工作人員允許靠近裝置的等待時間等提供重要的參考依據(jù)。

1 計算方法

1.1 總體設計

網(wǎng)格計數(shù)停機劑量率計算方法基于柵元計數(shù)嚴格二步法的思想，充分利用MCNP程序的網(wǎng)格計數(shù)能力。圖1描述了該方法的計算流程。

圖1 計算流程圖Fig.1 Calculation flow chart

該方法首先使用MCNP程序[6]，對已劃分了網(wǎng)格的計算模型進行網(wǎng)格計數(shù)中子輸運計算，獲得所有被計數(shù)網(wǎng)格的中子分群通量。然后調(diào)用MCNP程序?qū)γ恳粋€網(wǎng)格進行PTRAC計算，獲得每個網(wǎng)格的粒子分布情況，用于計算網(wǎng)格的平均材料。

然后，該方法根據(jù)上面得到的網(wǎng)格內(nèi)通量和材料，以及用戶提供的輻照方案，通過使用FISPACT程序[7]對每一個計數(shù)網(wǎng)格進行活化計算，得到不同冷卻時間下所有網(wǎng)格的衰變光子分布情況。之后，提取活化之后的光子譜及強度，生成光子源輸入文件。

最后，使用MCNP程序?qū)Σ煌瑫r間段下的衰變光子進行網(wǎng)格計數(shù)輸運計算，最終獲得不同時間段下計數(shù)網(wǎng)格的劑量率分布。

該方法通過接口程序，自動完成計算過程中信息的處理及傳遞，實現(xiàn)MCNP和FISPACT的雙向自動耦合，避免了其中耗時的、復雜困難易出錯的手工處理過程。

利用比濁法的原理，用光密度值(OD值)的變化來判斷菌株存活情況。根據(jù)朗伯比爾定律OD=log It IO=kbc，c表示樣品濃度，若樣品液厚度b一定，則OD值與樣品的濁度相關[11]。將活化好的菌種菌液以1%的接種比例，分別加入用1 mol/L HCl和1 mol/L NaOH調(diào)pH值(3，4，5，6，7)的滅菌YPD液體培養(yǎng)基中，28 ℃、180 r/min振蕩培養(yǎng)48 h后，將菌液稀釋至合適倍數(shù)(使吸光值在0.1～0.8之間)，在波長600 nm處測量吸光度值，每組試驗做3組平行。

1.2 網(wǎng)格材料計算

MCNP程序提供了一種名為粒子徑跡輸出卡(PTRAC: Particle Track Output)的計算功能。通過在MCNP輸入文件中填寫PTRAC卡，并調(diào)用MCNP程序進行輸運計算，得到名為ptrac的輸出文件，該文件中打印出了計數(shù)粒子信息。統(tǒng)計網(wǎng)格中這些計數(shù)粒子所在的柵元信息，可以得知網(wǎng)格所含柵元情況，從而計算出網(wǎng)格的平均材料。計算公式如下：

(1)

由于PTRAC卡及源描述不支持如圖2所示的圓柱幾何劃分出的扇形網(wǎng)格，本文研究并提出了一種方法來計算這類扇形幾何網(wǎng)格內(nèi)的材料：首先，該程序根據(jù)網(wǎng)格的邊界設定了立方體型的包圍體，在此包圍體內(nèi)進行PTRAC計算，然后去掉不在網(wǎng)格內(nèi)的源粒子，接著依照上文所述方法計算網(wǎng)格材料組成。

圖2 單元與包圍體Fig.2 Cell and bounding volume

1.3 源子程序

本文開發(fā)了一種基于網(wǎng)格的衰變光子取樣程序，并將該程序耦合到MCNP程序中，圖3給出了源子程序抽樣流程。

圖3 源子程序流程圖Fig.3 Source subroutine flowchart

源輸入文件中包含了網(wǎng)格的幾何坐標、體積、衰變光子強度、衰變光子譜、柵元號以及光子能群結構信息。源子程序使用隨機數(shù)函數(shù)rang()在包圍盒內(nèi)進行抽樣。

首先，源子程序按一定概率進行網(wǎng)格抽樣。對于不同網(wǎng)格，其體積和光子強度值越大被抽樣到的概率越大，抽樣概率計算公式如下：

(2)

其中，Vi表示第i個網(wǎng)格的體積，Di表示第i個網(wǎng)格的光子強度，Ri表示從n個網(wǎng)格中抽取的第i個網(wǎng)格的體積與光子強度值。

然后，在抽樣得到的網(wǎng)格內(nèi)進行均勻光子抽樣，通過公式(3)確定立方體網(wǎng)格內(nèi)衰變光子的幾何坐標(xxx,yyy,zzz)；通過公式(4)確定圓柱網(wǎng)格內(nèi)衰變光子的半徑r，然后通過坐標變換得到幾何坐標(xxx,yyy,zzz)。對于得到的衰變光子坐標，調(diào)用MCNP程序提供的庫函數(shù)確定該光子所在的柵元號icl，本文支持對重復結構的柵元處理。

x=(xmax-xmin)·rang()+xmin

(3)

(4)

最后，利用網(wǎng)格的光子譜和光子強度信息，計算出光子的能量erg。另外本文還修改了MCNP部分模塊，使其能夠支持光子輸運計算的并行化。

1.4 并行計算

精確計算具有復雜結構的裝置的停機劑量率，需要對其進行數(shù)量龐大的網(wǎng)格劃分。為了解決計算耗時問題，本文采用多節(jié)點和多線程計算技術。

程序設計過程中，以其中一個節(jié)點作為主節(jié)點運行網(wǎng)格停機劑量率計算程序，負責計算文件的發(fā)送與結果文件的接收。其余的節(jié)點負責調(diào)用FISPACT進行活化計算，并將計算輸出文件發(fā)送給主節(jié)點。主節(jié)點上的主程序完成計算輸入文件準備之后，啟動兩個線程，一個用于將計算文件平均發(fā)送給計算節(jié)點進行計算；另一個負責監(jiān)聽計算節(jié)點傳遞結果文件的請求，一旦接收到計算節(jié)點的請求，主程序啟動一個新的線程接收這個計算節(jié)點發(fā)回的計算結果文件。所有計算節(jié)點一直處于監(jiān)聽狀態(tài)，一旦監(jiān)聽到接收文件的請求，開始接收文件。當所有文件接收完畢，關閉通信線路，并開始調(diào)用FISPACT計算。當所有文件計算完畢，計算節(jié)點再次連接主節(jié)點，并請求發(fā)送文件。

其次，在計算節(jié)點該程序創(chuàng)建并啟動多個子線程，每個線程函數(shù)調(diào)用FISPACT進行計算，計算完成后子線程自動結束。主線程和子線程之間通過全局變量實現(xiàn)通信，主線程必須等待所有子線程結束后，再進行后續(xù)工作。各個子線程之間沒有信息傳遞，從而簡化了實現(xiàn)步驟。采用多線程的方式提高了CPU的利用率，減少了計算時間。

2 ITER-T426實驗例題測試

在意大利開展的ITER-T426停機劑量率實驗使用中子產(chǎn)生器FNG(Frascati Neutron Generator)產(chǎn)生的D-T中子來照射實驗裝置，然后使用探測器測量不同冷卻時間該裝置中的劑量率，將實驗結果與程序計算結果進行比較分析，以校核相關計算程序、方法及核數(shù)據(jù)庫的可靠性[8-10]。

ITER-T426模型由兩部分組成，中子產(chǎn)生器模型及立體屏蔽層模型，本文使用超級蒙卡核模擬軟件系統(tǒng)SuperMC[11-14]建立了該裝置模型，如圖4所示。

圖4 ITER-T426模型Fig.4 ITER-T426 model

在網(wǎng)格計數(shù)停機劑量率計算中，對于中子輸運計算，本文使用HENDL3.0[15-16]數(shù)據(jù)庫計算計數(shù)網(wǎng)格的中子通量分布。FISPACT程序活化計算采用活化庫EAF-2007[7]。中子及衰變光子輸運計算使用相同的計算模型。本文在圓柱坐標系下在模型R、θ、Z方向分別劃分43×8×66計數(shù)網(wǎng)格。在衰變光子的輸運計算中，本文使用衰變光子取樣程序進行衰變光子源抽樣。在光子輸運計算時，使用ITER推薦的通量劑量轉(zhuǎn)換因子[16]計算出模型在不同停機時間的劑量率。計算結果的可視化分析使用FDS團隊自主開發(fā)的科學計算可視化分析軟件SVIP[18-20]。圖5是該裝置的中子通量分布的剖面圖；圖6顯示了停堆12天后的光子源分布情況，可以看出右側(cè)的不銹鋼部分活化比較嚴重；圖7展示了在圖6的衰變光子源輸運之后的分布情況。

圖5 圓柱剖面的中子通量場(cm-2s-1)Fig.5 The profile of the neutron flux map (cm-2s-1)

圖6 停機12天后中心剖面衰變光子強度分布(cm-3s-1)Fig.6 The profile of intensity of decay gamma sources 12 days after shutdown (cm-3s-1)

圖7 停機12天后y-z平面光子場分布(cm-2s-1)Fig.7 The profile of photon flux map 12 days after shutdown (cm-2s-1)

圖8 計算值與實驗值的比值Fig.8 Comparison between calculated and experiment dose rate in the cavity center.

圖8給出了網(wǎng)格計數(shù)停機劑量率程序計算結果與實驗測量結果隨冷卻時間的變化。從該圖中可以看出，當冷卻時間比較短時，程序計算結果比實驗結果最大處約低估25%。這是因為計算使用的輻照方案是簡化的，不能夠完全模擬實際的實驗情況，并且有些對劑量貢獻較多的核素沒有被考慮。而冷卻時間比較長時，程序計算結果最高比實驗結果高出10%，此時的誤差的主要來自于網(wǎng)格劃分不夠細。本程序計算結果與實驗測量結果以及其他程序計算結果趨勢相吻合。

3 結束語

本文根據(jù)嚴格二步法的計算理論，研究了基于網(wǎng)格計數(shù)的停機劑量率計算方法，設計并實現(xiàn)了基于網(wǎng)格計數(shù)的停機劑量率計算程序。

為了驗證方法和程序的正確性，利用國際熱核聚變實驗堆停堆劑量基準實驗ITER-T426實驗進行測試，計算結果證明了該方法和程序的正確性以及處理復雜幾何裝置的適用性。

致謝

本工作得到中科院核能安全技術研究所·FDS團隊其他成員的大力幫助和支持，在此深表感謝！

[1] Davide Valenza, Hiromasa Iida, Romano Plenteda, et al, Proposal of Shutdown Dose Rate Estimation Method by Monte Carlo Code[J]. Fusion Engineering and Design. 2001, 55: 411-418.

[2] Y. Chen, Y. Wu, U. Fischer, The rigorous 2-step calculation of shutdown dose rates for nuclear fusion devices[J]. Nuclear technology. 2003, 26: 763-764.

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

[4] L. L. Wu, Q. Yang, J. Zou, et al. Shutdown dose rate calculation code system and its application to EAST tokamak. Fusion Engineering and Design, 2012, 87: 1315-1318.

[5] 吳亮亮，應棟川，邱岳峰，等. 三維停堆劑量率計算程序研發(fā)及其在EAST上的應用[J]. 核科學與工程，2011, 31(1): 80-85.

[6] MCNP Version 5, Diagnostic applications group(X-5)[J]. Los Alamos National Laboratory.5,2003:334-567.

[7] R.A. Forrest, FISPACT2007: User Manual[M]. UKAEA Fusion. 534, 2007:12-22.

[8] P. Batistoni, M. Angelone, L. Petrizzi, et al, Experi-mental Validation of Shut-Down Dose Rates, Final Report, June 2001.

[9] P. Batistoni, M. Angelone, L. Petrizzi, et al, “Bench-mark Experiment for the validation of shut down activa-tion and dose calculation in a fusion device”, Journal of Nuclear Science and Technology, Sup. 2, p. 974-977 (August 2001), ND2001.

[10] P. Batistoni, L. Petrizzi, Task T426- Neutronics Experiments, Experimental Validation of Shut Down Dose Rates, EFF-Doc-726, March 2000.

[11] 吳宜燦，李瑩，盧磊，等. 蒙特卡羅粒子輸運計算自動建模程序系統(tǒng)的研究與發(fā)展[J]. 核科學與工程，2006, 26(1):20-27.

[12] 曾勤，盧磊，李瑩，等. 蒙特卡羅粒子輸運計算自動建模程序MCAM在ITER核分析建模中的應用[J]. 原子核物理評論，2006, 23(2):138-141.

[13] Y. Wu, FDS Team. CAD-based Interface Programs for Fusion Neutron Transport Simulation. Fusion Engineering and Design, 2009, 84: 1987-1992.

[14] 吳宜燦，胡麗琴，龍鵬程，等. 先進核能系統(tǒng)設計分析軟件與數(shù)據(jù)庫研發(fā)進展[J]. 核科學與工程，2010, 30(1): 55-64.

[15] J. Zou, Z. He, Q. Zeng, et al. Development and testing of multigroup library with correction of self-shielding effects in fusion-fission hybrid reactor. Fusion Engineering and Design, 2010,85: 1587-1590.

[16] C.Cheng, J. Zou, D. Xu, et al. Development and benchmark of high energy continuous-energy neutron cross section library HENDL-ADS/MC. Nuclear Science and Engineering, 2012, 32(4): 74-79.

[17] M.J. Loughlin. Recommendations on Computation of Dose from Flux Estimates. ICRP 74, 2008.

[18] Y. T. Luo, P. C. Long, G. Y. Wu, et al. SVIP-N 1.0: An integrated visualization platform for neutronics analysis. Fusion Engineering and Design, 2010, 85(7-9): 1527-1530.

[19] P. C. Long, Q. Zeng, T. He, et al. Development of a geometry-coupled visual analysis system for MCNP. Progress in nuclear science and technology, 2010, 2: 280-283.

[20] 龍鵬程, 羅月童, 鄒俊, 等. 基于可編程圖形處理器的可視化技術在中子學分析中的應用研究 [J]. 核電子學與探測技術, 2010, 30(8): 1042-1045.

Rigorous Two-step Shutdown Dose rate Calculation Method Based on Mesh Tally

ZHENG Jian, LI Bin, ZOU Jun, YANG Qi

(Key Laboratory of Neutronics and Radiation Safety, Institute of Nuclear Energy Safety Technology, Chinese Academy of Sciences, Hefei, Anhui, 230031, China)

Based on the computation theory of rigorous two-step, a shutdown dose rate calculation code utilizing mesh tally technique was developed. This code supported cylindrical mesh. It integrated a source subroutine which could provide sources sampling from repeated structures geometry. Additionally, based on the platform of multi-core blade server, multi-node and multi-thread technique were adopted to conduct the mesh material and inventory calculations. The verification test was conducted on ITER-T426 experimental benchmark. The results were in good agreement with the experimental values, which proved the validity and reliability of this method and code.

Mesh tally; Shutdown dose rates; Rigorous two-step method; Automatic coupling

2016-11-21

國家ITER 973計劃(2011GB11306)；中科院戰(zhàn)略性先導科技專項(XDA03040000)；中國科學院知識創(chuàng)新工程重要方向項目(095CF2R211、KJCX2-YW-N35)

鄭 劍(1988—)，男，湖北黃岡，碩士研究生，現(xiàn)主要從事中子學計算研究工作

楊 琪:qi.yang@fds.org.cn

TL63

A

0258-0918(2016)06-0784-06