全局減方差方法在乏燃料干式貯存容器屏蔽計算中的應用

王夢琪，鄭 征，梅其良，黎 輝，程湯培

(1.上海核工程研究設計院有限公司，上海 200233；2.北京應用物理與計算數(shù)學研究所，北京 100190)

乏燃料干式貯存容器屏蔽設計，由于其不規(guī)則的通風、力學及屏蔽結(jié)構(gòu)設計，以及乏燃料組件和結(jié)構(gòu)部件復雜的光中子源分布，是典型的不規(guī)則幾何和中光子耦合源項問題。蒙特卡羅(MC)方法由于其數(shù)據(jù)庫精度高、幾何處理能力強、適用于中光子耦合輸運計算等特點，是乏燃料干式貯存容器屏蔽計算的主要方法之一。乏燃料干式貯存容器要將放射性極高的乏燃料組件輻射水平降至工作人員能常規(guī)進入的維護水平，是典型的深穿透屏蔽計算問題。對于深穿透問題，MC方法存在計算效率低、不易收斂、難以獲得可信結(jié)果等問題。

為解決MC方法在計算過程中遇到的深穿透和難收斂問題，提出了基于分段續(xù)算思想的方法，如面源續(xù)算、離散縱標(SN)方法與MC方法耦合計算等[1-3]。此外Wagner等[4]提出了CADIS(Consistent Adjoint Driven Importance Sampling)方法，該方法基于SN方法計算得到的共軛注量率生成MC方法的源偏倚和權(quán)窗參數(shù)，可實現(xiàn)對MC計算的局部區(qū)域加速。在CADIS方法的基礎上，提出了FW-CADIS(Forward-CADIS)方法[5]，F(xiàn)W-CADIS方法首先進行正向SN計算，基于正向SN計算結(jié)果對共軛SN計算源項進行權(quán)重(或偏倚)，根據(jù)共軛SN計算得到的注量率生成MC方法的源偏倚和權(quán)窗參數(shù)，然后進行正向MC計算。相對于CADIS方法，F(xiàn)W-CADIS方法充分利用了正向SN計算得到的全局信息，能實現(xiàn)全局/半全局分布量的加速，目前在堆腔漏束、乏燃料干式貯存、中子測井等深穿透屏蔽計算中得到了應用[6-8]，計算效率可提高2～4個數(shù)量級[9]。

本文基于FW-CADIS方法開發(fā)適用于中光子耦合源的全局減方差方法程序進行乏燃料干式貯存容器屏蔽計算。

1 理論模型

1.1 全局減方差方法

考慮到共軛注量率表示粒子對探測器響應的貢獻，可利用共軛注量率對正向源進行偏倚，得到偏倚源分布，從而在重要的相空間內(nèi)抽樣產(chǎn)生更多的對探測器響應貢獻較大的源粒子。權(quán)窗技巧根據(jù)空間或能量相關(guān)的重要性對粒子進行分裂和輪盤賭，對輸運過程進行偏倚，增加到達重要相空間的粒子數(shù)。權(quán)窗下限和共軛注量率呈反比關(guān)系。

根據(jù)文獻[10]，為了進行全局分布量的加速，可推導出用于探測器響應計算的共軛源為：

(1)

式中：S+(r,E)為共軛源；r為空間位置；E和E′為能量；φ為正向注量率；Σd為探測器響應函數(shù)。

從式(1)可看出，正向注量率與共軛源強呈反比關(guān)系。偏倚后的共軛源使得MC粒子在探測器內(nèi)均勻分布，能在探測器內(nèi)得到近似均勻的計數(shù)誤差。共軛源S+(r,E)不局限于某類型粒子，可以是中子、光子或中光子耦合源。但當光子源和中子源源強相差很大時，建議采用如下中子和初級光子分別計算的方式：1) 對于中子源，在正向SN計算時，計算得到共軛中子源和次級光子源，然后進行共軛SN計算，得到共軛中子注量率和次級光子注量率；2) 對于初級光子源，通過正向計算對正向光子源進行偏倚，然后共軛計算得到共軛光子注量率。

根據(jù)共軛注量率φ+(q,g,p)可計算出能量、空間和粒子類型源偏倚因子：

(2)

(3)

(4)

式中：q為網(wǎng)格編號；g為能群編號；p為粒子類型編號；T、B、R分別為能量、空間、粒子類型源偏倚因子；Q為源區(qū)網(wǎng)格總數(shù)；G為能群數(shù)。

