高 產(chǎn)

(1.中國工程物理研究院 核物理與化學(xué)研究所，四川 綿陽 621900；2.中國工程物理研究院 研究生部，北京 100088)

中子輸運方程是核反應(yīng)過程數(shù)值模擬的基本方程，是一退化雙曲型的積分-微分方程，只在極簡單的情況下才能得到解析解，對于一般實際問題，必須采用一些近似方法求解。中子輸運方程的數(shù)值求解方法[1-3]主要包括離散縱標(biāo)法、球諧函數(shù)法和穿透概率法等。由于求解中子輸運方程的離散縱標(biāo)法實現(xiàn)簡單，具有易處理真空邊界及計算速度較快的優(yōu)點，已被廣泛應(yīng)用于核裝置以及醫(yī)學(xué)等領(lǐng)域中輸運問題的數(shù)值求解。離散縱標(biāo)法對方向變量采用離散縱標(biāo)離散，對空間變量則采用節(jié)塊法、有限差分法和有限元法等離散。有限差分法用于較規(guī)則的幾何問題，不能求解非結(jié)構(gòu)網(wǎng)格問題，三角形節(jié)塊法雖可求解非結(jié)構(gòu)網(wǎng)格問題，但計算精度不高，而有限元方法可采用任意形狀的網(wǎng)格剖分求解區(qū)域，對區(qū)域的形狀有很大的適應(yīng)性，因此，可應(yīng)用于非結(jié)構(gòu)網(wǎng)格問題的求解。

隨著核能及核技術(shù)的發(fā)展，反應(yīng)堆物理計算的發(fā)展趨勢是精確求解三維任意幾何結(jié)構(gòu)的中子輸運問題，避免任何均勻化近似。離散縱標(biāo)-間斷有限元方法(DFEM)因可求解任意幾何結(jié)構(gòu)的中子輸運問題而成為當(dāng)今研究熱點。中子輸運方程的有限元求解方法主要有傳統(tǒng)的Galerkin法、最小二乘有限元法和間斷有限元法。傳統(tǒng)的Galerkin法形成的剛度矩陣非對稱，不利于方程組的快速迭代求解；最小二乘有限元法未對輸運方程做精確處理，難以保證具有較高的計算精度；而間斷有限元法在非均勻介質(zhì)上對輸運方程做了精確處理，因而具有較高的計算精度。美國阿拉莫斯國家實驗室[4-6]研制了ATTILA程序，基準(zhǔn)實驗結(jié)果表明，ATTILA具有很好的精度和計算速度，并能處理幾何復(fù)雜的中子輸運問題。美國德克薩斯農(nóng)工大學(xué)核能與計算科學(xué)系[7]研制了基于非結(jié)構(gòu)網(wǎng)格的并行輸運程序PDT，它對輸運掃描使用全新的迭代方法，并支持多物理應(yīng)用。北京應(yīng)用物理與計算數(shù)學(xué)研究所研制了二維柱幾何Lagrange網(wǎng)格上的非定態(tài)中子光子輸運程序2DSnDFE，并對它進(jìn)行了并行化[8-11]。洪振英等[12]將空間線性間斷有限元法應(yīng)用于一維球幾何動態(tài)粒子輸運方程的求解，數(shù)值算例表明，空間線性間斷有限元法在網(wǎng)格邊界的數(shù)值精度方面明顯高于有限差分法的指數(shù)格式和菱形格式，且通量在時間上的微分曲線相對光滑。巨海濤等[13-14]研制了二維和三維非結(jié)構(gòu)網(wǎng)格下求解中子輸運方程的程序LESFES，數(shù)值計算結(jié)果表明，該方法能用于非結(jié)構(gòu)網(wǎng)格，并具有較高的計算精度。

本文研究基于非結(jié)構(gòu)網(wǎng)格的三維中子輸運問題的離散縱標(biāo)-間斷有限元計算方法，并開發(fā)相應(yīng)的中子輸運計算程序。重點研究每個離散方向的非結(jié)構(gòu)網(wǎng)格的求解順序的方法以及中子輸運DFEM方程中幾個矩陣的矩陣元的精確求解方法，并對一系列基準(zhǔn)例題進(jìn)行計算。

