黃金龍,曹良志,賀清明,秦 帥,吳宏春

(西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院,陜西 西安 710049)

核數(shù)據(jù)是核反應(yīng)堆中子學(xué)計(jì)算的重要輸入?yún)?shù),其敏感性分析結(jié)果可反映對(duì)響應(yīng)參數(shù)敏感的核數(shù)據(jù),進(jìn)而對(duì)核反應(yīng)堆堆芯設(shè)計(jì)和不確定度量化[1]等有重要的指導(dǎo)意義。反復(fù)裂變幾率(IFP)方法[2]被眾多蒙特卡羅程序用于k特征值對(duì)核數(shù)據(jù)的敏感性分析[3-6],該方法的主要缺點(diǎn)是,內(nèi)存占用與每代模擬粒子數(shù)呈正比,當(dāng)每代模擬粒子數(shù)較大時(shí),內(nèi)存占用甚至不可接受。Kiedrowski在MCNP中提出稀疏矩陣[7]存儲(chǔ)方式,以降低IFP方法的內(nèi)存使用。數(shù)據(jù)表明,采用稀疏矩陣方法可降低10～100倍的內(nèi)存,但內(nèi)存占用量依舊較大。Perfetti[8]提出CLUTCH方法[9-10]可在不損失結(jié)果精度的前提下,顯著降低內(nèi)存占用。CLUTCH方法需要?jiǎng)澐志W(wǎng)格來(lái)統(tǒng)計(jì)權(quán)重函數(shù),且網(wǎng)格的大小[8]推薦為1 cm3,對(duì)于大規(guī)模問(wèn)題(如壓水堆全堆問(wèn)題),需要?jiǎng)澐执罅烤W(wǎng)格來(lái)統(tǒng)計(jì)權(quán)重函數(shù),在網(wǎng)格數(shù)量較多的情況下,保證權(quán)重函數(shù)收斂所需的粒子數(shù)是不可接受的。因此,本文將使用函數(shù)展開(kāi)計(jì)數(shù)(FET)方法[11]對(duì)權(quán)重函數(shù)進(jìn)行統(tǒng)計(jì),提高統(tǒng)計(jì)權(quán)重函數(shù)的精度和效率,從而實(shí)現(xiàn)對(duì)大規(guī)模問(wèn)題的高效率、高精度的連續(xù)能量核數(shù)據(jù)敏感性分析。

本文在蒙特卡羅程序NECP-MCX[12]上實(shí)現(xiàn)基于網(wǎng)格和函數(shù)展開(kāi)計(jì)數(shù)方法的CLUTCH方法(CLUTCH-Mesh和CLUTCH-FET),以進(jìn)行k特征值對(duì)核數(shù)據(jù)敏感性分析。以IFP方法[13]的計(jì)算結(jié)果為參考解,在Godiva、Flattop和AP1000全堆問(wèn)題上對(duì)CLUTCH-Mesh和CLUTCH-FET方法進(jìn)行驗(yàn)證,并比較兩種方法的計(jì)算精度和效率。

1 理論方法

基于一階微擾理論,k特征值對(duì)核數(shù)據(jù)的靈敏度系數(shù)S可表示為:

Sk,Σx(r,E)=

(1)

其中:Σx(r,E)為x反應(yīng)道在位置r處能量為E的宏觀截面;F為裂變項(xiàng)算子;S為散射項(xiàng)算子;T為碰撞項(xiàng)算子;ψ為中子通量;ψ*為共軛通量;〈,〉表示在相空間中積分。

由式(1)可看出,靈敏度系數(shù)是反應(yīng)率和共軛通量的乘積,反應(yīng)率在蒙特卡羅計(jì)數(shù)中容易統(tǒng)計(jì),共軛通量可用IFP和CLUTCH方法求解。

1.1 CLUTCH方法

基于Contributon理論[14],相空間中某點(diǎn)P0處的共軛通量表達(dá)式為:

(2)

