適用于任意幾何的特征線邊界條件處理方法

楊 睿，胡 赟，單浩棟，徐 李

(中國原子能科學(xué)研究院 反應(yīng)堆工程技術(shù)研究部，北京 102413)

特征線方法(MOC)具有強(qiáng)大的幾何處理能力和天然的可并行性，是實(shí)現(xiàn)下一代三維全堆芯高精度、高保真物理計算的優(yōu)秀方案[1]。當(dāng)前，直接三維全堆MOC計算程序主要有MPACT3D[2]、OpenMOC3D[3]和APOLLO3中的TDT[4]等。出于減少存儲和便于處理邊界條件的目的，這些程序均進(jìn)行了一定程度幾何簡化，并沒有真正適用于任意三維幾何。

為向任意幾何拓展，需研究任意幾何下的邊界條件處理問題。當(dāng)前主流的邊界計算方法包括循環(huán)特征線法[5]、插值方法[6-7]、打混方法[7]等。在循環(huán)特征線法中，特征線能首尾相連，處理簡單，不引入邊界誤差，但只能處理規(guī)則幾何，對幾何有很大限制。插值方法理論上可處理任意幾何，但在實(shí)踐中實(shí)現(xiàn)起來較麻煩，另外插值也會引入邊界誤差。打混方法將邊界劃分為若干個面，認(rèn)為每個面上的入射通量是均勻的，計算時將各方向的所有出射通量收集起來進(jìn)行反射，處理相對簡單，但引入誤差最多，不僅會引入空間上的模糊，對于曲面也會引入角度的模糊。上述方法均要求為邊界條件進(jìn)行迭代求解，增加了迭代次數(shù)。插值和打混方法還要求收集邊界通量用于下一次計算，增加了存儲要求。

為更好地解決邊界條件處理問題，本文提出一種將角通量場拆成源場和邊界場分離計算的方法，可同時保留循環(huán)特征線方法的連續(xù)性和插值方法的幾何任意性。

1 特征線方法

離散能群后的穩(wěn)態(tài)中子輸運(yùn)方程可寫為：

(1)

其中：φg(Ω，r)為角通量；Σt，g(r)為總截面；Qg(Ω，r)為總角度源；Ω為立體角；r為空間位置；g為能群標(biāo)號。

反照、反射、周期等邊界條件可寫成統(tǒng)一形式：

φg(r0,Ω)=α(r′)φg(r′out,Ω′)

(2)

其中：r0和Ω為入射位置和角度；r′out和Ω′為出射位置和角度；α(r′)為反照率，滿足0≤α≤1。式(2)表示某處某角度的入射角通量是另一位置和角度下出射角通量乘以反照率。

取r=r0+sΩ，代入式(1)可得：

Qg(r0,Ω,s)

(3)

其中，s為從r0出發(fā)沿Ω方向移動的距離。

截面參數(shù)中省去r0，通過常數(shù)變易法解得特征線方程：

(4)

式(4)中截面與空間有關(guān)，源項與空間和角度有關(guān)，需進(jìn)行相關(guān)處理。

1.1 空間和角度處理

類似于SN方法，取求積組為角度{Ωm}和求積系數(shù){ωm}，每個角度的角通量均滿足式(4)，角度矩的計算使用數(shù)值積分表示；認(rèn)為空間中材料是由若干等材料區(qū)域組成，等材料區(qū)域內(nèi)截面不變；將等材料區(qū)域中源項的空間變量以0階基函數(shù)(平源)展開。經(jīng)上述處理，在每個平源近似區(qū)F中的式(4)變?yōu)椋?/p>

(5)

其中，m為角度的標(biāo)號。

平角度源項可寫為：

(6)

(7)

其中：VF為體積；M為總的離散角度數(shù)。

1.2 空間矩計算和源更新

(8)

圖1 空間矩積分Fig.1 Integration of space moment

大部分特征線方法空間矩中的積分計算使用了控制容積法，本文使用數(shù)值積分法，將式(8)中的面積分離散為若干個點(diǎn)來求解積分，即：

(9)

(10)

代入式(8)和(9)中可得到：

(11)

如果近似認(rèn)為：