1 中子輸運方程

考慮區(qū)域V、邊界δV，各向同性固定源的多群SN方程[5]為：

Sm,g(r)+Qm,g(r)

(1)

(2)

ψm,g(r)=Fm,g(r)n·Ωm<0

(3)

其中：m、g為角度和能群指標(biāo)；ψ為離散縱標(biāo)角通量；F為入射角通量；σt為宏觀總截面；σs為宏觀散射截面；S為用球諧函數(shù)Yl,n展開表示的散射源；Q為固定源；n為δV的外方向單位法線矢量；Ω為中子運動方向。

(4)

用線性無關(guān)的一組基函數(shù)[γp(r),1≤p≤Pk](p為頂點指標(biāo)；Pk為單元k的頂點數(shù))，定義單元k內(nèi)的近似角通量和標(biāo)通量。這些近似函數(shù)為：

ψ(r)?Ψ(r)=ΘT(r)ψ

(5)

φ(r)?Φ(r)=ΘT(r)φ

(6)

其中：Θ為單元k內(nèi)的基函數(shù)；ψ和φ分別為單元k內(nèi)的角通量和標(biāo)通量。

Θ=[γ1(r),…,γPk(r)]T

(7)

ψ=[ψ1,…,ψPk]T

(8)

φ=[φ1,…,φPk]T

(9)

應(yīng)用Galerkin方法，將式(5)、(6)代入式(4)中，用Θ乘以式(4)，并進(jìn)行Vk的積分。用高斯散度定理，得到：

(10)

其中：ni為δV的外方向單位法線矢量的第i分量；ψs(r)為單元邊界上的角通量。對于每個單元，有：

(11)

其中：l為表面數(shù)；Nf為圍住單元k的表面總數(shù)，如四面體單元的Nf=4，六面體單元的Nf=6。

指定邊界的向上側(cè)面的角通量為：

ψs(r)=ΘT(r)ψs,l

(12)

其中，對于第l表面，有：

(13)

其中：ψinc為共用單元k的第l表面的單元的頂點角通量的相應(yīng)列向量；nav,l為第l表面的平均單位法線矢量，定義為：

(14)

(15)

則每個單元的DFEM方程為：

(16)

其中，

(17)

(18)

(19)

(20)

這些方法的執(zhí)行僅用到四面體單元的線性基函數(shù)。對于四面體網(wǎng)格，對這4個矩陣中的體積積分使用體積坐標(biāo)，面積積分使用面積坐標(biāo)，因而這4個矩陣可解析地精確求解，提高計算的速度和精度。

2 中子輸運方程的求解方法

對于中子輸運方程，遞推計算必須沿中子的運動方向進(jìn)行，這樣，在遞推過程中邊界角通量引進(jìn)的誤差是逐步衰減的，因此，計算格式是穩(wěn)定的。對于中子輸運計算，為加快收斂速度，采用從真空邊界向內(nèi)掃描的方式，本文從第Ⅶ卦限開始計算，計算次序依次是Ⅶ、Ⅷ、Ⅵ、Ⅲ、Ⅴ、Ⅳ、Ⅱ、Ⅰ(表1)。與結(jié)構(gòu)網(wǎng)格求解過程相比，對于某個給定的離散方向Ωm，空間網(wǎng)格的求解順序不再是按幾何位置的順序依次進(jìn)行，而是根據(jù)Ωm方向通過四面體網(wǎng)格邊界的情況決定；空間網(wǎng)格的求解順序?qū)ν回韵薜乃须x散方向不再一樣，不同的離散方向可能有不同的求解順序。

對于某離散方向Ωm，在確定空間網(wǎng)格的求解順序時，需檢查每個網(wǎng)格邊界的中子流向。對于一個平面，它的外法線矢量方向為n=(a×b)/(|a|·|b|)，據(jù)此可得到它的外法線矢量方向的分量nx、ny、nz。通過檢驗量Ωm·n的符號，就能確定離散方向Ωm通過網(wǎng)格邊界的流向。當(dāng)Ωm·n>0時，這個面是出射面，出射面上的角通量未知；當(dāng)Ωm·n<0時，這個面是入射面，入射面上的角通量必須已知，否則無法繼續(xù)求解，通過深度優(yōu)先搜索和廣度優(yōu)先搜索相結(jié)合的方式進(jìn)行網(wǎng)格求解排序，網(wǎng)格求解排序一次建立并儲存起來供以后的掃描使用。

