點(diǎn)探測器敏感性關(guān)聯(lián)抽樣計算方法

李 瑞,付元光,鄧 力,許海波

(1.中物院高性能數(shù)值模擬軟件中心,北京 100088;2.北京應(yīng)用物理與計算數(shù)學(xué)研究所,北京 100094)

在反應(yīng)堆輻射屏蔽設(shè)計中,常需要通過改變屏蔽體尺寸、材料密度與組分實(shí)現(xiàn)屏蔽方案的優(yōu)化。在此過程中,掌握這些因素導(dǎo)致的輻射劑量率變化規(guī)律可有效提高屏蔽方案優(yōu)化效率。另外,核數(shù)據(jù)本身的不確定性會對輻射屏蔽模擬結(jié)果的準(zhǔn)確性造成影響,在屏蔽方案的評估中也需要加以考慮。當(dāng)不同屏蔽方案的差異較小時,上述計算輻射劑量率變化的問題可視為微擾問題。

反應(yīng)堆輻射屏蔽模擬一般采用離散縱標(biāo)(SN)方法與蒙特卡羅(MC)方法。其中SN方法屬于確定論方法,在計算微擾問題時可直接采用兩次獨(dú)立計算并做差的方式進(jìn)行。但是,由于SN方法對于復(fù)雜幾何處理能力有限,并存在射線效應(yīng)等問題,其計算精度一般低于MC方法。與SN方法相比,MC方法可處理比較復(fù)雜的幾何體,采用的連續(xù)能量核反應(yīng)截面也更適用于核數(shù)據(jù)擾動分析。但MC方法的計算量巨大,特別是對于輻射屏蔽這一類深穿透問題,需要使用降方差技巧才能得到統(tǒng)計收斂結(jié)果。然而,在微擾問題計算中,由于輻射劑量率本身的改變就很小,如果采用直接MC模擬方法計算,則對于其統(tǒng)計收斂性要求更加苛刻,導(dǎo)致時間開銷進(jìn)一步增加。因此,為排除隨機(jī)漲落的影響,MC方法一般在與參考系統(tǒng)相同的隨機(jī)序列下完成微擾系統(tǒng)的模擬,常見的方法包括關(guān)聯(lián)抽樣[1]與微分算符[2]。

近年來,采用伴隨通量作為重要性的降方差技巧在MC輻射屏蔽模擬中被廣泛研究,如針對單一目標(biāo)的CADIS算法[3]以及多目標(biāo)的FW-CADIS算法[4]。這些算法的實(shí)質(zhì)是基于粒子權(quán)重的賭/分裂以及偏倚抽樣的組合。另一方面,在MC輻射屏蔽模擬中,常采用點(diǎn)探測器計數(shù)方法提高計數(shù)率。

MCNP軟件[5]基于微分算符法,提供了體平均通量及其衍生計數(shù)的敏感性計算功能。Perel等[6]在此基礎(chǔ)上研究了點(diǎn)探測器計數(shù)的敏感性計算方法。與體平均通量這一類徑跡長度估計方法不同,點(diǎn)探測器計數(shù)采用了下次事件估計方法[7],需要額外考慮虛粒子的模擬過程。在Perel等[6]展示的算法中,需多次遍歷粒子的徑跡歷史貢獻(xiàn)。在采用權(quán)重賭/分裂技巧的屏蔽問題模擬中,次級粒子的數(shù)目很多且平均粒子隨機(jī)游動鏈較長,這種算法的性能開銷會很大。本工作基于關(guān)聯(lián)抽樣算法架構(gòu),利用單因素擾動系統(tǒng)將系統(tǒng)擾動因素局限在可通過一個擾動權(quán)重描述的簡單情形,在此基礎(chǔ)上以隨輸運(yùn)過程實(shí)時更新的方式代替頻繁的粒子歷史遍歷。另外,將擾動權(quán)作為粒子擴(kuò)展屬性也可避免對次級粒子徑跡的回溯需求,從而降低算法的實(shí)現(xiàn)難度。

