線性源近似的中子輸運方程特征線解法

湯春桃

（上海核工程研究設計院，上海 200233）

中子輸運方程的特征線解法是從中子輸運方程的一階微分形式出發，沿特征線進行積分求解的方法，理論上該方法適用于任意復雜幾何輸運問題的求解。同時，該方法在求解輸運方程的過程中無需保存中子角通量和大規模耦合系數矩陣，所以十分節省計算機內存。特征線法早在20世紀50年代即被提出，但直至20世紀70年代才被運用于簡單幾何的反應堆物理計算中。最近十多年由于計算機硬件水平飛速發展以及工程計算中對復雜幾何處理能力的需求，該方法被國際上眾多機構廣泛研究。近期開發或維護的組件計算程序均采用了這種算法，如 CASMO－4E［1］、HELIOS－2［2］、WIMS［3］、DRAGON［4］等。

但現有的大部分特征線法輸運程序均是基于平源近似的步特征線法（SC）模型開發的，平源近似是特征線法理論模型中除角度變量直接離散外又一基本假定。為取得足夠的計算精度，該方法須采用足夠致密的特征線布置，同時平源近似的區域劃分也須足夠小。增加平源近似區的數目，一方面會給特征線的布置帶來更高的要求，另一方面會增加幾何以及相關物理量的存儲空間，還會極大地增加相應的跟蹤計算時間。

本工作在平源近似的步特征線法理論模型的基礎上，提出基于線性源近似的特征線法（LS）改進跟蹤計算模型。線性源近似是平源近似的拓展，通過在輸運方程右端中子源項部分引入線性項表達式實現。與平源近似的步特征線法相比，它可在采用較大網格剖分的前提下，獲得很高的計算精度。其中，中子源項表達式中線性斜率的計算是模型的核心，本文研究該線性斜率的解析表達式，完善坐標投影法的理論模型。此外，在探討這種改進的數值計算模型的過程中，提出相關負中子源分布的修正方法，使數值收斂過程穩定。

1 理論模型

與SN方法相同，特征線法也是由經角度變量直接離散后的中子輸運方程出發的。中子沿某一方向飛行滿足穩態多群中子輸運方程：

式中：θn為中子飛行方向Ωmn（m為方位角、n為極角）的極角；s為中子運行軌跡在x－y平面的投影；ψg（s）為中子角通量；Qg（s）為該處的中子源項；Σg（s）為宏觀總截面；下標g為能群。

在不考慮散射各向異性的情況下，式（1）右端源項可表示為：

式中：χg為裂變份額；keff為有效增殖因數；ν為每次裂變中子數；Σf，g′為裂變截面；Σs，g′→g為散射截面。

從式（2）可看出，Qg（s）應具備與ψ（s）類似的形狀分布。實際上，在通常的步特征線法計算模型中，在同一網格i內均未考慮Qg（s）的空間分布，而是將其看作一常數。這是特征線法理論模型中，除角度變量直接離散外又一基本假定。欲在這種平源近似的假定下獲得較高的計算精度，特征線法的空間網格離散（即平源近似區）的尺寸須劃分得足夠小。增加平源近似區的數目、減小平源近似區的尺寸，一方面會增加幾何及相關物理量的存儲，另一方面也會增加相應的處理時間。同時，還會給特征線的布置提出更高要求。

圖1a示出實際壓水堆燃料柵元常用的三區模型，欲獲得準確的計算結果，在平源近似的特征線法程序中需建立如圖1b所示的空間網格劃分。基于線性源近似的特征線法會緩解這一問題，它可在采用較大空間網格劃分的情況下，依然獲得良好的數值計算精度。

圖1 燃料柵元三區模型（a）與平源近似區網格劃分（b）Fig.1 Fuel cell model with three region（a）and mesh division for flat source approximation（b）

1.1 步特征線法

根據式（1），從中子輸運方程的一階微分形式出發，在假設網格i內宏觀截面和中子源強為常數的前提下，穿過該網格的第k條特征線的出射中子角通量可表示為：

由此，網格i內的平均標量中子通量可通過式（3）在所有網格i內的特征線及所有離散方向上積分獲得：

式中：ωm和ωn分別為對應離散方向Ωmn的輻角和極角權重；δAm為特征線間的投影間隔；Ai為網格i的面積。

為保證實際區域i的面積守恒，幾何掃描產生特征線時，投影長度sm，i，k需經式（5）修正：

1.2 線性源近似特征線法

在笛卡爾坐標（x－y）系統內，本工作探討的線性源特征線法模型的右端源分布采用如下表達形式：

為源分布的線性斜率，可表示為：

圖2示出坐標投影法示意圖。

在假設入射中子角通量、網格i的平均源強Qi和線性斜率已知的情況下，將式（1）沿特征線k積分，即可求得出射角中子通量的表達式：

圖2 坐標投影法示意圖Fig.2 Sketch map for projection method

其中：τ＝Σism，i，k／sinθn，為三維空間的光學距離；Ei函數可按下式定義：

