二維中子輸運問題的離散縱標短特征線方法研究

劉 聰，張 斌，張 亮，包博宇，趙好強，陳義學

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

離散縱標法(SN)是粒子輸運數值計算常用的確定論方法之一[1]，廣泛應用于堆芯物理分析和核設施屏蔽計算。空間變量離散是SN方程數值求解的重要組成部分，常見的離散方法包括有限差分法、有限元法等，不同離散方法各有優劣。空間離散方法的正定性、截斷誤差、網格步長敏感性等影響著輸運計算的可靠性。有限差分方法具有成熟的數學基礎且是各類SN程序普遍應用的空間離散方式，菱形差分(DD)[2]是最具代表性的低階差分方法，DD在一維平板幾何中具有二階收斂精度，但有負通量問題且在粗網格下計算精度較低。由DD發展出多種非線性低階差分方法，包括負通量置零修正菱形差分(DZ)、帶權重差分(WD)[3]、定向θ權重差分(DTW)[4]、指數定向差分(EDW)[5]和指數迭代差分(EDI)[6]等，這些差分方法通過引入修正或非線性假設保證角通量密度非負，在一定程度上改善了DD非物理振蕩和粗網精度低的問題。有限元法[7-8]利用基函數在網格內將角通量密度進行展開，求解離散后的弱形式輸運方程，具有較高計算精度和較好數值擴散特性，但多維情況下隨著展開階數升高，負通量抑制技術[9]并不能嚴格保證全局角通量密度非負。Lathrop[10]提出SN方程的常數短特征線離散方法，假設網格內源強分布構造半解析的特征線解，通過輸運方程的空間矩守恒推導網格角通量密度分布。Dehart等[11]將其推廣至二維任意形狀網格，開發了NEWT輸運計算程序。Larsen和Alcouffe[12]以及Mathews等[13]分別基于結構和非結構網格研究了線性短特征線方法，后者對負通量不進行修正。按照對網格源強分布處理方式的不同，自1990年至今陸續有多種基于短特征線離散的數值方法被提出。

相比差分方法僅考慮中子角通量密度沿單個空間變量的變化，短特征線方法基于半解析的中子衰減解，根據離散角度考慮網格出射邊界與入射的貢獻關系，減弱了因空間離散誤差造成的非物理振蕩。常數短特征線格式天然非負但計算精度不足，本文重點研究SN方程的線性短特征線空間離散方法，改善粗網精度、提高計算效率，通過對角通量密度分布的線性斜率進行旋轉修正，嚴格保證源強和通量密度非負，避免非物理結果。

1 短特征線空間離散方法

1.1 SN短特征線方法與長特征線方法的比較

圖1 SN短特征線和MOC方法的掃描示意圖Fig.1 Sketch of mesh sweep for SN short characteristic and MOC methods

短特征線方法是SN方程的空間變量離散方法，區別于組件、堆芯分析中常用的特征線方法(MOC方法，也稱為長特征線方法)。SN短特征線和MOC方法均是基于微分-積分形式中子輸運方程的確定論求解方法，方程中連續角度變量直接離散為1組有限離散方向，兩者最主要的理論差異體現在空間變量離散和網格掃描求解上。SN短特征線方法求解基于中子守恒的網格平衡方程，根據空間離散格式的理論假設，由網格入射角通量密度、網格源強求解得到網格平均和出射角通量密度，按離散方向進行角度-網格輸運掃描。短特征線離散根據特征線解與SN方程的空間矩守恒，推導網格入射與出射角通量密度的數值關系。MOC方法由平源或線性源近似對輸運方程進行積分可得到特征線解。每個離散方向下劃分1組平行射線穿過計算系統，幾何運算求出射線與網格交點、線段長度，利用射線追蹤技術進行掃描計算，求解各條射線在網格內入射、出射和該段平均角通量密度。SN短特征線和MOC方法按圖1所示路徑對計算模型進行掃描。短特征線方法基于SN網格平衡方程，針對局部網格，而MOC方法是沿穿過整個計算系統的特征線進行求解，針對網格內各段子射線。此外，SN短特征線和MOC方法在網格平均通量密度計算、邊界處理等方面存在差異。由于計算問題的不同，對角度求積組、各向異性散射展開階數的選擇也不同。

1.2 常數短特征線空間離散方法

穩態中子輸運方程中的角度變量采用SN方法直接離散，能量變量分群處理，二維xy幾何下的SN方程如下。

Σt,gψm,g(x,y)=Qm,g(x,y)

(1)

