尚 文，莊 坤,2,*，李 婷，湯曉斌,2，于東森

(1.南京航空航天大學 材料科學與技術(shù)學院，江蘇 南京 211106；2.南京航空航天大學 空間核技術(shù)應用與輻射防護工業(yè)和信息化部重點實驗室，江蘇 南京 211106)

變分節(jié)塊法(VNM)無需使用橫向積分技術(shù)即可將節(jié)塊的體積通量矩和偏中子流矩擴展為正交基函數(shù)之和，是核反應堆堆芯中子學計算最成功的節(jié)塊法之一，由西北大學和阿貢國家實驗室(ANL)首次提出，用于求解多群穩(wěn)態(tài)中子擴散和輸運方程。20世紀90年代中期，ANL開發(fā)的VARIANT是第一個基于VNM的程序[1]，它被用于ANL(如REBUS代碼)[2]和歐洲(如ERANOS代碼)[3-5]的快堆設計。2007年，ANL又開發(fā)了名為NODAL的VARIANT程序新版本，作為UNC軟件包中的求解器之一[6-7]。2011年，愛達荷州國家實驗室(INL)的INSTANT程序也采用了變分節(jié)塊法[8]。1995年，在DIF3D代碼中增加了VARIANT用以提高快堆的通量解[9]。2014年，西安交通大學的李云召等開發(fā)了適用于三維矩形幾何的反應堆堆芯計算程序VIOLET[10]，后被用于壓水堆堆芯計算的Bamboo-Core代碼[11-13]的中子擴散模塊。此外，VIOLET代碼得到擴展以求解六角形幾何的中子輸運方程，被用作快堆計算代碼系統(tǒng)NECP-SARAX的解算裝置[14]。2018年，Zhang等[15]提出了基于矩形網(wǎng)格的改進變分節(jié)塊法來求解三維穩(wěn)態(tài)多群中子輸運方程；2019年又將積分法應用于具有六角形組件的反應堆，提出了基于六角形網(wǎng)格的改進變分節(jié)塊法[16]。

有限元法和基于任意三角形網(wǎng)格的變分節(jié)塊法均采用Galerkin變分技術(shù)、泛函概念和非結(jié)構(gòu)網(wǎng)格。其不同之處在于：1) 有限元法一般采用小尺寸網(wǎng)格劃分和低階試驗函數(shù)，而變分節(jié)塊法可使用高階(>4)多項式函數(shù)來展開任意三角形網(wǎng)格內(nèi)的空間變量；2) 有限元法計算有限元節(jié)塊的展開系數(shù)，通過傳遞有限元節(jié)塊值實現(xiàn)耦合，變分節(jié)點法則利用三維高階(>4)多項式函數(shù)直接展開任意三角形網(wǎng)格內(nèi)的空間變量，通過掃描并利用節(jié)塊網(wǎng)格界面上凈中子流和中子通量的連續(xù)條件實現(xiàn)耦合；3) 對于變分節(jié)塊法，可直接由節(jié)塊內(nèi)中子通量的展開系數(shù)計算節(jié)塊的精細功率分布，而有限元法則需根據(jù)有限元節(jié)塊的展開系數(shù)進一步計算。

目前大多數(shù)變分節(jié)塊法基于六角形或矩形結(jié)構(gòu)幾何網(wǎng)格，然而隨著核能與核技術(shù)的發(fā)展，新概念堆型被不斷提出，這些新型反應堆的一個重要特點是其組件設計和堆芯布置不再采用單一常規(guī)的幾何形狀和堆芯結(jié)構(gòu)，因此基于矩形和六角形幾何的變分節(jié)塊方法不能準確計算新型反應堆堆芯中子學?？紤]到三角形網(wǎng)格對曲線或多角形邊界有很好的擬合性，理論上可用來逼近任意幾何形狀。另外，采用三角形網(wǎng)格進行計算時，可進行局部的網(wǎng)格加密，從而提高計算結(jié)果的準確性和可靠性。因此，本文擬研究基于任意三角形網(wǎng)格的多群中子擴散變分節(jié)塊法。首先，對計算區(qū)域進行三角形網(wǎng)格剖分，將任意三角形變換為正三角形；其次，建立泛函并利用正三角形內(nèi)正交基函數(shù)展開節(jié)塊內(nèi)參量；再次，利用變分原理獲得中子通量密度與節(jié)塊邊界上分中子流的響應關(guān)系；最后，基于傳統(tǒng)源迭代法對其進行求解?；谏鲜隼碚撃Ｐ烷_發(fā)程序TriVNM，并采用基準題驗證程序TriVNM的可靠性與精確性，以及對具有復雜幾何組件堆芯的適用性。