表1 輸運掃描次序

3 基準(zhǔn)校核

基準(zhǔn)問題1[14]是經(jīng)濟(jì)合作與發(fā)展組織/核能源機(jī)構(gòu)委員會的反應(yīng)堆物理的一個小壓水堆(LWR)的基準(zhǔn)問題?？紤]兩種情況：1) Case1，控制棒位置是空的(真空)；2) Case2，插入控制棒。

基準(zhǔn)問題2[14]是經(jīng)濟(jì)合作與發(fā)展組織/核能源機(jī)構(gòu)委員會的反應(yīng)堆物理的一個小快增殖堆(FBR)的基準(zhǔn)問題?？紤]兩種情況：1) Case1，控制棒提出(控制棒位置填充Na)；2) Case2，控制棒插進(jìn)去一半。

基準(zhǔn)問題3[15]是為驗證應(yīng)用有限元方法計算非結(jié)構(gòu)網(wǎng)格輸運問題所構(gòu)造的問題，它由一正方形中間套一圓組成。

散方向采用S4求積組；3) 采用的收斂標(biāo)準(zhǔn)為內(nèi)迭代，Δφ/φ<5×10-4，外迭代，Δkeff/keff<5×10-5。

控制棒價值定義為：

對于基準(zhǔn)問題1，計算了兩種情況的keff及區(qū)域平均中子注量率，SN參考解及開發(fā)的TetTran1.0程序的計算結(jié)果分別列于表2、3。

表2 基準(zhǔn)問題1的keff和控制棒價值

表3 基準(zhǔn)問題1的區(qū)域平均中子注量率

從表2可看到，TetTran1.0的keff與參考解的相對偏差在0.05%以內(nèi)，控制棒價值與參考解的相對偏差在0.65%左右。從表3可看到，對于堆芯區(qū)域的平均中子注量率，Tet-Tran1.0與參考解的相對偏差在0.022%以內(nèi)；對于反射層和控制棒/空隙區(qū)域的平均中子注量率，TetTran1.0與參考解的相對偏差在1.3%以內(nèi)。

對于基準(zhǔn)問題2，計算了兩種情況的keff及區(qū)域平均中子注量率，計算結(jié)果分別列于表4、5。

表4 基準(zhǔn)問題2的keff和控制棒價值

表5 基準(zhǔn)問題2的區(qū)域平均中子注量率

從表4可看到，TetTran1.0的keff與參考解的相對偏差在0.135%以內(nèi)，控制棒價值與參考解的相對偏差在2.6%左右。從表5可看到，對于堆芯區(qū)域的平均中子注量率，Tet-Tran1.0與參考解的相對偏差在0.16%以內(nèi)；對于反射層(軸向包層和徑向包層)和控制棒區(qū)域的平均中子注量率，TetTran1.0與參考解的相對偏差在2.6%以內(nèi)。

對于基準(zhǔn)問題3，計算了keff和每群的區(qū)域平均中子注量率，計算結(jié)果列于表6、7，其中，以MG-MCNP3B的計算結(jié)果作為參考解。

表6 基準(zhǔn)問題3的keff

表7 基準(zhǔn)問題3的區(qū)域平均中子注量率

從表6可看到，TetTran1.0的keff與參考解的相對偏差在0.1%以內(nèi)，而LESFES的keff與參考解的相對偏差在0.34%。從表7可看到，對于燃料區(qū)的平均中子注量率，Tet-Tran1.0與參考解的相對偏差在0.17%以內(nèi)；對于慢化區(qū)的平均中子注量率，TetTran1.0與參考解的相對偏差在1.72%以內(nèi)。

4 結(jié)論