為了使得結(jié)果無偏，對源粒子權(quán)重進行修正得到偏倚后的權(quán)重：

(5)

通過上述公式可得到中子、次級光子、初級光子各自對應的源偏倚和權(quán)窗參數(shù)，最終利用源偏倚和權(quán)窗參數(shù)進行MC計算，得到中子劑量率、次級光子劑量率和初級光子劑量率。

1.2 計算流程

全局減方差方法計算流程如圖1所示。全局減方差方法需進行1次正向SN計算、1次共軛SN計算和1次正向MC計算。主要計算步驟包括：1) 建立三維正向SN計算模型，計算得到三維多群正向注量率分布，利用正向注量率計算共軛源強；2) 建立三維共軛SN計算模型，計算得到三維多群共軛注量率分布，共軛計算源項空間分布在探測器位置，能譜為響應函數(shù)能譜，且能群與正向相反；3) 讀取三維多群共軛注量率，根據(jù)共軛注量率計算源偏倚因子，計算源粒子能量、空間和粒子類型累積分布函數(shù)，根據(jù)共軛注量率計算MC計算所需的權(quán)窗下限；4) 在MC輸入文件中增加權(quán)窗參數(shù)卡，進行正向輸運計算。

圖1 全局減方差方法計算流程Fig.1 Global variance reduction method calculation flowchart

2 乏燃料干式貯存容器計算過程

2.1 問題介紹

乏燃料干式貯存容器內(nèi)層為密封容器，裝有21盒乏燃料組件，每盒乏燃料組件內(nèi)有204個燃料元件，燃料元件布置方式為15×15。密封容器外圍結(jié)構(gòu)包括貯存容器屏蔽層、底板、頂蓋、通風口等，其結(jié)構(gòu)如圖2所示，其中主屏蔽體為1 m厚的混凝土，內(nèi)外表面為2.5 cm厚的鋼敷面。乏燃料組件源項是光中子耦合源，中子源主要來源于錒系核素的自發(fā)裂變以及燃料區(qū)氧化物的(α，n)反應。光子源則主要是由裂變產(chǎn)物、錒系核素及其子核衰變產(chǎn)生的，以及由結(jié)構(gòu)材料中由于中子活化產(chǎn)生的60Co衰變導致的。本文所計算的乏燃料組件源項如圖3所示，單盒乏燃料組件光子源強為3.39×1015s-1，中子源強為7.27×108s-1，光子源強比中子源強大近7個數(shù)量級。雖然中子源強較低，但由于其穿透力強，在接近容器外圍區(qū)域時，中子以及由于中子導致的次級光子貢獻并不可忽略。

圖2 乏燃料干式貯存容器軸向及徑向剖面示意圖Fig.2 Axial and radial schematic diagram of spent fuel dry storage cask

圖3 乏燃料組件源項Fig.3 Source term of spent fuel assembly

2.2 計算過程

全局減方差方法可適用于中光子耦合源問題，當光子源強和中子源強量級相差不懸殊時，通常采用中光子耦合源一步計算。本文光子源強比中子源強大近7個數(shù)量級，為了使得由于中子導致的次級光子和乏燃料組件本身的初級光子均能收斂，采用中子源和光子源分別計算的方案。

1) 中子源計算

正向SN計算能譜為乏燃料中子源譜，共軛SN計算采用中子通量-劑量率轉(zhuǎn)換因子作為共軛源的能譜，正向SN網(wǎng)格劃分與共軛SN計算相同。SN計算網(wǎng)格數(shù)為R×θ×Z=70×57×89，散射截面展開階數(shù)為P3，求積組階數(shù)為S8，數(shù)據(jù)庫為47群中子。因為SN計算只是為了獲得粒子輸運的趨勢，因此網(wǎng)格可不用劃分太細，以節(jié)省計算時間。

計算目標為方位角在0°～45°的1/8模型全區(qū)域。圖4示出中子源偏倚后前10萬抽樣粒子空間分布及各能量段的概率密度函數(shù)(PDF)。由圖4可見：粒子在空間上明顯向共軛源區(qū)域附近偏倚；源中子能量向高能偏倚，但低能中子仍占有不可忽略的份額，因為中子(包括低能中子)輸運過程中會與材料發(fā)生俘獲反應產(chǎn)生次級光子，故從圖4b中各能量段中子對計數(shù)區(qū)域的貢獻份額可見，低能中子的貢獻仍不可忽略。

