王新宇，張 斌，陳義學(xué)

(華北電力大學(xué) 核科學(xué)與工程學(xué)院，北京 102206)

離散縱標(biāo)(SN)方法是屏蔽計(jì)算中常用的算法之一[1]。在離散縱標(biāo)屏蔽計(jì)算中，因?qū)αＷ虞斶\(yùn)方程中空間變量和角度變量的近似處理，不可避免地存在離散誤差，特別是屏蔽計(jì)算中角通量密度在空間或角度上的分布不光滑甚至間斷的情況。射線效應(yīng)是角度離散化造成的標(biāo)通量密度呈現(xiàn)空間震蕩分布的現(xiàn)象，在帶有孔道、孤立點(diǎn)源的屏蔽問題中尤為明顯[2-3]。首次碰撞源方法[4-7]通過解析方法或高階方法求解點(diǎn)源產(chǎn)生的未碰撞角通量密度，計(jì)算出每個(gè)空間網(wǎng)格在離散方向上的首次碰撞源作為SN計(jì)算中的分布源。這樣即使在弱散射情況下，粒子也能運(yùn)動到整個(gè)求解空間從而減小角度離散誤差。首次碰撞源方法雖能在一定程度上緩解射線效應(yīng)，但無法解決空間離散導(dǎo)致的數(shù)值擴(kuò)散問題，同時(shí)首次碰撞源方法難以處理二次射線效應(yīng)[8]。對空間變量的離散易導(dǎo)致未碰撞通量密度的數(shù)值解出現(xiàn)數(shù)值擴(kuò)散效應(yīng)，具體表現(xiàn)為數(shù)值模擬中粒子的傳遞偏離真實(shí)的方向或粒子在傳遞過程中不能正確地衰減[9]。離散縱標(biāo)屏蔽計(jì)算中的空間變量離散格式主要可分為有限差分格式[10]、短特征線格式[11]和間斷有限元格式[12]等，其中一些空間離散格式雖具備良好的魯棒性，但不能緩解因角度離散導(dǎo)致的射線效應(yīng)。

本文基于射線追蹤研究一種半解析的多次碰撞源方法，以緩解射線效應(yīng)和數(shù)值擴(kuò)散對屏蔽計(jì)算精度的影響。通過計(jì)算在選定區(qū)域內(nèi)粒子發(fā)生多次碰撞的通量密度，將粒子每次碰撞看作一個(gè)散射的過程。由每次碰撞通量密度(包括未碰撞通量密度)和散射截面計(jì)算出選定區(qū)域網(wǎng)格的多次碰撞源，將此分布源作為固定源對整個(gè)模型進(jìn)行離散縱標(biāo)輸運(yùn)計(jì)算。然后將每次碰撞通量密度(包括未碰撞通量密度)和輸運(yùn)求解的碰撞通量密度求和，得到整個(gè)計(jì)算模型的通量密度分布。

1 理論方法

1.1 多次碰撞源方法

穩(wěn)態(tài)單群固定源無裂變輸運(yùn)方程如下：

(1)

Lψ=Sψ+Q

(2)

L為輸運(yùn)算子：

(3)

S為散射矩陣：

(4)

將整個(gè)模型分為區(qū)域A和區(qū)域B兩部分，假設(shè)選定計(jì)算粒子在區(qū)域A中碰撞多次的通量密度，區(qū)域A需包含源區(qū)。由此可得：

L(ψA+ψB)=(SA+SB)(ψA+ψB)+Q

(5)

其中：ψA和ψB分別為區(qū)域A和區(qū)域B中的通量密度，兩者之和表示整個(gè)模型內(nèi)的角通量密度分布情況；SA和SB分別為區(qū)域A和區(qū)域B中的散射矩陣，兩者之和表示整個(gè)模型內(nèi)的散射矩陣的分布。在區(qū)域A內(nèi)的ψB根據(jù)其物理意義均等于0，此時(shí)源區(qū)對區(qū)域A的通量密度占主要貢獻(xiàn)，將式(5)進(jìn)行化簡后可得到：