式中：ψ為中子角通量密度；Σt為宏觀總截面；x和y為空間變量；Q為源強，包含散射源、裂變源及外源；μ和η為角度向量與坐標軸夾角余弦；m和g分別為離散方向編號和能群編號，后文為簡潔表示，省略離散方向和能群編號下標。

根據離散方向導數關系將式(1)改寫為射線方程(式(2))，該方程具有式(3)形式的解。

(2)

(3)

常數短特征線方法也稱步特征線(SC)方法，假設網格源強和邊界入射角通量密度為常數分布。源強由初值或前次迭代結果給出，入射角通量密度由邊界條件或上風向網格出射值給出。建立計算網格的局部坐標系，二維幾何下角特征線與網格的相交情況如圖2所示，角特征線上、下分別由左邊界和下邊界入射貢獻。

a——交于上邊界；b——交于右邊界圖2 角特征線與網格關系圖Fig.2 Relationship between angular characteristic line and mesh

(4)

通過特征線解與目標函數的矩守恒可推導SC方法的邊界出射值與網格平均值，其中目標函數為常數分布，僅需零階矩守恒(式(5))。

(5)

(6)

(7)

(8)

計算網格的出射值作為下風向網格入射值，完成角度-網格掃描，更新網格通量矩和散射源，采用源迭代方法進行求解。

1.3 線性短特征線空間離散方法

SC方法具有正定性，在源強和入射值非負的條件下可保證角通量密度計算值非負，但由于基于常數假設，在粗網格下計算精度較低。線性短特征線(LC)方法基于線性假設改善了粗網格的計算精度，降低了輸運計算對網格步長的敏感性。計算網格局部坐標系下假設網格源強、入射角通量密度為線性分布，式(2)的SN方程改寫為式(9)形式。

0

(9)

(10)

Ψabove(x,y)=

Ψbelow(x,y)=

需要基于特征線解構造線性分布的出射ψi,out、ψj,out和網格ψcell角通量密度分布函數。

(11)

式中，ψx和ψy為角通量密度在網格內分布函數的線性斜率。

(12)

圖3 網格角通量密度求解示意圖Fig.3 Sketch for solving mesh angular flux

基于零階矩守恒可求解式(11)和(12)中邊界角通量密度平均值。

(13)

(14)

(15)

(16)

目標函數為線性分布，邊界線性斜率按式(17)根據差商近似求解。

(17)

對于相交邊界分別計算各段斜率。

(18)

根據一階矩守恒，將式(12)代入構造邊界的整體線性分布函數。

(19)

最終得到相交邊界的角通量密度分布函數的線性斜率。

(20)

(21)

(22)

(23)

按標準球諧函數展開法計算散射源，更新網格源強進行迭代計算。LC方法不是天然非負的，當網格步長過大時角通量密度分布可能在邊界或網格內出現負值，由于負通量是非物理的同時可能導致網格負散射源，因此進行旋轉斜率修正強制邊界或網格線性分布最小值為零，該修正方法破壞了網格一階矩守恒關系，但避免了數值求解的非物理結果，減弱了負通量對迭代的擾動。短特征線離散格式SC和LC已嵌入多維SN輸運計算程序ARES[14]中的二維求解模塊DONTRAN2D。

2 數值驗證

2.1 Gelbard固定源問題

該測試例題是雙群均勻介質固定源問題，取自美國阿貢國家實驗室基準題手冊ANL-7416[15]，四周為反射邊界，文獻給出P19階球諧函數解作為參考基準。模型幾何如圖4所示，材料截面和源強信息列于表1，表中，ν為每次裂變釋放的中子數，Σf為裂變截面，Σs為散射截面。

圖4 Gelbard固定源問題幾何示意圖Fig.4 Geometric model of Gelbard fixed-source problem

gνΣf,g/cm-1Σs,g-1/cm-1Σs,g-2/cm-1Σt,g/cm-1源強/(cm-3·s-1)10.06.947×10-32.343 4×10-29.210 4×10-26.546 0×10-320.00.04.850×10-31.008 77×10-11.770 1×10-2

網格尺寸設為1 cm，角度求積組選取S16階EON求積組[16]，通量密度收斂準則為5×10-5。沿y=80和y=139的計算結果示于圖5，SC和LC的計算結果與參考解吻合良好。強吸收介質引起SN方法的射線效應，造成沿y=139處的通量密度結果一定程度的振蕩。

圖5 各能群中子通量密度計算結果Fig.5 Neutron flux results for each group

2.2 自設固定源問題

圖6 單群固定源問題幾何模型Fig.6 Geometric model of one-group fixed-source problem

