有限元方法應(yīng)用于一維屏蔽計(jì)算的研究

趙京昌，張斌，陳義學(xué)

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

?

有限元方法應(yīng)用于一維屏蔽計(jì)算的研究

趙京昌，張斌，陳義學(xué)

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

摘要:應(yīng)用確定論方法，對(duì)一維平板及一維球的屏蔽計(jì)算進(jìn)行了研究，采用離散SN－有限元方法開(kāi)發(fā)出一維穩(wěn)態(tài)中子輸運(yùn)計(jì)算程序DONTRAN1D。在處理模型邊界條件時(shí)，提出設(shè)置虛擬點(diǎn)的方法。對(duì)負(fù)通量除采用置零修正外，還使用輸運(yùn)方程進(jìn)行中子平衡計(jì)算。最后應(yīng)用該程序進(jìn)行算例驗(yàn)證，并與國(guó)際上成熟的程序進(jìn)行比較分析，求解精度能得到很好的保證。同時(shí)，程序中還保留了大量功能模塊接口，便于以后繼續(xù)開(kāi)發(fā)之用。

關(guān)鍵詞:一維;屏蔽計(jì)算;有限元;中子輸運(yùn)

0　引言

隨著我國(guó)核電事業(yè)的迅速發(fā)展，公眾對(duì)核電安全性要求的呼聲越來(lái)越高。合理有效的屏蔽設(shè)計(jì)可以確保輻射劑量安全，也能保證核動(dòng)力裝置的正常運(yùn)行。屏蔽計(jì)算的主要任務(wù)是給出輻射粒子的空間精細(xì)分布，因此，在數(shù)值求解中對(duì)空間變量的離散處理至關(guān)重要。20世紀(jì)50年代，有限元方法開(kāi)始在工程問(wèn)題的數(shù)值求解中廣泛應(yīng)用。20世紀(jì)70年代初，有限元方法開(kāi)始應(yīng)用在反應(yīng)堆數(shù)值計(jì)算領(lǐng)域［1］，利用有限元方法求解一維單能穩(wěn)態(tài)中子輸運(yùn)問(wèn)題，顯示出了這一方法的巨大優(yōu)越性。本文首先簡(jiǎn)要介紹粒子輸運(yùn)的概念，然后給出基于有限元方法的一維平板和一維球的屏蔽計(jì)算迭代解法。在處理反射邊界處的通量時(shí)，在幾何模型外部設(shè)置了虛擬點(diǎn)，以適應(yīng)輸運(yùn)方程的求解。對(duì)負(fù)通量除采用置零修正外，還使用輸運(yùn)方程進(jìn)行中子平衡計(jì)算。通過(guò)算例，將自主開(kāi)發(fā)的計(jì)算程序DONTRAN1D與國(guó)際上先進(jìn)的程序進(jìn)行比較。

1　中子輸運(yùn)方程

1.1輸運(yùn)方程的建立

研究中子輸運(yùn)過(guò)程應(yīng)用的一條基本原理是“中子數(shù)守恒”，即在一定體積內(nèi)，中子總數(shù)對(duì)時(shí)間的變化率應(yīng)該等于該體積內(nèi)中子的產(chǎn)生率減去該體積內(nèi)中子的泄漏率和移出率。

式中: n為中子總數(shù); t為時(shí)間; Q為中子產(chǎn)生率; L為中子泄漏率; R為中子移出率。

本文研究屏蔽計(jì)算方法，只討論平衡狀態(tài)，也就是中子數(shù)變化率為零的情形。此時(shí)，輸運(yùn)平衡方程可寫(xiě)為

方程中的3項(xiàng)分別為泄漏項(xiàng)、移出項(xiàng)和產(chǎn)生項(xiàng)。產(chǎn)生項(xiàng)中的中子來(lái)源有散射源、裂變?cè)春酮?dú)立的外中子源。

1.2定解條件

上文導(dǎo)出的中子輸運(yùn)方程是一個(gè)微分積分方程，它只表示中子守恒規(guī)律的數(shù)學(xué)形式。為了封閉這一物理過(guò)程的數(shù)學(xué)描述，還必須確定初始條件和邊界條件，才能使方程的解有定值。從中子角通量的物理意義知道，在方程所適用的區(qū)域內(nèi)，它應(yīng)當(dāng)是連續(xù)的、有限非負(fù)實(shí)數(shù)。

邊界條件比較復(fù)雜，它依賴于所研究問(wèn)題的性質(zhì)。解中子輸運(yùn)方程常遇到的邊界條件如下。