其中:Qs(P0)為粒子進(jìn)入P0時(shí)的權(quán)重;χ(r,E)為r處出射能量為E的裂變譜;ν(E′)為入射能量為E′的裂變中子數(shù);Σf(r,E′)為位置r處能量為E′的裂變截面;Ω′和Ω分別為中子的入射方向和出射方向;ψ(r,Ω′,E′|P0)為從P0點(diǎn)產(chǎn)生的中子在(r,Ω′,E′)處產(chǎn)生的中子通量;G(P0→r)為轉(zhuǎn)移函數(shù),表示從P0處轉(zhuǎn)移到r處的中子通量;I*(r)為r處的權(quán)重函數(shù)。

權(quán)重函數(shù)I*(r)可根據(jù)未歸一的裂變譜的靈敏度系數(shù)求解,表達(dá)式為:

(3)

其中:D為式(1)的分母項(xiàng);式(3)的分母表示初始代r處產(chǎn)生的裂變中子數(shù),分子表示初始代由r處產(chǎn)生的裂變中子,經(jīng)過(guò)幾個(gè)蒙特卡羅代的模擬后,在漸近代產(chǎn)生的中子數(shù),因此,I*(r)為初始代r處平均產(chǎn)生的1個(gè)裂變中子在漸近代所產(chǎn)生的中子數(shù)。由于I*(r)只與空間位置有關(guān),因此,通?；诳臻g網(wǎng)格對(duì)I*(r)進(jìn)行統(tǒng)計(jì)且可在非活躍代統(tǒng)計(jì)得到。空間網(wǎng)格需要足夠精細(xì),以準(zhǔn)確描述權(quán)重函數(shù)的空間分布,同時(shí)需要模擬足夠多的粒子,以保證權(quán)重函數(shù)收斂。Perfetti[8]推薦,網(wǎng)格大小約為1 cm3,每個(gè)包含燃料的網(wǎng)格有1 000個(gè)非活躍代中子歷史,可保證權(quán)重函數(shù)足夠精細(xì)和收斂。以AP1000全堆問(wèn)題為例,為保證網(wǎng)格大小約為1 cm3,需采用320×320×400的網(wǎng)格,為保證每個(gè)包含燃料的網(wǎng)格有1 000個(gè)非活躍代中子歷史,若每代模擬的中子數(shù)為106,則需要模擬的非活躍代數(shù)為28 581代,計(jì)算代價(jià)很大。因此,本文將使用函數(shù)展開(kāi)計(jì)數(shù)方法對(duì)權(quán)重函數(shù)進(jìn)行統(tǒng)計(jì)。

1.2 函數(shù)展開(kāi)計(jì)數(shù)方法

1) FET系數(shù)求解

蒙特卡羅統(tǒng)計(jì)目標(biāo)量φ(ξ)可用任意正交完備的基函數(shù)展開(kāi)[11],表示為:

(4)

其中:ξ為歸一化后的相空間位置;n為函數(shù)展開(kāi)階數(shù);an為n階函數(shù)展開(kāi)系數(shù);φn為n階展開(kāi)函數(shù);ρ(ξ)為ξ處的權(quán)重函數(shù)。

根據(jù)基函數(shù)的正交性,可求解an:

(5)

其中,Γ為正交區(qū)域。

基于碰撞估計(jì)器,可得到an的估計(jì)表達(dá)式:

(6)

基于函數(shù)展開(kāi)計(jì)數(shù)方法對(duì)式(3)進(jìn)行統(tǒng)計(jì)時(shí)需分兩步:1) 初始代用函數(shù)展開(kāi)計(jì)數(shù)方法統(tǒng)計(jì)裂變中子的分布,即統(tǒng)計(jì)分母,并保存函數(shù)展開(kāi)系數(shù);2) 由于初始代的裂變中子分布已得到,即式(3)的分母已知,再用函數(shù)展開(kāi)計(jì)數(shù)方法對(duì)式(3)進(jìn)行統(tǒng)計(jì)。

2) 勒讓德多項(xiàng)式

勒讓德多項(xiàng)式[15]是在區(qū)間[-1,1]上的正交多項(xiàng)式,n階勒讓德多項(xiàng)式可表示為:

(7)

其在區(qū)間[-1,1]上的正交性表示為:

(8)

以基于勒讓德函數(shù)對(duì)三維空間變量φ(x,y,z)的函數(shù)展開(kāi)為例進(jìn)行說(shuō)明。φ(x,y,z)可表示為:

