劉 聰，張 斌，張 亮，包博宇，趙好強(qiáng)，陳義學(xué)

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

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

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

1 短特征線空間離散方法

1.1 SN短特征線方法與長(zhǎng)特征線方法的比較

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

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

1.2 常數(shù)短特征線空間離散方法

穩(wěn)態(tài)中子輸運(yùn)方程中的角度變量采用SN方法直接離散，能量變量分群處理，二維xy幾何下的SN方程如下。

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

(1)

式中：ψ為中子角通量密度；Σt為宏觀總截面；x和y為空間變量；Q為源強(qiáng)，包含散射源、裂變?cè)醇巴庠?；μ和η為角度向量與坐標(biāo)軸夾角余弦；m和g分別為離散方向編號(hào)和能群編號(hào)，后文為簡(jiǎn)潔表示，省略離散方向和能群編號(hào)下標(biāo)。

根據(jù)離散方向?qū)?shù)關(guān)系將式(1)改寫(xiě)為射線方程(式(2))，該方程具有式(3)形式的解。

(2)

(3)

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

a——交于上邊界；b——交于右邊界圖2 角特征線與網(wǎng)格關(guān)系圖Fig.2 Relationship between angular characteristic line and mesh

(4)

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

(5)

(6)

(7)

(8)

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

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

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

0

(9)

(10)

Ψabove(x,y)=

Ψbelow(x,y)=

需要基于特征線解構(gòu)造線性分布的出射ψi,out、ψj,out和網(wǎng)格ψcell角通量密度分布函數(shù)。

(11)

式中，ψx和ψy為角通量密度在網(wǎng)格內(nèi)分布函數(shù)的線性斜率。

(12)

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

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

(13)

(14)

(15)

(16)

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

(17)

對(duì)于相交邊界分別計(jì)算各段斜率。

(18)

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

(19)

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

(20)

(21)

(22)

(23)

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

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

2.1 Gelbard固定源問(wèn)題

該測(cè)試?yán)}是雙群均勻介質(zhì)固定源問(wèn)題，取自美國(guó)阿貢國(guó)家實(shí)驗(yàn)室基準(zhǔn)題手冊(cè)ANL-7416[15]，四周為反射邊界，文獻(xiàn)給出P19階球諧函數(shù)解作為參考基準(zhǔn)。模型幾何如圖4所示，材料截面和源強(qiáng)信息列于表1，表中，ν為每次裂變釋放的中子數(shù)，Σf為裂變截面，Σs為散射截面。

圖4 Gelbard固定源問(wèn)題幾何示意圖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源強(qiáng)/(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

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

圖5 各能群中子通量密度計(jì)算結(jié)果Fig.5 Neutron flux results for each group

2.2 自設(shè)固定源問(wèn)題

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

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

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

圖7示出了在網(wǎng)格步長(zhǎng)為0.2 cm條件下SC和LC的計(jì)算結(jié)果以及MCMG各計(jì)數(shù)器的統(tǒng)計(jì)誤差?？煽闯鯨C結(jié)果與MCMG基準(zhǔn)的相對(duì)偏差小于1.4%，SC計(jì)算精度不及LC方法。

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

(24)

式中，φn,calc和φn,ref分別為第n個(gè)區(qū)域的平均中子通量密度計(jì)算值和基準(zhǔn)值。

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

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

2.3 Stepanek沸水堆柵元臨界基準(zhǔn)題

Stepanek等[18]提出了輕水堆基準(zhǔn)題，用于比較不同確定論方法的計(jì)算精度與效率，本文選取其中沸水堆(BWR)柵元問(wèn)題進(jìn)行計(jì)算分析。該模型為雙群臨界問(wèn)題，內(nèi)部均勻化燃料棒由輕水包圍，四周為反射邊界，幾何模型如圖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

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

改變網(wǎng)格步長(zhǎng)進(jìn)行網(wǎng)格敏感性的測(cè)試，隨網(wǎng)格數(shù)量和計(jì)算時(shí)間改變的特征值結(jié)果變化示于圖10。可見(jiàn)LC在各網(wǎng)格劃分下的結(jié)果均優(yōu)于其他離散格式，6×6網(wǎng)格劃分時(shí)偏差也僅為58 pcm。幾種差分離散方法對(duì)于該問(wèn)題的計(jì)算效率差異較小，DZ與WD基本一致且略優(yōu)于EDW和DTW。SC隨著網(wǎng)格步長(zhǎng)的增大計(jì)算結(jié)果出現(xiàn)較大偏差，計(jì)算精度與網(wǎng)格步長(zhǎng)存在較強(qiáng)關(guān)聯(lián)。

表3 沸水堆柵元問(wèn)題的計(jì)算結(jié)果Table 3 Result of BWR cell problem

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

3 結(jié)論

本文采用離散縱標(biāo)短特征線方法求解二維中子輸運(yùn)方程，研究了基于常數(shù)和線性分布的短特征線空間離散方法，并將該方法應(yīng)用于多維SN輸運(yùn)計(jì)算程序ARES中的二維求解模塊DONTRAN2D中。固定源和臨界測(cè)試?yán)}的結(jié)果表明，線性短特征線方法對(duì)網(wǎng)格步長(zhǎng)敏感性較低，相比常數(shù)短特征線和低階差分離散具有更高的粗網(wǎng)精度。未來(lái)將擴(kuò)展至三維幾何并在實(shí)際復(fù)雜工程問(wèn)題上進(jìn)行應(yīng)用。