三維球床幾何稀疏條數長特征線加速方法

郭 建,郭 炯,李 富,*,嚴 睿,鄒 楊

(1.中國科學院 上海應用物理研究所,上海 201800;2.清華大學 核能與新能源技術研究院先進核能技術協同創新中心 先進反應堆工程與安全教育部重點實驗室,北京 100084)

特征線法是實現球床高溫氣冷堆復雜幾何全堆芯中子輸運最有前景的方法[1]。在數學上,特征線法是一種廣泛應用的偏微分方程解析解法,20世紀70年代,Askew[2]將偏微分方程特征線法和中子輸運方程的離散縱標法、碰撞概率法結合起來,發展出了求解中子輸運方程的特征線法。中子輸運特征線法在求解區域按照離散縱標法產生的角度方向稠密布置大量的特征線,角通量在特征線上逐段解析求解,并通過積分得到網格標通量。中子輸運方程的求解在特征線上完成,因而特征線法對求解區域幾何和材料分布沒有任何要求,非常適合于具有隨機分布燃料球的高溫氣冷堆中子輸運求解。

迭代次數多和計算速度慢是限制特征線法應用的主要因素。因此,特征線法通常作為組件計算中的二維輸運計算方法或者全堆芯計算2-D/1-D方法[3]的徑向二維輸運計算方法。WIMS[4]、CASMO[5]和HELIOS[6]等組件計算軟件均采用特征線法作為非均勻幾何中子輸運計算方法。為了解決特征線法計算速度慢的問題,CASMO、OpenMOC[7]和DeCART[8]等程序采用粗網有限差分(CMFD)方法降低迭代次數,可實現具有宏觀方形規則區域的加速,顯著減少了計算時間。基于2-D/1-D方法的MPACT[9]和NECP-X[10]等程序除了采用CMFD方法進行數值加速外,還實現了在集群計算機上的大規模并行,能夠在分鐘量級給出單個算例的結果。

為了突破數值加速方法中宏觀規則幾何的局限,Yamamoto[11]提出了廣義粗網再平衡(GCMR)方法,從理論上指出GCMR方法可應用于非規則幾何的粗網格,然而沒有給出具體可行的方法。柴曉明等[12]提出了廣義粗網有限差分(GCMFD)方法并引入了用于計算非規則幾何耦合系數的寬度因子,但寬度因子無法從理論上直接給出而只能在迭代過程中逐漸逼近,此外,非規則幾何區域網格邊界面積和體積通常只能近似計算,這都減弱了GCMFD方法的加速效果。Smith[13]提出了長特征線加速方法(macro-track transport acceleration method),不需要計算網格界面的中子流,具有和特征線法同樣的幾何處理能力,可真正實現任意幾何特征線法加速。

本文在長特征線加速方法中引入稀疏條數長特征線方法,可在不降低加速方程精度的前提下極大地縮短計算時間,解決三維幾何加速方法計算量大的問題。

1 傳統長特征線加速方法

長特征線加速方法將連續幾個線段合并為一個長特征線段,如圖1所示,實線段表示長特征線。根據特征線法中的角通量,見式(1),可以得到長特征線上的源,見式(2):

a——網格標通量計算方法示意圖;b——長特征線源計算方法示意圖圖1 長特征線加速方法網格標通量和長特征線源計算方法Fig.1 Calculation method of grid scalar flux and macro-track source for macro-track transport acceleration method

Ψi,g,m,k(s0,out)=Ψi,g,m,k(s0,in)e-τi,g,m,k+

(1)

Ψip,g,m,k(sp,in)

(2)

(3)

式中:s為特征線段編號;l為特征線段長度;s=p,p+1,…,q為長特線從起點p到終點q穿過的所有特征線段。

長特征線上的源可在特征線掃描過程中計算得到,然后可構造低階的基于特征線法的加速方程。長特征線上的角通量為:

(4)

網格標通量用穿過網格的所有長特征線平均角通量修正為:

(5)

式中:φ為網格標通量,上角括號為迭代次數,0為啟動長特征線加速的初始值,n為跳出加速時的迭代結果;I(M,K)為穿過網格i的所有角度方向的特征線段。

長特征線源用其穿過的所有網格源來修正為:

(6)

