基于最小二乘的線性源特征線方法研究

陳 瑩，梁 亮，張 乾，趙 強(qiáng)

(哈爾濱工程大學(xué) 核安全與仿真技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，黑龍江 哈爾濱 150001)

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

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

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

本文提出一種采用最小二乘的線性源特征線方法，通過(guò)壓水堆柵元及C5G7基準(zhǔn)題的計(jì)算，與坐標(biāo)投影線性源特征線方法和CASMO5中運(yùn)用的線性源特征線方法進(jìn)行對(duì)比分析。

1 理論推導(dǎo)

1.1 線性源特征線方法

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

(1)

(2)

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

(3)

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

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

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

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

(4)

(5)

(6)

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

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

(7)

(8)

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

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

(9)

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

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

(10)

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

(11)

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

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

(12)

(13)

(14)

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

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

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

2.1 燃料柵元算例

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

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

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

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

a——細(xì)網(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和計(jì)算時(shí)間Table 1 keff and calculation time of fuel cell problem

2.2 17×17燃料組件算例

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

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

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

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

a——細(xì)網(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和計(jì)算時(shí)間Table 2 keff and calculation time of fuel assembly problem

2.3 C5G7-MOX基準(zhǔn)題

C5G7-MOX基準(zhǔn)題[20]是由二氧化鈾燃料組件和MOX燃料組件混合裝載，共計(jì)16盒燃料組件，呈1/8對(duì)稱。1/4堆芯布置如圖8所示。每個(gè)組件的幾何尺寸為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基準(zhǔn)題keff和計(jì)算時(shí)間的計(jì)算結(jié)果列于表3，功率計(jì)算結(jié)果列于表4。

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

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

表3 C5G7-MOX基準(zhǔn)題的keff和計(jì)算時(shí)間Table 3 keff and calculation time for C5G7-MOX benchmark

表4 C5G7-MOX基準(zhǔn)題功率計(jì)算結(jié)果Table 4 Power calculation result of C5G7-MOX benchmark

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

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

3 小結(jié)

本文提出一種最小二乘線性源特征線方法，通過(guò)計(jì)算壓水堆柵元及C5G7基準(zhǔn)題等算例，得到以下結(jié)論。

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

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

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