王 杰，石志勇，宋 婧，郝麗娟，胡麗琴,1，葛 鵬

(1. 中國科學(xué)技術(shù)大學(xué)，安徽合肥230027；2.中國科學(xué)院核能安全技術(shù)研究所，中國科學(xué)院中子輸運理論與輻射安全重點實驗室，安徽合肥230031；3. 火箭軍裝備研究院，北京市100094)

基于CAD的點核劑量計算方法研究及初步應(yīng)用

王 杰1,2，石志勇3，宋 婧2，郝麗娟2，胡麗琴2,1，葛 鵬2

(1. 中國科學(xué)技術(shù)大學(xué)，安徽合肥230027；2.中國科學(xué)院核能安全技術(shù)研究所，中國科學(xué)院中子輸運理論與輻射安全重點實驗室，安徽合肥230031；3. 火箭軍裝備研究院，北京市100094)

點核積分方法是輻射劑量計算的基本方法之一，廣泛應(yīng)用于輻射防護(hù)領(lǐng)域。傳統(tǒng)的點核劑量計算中采用文本方式描述計算模型，存在難以描述復(fù)雜幾何、易出錯且耗時的問題。針對該問題，本文基于FDS團(tuán)隊自主研發(fā)的超級蒙卡核模擬軟件系統(tǒng)SuperMC，進(jìn)行了基于CAD的點核劑量計算方法研究與程序開發(fā)，可基于實際問題的CAD模型直接進(jìn)行支持多重源的光子點核劑量計算，提高了程序?qū)?fù)雜三維幾何問題的處理能力，并包含較為完備的核數(shù)據(jù)庫。使用ANSI/ANS6.6.1、ESIS和VisiPlan的基準(zhǔn)例題對程序進(jìn)行了測試驗證，測試結(jié)果與VisiPlan4.0對比吻合良好。同時將該方法初步應(yīng)用于ITER熱室屏蔽的設(shè)計中，說明了本方法及程序處理復(fù)雜場景問題的能力。

劑量率；點核；CAD；ITER；SuperMC

點核積分法是一種格林函數(shù)積分法，適用于計算和處理γ射線在復(fù)雜幾何空間的輻射屏蔽問題[1]。作為輻射劑量計算的重要方法之一，點核積分法不僅不受空間尺寸和屏蔽體厚度的限制，而且計算速度快，具有較高的計算效率。該方法采用線性衰減方法計算γ射線在物質(zhì)中的衰減過程，考慮到散射通量的影響，引入積累因子來計算空間某一點散射光子所造成的影響的積累或增加[2]。

能量為E(MeV)，源強為1photon/sec的γ點源在空間R處產(chǎn)生的劑量率為:

(1)

式中:F(E)為通量劑量轉(zhuǎn)換系數(shù);B(E，Ld)為劑量積累因子;R為γ點源到劑量計算點處的直線距離;Ld(E)為光學(xué)距離：

(2)

式中：n是γ射線穿過的材料種類數(shù)；μi為第i種材料的線衰減系數(shù)；ti為射線在第i種材料中的穿透距離，光學(xué)距離的單位以平均自由程(mfp)表示。所謂的點核積分方法，就是在輻射屏蔽幾何空間中計算以下積分：

D(r)=∫E?VS·D(E,Ld,r)dEdV

(3)

式中：V為輻射源所在的幾何區(qū)域；E為輻射源的能譜分布；S為輻射源的強度。由于幾何區(qū)域和輻射源的復(fù)雜性，(3)式中的積分一般不能用解析方法求出，而是通過對輻射源進(jìn)行空間和能譜離散使其成為點源，從而實現(xiàn)積分。

目前國內(nèi)外廣泛應(yīng)用于屏蔽設(shè)計的點核積分程序有QAD[3]、PUTZ、VisiPlan[4]等。然而，隨著先進(jìn)反應(yīng)堆[5-7]的發(fā)展，核設(shè)施的幾何結(jié)構(gòu)更加復(fù)雜，空間尺度也更大。QAD、PUTZ等一些采用文本形式描述計算模型的點核積分程序，難以方便地構(gòu)建復(fù)雜的三維幾何模型；而VisiPlan、Microshield等程序雖然增加了3D可視化建模功能，但仍需要對現(xiàn)實場景進(jìn)行很大程度的簡化才能計算，并且這些程序所使用的核數(shù)據(jù)庫并不完備，缺少一些常用的核數(shù)據(jù)，不能很好地滿足屏蔽設(shè)計的需求。為此，本文基于中國科學(xué)院核能安全技術(shù)研究所·FDS團(tuán)隊自主研發(fā)的超級蒙卡核模擬軟件系統(tǒng)SuperMC[8-11]，進(jìn)行了基于CAD點核劑量計算方法的研究和程序開發(fā)。