(12)

其中，Δsm為射線密度。則式(12)與控制容積法得到的結(jié)果是一致的。數(shù)值積分法可避免討論由于邊界不一致引入的誤差，式(11)的代數(shù)精度為0階，誤差由積分的一階項決定。

如果各方向的射線密度相同，可近似?。?/p>

(13)

其中，LF為各方向的總長度，與角度無關(guān)。

實(shí)際上對于任意三維幾何，由于不能使用循環(huán)特征線，所以任意方向的射線密度可以也應(yīng)當(dāng)取為相同。通過式(11)可對式(6)進(jìn)行更新。在源迭代中，源更新的同時還需進(jìn)行keff的更新，和一般迭代相同，使用：

(14)

式(6)和(14)組成的源更新在每次內(nèi)迭代計算完成后進(jìn)行。

2 邊界分離計算

內(nèi)迭代處理的是邊界條件(式(2))和特征線方程式(式(4))聯(lián)立的方程，式(4)的計算可通過平源近似(式(5))處理，還需要解決的是邊界條件。

2.1 源和邊界分離

將式(4)分解為來自源項和邊界的兩部分：

(15)

其中：

(16)

(17)

(18)

(19)

(20)

其中，n(σ)為σ上方向向外的法向量。

(21)

又0≤α(r(s))≤1，故：

(22)

式(22)和(20)矛盾，故不存在齊次(無源)方程成立的解，所以邊界計算有唯一解，它由源場決定。下面求解這個解。首先考慮源場僅在某處和某角度有出射角通量，其他位置為0：

(23)

圖2 邊界追蹤示意圖Fig.2 Boundary tracking diagram

在追蹤過程中自動滿足了反射邊界式(2)和(17)，而小于最小角通量后的追蹤可近似看作是0解，也成立，所以它是1組可行解，所產(chǎn)生的分布即為式(23)的解。由于反照率小于等于1，且式(17)只會發(fā)生衰減，所以整個邊界追蹤中角通量一直都在縮小，故必然可在有限長度下完成計算。

對于任意源場產(chǎn)生的邊界分布可寫為：

(24)

其中，δ為δ函數(shù)。

容易知道式(17)和(18)的對真空出射通量具有線性疊加性，所以邊界場為每個源場出射通量經(jīng)過邊界追蹤后的結(jié)果的疊加。

圖3 分離算法的解釋Fig.3 Explanation of separation algorithm

這里對分離算法做一解釋，如圖3所示，黃色為計算區(qū)域，上側(cè)和右側(cè)為反射邊界條件。它實(shí)際可等效為1個通過翻轉(zhuǎn)生成的4個原圖形的組合，翻轉(zhuǎn)時的方向如箭頭所示。對于從真空邊界出發(fā)到達(dá)另一真空邊界的綠→藍(lán)→紅射線，它等效于綠虛線→藍(lán)虛線→紅線。特征線方程類似于積分輸運(yùn)方程，都使用了首次碰撞的思想，中子生成后一旦發(fā)生包括散射在內(nèi)的首次碰撞就認(rèn)為消失。這樣綠虛線上生成的中子只會按原方向沿藍(lán)虛線、紅線組成的射線向前運(yùn)動，運(yùn)動時只會發(fā)生衰減，即綠線(源場)中生成中子，在藍(lán)、紅線(邊界場)中衰減。同樣，藍(lán)虛線上也會生成中子，在紅線中衰減。綠虛線、藍(lán)虛線和它們的后續(xù)追蹤，加上紅線部分就構(gòu)成了整條線上的角通量分布，也即分離計算的流程。

在實(shí)際計算時，認(rèn)為源場產(chǎn)生的邊界分布的離散方式與內(nèi)部追蹤的離散方式一致。這樣源場中的1條射線追蹤完成后，可不進(jìn)行下一條源場射線追蹤，而是接著進(jìn)行邊界場中的追蹤。源場和邊界場的射線可實(shí)現(xiàn)完全的首尾相連，不需進(jìn)行插值方法的復(fù)雜處理。同時由于邊界場天然有限，不要求和循環(huán)特征線一樣返回原點(diǎn)，解除了幾何限制。

