莊 坤，王連杰，劉 琨，顏江濤，尚 文

(1.南京航空航天大學(xué) 材料科學(xué)與技術(shù)學(xué)院，江蘇 南京 211106；2.中國核動力研究設(shè)計院 核反應(yīng)堆系統(tǒng)設(shè)計技術(shù)重點實驗室，四川 成都 610041)

隨著核能的快速發(fā)展，基于六角形組件的先進(jìn)堆型不斷涌現(xiàn)，如快中子反應(yīng)堆[1]、可再生沸水堆[2]、熔鹽堆[1]等，與壓水堆相比，這些新堆型往往采用更加先進(jìn)的燃料類型、更加不均勻的燃料布置，其中子能譜與壓水堆有很大區(qū)別，中子各向異性程度更強(qiáng)，傳統(tǒng)輕水堆計算分析方法中的擴(kuò)散近似假設(shè)不再適用。而堆芯燃料管理計算、中子動力學(xué)計算是先進(jìn)堆中子學(xué)設(shè)計過程中重要的內(nèi)容，在先進(jìn)堆的設(shè)計階段，需要反復(fù)大量的計算堆芯的功率分布等參數(shù)。如果利用中子輸運理論，雖然能獲得較好的結(jié)果但計算效率較低，如果采用傳統(tǒng)的擴(kuò)散理論，雖然可提高計算效率，但堆芯中復(fù)雜的中子能譜和強(qiáng)的中子各向異性等特性使得擴(kuò)散計算的精度較低。因此有必要研究一種既能滿足精度需求，又能滿足計算效率的先進(jìn)堆芯中子學(xué)計算方法。

Quasi-diffusion是一種低階中子輸運方程，于1964年首次被提出[3]，其形式與傳統(tǒng)中子擴(kuò)散方程相近，不同的是泄漏項表達(dá)方式。傳統(tǒng)擴(kuò)散理論建立的基本假定之一是中子的角通量是關(guān)于角度的一階函數(shù)，而Quasi-diffusion方程則不采用這一假設(shè)，用一個艾丁頓因子表示中子角通量與角度向量的乘積的積分項。艾丁頓因子使得Quasi-diffusion方程具有一些中子流各向異性的輸運特性，這是傳統(tǒng)擴(kuò)散理論中所忽略的，由此也可以推斷，Quasi-diffusion方程的解的精度要高于傳統(tǒng)擴(kuò)散方程。近年來，國內(nèi)外許多學(xué)者展開了更深入的研究，并探索其在工程中的應(yīng)用。2009年，William[4]基于2D Quasi-diffusion方程求解2D笛卡爾幾何中非結(jié)構(gòu)網(wǎng)格中的輸運問題，并給出了一維幾何下的艾丁頓因子解析解。2013年，Aristova、Baydin等[5-6]基于Quasi-diffusion方程求解多群輸運問題，并計算了三維六邊形組件快中子堆臨界參數(shù)。2013年，Razvan[7]基于節(jié)塊法和有限體積方法研究了矩形幾何和六角形幾何下的2維2群Quasi-diffusion方程的計算方法，并將其應(yīng)用于高溫氣冷堆中子通量和keff計算，其中艾丁頓因子的計算通過中子角通量的球諧函數(shù)展開進(jìn)行。2015年，Hall[8]利用修改后的SERPENT程序產(chǎn)生組件均勻化截面和艾丁頓因子，基于節(jié)塊展開法研究了一維情況下的Quasi-diffusion方程的數(shù)值計算方法，并推導(dǎo)了節(jié)塊不連續(xù)因子。2018年，Anistratov等[9]基于間斷有限元方法和Quasi-diffusion方程研究了一維平板幾何中輸運問題求解。2019年，Ghassemi等[10]基于Quasi-diffusion方程研究了2D矩形幾何下多群熱輻射傳輸問題求解方法。2019年，Jones等[11]基于局部網(wǎng)格離散化技術(shù)提出了一種2DR-z坐標(biāo)系下采用Quasi-diffusion方程求解穩(wěn)態(tài)非結(jié)構(gòu)網(wǎng)格輸運問題。2019年，Xu等[12]基于三維截面和Quasi-diffusion方程對TREAT堆芯進(jìn)行分析。2019年，Olivier等[13]研究了利用Quasi-diffusion加速SN高階輸運方程的計算方法。