LψA=SψA+Q

(6)

式(6)描述粒子在區(qū)域A運(yùn)動的過程，區(qū)域B的通量密度來自區(qū)域A的通量密度的貢獻(xiàn)，即在后續(xù)計(jì)算中將區(qū)域A看作模型的源區(qū)，式(5)、(6)聯(lián)立后可得式(7)：

LψB=SBψA+(SA+SB)ψB=SBψA+SψB

(7)

由式(6)可得：

(8)

(9)

LψSN=SψSN+Q(n+1)

(10)

最終，由兩步輸運(yùn)結(jié)果之和得到總通量密度，即：

(11)

多次碰撞源方法的主要思想在于將粒子每次發(fā)生碰撞的過程獨(dú)立處理，與標(biāo)準(zhǔn)的SN方法的源迭代求解流程有所區(qū)別。在計(jì)算粒子每次發(fā)生碰撞的過程時(shí)采用半解析的射線追蹤方法[13]，盡量減少因角度離散和空間離散引起的離散誤差，更準(zhǔn)確地描述粒子在模型內(nèi)的輸運(yùn)行為。多次碰撞源方法具有一些積分輸運(yùn)方法的性質(zhì)。當(dāng)選定的計(jì)算區(qū)域和碰撞次數(shù)均取極限值，即是計(jì)算整個(gè)模型內(nèi)粒子碰撞無窮次后的通量密度分布情況，與標(biāo)準(zhǔn)的SN方法的源迭代求解流程所描述的物理過程一致。

1.2 求解流程

標(biāo)準(zhǔn)的SN方法求解方法是對群內(nèi)散射源進(jìn)行迭代計(jì)算，迭代方程如下：

(12)

(13)

其中：D為離散角度轉(zhuǎn)化為矩的算符；φ(l)為通量矩；ε為迭代的收斂準(zhǔn)則。令迭代初始值φ(0)=0，ψ(n)即為至多經(jīng)歷n次碰撞的粒子角通量密度。在源迭代過程中得到的ψ(0)，ψ(1)，…，ψ(n)分別為粒子未經(jīng)碰撞的角通量密度，經(jīng)1次碰撞的角通量密度，…，經(jīng)n次碰撞的角通量密度。

選取Kobayashi基準(zhǔn)題模型Ⅱ[14]分析迭代過程中模型內(nèi)的通量密度分布情況。Kobayashi基準(zhǔn)題模型Ⅱ的幾何結(jié)構(gòu)如圖1所示。源強(qiáng)S和截面信息列于表1。

圖1 Kobayashi基準(zhǔn)題模型Ⅱ幾何示意圖Fig.1 Geometry of Kobayashi benchmark problem model Ⅱ

表1 源強(qiáng)與截面信息
Table 1 Source strength and cross section

區(qū)域S/(cm-3·s-1)Σt/cm-1Σs/cm-1源區(qū)10.10.05真空區(qū)010-40.5×10-4屏蔽區(qū)00.10.05

為盡可能減少離散誤差的影響，將模型劃分為步長為1 cm的60×100×60的空間網(wǎng)格集合，空間離散格式選擇置零修正的菱形差分格式，采用100階勒讓德-切比雪夫求積組進(jìn)行計(jì)算，圖2為當(dāng)源迭代過程進(jìn)行到第n次時(shí)，模型內(nèi)ψ(n)的分布情況。

從分布形狀上看，在初始幾次迭代中標(biāo)通量密度在空間上分布呈現(xiàn)非物理振蕩，而隨后的散射過程會使整個(gè)模型空間內(nèi)的標(biāo)通量密度的分布趨于光滑，標(biāo)通量密度的分布是隨著迭代的進(jìn)行越來越廣且各向同性越來越強(qiáng)的，其非均勻性逐漸減弱。數(shù)值結(jié)果表明，對于帶有孔道的屏蔽問題，前幾次迭代的標(biāo)通量密度數(shù)值上相對較大且主要分布在孔道內(nèi)，隨著迭代的進(jìn)行標(biāo)通量密度逐漸減小。因此僅在計(jì)算未碰撞通量密度和前n次碰撞通量密度時(shí)采用半解析的射線追蹤方法，后續(xù)的碰撞通量密度采用SN方法求解：