本文首先介紹了基于CAD點核劑量計算的關(guān)鍵方法，其次通過ANSI/ANS-6.6.1[12]、ESIS[13]和VisiPlan中的基準(zhǔn)例題對程序進(jìn)行了測試，最后將開發(fā)的程序應(yīng)用于ITER熱室屏蔽設(shè)計中。

1 基于CAD的點核劑量計算關(guān)鍵方法

1.1 CAD幾何建模

由(2)式可以看出，光學(xué)距離Ld的計算需要預(yù)先計算出γ射線穿透的材料種類和距離。因此在使用點核積分法進(jìn)行劑量計算之前，需要以構(gòu)造實體幾何法(CSG)形成相應(yīng)的幾何空間，采用射線追蹤技術(shù)計算γ射線在幾何空間中的穿透距離，然而用文本方式描述CSG模型存在著抽象易錯并且耗時的問題。針對該問題，本文在SuperMC已有版本基礎(chǔ)上，基于實際問題的CAD模型直接進(jìn)行點核劑量計算。

基于邊界表示法(BREP)描述的CAD幾何模型可以在SuperMC中自動轉(zhuǎn)換為基于構(gòu)造實體幾何法(CSG)表述的劑量計算模型[14]，算法的實質(zhì)是把復(fù)雜模型分解為簡單形體，然后基于基本的面(半空間)描述簡單形體，用簡單形體的布爾組合來完成對復(fù)雜形體的描述[15-16]。每個幾何區(qū)域都被賦予了材料信息，并轉(zhuǎn)換到CSG模型中用于點核劑量計算。SuperMC中采用層次樹結(jié)構(gòu)來描述幾何和材料信息(組成成分、密度等)，方便了在射線追蹤過程中對幾何體的查找，層次樹的根節(jié)點作為世界體，包含所有的幾何區(qū)域，避免了出現(xiàn)某些區(qū)域未定義的情況。借助于CAD建模工具(CATIA、AutoCAD等)可以描述任意幾何，從而能夠方便的構(gòu)建復(fù)雜的三維幾何模型[17]。

1.2 輻射源

點核積分方法中對于輻射源的描述參數(shù)主要分有能譜、幾何和源強。本文將輻射源的能譜離散為25個能群，能量范圍為15keV～15MeV。每個放射源的能譜信息包含不同的能量值及該能量光子在所有出射光子中所占的比率。對于給定的能譜分布，(1-3)中對于能量的積分轉(zhuǎn)換為：

(4)

式中：i為能群標(biāo)號；Ei為第i個能群的能量值；FEi為Ei能量的光子在所有出射光子中所占的比率。

程序中支持的輻射源類型包括點源、線源(直線、圓圈)、面源(球面、圓柱面、圓盤、矩形面)和體源(球、圓柱、立方體等)。計算輻射源的體積分時，采用了基于統(tǒng)計誤差的蒙卡積分方法：在抽樣的過程中同時計算統(tǒng)計誤差，并判斷是否達(dá)到指定的值，從而自適應(yīng)地確定蒙卡積分所需要的抽樣點數(shù)，從而將線源、面源和體源離散為N個點源。對于均勻分布的源，在抽樣點數(shù)為N的情況下，式(5)中的積分變?yōu)椋?/p>

(5)

統(tǒng)計誤差：

(6)

對于多重源的情況，則依次計算每個源產(chǎn)生的影響，所有輻射源產(chǎn)生的劑量總和作為該點的劑量率。

1.3 數(shù)據(jù)庫

從點核積分公式(1)和公式(2)可以看出，點核劑量計算中主要需要三個參數(shù)：材料線衰減系數(shù)μ、劑量積累因子B和通量劑量轉(zhuǎn)換系數(shù)F。線性衰減系數(shù)μ可表示為μm/ρ，μm為質(zhì)量衰減系數(shù)。因此，進(jìn)行點核劑量計算需要μm、B和F三個數(shù)據(jù)庫。程序中內(nèi)置了15keV～15MeV范圍內(nèi)25個能群的光子通量劑量轉(zhuǎn)換系數(shù)F。

1.3.1 質(zhì)量衰減系數(shù)

