基于先進(jìn)節(jié)塊法的六邊形三維中子物理程序的開發(fā)與驗(yàn)證

陸道綱，呂思宇,*，隋丹婷，郭勁松

(1.華北電力大學(xué) 核科學(xué)與工程學(xué)院，北京 102206；2.非能動核能安全技術(shù)北京重點(diǎn)實(shí)驗(yàn)室，北京 102206)

三維節(jié)塊擴(kuò)散理論作為準(zhǔn)確且有效的求解全堆芯耦合擴(kuò)散方程的方法被核工程領(lǐng)域廣泛接受，并應(yīng)用于輕水堆(LWR)與快堆(FBR)堆芯三維中子通量、功率分布的仿真計算中。由于目前世界上大部分反應(yīng)堆燃料組件為方形燃料組件，先進(jìn)的三維節(jié)塊擴(kuò)散理論往往是率先應(yīng)用于笛卡爾坐標(biāo)系下的中子擴(kuò)散方程求解。為使這些先進(jìn)節(jié)塊擴(kuò)散理論能應(yīng)用于具有六邊形燃料組件的堆型(如鈉冷快堆)中，各科研院所提出了很多方法。Wagner等[1-2]提出的應(yīng)用于DIF3D程序中的基于橫向積分技術(shù)的六邊形節(jié)塊展開法，Koyanma等[3-4]提出的通過邊界擬合坐標(biāo)轉(zhuǎn)換法將六邊形轉(zhuǎn)變?yōu)榈芽栕鴺?biāo)系的節(jié)塊展開法，和Chao等[5-7]提出的基于共形映射的六邊形節(jié)塊法等。在前期的工作中，華北電力大學(xué)提出了可應(yīng)用于六邊形組件的高階節(jié)塊展開法，并基于該方法開發(fā)了三維中子物理程序HNHEX[8-9]，該程序以六邊形組件為徑向計算節(jié)點(diǎn)，并根據(jù)需要對軸向進(jìn)行劃分，可計算三維中子通量分布，通過基準(zhǔn)題對該程序進(jìn)行了驗(yàn)證并得到了較好的仿真結(jié)果。但不可忽視的是，該方法是基于橫向積分將六邊形節(jié)塊內(nèi)部的中子擴(kuò)散方程通過3個徑向方向、1個軸向方向并通過泄漏項(xiàng)耦合的一維中子擴(kuò)散方程組。這種處理方法會在對面平均通量求二階導(dǎo)時產(chǎn)生奇異項(xiàng)。六邊形高階節(jié)塊展開法通過采用Wagner在文獻(xiàn)[1]中提出的方法處理奇異項(xiàng)，即通過在一維擴(kuò)散方程中忽略奇異項(xiàng)，在節(jié)點(diǎn)耦合方程中引入補(bǔ)償項(xiàng)來修復(fù)節(jié)點(diǎn)平衡關(guān)系。相對于笛卡爾坐標(biāo)系下的節(jié)塊展開法，這種近似損失了計算精度。而通過共形映射先將六邊形各邊映射到長方形上，再進(jìn)行橫向積分得到各方向上的橫向一維中子擴(kuò)散方程，則避免了奇異項(xiàng)的產(chǎn)生[5,10]，這種處理方法被證明可顯著減小計算誤差。

由于高階節(jié)塊展開法具有較快的計算速度，且目前反應(yīng)堆中材料不均勻性主要來自于徑向，在軸向不均勻性并不顯著。本文通過在軸向上使用高階節(jié)塊展開法，并將HNHEX程序中徑向上使用的展開法替換為基于共形映射的半解析節(jié)塊法，通過泄漏項(xiàng)進(jìn)行耦合開發(fā)SA-HNHEX(semi analytic-high order nodal expansion method for hexagonal-z)程序。對VVER-440的二維、三維基準(zhǔn)題分別進(jìn)行建模計算，初步驗(yàn)證該方法與程序的有效性和正確性。

1 求解方法