(1)兩種不同介質(zhì)的分界面。若兩種介質(zhì)直接接觸，其間沒(méi)有源和其他物質(zhì)插入，那么根據(jù)連續(xù)性條件，在分界面上應(yīng)滿足=(r，E，Ω)在交界面上沿Ω方向是r的連續(xù)函數(shù)。若交界面中插入第3種介質(zhì)或源，就必須考慮中子穿過(guò)這層介質(zhì)的效應(yīng)，這時(shí)就必須對(duì)上述邊界條件加以修正。

(2)自由表面。自由表面指粒子只能從系統(tǒng)中逃出，而不能再進(jìn)入這一系統(tǒng)。假定中子輸運(yùn)過(guò)程發(fā)生的區(qū)域由凸的且分塊光滑的曲面圍成，中子自表面逸出后就不可能再返回到域中來(lái)。這樣，與真空交界的凸表面就是自由表面。

(3)反射邊界。在反射邊界上，某個(gè)方向的入射中子角通量密度等于與之對(duì)應(yīng)的反射方向的出射中子通量密度。反射邊界條件往往用來(lái)表示解的對(duì)稱性質(zhì)，例如，在系統(tǒng)的對(duì)稱點(diǎn)、對(duì)稱軸或?qū)ΨQ面上經(jīng)常使用反射邊界條件來(lái)減少計(jì)算量。

(4)反照邊界。反照邊界與反射邊界類(lèi)似，但入射角通量密度或中子流與出射角通量密度或中子流不相等，兩者的比值為一個(gè)小于1的常數(shù)α，通常將α稱為反照系數(shù)。

2　有限元方法應(yīng)用于一維屏蔽計(jì)算

輸運(yùn)方程中，中子通量密度是由空間位置、能量和角度共同約束的。進(jìn)行數(shù)值計(jì)算時(shí)，需要對(duì)這些變量做相應(yīng)的離散化處理。

2.1角度及能量變量的近似

對(duì)角度變量的處理采用傳統(tǒng)的離散縱標(biāo)SN方法，但僅限于中子通量密度，源項(xiàng)仍設(shè)置為各向同性源。離散縱標(biāo)SN方法就是利用特別選定的一組方向余弦值{μm} (m = 1，2，…，N)把輸運(yùn)方程離散化。用一組離散的角度坐標(biāo)變量Ωm(m = 1，2，…，N)來(lái)代替連續(xù)的坐標(biāo)變量Ω，在這些特定的方向上對(duì)輸運(yùn)方程進(jìn)行數(shù)值求解，再用數(shù)值求積公式近似表示函數(shù)對(duì)Ω的積分［2］。對(duì)中子角通量密度的積分可表示為

式中: N為總的離散方向數(shù);ωm為每個(gè)方向的求積權(quán)重。為表示方便，這里略去了能量變量下標(biāo)。

一維平板和一維球的輸運(yùn)方程經(jīng)角度離散后變?yōu)?/p>

本文對(duì)能量變量的近似處理采用的是數(shù)值計(jì)算中最常用的分群方法。角通量的能量關(guān)聯(lián)問(wèn)題，與坐標(biāo)系沒(méi)有多大關(guān)系，因此，為不失一般性，以一維平板輸運(yùn)方程為例，它的多群形式可寫(xiě)為

式中:下標(biāo)g為能群編號(hào); Qs和Qf分別為本能群的總散射源和總裂變?cè)础?/p>

2.2空間變量的近似

在這一部分，本文不再使用傳統(tǒng)的菱形差分，而是將線性間斷有限元方法用于空間變量的離散。采用這種非連續(xù)的方法，可得到非常精確和穩(wěn)定的差分組合。

每個(gè)空間網(wǎng)格都有兩個(gè)邊界通量值，但兩個(gè)相鄰網(wǎng)格有一條公共邊界，因此每條公共邊界有兩個(gè)邊界通量值，將這兩個(gè)通量值看作邊界兩旁接近邊界處的通量值。這兩個(gè)值可能不等，即角通量在網(wǎng)格內(nèi)部連續(xù)，但在網(wǎng)格邊界是間斷的［3］。

空間變量離散后，對(duì)輸運(yùn)方程的求解轉(zhuǎn)變?yōu)閷?duì)一系列線性方程組的求解。以?shī)A角余弦為正的情形為例，對(duì)于平板的第k個(gè)網(wǎng)格和第m個(gè)離散方向有