(9)

其中:i,j,k分別為x,y,z方向展開(kāi)系數(shù)索引;I,J,K分別為x,y,z方向展開(kāi)階數(shù);fi,j,k為函數(shù)展開(kāi)系數(shù),根據(jù)式(8)所示勒讓德函數(shù)的正交性,fi,j,k可根據(jù)式(10)求得。

(10)

基于碰撞估計(jì)器,可得到fi,j,k的估計(jì)式:

(11)

在非活躍代基于勒讓德多項(xiàng)式對(duì)權(quán)重函數(shù)進(jìn)行函數(shù)展開(kāi)計(jì)數(shù),將得到的各階展開(kāi)系數(shù)代入式(9),在活躍代即可得到空間各位置處的權(quán)重函數(shù),將得到的權(quán)重函數(shù)代入式(2),用于共軛通量的求解,進(jìn)而由式(1)計(jì)算得到靈敏度系數(shù)。

2 數(shù)值結(jié)果

基于上述理論方法在NECP-MCX中實(shí)現(xiàn)基于網(wǎng)格和函數(shù)展開(kāi)計(jì)數(shù)的CLUTCH方法,用于k特征值對(duì)核數(shù)據(jù)的敏感性分析,并與IFP方法的計(jì)算結(jié)果比較,對(duì)比CLUTCH-Mesh和CLUTCH-FET方法的計(jì)算精度以及效率,其中計(jì)算效率用品質(zhì)因子(FOM)衡量,品質(zhì)因子的定義如式(12)所示。

(12)

其中:R為蒙特卡羅統(tǒng)計(jì)量的相對(duì)統(tǒng)計(jì)標(biāo)準(zhǔn)差;T為蒙特卡羅計(jì)算時(shí)間。品質(zhì)因子越大,則計(jì)算效率越高。本文選取Godiva、Flattop和AP1000全堆問(wèn)題進(jìn)行驗(yàn)證,使用ACE格式的ENDF/BⅦ.0數(shù)據(jù)庫(kù)。

2.1 Godiva問(wèn)題

Godiva基準(zhǔn)題[16]是一個(gè)均勻裸球裝置,材料為高濃鈾,半徑為8.741 cm,如圖1所示。靈敏度系數(shù)計(jì)算采用2 000個(gè)蒙特卡羅代,其中活躍代設(shè)為1 000代,每代模擬粒子數(shù)為105,塊的大小設(shè)置為10代。由于Godiva屬于一維問(wèn)題,因此只需在半徑方向上對(duì)權(quán)重函數(shù)進(jìn)行統(tǒng)計(jì)。

圖1 Godiva基準(zhǔn)題Fig.1 Model for Godiva problem

選用不同展開(kāi)階數(shù)的勒讓德多項(xiàng)式對(duì)權(quán)重函數(shù)進(jìn)行統(tǒng)計(jì),得到不同展開(kāi)階數(shù)下CLUTCH-FET方法計(jì)算的靈敏度系數(shù),如表1所列,其中Ave和RSD分別為靈敏度系數(shù)均值和相對(duì)統(tǒng)計(jì)漲落。

表1 Godiva問(wèn)題中不同展開(kāi)階數(shù)下CLUTCH-FET方法計(jì)算的靈敏度系數(shù)Table 1 Sensitivity coefficient of CLUTCH-FET method under different expansion orders for Godiva problem

(13)

其中:m為反應(yīng)道序號(hào);M為反應(yīng)道總數(shù),即M=10;dn,m為勒讓德階數(shù)取n階時(shí)第m個(gè)反應(yīng)道的靈敏度系數(shù)計(jì)算結(jié)果與參考解之間的差異,表示為:

(14)

其中:Aven,m和RSDn,m分別為展開(kāi)系數(shù)階數(shù)取n的第m個(gè)反應(yīng)道的靈敏度系數(shù)的均值和相對(duì)統(tǒng)計(jì)漲落;AveIFP,m和RSDIFP,m分別為IFP方法的第m個(gè)反應(yīng)道的靈敏度系數(shù)的均值和相對(duì)統(tǒng)計(jì)漲落。