可以看出，目前大多數(shù)研究將Quasi-diffusion方程應(yīng)用于加速中子輸運計算，艾丁頓因子通過高階中子輸運方程解進(jìn)行計算，而將Quasi-diffusion方程作為一種像傳統(tǒng)中子擴(kuò)散方程一樣的獨立方程應(yīng)用于反應(yīng)堆計算中的研究還相對較少，主要原因是目前組件程序還不能提供艾丁頓因子。同時，國內(nèi)外對艾丁頓因子的計算主要采用近似方法，并針對1D和2D矩形幾何開展相關(guān)研究，而對三維六角形幾何的Quasi-diffusion方程的求解方法研究還有待進(jìn)一步加強(qiáng)。針對艾丁頓因子計算的非線性特點，本文基于兩步法的思想，將艾丁頓因子看作是一個特殊的少群參數(shù)，采用中子能譜適應(yīng)性更好的蒙特卡羅程序SERPENT計算?？紤]到變分節(jié)塊法的優(yōu)勢，利用三維節(jié)塊內(nèi)的正交基函數(shù)對節(jié)塊內(nèi)相關(guān)變量進(jìn)行展開，避免傳統(tǒng)橫向積分技術(shù)在處理六角形Quasi-diffusion方程時橫向泄漏產(chǎn)生奇異項[14]。

本文對傳統(tǒng)擴(kuò)散變分節(jié)塊法進(jìn)行拓展，開發(fā)六角形組件堆芯計算程序VNMQD，并采用3D VVER1000基準(zhǔn)題、1D RBWR問題、3D BN600簡化模型對程序進(jìn)行驗證。

1 理論模型與數(shù)值方法

1.1 Quasi-diffusion方程及邊界條件

艾丁頓因子張量E表示如下：

(1)

(2)

Quasi-diffusion方程中的連續(xù)條件是基于中子角通量ψ在兩種介質(zhì)交界面上都是空間r的連續(xù)函數(shù)得到的，可用下面1組近似條件表示：

(3)

式中：Yn，m(Ω)為球諧函數(shù)；n為展開階數(shù)，m=-n,…,n；n為交界面法線方向向量；ψ(r，Ω，E)為r處、Ω方向、能量為E的中子角通量密度。

經(jīng)過化簡即可得到擴(kuò)散近似方程的連續(xù)條件。

流連續(xù)邊界條件為：

(4)

通量連續(xù)條件為：

n·Ex(r,E)φ(r,E)

n·Ey(r,E)φ(r,E)

n·Ez(r,E)φ(r,E)

(5)

式中：J(r，E)為r處能量為E的中子流密度；φ(r，E)為r處能量為E的中子標(biāo)通量。

Quasi-diffusion方程中除包含傳統(tǒng)的宏觀截面(如總截面、裂變截面、散射截面)外，還包含一個重要的量——艾丁頓因子。艾丁頓因子具有以下特點：式(1)中所有元素都在0和1之間，對角元素在1/3附近，非對角元素約等于0，對角元與非對角元有幾個數(shù)量級的差別。值得指出的是，如果艾丁頓因子在每個區(qū)域、每個能群及每個方向上均為1/3，則Quasi-diffusion方程退化為擴(kuò)散方程，此時意味著中子流是各向同性的。

目前組件程序還不能計算艾丁頓因子，考慮到傳統(tǒng)熱中子反應(yīng)堆組件程序在復(fù)雜能譜反應(yīng)堆中應(yīng)用的限制，以及蒙特卡羅方法具有幾何和中子能譜適應(yīng)性強(qiáng)的特點，本文基于SERPENT進(jìn)行組件均勻化計算。原版SERPENT不具備艾丁頓因子的計算功能，經(jīng)過二次開發(fā)后可實現(xiàn)艾丁頓因子的計算，并得到了成功的應(yīng)用[8]。由式(2)可看出，計算艾丁頓因子首先需要知道中子角通量分布，而中子角通量正是需要求解的，因此這是一個非線性問題。通過分析發(fā)現(xiàn)，艾丁頓因子與宏觀截面具有類似的非線性特點，本文基于兩步法的思想，將艾丁頓因子看作是一個特殊的宏觀參數(shù)，在SERPENT中通過修改其源代碼統(tǒng)計每個粒子的飛行方向Ωu、Ωv、Ωw以及沿這個方向的中子角通量變量，并基于式(5)分別計算分子、分母最終得到艾丁頓因子，然后通過并群并區(qū)計算得到均勻化區(qū)域的少群艾丁頓因子。

