仲港其 徐 坤 陸玉東 毛世峰 葉民友

1（中國科學(xué)技術(shù)大學(xué)核科學(xué)技術(shù)學(xué)院 合肥 230026）

2（中國科學(xué)院等離子體物理研究所 合肥 230031）

聚變中子學(xué)分析對聚變堆的設(shè)計與建造、安全運行與維護具有重要意義。聚變中子學(xué)分析即應(yīng)用中子學(xué)理論、輸運計算程序以及核數(shù)據(jù)庫研究中子與結(jié)構(gòu)材料的相互作用，聚變中子由等離子體區(qū)域向外部擴散的過程中會與包層、偏濾器、真空室、超導(dǎo)磁體等部件的材料發(fā)生作用，生成中子、光子和電子等粒子，進而引發(fā)嬗變、活化以及輻照損傷等問題。目前，中子學(xué)分析主要有兩種方法：確定論方法和概率論方法（蒙特卡羅方法）［1］。由于托卡馬克模型的復(fù)雜性，聚變中子學(xué)分析采用的主要模擬方法是蒙特卡羅方法與程序，其中應(yīng)用最廣泛的是美國洛斯阿拉莫斯國家實驗室（Los Alamos National Laboratory，LANL）開發(fā)的三維蒙特卡羅程序MCNP（Monte Carlo N Particle Transport Code）［2］，國內(nèi)外相關(guān)研究機構(gòu)也開展了其他三維蒙特卡羅程序的應(yīng)用研究，如Serpent［3-4］、Tripoli［5］、Geant4［6］和OpenMC等。OpenMC是麻省理工學(xué)院于2011年開始開發(fā)的面向現(xiàn)代化計算機集群并行高效計算的開源蒙特卡羅粒子輸運軟件，采用高性能的并行算法，提高了并行擴展性，減少了進程通信時間，計算效率更為高效［7］。英 國 原 子 能 管 理 局（UK Atomic Energy Authority，UKAEA）、德國卡爾斯魯厄理工學(xué)院（Karlsruhe Institute of Technology，KIT）以及威斯康星大學(xué)麥迪遜分校等研究機構(gòu)不斷推進OpenMC在聚變中子學(xué)分析中的應(yīng)用研究，包括三維CAD建模功能、減方差方法的應(yīng)用以及基準測試等［8-9］。KIT的研究人員在OpenMC中添加了McDeLicious計算代碼，并且基于IFMIF-DONES裝置模型對比了OpenMC-McDeLicious與MCNP-McDeLicious計 算結(jié)果，獲得了良好的一致性結(jié)果［10］。由于聚變堆工程模型極其復(fù)雜，使用傳統(tǒng)的構(gòu)造實體幾何（Constructive Solid Geometry，CSG）方法建模十分繁瑣和耗時，目前主要有兩種方式解決該問題：CAD模型到CSG模型的轉(zhuǎn)換，典型的模型轉(zhuǎn)換程序包括CMGC［11］和McCAD［12］等；直接基于CAD模型進行粒子輸運模擬。DAGMC（Direct Accelerated Geometry Monte Carlo）是威斯康星大學(xué)麥迪遜分校開發(fā)的開源工具包，通過在蒙特卡羅軟件中插入核心幾何庫，實現(xiàn)了基于CAD模型的直接粒子輸運模擬［13］，避免了通過構(gòu)造實體幾何方法建模。DAGMC已經(jīng)被集成到多個蒙特卡羅程序中，如DAG-MCNP、DAG-Geant4和DAG-OpenMC等。DAG-OpenMC繼承了DAGMC與OpenMC的優(yōu)點：基于CAD建模并直接在CAD模型上執(zhí)行粒子輸運模擬；高性能并行計算。這些優(yōu)點在復(fù)雜聚變工程中子學(xué)分析中具有極大優(yōu)勢，被認為是解決聚變中子學(xué)分析問題的潛在方法。威斯康辛大學(xué)的相關(guān)人員通過耦合DAG-OpenMC與PyNE R2S，對FNGITER停堆劑量率基準題展開測試，通過對比MCNP、DAG-OpenMC以及實驗結(jié)果證明了DAGOpenMC和PyNE R2S耦合可以有效計算聚變系統(tǒng)的停堆劑量率［14］。由于聚變裝置工程模型非常復(fù)雜，建模十分煩瑣和耗時，并且等離子體區(qū)域中子源的空間分布極其復(fù)雜，難以準確定義。針對上述問題，本文開展了OpenMC在聚變堆中子學(xué)工程分析中的應(yīng)用研究，探究了OpenMC在復(fù)雜工程結(jié)構(gòu)條件下的聚變中子學(xué)分析中的準確性與適用性。