1 理論模型

首先，從多群中子擴散方程出發(fā)，有：

Σs,gΦg(r)+Sg

(1)

其中：g為能群編號；Φg(r)為中子標通量密度，cm-2·s-1；Dg為中子擴散系數(shù)，cm；Σt,g為中子宏觀總截面，cm-1；Σs,g為群內(nèi)宏觀散射截面，cm-1；Sg為中子源項，cm-3·s-1，包括散射源項和裂變源項。

(2)

其中，g′為能群編號。

其次，利用Galerkin變分技術(shù)，對中子擴散方程在整個求解域上建立一個包含三角形節(jié)塊中子平衡方程的泛函F：

(3)

其中，節(jié)塊v的貢獻為：

(4)

χγ=J·nγ

(5)

(6)

其中：nγ為邊界γ的外法線方向向量；χγ為節(jié)塊邊界上沿外法線方向的凈中子流密度。

再次，采用坐標變換將ANSYS軟件對計算區(qū)域剖分的任意三角形節(jié)塊變換為成正三角形節(jié)塊，如圖1所示。

圖1 任意三角形坐標變換示意圖Fig.1 Arbitrary trianglar coordinate transformation

任意三角形和正三角形的坐標轉(zhuǎn)換關(guān)系為：

(7)

坐標變換后進行基函數(shù)的構(gòu)造。線性無關(guān)的函數(shù){1，x，y，z，x2，xy，xz，y2，…}構(gòu)成函數(shù)向量g(r)：

g(r)=[1，x，y，z，x2，…]T

(8)

顯然，該函數(shù)系在積分區(qū)域內(nèi)不是正交歸一的。而正定對稱矩陣G一定與單位矩陣相似，即存在滿秩矩陣Q并滿足下式：

(9)

QGQT=I

(10)

令f=Qg，滿秩矩陣Q保證g(r)和G兩組基函數(shù)所張成的函數(shù)空間相同，即：

span{f1,f2,…}=span{g1,g2,…}

(11)

(12)

矩陣Q可用Gram-Schmidt算法求解。通過上述坐標變換，即可構(gòu)造出正交歸一的基函數(shù)。

將式(4)中節(jié)塊內(nèi)中子標通量密度Φg(r)、中子源項Sg和節(jié)塊邊界上的凈中子流密度χγ分別按照正三角形的基函數(shù)展開，其中，空間基函數(shù)fi和hkγ為完全的正交多項式。

(13)

(14)

于是，中子標通量密度和源項的展開式系數(shù)之間滿足以下關(guān)系式：

(15)

將式(13)代入式(4)，得到：

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

(16)

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

Aij=(3Σtr)-1Pij+δijVvΣr

(17)

(18)

M=[M1M2…Mγ…]

(19)

(20)

分別令關(guān)于φT和χγ的一階變分形式為0，可得：

φ=A-1s-A-1Mχ

(21)

(22)

Ψ=MTA-1s-MTA-1Mχ

(23)

為將響應矩陣表達成通用形式，定義邊界上的分中子流密度為：

(24)

Ψ=2(j++j-)

(25)

χ=j+-j-

(26)

將式(25)、(26)代入式(24)，得到離散后的中子擴散方程：

j+=Bs+Rj-

(27)

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

(28)

式中：j+和j-分別為節(jié)塊邊界上的出射流中子密度矩向量和入射流中子密度矩向量；B、C、H和R為響應矩陣。

最后，根據(jù)上述推導，得到變分節(jié)塊法中的節(jié)塊響應關(guān)系：