圖4 偏倚中子源空間(a)及能量概率密度函數(shù)(b)分布Fig.4 Space (a) and PDF (b) distributions for biased neutron source

權(quán)窗下限是與能量和空間相關(guān)的函數(shù)，在輸運計算過程中每群中子和光子均有1套全區(qū)域的權(quán)窗下限，圖5示出3.67～4.97 MeV的中子及1.0～1.5 MeV的次級光子的權(quán)窗下限在空間的分布。由圖5可見，權(quán)窗下限與共軛注量率呈反比關(guān)系，收斂難度越大的計數(shù)區(qū)域權(quán)窗下限越低，在整個模型區(qū)域內(nèi)權(quán)窗下限相差近10個數(shù)量級。

2) 光子源計算

正向SN計算能譜為乏燃料初級光子源譜，共軛SN計算采用光子通量-劑量率轉(zhuǎn)換因子作為共軛源的能譜，正向及共軛SN計算網(wǎng)格數(shù)為R×θ×Z=70×57×89，散射截面展開階數(shù)為P3，求積組階數(shù)為S8，數(shù)據(jù)庫為20群光子。

由于光子在混凝土和鋼等屏蔽材料中相對中子衰減更快、自由程更短，導致由于初級光子在各區(qū)域產(chǎn)生的劑量場更加難收斂。通過分析發(fā)現(xiàn)，目前即便采用全局減方差方法想使得整個徑向和軸向范圍內(nèi)的初級光子劑量場分布收斂仍十分困難。本文選取乏燃料貯存容器關(guān)鍵區(qū)域(貯存容器外表面的活性區(qū)中平面區(qū)域、底部通風口區(qū)域、頂部通風口區(qū)域)進行分析，故將共軛源放置在上述局部區(qū)域，具體位置如圖6a所示。

圖6示出初級光子源偏倚后前50萬偏粒子抽樣的空間分布及各能量段的概率密度函數(shù)。由圖6可見：粒子在空間上明顯向共軛源區(qū)域附近偏倚；初級光子很難穿透到容器外，能量小于1.0 MeV的光子對容器外表面的貢獻小于1.0×10-5，因此能量上得到很大提高，無偏初級光子源的能量主要集中在0.6～0.7 MeV，偏倚后初級光子能量主要在1.5 MeV以上。

圖5 中子(a)和次級光子(b)的權(quán)窗下限分布Fig.5 Weight-window lower bound for neutron (a) and secondary photon (b)

圖6 偏倚光子源空間(a)及能量概率密度函數(shù)(b)分布Fig.6 Space (a) and PDF (b) distributions for biased photon source

能量為1.0～1.5 MeV的初級光子，其權(quán)窗下限在空間分布的剖面如圖7所示。

圖7 初級光子的權(quán)窗下限分布Fig.7 Weight-window lower bound for primary photon

3 數(shù)值結(jié)果及分析

3.1 基于乏燃料中子源的結(jié)果分析

中子及次級光子劑量率計算時各程序所使用的進程數(shù)和并行計算時間列于表1。計算平臺硬件參數(shù)為：64位SUSE Linux Enterprise Server 11系統(tǒng)，644核Intel Xeon CPU E7-8830處理器，主頻為2.4 GHz。

圖8示出劑量率在0°～45°角全區(qū)域分布的剖面圖。圖9示出在相同計算條件下采用與不采用全局減方差方法區(qū)域內(nèi)所有計數(shù)網(wǎng)格收斂計數(shù)誤差分布。由圖9可知：采用全局減方差方法后有92.3%網(wǎng)格計數(shù)相對誤差小于10.0%，有98.8%的網(wǎng)格計數(shù)相對誤差在20.0%以下；不采用全局減方差方法只有少量(<49.7%)靠近源區(qū)的網(wǎng)格能收斂到相對誤差在10.0%以下，有58.2%的網(wǎng)格計數(shù)相對誤差在20.0%以下。參照文獻[11]中的方法，對所有計算區(qū)域進行品質(zhì)因子(FOM)統(tǒng)計，采用MC直接計算后中子劑量率的FOM為1.93×10-1，次級光子劑量率的FOM為2.90×10-4；采用全局減方差方法后中子劑量率的FOM為2.16×10-3，次級光子劑量率的FOM為6.14×10-2。可見全局減方差方法能使得中子收斂效率提高約89倍，次級光子收斂效率提高212倍。