中國聚變工程試驗堆（China Fusion Engineering Test Reactor，CFETR）是中國磁約束聚變能計劃的下一步試驗裝置，旨在為未來建造的示范堆提供參考［15］。本文首先基于CFETR一維柱殼中子學(xué)模型驗證了OpenMC與MCNP程序關(guān)于聚變中子學(xué)典型物理量計算結(jié)果的一致性，表明了OpenMC可以準確計算聚變中子學(xué)物理量。進一步探究了DAG-OpenMC基于CAD建立復(fù)雜聚變堆工程模型的功能；基于平衡等離子體密度和溫度分布計算得到等離子體中子輻射系數(shù)矩陣，根據(jù)該系數(shù)矩陣編寫C++源文件定義CFETR三維中子學(xué)分析模型的D-T中子源分布，解決了OpenMC自帶的源定義功能無法描述復(fù)雜的聚變中子源分布問題。將DAG-OpenMC應(yīng)用于CFETR三維中子學(xué)模型，計算得到了中子壁負載分布、包層氚增殖率以及核熱沉積等。

1 方法

1.1 基于CAD建模

DAGMC采用實體建模軟件Trelis作為前處理模塊，實現(xiàn)對CAD模型的預(yù)處理并轉(zhuǎn)化為面網(wǎng)格模型?；赥relis預(yù)處理模型的過程如圖1所示，首先利用建模軟件（如CATIA）創(chuàng)建托卡馬克的復(fù)雜CAD模型，然后利用SpaceClaim對CAD模型進行簡化，簡化后將模型導(dǎo)入Trelis進行預(yù)處理。DAGOpenMC在粒子輸運過程中執(zhí)行嚴格的模型檢查，即模型不允許出現(xiàn)幾何干涉，因此要在Trelis中進行檢測并修復(fù)模型。模型檢測通過后，對模型所有表面進行“壓印”（Imprint），確保面的準確定義。表面“壓印”成功后，合并模型并再次檢測模型，此時檢測幾何干涉和表面重疊。重復(fù)上述步驟直至Trelis對模型預(yù)處理成功，即模型不存在任何錯誤。

圖1 DAGMC建模過程[16]Fig.1 Modeling process of DAGMC[16]

模型預(yù)處理結(jié)束后，利用Svalinn提供的插件［17］在Trelis程序中為模型的幾何表面賦予邊界條件以及為幾何體賦予材料屬性，設(shè)定粒子輸運邊界限制粒子運動區(qū)域。在Trelis中處理好模型后，導(dǎo)出為dagmc.h5m文件作為OpenMC程序的幾何輸入。

1.2 等離子體源

托卡馬克裝置中復(fù)雜的等離子體物理，決定了中子源的復(fù)雜空間分布。MCNP采用Fortran source subroutine或者復(fù)雜sdef定義的格式，實現(xiàn)了復(fù)雜中子源的準確定義。而OpenMC提供的源定義功能較弱，在Fixed source運行模式下，僅支持定義點源、環(huán)形線源、球形源和立方體源。為了改善這一情況以及更精確的定義等離子體源，本文基于OpenMC的自定義源擴展接口，編寫C++源文件定義復(fù)雜的聚變中子源。源文件中包含一個自定義的源類和在該類中生成源中子的函數(shù)，在該函數(shù)中定義源中子的空間分布、權(quán)重、方向和能量等信息，隨后編譯成動態(tài)鏈接庫，在OpenMC程序運行時抽樣源中子。

2 基于CFETR一維柱殼模型的聚變中子學(xué)參數(shù)計算與對比

在中子學(xué)分析中，常采用一維模型為裝置的尺寸設(shè)計提供指導(dǎo)和優(yōu)化。本節(jié)基于CFETR一維柱殼模型，選用MCNP5作為對比程序驗證OpenMC在聚變中子學(xué)分析中的計算準確性，其中，MCNP版本為5.1.60，OpenMC版本為12.1。對比計算中，采用的中子截面庫均為FENDL3.1d［18］，該截面庫從ACE格式轉(zhuǎn)換為HDF5格式供OpenMC使用。OpenMC與MCNP模擬的粒子數(shù)量均為1×108個。CFETR功率為200 MW，其中聚變中子功率為160 MW，計算得出單位時間產(chǎn)生的中子數(shù)為7.11×1019。對比參數(shù)包括各部件中子通量與光子通量，包層中子能譜分布與光子能譜分布，核熱沉積以及包層氚增殖率等。

2.1 CFETR一維柱殼模型