(29)

(Ij-RgП)jg=Bgsg

(30)

φg=Hgsg-Cg(Ij-Π)jg

(31)

式中：jg為節(jié)塊的出射中子流密度向量(省去了上標“+”)，cm-2·s-1；Bg、Cg、Hg和Rg為由幾何與材料共同決定的響應矩陣；Ij為大小與分中子流密度向量相對應的單位矩陣；Π為節(jié)塊分中子流密度之間的聯(lián)系矩陣，包含了內(nèi)部邊界上相鄰節(jié)塊間的互為出入射關(guān)系和外邊界上的邊界條件。

將式(29)代入式(30)、(31)，可得：

(32)

(33)

將所有能群合并，可得到算子形式的變分節(jié)塊法：

(34)

(35)

Fgg′=χgνΣf,g′

(36)

所有三角形節(jié)塊通過中子流相互耦合，然后利用傳統(tǒng)的源迭代法對其進行求解。含特征值的裂變源問題采用冪法迭代求解，稱為裂變源迭代。在每次裂變源迭代中，需在已知裂變源f的情況下，進行多群迭代。

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

通過以上理論模型，本文基于FORTRAN90語言開發(fā)了基于任意三角形網(wǎng)格的變分節(jié)塊法中子擴散計算程序TriVNM。分別采用2D/3D-IAEA、2D/3D-LRA、2D/3D-VVER440、2D/3D-VVER1000、不規(guī)則幾何基準題驗證程序TriVNM正確性。本文只選擇矩形幾何組件的2D-IAEA、六角形幾何組件不帶反射層的3D-VVER1000和非結(jié)構(gòu)幾何組件這3個具有代表性的基準題進行說明。以下所有計算均基于Intel(R) Core(TM) i7-7700 CPU@3.60 GHz。

2.1 2D-IAEA基準題

2D-IAEA基準問題[17]是一個簡化的兩群PWR基準問題，堆芯共有177盒燃料組件，組件的幾何尺寸為20 cm×20 cm，堆芯1/8對稱分布，堆芯布置和網(wǎng)格剖分如圖2所示。堆芯按雙區(qū)布料方案布置，徑向有20 cm厚的水反射層，堆芯外邊界條件為真空邊界(即入射中子流密度為0)，所有組件均采用等效均勻化的參數(shù)。由于此基準問題堆芯布置了控制棒，而且在堆芯和反射層交界面上，熱中子通量梯度很大，因此廣泛用于校核雙群中子擴散方程的數(shù)值計算模型及其計算精度。

圖2 2D-IAEA基準題1/4堆芯布置示意圖(a)和網(wǎng)格剖分示意圖(b)Fig.2 Configuration of 2D-IAEA benchmark (a) and its meshing by ANSYS code (b)

TriVNM計算的有效增殖因數(shù)keff為1.029 560 7，2D-IAEA基準題給出的keff參考值為1.029 585?？煽闯觯c參考值相比，有效增殖因數(shù)的計算偏差僅為2.360 pcm。圖3為用程序TriVNM計算的歸一化功率分布及與參考值的比較，可看出，歸一化功率的最大偏差出現(xiàn)在堆芯最內(nèi)圈的組件上，約為0.20%。表1為TriVNM計算結(jié)果與其他程序計算結(jié)果的比較，可看出，TriVNM計算結(jié)果符合參考解且較其他程序的偏差小。

表1 2D-IAEA基準題在不同程序下的計算結(jié)果Table 1 Calculation results of 2D-IAEA benchmark problem with different programs

圖3 歸一化功率分布計算結(jié)果及與參考解的對比Fig.3 Calculation result of normalized power distribution and comparison with reference result

2.2 3D-VVER1000基準題

3D-VVER1000基準題是不帶反射層的2D-VVER1000基準問題[17]的擴展，徑向上有8圈燃料組件，全堆芯共有25根控制棒，堆芯布置呈1/6旋轉(zhuǎn)對稱，組件的對邊距為23.6 cm。堆芯軸向高度為200 cm，25根控制棒中的6根提到堆芯中部。堆芯的徑向外邊界反照率為0.125，軸向反照率為0.15。TriVNM計算時，徑向上將每個六角形節(jié)塊劃分為8個三角形，軸向上每10 cm劃分為1層，參考解仍是由細網(wǎng)差分程序DIF3D-FD在徑向上將每個六角形組件劃分為384和486個三角形、軸向上每5.0 cm為1層的計算結(jié)果外推得到。