圖2 Godiva問(wèn)題中不同展開(kāi)階數(shù)下CLUTCH-FET方法計(jì)算的靈敏度系數(shù)與參考解間的差異Fig.2 Difference between sensitivity coefficients calculated by CLUTCH-FET under different expansion orders and reference for Godiva problem

CLUTCH-Mesh方法在半徑方向上劃分16等分作為統(tǒng)計(jì)權(quán)重函數(shù)的網(wǎng)格。CLUTCH-FET方法采用4階勒讓德多項(xiàng)式對(duì)權(quán)重函數(shù)在半徑方向上的分布進(jìn)行統(tǒng)計(jì),兩種方法求得的權(quán)重函數(shù)分布如圖3所示,IFP、CLUTCH-Mesh和CLUTCH-FET方法計(jì)算得到的靈敏度系數(shù)列于表2,FOM值示于圖4。

表2 Godiva問(wèn)題靈敏度系數(shù)計(jì)算結(jié)果Table 2 Sensitivity coefficient result for Godiva problem

圖3 Godiva問(wèn)題權(quán)重函數(shù)分布Fig.3 Distribution of importance weighting functions for Godiva problem

圖4 Godiva問(wèn)題FOM值對(duì)比Fig.4 Comparison of FOM values for Godiva problem

圖3表明,基于網(wǎng)格和基于函數(shù)展開(kāi)計(jì)數(shù)方法計(jì)算得到的權(quán)重函數(shù)符合得較好,隨著半徑的增大,權(quán)重函數(shù)的數(shù)值逐漸減小。表2顯示,CLUTCH-Mesh和CLUTCH-FET方法的計(jì)算結(jié)果與IFP方法計(jì)算結(jié)果的相對(duì)偏差基本小于1%,表明了兩種方法的準(zhǔn)確性。通過(guò)對(duì)比圖4所示FOM值,可看出, CLUTCH-Mesh和CLUTCH-FET方法計(jì)算得到的靈敏度系數(shù)的FOM值相近,CLUTCH-Mesh方法計(jì)算得到的FOM值略高于CLUTCH-FET方法,兩種CLUTCH方法的FOM值普遍高于IFP方法,表明CLUTCH-Mesh和CLUTCH-FET方法相比于IFP方法計(jì)算效率有較大提升。同一方法、不同核數(shù)據(jù)的FOM值也有較大差異,這與不同反應(yīng)道發(fā)生的概率有關(guān)。

2.2 Flattop問(wèn)題

Flattop基準(zhǔn)題[16]布置如圖5所示,裝置中心材料是钚合金,反射層材料是鈾,堆芯半徑為4.533 2 cm,反射層外徑為24.142 cm。CLUTCH-Mesh方法采用半徑方向均分48份一維球狀網(wǎng)格,CLUTCH-FET方法在半徑方向上全局和分段進(jìn)行函數(shù)展開(kāi)計(jì)數(shù),其中全局函數(shù)展開(kāi)計(jì)數(shù)表示在半徑方向上對(duì)[0 cm,24.142 cm]進(jìn)行勒讓德函數(shù)展開(kāi),展開(kāi)階數(shù)分別取5、10、15、20階,分段函數(shù)展開(kāi)計(jì)數(shù)表示在半徑方向上分別對(duì)[0 cm,4.533 2 cm]和[4.533 2 cm,24.142 cm]進(jìn)行勒讓德函數(shù)展開(kāi),得到兩套函數(shù)展開(kāi)系數(shù),展開(kāi)階數(shù)均為4階,其他計(jì)算參數(shù)與Godiva問(wèn)題相同。得到的權(quán)重函數(shù)沿半徑的分布如圖6所示,3種方法計(jì)算的靈敏度系數(shù)列于表3,FOM值示于圖7。

表3 Flattop問(wèn)題靈敏度系數(shù)計(jì)算結(jié)果Table 3 Sensitivity coefficient result for Flattop problem

圖5 Flattop基準(zhǔn)題Fig.5 Model for Flattop problem

圖6 Flattop問(wèn)題權(quán)重函數(shù)分布Fig.6 Distribution of importance weighting functions for Flattop problem

圖7 Flattop問(wèn)題FOM值對(duì)比Fig.7 Comparison of FOM values for Flattop problem

