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

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

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

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

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

1 理論模型

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

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

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

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

m=1,2,3,4

(1)

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

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

(2)

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

(3)

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

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

(13)

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

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

圖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)驗(yàn)證

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

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

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

(14)

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

(15)

由表1可看出，ICLDFE求積組在多項(xiàng)式(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對于不同多項(xiàng)式階數(shù)的數(shù)值積分相對誤差均為0.00%。數(shù)據(jù)表明該求積組在局部區(qū)域內(nèi)具有較高的數(shù)值積分精度，隨著求積組階數(shù)的增加，數(shù)值積分精度逐漸提高。圖7表明，平均角度網(wǎng)格步長每減少約1個(gè)量級，5階及5階以下多項(xiàng)式數(shù)值積分的相對誤差均下降約4個(gè)量級，ICLDFE求積組具有4階收斂性。

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

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

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

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

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

2) 結(jié)果分析

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

(16)

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

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

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

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

3 結(jié)論

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