表1 中子及次級光子劑量率計算所使用的進程數(shù)和計算時間Table 1 Number of process and time for calculation of neutron and secondary photon dose rate

圖8 中子劑量率(a)和次級光子劑量率(b)分布Fig.8 Dose rate distribution for neutron (a) and secondary photon (b)

圖9 網(wǎng)格計數(shù)收斂比例分布Fig.9 Convergent proportional distribution of mesh tally

表2列出幾處關(guān)鍵位置直接計算、全局減方差方法和三維SN程序JSNT的計算結(jié)果比較。由表2可知：直接計算源區(qū)收斂效果很好，但容器外計數(shù)無法收斂；全局減方差方法乏燃料容器外中子和次級光子結(jié)果均已收斂。同已收斂的直接計算結(jié)果(源區(qū))相比，中子劑量率相對偏差為-0.67%，次級光子劑量率相對偏差為0.40%，結(jié)果吻合良好。與SN的計算結(jié)果比較可知：中子劑量率結(jié)果吻合良好，其中最大相對偏差為-20.1%；次級光子劑量率在源區(qū)相對偏差為-2.0%，吻合良好，容器外表面及通風口處SN結(jié)果相對偏小，相對偏差大于20%，其偏差主要來自于SN計算時所使用的Bugle-96數(shù)據(jù)庫與MC連續(xù)能量截面庫的差別[3,7]。由于SN計算所使用的數(shù)據(jù)庫Bugle-96熱中子群(E<0.625 eV)只有兩群半，且熱中子截面不保守，導致容器外熱中子注量率水平和由于(n，γ)反應產(chǎn)生的次級光子偏小。

表2 乏燃料干式貯存容器中子及次級光子的劑量率Table 2 Dose rate of neutron and secondary photon for spent fuel dry storage cask

注：括號內(nèi)為MC計算結(jié)果的統(tǒng)計偏差

3.2 基于乏燃料光子源的結(jié)果分析

初級光子劑量率計算時各程序所使用的進程數(shù)和并行計算時間列于表3。

表3 初級光子劑量率計算所使用的進程數(shù)和計算時間Table 3 Number of process and time for calculation of primary photon dose rate

圖10示出初級光子劑量率分布。由圖10a可見，MC直接計算大部分區(qū)域沒有計數(shù)，采用全局減方差方法后，粒子向關(guān)注區(qū)域偏倚，探測器外關(guān)注區(qū)域基本均有計數(shù)。針對計算區(qū)域的FOM進行統(tǒng)計分析：對于整個模型區(qū)域，MC直接計算的FOM為2.20×10-5，采用全局減方差方法后的FOM為9.14×10-5，可見全局減方差方法能使光子源全局收斂效率提高約4倍，收斂效果提升并不明顯，這是因為初級光子計算時共軛源放在重點關(guān)注的局部區(qū)域(貯存容器外表面的活性區(qū)中平面區(qū)域、底部通風口區(qū)域、頂部通風口區(qū)域)，而其他未關(guān)注區(qū)域的收斂效果并未得到改善。表4列出乏燃料干式貯存容器的初級光子劑量率。對表4所列幾處關(guān)鍵位置的FOM進行分析，其中源區(qū)MC直接計算的FOM為2.69，采用全局減方差方法后的FOM為1.23，可見對于源區(qū)計數(shù)，采用全局減方差方法后收斂效率反而下降。這是因為采用全局減方差技巧后，源粒子在空間上明顯向共軛源區(qū)域附近偏倚，源內(nèi)部大部分區(qū)域很難抽到粒子；對于貯存容器外表面的活性區(qū)中平面區(qū)域、底部通風口區(qū)域、頂部通風口區(qū)域，MC直接計算均無計數(shù)，F(xiàn)OM為1.42×10-5，采用全局減方差方法后，上述3個區(qū)域的FOM分別為5.04×10-4、1.88×10-3、1.75×10-4，收斂效率是MC直接計算的12～132倍，收斂效率提高1～2個數(shù)量級。