堆芯以每盒組件為1個徑向節(jié)點(diǎn)，軸向節(jié)點(diǎn)數(shù)量可根據(jù)需要自行設(shè)定，通常為10 cm至25 cm高度。多群擴(kuò)散方程可寫成如下形式：

(1)

(2)

對式(1)的求解需要建立額外的面平均泄漏和節(jié)點(diǎn)中子通量之間的關(guān)系。將式(1)分別對x-y平面與z方向進(jìn)行平面積分，可得到通過泄漏項(xiàng)耦合的徑向x-y平面與軸向z方向上的中子平衡方程：

w=x,y

(3)

k=1,…,K；g=1,…,G

(4)

由于材料不均勻多產(chǎn)生于徑向，所以在徑向上擴(kuò)散方程的求解采用基于共形映射的六邊形半解析節(jié)塊法，軸向上高階節(jié)塊展開法已經(jīng)能夠得到較高的計算精度且計算效率相對于半解析方法要高，所以在軸向上沿用高階節(jié)塊展開法。

1.1 徑向基于共形映射的六邊形半解析節(jié)塊法

在應(yīng)用于中子物理計算之前，共形映射被應(yīng)用于靜電學(xué)和流體力學(xué)將等勢曲線和場曲線從一種幾何圖形映射到另一種幾何圖形上。這種映射方法保留了原來坐標(biāo)下的角度信息，因此在新的拉普拉斯型微分方程中的拉普拉斯算子與共形映射前保持不變。這正是共形映射在靜電學(xué)與流體力學(xué)中有效應(yīng)用的原因，類似的，中子擴(kuò)散方程也是拉普拉斯微分方程，中子流是中子通量的矢量梯度，正如電場是電勢的矢量梯度一樣。中子通量和中子流之間的正交性在共形映射后也可得到保留。圖1為六邊形至矩形共形映射示意圖[5,10]。圖中：JLT為六邊形節(jié)塊徑向左側(cè)上表面中子凈流；JLB為六邊形節(jié)塊徑向左側(cè)下表面中子凈流；JRT為六邊形節(jié)塊徑向右側(cè)上表面中子凈流；JRB為六邊形節(jié)塊徑向右側(cè)下表面中子凈流。

在圖1左側(cè)復(fù)平面Z上，二維單群中子擴(kuò)散方程可表示為：

(5)

如圖1所示，復(fù)平面Z上的1組正交曲線坐標(biāo)(x,y)通過共形映射函數(shù)w=f(z)變換為復(fù)平面W上的1組正交曲線坐標(biāo)(u,v)，其中w=u+iv、z=x+iy為復(fù)變函數(shù)。z的實(shí)部和虛部與w的實(shí)部和虛部分別在復(fù)平面Z和W上組成直角坐標(biāo)系。整個x-y平面上的六邊形內(nèi)部區(qū)域被映射到u-v平面上的矩形內(nèi)部。矩形內(nèi)部相互垂直的直線U和V是六邊形內(nèi)部的相互垂直的曲線U和V在復(fù)平面W上的映射。式(5)中的拉普拉斯算子通過映射得到：

圖1 六邊形至矩形共形映射示意圖

(6)

其中，g(u,v)為共形映射的線型比例函數(shù)，其平方為映射區(qū)域面積比例函數(shù)，映射線型比例函數(shù)在映射區(qū)域內(nèi)始終為正，表示由于映射而產(chǎn)生的局部線性比例變化。值得注意的是，保形映射的局部線性尺度變化與方向無關(guān)，即映射在所有方向上對局部區(qū)域均是各向同性縮放。從式(6)可知，當(dāng)φ映射到W平面時，其空間梯度將被在所有方向上按1/g因子縮放?？芍?，共形映射后的中子擴(kuò)散方程形式如下：

(7)

