陳 瑩，梁 亮，張 乾，趙 強

(哈爾濱工程大學 核安全與仿真技術國防重點學科實驗室，黑龍江 哈爾濱 150001)

隨著反應堆物理計算精度要求的提高，基于特征線方法(MOC)的中子輸運方法得到迅速發(fā)展。該方法是對中子輸運方程的微分形式沿中子飛行軌跡進行積分的一種求解方法。理論上，該方法可求解任意復雜的幾何輸運問題[1]。近年來，國內(nèi)外的一些組件計算程序，如CASMO[2]、DRAGON[3]、WIMS[4]等，一步法堆芯計算程序，如CRX-3D[5]、DeCART[6]、nTRACER[7]、MPACT[8]、NECP-X[9]等，均使用了特征線方法作為輸運求解模塊。除能量、時間和空間變量離散外，在特征線方法中，平源近似是重要的基本假設，即同一網(wǎng)格內(nèi)源項相同。平源近似下，精確的計算結果需要較小的網(wǎng)格劃分和較小的特征線線寬。因此平源近似需要的網(wǎng)格數(shù)量較大，導致幾何信息和相關物理量的存儲空間較大，影響了總的計算效率。

為改進平源近似帶來的問題，Santandrea等[10]嘗試將源項采用線性表示。線性源是將中子輸運方程的右端項由常量替換為線性表達式。與平源近似相比，線性源可在較大網(wǎng)格劃分時得到良好的計算結果。

Petkov等[11]在MARIKO程序中使用線性源特征線方法。Santandrea等[12]提出正定的表面特征線格式及改進的非線性表面特征線格式的線性源方法。由于該方法的非線性導致很多加速方法無法使用，Le Tellier等[13]提出一種簡化線性特征線格式線性源方法，雖然這種計算方式在計算過程中使用了線性函數(shù)表示源項，但它未將通量表示為線性函數(shù)。湯春桃[14]在PEACH程序中采用坐標投影線性源方法計算源項斜率，并提出將源項在x、y方向投影的斜率計算公式，該方法的通量同樣是平通量假設，因此對計算精度的提高有限。Chai等[15]提出的線性源方法列出每條特征線起始點、中點和終點3點的平衡方程，得到表示源和通量的斜率系數(shù)的計算方式，并開發(fā)了三維線性源特征線方法程序TCM_L，該方法給出的是三維幾何模型的測試結果。Rodolfo等[16]提出一種線性源方法并應用在CASMO5中，該方法在使用線性函數(shù)表示源項的同時對通量也使用線性假設，因而理論上其對計算使用的網(wǎng)格數(shù)量不敏感，即使在大網(wǎng)格條件下也可得到精確的計算結果，但該方法在計算帶氣隙問題時需引入近似處理。朱成林[17]基于坐標投影線性源方法使用Dragon程序?qū)崿F(xiàn)了坐標投影線性源方法在三維上的應用。鄭勇[18]基于簡化線性特征線格式線性源方法開發(fā)了自適應性線性源方法，將線性源應用到矩陣特征線方法中。朱成林和鄭勇是基于上述幾種線性源方法挑選出其中一種將其應用到自己使用的程序中以達到完善程序完整性的目的，并不是方法研究。

本文提出一種采用最小二乘的線性源特征線方法，通過壓水堆柵元及C5G7基準題的計算，與坐標投影線性源特征線方法和CASMO5中運用的線性源特征線方法進行對比分析。

1 理論推導

1.1 線性源特征線方法

在笛卡爾坐標系中，線性源特征線輸運方程如式(1)所示：

(1)

(2)

式中：τ為光學長度，τ=Σtsm,i,k/sinθn；Γ1(τ)=1-e-τ；Γ3(τ)=2(τ-(1-e-τ))-τ(1-e-τ)。

(3)

式中：ωn和ωm分別為中子方向Ωmn的極角權重和方位角權重；δAm為特征線密度；Vi為平源區(qū)i的體積。

1.2 最小二乘線性源特征線方法

本文提出了一種基于最小二乘的線性源特征線方法。首先假設源項分布為：

Q(x,y)=Qi+Qi,xx+Qi,yy

(4)

(5)

(6)

式中：Qi，x與Qi，y分別為平源區(qū)i內(nèi)x、y方向的源項分布系數(shù)；(xc，yc)為特征線段的中點坐標；φm為方位角角度。

在各向同性條件下，根據(jù)中子源強的表達式可得到Qi，x與Qi，y：

(7)

(8)

假設中子角通量密度Ψ(x，y)為：

Ψ(x,y)=Ψ0+Ψxx+Ψyy

(9)

式中：Ψ0為某一網(wǎng)格內(nèi)平均中子角通量密度；Ψx和Ψy分別為x和y方向的角通量密度分布系數(shù)。

