研究堆退役中的快速輻射傳輸計(jì)算方法研究

畢遠(yuǎn)杰，楊宏偉，郭金森，駱志平，李傳龍，陳 凌

（中國原子能科學(xué)研究院 輻射安全研究所，北京 102413）

研究堆退役是核燃料循環(huán)的重要環(huán)節(jié)之一，目前我國的研究堆已逐步進(jìn)入退役階段。然而，現(xiàn)有的退役工作通常是基于初步的源項(xiàng)分析和輻射監(jiān)測(cè)數(shù)據(jù)進(jìn)行的。這樣的操作流程存在的問題是：隨著設(shè)備的拆除和去污工作的進(jìn)行，輻射源項(xiàng)和屏蔽條件均在不停地發(fā)生變化，而輻射監(jiān)測(cè)數(shù)據(jù)不能及時(shí)地跟進(jìn)。同時(shí)，在工作過程中，即使選取了參考點(diǎn)進(jìn)行輻射劑量的測(cè)量，由于輻射場分布的差別，參考點(diǎn)也不一定能反映人員所在位置的實(shí)際劑量。因此，操作人員可能在實(shí)際的工作過程中受到額外的照射，對(duì)他們的身體造成傷害。

隨著計(jì)算機(jī)技術(shù)的發(fā)展，目前在大型核設(shè)施的退役過程中，人們更希望能通過計(jì)算機(jī)模擬給出整個(gè)廠房內(nèi)部的三維劑量場分布及其動(dòng)態(tài)變化情況，為工作人員在放射性工作場所內(nèi)進(jìn)行放射性操作時(shí)的工作方案設(shè)計(jì)提供保障。根據(jù)計(jì)算機(jī)的模擬結(jié)果，工作人員可對(duì)不同的退役方案進(jìn)行評(píng)估，選擇最優(yōu)化的方案進(jìn)行實(shí)施。因而，進(jìn)行三維輻射場快速計(jì)算成為時(shí)代的需求。

在進(jìn)行三維輻射場快速計(jì)算的過程中存在兩個(gè)問題：1）空間，需進(jìn)行整個(gè)廠房內(nèi)全空間的輻射場計(jì)算；2）時(shí)間，伴隨著設(shè)備的拆卸需在盡量短的時(shí)間內(nèi)給出新的劑量分布情況。目前傳統(tǒng)的計(jì)算方法很難同時(shí)解決這兩個(gè)問題，因此，本文擬制定一種新的計(jì)算方案，可在較短的時(shí)間內(nèi)給出全空間的輻射場分布情況。

1 輻射場計(jì)算方法現(xiàn)狀

1.1 傳統(tǒng)計(jì)算方案

解決輻射傳輸問題的基本工具是Boltzmann方程［1-2］。求解Boltzmann 方程的方法分為兩類：1）確定論方法，目前最常用的是離散縱標(biāo)方法，代 表 性 的 程 序 有ANISN［3］、DORT［4］、TORT［5］等；2）蒙特卡羅方法（簡稱蒙卡方法），蒙卡方法是一種隨機(jī)性方法，代表性的程序有MCNP［6］、FLUKA［7-8］、GEANT4［9-10］、MARS［11］、PHITS［12］等。在實(shí)際的工程應(yīng)用中，還常用到Boltzmann方程的近似解法，如在屏蔽計(jì)算中用點(diǎn)核方法進(jìn)行估算，常用的點(diǎn)核程序有MICROSHIELD［13］、QAD［14］、MERCURE［15］等。

點(diǎn)核方法的主要優(yōu)點(diǎn)是快速。除廣泛用于屏蔽設(shè)計(jì)的估算外，近年來新興的用于輻射防護(hù)最優(yōu)化（ALARA）的輻射場快速計(jì)算軟件也多基于點(diǎn)核方法編寫。國際上已有一些相關(guān)軟件，如比利時(shí)核研究中心開發(fā)的VISIPLAN 軟件［16］；法 國 的CEA-LIST 在 點(diǎn) 核 程 序MERCURE 5的基礎(chǔ)上開發(fā)的NARVEOS軟件［17］；日本核燃料循環(huán)開發(fā)署與日本原子能研究機(jī)構(gòu)、挪威OECD／NEA Halden的反應(yīng)堆項(xiàng)目聯(lián)合開發(fā)的退役工程支持系統(tǒng)DEXUS［18］；美國電力研究協(xié)會(huì)EPRI 開發(fā)的三維可視化系統(tǒng)［19］等。國內(nèi)方面如中國科學(xué)院核能安全技術(shù)研究所FDS團(tuán)隊(duì)基于國家大科學(xué)工程全超導(dǎo)托卡馬克核聚變實(shí)驗(yàn)裝置EAST 建立的輻射虛擬仿真系統(tǒng)RVIS［20］，也將建立基于點(diǎn)核方法的快速計(jì)算模塊。以上基于點(diǎn)核方法的虛擬仿真軟件雖快速，但存在源項(xiàng)周圍劑量估算不準(zhǔn)確、不能處理復(fù)雜的幾何結(jié)構(gòu)、可靠性低等問題。

