DAG-OpenMC在聚變中子學(xué)分析中的應(yīng)用研究

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

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

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

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

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

1 方法

1.1 基于CAD建模

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

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

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

1.2 等離子體源

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

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

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

2.1 CFETR一維柱殼模型

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

圖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%，冷屏之外離等離子體較遠(yuǎn)的區(qū)域核熱沉積的偏差達(dá)5%～10%，但基本不超過兩倍OpenMC的統(tǒng)計(jì)誤差（2σ）。MCNP計(jì)算的總核熱為217.214 MW，其中中子總核熱為160.839 MW，光子總核熱為56.375 MW；OpenMC計(jì)算的總核熱為217.208 MW，其中中子總核熱為160.717 MW，光子總核熱為56.490 MW，OpenMC計(jì)算的中子總核熱相較于MCNP低0.076%，光子總核熱相較于MCNP則高0.204%，最終的總核熱相較于MCNP低0.002%?？傮w來說，OpenMC與MCNP的計(jì)算結(jié)果一致。

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

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

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

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

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

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

3.1 CFETR建模及中子源

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

圖6 單個(gè)扇段包層模塊分布(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聚變譜分布。空間分布大致呈D形，等離子體中心中子源密度最高，并向等離子體的邊緣逐漸減少。基于平衡等離子體密度和溫度分布計(jì)算得到CFETR等離子體中子輻射系數(shù)矩陣［22］，根據(jù)該系數(shù)矩陣以及參考MCNP網(wǎng)格離散體源的定義方式，將中子源離散成多個(gè)各向同性的體源，每個(gè)體源權(quán)重不同，在同一體源內(nèi)源中子均勻抽樣?；贠penMC提供的源擴(kuò)展接口，使用C++實(shí)現(xiàn)了二維（R-Z）各向同性的40×40網(wǎng)格離散體源來描述CFETR模型中D-T中子源的空間分布，并將結(jié)果與MCNP-sdef定義的中子源作對比，如圖8所示?？梢园l(fā)現(xiàn)，兩者定義的中子源在空間分布上具有極好的一致性（在等離子體邊緣區(qū)域，源概率低于10-8，抽樣統(tǒng)計(jì)誤差較大，導(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 中子壁負(fù)載分布

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

圖9 中子壁負(fù)載分布Fig.9 Neutron wall load distribution

3.3 包層氚增殖率

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

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

3.4 核熱沉積

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

表2列出了CFETR各部件的核熱沉積和占比情況，全堆產(chǎn)生的核熱共計(jì)211.73 MW，全局能量倍增因子為1.32。其中包層系統(tǒng)產(chǎn)生的核熱達(dá)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é)分析提供了一種選擇，該程序易于獲取和進(jìn)行二次開發(fā)。本文基于CFETR一維柱殼模型對比分析了OpenMC與MCNP關(guān)于中子/光子通量、中子/光子核熱沉積、中子/光子能譜以及包層氚增殖率等聚變中子學(xué)參數(shù)的計(jì)算結(jié)果，結(jié)果顯示，OpenMC與MCNP在中子學(xué)分析計(jì)算方面具有良好的一致性，表明OpenMC可以應(yīng)用于聚變中子學(xué)分析。

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

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