線性源近似的中子輸運(yùn)方程特征線解法

湯春桃

（上海核工程研究設(shè)計(jì)院，上海 200233）

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

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

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

1 理論模型

與SN方法相同，特征線法也是由經(jīng)角度變量直接離散后的中子輸運(yùn)方程出發(fā)的。中子沿某一方向飛行滿足穩(wěn)態(tài)多群中子輸運(yùn)方程：

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

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

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

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

圖1a示出實(shí)際壓水堆燃料柵元常用的三區(qū)模型，欲獲得準(zhǔn)確的計(jì)算結(jié)果，在平源近似的特征線法程序中需建立如圖1b所示的空間網(wǎng)格劃分?；诰€性源近似的特征線法會(huì)緩解這一問(wèn)題，它可在采用較大空間網(wǎng)格劃分的情況下，依然獲得良好的數(shù)值計(jì)算精度。

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

1.1 步特征線法

根據(jù)式（1），從中子輸運(yùn)方程的一階微分形式出發(fā)，在假設(shè)網(wǎng)格i內(nèi)宏觀截面和中子源強(qiáng)為常數(shù)的前提下，穿過(guò)該網(wǎng)格的第k條特征線的出射中子角通量可表示為：

由此，網(wǎng)格i內(nèi)的平均標(biāo)量中子通量可通過(guò)式（3）在所有網(wǎng)格i內(nèi)的特征線及所有離散方向上積分獲得：

式中：ωm和ωn分別為對(duì)應(yīng)離散方向Ωmn的輻角和極角權(quán)重；δAm為特征線間的投影間隔；Ai為網(wǎng)格i的面積。

為保證實(shí)際區(qū)域i的面積守恒，幾何掃描產(chǎn)生特征線時(shí)，投影長(zhǎng)度sm，i，k需經(jīng)式（5）修正：

1.2 線性源近似特征線法

在笛卡爾坐標(biāo)（x－y）系統(tǒng)內(nèi)，本工作探討的線性源特征線法模型的右端源分布采用如下表達(dá)形式：

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

圖2示出坐標(biāo)投影法示意圖。

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

圖2 坐標(biāo)投影法示意圖Fig.2 Sketch map for projection method

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

根據(jù)式（9）求得出射角中子通量后，更新網(wǎng)格平均標(biāo)量中子通量，進(jìn)而更新散射源、裂變?cè)吹冗^(guò)程與平源近似的步特征線模型相同。需注意，每完成1次外迭代均需根據(jù)出射角中子通量更新線性斜率的貢獻(xiàn)項(xiàng)、：

式中：Km為經(jīng)過(guò)區(qū)域i、方向?yàn)閙的特征線總數(shù)目；Nθ、Nφ分別為方向離散輻角和極角的總數(shù)目。

在線性源近似模型中，中子源強(qiáng)分布函數(shù)非負(fù)是數(shù)值計(jì)算過(guò)程中對(duì)物理量的基本要求。式（6）右端項(xiàng)非負(fù)的充要條件是：

其中，（sm，i，k）max是特征線段sm，i，k在 網(wǎng)格i內(nèi)的最大值。

在投影法求線性斜率的過(guò)程中，式（12）并不能自動(dòng)滿足，這意味著若不做相應(yīng)的修正，很有可能出現(xiàn)負(fù)中子源強(qiáng)，這是不允許發(fā)生的。

在迭代求解過(guò)程中，一旦式（12）得不到滿足，程序?qū)?huì)在保留原斜率符號(hào)的前提下自動(dòng)調(diào)整它的數(shù)值，即：

在實(shí)際計(jì)算過(guò)程中，不滿足式（12）的情況絕大部分只發(fā)生在迭代初期。經(jīng)基準(zhǔn)題的數(shù)值檢驗(yàn)，式（13）的修正方式可使整個(gè)收斂過(guò)程更穩(wěn)定。另外，它既不影響網(wǎng)格內(nèi)的中子平衡，又在物理量可允許的范圍內(nèi)最大限度地考慮到了中子源強(qiáng)的線性分布，是一種非常有效的修正方法。

2 數(shù)值檢驗(yàn)

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

2.1 C5G7－MOX 2D基準(zhǔn)問(wèn)題

該基準(zhǔn)題堆芯由UO2燃料組件和MOX燃料組件混合裝載，共計(jì)16盒燃料組件，呈1／8對(duì)稱。由于它具有強(qiáng)泄漏、組件間能譜差異大、非均勻性強(qiáng)等特點(diǎn)，目前被美、日、韓等國(guó)家的研究機(jī)構(gòu)廣泛用于新一代堆芯物理分析。關(guān)于該基準(zhǔn)題的具體幾何結(jié)構(gòu)和各種材料的宏觀截面等參數(shù)參考文獻(xiàn)［7］。