圖2 Kobayashi基準(zhǔn)題模型Ⅱ內(nèi)通量密度分布示意圖Fig.2 Flux density distribution in Kobayashi benchmark problem model Ⅱ

(14)

其中：ψ(i)為粒子第i次碰撞的通量密度；Q(i)為粒子第i次碰撞時(shí)的源項(xiàng)。射線追蹤方法主要在于計(jì)算粒子在幾何空間中從ro點(diǎn)穿行到rp點(diǎn)的徑跡長度，即光學(xué)距離τ(ro,rp)，其表達(dá)式如下：

(15)

其中：μ為粒子運(yùn)動方向的角余弦；Σt(ro+μΩrp)為ro+μΩo→p位置點(diǎn)的總截面，Ωo→p是從ro點(diǎn)穿行到rp點(diǎn)的單位向量。根據(jù)角通量密度的球諧矩定義：

(16)

碰撞角通量密度的球諧矩可由式(14)和式(16)確定：

(17)

采用多次碰撞源方法計(jì)算粒子在整個(gè)模型空間內(nèi)碰撞無窮次的情況時(shí)，總通量密度ψtotal的迭代求解流程如式(18)所示，表現(xiàn)出多次碰撞源方法的積分性質(zhì)。

Lψ(1)=Q,Sψ(2)=Sψ(1),…,Lψ(n)=Sψ(n-1)

(18)

2 基準(zhǔn)驗(yàn)證

2.1 自設(shè)屏蔽問題

圖3 自設(shè)屏蔽問題幾何模型Fig.3 Geometric model of self-designed shielding problem

該自設(shè)屏蔽問題幾何模型如圖3所示，尺寸為50 cm×50 cm×50 cm，邊界條件均為真空邊界；源區(qū)位于模型中心邊長為2 cm的立方體，中子源強(qiáng)為1 cm-3·s-1；源區(qū)外圍是一層尺寸為4 cm×4 cm×4 cm的高散射比材料；屏蔽層內(nèi)分布著4個(gè)總截面較大的立方體塊，尺寸為1 cm×1 cm×1 cm。具體的截面信息列于表2。

表2 自設(shè)屏蔽問題的截面信息Table 2 Cross section of self-designed shielding problem

該問題源區(qū)的四周分布有總截面相對較大的高散射比單元，因射線效應(yīng)形成的粒子束穿過這些單元時(shí)可能會導(dǎo)致新的或次要的射線效應(yīng)，采用首次碰撞源方法修正不能消除由這些單元產(chǎn)生的二次射線效應(yīng)，需更高次的碰撞源修正方法。

整個(gè)模型被劃分為步長1 cm的50×50×50空間網(wǎng)格集合，采用16階勒讓德-切比雪夫求積組，置零修正的菱形差分格式。分別使用標(biāo)準(zhǔn)SN方法、首次碰撞源方法和多次碰撞源方法計(jì)算通量密度分布，其中多次碰撞源方法的計(jì)算區(qū)域?yàn)樵磪^(qū)和區(qū)域1，碰撞次數(shù)為2次。3種方法的計(jì)算時(shí)間分別為26、25和26 s，計(jì)算結(jié)果如圖4～6所示，首次碰撞源方法在一定程度上緩解了離散誤差，但射線效應(yīng)依然十分明顯，而多次碰撞源方法減弱二次射線效應(yīng)的效果更加顯著。

2.2 Kobayashi基準(zhǔn)題模型Ⅲ

Kobayashi基準(zhǔn)題是經(jīng)濟(jì)合作和發(fā)展組織核能機(jī)構(gòu)(OECD/NEA)輻射輸運(yùn)專家組提出的用于驗(yàn)證輸運(yùn)程序屏蔽計(jì)算能力的基準(zhǔn)題。選取Kobayashi基準(zhǔn)題模型Ⅲ，用標(biāo)準(zhǔn)SN方法、首次碰撞源方法和多次碰撞源方法分別進(jìn)行計(jì)算，分析多次碰撞源方法的計(jì)算精度、計(jì)算效率及離散誤差的消除效果。