1.2 輻射場耦合計(jì)算

蒙卡、離散縱標(biāo)和點(diǎn)核方法各有優(yōu)缺點(diǎn)（表1）。單一的計(jì)算方法均存在各自的不足，蒙卡方法對(duì)于大空間厚屏蔽問題無法在合理的時(shí)間內(nèi)給出可靠的計(jì)算結(jié)果，確定論方法則不適合復(fù)雜源項(xiàng)和幾何區(qū)域。耦合計(jì)算則可將整個(gè)問題劃分為兩個(gè)區(qū)域進(jìn)行求解，在復(fù)雜源項(xiàng)和幾何區(qū)域采用蒙卡方法進(jìn)行模擬，在大空間厚屏蔽區(qū)域采用確定論方法求解深穿透問題。

1）蒙卡-離散縱標(biāo)耦合方法

目前存在的蒙卡-離散縱標(biāo)耦合方法的研究，主要是針對(duì)MCNP和DOORS代碼包中程序的耦合。如MCNP 與一維離散縱標(biāo)程序ANISN 的 耦 合［21-23］，MCNP 與 二 維 離 散 縱 標(biāo)程 序DORT 的 耦 合［24-25］，MCNP 與 三 維 離 散縱標(biāo)程序TORT 的耦合［26-27］。蒙卡-離散縱標(biāo)耦合的特色在于進(jìn)行全空間的趨于精確的計(jì)算。但在離散縱標(biāo)方法中，考慮到能量多群、方位角離散和空間網(wǎng)格劃分，得到的是一系列高維積分微分方程組，需各種技巧迭代至收斂。尤其是三維離散縱標(biāo)程序，計(jì)算量龐大。同時(shí)，在蒙卡和離散縱標(biāo)的耦合過程中，若連接面上通量的空間分布和方位角分布缺乏對(duì)稱性，則需調(diào)整方位角離散和空間網(wǎng)格劃分的大小來保證離散縱標(biāo)方法求解的精度，這會(huì)大幅增加計(jì)算量。所以，蒙卡-離散縱標(biāo)耦合方法并不適用于快速計(jì)算。

表1 3種計(jì)算方法的優(yōu)缺點(diǎn)Table 1 Comparison of advantage and disadvantage of three methods

2）蒙卡-點(diǎn)核耦合方法

蒙卡與點(diǎn)核耦合進(jìn)行全空間輻射場快速計(jì)算還未見詳細(xì)報(bào)道，蒙卡-點(diǎn)核耦合方法的特點(diǎn)為：（1）可給出全空間的三維輻射場，解決點(diǎn)核方法在復(fù)雜源項(xiàng)和幾何區(qū)域存在的問題；（2）相對(duì)于蒙卡-離散縱標(biāo)耦合方法的趨于精確計(jì)算，蒙卡-點(diǎn)核耦合方法更加快速，可隨著設(shè)備的拆卸在盡量短的時(shí)間內(nèi)給出新的劑量場分布。

本文基于蒙卡-點(diǎn)核耦合的計(jì)算方案，編寫耦合接口模塊，用于輻射場的三維快速計(jì)算。

2 耦合計(jì)算系統(tǒng)

圖1為耦合計(jì)算程序系統(tǒng)示意圖。在耦合計(jì)算中，將整個(gè)問題劃分為兩個(gè)區(qū)域：復(fù)雜源項(xiàng)和幾何區(qū)域及簡單幾何大屏蔽區(qū)域。復(fù)雜源項(xiàng)和幾何區(qū)域的結(jié)構(gòu)通常較復(fù)雜，所以在這種區(qū)域采用蒙卡方法進(jìn)行粒子輸運(yùn)模擬。而在大屏蔽部分，幾何通常較簡單，同時(shí)屬于粒子深穿透問題，采用點(diǎn)核方法進(jìn)行計(jì)算。結(jié)合蒙卡程序和點(diǎn)核計(jì)算的結(jié)果，得到三維輻射場快速計(jì)算軟件，并對(duì)軟件進(jìn)行基準(zhǔn)校驗(yàn)和對(duì)比分析。

圖1 耦合計(jì)算程序系統(tǒng)示意圖Fig.1 Schematic of coupled calculation program