由于艾丁頓因子張量對角元與非對角元有幾個數(shù)量級的差別，本文在實際應(yīng)用中忽略了非對角元素，因此三維多群Quasi-diffusion方程表示如下：

Σsgφg+Sg

(6)

式中：g為能群編號；Σtr,g為能群g的輸運截面；φg為中子標(biāo)通量密度，cm-2·s-1；Σtg和Σsg分別為中子宏觀總截面和群內(nèi)宏觀散射截面，cm-1；Euug為少群艾丁頓因子；Sg為中子源項，cm-3·s-1，包括散射源項和裂變源項：

(7)

Quasi-diffusion方程中的中子流Jg表達(dá)式為：

(8)

式中，Exxg、Eyyg、Ezzg為能群g的艾丁頓因子張量對角線元素。

1.2 六角形幾何Quasi-diffusion方程變分節(jié)塊法求解

本文從三維多群Quasi-diffusion方程出發(fā)，利用Galerkin變分和Lagrange乘子法在整個求解域上建立一個包含節(jié)塊平衡關(guān)系和邊界條件的泛函，再采用Ritz離散方法將這個泛函進(jìn)行展開，相應(yīng)的基函數(shù)為空間上的正交多項式函數(shù)，得到耦合節(jié)塊體內(nèi)的中子通量密度和節(jié)塊邊界上的分中子流密度的節(jié)塊響應(yīng)矩陣，最后通過迭代計算實現(xiàn)對方程的數(shù)值求解。

根據(jù)Galerkin變分原理，Quasi-diffusion方程可在由若干節(jié)塊組成的整個求解域及其邊界上建立全局泛函F：

(9)

式中，節(jié)塊v的貢獻(xiàn)為：

(10)

式中：dv表示對節(jié)塊v的空間積分；γ為邊界面；χγg為節(jié)塊邊界上的凈中子流密度。

對式(10)中的中子標(biāo)通量密度φg、中子源項Sg和節(jié)塊邊界上的凈中子流密度χg分別利用六角形節(jié)塊內(nèi)的正交基函數(shù)進(jìn)行如下展開，為簡化省去了能群編號g。

(11)

式中，空間基函數(shù)fi、hkγ為完全的正交多項式：

(12)

式中，δij為克羅內(nèi)克符號。

將式(11)代入式(10)可得到節(jié)塊泛函與中子通量密度、中子源和凈中子流密度的展開系數(shù)的關(guān)系。

Fv[φ,χ]=φTAφ-2φTs+2φTMχ

(13)

式中：φ、s和χ為中子通量密度、中子源和凈中子流密度的展開系數(shù)組成的向量；A和M為矩陣，其計算公式為：

Aij=Pij+δijVvΣr

M=[M1,M2,M3,…]

(14)

隨后采用與文獻(xiàn)[14]相同的處理方式可得到節(jié)塊分中子流J+和J-、節(jié)塊平均中子通量密度φ之間的響應(yīng)關(guān)系：

J+=Bs+RJ-

φ=Hs-C(J+-J-)

(15)

式中，B、R、H、C為幾何與材料共同決定的響應(yīng)矩陣。

式(15)表示六角形節(jié)塊之間是通過表面的分中子流耦合的，在Quasi-diffusion方程處理時應(yīng)充分考慮式(4)和(5)的連續(xù)性條件。在如圖1所示的六角形節(jié)塊中，對于垂直于坐標(biāo)軸的面1、面2，其法線方向的凈中子流密度由式(16)表示，對于非垂直于坐標(biāo)軸的面3、4、5、6，其法線方向的凈中子流密度由式(17)表示，相應(yīng)的通量密度連續(xù)性條件如(18)、(19)所示。

(16)

(17)

(18)

(19)

圖1 六角形節(jié)塊示意圖Fig.1 Scheme of hexagonal nodal

通過式(18)、(19)及分中子流與凈中子流的關(guān)系，很容易得到相鄰節(jié)塊在交界面處分中子流的修正表達(dá)式：

(20)

fu=Exx對于1,2面

fu=Exxcos2α+Eyycos2β

對于 3, 4, 5, 6面

2 數(shù)值結(jié)果

