基于任意三角形網(wǎng)格的中子擴(kuò)散變分節(jié)塊法

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

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

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

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

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

1 理論模型

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

Σs,gΦg(r)+Sg

(1)

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

(2)

其中，g′為能群編號(hào)。

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

(3)

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

(4)

χγ=J·nγ

(5)

(6)

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

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

圖1 任意三角形坐標(biāo)變換示意圖Fig.1 Arbitrary trianglar coordinate transformation

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

(7)

坐標(biāo)變換后進(jìn)行基函數(shù)的構(gòu)造。線(xiàn)性無(wú)關(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)不是正交歸一的。而正定對(duì)稱(chēng)矩陣G一定與單位矩陣相似，即存在滿(mǎn)秩矩陣Q并滿(mǎn)足下式：

(9)

QGQT=I

(10)

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

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

(11)

(12)

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

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

(13)

(14)

于是，中子標(biāo)通量密度和源項(xiàng)的展開(kāi)式系數(shù)之間滿(mǎn)足以下關(guān)系式：

(15)

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

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

(16)

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

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)

為將響應(yīng)矩陣表達(dá)成通用形式，定義邊界上的分中子流密度為：

(24)

Ψ=2(j++j-)

(25)

χ=j+-j-

(26)

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

j+=Bs+Rj-

(27)

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

(28)

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

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

(29)

(Ij-RgП)jg=Bgsg

(30)

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

(31)

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

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

(32)

(33)

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

(34)

(35)

Fgg′=χgνΣf,g′

(36)

所有三角形節(jié)塊通過(guò)中子流相互耦合，然后利用傳統(tǒng)的源迭代法對(duì)其進(jìn)行求解。含特征值的裂變?cè)磫?wèn)題采用冪法迭代求解，稱(chēng)為裂變?cè)吹?。在每次裂變?cè)吹?，需在已知裂變?cè)磃的情況下，進(jìn)行多群迭代。

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

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

2.1 2D-IAEA基準(zhǔn)題

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

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

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

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

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

2.2 3D-VVER1000基準(zhǔn)題

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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