CFETR一維柱殼模型是一個多層同心圓柱體，徑向上以多層圓柱殼結(jié)構(gòu)表示CFETR中平面處徑向各部件，以等離子體區(qū)域為界分為內(nèi)側(cè)和外側(cè)兩部分，上下平面設(shè)為反射面，中子學(xué)模型如圖2所示。高場側(cè)從等離子體區(qū)域至中心螺線管之間的結(jié)構(gòu)分別為：第一壁鎢鎧甲、第一壁、增殖包層、氦氣分配管組、屏蔽包層、真空室、冷屏、環(huán)向場線圈和中心螺線管；低場側(cè)從等離子體區(qū)域至中心螺線管之間的結(jié)構(gòu)分別為：第一壁鎢鎧甲、第一壁、增殖包層、氦氣分配管組、屏蔽包層、真空室、冷屏、環(huán)向場線圈、杜瓦和生物屏蔽墻。各部件使用均質(zhì)材料填充。各向同性的中子源均勻分布在等離子體區(qū)域內(nèi)，能量分布遵從Gauss聚變譜分布，平均能量為14.079 1 MeV。

圖2 CFETR一維中子學(xué)模型Fig.2 One dimensional neutronics model of CFETR

2.2 結(jié)果對比與分析

基于CFETR一維柱殼模型計算得到了在不同位置處的中子通量和光子通量（圖3）。OpenMC與MCNP計算的第一壁、包層、屏蔽層以及真空室等靠近等離子體區(qū)域的中子通量與光子通量偏差不超過0.1%，具備良好的一致性。表明即使OpenMC處理光子輸運時模擬的物理過程相較于MCNP存在部分不同，OpenMC的光子截面庫相比于MCNP額外包含Geant4 EMLOW數(shù)據(jù)庫中的Seltzer和Berger的軔致輻射截面、NIST ESTAR數(shù)據(jù)庫的平均激發(fā)能量以及Biggs等計算的康普頓截面數(shù)據(jù)［19］，但對結(jié)果并無太大影響。冷屏之外離等離子體較遠的區(qū)域，中子通量與光子通量偏差較大，可達5%左右，但偏差基本不超過兩倍OpenMC的統(tǒng)計誤差（2σ），推測是通量統(tǒng)計誤差造成的。

圖3 中子通量(a)和OpenMC/MCNP中子通量比率(b)，光子通量(c)和OpenMC/MCNP光子通量比率(d)Fig.3 Neutron flux(a)and OpenMC/MCNP neutron flux ratio(b),photon flux(c)and OpenMC/MCNP photon flux ratio(d)

聚變堆的核熱沉積是等離子體產(chǎn)生的中子能量沉積和光子能量沉積的總和。中子核熱沉積和光子核熱沉積如圖4所示，在不同位置處的OpenMC/MCNP核熱沉積比值呈現(xiàn)了與通量比值相似的情況，第一壁、包層、屏蔽層以及真空室等靠近等離子體的區(qū)域核熱沉積偏差基本不超過1%，冷屏之外離等離子體較遠的區(qū)域核熱沉積的偏差達5%～10%，但基本不超過兩倍OpenMC的統(tǒng)計誤差（2σ）。MCNP計算的總核熱為217.214 MW，其中中子總核熱為160.839 MW，光子總核熱為56.375 MW；OpenMC計算的總核熱為217.208 MW，其中中子總核熱為160.717 MW，光子總核熱為56.490 MW，OpenMC計算的中子總核熱相較于MCNP低0.076%，光子總核熱相較于MCNP則高0.204%，最終的總核熱相較于MCNP低0.002%?？傮w來說，OpenMC與MCNP的計算結(jié)果一致。

中子能譜和光子能譜分析選取低場側(cè)包層，能群劃分選擇Vitamin-j結(jié)構(gòu)，該結(jié)構(gòu)中子能譜的統(tǒng)計范圍為10-7～20 MeV，共175個能量箱，光子能譜的統(tǒng)計范圍為10-2～20 MeV，共125個能量箱。計算得到了OpenMC與MCNP的中子能譜圖及光子能譜圖（圖5）。當(dāng)中子能量大于10-5MeV時，OpenMC與MCNP計算出的中子通量偏差不超過1%。當(dāng)中子能量低于10-5MeV時，OpenMC與MCNP的中子通量偏差變大，可達5%左右，但仍不超過兩倍的OpenMC統(tǒng)計誤差（2σ）。當(dāng)光子能量位于10-2～10 MeV范圍內(nèi)時，OpenMC與MCNP計算出的光子通量偏差不超過1%。當(dāng)光子能量高于10 MeV時，光子通量偏差變大，達5%左右，但基本不超過兩倍的OpenMC統(tǒng)計誤差（2σ）。