a——整體模型；b——y=25 cm平面圖4 SN方法計(jì)算的通量密度分布Fig.4 Flux density distribution calculated by SN method

a——整體模型；b——y=25 cm平面圖5 首次碰撞源方法計(jì)算的通量密度分布Fig.5 Flux density distribution calculated by first collision source method

a——整體模型；b——y=25 cm平面圖6 多次碰撞源方法計(jì)算的通量密度分布Fig.6 Flux density distribution calculated by n’th collision source method

該模型是一個(gè)具有折線形孔道的深穿透屏蔽問題，幾何模型如圖7所示。孔道內(nèi)的角通量密度分布呈現(xiàn)強(qiáng)各向異性，且孔道與屏蔽區(qū)交界面處的通量密度梯度極大，存在較大的離散誤差。源強(qiáng)與截面信息列于表1，選取基準(zhǔn)報(bào)告中提供的沿孔道的22個(gè)坐標(biāo)點(diǎn)作為關(guān)鍵點(diǎn)輸出。

圖7 Kobayashi基準(zhǔn)題模型Ⅲ幾何示意圖Fig.7 Geometry of Kobayashi benchmark problem model Ⅲ

圖8 Kobayashi基準(zhǔn)題模型Ⅲ關(guān)鍵點(diǎn)通量密度計(jì)算值與基準(zhǔn)值之比Fig.8 Ratio of calculated value and benchmark value in Kobayashi benchmark problem model Ⅲ

模型劃分為步長2 cm的30×50×30空間網(wǎng)格集合，選取置零修正的菱形差分格式，迭代收斂準(zhǔn)則為1×10-4。以基準(zhǔn)報(bào)告中提供的蒙特卡羅程序計(jì)算結(jié)果作為參考值，分別采用SN方法、首次碰撞源方法和多次碰撞源方法的計(jì)算值與參考值進(jìn)行比較，其中多次碰撞源方法的計(jì)算區(qū)域?yàn)樵磪^(qū)和屏蔽區(qū)，碰撞次數(shù)為2次，沿y軸關(guān)鍵點(diǎn)的相對誤差變化曲線如圖8所示。數(shù)值結(jié)果表明，孔道的存在導(dǎo)致角通量密度的各向異性程度隨距離增加而增強(qiáng)。當(dāng)采用16階勒讓德-切比雪夫求積組計(jì)算時(shí)，孔道內(nèi)的中子通量密度非物理震蕩強(qiáng)烈；當(dāng)求積組的階數(shù)提高，采用100階PNTN求積組進(jìn)行計(jì)算，通量密度的相對誤差隨距離的增加仍呈現(xiàn)逐漸上升的趨勢。采用首次碰撞源方法計(jì)算時(shí)，同樣使用16階PNTN求積組計(jì)算得到的相對誤差整體小于10%，但離散誤差的影響并未完全消除。另一方面，由于射線追蹤的性質(zhì)，首次碰撞源方法低估了源區(qū)附近的未碰撞通量密度。采用多次碰撞源方法計(jì)算時(shí)，使用16階PNTN求積組計(jì)算結(jié)果的相對誤差整體小于5%，相對誤差的均方根和標(biāo)準(zhǔn)差均有所降低，多次碰撞源方法也存在低估源區(qū)附近通量密度的問題。

表3列出了16階PNTN求積組、首次碰撞源方法和多次碰撞源方法的計(jì)算值與參考值的相對誤差和計(jì)算時(shí)間。數(shù)值結(jié)果表明，多次碰撞源方法的相對誤差均方根小于3%，相對誤差的標(biāo)準(zhǔn)差為1.025 2×10-2，消除離散誤差的能力優(yōu)于另外兩種計(jì)算方法，但計(jì)算效率方面，相同情況下多次碰撞源方法的計(jì)算耗時(shí)較長，原因主要是多次碰撞源方法計(jì)算多次碰撞通量密度時(shí)，需對所選區(qū)域的每個(gè)網(wǎng)格進(jìn)行射線追蹤，因此多次碰撞源方法的計(jì)算量與所選計(jì)算區(qū)域的網(wǎng)格數(shù)密切相關(guān)。