式中:Δx為網(wǎng)格長(zhǎng)度;上標(biāo)b表示該變量來(lái)自相鄰網(wǎng)格。

2.3負(fù)通量修正

在網(wǎng)格邊界處求得的角通量可能為負(fù)，尤其是在無(wú)源粗網(wǎng)格中。由于負(fù)中子通量沒(méi)有物理意義，并且它還會(huì)在很大程度上影響計(jì)算的穩(wěn)定與精確性，應(yīng)予以修正。

本文在置零修正法的基礎(chǔ)上又應(yīng)用了一個(gè)加強(qiáng)條件。對(duì)一個(gè)給定空間網(wǎng)格的求解，首先使用標(biāo)準(zhǔn)的線性間斷有限元方法進(jìn)行計(jì)算。當(dāng)負(fù)通量出現(xiàn)后，立即將其置零，網(wǎng)格內(nèi)其余點(diǎn)的通量重新計(jì)算。置零后，網(wǎng)格內(nèi)有一個(gè)點(diǎn)的通量變?yōu)橐阎?，這時(shí)只需求解一個(gè)方程便能確定另一點(diǎn)唯一的通量值。因?yàn)樗说臋?quán)函數(shù)值為1，所以形式上就是離散化了的輸運(yùn)方程。

3　程序開(kāi)發(fā)與驗(yàn)證

3.1功能及特點(diǎn)

DONTRAN是華北電力大學(xué)核反應(yīng)堆物理與屏蔽研究所正在開(kāi)發(fā)的大型三維輸運(yùn)屏蔽計(jì)算程序，整體采用離散縱標(biāo)SN方法，DONTRAN1D是對(duì)其中的一維問(wèn)題采用更有優(yōu)勢(shì)的有限元方法求解的程序。DONTRAN1D采用確定論方法對(duì)中子輸運(yùn)方程進(jìn)行數(shù)值迭代求解。方向變量的離散使用離散縱標(biāo)法，空間變量離散使用有限元方法，能量變量則用分群近似處理。采用外－內(nèi)迭代方法，在所有幾何、能群、方向上對(duì)離散方程組進(jìn)行迭代求解，獲得中子通量密度的分布。本程序可對(duì)一維平板以及一維球模型的多群中子輸運(yùn)方程求解，主要用于核系統(tǒng)輻射屏蔽計(jì)算，也可用于核裝置的設(shè)計(jì)分析。

3.2測(cè)試算例

本文針對(duì)DONTRAN1D的功能設(shè)計(jì)了多個(gè)算例。屏蔽中最常使用的材料是不銹鋼，設(shè)計(jì)算例時(shí)將其簡(jiǎn)化，只取56Fe一種核素，核子密度采用SS316不銹鋼中鐵的核子密度，溫度為常溫。使用TRANSX程序計(jì)算出反應(yīng)截面數(shù)據(jù)。求積組統(tǒng)一取S16高斯－勒讓德求積組，內(nèi)迭代收斂精度為0.0001，網(wǎng)格均勻劃分。驗(yàn)證單群計(jì)算功能的算例有4個(gè)，分別測(cè)試不同的幾何形狀和固定源分布，計(jì)算參數(shù)見(jiàn)表1。

表1　單群計(jì)算參數(shù)

多群中子的測(cè)試只計(jì)算三群中子的情況，因?yàn)槎嘤谌褐凶拥牡惴ㄅc三群中子的算法是相同的。由于單群測(cè)試時(shí)已經(jīng)驗(yàn)證過(guò)，所以也不再考慮固定源分布的差異。多群計(jì)算參數(shù)見(jiàn)表2。

3.3計(jì)算結(jié)果與分析

將DONTRAN1D的計(jì)算結(jié)果與ANISN的計(jì)算結(jié)果作對(duì)比分析。ANISN是美國(guó)橡樹(shù)嶺國(guó)家實(shí)驗(yàn)室研發(fā)的一維輸運(yùn)程序，是公認(rèn)的精度高和穩(wěn)定性好的計(jì)算程序。對(duì)角度變量的離散同樣使用離散縱標(biāo)法，而對(duì)空間變量的離散則采用傳統(tǒng)的有限差分法。