1 算法介紹

關(guān)聯(lián)抽樣算法的重點(diǎn)在于采用與原系統(tǒng)相同的隨機(jī)序列模擬擾動系統(tǒng)的輸運(yùn)與計數(shù)統(tǒng)計過程。在形式上,擾動系統(tǒng)粒子隨機(jī)游動過程與原系統(tǒng)是完全相同的,統(tǒng)計差異僅體現(xiàn)在粒子權(quán)重上。本工作基于三維粒子輸運(yùn)模擬軟件JMCT開展算法研究[8]。

對于輻射屏蔽問題,點(diǎn)探測器計數(shù)估計量[6]可表示為:

(1)

(2)

式中:δ為克羅內(nèi)克函數(shù),表示僅當(dāng)虛粒子飛行方向朝向點(diǎn)探測器時觸發(fā)計數(shù);Ωd為虛粒子飛行方向矢量;Rd為點(diǎn)探測器的位置矢量;R為隨機(jī)游動鏈中粒子的位置矢量;Ps為粒子發(fā)生反應(yīng)時在點(diǎn)探測器方向發(fā)射粒子的概率,即指向概率;Σt為總反應(yīng)截面;s為虛粒子到達(dá)探測器過程中經(jīng)歷的徑跡段長度。

式(1)中的Pik(Ek)為粒子隨機(jī)游動鏈中徑跡段的產(chǎn)生概率,按照Perel等[6]對于徑跡段指標(biāo)的約定,每個徑跡段以碰撞或源抽樣事件開始,因此點(diǎn)探測器計數(shù)是在徑跡段開始處觸發(fā)的。由于虛粒子與實(shí)際碰撞或源抽樣過程是獨(dú)立的,因此影響第k個徑跡段點(diǎn)探測器計數(shù)的所有徑跡段k∈lin(j)并不包括k徑跡段本身[6]。因此,對于首個點(diǎn)探測器計數(shù),僅依賴徑跡段0,為保證形式上的一致性,將徑跡段指標(biāo)定義進(jìn)行簡單擴(kuò)展,擴(kuò)展之后的徑跡段概率為:

(3)

Q為輻射屏蔽計算中的外源分布,源粒子的位置屬性指標(biāo)為0,能量與方向指標(biāo)為1。由于粒子從R0飛行到R1過程中,其能量與方向是不變的,因此E1與Ω1分別對應(yīng)了第2個徑跡段開始之前粒子的能量與飛行方向。對于k>1的情況,位置屬性指標(biāo)與徑跡段的關(guān)系以此類推。

T為輸運(yùn)核,表示粒子飛行至下一次碰撞的概率,其形式為:

T(Rk,Ek,Ωk;Rk-1)=

(4)

Cx為各類反應(yīng)的碰撞核,如式(5)所示。

Cx(Rk-1,Ek,Ωk;Ek-1,Ωk-1)=

(5)

式中:Σt為宏觀總截面;Σx為宏觀反應(yīng)截面;νx為反應(yīng)發(fā)生后產(chǎn)生的次級粒子數(shù)目,主要用于中子誘發(fā)裂變、(n,xn)以及中子產(chǎn)光等反應(yīng)過程。

(6)

(7)

式(6)所示為權(quán)重技巧下的點(diǎn)探測器計數(shù)形式。當(dāng)系統(tǒng)存在擾動時,基于關(guān)聯(lián)抽樣計算點(diǎn)探測器計數(shù)對應(yīng)變化量Δ(Td)的過程為:

Δ(Td)=Est*[Td]-Est[Td]=

(8)

式中,上標(biāo)*表示擾動系統(tǒng)相關(guān)項。

式(8)對應(yīng)的隨機(jī)游動過程與式(6)類似,唯一的區(qū)別就是引入了擾動權(quán)乘子

(9)

(10)

(11)

根據(jù)式(8)～(11)可構(gòu)建點(diǎn)探測器擾動量的關(guān)聯(lián)抽樣算法。僅需要在原系統(tǒng)隨機(jī)游動模擬過程中更新并統(tǒng)計擾動權(quán)對應(yīng)的計數(shù)即可,算法流程如圖1所示。