將每條特征線的起點、終點及其平均通量代入式(9)，每條特征線段可得到3個方程：

(10)

根據(jù)上述公式可得到3倍特征線段數(shù)量的方程組：

(11)

式中：(xin，yin)為某條特征線段的起點；(xout，yout)為某條特征線段的終點；Ψin、Ψc和Ψout分別為某條特征線段入射、平均和出射中子角通量密度；n為特征線段數(shù)量。

使用最小二乘法求解式(11)，即式(11)兩端分別乘以AT得到正定方程組：

(12)

(13)

(14)

給定源項初值計算得到平源區(qū)內(nèi)的標通量及標通量分布系數(shù)，再根據(jù)平源區(qū)內(nèi)的標通量及標通量分布系數(shù)得到平源區(qū)內(nèi)每條特征線段上的源強和斜率。

2 數(shù)值驗證

通過壓水堆柵元、組件問題及C5G7-MOX基準題驗證最小二乘線性源特征線方法的計算精度，作為對比，同時給出平源近似特征線方法、坐標投影線性源特征線方法、CASMO5中線性源特征線方法的計算精度。本文給出的計算時間是在CPU為Intel(R) Core(TM) i7-8550U CPU@1.80 GHz、8G內(nèi)存微機上運行結果，中子通量和keff的收斂判據(jù)分別為1.0×10-6、1.0×10-7。

2.1 燃料柵元算例

以單柵元69群算例進行驗證，柵元幾何模型包括4個區(qū)域，由內(nèi)至外分別為燃料區(qū)、氣隙、包殼及慢化劑區(qū)，幾何尺寸如圖1所示。多群截面通過蒙特卡羅程序統(tǒng)計得到。柵元的邊界條件均為反射邊界。柵元網(wǎng)格劃分采用如圖2所示的細網(wǎng)格和粗網(wǎng)格兩種。

圖3示出有效增殖因數(shù)的偏差Δkeff隨網(wǎng)格數(shù)量的變化，計算條件為：特征線密度為0.01 cm，0°～360°內(nèi)方位角數(shù)量為80個，Y-T求積組[19]極角數(shù)量為3個。選取48個網(wǎng)格下特征線線寬為0.001 cm、方位角數(shù)量為128的平源近似結果為基準解。由圖3可看出，計算時網(wǎng)格數(shù)量越多計算得到的結果越精確，且隨網(wǎng)格數(shù)量的逐漸增多，坐標投影線性源特征線方法和最小二乘線性源特征線方法的網(wǎng)格敏感性低于平源近似模型的特征線方法。最小二乘線性源特征線方法的Δkeff穩(wěn)定在10～15 pcm范圍內(nèi)，比平源近似特征線方法更加準確。

圖1 燃料柵元幾何尺寸示意圖Fig.1 Schematic of fuel cell geometry

表1列出燃料柵元算例的計算結果。由表1可見，粗網(wǎng)格(8個網(wǎng)格)條件下，線性源特征線方法與平源近似特征線方法相比計算精度更高。線性源方法中，最小二乘線性源特征線方法的計算精度高于坐標投影線性源特征線方法。計算時間方面，在精度相當?shù)募毦W(wǎng)格下，最小二乘線性源特征線方法的計算時間比平源近似特征線方法的更少。線性源特征線方法中，CASMO5在計算氣隙時無法收斂。

a——細網(wǎng)格；b——粗網(wǎng)格圖2 燃料柵元網(wǎng)格劃分示意圖Fig.2 Schematic of fuel cell grid division

圖3 燃料柵元算例Δkeff隨網(wǎng)格數(shù)量的變化Fig.3 Δkeff change with mesh number in fuel cell problem

表1 燃料柵元算例的keff和計算時間Table 1 keff and calculation time of fuel cell problem

2.2 17×17燃料組件算例

以C5G7基準題[20]中燃料組件為計算對象，燃料組件是由導向管、UO2燃料柵元和裂變室組成，方形排列17×17元件的組件。組件的幾何尺寸為21.42 cm×21.42 cm，組件中柵元布置方式如圖4所示。柵元的具體幾何信息如圖5所示。網(wǎng)格劃分采用如圖6所示的細網(wǎng)格和粗網(wǎng)格兩種。圖7示出隨網(wǎng)格數(shù)量逐漸增加不同方法進行輸運計算得到的Δkeff，計算條件為：特征線密度為0.02 cm，0°～360°內(nèi)方位角數(shù)量為96個，Y-T求積組[19]極角數(shù)量為3個。由圖7可見，在8個網(wǎng)格時線性源特征線方法計算得到的精度達到了平源近似特征線方法使用40個網(wǎng)格數(shù)量計算精度。由于算例中不存在氣隙，CASMO5中的線性源特征線方法計算過程中可收斂。隨著網(wǎng)格數(shù)量的增加，4種方法計算得到的結果更加精確，且線性源特征線方法的網(wǎng)格敏感性低于平源近似特征線方法。