質(zhì)量衰減系數(shù)表征了射線在材料中與物質(zhì)的相互作用，對于γ射線，主要為光電效應(yīng)、康普頓散射和電子對效應(yīng)。本文在丹麥工業(yè)大學(xué)開發(fā)的WinXCom[18]程序和美國國家標(biāo)準(zhǔn)ANSI/ANS-6.4.3[19]的基礎(chǔ)上設(shè)計制作了衰減系數(shù)數(shù)據(jù)庫，包含了1-100號元素和水、空氣、標(biāo)準(zhǔn)混凝土(NBS混凝土)的質(zhì)量衰減系數(shù)。對于水、空氣、NBS混凝土之外的混合物和化合物的質(zhì)量衰減系數(shù)，本文采用下式計算：

(7)

式中：ωi為元素i的質(zhì)量分?jǐn)?shù)；μmi為元素i的質(zhì)量衰減系數(shù)。程序?qū)⒏鶕?jù)材料的組成成分，自動計算該材料的質(zhì)量衰減系數(shù)。

1.3.2 積累因子

積累因子是進(jìn)行γ射線屏蔽分析時所使用的一切點核公式的基礎(chǔ)。常用的積累因子分為劑量積累因子和能量吸收積累因子，劑量積累因子又稱為照射量積累因子，是劑量計算點實際劑量與未經(jīng)碰撞粒子產(chǎn)生的劑量的比值[2]。常用的積累因子計算公式有G-P公式、泰勒公式、凱波多項式等。本文采用泰勒公式計算γ射線的劑量積累因子:

B(E,Ld)=A1e-α1Ld+(1-A1)e-α2Ld

(8)

式中參數(shù)A1、α1、α2是通過在矩方法解得的各向同性點源在無線均勻介質(zhì)中輸運結(jié)果的基礎(chǔ)上擬合得到的，與材料和射線的能量相關(guān)。本文所使用的積累因子數(shù)據(jù)庫包含了15種參考材料的泰勒公式參數(shù)，數(shù)據(jù)來源于ANSI/ANS-6.4.3[19]。對于不同的劑量計算點，程序支持選取不同的積累因子，一次計算中可選擇多種積累因子。

2 基準(zhǔn)驗證

為了驗證程序的正確性，本文采用ANSI[12]、ESIS[13]和VisiPlan的相關(guān)例題對程序進(jìn)行了測試，測試結(jié)果與VisiPlan4.0[4]進(jìn)行了對比。

2.1 ANSI/ANS例題測試

本文采用了ANSI/ANS 6.6.1[12]中的點源例題I.1、I.2和體源例題II.1、II.2對程序進(jìn)行了測試，測試結(jié)果如圖1所示。點源例題計算結(jié)果與VisiPlan4.0一致，體源例題所有計算點的計算結(jié)果與VisiPlan相比，偏差均在4.0%以內(nèi)。

圖1 ANSI/ANS例題測試結(jié)果Fig.1 The test results of ANSI/ANS problem

2.2 ESIS例題測試

本文采用ESIS發(fā)布的屏蔽例題對程序進(jìn)行測試驗證，其CAD模型如圖2所示，詳見文獻(xiàn)[13]。測試結(jié)果如表1所示，和VisiPlan的相對偏差分別為0.87%和1.30%。

圖2 ESIS例題的CAD模型Fig.2 The CAD model of ESIS problem表1 ESIS測試結(jié)果Table 1 The test results of ESIS

VisiPlan/(mSv/h)SuperMC/(mSv/h)deviationD149102×102486716×102-087%D283205×10-3842903×10-3130%

2.3 VisiPlan應(yīng)用例題測試

本文采用了VisiPlan中的蒸汽發(fā)生器例題進(jìn)行了測試，如圖3。一回路和二回路中充滿水，一回路傳熱管束簡化為由鐵和水的混合物組成的圓柱體，放射源在其中均勻分布。部分計算結(jié)果如圖4所示，與VisiPlan相比，有兩個點的相對偏差達(dá)到8%，其他均保持在5%以內(nèi)。

圖3 蒸汽發(fā)生器CAD模型Fig.3 The CAD model of steam generator

圖4 蒸汽發(fā)生器例題計算結(jié)果Fig.4 The test results of steam generator

通過ANSI、ESIS和VisiPlan中基準(zhǔn)例題的測試，本文計算結(jié)果與VisiPlan相比最大偏差為8%，出現(xiàn)偏差的主要原因是由于兩個程序采用了不同的質(zhì)量衰減系數(shù)數(shù)據(jù)庫，并且對積分點的抽取方式不同而導(dǎo)致的。

3 初步應(yīng)用：ITER熱室屏蔽設(shè)計