為分析各空間離散方法對屏蔽問題的計算精度，自設單群均勻介質固定源問題。模型如圖6所示，左邊界和下邊界為反射邊界，其余為真空邊界，材料總截面為0.5 cm-1，散射比為0.5，該模型長度為15倍中子自由程。使用多群蒙特卡羅程序MCMG[17]對該問題進行模擬，沿下邊界設置6個計數器，每個邊長為5 cm，統計區域平均中子通量密度作為參考解。

圖7 單群固定源問題的各區域平均中子通量密度和相對偏差Fig.7 Regional average neutron flux and relative deviation for one-group fixed-source problem

圖7示出了在網格步長為0.2 cm條件下SC和LC的計算結果以及MCMG各計數器的統計誤差。可看出LC結果與MCMG基準的相對偏差小于1.4%，SC計算精度不及LC方法。

網格步長分別設為0.1、0.2、0.5、1、2.5和5 cm，對短特征線、差分空間離散方法進行測試，由各區域相對偏差按式(24)計算離散方法的均方根(RMS)偏差。

(24)

式中，φn,calc和φn,ref分別為第n個區域的平均中子通量密度計算值和基準值。

圖8示出了各方法均方根偏差關于網格光學厚度(ΣtΔl，Σt為材料總截面，Δl為網格邊長)和計算時間的變化情況。可見，LC和EDW對于該問題的計算效率最高，在相近時間內可獲得更高的計算精度，對網格步長敏感性較低，網格步長為2.5倍自由程時的RMS偏差約為10%。但EDW在網格步長為0.05和0.1倍自由程時計算未收斂，由于該方法假設通量密度在網格內為指數函數分布，網格劃分過細時網格內的通量密度衰減很小，導致指數因子擬合出現問題。SC方法基于常數假設，計算結果對于網格步長十分敏感，粗網下精度不足。

圖8 各空間離散方法的RMS偏差分布Fig.8 Distribution of RMS deviation for each spatial discretization method

2.3 Stepanek沸水堆柵元臨界基準題

Stepanek等[18]提出了輕水堆基準題，用于比較不同確定論方法的計算精度與效率，本文選取其中沸水堆(BWR)柵元問題進行計算分析。該模型為雙群臨界問題，內部均勻化燃料棒由輕水包圍，四周為反射邊界，幾何模型如圖9所示，材料截面列于表2，裂變譜快群為1.0、熱群為0.0。

圖9 沸水堆柵元幾何示意圖Fig.9 Geometric model of BWR cell problem

材料gνΣf,g/cm-1Σs,g-1/cm-1Σs,g-2/cm-1Σt,g/cm-1116.203×10-31.780×10-11.002×10-21.966 47×10-121.101×10-11.089×10-35.255×10-15.961 59×10-1210.01.995×10-12.188×10-22.220 64×10-120.01.558×10-38.783×10-18.878 74×10-1

網格步長約0.2 cm，特征值收斂準則設為10-6，參考解來自MCNP的蒙特卡羅計算結果[19]。特征值Kinf和區域平均通量密度的計算結果列于表3，同時給出了ARES的有限差分離散的計算結果，區域平均通量密度按燃料區快群結果進行了歸一。該網格劃分下幾種空間離散格式的特征值結果與參考值的偏差均小于100 pcm。其中LC結果與參考值最接近，僅偏大16 pcm，其次是DZ和WD格式較準確。

改變網格步長進行網格敏感性的測試，隨網格數量和計算時間改變的特征值結果變化示于圖10。可見LC在各網格劃分下的結果均優于其他離散格式，6×6網格劃分時偏差也僅為58 pcm。幾種差分離散方法對于該問題的計算效率差異較小，DZ與WD基本一致且略優于EDW和DTW。SC隨著網格步長的增大計算結果出現較大偏差，計算精度與網格步長存在較強關聯。

表3 沸水堆柵元問題的計算結果Table 3 Result of BWR cell problem

圖10 沸水堆柵元問題特征值Fig.10 Eigenvalue result of BWR cell problem

3 結論

本文采用離散縱標短特征線方法求解二維中子輸運方程，研究了基于常數和線性分布的短特征線空間離散方法，并將該方法應用于多維SN輸運計算程序ARES中的二維求解模塊DONTRAN2D中。固定源和臨界測試例題的結果表明，線性短特征線方法對網格步長敏感性較低，相比常數短特征線和低階差分離散具有更高的粗網精度。未來將擴展至三維幾何并在實際復雜工程問題上進行應用。