耦合接口模塊的編制是計(jì)算程序系統(tǒng)的重點(diǎn)，耦合接口模塊需解決的主要問題是如何將蒙卡輸出的粒子軌跡轉(zhuǎn)換為點(diǎn)核計(jì)算的源項(xiàng)。首先需從物理上建立兩者之間的聯(lián)系，即獲得等效面源。通過分析Boltzmann方程可得出：在一個(gè)封閉曲面內(nèi)的源和物質(zhì)，可由這個(gè)封閉曲面上的粒子流量率來代替，這兩種情況對(duì)探測(cè)器的響應(yīng)相同［1］。在這個(gè)等效問題中，通過連接曲面的流量率相當(dāng)于等效面源的源強(qiáng)。

在計(jì)算過程中，首先對(duì)交界面進(jìn)行空間上的網(wǎng)格劃分得到面元dA，接著對(duì)dA 出射的粒子進(jìn)行方位角離散和能量分群，最終得到（dA，dE，dΩ）。設(shè)dA 的中心點(diǎn)為A′、外法向?yàn)閚A，則位于P 點(diǎn)的探測(cè)器和dA 的外法線方向的夾角，因 此，dA 對(duì)P點(diǎn)探測(cè)器的貢獻(xiàn)為：

其中：d（θ）為沿θ方向穿過屏蔽的厚度；r為點(diǎn)到探測(cè)器的距離；Hθ為無屏蔽的情況下在θ 方向與源單位距離處的劑量率；μ 為γ射線的線衰減系數(shù)；B（d）為累積因子。

3 退役算例

本文對(duì)退役過程中放射性廢物儲(chǔ)存容器所在的大尺度厚屏蔽房間的輻射場進(jìn)行了計(jì)算，計(jì)算算例幾何和參考點(diǎn)如圖2所示。圓柱形放射性廢物儲(chǔ)存罐位于x＝-290cm 處，內(nèi)部裝有含放射性核素137Cs的淤泥，容器壁為1cm厚不銹鋼。放射性核素濃度在z方向呈指數(shù)分布，抽樣函數(shù)為：其中：淤泥在z 方向的范圍為從z1到z2；ξ 為［0，1］上均勻分布的隨機(jī)變量；q為系數(shù)。

圖2 計(jì)算算例幾何和參考點(diǎn)Fig.2 Geometry of model and reference points

在x＝-683cm～x＝-663cm 放置一短混凝土墻，在x＝-1 283cm～x＝-1 263cm放置一長混凝土墻。

計(jì)算結(jié)果采用通用蒙卡軟件FLUKA 進(jìn)行校準(zhǔn)，表2比較了FLUKA 和耦合程序的計(jì)算結(jié)果。對(duì)于選擇的參考點(diǎn)，兩個(gè)計(jì)算程序的相對(duì)偏差小于10%。在計(jì)算時(shí)間方面，將整個(gè)房間劃分為105個(gè)網(wǎng)格進(jìn)行三維輻射場計(jì)算，在3.4GHz、4GB內(nèi)存的個(gè)人計(jì)算機(jī)上進(jìn)行測(cè)試，F(xiàn)LUKA 達(dá)到表2 中的計(jì)算誤差所需的計(jì)算時(shí)間為12h，耦合程序用20min得到相對(duì)誤差1%的等效面源，18min進(jìn)行點(diǎn)核計(jì)算。

表2 參考點(diǎn)周圍劑量當(dāng)量率比對(duì)Table 2 Comparison of ambient dose equivalent in reference points

4 結(jié)論

本文基于蒙卡與點(diǎn)核耦合的輻射場計(jì)算方法，編制了相應(yīng)的耦合計(jì)算程序，用于研究堆退役過程中的輻射場快速計(jì)算。蒙卡-點(diǎn)核耦合方法結(jié)合了蒙卡和點(diǎn)核計(jì)算各自的優(yōu)點(diǎn)，不僅適合于復(fù)雜源項(xiàng)和屏蔽體區(qū)域的計(jì)算，同時(shí)也適合于大屏蔽區(qū)域的輻射場快速計(jì)算。作為算例，本文采用編制的耦合程序?qū)ρ芯慷淹艘圻^程中放射性廢物儲(chǔ)存房間的輻射場進(jìn)行了計(jì)算，計(jì)算結(jié)果驗(yàn)證了程序的可靠性和適用性。

［1］ CHILTON B，SHULTIS J K，F(xiàn)AW R E.Principles of radiation shielding［M］.USA：Prentice-Hall，Inc.，1984.

［2］ SHULTIS J K，F(xiàn)AW R E.Radiation shielding［M］.USA：American Nuclear Society，Inc.，2000.

［3］ ENGLE W W，Jr.ANISN，a one-dimensional discrete ordinates transport code with anisotropic scattering，K-1693［R］.USA：ORNL，1967.