2.2 空間矩計算

圖4 反射段無法和已有段對齊示意圖Fig.4 Scheme of unmatched reflective segment

為此通過插值方法計算該段的角度和空間求積系數(shù)。先進(jìn)行角度插值，由于絕大多數(shù)求積組的方位角是均勻分布的，故每一極角層中各方位角的權(quán)重一致。假設(shè)Ω的極角余弦為μ，它介于μn和μn+1兩層之間，每層求積系數(shù)分別是ωn和ωn+1，則可?。?/p>

μn<μ<μn+1

(25)

其中，φ為方位角。空間求積系數(shù)的計算為：

(26)

角度的空間求積系數(shù)w(μn,φi)按式(10)計算。若各方向射線密度相同，利用式(13)，則各方向空間求積系數(shù)相同，不需要計算式(26)。

權(quán)重計算完成后，空間矩計算可寫為：

(27)

2.3 迭代計算流程

迭代計算流程如圖5所示。這里雖保留了內(nèi)迭代的說法，但內(nèi)迭代實(shí)際只是角度、真空射線、能群的循環(huán)，不存在迭代。

在實(shí)際實(shí)現(xiàn)時，源場中射線各段信息是按照預(yù)先設(shè)定的射線密度提前追蹤完成。但反射場中的射線由于開始時并不知道要追蹤多長，因此采用了一種“需要時生成并保存”的方式(圖5中虛線框)。當(dāng)射線到達(dá)邊界且出射通量沒有達(dá)到最小通量時，那么檢測反射射線是否已生成，若未生成，生成并保存；若已生成，則沿著先前保存的射線信息計算。另外，由于權(quán)重計算只與角度有關(guān)而與具體追蹤過程無關(guān)，且計算成本很小，可以在發(fā)生反照或反射時進(jìn)行計算。

每條源場射線的初始通量為0，先按式(16)進(jìn)行有源計算，然后通過邊界處理式(2)后開始邊界追蹤。虛線框下的通量過小判斷實(shí)際是在計算式(17)的同時進(jìn)行，對于式(17)中的每一段均會檢測，如果出現(xiàn)通量小于限定值則進(jìn)入下一源場射線追蹤。如果該邊界線追蹤完成后還沒有過小則繼續(xù)追蹤下一條邊界線。因為邊界追蹤是持續(xù)衰減的，所以經(jīng)過有限長度后必然可以小于給出的限定值。

圖5 迭代計算示意圖Fig.5 Scheme of iteration calculation

邊界分離計算是指源場射線和它的反射射線采用有源和無源兩種方法，是計算方法的分離，但幾何和通量傳遞上是連續(xù)的。與循環(huán)特征線方法相比，這種方法相當(dāng)于一種長度有限的循環(huán)特征線方法，特征線首尾相連。但由于邊界自然衰減至0，并不要求特征線需經(jīng)過若干次反射回歸原點(diǎn)，避免了繁瑣的討論及幾何限制。相比于插值方法，該方法相當(dāng)于將邊界上的空間角度插值轉(zhuǎn)移到角度權(quán)重插值，而特征線中的求積組又是固定的，插值簡單且一次性。

此外該方法還有兩個優(yōu)點(diǎn)：1) 單次特征線循環(huán)就可保證邊界條件收斂，不需存儲邊界出射角通量，也不需邊界迭代；2) 每條特征線均是獨(dú)立計算，與其他特征線的角通量無關(guān)，初始和結(jié)束都為0，不需保存或傳遞中間通量，便于實(shí)現(xiàn)特征線并行計算。

3 數(shù)值驗證

3.1 Takeda算例

采用Takeda算例[8]的問題1驗證本文方法的正確性。Takeda算例是1個三維簡化堆芯模型，幾何布置如圖6所示。

圖6 Takeda算例問題1幾何模型Fig.6 Calculation geometry of Takeda problem 1

圖6中x=0、y=0和z=0平面為反射邊界條件，其余為真空邊界條件；能群數(shù)為2。case1中紅色部分填充材料為空塊，case2中紅色部分填充材料為控制棒。