氚增殖率（Tritium Breeding Ratio，TBR）是衡量聚變堆包層氚增殖能力的參數(shù)，表示單位時間內(nèi)包層中產(chǎn)生的氚原子數(shù)與等離子體中產(chǎn)生的聚變中子數(shù)的比值［20］。在CFETR中，主要產(chǎn)氚核素為6Li，同時，7Li和9Be也可產(chǎn)生極少量的氚。使用OpenMC與MCNP分別計算了6Li、7Li和9Be核素的氚增殖率，結(jié)果如表1所示。對于這三種核素，OpenMC計算的氚增殖率略低于MCNP。OpenMC的總氚增殖率為1.147 7，相較于MCNP的1.152 1低0.38%。總的來說，在計算氚增殖率時，OpenMC和MCNP的計算結(jié)果一致。

表1 OpenMC與MCNP包層6Li、7Li和9Be的TBR對比Table 1 Comparison of TBR of6Li,7Li and9Be between OpenMC and MCNP

總體來說，CFETR一維柱殼模型對比分析的結(jié)果表明OpenMC在中子學(xué)分析計算（包括通量、核熱、能譜及TBR等）方面與MCNP一致，因此，OpenMC可以應(yīng)用于聚變中子學(xué)分析。

3 面向CFETR工程模型的三維中子學(xué)分析應(yīng)用

CFETR工程結(jié)構(gòu)十分復(fù)雜，相較于一維柱殼模型，三維模型建模更加困難和繁雜，因此，減少建模工程的工作量尤為重要，而DAG-OpenMC的CAD幾何功能在處理聚變堆復(fù)雜模型時可以顯著提高建模效率，減少模型錯誤。聚變中子源對計算結(jié)果的可靠性具有直接影響，因此，在三維中子學(xué)分析中，使用符合實際情況的中子源來提高計算的準確性。

3.1 CFETR建模及中子源

CFETR第二階段的設(shè)計大半徑R=7.2 m，小半徑a=2.2 m。主機由16個扇段組成，每個扇段包括2個內(nèi)包層段和3個外包層段。典型的CFETR主機模型包括包層、真空室、偏濾器、冷屏、超導(dǎo)縱場線圈、中心螺線管和極向場線圈等部件。氦冷陶瓷增殖 包 層（Helium Cooled Ceramic Breeder blanket，HCCB）是CFETR的候選包層概念之一，該包層概念由核工業(yè)西南物理研究院提出，采用多模塊化設(shè)計方式，在極向上分布有11個包層模塊，內(nèi)包層有6塊（1#～6#），外包層有5塊（7#～11#）。HCCB包層選用6Li富集度為90%的Li4SiO4小球作為氚增殖劑，金屬Be小球作為中子倍增劑，Li4SiO4小球和Be小球以球床填充的方式交替分區(qū)布置在包層中［21］。HCCB包層單個扇段的布局和外側(cè)中平面處增殖包層模塊分別如圖6（a）和（b）所示?？紤]到CFETR模型的對稱性，進行三維中子學(xué)分析時，在環(huán)向方向上選取11.25°的扇區(qū)作為分析模型，并采用反射邊界條件模擬360°全堆模型。圖7展示了CFETR HCCB包層的中子學(xué)計算模型。

圖6 單個扇段包層模塊分布(a)和外側(cè)中平面包層模塊(b)Fig.6 Module distribution of single sector blanket(a)and outer mid-plane blanket module(b)

圖7 CFETR中子學(xué)模型的縱向剖面(a)和內(nèi)、外側(cè)結(jié)構(gòu)的水平剖面(b)Fig.7 Vertical cross section(a)and horizontal cross section of the CFETR neutronics model(b)

托卡馬克裝置發(fā)生D-T聚變反應(yīng)產(chǎn)生的中子平均能量為14.079 1 MeV，遵從Gauss聚變譜分布?？臻g分布大致呈D形，等離子體中心中子源密度最高，并向等離子體的邊緣逐漸減少?；谄胶獾入x子體密度和溫度分布計算得到CFETR等離子體中子輻射系數(shù)矩陣［22］，根據(jù)該系數(shù)矩陣以及參考MCNP網(wǎng)格離散體源的定義方式，將中子源離散成多個各向同性的體源，每個體源權(quán)重不同，在同一體源內(nèi)源中子均勻抽樣?；贠penMC提供的源擴展接口，使用C++實現(xiàn)了二維（R-Z）各向同性的40×40網(wǎng)格離散體源來描述CFETR模型中D-T中子源的空間分布，并將結(jié)果與MCNP-sdef定義的中子源作對比，如圖8所示。可以發(fā)現(xiàn)，兩者定義的中子源在空間分布上具有極好的一致性（在等離子體邊緣區(qū)域，源概率低于10-8，抽樣統(tǒng)計誤差較大，導(dǎo)致兩者空間分布在此處存在偏差）。