式中,Q為網格內的源。長特征線加速方法的計算流程如圖2所示。Smith[13]基于輕水堆1/4堆芯二維問題研究了長特征線方法的加速效果,在迭代步上得到了20～50倍的加速效果。長特征線加速方法是一種非線性加速方法,本質上是對輸運方程特征法的低階近似,具有良好的穩定性,因此可以用于擴散非線性加速方法發散時問題的加速。長特征線加速方法的缺點是計算量要明顯高于其他加速方法,加速方程求解過于耗時導致其在時間上的加速效果弱于迭代步上的加速效果。

圖2 長特征線加速方法計算流程Fig.2 Flow chart of macro-track transport acceleration method

2 稀疏條數長特征線加速方法

為了解決長特征線加速方法計算量大的問題,二維方法中通常采用稀疏角度的方式。二維特征線法將角度方向劃分為M個方位角和P個極角,總的方向數為MP,方位角按照等角度間隔的方式劃分,極角有多種劃分方式,包括:等角度、等權重和TY求積組[14]等。三維特征線法一般采用旋轉對稱的求積組,常用的求積組有Carlson求積組和Lee求機組等[15]。圖3對比了二維和三維特征線法求積組的不同,可以看出二維方法中求積方向在兩極相對稠密,因此二維特征線法可以選擇具有反射對稱性的少量幾個方位角和1個極角作為加速方程的求解方向,此時的網格標通量計算公式為:

a——二維特征線法求積組,16個方位角,6個極角;b——三維特征線法旋轉對稱求積組,S8圖3 二維和三維特征線法的求積組對比Fig.3 Comparison of quadrature set for 2-D and 3-D method of characteristics

(7)

三維問題相比二維問題特征線數量高一個數量級,本文提出了稀疏條數長特征線加速方法來降低三維加速方法的計算量,并引入了條數稀疏度概念,如圖4所示。三維幾何中,按照每t個相鄰特征線取一的方式選擇用于求解加速方程的特征線,這里的t即為條數稀疏度,此時網格標通量的計算公式為:

a——特征線方法,M=12;b——稀疏角度長特征線加速方法,M=4;c——稀疏條數長特征線加速方法,M=12,t=3圖4 長特征線布置方法Fig.4 Macro-track arrangement method

(8)

采用稀疏條數長特征線加速方法相比稀疏角度的優點是不會降低加速方程在角度離散上的精度。一維離散縱標法可以證明,SN階方程的精度相當于PN-1階,對于二維和三維問題雖然不能得到相同的結論,但求積方向越多方程的精度越高仍然成立。稀疏條數長特征線加速方法不降低特征線方向數,只在空間上減少特征線的數量,由于三維幾何特征線布置相對二維幾何是非常稠密的,減少特征線數量不會明顯影響式(8)特征線數量積分的準確性。

3 加速方法參數選取

條數稀疏度、長特征線長度和長特征線迭代步是影響加速效果的主要因素。本文基于小型輕水堆基準題對長特征線加速方法的加速效果進行了分析。小型輕水堆基準題[16]是一個邊長為50 cm的正方體,計算模型劃分為50×50×50個邊長為1 cm的正方體小網格,三維特征線間距為0.2 cm,長特征線最大長度限制為2 cm,每次啟動加速方法迭代50步后退出。在進行加速參數研究前對特征線間距進行敏感性分析以避免特征線布置的疏密影響加速效果,計算結果如圖5所示,有效增殖因數keff偏差為相對特征線間距為0.05 cm時的keff計算結果,可以看出keff偏差隨著特征線間距降低而趨于穩定,本文特征線間距是根據keff偏差小于5 pcm的原則選取的。

圖5 keff偏差隨特征線間距變化Fig.5 keff difference variation with track spacing

條數稀疏度對加速效果的影響如表1所列,條數稀疏度為1時的總掃描時間(包括特征線掃描和長特征線掃描)增加了接近1倍,簡單計算可以得到長特征線每次掃描所需時間與特征線掃描所需時間相當。條數稀疏度不大于5時,隨著條數稀疏度增加總掃描時間快速降低,迭代步和加速迭代步均不變,條數稀疏度進一步變大時加速效果開始出現明顯減弱,迭代步和加速迭代步也大幅度增加。說明條數稀疏度大于5時,由于長特征線數量減少,式(8)的計算準確度會嚴重偏離真實情況。