圖1～圖4為4個(gè)單群中子算例的計(jì)算結(jié)果，從圖中曲線可見(jiàn)，DONTRAN1D與ANISN的計(jì)算結(jié)果吻合度極高。4個(gè)算例所有網(wǎng)格點(diǎn)通量值的最大相對(duì)誤差分別為－0.33％，－0.15％，－9.34％和2.68％，算例3和算例4的最大相對(duì)誤差分別出現(xiàn)在第3個(gè)網(wǎng)格和第1個(gè)網(wǎng)格，在離源較遠(yuǎn)的區(qū)域，相對(duì)誤差均下降到0.10％以下。

表2　多群計(jì)算參數(shù)

圖1　算例1計(jì)算結(jié)果

圖2　算例2計(jì)算結(jié)果

圖3　算例3計(jì)算結(jié)果

圖4　算例4計(jì)算結(jié)果

經(jīng)分析，產(chǎn)生上述問(wèn)題的原因在于，DONTRAN1D處理反射邊界條件時(shí)，在幾何模型外部設(shè)置了一個(gè)等值虛擬點(diǎn)，然后套用輸運(yùn)方程求解。這里只使用了輸運(yùn)方程組中的一個(gè)方程而非整個(gè)方程組，權(quán)函數(shù)的作用沒(méi)有完全體現(xiàn)出來(lái)，由此求出的通量值對(duì)后續(xù)網(wǎng)格的計(jì)算勢(shì)必產(chǎn)生影響。

三群中子算例的計(jì)算結(jié)果如圖5、圖6所示，從圖中曲線可以看出，兩個(gè)程序的計(jì)算結(jié)果吻合度也是非常高的。

圖5　算例5計(jì)算結(jié)果

圖6　算例6計(jì)算結(jié)果

算例5(幾何為平板)的最大相對(duì)誤差出現(xiàn)在第3群第2個(gè)網(wǎng)格，為－0.43％;算例6(幾何為球)的最大相對(duì)誤差出現(xiàn)在第3群第3個(gè)網(wǎng)格，為－9.70％。誤差出現(xiàn)的原因，除上述提到的虛擬點(diǎn)外，還有群間散射源的影響。在計(jì)算第1群散射到第3群以及第2群散射到第3群的散射源時(shí)，要用到第1群和第2群的中子通量密度，所以前兩群計(jì)算結(jié)果的誤差將會(huì)累積到第3群中，使得第3群的相對(duì)誤差最大。

4　結(jié)束語(yǔ)

本文對(duì)有限元方法應(yīng)用于一維屏蔽問(wèn)題的數(shù)值計(jì)算進(jìn)行了初步研究，介紹了自主開(kāi)發(fā)的計(jì)算程序DONTRAN1D，并通過(guò)算例與國(guó)際上的先進(jìn)程序進(jìn)行對(duì)比，計(jì)算精度總體上令人滿意。在處理模型邊界條件時(shí)，提出設(shè)置虛擬點(diǎn)的方法。對(duì)負(fù)通量的處理除采用置零修正外，還再次利用輸運(yùn)方程進(jìn)行中子平衡計(jì)算，效果十分明顯。

參考文獻(xiàn):

［1］MACHORRO E.Discontinuous Galerkin finite element method applied to the 1－D spherical neutron transport equation［J］.Journal of computational physics，2007，223 (1) : 67－81.

［2］HILL T R.ONETRAN: a discrete ordinates finite element code for the solution of the one-dimensional multigroup transport equation［R］.New Mexico: Los Alamos Scientific Laboratory，1975.

［3］杜書(shū)華.輸運(yùn)問(wèn)題的計(jì)算機(jī)模擬［M］.長(zhǎng)沙:湖南科學(xué)技術(shù)出版社，1989.

［4］ADAMS M L.Discontinuous finite-element transport solutions in the thick diffusion limit in Cartesian geometry［R］.Califor-nia: Lawrence Livermore National Laboratory，1991.

［5］KUZMIN D.A guide to numerical methods for transport equa-tions［D］.Freistaat Bayern: Friedrich-Alexander-Universit?t Erlangen－Nürnberg，2010.

(本文責(zé)編:劉芳)

趙京昌(1990—)，男，山東濰坊人，在讀碩士研究生，從事反應(yīng)堆輸運(yùn)屏蔽方面的研究工作(E-mail: zhaojczhaojc@ 163.com)。

作者簡(jiǎn)介:

基金項(xiàng)目:國(guó)家自然科學(xué)基金(11575061)

收稿日期:2015－10－15;修回日期:2015－12－20

中圖分類(lèi)號(hào):TL 328

文獻(xiàn)標(biāo)志碼:A

文章編號(hào):1674－1951(2016)01－0001－04