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

畢遠杰，楊宏偉，郭金森，駱志平，李傳龍，陳 凌

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

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

隨著計算機技術的發展，目前在大型核設施的退役過程中，人們更希望能通過計算機模擬給出整個廠房內部的三維劑量場分布及其動態變化情況，為工作人員在放射性工作場所內進行放射性操作時的工作方案設計提供保障。根據計算機的模擬結果，工作人員可對不同的退役方案進行評估，選擇最優化的方案進行實施。因而，進行三維輻射場快速計算成為時代的需求。

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

1 輻射場計算方法現狀

1.1 傳統計算方案

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

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

1.2 輻射場耦合計算

蒙卡、離散縱標和點核方法各有優缺點（表1）。單一的計算方法均存在各自的不足，蒙卡方法對于大空間厚屏蔽問題無法在合理的時間內給出可靠的計算結果，確定論方法則不適合復雜源項和幾何區域。耦合計算則可將整個問題劃分為兩個區域進行求解，在復雜源項和幾何區域采用蒙卡方法進行模擬，在大空間厚屏蔽區域采用確定論方法求解深穿透問題。

1）蒙卡-離散縱標耦合方法

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

表1 3種計算方法的優缺點Table 1 Comparison of advantage and disadvantage of three methods

2）蒙卡-點核耦合方法

蒙卡與點核耦合進行全空間輻射場快速計算還未見詳細報道，蒙卡-點核耦合方法的特點為：（1）可給出全空間的三維輻射場，解決點核方法在復雜源項和幾何區域存在的問題；（2）相對于蒙卡-離散縱標耦合方法的趨于精確計算，蒙卡-點核耦合方法更加快速，可隨著設備的拆卸在盡量短的時間內給出新的劑量場分布。

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

2 耦合計算系統

圖1為耦合計算程序系統示意圖。在耦合計算中，將整個問題劃分為兩個區域：復雜源項和幾何區域及簡單幾何大屏蔽區域。復雜源項和幾何區域的結構通常較復雜，所以在這種區域采用蒙卡方法進行粒子輸運模擬。而在大屏蔽部分，幾何通常較簡單，同時屬于粒子深穿透問題，采用點核方法進行計算。結合蒙卡程序和點核計算的結果，得到三維輻射場快速計算軟件，并對軟件進行基準校驗和對比分析。

圖1 耦合計算程序系統示意圖Fig.1 Schematic of coupled calculation program

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

在計算過程中，首先對交界面進行空間上的網格劃分得到面元dA，接著對dA 出射的粒子進行方位角離散和能量分群，最終得到（dA，dE，dΩ）。設dA 的中心點為A′、外法向為nA，則位于P 點的探測器和dA 的外法線方向的夾角，因 此，dA 對P點探測器的貢獻為：

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

3 退役算例

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

圖2 計算算例幾何和參考點Fig.2 Geometry of model and reference points

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

計算結果采用通用蒙卡軟件FLUKA 進行校準，表2比較了FLUKA 和耦合程序的計算結果。對于選擇的參考點，兩個計算程序的相對偏差小于10%。在計算時間方面，將整個房間劃分為105個網格進行三維輻射場計算，在3.4GHz、4GB內存的個人計算機上進行測試，FLUKA 達到表2 中的計算誤差所需的計算時間為12h，耦合程序用20min得到相對誤差1%的等效面源，18min進行點核計算。

表2 參考點周圍劑量當量率比對Table 2 Comparison of ambient dose equivalent in reference points

4 結論

本文基于蒙卡與點核耦合的輻射場計算方法，編制了相應的耦合計算程序，用于研究堆退役過程中的輻射場快速計算。蒙卡-點核耦合方法結合了蒙卡和點核計算各自的優點，不僅適合于復雜源項和屏蔽體區域的計算，同時也適合于大屏蔽區域的輻射場快速計算。作為算例，本文采用編制的耦合程序對研究堆退役過程中放射性廢物儲存房間的輻射場進行了計算，計算結果驗證了程序的可靠性和適用性。

［1］ CHILTON B，SHULTIS J K，FAW R E.Principles of radiation shielding［M］.USA：Prentice-Hall，Inc.，1984.

［2］ SHULTIS J K，FAW 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，Feature，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，FERRARI 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，FILGES 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］韓靜茹.三維蒙特卡羅-離散縱標雙向耦合屏蔽計算方法研究［D］.北京：華北電力大學，2012.