表1 長特征線加速方法加速效果隨條數稀疏度的變化Table 1 Macro-track transport acceleration performance with respect to sparsity

長特征線長度對加速效果的影響如表2所列,計算過程過中條數稀疏度為3,長特征線長度越長單次掃描所需時間越短,長特征線長度越長則低階加速方程相比特征線方程近似越多,需要的迭代步越多,總的加速效果由以上兩個方面因素共同決定,表現為總掃描時間先縮短后增加。

表2 長特征線加速方法加速效果隨長特征線長度的變化Table 2 Macro-track transport acceleration performance with respect to length of macro-track

長特征線迭代步用來控制每次啟動加速方法后的迭代次數。長特征線迭代步對加速效果的影響如表3所列,計算過程中條數稀疏度為3,長特征線長度為2.0 cm。長特征線迭代步不大于60步時,迭代步隨長特征線迭代步增加而降低,之后即使進一步增加長特征線迭代步,迭代步也保持不變。總掃描時間的最小值未出現在長特征線迭代步為60步時,而是出現在20步時,原因是總的掃描時間中長特征線掃描占很大一部分,長特征線迭代步的增多抵消了迭代步的降低。

表3 長特征線加速方法加速效果隨長特征線迭代步的變化Table 3 Macro-track transport acceleration performance with respect to macro-track iteration step

綜合以上分析,條數稀疏度取3～5,長特征線長度取2.0 cm左右,長特征線迭代步取20～60步時可以獲得比較好的加速效果。使用最佳參數時,三維長特征線加速方法在計算時間上可以獲得7倍左右的加速效果,此時在迭代步上可以獲得20倍左右的加速效果;只考慮迭代步最多時可以得到25倍的加速效果,與二維方法迭代步加速效果基本一致。

4 球床幾何加速效果

使用球床幾何模型對稀疏條數長特征線加速方法的加速效果進行了驗證,受限于服務器的計算能力,構造了1個由172個面心立方堆積的燃料球的簡化模型,如圖6所示,模型邊長為60 cm,中心有一個邊長為31.46 cm的正方形空腔區域,燃料球半徑為3 cm,球床區域填充率為0.625,三維特征線間距為0.1 cm。燃料球使用小型輕水堆基準題堆芯截面[16],反射層使用小型輕水堆基準題反射層截面,燃料球間隙為空腔,蒙特卡羅程序計算得到的keff=0.843 10(0.00 016)。

a——正方形球床幾何模型;b——球床區域燃料球堆積模型圖6 三維面心立方規則堆積球床計算模型Fig.6 Calculation model of 3-D pebble-bed geometry with face-centered regularly packing pebble

計算結果如表4所列,計算過程中長特征線迭代步固定為50步,重點研究了條數稀疏度和長特征線長度對加速效果的影響。從計算結果可以看出,使用稀疏條數長特征線加速后,迭代步可以獲得最多30多倍的加速效果,特征線追蹤時間可以獲得最多7倍左右的加速效果,此時迭代步加速效果為20倍左右。三維面心立方堆積球床的計算證明稀疏條數長特征線加速方法在包含空腔的非規則幾何中同樣可以獲得比較好的加速效果。

表4 長特征線加速方法對三維球床幾何的加速效果Table 4 Macro-track transport acceleration performance in 3-D pebble-bed geometry

5 結論

為解決球床高溫氣冷堆三維特征線全堆芯輸運問題的幾何復雜性和迭代速度慢的問題,本文采用長特征線加速方法進行了加速研究,主要結論如下。

1) 將長特征線加速方法應用于三維非規則球床幾何并提出了稀疏條數長特征線加速方法,極大地降低了加速方法的計算時間,顯著改善了加速效果。

2) 稀疏條數長特征線加速方法在三維非規則球床幾何掃描時間上可以獲得7倍左右的加速,此時迭代步上可以獲得20倍左右的加速,僅考慮迭代步時最多可以得到30多倍的加速效果,總體來說迭代步加速效果與二維特征線法相近。

3) 通過基準題研究,得到了所測試球床問題加速參數的選取范圍。條數稀疏度取3～5、長特征線長度取2.0 cm左右、長特征線迭代步取20～60步時可以獲得比較好的加速效果。