基于二十面體間斷有限元離散求積組的IRI-TUB基準(zhǔn)題驗證

代 妮，張 斌，陳義學(xué)

(華北電力大學(xué) 核科學(xué)與工程學(xué)院，北京 102206)

粒子輸運計算在核裝置的屏蔽設(shè)計中起到至關(guān)重要的作用，屏蔽計算的核心是粒子輸運方程的求解。穩(wěn)態(tài)輸運方程是一個包含能量、空間和角度共6個自變量的線性微分-積分方程，通常采用數(shù)值離散的方法近似求解[1-2]。離散縱標(biāo)法(SN)作為應(yīng)用最廣泛的確定論輸運計算方法之一，具有計算速度較快、精度較高、適合求解深穿透輸運問題等優(yōu)點[3]。SN對輸運方程的角度變量直接離散，將原來連續(xù)的方向轉(zhuǎn)化為有限個離散方向進行求解。對于各向異性較強的深穿透屏蔽問題，通量密度的角度分布不光滑甚至出現(xiàn)間斷，傳統(tǒng)求積組無法準(zhǔn)確對高階球諧函數(shù)積分，導(dǎo)致離散誤差較大，嚴重影響了屏蔽計算的精度與可靠性。因此，研究具備積分分段間斷函數(shù)能力的求積組是屏蔽計算中的一個重要研究方向。

在SN的求積組研究領(lǐng)域，1964年Carlson和Lathrop[4]提出層對稱求積組(LS)。該求積組按層劃分，層的交點即為求積點，根據(jù)矩條件求解方向余弦與權(quán)重系數(shù)。LS是目前應(yīng)用最廣泛的求積組，但當(dāng)離散階數(shù)大于20時，其權(quán)重系數(shù)會為負。隨后Carlson[5]又提出等權(quán)重求積組(EQN)，人為地使各離散方向的權(quán)重系數(shù)相等以避免權(quán)重為負，但不便于構(gòu)造高階求積組。1987年Walters[6]提出勒讓德-切比雪夫求積組(PNTN)。該求積組具有積分精度高、易于構(gòu)造大量離散方向等優(yōu)勢。1995年Thurgood等[7]提出全對稱TN求積組，采用八面體向球面投影生成求積點，對應(yīng)球面三角形面積為權(quán)重系數(shù)，相比于LS精度更高且保證權(quán)重恒為正。2007年，Endo和Yamamoto[8]研究了偶階奇階求積組(EON)，使用單動量偶階矩和奇階矩求解求積組參數(shù)。該求積組能同時滿足更多矩條件，但仍存在階數(shù)限制，大于16階會出現(xiàn)負權(quán)重。2010年Jarrell等[9]開發(fā)了線性間斷有限元求積組，根據(jù)八面體三角形幾何特性定義離散方向，通過對線性間斷有限元基函數(shù)積分得到相應(yīng)權(quán)重。該求積組可在1/8球面上準(zhǔn)確積分0階和1階球諧函數(shù)，能產(chǎn)生大量離散方向且易于實現(xiàn)局部細化。2012年Ahrens[10]根據(jù)旋轉(zhuǎn)不變性提出了二十面體求積組，該求積組能準(zhǔn)確積分高階球諧函數(shù)，但由于其不具有八面體對稱性而難以處理反射邊界條件。2015年Lau和Adams[11]基于球面四邊形研究了線性和二次間斷有限元方法及相適應(yīng)的局部細化求積組。

本文結(jié)合間斷有限元思想，構(gòu)造二十面體線性和二次間斷有限元求積組，以適應(yīng)角通量密度各向異性較強的屏蔽問題，減小角度離散誤差，提高屏蔽計算方法的計算精度與可靠性。

1 理論模型

1.1 間斷有限元求積組構(gòu)造流程

在離散縱標(biāo)法中，離散方向及其對應(yīng)權(quán)重系數(shù)的集合稱為離散求積組。一般來說，采用離散方向的極角θ及方位角γ的方向余弦(μ,η,ξ)定義坐標(biāo)系。間斷有限元求積組的構(gòu)造流程如圖1所示。首先將正二十面體嵌入單位球內(nèi)，其每個面分別對應(yīng)1個角度區(qū)域。在每個平面三角形內(nèi)建立坐標(biāo)系(u,v)并定義初始離散方向，將平面上的離散點投影到單位球表面，即可得到相應(yīng)且唯一的求積點，如圖2所示。其次根據(jù)初始離散方向求解間斷有限元基函數(shù)，在每個角度區(qū)域內(nèi)對間斷有限元基函數(shù)積分即可得到初始離散方向?qū)?yīng)的求積權(quán)重。隨著離散方向數(shù)的增加，權(quán)重系數(shù)可能出現(xiàn)負值，通過調(diào)整離散方向位置迭代優(yōu)化求積組以保證權(quán)重非負。最后根據(jù)二十面體的旋轉(zhuǎn)與對稱特性產(chǎn)生全角度區(qū)域內(nèi)求積組。