圖6表明,權(quán)重函數(shù)在半徑方向上的變化達(dá)到上百倍,若不進(jìn)行分段函數(shù)展開(kāi)計(jì)數(shù),展開(kāi)階數(shù)取20階,仍無(wú)法準(zhǔn)確統(tǒng)計(jì)出權(quán)重函數(shù)的分布,因此需要進(jìn)行分段函數(shù)展開(kāi)計(jì)數(shù),才能較好地?cái)M合權(quán)重函數(shù)的分布。結(jié)果顯示,分兩段進(jìn)行函數(shù)展開(kāi)計(jì)數(shù),展開(kāi)階數(shù)取4階時(shí),便可得到較好的權(quán)重函數(shù)分布結(jié)果。因此,當(dāng)權(quán)重函數(shù)變化較大時(shí),分段進(jìn)行函數(shù)展開(kāi)計(jì)數(shù)較全局函數(shù)展開(kāi)計(jì)數(shù)的效果更好。一般函數(shù)展開(kāi)計(jì)數(shù)的分段選在材料分界處,保證每段函數(shù)展開(kāi)計(jì)數(shù)區(qū)域內(nèi)權(quán)重函數(shù)的梯度不大,以得到更精確的權(quán)重函數(shù)。圖6還表明,分段進(jìn)行權(quán)重函數(shù)展開(kāi)計(jì)數(shù)時(shí),在區(qū)間邊界處會(huì)出現(xiàn)間斷,但由于權(quán)重函數(shù)僅在燃料區(qū)域有意義,且研究結(jié)果[8]表明靈敏度系數(shù)的計(jì)算對(duì)權(quán)重函數(shù)的精度要求為10%～20%,因此權(quán)重函數(shù)在區(qū)間邊界處的間斷性對(duì)靈敏度系數(shù)的計(jì)算影響較小。由表3可看出,基于CLUTCH-Mesh方法和CLUTCH-FET方法計(jì)算的靈敏度系數(shù)與IFP結(jié)果的相對(duì)偏差均小于1%,且兩種方法的計(jì)算偏差數(shù)值相近,表明兩種方法具有較高的精度且精度相近。圖7表明,基于CLUTCH-Mesh方法和CLUTCH-FET方法計(jì)算得到的靈敏度系數(shù)的FOM值相近,且普遍高于IFP方法的FOM值,表明兩種CLUTCH方法的計(jì)算效率比IFP方法更高。

2.3 AP1000全堆問(wèn)題

AP1000堆芯布置如圖8所示,采用15×15的燃料布置,共157個(gè)燃料組件,控制棒處于全提狀態(tài)。在統(tǒng)計(jì)權(quán)重函數(shù)時(shí),CLUTCH-Mesh方法采用340×340×400的立方體網(wǎng)格,以保證網(wǎng)格的大小約為1 cm3,CLUTCH-FET方法對(duì)每個(gè)燃料組件分別進(jìn)行函數(shù)展開(kāi)計(jì)數(shù),即采用157套函數(shù)展開(kāi)計(jì)數(shù),每個(gè)燃料組件內(nèi),x、y方向各選取4階勒讓德函數(shù)展開(kāi)計(jì)數(shù),由于燃料區(qū)域的材料在z方向沒(méi)有變化,因此在z方向不分段,由于z方向燃料區(qū)域較長(zhǎng),適當(dāng)選取高階展開(kāi)系數(shù),選取10階勒讓德函數(shù)展開(kāi)計(jì)數(shù)。每代模擬的粒子數(shù)為106,其他計(jì)算參數(shù)同Godiva問(wèn)題。IFP、CLUTCH-Mesh和CLUTCH-FET方法計(jì)算的靈敏度系數(shù)列于表4,FOM值對(duì)比示于圖9。

表4 AP1000全堆問(wèn)題靈敏度系數(shù)計(jì)算結(jié)果Table 4 Sensitivity coefficient result for AP1000 problem

圖8 AP1000堆芯布置Fig.8 AP1000 core layout

圖9 AP1000問(wèn)題FOM值對(duì)比Fig.9 Comparison of FOM values for AP1000 problem