所有節(jié)塊通過式(15)耦合，通過紅黑掃描策略基于傳統(tǒng)的源迭代方法可對整個系統(tǒng)進(jìn)行求解計算?；谏鲜隼碚撃Ｐ?，本文采用FORTRAN90語言開發(fā)了六角形變分節(jié)塊法Quasi-diffusion程序VNMQD，程序計算流程如圖2所示。隨后采用六角形純擴(kuò)散基準(zhǔn)題如2D/3D VVER系列問題對VNMQD進(jìn)行驗證，由于該問題為宏觀問題，在VNMQD進(jìn)行計算時設(shè)置艾丁頓因子Exx=Eyy=Ezz=1/3。

圖2 計算流程圖Fig.2 Flowchart of calculation

為簡化篇幅，本文只給出不帶反射層3D VVER1000基準(zhǔn)題計算結(jié)果。該基準(zhǔn)題堆芯軸向高度為200 cm，包含8圈燃料組件(其中包括25根控制棒)，呈1/6旋轉(zhuǎn)對稱，徑向和軸向反照率分別為0.125、0.15，組件對邊距為23.6 cm。VNMQD計算時軸向高度劃分為10 cm 1個網(wǎng)格，詳細(xì)的基準(zhǔn)題幾何及宏觀截面數(shù)據(jù)參考文獻(xiàn)[15-16]。VNMQD計算的功率分布及堆芯有效增殖因數(shù)與參考解的比較如圖3、表1所示，其中參考解取自文獻(xiàn)[15]，對比結(jié)果均來自擴(kuò)散，計算采用0階散射截面。對比結(jié)果可以發(fā)現(xiàn)，VNMQD組件功率計算結(jié)果與參考解吻合較好，最大相對功率偏差出現(xiàn)在最外層組件為0.3%；VNMQD計算的keff與參考解符合很好，誤差為-11 pcm，且與國際上同類擴(kuò)散程序AFEN、ANC-H、GTDIF-H、NLSANM具有相同的計算精度。當(dāng)艾丁頓因子為1/3，Quasi-diffusion方程退化為擴(kuò)散理論，該系列基準(zhǔn)題驗證表明VNMQD程序在艾丁頓因子為1/3時編寫正確。

圖3 VVER1000基準(zhǔn)題組件歸一化功率分布Fig.3 Assembly normalized power distribution of VVER1000 benchmark

表1 VVER1000基準(zhǔn)題有效增殖因數(shù)Table 1 Effective multiplication factor of VVER1000 benchmark

為進(jìn)一步驗證VNMQD的正確性，本文設(shè)計了一個1D RBWR單組件問題和3D BN600簡化堆芯問題。為盡可能避免由于均勻化模型帶來的均勻化誤差，特別是艾丁頓因子的均勻化誤差，本文采用了全堆芯計算的均勻化模型獲得相關(guān)均勻化截面及艾丁頓因子。實際上，選取合理的均勻化模型、邊界條件等可以避免全堆芯均勻化模型的使用，以增加實用，在這里不作討論。RBWR單組件問題來源于日本可再生沸水堆設(shè)計，該問題具有較強(qiáng)的非均勻性[17]，相關(guān)設(shè)計參數(shù)及材料組成參考文獻(xiàn)[8]，如圖4所示，外邊為全反射邊界。由于SERPENT自身畫圖原因，圖4描述的是基于無限組件排布的單組件矩形圖。該問題中組件對邊距為19.766 8 cm，徑向全反射軸向為真空邊界條件，軸向高度為134.3 cm，VNMQD在計算中將軸向從下到上分為34層，高度分別為5×5.6 cm、8×2.412 5 cm、8×6.5 cm、8×3.5 cm、5×1.4 cm。采用SERPENT產(chǎn)生12群少群均勻化截面及艾丁頓因子，能群劃分列于表2。每層網(wǎng)格產(chǎn)生1套均勻化參數(shù)，即34種均勻化截面，VNMQD計算采用0階散射截面，并對總截面進(jìn)行1階散射修正。SERPENT的非均勻計算結(jié)果作為參考解，對比包括堆芯有效增殖因數(shù)和節(jié)塊平均中子通量密度。

表2 RBWR問題12群能群結(jié)構(gòu)Table 2 Energy structure of 12 energy groups for RBWR benchmark