計算時的參數(shù)選擇為：網(wǎng)格尺寸取為0.5 cm，射線密度取為0.05 cm×0.05 cm，求積組選擇勒讓德-切比雪夫求積組，8極角和16方位角。keff計算結(jié)果列于表1，循環(huán)特征線的計算結(jié)果來自MPACT3D[2]。由表1可看出，keff計算結(jié)果誤差很小。比較了3種材料中各能群平均通量的誤差，結(jié)果列于表2。

表1 keff計算結(jié)果Table 1 Calculation result of keff

表2 平均通量誤差的計算結(jié)果Table 2 Calculation result of error of average flux

由表2可看出，大部分區(qū)域計算結(jié)果與參考解的誤差總體很小。然而無論是case1還是case2，空塊和反射層中第2群的計算結(jié)果的誤差均較大。這一問題同樣在MPACT3D中出現(xiàn)，初步分析主要原因可能是低階散射帶來的影響。為排除這一因素，使用相同截面單獨(dú)與MAPCT3D的結(jié)果進(jìn)行對比，結(jié)果列于表3。由表3可看出，循環(huán)特征線方法與邊界分離計算的平均通量結(jié)果相差很小，均在0.5%以內(nèi)。

表3 與循環(huán)特征線通量結(jié)果對比Table 3 Flux comparison with cyclic characteristics method

3.2 C5G7算例

本文方法可用于1/8、1/6或1/12對稱的堆芯。以C5G7-2D算例[9]為例，在上、下表面設(shè)置反射邊界條件以構(gòu)成三維算例。根據(jù)對稱性，該問題是1/8對稱，取分界面如圖7a對角線所示，每個組件均會被一分為二，存在半柵元(圖7b)。該分界面被設(shè)置為全反射邊界條件。

分別計算了沒有分界面的1/4堆芯和有分界面的1/8堆芯，柵元內(nèi)的網(wǎng)格劃分為兩環(huán)和8個扇區(qū)，計算結(jié)果列于表4。由表4可見，1/4和1/8堆芯的計算結(jié)果相差很小，但1/8堆芯可節(jié)約48.27%的時間。此類問題需將組件和柵元進(jìn)行切割，在常規(guī)確定論計算中較少見，借助于邊界分離算法可進(jìn)行求解。

圖7 1/8的C5G7問題Fig.7 Eighth C5G7 problem

表4 C5G7問題計算結(jié)果Table 4 Calculation result of C5G7 problem

3.3 單鈾球水腔模型

為進(jìn)一步驗證邊界的任意性，采用文獻(xiàn)[10]中單鈾球水腔模型進(jìn)行驗證。該問題的最外層為球形反射邊界條件，內(nèi)部為輕水，外部為燃料U-10Zr，模型構(gòu)型如圖8所示。

因文獻(xiàn)[10]中未給出具體溫度，這里溫度取為297 K，U-10Zr的密度取為15.92 g/cm3。利用Scale程序進(jìn)行并群，表5列出MOC程序和Scale程序的計算結(jié)果。

計算時采用了0階散射，對于這種各向異性較強(qiáng)的問題，從計算結(jié)果看誤差相對于前兩個算例較大，但仍在可接受范圍內(nèi)。

圖8 單鈾球水腔模型構(gòu)型Fig.8 Profile of single uranium sphere model with water cavity

表5 單鈾球水腔模型的keff計算結(jié)果Table 5 keff calculation result of single uranium sphere model with water cavity

4 結(jié)論

本文提出了一種適用于任意幾何的特征線邊界計算方法，將角通量場分離成源場和邊界場處理，基于數(shù)值積分和權(quán)重插值給出了迭代計算流程，兼具了插值方法的幾何適應(yīng)性和循環(huán)特征線的首尾相連性。經(jīng)Takeda算例、單鈾球水腔模型和C5G7算例驗證，與參考解keff的最大誤差分別為21、319和138.8 pcm，表明在多種邊界條件下計算結(jié)果仍然可靠。在該方法下，邊界嚴(yán)格對齊且不需存儲邊界通量，避免了邊界條件的迭代；每條特征線可完全孤立計算，適用于特征線并行計算。后續(xù)正在開發(fā)全堆任意幾何的特征線并行計算程序。