圖1 間斷有限元求積組構(gòu)造流程Fig.1 Construction of discontinuous finite element quadrature sets

圖2 平面三角形到球面三角形的投影Fig.2 Projection from flat triangle to spherical triangle

1.2 初始離散方向與基函數(shù)

a——線性間斷有限元求積組；b——二次間斷有限元求積組圖3 各次間斷有限元求積組Fig.3 Discontinuous finite element quadrature sets of different orders

二十面體映射方法通常選取子三角形區(qū)域的重心作為初始離散方向，如圖3所示。間斷有限元求積組通過細化三角形提高求積組的階數(shù)，對于1階間斷有限元求積組，線性間斷有限元求積組(ICLDFE)每組有4個離散方向，二次間斷有限元求積組(ICQDFE)每組有9個離散方向。N階線性間斷有限元求積組和二次間斷有限元求積組在單個三角形內(nèi)分別含有4N和9N個離散方向。離散方向數(shù)隨求積組的階數(shù)呈指數(shù)增長，能構(gòu)造大量離散方向。

選定初始離散方向后，對各已知離散求積點(μm,ηm,ξm)建立基函數(shù)bm，其中am,n為基函數(shù)系數(shù)，m為初始離散方向編號，n為基函數(shù)的自由度，與間斷有限元基函數(shù)的次數(shù)有關(guān)。線性間斷有限元求積組n=4，基函數(shù)表達式為：

bm(Ω)=am,1+am,2μm+am,3ηm+am,4ξm

m=1,2,3,4

(1)

二次間斷有限元求積組n=9，基函數(shù)表達式為：

bm(Ω)=am,1+am,2μm+am,3ηm+am,4ξm+

(2)

在離散角度區(qū)域內(nèi)，基函數(shù)在其對應(yīng)離散方向上的值為1，在其余方向上值均為0。以二次間斷有限元基函數(shù)建立矩陣方程(3)，通過求逆運算即可求解間斷有限元基函數(shù)中的未知系數(shù)。

(3)

1.3 求積權(quán)重及優(yōu)化策略

對每個離散方向上的基函數(shù)在其相應(yīng)離散角度區(qū)域ΔΩ內(nèi)積分得到求積權(quán)重：

(4)

由于上述積分邊界是曲線，數(shù)值上難以求解。因此，利用雅各比行列式將球面積分(μ,η,ξ)轉(zhuǎn)變?yōu)槠矫娣e分(u,v)。根據(jù)雅各比行列式定義，得：

(5)

其中，|J|為雅各比行列式，即：

(6)

由于坐標(biāo)系選擇的特殊性：

(7)

式(6)可簡化為：

(8)

由二十面體的幾何特性可得：

(9)

(10)

其中，h為平面三角形的高，由二十面體的幾何性質(zhì)決定：

(11)

將式(9)和(10)代入式(8)，并化簡得：

(12)

實際上，隨著求積組階數(shù)的增加，中心三角形的權(quán)重逐漸減小，從而導(dǎo)致權(quán)重出現(xiàn)負值。根據(jù)4個初始離散方向得到的線性基函數(shù)如圖4所示，可看出基函數(shù)在相應(yīng)離散區(qū)域并不完全為正，因此需通過調(diào)整離散方向的位置以保證權(quán)重非負。如圖5所示，規(guī)定每個離散方向僅在其與平面三角形中心點的連線上移動，且不超過其所在子三角形。定義每個離散方向的比例系數(shù)：

(13)

其中：d1為三角形中心到頂點的距離；d2為三角形中心到離散方向點的距離。通過迭代計算調(diào)整比例系數(shù)使中心三角形的權(quán)重等于其對應(yīng)球面三角形的面積。

已知一個三角形上的求積點之后，根據(jù)二十面體的幾何特性通過旋轉(zhuǎn)與對稱即可得到其余19個面的求積點及權(quán)重。圖6為2階線性間斷有限元求積組分布圖，共有320個離散方向均勻分布球表面，最大權(quán)重與最小權(quán)重之比為1.58。另外，二十面體的不同三角形面可采用不同階數(shù)或不同次數(shù)的求積組，易于實現(xiàn)局部角度細化。