為了說明程序?qū)?fù)雜場景問題的處理能力，本文對ITER熱室屏蔽問題進(jìn)行了初步分析，給出了熱室中門7的設(shè)計方案。

ITER熱室建筑位于托卡馬克大廳[20]北邊，如圖5(a)所示。熱室中存在7個混凝土(2.6g/cm3)門，位置如圖5(b)所示，當(dāng)4個被活化的第一壁組件被運輸?shù)綗崾抑泻螅箝T2、3、4、7后的劑量應(yīng)當(dāng)?shù)陀?0μSv/h，門5、6后的劑量應(yīng)在1000μSv/h以下[21]。為此，需要計算出各個門的最佳厚度。

本文僅以門7為例，采用開發(fā)的程序?qū)﹂T7的厚度進(jìn)行初步設(shè)計。為了簡化計算過程，把每個第一壁組件等效為一個平板體源和圓柱體源的組合體。采用等效后的輻射源，可以計算出在門7前的γ劑量率為0.085Sv/h，與采用原始模型MCNP計算的結(jié)果(0.08Sv/h)符合很好，因此這樣的簡化是可行的。

圖5 ITER HCB幾何模型Fig.5 The geometry of ITER HCB (a) 熱室CAD模型;(b) 熱室中7個門的位置

在等效源的基礎(chǔ)上，本文計算了在不同厚度下，門7后的劑量率。隨著厚度的增加，γ劑量率的衰減如圖6所示。

圖6 不同厚度的門7對γ劑量率的衰減Fig.6 The attenuation of dose rate with different thickness

混凝土門的厚度X與衰減后的γ劑量率D可用下式表示：

(9)

b為常數(shù)，表征γ射線在混凝土中的衰減效果。根據(jù)圖6，可以計算出b≈0.0585。

因此，為了使門7后的劑量率小于10μSv/h，門7的厚度至少為67cm。表1給出了在不同安全系數(shù)下，門7所需的厚度以及使用MCNP設(shè)計出的厚度。

表1 不同安全系數(shù)所要求的門7厚度Table 1 Minimum thickness of door 7

從表中看出，本文計算結(jié)果比MCNP結(jié)果大了5～7cm，這是因為點核積分法使用的是無限介質(zhì)積累因子，包含了在邊界處的反射貢獻(xiàn)，導(dǎo)致計算結(jié)果偏高，因此本文得出的是一個偏保守的設(shè)計方案，能夠確保實際的劑量率滿足劑量限制。在計算時間上，使用MCNP進(jìn)行一次計算需要數(shù)周的時間，而采用本文所提出的方法僅需要幾個小時，大大縮減了計算時間。

4 結(jié)論

本文基于實際問題的CAD模型直接進(jìn)行點核劑量計算，提高了程序?qū)?fù)雜三維幾何問題的處理能力，可進(jìn)行支持多重源的光子劑量計算，并包含較為完備的核數(shù)據(jù)庫。采用ANSI/ANS-6.6.1、ESIS和VisiPlan中的基準(zhǔn)例題對程序進(jìn)行了測試驗證，計算結(jié)果與VisiPlan4.0的計算結(jié)果吻合很好，證明了本文方法和程序的正確性。將開發(fā)的程序初步應(yīng)用于ITER熱室屏蔽問題中，給出了對于門7的初步設(shè)計方案，顯示程序處理復(fù)雜場景問題的能力。在后續(xù)的研究中，將會結(jié)合實際反應(yīng)堆應(yīng)用場景開展應(yīng)用。

致謝

本工作得到中科院核能安全技術(shù)研究所·FDS團(tuán)隊各位老師的幫助與指導(dǎo)，在此深表感謝。

[1] 萬海霞. 輻射場可視化平臺劑量分布快速計算程序開發(fā)[D]. 華北電力大學(xué)，2012.

[2] Schaeffer N M，Reactor Shielding for Nuclear Engineer[R]. Radiation Research Associates，Inc. Fort Worth，Tex. (USA),1973.

[3] Cain V R. QAD-CG: A Combinatorial Geometry Version of QAD-P5A Point-Kernel Code for Neutron and Gamma-ray Shielding Calculation[Z]. 1977，ORNL-CCC307.

[4] Vermeersch F，Van Bosstraeten C. Development of the VisiPlan ALARA planning tool[J]. Proceeding of the International Conference on Topical issues in Nuclear Radiation and Radioactive Waste Safety，Vienna Austria，August 31 to September 4，1998.

[5] Y. Wu，H. chen，S. Liu，et al. Fusion-Based Hydrogen Production Reactor and Its Material Selection[J]. Journal of Nuclear Materials，2009，386-388:122-126.