圖1 點(diǎn)探測器擾動關(guān)聯(lián)抽樣計算流程Fig.1 Flowchart of point detector perturbation based on correlated sampling

另一方面,為完成點(diǎn)探測器擾動的模擬,對于每一個擾動因素,都需要建立與原系統(tǒng)結(jié)構(gòu)相同的統(tǒng)計計數(shù)。結(jié)構(gòu)相同是指點(diǎn)探測器的數(shù)目、能量區(qū)間劃分、時間區(qū)間劃分都相同。具有相同結(jié)構(gòu)的擾動系統(tǒng)計數(shù)可通過對原系統(tǒng)計數(shù)進(jìn)行克隆得到。所有擾動系統(tǒng)克隆計數(shù)的總內(nèi)存使用量M如式(12)所示。

(12)

式中:Np為擾動系統(tǒng)總數(shù);Nd為點(diǎn)探測器計數(shù)數(shù)目;Mi為每個點(diǎn)探測器計數(shù)的內(nèi)存使用量。

JMCT中允許在一次模擬中同時分析多個擾動因素,這些擾動因素包括核素截面與材料密度的變化,以及針對擾動因素的能量區(qū)間劃分。在本工作中,將某一能量區(qū)間內(nèi)具有單一核素與反應(yīng)道擾動或材料密度擾動的系統(tǒng)定義為單因素擾動系統(tǒng)。因此,一次MC模擬過程中,可存在多個單因素擾動系統(tǒng)。擾動系統(tǒng)對應(yīng)的計數(shù)內(nèi)存與內(nèi)部狀態(tài)是獨(dú)占的,有效地避免了不同擾動因素在計算流程上相互干擾。對于每個單因素擾動系統(tǒng),可通過擴(kuò)展固定數(shù)目為1的粒子屬性即可完成輸運(yùn)擾動權(quán)乘子的管理。此外,JMCT允許定義擾動所在幾何體與材料,從而限制擾動的空間范圍。

2 計算結(jié)果

利用勞倫斯利弗莫爾國家實(shí)驗(yàn)室(Lawrence Livermore National Laboratory, LLNL)的脈沖球形實(shí)驗(yàn)裝置對算法進(jìn)行正確性檢驗(yàn)。該實(shí)驗(yàn)探測D-T脈沖中子在半徑為20.96 cm的石墨球中輸運(yùn)后的時間泄漏譜,裝置如圖2所示。其中D-T中子源發(fā)生器通過孔道經(jīng)-x方向放置于x=-0.131 cm處。D-T中子能量介于13.20～15.11 MeV區(qū)間內(nèi),具體分布依賴于出射角度。實(shí)驗(yàn)中,NE-213探測器與-x軸呈30°角,距離石墨球心776 cm。文獻(xiàn)[9]給出了該模型相應(yīng)的MCNP程序[5]輸入文件。

圖2 脈沖球形實(shí)驗(yàn)裝置模型Fig.2 Model of pulsed sphere experiment device

JMCT首先基于該模型進(jìn)行了泄漏譜的計算,作為對比的實(shí)驗(yàn)值來自文獻(xiàn)[10]。JMCT整體計算結(jié)果與Perel等[6]工作符合良好,如圖3所示。

圖3 脈沖球形實(shí)驗(yàn)裝置時間泄漏譜計算結(jié)果Fig.3 Calculation result of time leakage spectrum for pulsed sphere experiment device