圖4 線性間斷有限元基函數(shù)圖像Fig.4 Linear discontinuous finite element basis function

圖5 線性間斷有限元求積組優(yōu)化策略Fig.5 Optimization strategy of linear discontinuous finite element quadrature sets

圖6 2階二十面體線性間斷有限元求積組Fig.6 Second-order discontinuous finite element quadrature sets on icosahedron

2 精度分析及IRI-TUB基準(zhǔn)驗證

2.1 球諧函數(shù)積分精度

為保證粒子平衡，要求離散求積組在全角度區(qū)域內(nèi)盡可能準(zhǔn)確積分高階球諧函數(shù)。但實際屏蔽問題中，局部角度區(qū)域內(nèi)的角通量密度通常是劇烈變化甚至間斷的。為研究二十面體間斷有限元求積組在局部角度區(qū)域內(nèi)的球諧函數(shù)積分精度，用如下公式計算各求積組在1/20

球面上的多項式數(shù)值積分：

(14)

其中，a、b、c為多項式階數(shù)，非負整數(shù)。表1列出1/20球面內(nèi)各求積組數(shù)值積分結(jié)果與數(shù)值解的相對偏差。其中ICLDFE-SN為N階二十面體線性間斷有限元求積組；ICQDFE-SN為N階二十面體二次間斷有限元求積組。圖7示出了分別在不同多項式階數(shù)下ICLDFE求積組的收斂情況，其中縱坐標(biāo)為相對誤差的絕對值，定義橫坐標(biāo)為平均角度網(wǎng)格步長h：

(15)

由表1可看出，ICLDFE求積組在多項式(0,0,0)與(1,0,0)處的相對誤差均為0.00%，ICLDFE求積組能在1/20球面內(nèi)準(zhǔn)確積分0階以及1階的球諧函數(shù)。同理，ICQDFE求積組則可準(zhǔn)確積分2階及2階以下的球諧函數(shù)。ICLDFE-S4及ICLDFE-S5對于不同多項式階數(shù)的數(shù)值積分相對誤差均為0.00%。數(shù)據(jù)表明該求積組在局部區(qū)域內(nèi)具有較高的數(shù)值積分精度，隨著求積組階數(shù)的增加，數(shù)值積分精度逐漸提高。圖7表明，平均角度網(wǎng)格步長每減少約1個量級，5階及5階以下多項式數(shù)值積分的相對誤差均下降約4個量級，ICLDFE求積組具有4階收斂性。

2.2 IRI-TUB基準(zhǔn)題

1) 基準(zhǔn)題描述

IRI-TUB實驗孔道模型為布達佩斯技術(shù)大學(xué)功率為100 kW的研究反應(yīng)堆[12-13]，其目的是研究粒子在孔道內(nèi)的輸運行為及驗證粒子輸運計算方法。由于孔道的存在會使某些區(qū)域的角通量密度各向異性分布趨勢加劇，當(dāng)離散求積組數(shù)值積分角通量密度精度較低時會導(dǎo)致射線效應(yīng)問題[14]，極大地影響了屏蔽計算的精度及輻射屏蔽設(shè)計的可靠性。因此，采用二十面體線性間斷有限元離散求積組進行計算以驗證該求積組的數(shù)值積分精度及對于實際屏蔽模型的適應(yīng)性。

表1 1/20球面內(nèi)各求積組數(shù)值積分相對誤差Table 1 Relative error of numerical integration for quadrature sets in one-twentieth sphere

圖7 不同多項式階數(shù)下數(shù)值積分相對誤差隨平均角度網(wǎng)格步長的變化Fig.7 Relative error of different order polynomials integration versus angular mesh length

圖8 IRI-TUB基準(zhǔn)題幾何結(jié)構(gòu)Fig.8 Geometry structure of IRI-TUB benchmark problem

幾何模型如圖8所示。堆芯由燃料組件、水及石墨組成。堆芯活性區(qū)由24個7.2 cm×7.2 cm的燃料組件組成，高度為50 cm。堆芯外側(cè)配備大型輻射通道，通道內(nèi)由混凝土填充，混凝土塊內(nèi)有一內(nèi)徑為11.8 cm的圓柱形空腔孔道，孔道外側(cè)包圍著厚度為4.5 mm的不銹鋼層，總長度為187 cm。沿孔道中心距離孔道入口分別為0、67、121、148、175 cm處裝有探測器，可測得孔道內(nèi)不同位置中子反應(yīng)率的實驗值。

