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

鄭 劍，李 斌，鄒 俊，楊 琪

(中國(guó)科學(xué)院核能安全技術(shù)研究所，中國(guó)科學(xué)院中子輸運(yùn)理論與輻射安全重點(diǎn)實(shí)驗(yàn)室，合肥230031)

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

鄭 劍，李 斌，鄒 俊，楊 琪

(中國(guó)科學(xué)院核能安全技術(shù)研究所，中國(guó)科學(xué)院中子輸運(yùn)理論與輻射安全重點(diǎn)實(shí)驗(yàn)室，合肥230031)

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

網(wǎng)格計(jì)數(shù)；停機(jī)劑量率；嚴(yán)格二步法；自動(dòng)耦合

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

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

1 計(jì)算方法

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

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

圖1 計(jì)算流程圖Fig.1 Calculation flow chart

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

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

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

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

利用比濁法的原理，用光密度值(OD值)的變化來判斷菌株存活情況。根據(jù)朗伯比爾定律OD=log It IO=kbc，c表示樣品濃度，若樣品液厚度b一定，則OD值與樣品的濁度相關(guān)[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之間)，在波長(zhǎng)600 nm處測(cè)量吸光度值，每組試驗(yàn)做3組平行。

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

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

(1)

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

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

1.3 源子程序

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

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

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

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

(2)

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

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

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

(3)

(4)

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

1.4 并行計(jì)算

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

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

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

2 ITER-T426實(shí)驗(yàn)例題測(cè)試

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

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

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

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

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

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

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

圖8 計(jì)算值與實(shí)驗(yàn)值的比值Fig.8 Comparison between calculated and experiment dose rate in the cavity center.

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

3 結(jié)束語

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

為了驗(yàn)證方法和程序的正確性，利用國(guó)際熱核聚變實(shí)驗(yàn)堆停堆劑量基準(zhǔn)實(shí)驗(yàn)ITER-T426實(shí)驗(yàn)進(jìn)行測(cè)試，計(jì)算結(jié)果證明了該方法和程序的正確性以及處理復(fù)雜幾何裝置的適用性。

致謝

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

[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] 吳宜燦，李靜驚，李瑩，等. 大型集成多功能中子學(xué)計(jì)算與分析系統(tǒng)VisualBUS的研究與發(fā)展[J]. 核科學(xué)與工程，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] 吳亮亮，應(yīng)棟川，邱岳峰，等. 三維停堆劑量率計(jì)算程序研發(fā)及其在EAST上的應(yīng)用[J]. 核科學(xué)與工程，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] 吳宜燦，李瑩，盧磊，等. 蒙特卡羅粒子輸運(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] Y. Wu, FDS Team. CAD-based Interface Programs for Fusion Neutron Transport Simulation. Fusion Engineering and Design, 2009, 84: 1987-1992.

[14] 吳宜燦，胡麗琴，龍鵬程，等. 先進(jìn)核能系統(tǒng)設(shè)計(jì)分析軟件與數(shù)據(jù)庫研發(fā)進(jìn)展[J]. 核科學(xué)與工程，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] 龍鵬程, 羅月童, 鄒俊, 等. 基于可編程圖形處理器的可視化技術(shù)在中子學(xué)分析中的應(yīng)用研究 [J]. 核電子學(xué)與探測(cè)技術(shù), 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

國(guó)家ITER 973計(jì)劃(2011GB11306)；中科院戰(zhàn)略性先導(dǎo)科技專項(xiàng)(XDA03040000)；中國(guó)科學(xué)院知識(shí)創(chuàng)新工程重要方向項(xiàng)目(095CF2R211、KJCX2-YW-N35)

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

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

TL63

A

0258-0918(2016)06-0784-06