在此基礎(chǔ)上,對14.918～15.296 MeV區(qū)間內(nèi)的(n,n′1)反應(yīng)截面引入+1%的擾動,計算時間譜的變化。由于該實(shí)驗(yàn)裝置中D-T中子在石墨球中平均運(yùn)行2.9個自由程,因此擾動引起的時間泄漏譜變化主要集中在0～2次散射對應(yīng)的時間區(qū)間上。JMCT計算結(jié)果如圖4所示。其中由于(n,n′1)反應(yīng)截面的增加導(dǎo)致擾動系統(tǒng)的總截面變大,進(jìn)而造成點(diǎn)探測器直穿項計數(shù)降低。另外,散射截面的貢獻(xiàn)會由于散射截面的增加而變大,其中一次散射的貢獻(xiàn)最為顯著。將JMCT點(diǎn)探測器直穿項和一次散射峰的變化與Perel等[6]工作中微分算符計算結(jié)果進(jìn)行定量對比,相對偏差均在5%以內(nèi),如表1所列。

表1 時間泄漏譜擾動結(jié)果對比Table 1 Comparison of time leakage spectrum perturbation

圖4 脈沖球形實(shí)驗(yàn)裝置時間泄漏譜擾動計算結(jié)果Fig.4 Calculation result of time leakage spectrum perturbation for pulsed sphere experiment device

針對反應(yīng)堆輻射屏蔽問題,采用NUREG/CR-6115屏蔽基準(zhǔn)題[11]驗(yàn)證使用降方差技巧的點(diǎn)探測器敏感性計算結(jié)果。通過在壓力容器材料中引入+1%的核子密度擾動,計算壓力容器外中子能譜的擾動。為保證計算收斂,該模擬采用了CADIS降方差方法[3],其中伴隨通量來自JSNT[12]計算。在4×108粒子歷史規(guī)模下,JMCT壓力容器外表面中子能譜與文獻(xiàn)[11]中的DORT結(jié)果對比如圖5所示。

圖5 壓力容器外表面中子能譜計算結(jié)果Fig.5 Calculation result of neutron energy spectrum on outer surface of pressure vessel

這里采用直接MC模擬作為參考對比,關(guān)聯(lián)抽樣的計算結(jié)果如圖6所示。結(jié)果表明,兩者在趨勢上是符合的,即1～2 MeV能量區(qū)間的中子通量下降,而1 MeV能量以下的中子通量增加。導(dǎo)致這一現(xiàn)象的原因主要有2個:一方面是由于壓力容器材料對快中子的慢化作用導(dǎo)致高能中子慢化為低能中子;另一方面,壓力容器內(nèi)的Fe等核素對快中子的輻射俘獲較低能中子更加顯著。在該算例中,為避免兩次直接MC模擬的統(tǒng)計誤差掩蓋真實(shí)擾動,要求各能群統(tǒng)計標(biāo)準(zhǔn)差均低于0.02,兩次MC計算共耗時4 286.2CPU小時,關(guān)聯(lián)抽樣計算耗時4.1CPU小時。

圖6 壓力容器外表面中子能譜擾動計算結(jié)果Fig.6 Calculation result of neutron energy spectrum perturbation on outer surface of pressure vessel

3 結(jié)論

脈沖球形實(shí)驗(yàn)裝置計算結(jié)果表明,基于關(guān)聯(lián)抽樣的點(diǎn)探測器敏感性計算方法可達(dá)到與Perel等[6]工作中微分算符方法相當(dāng)?shù)挠嬎憔?說明了算法對輸運(yùn)與計數(shù)過程的微擾描述是正確的。在NUREG/CR-6115屏蔽基準(zhǔn)題模擬中,CADIS降方差技巧的使用一方面對源抽樣與輸運(yùn)過程進(jìn)行了偏倚,另一方面也會在權(quán)重輪盤賭過程中產(chǎn)生大量的次級粒子。本工作利用單因素擾動系統(tǒng)的概念將擾動因素進(jìn)行了規(guī)范化,在每個單因素擾動系統(tǒng)內(nèi),以粒子擴(kuò)展屬性的方式實(shí)現(xiàn)了擾動權(quán)乘子的更新。這一方法解決了權(quán)重技巧下輸運(yùn)與計數(shù)過程的微擾描述問題,并且利用次級粒子的存庫與獲取順序避免了敏感性統(tǒng)計中對粒子歷史進(jìn)行頻繁的回溯,從而降低了算法的運(yùn)行時間開銷。