網(wǎng)格劃分為3 cm×3 cm×3 cm，差分方式選擇菱形差分置零修正方法，多群截面數(shù)據(jù)由處理172群的MUSE1.0截面數(shù)據(jù)庫得到，考慮各向異性P3階勒讓德多項式近似，求積組分別選用具有相近方向數(shù)的PNTN-S60求積組(3 720個方向)與線性間斷有限元局部細化求積組(3 584個方向)，求積組分布如圖9所示。邊界條件均為真空，內(nèi)迭代收斂準(zhǔn)則5×10-4。采用三維離散縱標(biāo)屏蔽程序ARES[15]直接對該模型進行三維輸運計算。

圖9 二十面體線性間斷有限元局部細化求積組Fig.9 Discontinuous finite element local refined quadrature sets on icosahedron

2) 結(jié)果分析

115In(n,n′)115Inm反應(yīng)率的實驗值與計算值比較列于表2。其中C/E(ICLDFE)為采用二十面體間斷有限元局部細化求積組的計算結(jié)果與實驗測量值的比值，C/E(PNTN)為采用60階勒讓德-切比雪夫求積組的計算結(jié)果與實驗測量值的比值，C/E(文獻)為文獻中的計算值與實驗測量值的比值。文獻中的計算值是采用蒙特卡羅程序MORSE-SGC/S接續(xù)二維離散縱標(biāo)程序DOT3.5計算得到。反應(yīng)率計算公式為：

(16)

表2 IRI-TUB基準(zhǔn)題反應(yīng)率的實驗值與計算值比較Table 2 Comparison of measured and calculated reaction rates for IRI-TUB benchmark problem

由表2可看出，采用三維模型直接進行輸運計算得到的結(jié)果均優(yōu)于文獻中給出的計算結(jié)果。對于PNTN-S60求積組，計算結(jié)果與實驗測量值的相對偏差在35%以內(nèi)，且最大偏差位于孔道出口附近。這是由于在孔道內(nèi)部，中子通量密度的各向異性隨距孔道入口距離的增加而增強，PNTN-S60仍不能準(zhǔn)確描述孔道內(nèi)中子通量密度的各向異性分布情況，會嚴重低估中子通量密度。采用二十面體間斷有限元局部細化求積組的計算結(jié)果更接近實驗測量值，與實驗測量值最大相對偏差為25%。且隨著距孔道入口距離的增加，C/E由0.75增加到1.04。二十面體線性間斷有限元局部細化求積組對沿孔道方向的角度區(qū)域進行局部細化，提高局部求積精度。因而孔道內(nèi)的計算值與實驗值偏差較小。

圖10 IRI-TUB基準(zhǔn)題中子能譜的實驗值與計算值比較Fig.10 Comparison of measured and calculated spectra for IRI-TUB benchmark problem

圖10示出了測量點1和測量點4處大于10 eV中子能譜的計算值與實驗值。其中測量點1為孔道入口處，測量點4為距孔道入口148 cm處。從圖10可看出，能量為10 eV以上的中子能譜計算值與實驗測量值呈現(xiàn)相同變化規(guī)律且符合較好。測量點1和4的能量10 eV以上的中子通量密度測量值與實驗值之比分別為0.910和0.867。隨著距離孔道入口的增加，中子通量密度逐漸衰減，從測量點1至4下降約兩個數(shù)量級。

3 結(jié)論

基于間斷有限元思想，構(gòu)造二十面體線性和二次間斷有限元求積組，并對求積方向及權(quán)重進行優(yōu)化。該求積組具有產(chǎn)生大量離散方向且易于局部細化和權(quán)重非負的特點。數(shù)值結(jié)果表明，在1/20球面內(nèi)，二十面體求積組能準(zhǔn)確積分對應(yīng)階數(shù)以下的所有階球諧函數(shù)，且具有4階收斂性；IRI-TUB基準(zhǔn)題中孔道內(nèi)各關(guān)鍵點反應(yīng)率的計算值與實驗值的相對偏差在25%以內(nèi)，中子能譜的計算值與實驗值符合較好。初步驗證了二十面體間斷有限元求積組對于實際深穿透孔道屏蔽問題的適應(yīng)性與可靠性。但由于二十面體的幾何特性，該求積組難以處理反射邊界條件。另外，為保證精度同時提高計算效率，生成局部細化求積組的自適應(yīng)算法有待進一步研究。