圖8 中子源空間分布(a)和OpenMC與MCNP的中子源空間分布概率比值(b)Fig.8 Spatial distribution of neutron sources(a)and Probability ratio of neutron source spatial distribution between OpenMC and MCNP(b)

3.2 中子壁負載分布

中子壁負載指的是托卡馬克裝置中D-T聚變反應(yīng)生成的14.08 MeV的中子直接作用到第一壁表面的中子能量通量密度。包層中子第一壁總負載約為144.4 MW，占160 MW聚變中子功率的90.25%。圖9為各個包層模塊的中子壁負載分布，內(nèi)包層中子壁負載峰值出現(xiàn)在3#包層，達到0.156 6 MW?m-2，外包層中子壁負載峰值出現(xiàn)在9#包層上，達到0.227 1 MW?m-2，中 子 壁 負 載 的 平 均 值 為0.169 4 MW?m-2。

圖9 中子壁負載分布Fig.9 Neutron wall load distribution

3.3 包層氚增殖率

計算得到了HCCB包層中各包層模塊在200 MW聚變功率下對氚增殖率的貢獻（圖10）。包層總TBR為1.141，其中外包層總TBR為0.798，內(nèi)包層總TBR為0.343。

圖10 各包層模塊對TBR的貢獻Fig.10 Contribution of each blanket module to TBR

3.4 核熱沉積

聚變堆的核熱沉積由中子核熱沉積和光子核熱沉積組成，全局能量倍增因子是全堆核熱總功率與聚變中子功率的比值［23］，是評估聚變堆核熱產(chǎn)出能力的重要參數(shù)之一。

表2列出了CFETR各部件的核熱沉積和占比情況，全堆產(chǎn)生的核熱共計211.73 MW，全局能量倍增因子為1.32。其中包層系統(tǒng)產(chǎn)生的核熱達185.64 MW，占總核熱的主要部分（占比87.68%），偏濾器產(chǎn)生的核熱為24.98 MW，占比11.80%。

表2 CFETR各部件核熱沉積Table 2 Nuclear heating of CFETR components

4 結(jié)語

開源代碼OpenMC為聚變中子學(xué)分析提供了一種選擇，該程序易于獲取和進行二次開發(fā)。本文基于CFETR一維柱殼模型對比分析了OpenMC與MCNP關(guān)于中子/光子通量、中子/光子核熱沉積、中子/光子能譜以及包層氚增殖率等聚變中子學(xué)參數(shù)的計算結(jié)果，結(jié)果顯示，OpenMC與MCNP在中子學(xué)分析計算方面具有良好的一致性，表明OpenMC可以應(yīng)用于聚變中子學(xué)分析。

在CFETR的三維中子學(xué)分析中，使用實體建模軟件Trelis對三維CFETR的CAD模型進行預(yù)處理，集成DAGMC實現(xiàn)OpenMC在該模型上直接進行粒子輸運模擬，實現(xiàn)了DAG-OpenMC中聚變堆的建模?；谟嬎愕玫降牡入x子體中子輻射系數(shù)矩陣，參考MCNP網(wǎng)格離散體源的定義方式，編寫C++源文件在OpenMC中實現(xiàn)了精確的D-T聚變中子源定義。因此，在聚變堆復(fù)雜工程結(jié)構(gòu)條件下，DAGOpenMC在建模及設(shè)置中子源方面均滿足要求。經(jīng)計算，CFETR中子壁負載峰值出現(xiàn)在中平面處，氦冷陶瓷增殖包層TBR為1.141。在聚變中子功率為160 MW時，產(chǎn)生的總核熱為211.73 MW，全局能量倍增因子為1.32，其中包層產(chǎn)生的總核熱為185.64 MW，占總核熱的87.68%。

作者貢獻聲明仲港其：實施研究，采集數(shù)據(jù)，分析/解釋數(shù)據(jù)，起草文章；徐坤：指導(dǎo)，統(tǒng)計分析，技術(shù)及軟件支持，對文章的知識型內(nèi)容作批評性審閱；陸玉東：采集數(shù)據(jù)；毛世峰：指導(dǎo)，支持性貢獻；葉民友：指導(dǎo)，對文章的知識型內(nèi)容作批評性審閱。