根據式（9）求得出射角中子通量后，更新網格平均標量中子通量，進而更新散射源、裂變源等過程與平源近似的步特征線模型相同。需注意，每完成1次外迭代均需根據出射角中子通量更新線性斜率的貢獻項、：

式中：Km為經過區域i、方向為m的特征線總數目；Nθ、Nφ分別為方向離散輻角和極角的總數目。

在線性源近似模型中，中子源強分布函數非負是數值計算過程中對物理量的基本要求。式（6）右端項非負的充要條件是：

其中，（sm，i，k）max是特征線段sm，i，k在 網格i內的最大值。

在投影法求線性斜率的過程中，式（12）并不能自動滿足，這意味著若不做相應的修正，很有可能出現負中子源強，這是不允許發生的。

在迭代求解過程中，一旦式（12）得不到滿足，程序將會在保留原斜率符號的前提下自動調整它的數值，即：

在實際計算過程中，不滿足式（12）的情況絕大部分只發生在迭代初期。經基準題的數值檢驗，式（13）的修正方式可使整個收斂過程更穩定。另外，它既不影響網格內的中子平衡，又在物理量可允許的范圍內最大限度地考慮到了中子源強的線性分布，是一種非常有效的修正方法。

2 數值檢驗

根據上述理論模型，在自行研制的特征線法程序PEACH［5－6］中加入了線性源近似的計算模型，該模型的引入并不影響程序中原有的幾何處理方法以及各種加速方法的實施。為了檢驗線性源近似模型的精度與速度，采用由OECD／NEA發布的廣泛用于檢驗輸運程序求解非均勻堆芯能力的C5G7－MOX 2D基準問題［7］和自定義沸水堆小堆芯問題。

2.1 C5G7－MOX 2D基準問題

該基準題堆芯由UO2燃料組件和MOX燃料組件混合裝載，共計16盒燃料組件，呈1／8對稱。由于它具有強泄漏、組件間能譜差異大、非均勻性強等特點，目前被美、日、韓等國家的研究機構廣泛用于新一代堆芯物理分析。關于該基準題的具體幾何結構和各種材料的宏觀截面等參數參考文獻［7］。

程序PEACH在計算C5G7－MOX 2D基準問題時，將特征線法的標準控制參數設定為：1／8卦限內輻角數目為10、極角數目為3，平均相鄰特征線間的間隔為0.02cm，每燃料柵元內的網格數為48（圖3），靠近燃料棒的4列反射層柵元采用6×6等分、再外面的反射層采用4×4等分。該標準參數下的數值結果由PEACH的平源近似步特征線法模型給出。

表1列出C5G7－MOX 2D問題特征線法平源近似和線性源近似結果比較。從表1可見，PEACH已達到文獻［7］公布的國際同類軟件的計算精度。

圖3 C5G7－MOX 2D燃料柵元的2種網格劃分Fig.3 Two types of fuel cell division employed for C5G7－MOX 2Dproblem

為展示本工作提出的線性源特征線法的優越性，將控制參數放松為：1／8卦限內輻角數目為8、極角數目為3，平均相鄰特征線間的間隔為0.04cm，每燃料柵元內的網格數為8（圖3），反射層柵元均采用4×4等分。由表1可知，一旦平源近似網格劃分較大，平源近似特征線法模型的精度將受到影響，keff的偏差為29pcm、棒功率的最大相對偏差可達4.18%。圖4分別示出這兩種計算模塊得出的精細棒功率相對偏差分布。由圖4可見，平源近似特征線法的棒功率分布明顯存在傾斜。這一現象是較易理解的，靠近反射層的燃料柵元內注量梯度很大，較大網格的平源近似模型很難體現源的梯度。

但在相同控制參數的情況下，線性源近似特征線法模塊可給出與標準參數下精度相當的計算結果。另外，從計算機的耗時和存儲量兩方面來看，線性源近似特征線法模型均是標準參數下平源近似特征線法模塊的1／2，可見，線性源近似特征線法具備較明顯的優勢。

2.2 自定義沸水堆小堆芯問題

該問題的1／4堆芯由9盒燃料組件組成，堆芯外圍是寬度為15.24cm的水反射層。堆芯布置如圖5所示，燃料組件的外形尺寸為15.24cm×15.24cm，分成兩種類型，一種是新料（組件編號為1、2、3），另一種是燃耗深度為20GW·d／tU的舊料（組件編號為4、5、6）。

燃料組件內部共含有6種富集度的UO2燃料和2種類型的含GD燃料棒。燃料柵元為三區的幾何結構，分別是燃料、包殼、冷卻劑。組件中心位置為約2×2柵元尺寸的水洞，每1／4組件由5×5的燃料柵元規則排列組成。燃料組件外圍由盒壁固定，盒壁與燃料間有很窄的間隙，盒壁外面分別是寬水隙和窄水隙，其中，寬水隙是為給控制棒提供足夠的插入空間。

