基于任意三角形網格的中子擴散變分節塊法

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

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

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

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

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

1 理論模型

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

Σs,gΦg(r)+Sg

(1)

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

(2)

其中，g′為能群編號。

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

(3)

其中，節塊v的貢獻為：

(4)

χγ=J·nγ

(5)

(6)

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

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

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

任意三角形和正三角形的坐標轉換關系為：

(7)

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

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

(8)

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

(9)

QGQT=I

(10)

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

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

(11)

(12)

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

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

(13)

(14)

于是，中子標通量密度和源項的展開式系數之間滿足以下關系式：

(15)

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

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

(16)

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

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

(17)

(18)

M=[M1M2…Mγ…]

(19)

(20)

分別令關于φ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-分別為節塊邊界上的出射流中子密度矩向量和入射流中子密度矩向量；B、C、H和R為響應矩陣。

最后，根據上述推導，得到變分節塊法中的節塊響應關系：

(29)

(Ij-RgП)jg=Bgsg

(30)

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

(31)

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

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

(32)

(33)

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

(34)

(35)

Fgg′=χgνΣf,g′

(36)

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

2 數值結果

通過以上理論模型，本文基于FORTRAN90語言開發了基于任意三角形網格的變分節塊法中子擴散計算程序TriVNM。分別采用2D/3D-IAEA、2D/3D-LRA、2D/3D-VVER440、2D/3D-VVER1000、不規則幾何基準題驗證程序TriVNM正確性。本文只選擇矩形幾何組件的2D-IAEA、六角形幾何組件不帶反射層的3D-VVER1000和非結構幾何組件這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對稱分布，堆芯布置和網格剖分如圖2所示。堆芯按雙區布料方案布置，徑向有20 cm厚的水反射層，堆芯外邊界條件為真空邊界(即入射中子流密度為0)，所有組件均采用等效均勻化的參數。由于此基準問題堆芯布置了控制棒，而且在堆芯和反射層交界面上，熱中子通量梯度很大，因此廣泛用于校核雙群中子擴散方程的數值計算模型及其計算精度。

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

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

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

圖3 歸一化功率分布計算結果及與參考解的對比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旋轉對稱，組件的對邊距為23.6 cm。堆芯軸向高度為200 cm，25根控制棒中的6根提到堆芯中部。堆芯的徑向外邊界反照率為0.125，軸向反照率為0.15。TriVNM計算時，徑向上將每個六角形節塊劃分為8個三角形，軸向上每10 cm劃分為1層，參考解仍是由細網差分程序DIF3D-FD在徑向上將每個六角形組件劃分為384和486個三角形、軸向上每5.0 cm為1層的計算結果外推得到。

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

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

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

2.3 帶無限大反射層的兩區堆芯例題

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

圖5 非結構幾何組件堆芯布置Fig.5 Layout of irregular geometric assembly problem

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

圖6 非結構幾何組件基準問題的網格剖分示意圖Fig.6 Mesh used in irregular geometric assembly benchmark problem

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

3 結論

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