通過對比式(5)與式(7)可發(fā)現(xiàn)，均勻六邊形節(jié)塊被轉(zhuǎn)化為均勻矩形節(jié)塊后，與裂變截面和吸收截面相乘的非均勻性因子g2是僅與幾何相關(guān)的映射面積比例函數(shù)。映射函數(shù)可通過Schwartz-Christoffel方法[11-12]得到，該映射是唯一的，且矩形形狀是固定的，具有特定的等高比。如圖2設(shè)映射后矩形節(jié)塊的高度為b，底部從-a/2到a/2。圖2中：JL為六邊形節(jié)塊徑向左側(cè)表面中子凈流；JR為六邊形節(jié)塊徑向右側(cè)表面中子凈流；JT為六邊形節(jié)塊徑向上表面中子凈流；JB為六邊形節(jié)塊徑向下表面中子凈流。

圖2 節(jié)點(diǎn)k內(nèi)部各方向中子流共形映射關(guān)系

將(7)沿v軸方向橫向積分，得到一維擴(kuò)散方程如下：

(8)

與式(3)類似，將式(8)左右兩端進(jìn)行處理，橫向積分方程可表示為：

(9)

其中：Σa0為體平均總移出截面；Qu(u)為有效中子源；Lu(u)為橫向泄漏項(xiàng)。

(10)

(11)

展開式中的5個展開系數(shù)可通過在u=±a/2處的邊界條件；橫向積分方程(式(9))滿足零次矩(W0)、一次矩(W1)、二次矩(W2)權(quán)重獲得。通過式(10)可得到在u=±a/2處的面平均中子凈流和源項(xiàng)、泄漏項(xiàng)的零次矩、一次矩、二次矩。對相鄰節(jié)塊重復(fù)同樣的展開方式可得到橫向一維凈流耦合方程的基本形式：

(12)

其中，Cn為群常數(shù)確定的耦合系數(shù)。將以上方法通過徑向掃描原六邊形節(jié)塊可得到3個橫向一維凈流耦合方程。對于二維徑向上總泄漏項(xiàng)Lu的近似需考慮其在徑向上其他方向的泄漏與軸向上的泄漏，即：

L(u)=Lxy(u)+Lz(u)

(13)

其中，徑向泄漏項(xiàng)Lxy(u)的近似采用階躍泄漏近似，如圖2所示：

(14)

本文軸向泄漏項(xiàng)Lz通過軸向使用的高階節(jié)塊展開法得到的軸向表面中子流計算得到。

1.2 軸向高階節(jié)塊展開法

軸向上偏中子通量采用高階節(jié)塊展開法處理，展開式與基函數(shù)如下：

(15)

(16)

(17)

(18)

(19)

圖3 徑向掃描方向

1.3 耦合方法

徑向的半解析方法與軸向使用的高階節(jié)塊展開法通過泄漏項(xiàng)進(jìn)行耦合，其耦合邏輯如圖4所示。圖4中：1)軸向求解一維擴(kuò)散方程得到軸向節(jié)塊表面中子流；2)徑向通過軸向表面中子流計算徑向擴(kuò)散方程中的軸向泄漏項(xiàng)；3)求解徑向中子凈流耦合的內(nèi)迭代；4)內(nèi)迭代收斂后更新軸向一維擴(kuò)散方程泄漏項(xiàng)中的3個方向上的徑向泄漏項(xiàng)；5)重復(fù)步驟1～4，直到徑向與軸向上中子泄漏耦合收斂后，更新徑向、軸向偏通量與相應(yīng)源項(xiàng)系數(shù)完成源迭代[13]。

圖4 耦合方法流程圖

2 基準(zhǔn)題驗(yàn)證結(jié)果

通過對原中子物理程序HNHEX求解器的修改并加入了中子泄漏項(xiàng)耦合模塊，開發(fā)了能用于求解六邊形三維中子擴(kuò)散方程的中子物理程序SA-HNHEX。本文采用VVER-440的二維和三維基準(zhǔn)題對該程序進(jìn)行驗(yàn)證，除基準(zhǔn)題參考解外，加入了使用節(jié)塊展開法的DIF3D-N與使用半解析節(jié)塊法的ANC-HM計算結(jié)果進(jìn)行對比。