表1 C5G7－MOX 2D問題特征線法平源近似和線性源近似結果比較Table 1 Performance comparisons between SC and LS schemes for C5G7－MOX 2Dproblem

圖4 粗控制參數下平源近似（a）和線性源近似（b）的精細棒功率相對偏差分布Fig.4 Comparison of pin power relative deviation distributions between SC scheme（a）and LS scheme（b）with coarse parameters

各種材料的69群宏觀截面由IAEA 69群截面庫經DRAGON程序［4］計算產生。該問題的參考解（keff和精細棒功率分布）由多群蒙特卡羅程序 MCMG［8］給出。為獲得可靠的計算結果，在MCMG計算過程中每代投入的粒子數為160 000，總共計算1 200代，其中前100代不參與統計。

在計算該問題時，將特征線法的標準控制參數設定為：1／8卦限內輻角數目為8、極角數目為3，平均相鄰特征線間的間隔為0.03cm，每燃料柵元內的網格數為40，靠近燃料棒的4列反射層柵元采用5×5等分、再外面的反射層用3×3等分。該標準參數下的數值結果由PEACH的平源近似步特征線法模型給出。

圖5 自定義沸水堆小堆芯問題布置Fig.5 Configuration of self－defined BWR mini－core problem

為了展示本工作提出的線性源特征線法的優越性，將控制參數放松為：1／8卦限內輻角數目為4、極角數目為2，平均相鄰特征線間的間隔為0.07cm，每燃料柵元內的網格數為12，反射層柵元網格劃分與標準參數時相同。

表2列出程序PEACH平源近似模塊和線性源近似模塊針對該問題與參考解相比較的偏差。由表2可知，一旦平源近似網格劃分較大時，平源近似特征線法模塊的精度將受影響，針對該問題，keff的偏差為－140pcm、最大棒功率的相對偏差可達3.54%。但在相同控制參數的情況下，線性源近似特征線法模塊可給出與標準參數下精度相當的計算結果。另外，從計算機的耗時和存儲量兩方面來看，線性源近似特征線法模塊不到標準參數下平源近似特征線法模塊的1／2，再次證明線性源近似特征線法具備較明顯的優勢。

表2 沸水堆小堆芯問題特征線法平源近似和線性源近似結果比較Table 2 Performance comparisons between SC and LS schemes for BWR mini－core problem

3 結論

本工作提出一種基于坐標投影法的線性源特征線法跟蹤計算模型，提出了迭代求解過程中相關負中子源分布的修正方法，并將該模型成功加入程序PEACH中。通過C5G7－MOX 2D基準題和自定義沸水堆小堆芯問題的數值驗算結果表明，本工作提出的線性源近似特征線法模型在相同計算精度的前提下，占用更少的系統內存和運行時間，值得在中子輸運方程的特征線解法領域中推廣。

本工作是在上海交通大學趙榮安教授和張少泓副教授指導下完成的，在此對二位導師表示最誠摯的感謝。

［1］SMITH K S，RHODES J D.Full－core 2－D LWR core calculations with CASMO－4E［C／CD］∥Proceedings of PHYSOR 2002.Seoul，Korea：［s.n.］，2002.

［2］WEMPLE C A， GHEORGHIU H N M，STAMM′LER R J J，et al.Recent advances in the HELIOS－2lattice physics code［C］∥Proceedings of PHYSOR 2008.Interlaken，Switzerland：［s.n.］，2008.

［3］NETWON T，HOSKING G，HUTTON L，et al.Developments within WIMS10［C］∥Proceedings of PHYSOR 2008.Interlaken，Switzerland：［s.n.］，2008.

［4］MARLEAU G，HéBERT A，ROY R.A user’s guide for DRAGON，IGE－174［R］.［S.l.］：Ecole Polytechnique de Montréal，2007.

［5］湯春桃，張少泓.以柵元為模塊進行特征線跟蹤的中子輸運方程解法［J］.核動力工程，2009，30（4）：32－36.TANG Chuntao，ZHANG Shaohong.Study on method of characteristics based on cell modular ray tracing［J］.Nuclear Power Engineering，2009，30（4）：32－36（in Chinese）.

［6］湯春桃，張少泓.CMFD加速在特征線法輸運計算中的應用［J］.核動力工程，2009，30（5）：8－12.TANG Chuntao，ZHANG Shaohong.Application of coarse－mesh finite difference acceleration in transportation calculation by method of characteristics［J］.Nuclear Power Engineering，2009，30（5）：8－12（in Chinese）.

［7］LEWIS E E，SMITH M A，TSOULFANIDIS N，et al.Benchmark specification for deterministic 2－D／3－D MOX fuel assembly transport calculations without spatial homogenisation （C5G7－MOX），NEA／NSC／DOC（2001）4［R］.USA：OECD／NEA，2001.

［8］DENG L，XIE Z S，ZHANG J M.A 3－D multigroup P－3Monte Carlo code and its benchmarks［J］.J Nucl Sci Technol，2000，37（7）：608－614.