［4］ RHOADES W A，CHILDS R L.The DORT two-dimensional discrete ordinates transport code［J］.Nuclear Science and Engineering，1988，99：88-89.

［5］ AMY Y Y.The three-dimensional，discrete ordinates neutral particle transport code TORT：An overview［C］∥OECD／NEA Meeting on 3D Deterministic Radiation Transport Computer Programs，F(xiàn)eature，Applications and Perspectives.Paris：［s.n.］，1996.

［6］ HENDRICKS J S，MCKINNEY G W，WATERS L S，et al.MCNPX user's manual：Ver-sion 2.5.0，LAUR-05-2675［R］.USA：LANL，2005.

［7］ BOHLEN T T，CERUTTI F，CHIN M P W，et al.The FLUKA code：Developments and challenges for high energy and medical applications［J］.Nuclear Data Sheets，2014，120：211-214.

［8］ FASSO A，F(xiàn)ERRARI A，RANFT J，et al.FLUKA： A multi-particle transport code，SLAC-R-773［R］.［S.l.］：［s.n.］，2005.

［9］ AGOSTINELLI S，ALLISON J，AMAKO K，et al.GEANT4：A simulation toolkit［J］.Nucl Instrum Phys Methods Res A，2003，506（3）：250-303.

［10］ALLISON J，AMAKO K，APOSTOLAKIS J，et al.GEANT4 developments and applications［J］.IEEE Trans Nucl Sci，2006，53（1）：270.

［11］MOKHOV N V，STRIGANOV S I.MARS15 overview［C］.［S.l.］：［s.n.］，2007.

［12］NIITA K，MATSUDA N，IWAMOTO Y，et al.PHITS overview［C］.［S.l.］：［s.n.］，2007.

［13］MICROSHIELD?user's manual［M］.［S.l.］：Grove Software，INC.，2009.

［14］CAIN V R.QAD-CG，a combinatorial geometry version of QAD-P5A，apoint kernel code for neutron and gamma-ray shielding calculations，ORNL-CCC307［R］.USA：ORNL，1977.

［15］DUPONT C，NIMAL J C.MERCURE-4：A three-dimensional Monte Carlo program for the integration of specific point attenuation kernels in a straight line，OLS-82-116［R］.［S.l.］：［s.n.］，1980.

［16］VERMEERSCH F.Dose assessment and dose optimization in decommissioning using the VISIPLAN 3D ALARA planning tool［C］∥Radiation Protection and Decommissioning Conference ABR／BVS.Brussels：［s.n.］，2003.

［17］THEVENON J B.Narveos：A new tool supporting ALARA studies［R］.French：Marcoule Cea Center，2011.

［18］IGUCHI Y，KANEHIRA Y，TACHIBANA M，et al.Development of decommissioning engineer-ing support system （DEXUS）of the Fugen Nuclear Power Station［J］.J Nucl Sci Technol，2004，41（3）：367-375.

［19］Work planning and 3D visualization tool［R］.USA：Electric Power Research Institute，2009.

［20］TANG Z D，LONG P C，HUANG S Q，et al.Real-time dose assessment and visualization of radiation field for EAST Tokamak［J］.Fusion Engineering and Design，2010，85（7-9）：1 591-1 599.

［21］GALLMEIER F X，PEVEY R E.Creation of a set of interface utilities to allow coupled Monte Carlo／discrete ordinate shielding analysis［C］∥The Third International Topical Meeting on Nuclear Applications of Accelerator Technology.USA：American Nuclear Society，1999：404-409.

［22］GABRIEL T A.CALOR：A Monte Carlo program package for the design and analysis of calorimeter systems，ORNL／TM-5619［R］.USA：ORNL，1977.

［23］CLOTH P，F(xiàn)ILGES D，NEEF R D，et al.HERMES，a Monte Carlo program system for beam-materials interaction studies［R］.Germany：KFA，1988.

［24］JOHNSON J O.Shielding design of the spallation neutron source，shielding aspects of accelerators，targets and irradiation facilities［C］∥SATIF 4 Workshop Proceedings.Knoxville：［s.n.］，1998：89-100.

［25］JOHNSON J O，SANTORO R T，LILLIE R A，et al.The SNS target station preliminary title I shielding analyses［J］.J Nucl Sci Technol，2000，37（Suppl.）：35-39.

［26］CHEN Y.Coupled Monte Carlo discrete ordinates computational scheme for three-dimensional shielding calculations of large and complex nuclear facilities［D］.Forschungszentrum Karlsruhe GmbH：Institute Für Reaktorsicherheit，2005.

［27］韓靜茹.三維蒙特卡羅-離散縱標(biāo)雙向耦合屏蔽計(jì)算方法研究［D］.北京：華北電力大學(xué)，2012.