本文使用離散縱標(biāo)-間斷有限元方法求解了三維中子輸運方程，它對能量變量采用多群近似離散，對方向變量采用離散縱標(biāo)法離散，對空間變量采用間斷有限元離散；給出了每個SN離散方向的空間網(wǎng)格的求解排序及中子輸運DFEM方程中幾個矩陣的矩陣元的精確求解方法，并據(jù)此開發(fā)了基于非結(jié)構(gòu)網(wǎng)格的三維輸運計算程序TetTran1.0?；鶞?zhǔn)例題校核表明，該程序具有很高的計算精度。

參考文獻(xiàn)：

[1] LEWIS E E, MILLER W F. Computational methods of neutron transport[M]. US: American Nuclear Society, 1993.

[2] 謝仲生，鄧力. 中子輸運理論數(shù)值計算方法[M]. 西安：西北工業(yè)大學(xué)出版社，2005.

[3] 杜書華，等. 輸運問題的計算機(jī)模擬[M]. 長沙：湖南科技出版社，1989.

[4] WAREING T A, MCGHEE J M, MOREL J E. ATTILA: A three-dimensional, unstructured tetrahedral mesh discrete ordinates transport code[J]. Trans Am Nucl Soc, 1996, 75: 146-147.

[5] WAREING T A, McGHEE J M, MOREL J E, et al. Discontinuous finite element SNmethods on three-dimensional unstructured grids[J]. Nucl Sci Eng, 2001, 138(3): 256-268.

[6] WARSA J S, WAREING T A, MOREL J E. Fully consistent diffusion synthetic acceleration of linear discontinuous SNtransport discretizations on unstructured tetrahedral meshes[J]. Nucl Sci Eng, 2002, 141(3): 236-251.

[7] ASCI project[M/OL]. http:∥parasol.tamu.edu/asci/web_page.

[8] 陽述林. 高維中子輸運方程的離散格式與并行算法研究[D]. 綿陽：中國工程物理研究院，2004.

[9] 莫則堯，傅連祥，陽述林. 非結(jié)構(gòu)網(wǎng)格上求解中子輸運方程的并行流水線SN掃描算法[J]. 計算機(jī)學(xué)報，2004，27(5)：587-595.

MO Zeyao, FU Lianxiang, YANG Shulin. Parallel pipelined SNsweeping algorithm for neutron transport on unstructured grid[J]. Chinese Journal of Computers, 2004, 27(5): 587-595(in Chinese).

[10] MO Zeyao, FU Lianxiang. Parallel flux sweep algorithm for neutron transport on unstructured grid[J]. The Journal of Supercomputing, 2004, 30(1): 5-17.

[11] 魏軍俠，陽述林，王雙虎，等. 非匹配網(wǎng)格上中子輸運方程的間斷有限元方法[J]. 核動力工程，2010，31(S2)：25-28.

WEI Junxia, YANG Shulin, WANG Shuanghu, et al. Discontinuous finite element method for neutron transport equations on no matching mesh[J]. Nuclear Power Engineering, 2010, 31(S2): 25-28(in Chinese).

[12] 洪振英，袁光偉. 粒子輸運方程的線性間斷有限元方法[J]. 計算物理，2009，26(3)：325-334.

HONG Zhenying, YUAN Guangwei. Linear discontinuous finite element method for particle transport equation[J]. Chinese Journal of Computational Physics, 2009, 26(3): 325-334(in Chinese).

[13] 巨海濤，吳宏春，曹良志，等. 二維中子輸運方程的非結(jié)構(gòu)網(wǎng)格離散縱標(biāo)數(shù)值解法[J]. 核動力工程，2006，27(3)：1-5.

JU Haitao, WU Hongchun, CAO Liangzhi, et al. Discrete ordinates method for two-dimensional neutron transport equation based on unstructured-meshes[J]. Nuclear Power Engineering, 2006, 27(3): 1-5(in Chinese).

[14] TAKEDA T, IKEDA H. 3D neutron transport benchmarks, technical report OECD/NEA committee on reactor physics (NEACRP-L-330)[R]. Osaka, Japan: Department of Nuclear Engineering, 1991.

[15] CAO Liangzhi, WU Hongchun. Spheriacal harmonics method for neutron transport equation based on unstructured-meshes[J]. Nuclear Sci Tech, 2004, 15(6): 335-339.