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

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

?

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

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

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

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

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

0　引言

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

1　中子輸運方程

1.1輸運方程的建立

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

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

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

方程中的3項分別為泄漏項、移出項和產(chǎn)生項。產(chǎn)生項中的中子來源有散射源、裂變源和獨立的外中子源。

1.2定解條件

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

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

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

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

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

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

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

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

2.1角度及能量變量的近似

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

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

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

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

式中:下標g為能群編號; Qs和Qf分別為本能群的總散射源和總裂變源。

2.2空間變量的近似

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

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

空間變量離散后，對輸運方程的求解轉(zhuǎn)變?yōu)閷σ幌盗芯€性方程組的求解。以夾角余弦為正的情形為例，對于平板的第k個網(wǎng)格和第m個離散方向有

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

2.3負通量修正

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

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

3　程序開發(fā)與驗證

3.1功能及特點

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

3.2測試算例

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

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

多群中子的測試只計算三群中子的情況，因為多于三群中子的迭代算法與三群中子的算法是相同的。由于單群測試時已經(jīng)驗證過，所以也不再考慮固定源分布的差異。多群計算參數(shù)見表2。

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

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

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

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

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

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

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

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

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

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

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

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

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

4　結(jié)束語

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

參考文獻:

［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］杜書華.輸運問題的計算機模擬［M］.長沙:湖南科學(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)堆輸運屏蔽方面的研究工作(E-mail: zhaojczhaojc@ 163.com)。

作者簡介:

基金項目:國家自然科學(xué)基金(11575061)

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

中圖分類號:TL 328

文獻標志碼:A

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