沈華韻,潘流俊,衷 斌
(北京應用物理與計算數(shù)學研究所,北京 100094)
中子輸運外源問題是核能開發(fā)利用中廣泛存在且重要的研究課題,如核裝置的輻射屏蔽研究、加速器驅動次臨界系統(tǒng)(ADS)和聚變裂變混合堆的研發(fā)、在役核電站裝換料階段的臨界安全分析[1]。與屏蔽問題相比,后兩類問題有如下共同特點:系統(tǒng)中的裂變裝置處于次臨界狀態(tài),中子注量率直接受源強影響。為此,可將其統(tǒng)稱為外源驅動次臨界系統(tǒng)。
對外源驅動次臨界系統(tǒng)的中子物理學問題的研究,同樣可采用確定論(如SN)方法或蒙特卡羅(MC)方法[2]。但由于外中子源通常只分布在局部空間,使得SN方法的射線效應嚴重影響此類問題的計算精度。以大亞灣核電站為例,一次源的半徑和高僅0.053cm 和4.00cm,而堆芯的半徑和高分別為161cm 和183cm[1]。另外,SN方法難以精確描述源的復雜幾何結構,在一定程度上也將影響計算精度。MC 方法雖不存在上述問題,但計算效率較低,影響了該方法的應用。
雖然以ADS和混合堆為代表的外源驅動次臨界系統(tǒng)是當前的研究熱點,國內外均已開展或正在開展大量的研究工作,但這些研究往往關注于裝置的物理設計,而未涉及計算方法[3-5]。為此,本文針對外源驅動次臨界系統(tǒng)內中子輸運過程的特點,提出以中子首次裂變?yōu)轳詈宵c的MC/SN耦合計算方法。通過充分發(fā)揮MC 方法和SN方法各自的優(yōu)點,實現(xiàn)外源驅動問題高效高精度的模擬。在算法研究的基礎上實現(xiàn)MC/SN耦合輸運計算,并進行測試模型的計算與分析。
中子輸運方程可寫成3種不同的形式:通量型、碰撞密度型和發(fā)射密度型[6]??紤]到算法討論的便利性,本文將采用發(fā)射密度型積分方程。對于外源問題,該方程可寫為:
其中:P 為相空間r×E×Ω 上任意點;χ(P)為P 點的中子產生率,包括外源中子和中子核反應產生的次級中子;S(P)為外源中子產生率;散射核Ks(P′→P)為1個P′的中子經過輸運并發(fā)生非裂變反應產生狀態(tài)為P 的中子數(shù)量的期望值;裂變核Kf(P′→P)為1個P′的中子經過輸運并發(fā)生裂變反應產生狀態(tài)為P 的中子數(shù)量的期望值。式(1)的物理含義非常清晰,即中子發(fā)射密度等于源中子發(fā)射密度、非裂變核反應次級中子發(fā)射密度和裂變反應中子發(fā)射密度之和。
定義:
其中:
在式(2)兩邊對n 求和,同時利用式(3)和(4),得到:
對比式(1)可知,外源問題的解可表示為級數(shù)形式:
由式(2)和(4)可知,χ′(P)滿足如下方程:
χ″(P)滿足如下方程:
即有
其中,根據(jù)式(3),有:
可見,由式(7)與(8)共同給出的定義完全等價于式(1)。對比式(7)與式(1)可知,χ′(P)與χ(P)對應的外源均為S(P),僅在中子裂變反應的處理方式上存在差異,計算χ′(P)時需將裂變處理為吸收。因此,物理上χ′(P)包含兩方面對系統(tǒng)內中子分布的貢獻——外源中子的直接貢獻和外源中子通過非裂變反應產生次級中子帶來的貢獻。
明確χ′(P)的物理含義后,由式(9)表示的S1(P)的物理含義已非常清晰,具有首次裂變中子源的意義,對應的包含兩種來源的裂變反應產生的中子——外源中子直接引起的裂變和由外源中子發(fā)生非裂變反應產生的次級中子引起的裂變。與實際外源S(P)相比,S1(P)具有如下特點:運動方向上各向同性;分布于整個裂變區(qū)等。
最后,對比式(8)與(1)可知,χ″(P)與χ(P)的差異僅表現(xiàn)在對應的源分布上。前者為首次裂變中子源,后者為實際的中子源。顯然,χ″(P)包含了首次裂變中子通過各種核反應引起的對系統(tǒng)內中子分布的貢獻。因此,物理上也解釋了式(6)的具體含義。
由上述分析可知,求解式(1)可轉化為耦合求解式(7)和(8),這兩個方程的特點決定了各自適宜采用的方法。對于式(7),由于源S(P)通常具有空間局部分布、幾何結構較復雜等特點,適合采用MC 方法以保證計算精度。同時,由于已將裂變處理為吸收,不存在中子增殖過程,使得需要跟蹤的中子歷史較短,不會帶來過多計算量。對于式(8),由于源S1(P)分布于整個裂變區(qū),且各向同性,可采用SN方法,在保證計算精度的同時提高計算效率。基于此,形成如下以中子首次裂變?yōu)轳詈宵c的MC/SN耦合輸運算法。
1)MC 方法求解中子輸運方程(式(7))。從實際源分布S(P)中抽取中子,將裂變處理為吸收,并通過計數(shù)獲得該中子輸運過程對應的中子注量率φ′(P)和首次裂變源S1(P)。
2)SN方法求解輸運方程(式(8))。此時求解的是對應于S1(P)的外源問題,并獲得中子注量率φ″(P)。
3)疊加獲得輸運方程(式(1))對應的中子注量率φ(P)=φ′(P)+φ″(P)。
可見,該耦合輸運算法的思想是:將1個難以在兼顧精度和效率條件下實現(xiàn)求解的外源中子輸運問題,轉化為另外兩個耦合的外源問題,然后分別采用MC和SN方法實現(xiàn)耦合問題的求解。最后,需特別指出的是,該算法中的MC與SN耦合是易于實現(xiàn)且嚴格精確的,這由如下兩方面保證:1)僅存在MC 到SN的單向單次數(shù)據(jù)傳輸;2)首次裂變中子源S′(P)是各向同性的,可精確實現(xiàn)角度的離散。
為驗證耦合算法的有效性,構造如下裝置:二維xy 幾何,尺寸為20cm×20cm,x=0和y=0處為反射邊界,其他為自由面,材料為富集度10%的濃縮鈾。SN計算時采用S8,共80個離散方向。另外,假定源中子的能量為14.1 MeV,尺寸為1cm×1cm,考慮不同空間分布的情況。在系列算例中,均以等值線的形式給出14.1 MeV 和2 MeV 附近的中子注量率,并對標注作如下約定:MCNP表示僅用MC方法計算的結果,作為基準值;MCS8表示采用耦合算法計算的結果;S8表示僅用SN方法計算的結果。
模型1為僅在中心分布著單一中子源的情況,圖1示出了3種計算方式的結果。
從圖1a可見,SN方法的射線效應非常明顯。在計算精度上,無論是源中子所在能量附近還是裂變中子最可幾能量(2 MeV)附近,耦合計算的結果均能與MC 符合很好,而SN與MC的結果有明顯差異。
在模型1的基礎上,將外源沿x方向移至x=10cm 處,得到模型2,圖2示出了3種計算方式的結果。
圖2a中同樣可見SN方法的射線效應對計算結果的影響。在計算精度上,與模型1一致,耦合計算的結果與MC 方法的符合得很好,而SN的結果與MC方法的有明顯差異。
在模型3中,將考慮同時存在兩個中子源的情況,即假定同時存在模型1和2的中子源,圖3示出了3種計算方式的結果。
模型3的對比結果與前兩個一致,耦合計算的結果與MC 方法的符合得很好,而SN方法的結果與MC方法的有較大差異。
圖1 模型1中子注量率等值線對比Fig.1 Comparison of neutron fluence rate contour line of model 1
圖2 模型2中子注量率等值線對比Fig.2 Comparison of neutron fluence rate contour line of model 2
圖3 模型3中子注量率等值線對比Fig.3 Comparison of neutron fluence rate contour line of model 3
根據(jù)外源驅動次臨界系統(tǒng)中子輸運的特點,本文提出了以中子首次裂變?yōu)轳詈宵c的MC/SN耦合算法。在該耦合算法中,MC 方法用于模擬外源中子的早期輸運過程,直至發(fā)生裂變;對首次裂變中子的模擬,則采用SN方法。由于首次裂變中子源分布于整個裂變區(qū),且各向同性,采用SN方法已能保證計算精度。
不同分布源模型的計算表明,該耦合算法的計算結果與作為基準的MC 方法結果具有很好的一致性,而采用單一SN方法計算的結果則存在較大的偏差,初步驗證了基于中子首次裂變的MC/SN耦合算法在模擬局部源問題時的有效性和正確性。
[1] 沈華韻,李樹,衷斌,等.堆腔注水改進項對堆外探測器的影響分析[R].北京:北京應用物理與計算數(shù)學研究所,2012.
[2] LEWIS E E,MILLER W F.Computational methods of neutron transport[M].US:American Nuclear Society,Inc.,1993:156-358.
[3] SOOBY E,BATY A,BENE?O,et al.Candidate molten salt investigation for an accelerator driven subcritical core[J].Journal of Nuclear Materials,2013,440:298-303.
[4] LIN Z,CHEN J,GUO W,et al.The conceptual design of electron-accelerator-driven subcritical thorium molten salt system[J].Energy Procedia,2013,39:267-274.
[5] SIDDIQUE M T,KIM M H.Conceptual design study of Hyb-WT as fusion-fission hybrid reactor for waste transmutation[J].Annals of Nuclear Energy,2014,65:299-306.
[6] 裴鹿成,張孝澤.蒙特卡羅方法及其在粒子輸運問題中的應用[M].北京:科學出版社,1980:213-219.