表4表明,在當(dāng)前計(jì)算條件下,CLUTCH-Mesh方法計(jì)算得到的238U-σγ、1H-σelas、10B-σt、235U-σγ和1H-σγ的靈敏度系數(shù)的計(jì)算相對(duì)偏差較大,均大于1%,1H-σγ的甚至達(dá)到-3.04%,這是由于AP1000全堆問(wèn)題規(guī)模較大,需要?jiǎng)澐肿銐蚓?xì)的網(wǎng)格才能準(zhǔn)確描述權(quán)重函數(shù)的空間分布,網(wǎng)格數(shù)的增多使得每個(gè)網(wǎng)格統(tǒng)計(jì)得到收斂的權(quán)重函數(shù)所需的粒子數(shù)增多,根據(jù)前文結(jié)果,在當(dāng)前網(wǎng)格設(shè)置下,每代模擬106個(gè)粒子,需要28 581代的非活躍代模擬才能使權(quán)重函數(shù)收斂,而本次計(jì)算僅采用1 000代非活躍代,權(quán)重函數(shù)未收斂,導(dǎo)致靈敏度系數(shù)的計(jì)算相對(duì)偏差較大。而在當(dāng)前計(jì)算條件設(shè)置下,CLUTCH-FET方法計(jì)算得到的靈敏度系數(shù)的相對(duì)偏差均小于1%,表明CLUTCH-FET方法具有較高的精度。對(duì)比圖9所示的FOM值,發(fā)現(xiàn)CLUTCH-FET方法的FOM值普遍高于CLUTCH-Mesh和IFP方法,CLUTCH-Mesh方法的FOM值最低,這是由于CLUTCH-Mesh方法中網(wǎng)格數(shù)達(dá)到千萬(wàn)級(jí)別,對(duì)于數(shù)據(jù)量在千萬(wàn)級(jí)別的數(shù)組進(jìn)行規(guī)約和運(yùn)算非常耗時(shí),導(dǎo)致CLUTCH-Mesh方法的FOM值比CLUTCH-FET方法小。相較于IFP方法和CLUTCH-Mesh方法,CLUTCH-FET方法的FOM值分別最大提高了5.2倍和6.0倍。

3 結(jié)論

本文分別基于網(wǎng)格和函數(shù)展開(kāi)計(jì)數(shù)方法對(duì)CLUTCH方法中的權(quán)重函數(shù)進(jìn)行了統(tǒng)計(jì),函數(shù)展開(kāi)計(jì)數(shù)選取的基函數(shù)為勒讓德多項(xiàng)式,在Godiva、Flattop和AP1000全堆問(wèn)題上進(jìn)行了驗(yàn)證,并將兩種CLUTCH方法(CLUTCH-Mesh和CLUTCH-FET)計(jì)算得到的靈敏度系數(shù)與IFP方法進(jìn)行了比較。結(jié)果表明,對(duì)于Godiva和Flattop問(wèn)題,CLUTCH-Mesh和CLUTCH-FET方法具有和IFP方法相當(dāng)?shù)木?且相較于IFP方法,兩種CLUTCH方法的計(jì)算效率有顯著提升。對(duì)于AP1000全堆問(wèn)題,因該問(wèn)題規(guī)模較大,統(tǒng)計(jì)權(quán)重函數(shù)所需的網(wǎng)格數(shù)較多,導(dǎo)致權(quán)重函數(shù)難收斂,從而使得CLUTCH-Mesh方法的計(jì)算精度下降,表明CLUTCH-Mesh方法不適用于大規(guī)模問(wèn)題的靈敏度系數(shù)求解。而在選取合適的函數(shù)展開(kāi)階數(shù)的條件下,CLUTCH-FET方法的計(jì)算精度與IFP方法相當(dāng),計(jì)算效率有所提升,其品質(zhì)因子較IFP方法最大提高了5.2倍,較CLUTCH-Mesh方法最大提高了6.0倍,表明CLUTCH-FET方法適用于大規(guī)模問(wèn)題的靈敏度系數(shù)求解。Flattop問(wèn)題的權(quán)重函數(shù)分布結(jié)果表明,分段進(jìn)行函數(shù)展開(kāi)計(jì)數(shù)的效果較不分段的效果好。