[6] Y. Wu，S. Zheng，X. Zhu，et al. Conceptual Design of the Fusion-driven Subcritical System FDS-I[J]. Fusion Engineering and Design，2006，81(8): 1305-1311.

[7] Y. Wu，F(xiàn)DS Team. Conceptual Design of the China Fusion Power Plant FDS-II[J]. Fusion Engineering and Design，2008，83: 1683-1689.

[8] J. Song，G. Sun，Z. Chen，et al. Benchmarking of CAD-based SuperMC with ITER Benchmark Model[J]. Fusion Engineering and Design，DOI: 10.1016/j. fusengdes.2014.05.003.

[9] 孫光耀，宋婧，鄭華慶，等. 超級蒙特卡羅計算軟件SuperMC2.0中子輸運計算校驗[J]. 原子能科學(xué)與技術(shù)，2013，第47卷增刊: 140-145.

[10] Y. Wu，J. Song and FDS Team. Development of Super Monte Carlo Calculation Program SuperMC 2.0. Proc. Int. Conf. ANS National Meeting-2013 ANS Winter Meeting and Technology Expo，American Nuclear Society，Washington D. C.，USA，November 10-14，2013.

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

[12] Calculation and Measurement of Direct and Scattered Gamma Radiation from LWR Nuclear Power Plants，ANSI/ANS-6.6.1-1991，Amer. Nat. Standards Int.，New York，1991.

[13] Specification for Gamma Ray Shield Benchmark Applicable to A Nuclear Radioactive Waste Facility. ESIS Newsletter 5 (April 1973).

[14] Y. Wu，J. Song，H. Zheng，et al. CAD-Based Monte Carlo Program for Integrated Simulation of Nuclear System SuperMC[J]. Annals of Nuclear Energy，2015，82:162-168.

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

[16] Y. Wu，F(xiàn)DS Team. CAD-based Interface Programs for Fusion Neutron Transport Simulation[J]. Fusion Engineering and Design，2009，84(7):1987-1992.

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

[18] Gerward L，Guilber N，Jensen K B，et al. WinXCom-A Program for Calculation X-Ray Attenuation Cofficients[J]. Radiation Physics and Chemistry，2004，7(3): 653-654.

[19] Gamma-Ray Attenuation Cofficient and Buildup Factor for Engineering Materials，ANSI/ANS-6.4.3-1991，Amer. Nat. Standards Int.，New York，1991.

[20] Y. Li，L. Lu，A. Ding，et al. Benchmarking of MCAM 4.0 with the ITER 3D Model[J].Fusion Engineering and Design，2007，82(15): 2861-2866.

[21] Q. Yang. Dose Rates in Port Cells and Shielding of Port Cell Door，Nov. 2013，ITER_D_ LY6K73 v1.0.

Research and Preliminary Application of CAD-Based Point-KernelIntegral Method for Dose Calculation

WANG Jie1，2，SHI Zhi-yong3，SONG Jing2，HAO Li-juan2，HU Li-qin2，1，GE Peng2

(1. University of Science and Technology of China，Hefei，Anhui，230027，China;2. Key Laboratory of Neutronics and Radiation Safety，Institute of Nuclear Energy Safety Technology，Chinese Academy of Sciences，Hefei，Anhui，230031，China;3. EquipmentAcademy of the Rocket Force)

Point-kernel integral method is widely used in radiation shielding field and taken as a basic method for radiation dose calculation. However，for traditional point-kernel codes，describing the geometry of problems using text files is difficult to describe complex geometry and error-prone. Based on Super Monte Carlo Simulation Program for Nuclear and Radiation Process SuperMC，developed by FDS Team，a CAD-based point-kernel integral method was proposed in this study and implemented. The CAD model of problems can be directly used for calculation，which improves the ability of handing complex geometry. Gamma dose calculation can be performed with multi-source and complete database based on point-kernel method. Benchmark cases from ANSI/ANS-6.6.1，ESIS and VisiPlan were used to validate the code. The results were in consistent with those of VisiPlan4.0. Using SuperMC，the preliminary design of shielding design of ITER hot cell building was performed，which demonstrated the ability of this method for handling complex problems.

Dose Rate; Point-Kernel; CAD; ITER; SuperMC

2016-07-21

國家自然科學(xué)基金91026004；國家ITER 973計劃2014GB112001

王 杰(1991—)，男，安徽人，碩士研究生，主要從事點核劑量計算方法研究工作

葛 鵬:peng.ge@fds.org.cn

329+.2

A

0258-0918(2016)05-0656-07