程序PEACH在計(jì)算C5G7－MOX 2D基準(zhǔn)問(wèn)題時(shí)，將特征線法的標(biāo)準(zhǔn)控制參數(shù)設(shè)定為：1／8卦限內(nèi)輻角數(shù)目為10、極角數(shù)目為3，平均相鄰特征線間的間隔為0.02cm，每燃料柵元內(nèi)的網(wǎng)格數(shù)為48（圖3），靠近燃料棒的4列反射層?xùn)旁捎?×6等分、再外面的反射層采用4×4等分。該標(biāo)準(zhǔn)參數(shù)下的數(shù)值結(jié)果由PEACH的平源近似步特征線法模型給出。

表1列出C5G7－MOX 2D問(wèn)題特征線法平源近似和線性源近似結(jié)果比較。從表1可見(jiàn)，PEACH已達(dá)到文獻(xiàn)［7］公布的國(guó)際同類軟件的計(jì)算精度。

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

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

但在相同控制參數(shù)的情況下，線性源近似特征線法模塊可給出與標(biāo)準(zhǔn)參數(shù)下精度相當(dāng)?shù)挠?jì)算結(jié)果。另外，從計(jì)算機(jī)的耗時(shí)和存儲(chǔ)量?jī)煞矫鎭?lái)看，線性源近似特征線法模型均是標(biāo)準(zhǔn)參數(shù)下平源近似特征線法模塊的1／2，可見(jiàn)，線性源近似特征線法具備較明顯的優(yōu)勢(shì)。

2.2 自定義沸水堆小堆芯問(wèn)題

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

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

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

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

各種材料的69群宏觀截面由IAEA 69群截面庫(kù)經(jīng)DRAGON程序［4］計(jì)算產(chǎn)生。該問(wèn)題的參考解（keff和精細(xì)棒功率分布）由多群蒙特卡羅程序 MCMG［8］給出。為獲得可靠的計(jì)算結(jié)果，在MCMG計(jì)算過(guò)程中每代投入的粒子數(shù)為160 000，總共計(jì)算1 200代，其中前100代不參與統(tǒng)計(jì)。

在計(jì)算該問(wèn)題時(shí)，將特征線法的標(biāo)準(zhǔn)控制參數(shù)設(shè)定為：1／8卦限內(nèi)輻角數(shù)目為8、極角數(shù)目為3，平均相鄰特征線間的間隔為0.03cm，每燃料柵元內(nèi)的網(wǎng)格數(shù)為40，靠近燃料棒的4列反射層?xùn)旁捎?×5等分、再外面的反射層用3×3等分。該標(biāo)準(zhǔn)參數(shù)下的數(shù)值結(jié)果由PEACH的平源近似步特征線法模型給出。

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

為了展示本工作提出的線性源特征線法的優(yōu)越性，將控制參數(shù)放松為：1／8卦限內(nèi)輻角數(shù)目為4、極角數(shù)目為2，平均相鄰特征線間的間隔為0.07cm，每燃料柵元內(nèi)的網(wǎng)格數(shù)為12，反射層?xùn)旁W(wǎng)格劃分與標(biāo)準(zhǔn)參數(shù)時(shí)相同。

表2列出程序PEACH平源近似模塊和線性源近似模塊針對(duì)該問(wèn)題與參考解相比較的偏差。由表2可知，一旦平源近似網(wǎng)格劃分較大時(shí)，平源近似特征線法模塊的精度將受影響，針對(duì)該問(wèn)題，keff的偏差為－140pcm、最大棒功率的相對(duì)偏差可達(dá)3.54%。但在相同控制參數(shù)的情況下，線性源近似特征線法模塊可給出與標(biāo)準(zhǔn)參數(shù)下精度相當(dāng)?shù)挠?jì)算結(jié)果。另外，從計(jì)算機(jī)的耗時(shí)和存儲(chǔ)量?jī)煞矫鎭?lái)看，線性源近似特征線法模塊不到標(biāo)準(zhǔn)參數(shù)下平源近似特征線法模塊的1／2，再次證明線性源近似特征線法具備較明顯的優(yōu)勢(shì)。

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

3 結(jié)論

本工作提出一種基于坐標(biāo)投影法的線性源特征線法跟蹤計(jì)算模型，提出了迭代求解過(guò)程中相關(guān)負(fù)中子源分布的修正方法，并將該模型成功加入程序PEACH中。通過(guò)C5G7－MOX 2D基準(zhǔn)題和自定義沸水堆小堆芯問(wèn)題的數(shù)值驗(yàn)算結(jié)果表明，本工作提出的線性源近似特征線法模型在相同計(jì)算精度的前提下，占用更少的系統(tǒng)內(nèi)存和運(yùn)行時(shí)間，值得在中子輸運(yùn)方程的特征線解法領(lǐng)域中推廣。

本工作是在上海交通大學(xué)趙榮安教授和張少泓副教授指導(dǎo)下完成的，在此對(duì)二位導(dǎo)師表示最誠(chéng)摯的感謝。

［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ìn)行特征線跟蹤的中子輸運(yùn)方程解法［J］.核動(dòng)力工程，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加速在特征線法輸運(yùn)計(jì)算中的應(yīng)用［J］.核動(dòng)力工程，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.