三維球床幾何稀疏條數(shù)長(zhǎng)特征線加速方法

郭 建,郭 炯,李 富,*,嚴(yán) 睿,鄒 楊

(1.中國(guó)科學(xué)院 上海應(yīng)用物理研究所,上海 201800;2.清華大學(xué) 核能與新能源技術(shù)研究院先進(jìn)核能技術(shù)協(xié)同創(chuàng)新中心 先進(jìn)反應(yīng)堆工程與安全教育部重點(diǎn)實(shí)驗(yàn)室,北京 100084)

特征線法是實(shí)現(xiàn)球床高溫氣冷堆復(fù)雜幾何全堆芯中子輸運(yùn)最有前景的方法[1]。在數(shù)學(xué)上,特征線法是一種廣泛應(yīng)用的偏微分方程解析解法,20世紀(jì)70年代,Askew[2]將偏微分方程特征線法和中子輸運(yùn)方程的離散縱標(biāo)法、碰撞概率法結(jié)合起來(lái),發(fā)展出了求解中子輸運(yùn)方程的特征線法。中子輸運(yùn)特征線法在求解區(qū)域按照離散縱標(biāo)法產(chǎn)生的角度方向稠密布置大量的特征線,角通量在特征線上逐段解析求解,并通過(guò)積分得到網(wǎng)格標(biāo)通量。中子輸運(yùn)方程的求解在特征線上完成,因而特征線法對(duì)求解區(qū)域幾何和材料分布沒(méi)有任何要求,非常適合于具有隨機(jī)分布燃料球的高溫氣冷堆中子輸運(yùn)求解。

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

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

本文在長(zhǎng)特征線加速方法中引入稀疏條數(shù)長(zhǎng)特征線方法,可在不降低加速方程精度的前提下極大地縮短計(jì)算時(shí)間,解決三維幾何加速方法計(jì)算量大的問(wèn)題。

1 傳統(tǒng)長(zhǎng)特征線加速方法

長(zhǎng)特征線加速方法將連續(xù)幾個(gè)線段合并為一個(gè)長(zhǎng)特征線段,如圖1所示,實(shí)線段表示長(zhǎng)特征線。根據(jù)特征線法中的角通量,見(jiàn)式(1),可以得到長(zhǎng)特征線上的源,見(jiàn)式(2):

a——網(wǎng)格標(biāo)通量計(jì)算方法示意圖;b——長(zhǎng)特征線源計(jì)算方法示意圖圖1 長(zhǎng)特征線加速方法網(wǎng)格標(biāo)通量和長(zhǎng)特征線源計(jì)算方法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為特征線段編號(hào);l為特征線段長(zhǎng)度;s=p,p+1,…,q為長(zhǎng)特線從起點(diǎn)p到終點(diǎn)q穿過(guò)的所有特征線段。

長(zhǎng)特征線上的源可在特征線掃描過(guò)程中計(jì)算得到,然后可構(gòu)造低階的基于特征線法的加速方程。長(zhǎng)特征線上的角通量為:

(4)

網(wǎng)格標(biāo)通量用穿過(guò)網(wǎng)格的所有長(zhǎng)特征線平均角通量修正為:

(5)

式中:φ為網(wǎng)格標(biāo)通量,上角括號(hào)為迭代次數(shù),0為啟動(dòng)長(zhǎng)特征線加速的初始值,n為跳出加速時(shí)的迭代結(jié)果;I(M,K)為穿過(guò)網(wǎng)格i的所有角度方向的特征線段。

長(zhǎng)特征線源用其穿過(guò)的所有網(wǎng)格源來(lái)修正為:

(6)

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

圖2 長(zhǎng)特征線加速方法計(jì)算流程Fig.2 Flow chart of macro-track transport acceleration method

2 稀疏條數(shù)長(zhǎng)特征線加速方法

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

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

(7)

三維問(wèn)題相比二維問(wèn)題特征線數(shù)量高一個(gè)數(shù)量級(jí),本文提出了稀疏條數(shù)長(zhǎng)特征線加速方法來(lái)降低三維加速方法的計(jì)算量,并引入了條數(shù)稀疏度概念,如圖4所示。三維幾何中,按照每t個(gè)相鄰特征線取一的方式選擇用于求解加速方程的特征線,這里的t即為條數(shù)稀疏度,此時(shí)網(wǎng)格標(biāo)通量的計(jì)算公式為:

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

(8)

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

3 加速方法參數(shù)選取

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

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

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

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

長(zhǎng)特征線長(zhǎng)度對(duì)加速效果的影響如表2所列,計(jì)算過(guò)程過(guò)中條數(shù)稀疏度為3,長(zhǎng)特征線長(zhǎng)度越長(zhǎng)單次掃描所需時(shí)間越短,長(zhǎng)特征線長(zhǎng)度越長(zhǎng)則低階加速方程相比特征線方程近似越多,需要的迭代步越多,總的加速效果由以上兩個(gè)方面因素共同決定,表現(xiàn)為總掃描時(shí)間先縮短后增加。

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

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

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

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

4 球床幾何加速效果

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

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

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

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

5 結(jié)論

為解決球床高溫氣冷堆三維特征線全堆芯輸運(yùn)問(wèn)題的幾何復(fù)雜性和迭代速度慢的問(wèn)題,本文采用長(zhǎng)特征線加速方法進(jìn)行了加速研究,主要結(jié)論如下。

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

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

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