3D-VVER1000全堆芯計算耗時446.81 s，TriVNM計算的keff為1.011 240 9，3D-VVER1000基準題給出的keff參考值為1.011 350，與參考值相比，有效增殖因數(shù)的計算偏差僅為-10.788 pcm。圖4為1/6堆芯示意圖及TriVNM計算的歸一化功率分布及與參考值的比較，歸一化功率的最大相對偏差出現(xiàn)在最外圈組件上，約為0.344%。為驗證程序的計算精度，對同類程序的計算結(jié)果進行對比，結(jié)果列于表2?？煽闯觯c其他程序相比，程序TriVNM計算結(jié)果與參考解的偏差較小。

表2 3D-VVER1000基準題在不同程序下的計算結(jié)果Table 2 Calculation results of 3D-VVER-1000 benchmark problem with different programs

圖4 歸一化功率分布計算結(jié)果及與參考解的對比Fig.4 Calculation result of normalized power distribution and comparison with reference result

2.3 帶無限大反射層的兩區(qū)堆芯例題

為檢驗程序在不規(guī)則非結(jié)構(gòu)幾何下的計算準確性，選擇文獻[18]中的例題3進行計算。非結(jié)構(gòu)幾何組件堆芯布置如圖5所示，兩區(qū)堆芯位于無限大的水池中。堆芯分Ⅰ、Ⅱ兩區(qū)排布，堆芯外圍為反射層(即水池)，取其半徑為45 cm來代替無限大反射層。

圖5 非結(jié)構(gòu)幾何組件堆芯布置Fig.5 Layout of irregular geometric assembly problem

基于三角形網(wǎng)格的程序TriVNM對復雜幾何形狀具有良好的適用性，因為三角形網(wǎng)格組合理論上可逼近任意幾何形狀。圖6為網(wǎng)格剖分示意圖，由TriVNM計算的歸一化功率和keff以及其他程序的計算結(jié)果列于表3，其中FEM2D和ABFEM-T均基于有限元方法。相比之下，由TriVNM計算的keff偏差為20.1 pcm，且由TriVNM計算的keff和歸一化功率均與其他程序較為吻合。驗證了基于任意三角形網(wǎng)格的程序TriVNM適用于復雜幾何堆芯。

圖6 非結(jié)構(gòu)幾何組件基準問題的網(wǎng)格剖分示意圖Fig.6 Mesh used in irregular geometric assembly benchmark problem

表3 不規(guī)則幾何組件基準題計算結(jié)果Table 3 Calculation result of irregular geometric assembly benchmark

3 結(jié)論

隨著新概念堆型的不斷提出，其組件設計和堆芯布置不再采用常規(guī)的幾何形狀，傳統(tǒng)的基于常規(guī)幾何的堆芯計算程序已不適用于具有復雜幾何的新型反應堆堆芯計算，因此，本文開展了基于任意三角形網(wǎng)格的多群中子擴散數(shù)值計算方法研究。首先采用ANSYS軟件對計算區(qū)域進行三角形網(wǎng)格剖分，利用坐標變換將任意三角形轉(zhuǎn)換為正三角形；其次采用Galerkin變分技術(shù)建立包含節(jié)塊中子平衡方程的泛函，將參量利用正三角形內(nèi)正交基函數(shù)進行展開；最后利用變分原理，獲得中子通量密度與節(jié)塊邊界上分中子流的相應關(guān)系，并基于源迭代法對其進行求解?；谏鲜隼碚撃Ｐ烷_發(fā)了基于任意三角形網(wǎng)格的變分節(jié)塊法中子擴散計算程序TriVNM，并采用基準題進行了驗證。結(jié)果表明，TriVNM具有較高的計算精度，并對復雜幾何的新型反應堆具有適用性。