基于網格計數的停機劑量率計算方法研究

鄭 劍，李 斌，鄒 俊，楊 琪

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

基于網格計數的停機劑量率計算方法研究

鄭 劍，李 斌，鄒 俊，楊 琪

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

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

網格計數；停機劑量率；嚴格二步法；自動耦合

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

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

1 計算方法

1.1 總體設計

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

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

該方法首先使用MCNP程序[6]，對已劃分了網格的計算模型進行網格計數中子輸運計算，獲得所有被計數網格的中子分群通量。然后調用MCNP程序對每一個網格進行PTRAC計算，獲得每個網格的粒子分布情況，用于計算網格的平均材料。

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

最后，使用MCNP程序對不同時間段下的衰變光子進行網格計數輸運計算，最終獲得不同時間段下計數網格的劑量率分布。

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

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

1.2 網格材料計算

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

(1)

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

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

1.3 源子程序

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

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

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

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

(2)

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

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

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

(3)

(4)

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

1.4 并行計算

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

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

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

2 ITER-T426實驗例題測試

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

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

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

在網格計數停機劑量率計算中，對于中子輸運計算，本文使用HENDL3.0[15-16]數據庫計算計數網格的中子通量分布。FISPACT程序活化計算采用活化庫EAF-2007[7]。中子及衰變光子輸運計算使用相同的計算模型。本文在圓柱坐標系下在模型R、θ、Z方向分別劃分43×8×66計數網格。在衰變光子的輸運計算中，本文使用衰變光子取樣程序進行衰變光子源抽樣。在光子輸運計算時，使用ITER推薦的通量劑量轉換因子[16]計算出模型在不同停機時間的劑量率。計算結果的可視化分析使用FDS團隊自主開發的科學計算可視化分析軟件SVIP[18-20]。圖5是該裝置的中子通量分布的剖面圖；圖6顯示了停堆12天后的光子源分布情況，可以看出右側的不銹鋼部分活化比較嚴重；圖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給出了網格計數停機劑量率程序計算結果與實驗測量結果隨冷卻時間的變化。從該圖中可以看出，當冷卻時間比較短時，程序計算結果比實驗結果最大處約低估25%。這是因為計算使用的輻照方案是簡化的，不能夠完全模擬實際的實驗情況，并且有些對劑量貢獻較多的核素沒有被考慮。而冷卻時間比較長時，程序計算結果最高比實驗結果高出10%，此時的誤差的主要來自于網格劃分不夠細。本程序計算結果與實驗測量結果以及其他程序計算結果趨勢相吻合。

3 結束語

本文根據嚴格二步法的計算理論，研究了基于網格計數的停機劑量率計算方法，設計并實現了基于網格計數的停機劑量率計算程序。

為了驗證方法和程序的正確性，利用國際熱核聚變實驗堆停堆劑量基準實驗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] 吳宜燦，李靜驚，李瑩，等. 大型集成多功能中子學計算與分析系統VisualBUS的研究與發展[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] 吳亮亮，應棟川，邱岳峰，等. 三維停堆劑量率計算程序研發及其在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] 吳宜燦，李瑩，盧磊，等. 蒙特卡羅粒子輸運計算自動建模程序系統的研究與發展[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] 吳宜燦，胡麗琴，龍鵬程，等. 先進核能系統設計分析軟件與數據庫研發進展[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)；中科院戰略性先導科技專項(XDA03040000)；中國科學院知識創新工程重要方向項目(095CF2R211、KJCX2-YW-N35)

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

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

TL63

A

0258-0918(2016)06-0784-06