圖10b為乏燃料核材料區(qū)中平面高度光子劑量率沿徑向分布及MC全局減方差結(jié)果與SN結(jié)果的相對偏差。由圖10b可知：兩種方法計算的光子劑量率吻合良好，平均相對偏差約為14.5%，在乏燃料容器外表面最大相對偏差為-18.7%。在屏蔽層中兩種方法計算的光子劑量率的相對偏差較大，在20%～40%之間，主要原因為：1) MC計數(shù)網(wǎng)格和SN計數(shù)網(wǎng)格不完全一致，本文采用中心點插值的方法計算同一位置處的劑量率水平，由于在屏蔽層中通量梯度很大，受周圍網(wǎng)格大小的影響更明顯；2) MC全局減方差計算在屏蔽層中未完全收斂，計數(shù)相對誤差在5%～15%之間，結(jié)果的統(tǒng)計誤差較大。

a——軸向剖面圖；b——核材料區(qū)中平面高度沿徑向分布圖10 初級光子劑量率分布Fig.10 Dose rate distribution for primary photon

表4 乏燃料干式貯存容器的初級光子劑量率Table 4 Primary photon dose rate for spent fuel dry storage cask

由表4可知：MC直接計算采用140核并行，共計算593 min，容器外無法獲得計數(shù)，這是由于初級光子能量主要集中在1.0 MeV以下，而能量小于1.0 MeV的光子幾乎無法穿出到達乏燃料容器外部，對容器外表面的貢獻小于1.0×10-5，源粒子抽樣時1.0 MeV以下的初級光子占源粒子總數(shù)的98%以上，若不進行源偏倚以及權(quán)窗，即便計算更長的時間容器外也無法獲得計數(shù)；MC全局減方差方法與MC直接計算源區(qū)已收斂結(jié)果比較可知，二者相對偏差為5.8%，吻合良好；MC全局減方差方法同SN的計算結(jié)果比較可知，MC全局減方差方法初級光子結(jié)果偏小，其中源區(qū)和側(cè)面筒體外表面均已收斂，相對偏差分別為0.9%和15.0%，吻合良好；采用全局減方差方法后，大部分區(qū)域獲得計數(shù)，但底部通風口和頂部通風口仍很難收斂，本文在采用140核并行、共計算593 min的情況下在底部和頂部通風口計數(shù)率相對誤差大于10%，同SN的計算結(jié)果相比相對偏差為25.0%～39.6%，這是因為底部通風口MC全局減方差方法計算未完全收斂，導致結(jié)果仍有較大不確定性。

4 結(jié)論

本文開展了基于SN方法的全局減方差方法在乏燃料干式貯存容器屏蔽計算中的應用研究，針對乏燃料干式貯存容器，分別建立了中子源及光子源MC直接計算模型、SN計算模型及全局減方差方法計算模型，通過和MC直接計算及SN計算結(jié)果的比較，可得到以下結(jié)論。

1) 針對其中子源部分進行計算時，在計算時間相當?shù)那闆r下，采用全局減方差方法后乏燃料容器外中子和次級光子結(jié)果均已收斂，整個計算區(qū)域超過92.3%網(wǎng)格計數(shù)相對誤差能小于10.0%，而MC直接計算只有少量(<49.7%)靠近源區(qū)的網(wǎng)格能收斂到相對誤差在10.0%以下，通過對所有計算區(qū)域進行FOM統(tǒng)計分析，可知全局減方差方法能使中子及次級光子收斂效率提高約2個數(shù)量級。

2) 當光子源和中子源源強相差懸殊時，基于中子源而生成的次級光子權(quán)窗參數(shù)起到了很好的效果，次級光子源伴隨中子源的收斂而收斂。

3) 針對初級光子的計算，無偏初級光子源的能量主要集中在0.6～0.7 MeV，由于其衰減系數(shù)大、自由程短，很難達到貯存容器外表面，導致MC直接計算無法獲得計數(shù)。全局減方差方法通過源偏倚使得抽樣源粒子平均能量達到1.5 MeV以上，結(jié)合相空間權(quán)窗設置，使得大部分區(qū)域初級光子均有計數(shù)，已收斂區(qū)域和MC直接計算及SN計算結(jié)果的相對偏差均在20%以內(nèi)，通過FOM統(tǒng)計分析可知，全局減方差方法可使乏燃料容器外區(qū)域初級光子收斂效率提高1～2個數(shù)量級。同時研究發(fā)現(xiàn)，初級光子的深穿透收斂較中子更加困難，實現(xiàn)大片區(qū)域或離散區(qū)域的計數(shù)同時收斂計算代價很大，后續(xù)將進一步開展初級光子共軛源位置選取及源強設置等研究工作。