VVER-440基準(zhǔn)題是由Seidel等于1984年提出的[6]，其堆芯直徑上有25盒組件，組件對邊距為14.7 cm。反射層組件外部邊界條件考慮為真空。整個堆芯有7盒控制棒組件，當(dāng)控制棒組件插入堆芯時，會將下部的燃料組件從堆芯下部推出，從而帶來更大的負(fù)反應(yīng)性反饋，在二維基準(zhǔn)題中，控制棒組件為全部插入堆芯，在三維基準(zhǔn)題中為插入1/2。表1列出VVER-440基準(zhǔn)題所使用的群常數(shù)，圖5為二維VVER-440基準(zhǔn)題中1/12堆芯布置情況及SA-HNHEX計算結(jié)果相對于參考值(DIF3D-FD采用600和864三角形/六邊形網(wǎng)格推測得到[6])的相對偏差。

圖5 VVER-440二維基準(zhǔn)題堆芯布置及真空邊界條件計算結(jié)果

表1 VVER-440 基準(zhǔn)題材料群常數(shù)

表2列出VVER-440二維基準(zhǔn)題計算結(jié)果對比。如表2所列，SA-HNHEX計算VVER-440二維基準(zhǔn)題本征值keff為1.016 01，與基準(zhǔn)參考值誤差為63 pcm，與同樣使用共形映射解析節(jié)塊方法的ANC-H相比，精度相當(dāng)。與使用高階節(jié)塊展開法的HNHEX相比，在節(jié)點(diǎn)功率最大誤差與本征值誤差上計算精度均有明顯提高，但在計算時間上有所增加。

表2 VVER-440二維基準(zhǔn)題計算結(jié)果對比

VVER-440三維基準(zhǔn)題使用的群常數(shù)及堆芯徑向布置方案均與二維基準(zhǔn)題相同，軸向上控制棒組件為半插入堆芯，如圖6所示。堆芯高度為250 cm，堆芯上下均增加了25 cm反射層，軸向、徑向上的邊界條件均考慮為真空條件。

圖6 VVER-440三維基準(zhǔn)題控制棒軸向位置

表3列出VVER-440三維基準(zhǔn)題計算結(jié)果對比。如表3所列，SA-HNHEX與HNHEX均采用了兩組軸向節(jié)點(diǎn)數(shù)(12節(jié)點(diǎn)與24節(jié)點(diǎn))對VVER-440三維基準(zhǔn)題進(jìn)行計算。由于SA-HNHEX在徑向上采用了更精確的半解析節(jié)塊法，其計算精度有顯著提升，在12個軸向計算節(jié)點(diǎn)的情況下已達(dá)到HNHEX使用24軸向計算節(jié)點(diǎn)的精度。在計算耗時上，由于SA-HNHEX在軸向上使用更快的高階節(jié)塊展開法，相較于ANC-H在保證計算精度相當(dāng)?shù)那闆r下有更好的計算速度表現(xiàn)；相較于HNHEX軸向24個計算節(jié)點(diǎn)，其計算耗時的小幅增加所帶來的精度上的顯著提高也體現(xiàn)了該方法較高的性價比。

表3 VVER-440三維基準(zhǔn)題計算結(jié)果對比

3 結(jié)論與展望

本文基于先進(jìn)節(jié)塊方法提出了一種綜合徑向共形映射半解析-軸向高階節(jié)塊展開法的三維六邊形中子擴(kuò)散方程求解方法，通過對原六邊形高階節(jié)塊展開法求解器HNHEX的修改開發(fā)了SA-HNHEX三維六邊形中子擴(kuò)散方程求解程序，通過使用VVER-440型反應(yīng)堆的二維、三維基準(zhǔn)題對該方法及程序進(jìn)行了初步驗(yàn)證。計算結(jié)果與參考值符合較好，在計算精度上與同類程序相當(dāng)，并在計算精度上相比原程序HNHEX有顯著提升，且在計算耗時的增加上較為友好。

在未來的工作中，將對SA-HNHEX的迭代進(jìn)行加速優(yōu)化，并將該求解器與快堆系統(tǒng)分析程序SAC[14-16]進(jìn)行耦合開發(fā)與驗(yàn)證。