圖4 柵元布置方式示意圖Fig.4 Schematic of cell layout

圖5 柵元幾何尺寸示意圖Fig.5 Schematic of cell geometry

表2列出燃料組件算例計算結果。由表2可見，在網(wǎng)格數(shù)量較少時，線性源特征線方法與平源近似特征線方法相比計算精度更高，且線性源特征線方法的計算精度與細網(wǎng)格下平源近似的計算精度相當。線性源特征線方法中，最小二乘與CASMO5的計算精度相當，且高于坐標投影線性源特征線方法。計算時間方面，最小二乘線性源特征線方法與細網(wǎng)格下的平源近似特征線方法相比，計算時間要更少。該算例驗證了組件規(guī)模問題下最小二乘線性源特征線方法準確性較好。

a——細網(wǎng)格；b——粗網(wǎng)格圖6 柵元網(wǎng)格劃分示意圖Fig.6 Schematic of cell grid division

圖7 燃料組件算例Δkeff隨網(wǎng)格數(shù)量的變化Fig.7 Δkeff change with mesh number in fuel assembly problem

表2 燃料組件算例的keff和計算時間Table 2 keff and calculation time of fuel assembly problem

2.3 C5G7-MOX基準題

C5G7-MOX基準題[20]是由二氧化鈾燃料組件和MOX燃料組件混合裝載，共計16盒燃料組件，呈1/8對稱。1/4堆芯布置如圖8所示。每個組件的幾何尺寸為21.42 cm×21.42 cm，反射層寬度為21.42 cm。二氧化鈾組件和MOX組件均為17×17的排列方式。二氧化鈾組件內(nèi)燃料棒富集度相同，MOX組件內(nèi)有3種含量不同的MOX燃料柵元，如圖9所示。柵元幾何尺寸與圖5一致，網(wǎng)格劃分方式與圖6一致。C5G7-MOX基準題keff和計算時間的計算結果列于表3，功率計算結果列于表4。

圖8 1/4堆芯布置示意圖Fig.8 Schematic of 1/4 core layout

圖9 C5G7-MOX基準題組件幾何布置示意圖Fig.9 Schematic of component geometry layout for C5G7-MOX benchmark

表3 C5G7-MOX基準題的keff和計算時間Table 3 keff and calculation time for C5G7-MOX benchmark

表4 C5G7-MOX基準題功率計算結果Table 4 Power calculation result of C5G7-MOX benchmark

由表3、4可見，用平源近似特征線方法進行計算，粗網(wǎng)格條件下的計算精度無法達到細網(wǎng)格條件下的結果。使用粗網(wǎng)格進行計算，線性源特征線方法比平源近似特征線方法計算精度更高，且線性源特征線方法的計算精度與40個網(wǎng)格下平源近似特征線方法的計算精度相當。線性源特征線方法中，最小二乘線性源特征線方法較坐標投影線性源特征線方法計算精度更高，且與CASMO5中的線性源特征線方法的精度相同。計算時間方面，兩種線性源特征線方法比40個網(wǎng)格下的平源近似特征線方法的更少。功率水平方面，粗網(wǎng)格下線性源特征線方法的計算結果高于平源近似特征線方法的計算結果，CASMO5的功率計算結果更準確，最小二乘線性源特征線方法與坐標投影線性源特征線方法的功率結果精度一致。

上述數(shù)值結果驗證了最小二乘線性源特征線方法的網(wǎng)格敏感性低且準確性較好。

3 小結

本文提出一種最小二乘線性源特征線方法，通過計算壓水堆柵元及C5G7基準題等算例，得到以下結論。

1) 最小二乘線性源特征線方法能在網(wǎng)格數(shù)量較少的計算條件下得到較精確的計算結果，且計算精度可達到細網(wǎng)格下平源近似特征線方法的計算精度。

2) 在網(wǎng)格數(shù)量較少且輸入條件相同時，線性源特征線方法比平源近似特征線方法更為準確。線性源特征線方法中，最小二乘線性源特征線方法比坐標投影線性源特征線方法更為準確，且與CASMO5中的線性源特征線方法的計算精度相當。

3) 最小二乘線性源特征線方法與相同精度下的平源近似特征線方法相比，計算時間更少。最小二乘線性源特征線方法與坐標投影線性源特征線方法相比，計算精度更高且網(wǎng)格敏感性更低；與CASMO5中的線性源特征線方法相比，最小二乘線性源特征線方法可直接處理帶氣隙的小網(wǎng)格問題，無需引入其他近似處理，適應性更好。