1D RBWR問題的VNMQD計算結(jié)果列于表3，其中no edd情況表示為在計算中艾丁頓因子為1/3，即轉(zhuǎn)化為標(biāo)準(zhǔn)中子擴(kuò)散計算，其他少群截面則與edd情況保持一致。由表3可見，相比于VNMQD(no edd)，VNMQD(edd)計算的堆芯keff與參考誤差更小，為-73 pcm，而VNMQD(no edd)誤差達(dá)到了-1 738 pcm，兩者計算時間相當(dāng)，約為6.2 s。同樣對比了各節(jié)塊中子通量密度分布及誤差，由于趨勢類似，本文只給出第1、4群的對比，如圖5所示。由圖5可見，VNMQD(edd)計算結(jié)果與參考值符合更好，第1、4群中子通量密度最大相對誤差不超過10%，而VNMQD(no edd)結(jié)果最大相對誤差可以達(dá)到160%。

表3 RBWR問題真空邊界條件下VNMQD計算的keff對比Table 3 keff comparison of RBWR benchmark under axial vaccum condition

a——第1群；b——第4群圖5 RBWR問題節(jié)塊中子通量密度對比Fig.5 Comparison of nodal neutron flux density for RBWR benchmark

為驗證VNMQD在中子各向異性較強(qiáng)的六角形堆芯中的應(yīng)用，本文基于BN600快堆基準(zhǔn)題[18]中的燃料組件(LEZ、HEZ)和控制組件(SHR、SCR)建立了一個簡化的具有12圈燃料構(gòu)成的堆芯模型，如圖6所示，其中組件內(nèi)的相關(guān)材料進(jìn)行了體積權(quán)重打混處理，軸向分層與基準(zhǔn)題所述一致，軸向及徑向為真空邊界條件。采用SERPENT進(jìn)行精細(xì)建模，基于11少群能群結(jié)構(gòu)產(chǎn)生每個節(jié)塊的均勻化截面及艾丁頓因子，VNMQD計算時軸向高度從底到頂分別為30、4.5、15、4.5、23、5.3、41.15、5.1、41.15、5.5、29.7、30 cm。

圖6 BN600簡化模型徑向和軸向布置Fig.6 Radial and axial configurations of BN600 simplified model

堆芯keff計算結(jié)果列于表4，SERPENT計算結(jié)果作為參考解，VNMQD計算采用0階散射截面，并對總截面進(jìn)行1階散射修正。由表4可見，與參考解相比，VNMQD(edd)計算結(jié)果要比VNMQD(no edd)誤差更小，僅為11 pcm，而VNMQD(no edd)誤差可以達(dá)到-353 pcm，而兩者的計算效率接近，計算時間約為3 s。同樣對比了不同方法計算的組件平均中子通量密度，如圖7所示。由圖7可見，VNMQD(edd)計算結(jié)果要比VNMQD(no edd)誤差更小，在靠近邊界處VNMQD(edd)計算精度更高。

表4 BN600問題真空邊界條件下VNMQD計算的keff對比Table 4 keff comparison of VNMQD result for BN600 benchmark at vaccum condition

a——第1群；b——第4群圖7 BN600問題節(jié)塊中子通量密度對比Fig.7 Comparison of nodal neutron flux density for BN600 benchmark

3 結(jié)論

隨著核能技術(shù)的快速發(fā)展，更加非均勻布置、更強(qiáng)的中子各向異性六角形組件堆型不斷涌現(xiàn)，Quasi-diffusion方程比傳統(tǒng)擴(kuò)散方程引入更少的近似，計算精度更高。本文基于兩步法的概念，將Quasi-diffusion方程中的艾丁頓因子看作是一個特殊的少群參數(shù)，采用幾何和中子能譜適應(yīng)性更好的蒙特卡羅程序SERPENT進(jìn)行計算，拓展傳統(tǒng)擴(kuò)散變分節(jié)塊法實現(xiàn)Quasi-diffusion方程的求解，并開發(fā)了相應(yīng)的堆芯中子通量密度計算程序VNMQD。采用基準(zhǔn)題3D VVER1000、RBWR單組件、3D BN600簡化模型對程序進(jìn)行了驗證，結(jié)果表明，VNMQD具有較高的計算精度，對于非均勻性較強(qiáng)的堆芯，Quasi-diffusion比傳統(tǒng)擴(kuò)散在堆芯有效增殖因數(shù)和中子通量密度方面計算精度均有巨大提升，而兩者的計算效率非常接近。VNMQD開發(fā)正確，實現(xiàn)了計算精度和計算效率的平衡。