表3 Kobayashi基準(zhǔn)題模型Ⅲ計(jì)數(shù)點(diǎn)位中子通量密度Table 3 Neutron flux density at key point for Kobayashi benchmark problem model Ⅲ

注：括號內(nèi)為相對誤差

2.3 自設(shè)多群固定源問題

為了分析多次碰撞源方法對多群屏蔽問題的計(jì)算精度，自設(shè)均勻介質(zhì)固定源問題。該自設(shè)問題幾何模型如圖9所示，尺寸為10 cm×10 cm×10 cm，邊界條件均為真空邊界，材料為304不銹鋼[15]；源區(qū)位于模型中心為邊長1 cm的立方體，中子源強(qiáng)為1 cm-3·s-1，能量為14 MeV。計(jì)算采用30群的MATXS-10截面庫，各向異性散射P3階展開，網(wǎng)格尺寸為0.25 cm，100階PNTN求積組的計(jì)算結(jié)果作為參考解。

圖9 固定源問題幾何模型Fig.9 Geometric model of fixed-source problem

圖10 探測器中子能譜Fig.10 Neutron spectrum of detector

圖10為標(biāo)準(zhǔn)SN方法、首次碰撞源方法和多次碰撞源方法在探測器處的能譜結(jié)果，其中多次碰撞源方法的計(jì)算區(qū)域?yàn)檎麄€(gè)模型區(qū)域，碰撞次數(shù)為2次，且僅在第2群(13.5～15 MeV)使用多次碰撞源方法。在能量較低的能量區(qū)間內(nèi)各種方法的計(jì)算結(jié)果與參考解相近，但在0.1 eV～14 MeV能量區(qū)間內(nèi)僅有多次碰撞源方法與參考解吻合較好。這是由于高能部分的群內(nèi)散射截面較小同時(shí)各向異性散射較強(qiáng)，導(dǎo)致能量最高能群的角通量密度分布各向異性程度最強(qiáng)，第2群的角度離散誤差最大，射線效應(yīng)最嚴(yán)重。而隨著粒子的慢化，較低能量能群的角通量密度各向異性程度減弱。通過比對計(jì)算，驗(yàn)證了多次碰撞源方法可有效減弱離散誤差，具有較高的可靠性。

表4列出不同條件下的計(jì)算時(shí)間，多次碰撞源方法的計(jì)算時(shí)間與16階PNTN求積組和首次碰撞源方法的計(jì)算時(shí)間相比耗時(shí)較長。綜合計(jì)算精度考慮，多次碰撞源方法相較于100階PNTN求積組，兩者精度相近，計(jì)算時(shí)間的加速比為3.49。

表4 不同條件下的計(jì)算時(shí)間Table 4 Calculation time of different conditions

3 結(jié)論

針對屏蔽計(jì)算中強(qiáng)非均勻效應(yīng)造成離散誤差過大的問題，研究了基于射線追蹤技術(shù)的多次碰撞源方法。多次碰撞源方法通過射線追蹤計(jì)算指定區(qū)域網(wǎng)格內(nèi)粒子碰撞多次的通量密度，將孤立源等效為模型內(nèi)的分布源，在求積組階數(shù)較低時(shí)仍可得到準(zhǔn)確的結(jié)果，有效地消除了離散誤差對屏蔽計(jì)算精度的影響。數(shù)值計(jì)算結(jié)果表明，多次碰撞源方法相比首次碰撞源方法具有消除次級射線效應(yīng)的能力，且能緩解因離散誤差導(dǎo)致的射線效應(yīng)和數(shù)值擴(kuò)散現(xiàn)象。