SARAX程序在復(fù)雜非結(jié)構(gòu)幾何堆芯計(jì)算中的應(yīng)用研究

肖博文 鄭友琦 王永平 喬 梁 陶昱姍 劉嘯岳

1（西安交通大學(xué) 西安 710049）

2（中國(guó)核動(dòng)力研究設(shè)計(jì)院 成都 610041）

隨著核能的發(fā)展，核反應(yīng)堆的應(yīng)用領(lǐng)域也逐漸多元化。為了使得核反應(yīng)堆能夠滿足在特定場(chǎng)合的應(yīng)用需求，堆芯設(shè)計(jì)也趨于復(fù)雜化，不再僅僅局限于傳統(tǒng)的規(guī)則幾何排布。對(duì)于此類復(fù)雜幾何堆芯的物理計(jì)算，必須開(kāi)展相應(yīng)的新方法研究。目前，針對(duì)復(fù)雜幾何堆芯計(jì)算的方法主要分為兩類，即蒙特卡羅方法和非結(jié)構(gòu)幾何確定論方法，其中蒙特卡羅方法基于統(tǒng)計(jì)學(xué)原理，理論上不受幾何限制，但是由于在進(jìn)行堆芯物理計(jì)算時(shí)花費(fèi)時(shí)間過(guò)長(zhǎng)，此類方法目前仍難以用于工程設(shè)計(jì)計(jì)算，因此國(guó)內(nèi)外針對(duì)反應(yīng)堆設(shè)計(jì)計(jì)算，大多聚焦于非結(jié)構(gòu)幾何的確定論方法。

美國(guó)阿貢國(guó)家實(shí)驗(yàn)室在美國(guó)能源部NEAMS［1］項(xiàng)目的支持下，開(kāi)發(fā)了堆芯物理計(jì)算程序PROTEUS［2］，其 中PROTEUS-PN 和PROTEUSSN［3-4］都是三維非結(jié)構(gòu)幾何計(jì)算程序，二者在空間上都使用連續(xù)有限元離散，幾何適應(yīng)性強(qiáng)。但是由于基于二階形式的輸運(yùn)方程，其難以處理內(nèi)真空以及稀薄介質(zhì)的問(wèn)題。

美國(guó)阿拉莫斯國(guó)家實(shí)驗(yàn)室開(kāi)發(fā)了基于任意四面體或六面體的一階間斷有限元SN程序ATTILA［5-6］，在百核并行時(shí)擁有良好的并行性能。但是，其目前只用于聚變堆包層的計(jì)算，在反應(yīng)堆中應(yīng)用尚面臨困難。

此外，美國(guó)愛(ài)達(dá)荷國(guó)家實(shí)驗(yàn)室開(kāi)發(fā)了非結(jié)構(gòu)幾何輻射輸運(yùn)程序Rattlesnake［7］；法國(guó)原子能委員會(huì)開(kāi)發(fā)了一個(gè)任意六面體的間斷SN 有限元程序SNATCH［8］和MINARET［9］，可分別用于二維任意三角形或三維任意三棱柱幾何的計(jì)算。但由于公開(kāi)發(fā)表資料有限，上述方法和程序的性能與幾何處理能力尚不明確。

近年來(lái)，國(guó)內(nèi)多家單位也在開(kāi)展復(fù)雜幾何的輸運(yùn)算法研究。北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所傅連祥等開(kāi)發(fā)了二維SN 間斷有限元程序NTXY2D［10］，并對(duì)并行算法進(jìn)行了相應(yīng)的研究［11-13］，在百核時(shí)具有較好的并行效率。中國(guó)科學(xué)院計(jì)算機(jī)體系結(jié)構(gòu)國(guó)家重點(diǎn)實(shí)驗(yàn)室閆潔等開(kāi)發(fā)了非結(jié)構(gòu)網(wǎng)格的二維SN程序［14］，國(guó)防科技大學(xué)遲利華等［15］實(shí)現(xiàn)了對(duì)非結(jié)構(gòu)網(wǎng)格的SN 輸運(yùn)算法的多角度多能群同時(shí)并行的計(jì)算。中國(guó)核動(dòng)力院基于二維三角形幾何開(kāi)發(fā)了通過(guò)掃描求解的SN 連續(xù)有限元程序FEGT［16］。但上述方法大多面向輻射場(chǎng)模擬，尚未在反應(yīng)堆設(shè)計(jì)計(jì)算中有所應(yīng)用。

面向新型核反應(yīng)堆設(shè)計(jì)需求，本文采用了基于任意三棱柱網(wǎng)格的SN節(jié)塊法［17］，這種方法能夠適用于毫米至分米級(jí)分辨率的網(wǎng)格離散，可以針對(duì)反應(yīng)堆幾何特點(diǎn)大幅度降低計(jì)算的網(wǎng)格數(shù)，并極大地節(jié)省計(jì)算內(nèi)存，理論上可以兼顧堆芯物理計(jì)算對(duì)精度和效率的要求。

本文在現(xiàn)有的SARAX程序［18-19］系統(tǒng)的基礎(chǔ)上，提出了針對(duì)非結(jié)構(gòu)幾何堆芯特點(diǎn)的靈活建模和空間離散模型，實(shí)現(xiàn)了并行化的任意三棱柱網(wǎng)格SN節(jié)塊方法，并將其應(yīng)用于空間堆、熱管堆等新型復(fù)雜幾何的堆芯物理計(jì)算中，表現(xiàn)出了高的計(jì)算精度和計(jì)算效率。

1 計(jì)算方法

1.1 任意三棱柱網(wǎng)格的離散縱標(biāo)節(jié)塊法

針對(duì)復(fù)雜幾何堆芯，考慮其在徑向上的結(jié)構(gòu)較為復(fù)雜，軸向上不存在彎曲的結(jié)構(gòu)特點(diǎn)，同時(shí)保證計(jì)算的精度以及效率，采用了基于任意三棱柱幾何的離散縱標(biāo)節(jié)塊法［17］。

多群中子輸運(yùn)方程的形式如下所示［20］：

其中：

對(duì)于式（1）考慮各向同性散射，角度使用SN離散，在一三棱柱內(nèi)的多群中子輸運(yùn)方程形式可寫(xiě)為：

式中：m表示SN離散后的某一角度方向；μm、ηm、ξm表示角度方向m在x、y、z坐標(biāo)軸上的分量；g表示能量分群后的某一能群；Ψ m g(x，y，z)表示m方向上第g群的中子角通量密度，cm-2·s-1；hz為三棱柱高度，cm。

中子源項(xiàng)Q?g(x，y，z)包括裂變?cè)春蜕⑸湓?，cm-3·s-1，具體表達(dá)式可寫(xiě)為：

令：

為了使得上述三棱柱幾何的中子輸運(yùn)方程能夠進(jìn)行任意形狀的堆芯計(jì)算，將圖1 所示的任意三角形的坐標(biāo)（x，y）經(jīng)過(guò)有限元的坐標(biāo)變換后，可得如圖2所示計(jì)算坐標(biāo)系(x′，y′)下的中子輸運(yùn)方程：

圖1 直角坐標(biāo)系下任意三角形(x, y)Fig.1 Arbitrary triangular prism in (x, y)

圖2 直角坐標(biāo)系下等邊三角形(x′, y′)Fig.2 Equilateral triangular prism in (x′, y′)

式中：-2/3 ≤x′≤1/3，-ys(x′)≤y′≤ys(x′)，ys(x′)=(x′+ 2/3)/3，-1/2 ≤z≤1/2。

對(duì)于式（4），略去x′和y′的上標(biāo)，可寫(xiě)為：

對(duì)于式（5），在區(qū)間-1/2 ≤z≤1/2 和-ys(x)≤y≤ys(x)上積分，可得x方向的一維橫向積分方程：

然后對(duì)式（6）在-2/3 ≤x≤1/3 和-ys(x)≤y≤ys(x)區(qū)間上積分，可得z方向的一維橫向積分方程：

式中：

對(duì)于式（6）和式（7）進(jìn)行微分方程求解即可得到相應(yīng)方向的中子通量密度。

對(duì)式（6）在區(qū)間-2/3 ≤x≤1/3 上積分，可得三角形節(jié)塊的中子平衡方程：

聯(lián)立求解上述方程，即可得到每個(gè)三棱柱節(jié)塊的中子通量密度。

1.2 幾何建模

根據(jù)上述所采用的堆芯中子輸運(yùn)方法可知，對(duì)于非結(jié)構(gòu)堆芯的計(jì)算，需要將堆芯劃分為三棱柱網(wǎng)格才能夠進(jìn)行輸運(yùn)計(jì)算，因此，對(duì)于堆芯的幾何模型只需要將其在徑向上劃分為三角形網(wǎng)格。采用類似構(gòu)造實(shí)體幾何（Constructive Solid Geometry，CSG）的建模方式，通過(guò)線線組合直接對(duì)非結(jié)構(gòu)幾何的堆芯進(jìn)行真實(shí)建模，然后通過(guò)三角形網(wǎng)格剖分程序Triangle［21］生成三角形網(wǎng)格，從而進(jìn)行全堆的輸運(yùn)計(jì)算。

圖3展示了一些非結(jié)構(gòu)堆芯的徑向建模及網(wǎng)格剖分示意圖。

圖3 非結(jié)構(gòu)幾何堆芯建模網(wǎng)格剖分示意圖(a) 熱管堆，(b) 空間堆Fig.3 Mesh generation of unstructured geometry(a) Heat pipe reactor, (b) Space reactor

1.3 并行算法

由于堆芯計(jì)算采用的是中子輸運(yùn)方程，并且采用的是非結(jié)構(gòu)網(wǎng)格的精細(xì)建模，從而導(dǎo)致堆芯計(jì)算量較大，因此，需要采用并行計(jì)算以減少計(jì)算時(shí)間。輸運(yùn)方程的并行可以從能量、角度和空間三個(gè)維度進(jìn)行分析，首先，能群的并行由于能群數(shù)較少，導(dǎo)致計(jì)算迭代退化嚴(yán)重，并且CPU的數(shù)量受到能群數(shù)的限制；其次，角度上由于數(shù)目的限制以及需要在迭代過(guò)程中進(jìn)行全局的通信操作，因此，能群和角度不適合進(jìn)行并行操作；最后，空間上由于網(wǎng)格數(shù)量較多，能夠適用于分布式并行，分布式并行可以有效地減少內(nèi)存需求，因此本文主要討論在空間上的并行。

空間并行首先需要對(duì)于空間的網(wǎng)格進(jìn)行區(qū)域分解，區(qū)域分解方式采用開(kāi)源的圖形分割程序METIS［22］，該程序能夠?qū)Ψ墙Y(jié)構(gòu)網(wǎng)格進(jìn)行快速而高效地分解，以獲取通信比最小的劃分方式，這樣可以最大程度地降低通信的數(shù)據(jù)量。對(duì)于空間的并行算法，非結(jié)構(gòu)網(wǎng)格具有復(fù)雜的上下游依賴關(guān)系，嚴(yán)重制約了并行算法的效率，根據(jù)掃描算法，可知在不同CPU上進(jìn)行計(jì)算的網(wǎng)格，下游CPU的網(wǎng)格需要等到上游CPU邊界網(wǎng)格計(jì)算完成才能夠開(kāi)始進(jìn)行掃描，隨著CPU數(shù)量的增加，CPU的等待時(shí)間也會(huì)隨之增長(zhǎng)。為了解決這一問(wèn)題，本文采用了塊雅各比的并行算法［23］，其主要依賴于分塊矩陣迭代格式的退化。串行的掃描算法類似于高斯賽德?tīng)柕?，每一個(gè)網(wǎng)格在計(jì)算時(shí)都采用上游網(wǎng)格的最新迭代計(jì)算結(jié)果，該算法將每一個(gè)分塊矩陣交接處的迭代通量改為上一次迭代的通量，從而實(shí)現(xiàn)各個(gè)分塊矩陣的并行計(jì)算。該算法在分塊區(qū)域邊界上將迭代格式退化為雅可比迭代，而內(nèi)部仍然保持串行計(jì)算的格式，因此該方法被稱為塊雅可比算法。該方法的優(yōu)點(diǎn)在于實(shí)現(xiàn)簡(jiǎn)單，且僅需要在每一次迭代完成后進(jìn)行所有區(qū)域邊界面通量的交換，因此通信次數(shù)較少。

2 數(shù)值驗(yàn)證與應(yīng)用

2.1 轉(zhuǎn)鼓反應(yīng)堆計(jì)算驗(yàn)證

轉(zhuǎn)鼓控制反應(yīng)堆選擇一體化斯特林空間堆（Autonomous Circulation Miniature Integrated nuclear Reactor，ACMIR）［24］進(jìn)行計(jì)算。ACMIR為高溫氣冷堆，采用氦氣冷卻，UN 裝料。所有燃料棒被安置在圓柱形壓力容器中，堆芯外環(huán)繞BeO 反射層和6個(gè)控制轉(zhuǎn)鼓。堆芯橫截面如圖4所示。

圖4 ACMIR堆芯徑向示意圖Fig.4 Diagram of the radial core layout of ACMIR

對(duì)于ACMIR的建模方式，燃料組件采用六角形陣列模式，其余部分采用CSG 進(jìn)行真實(shí)建模，SARAX的徑向建模及網(wǎng)格剖分如圖5所示。

圖5 堆芯徑向網(wǎng)格剖分示意圖Fig.5 Diagram of mesh generation of radial core

使用截面產(chǎn)生程序TULIP［25］制作截面，然后將截面提供給堆芯計(jì)算程序進(jìn)行計(jì)算，同時(shí)將截面轉(zhuǎn)化為宏觀截面，能群數(shù)為33群，采用P1近似，使用蒙特卡羅程序NECP-MCX［10］進(jìn)行多群計(jì)算作為參考值，其中蒙特卡羅計(jì)算使用的粒子總數(shù)為1 000 000，計(jì)算總代數(shù)為600代，其中非活躍代數(shù)為200 代。SARAX 計(jì)算采用S4離散角度數(shù)，真空邊界條件，特征值結(jié)果如表1所示。

由表1可知，SARAX的特征值計(jì)算結(jié)果與多群蒙特卡羅的偏差為2.97×10-3。三角形網(wǎng)格的空間離散造成了計(jì)算偏差，而蒙特卡羅計(jì)算采用的是精確幾何模型，沒(méi)有引入網(wǎng)格的近似。之后統(tǒng)計(jì)了燃料組件的功率分布如圖6 所示，第一行為多群蒙特卡羅NECP-MCX 的計(jì)算結(jié)果，第二行為SARAX 的結(jié)果，第三行為兩者之間的相對(duì)偏差。由圖6 可以看出，統(tǒng)計(jì)的組件功率的偏差在-0.33%～0.17%，兩者之間的相對(duì)偏差較小。

圖6 堆芯燃料組件徑向功率分布Fig.6 Radial power distribution of fuel assembly

表 1 空間堆特征值計(jì)算結(jié)果Table 1 The keff results of ACRIM

2.2 熱管堆計(jì)算驗(yàn)證

靜默式鈾基可移動(dòng)反應(yīng)堆［11］是一種兆瓦級(jí)熱管核反應(yīng)堆。堆芯設(shè)計(jì)布局如圖7 所示，堆芯由燃料區(qū)、保溫層、固定式反射層、滑動(dòng)反射層、安全棒、應(yīng)急余排通道以及碳化硼部分組成。燃料區(qū)外設(shè)計(jì)有保溫層，保溫層外布置有大量的鈹反射層，并用碳化硼部分包圍堆芯以減少堆芯與外部環(huán)境之間的相互影響。

圖7 熱管堆布置圖Fig.7 Raidal layout of heat pipe core

堆芯中方形布置的燃料組件以及熱管，使用等效均勻化方法建模，堆芯外圍的控制棒及反射層進(jìn)行真實(shí)建模，徑向分布如圖8所示。

圖8 熱管堆堆芯徑向建模(a)及網(wǎng)格剖分(b)Fig.8 Radial modeling (a) and mesh generation (b) of heat pipe core

表 2 熱管堆特征值計(jì)算結(jié)果Table 2 The keff results

與轉(zhuǎn)鼓反應(yīng)堆計(jì)算類似，采用截面產(chǎn)生程序TULIP 制作截面，然后將截面提供給堆芯輸運(yùn)計(jì)算程序進(jìn)行計(jì)算，同時(shí)將截面轉(zhuǎn)化為宏觀截面，能群數(shù)為33 群，采用P1 近似，使用蒙特卡羅程序NECPMCX 進(jìn)行多群計(jì)算作為參考值，NECP-MCX 計(jì)算粒子數(shù)為1 000 000，計(jì)算總代數(shù)為600代，其中非活躍代數(shù)為200代。SARAX計(jì)算采用S4離散角度數(shù)，真空邊界條件。特征值結(jié)果如表2所示。

由表2可知，SARAX的特征值計(jì)算結(jié)果與多群蒙特卡羅的偏差為1.29×10-3。然后分別統(tǒng)計(jì)了SARAX和NECP-MCX計(jì)算的燃料組件徑向功率分布如圖9所示。由圖9可以看出，統(tǒng)計(jì)的組件功率的偏差在-0.73%～1.35%，兩者之間的相對(duì)偏差較小。

圖9 燃料組件徑向功率分布Fig.9 Radial power distribution of fuel assembly

2.3 并行效率測(cè)試

在天河二號(hào)服務(wù)器上對(duì)SARAX 的并行算法的強(qiáng)并行效率進(jìn)行了測(cè)試，即問(wèn)題規(guī)模不變，不斷增加并行核數(shù)。對(duì)于測(cè)試的問(wèn)題，選擇了網(wǎng)格數(shù)為5 617×14 的模型A 和網(wǎng)格數(shù)為7 132×9 的模型B 進(jìn)行測(cè)試。計(jì)算能群數(shù)均為33群，S4P1。模型A和模型B的并行計(jì)算結(jié)果如表3所示。

表3 模型A和模型B的并行性能Table 3 Parallel performance of model A and model B

從表3 可看出，模型A 和模型B 的并行效率在96 核時(shí)并行效率約為60%，在核數(shù)較多時(shí)并行效率開(kāi)始降低。表明基于塊雅克比并行的三維任意三棱柱節(jié)離散縱標(biāo)節(jié)塊法在百核時(shí)具有較高的并行效率。其在核數(shù)較多（百核以上）時(shí)并行效率降低，是由于塊雅克比的迭代退化從而導(dǎo)致并行效率下降。

3 結(jié)語(yǔ)

本文針對(duì)復(fù)雜非結(jié)構(gòu)幾何堆芯的特點(diǎn)，開(kāi)發(fā)了非結(jié)構(gòu)網(wǎng)格的堆芯物理計(jì)算程序SARAX，實(shí)現(xiàn)了復(fù)雜非結(jié)構(gòu)幾何堆芯的精確建模與網(wǎng)格剖分功能，實(shí)現(xiàn)了區(qū)域分解的塊雅各比并行算法，減少了程序在堆芯設(shè)計(jì)計(jì)算中的計(jì)算時(shí)間，并對(duì)程序在復(fù)雜堆芯的計(jì)算精度進(jìn)行了驗(yàn)證應(yīng)用。在空間堆以及熱管堆的應(yīng)用中，實(shí)現(xiàn)了堆芯的精細(xì)建模，對(duì)特征值與功率分布進(jìn)行了計(jì)算，與蒙特卡羅程序NECP-MCX的多群計(jì)算結(jié)果相比，特征值的計(jì)算偏差小于3.00×10-3，徑向功率分布的相對(duì)偏差小于1.5%，驗(yàn)證了SARAX 在非結(jié)構(gòu)幾何堆芯中計(jì)算中具有較高的精度。并對(duì)SARAX 程序的計(jì)算效率進(jìn)行了測(cè)試，在百核時(shí)的計(jì)算效率能夠達(dá)到60%左右。

作者貢獻(xiàn)聲明肖博文負(fù)責(zé)論文編寫(xiě)，程序設(shè)計(jì)開(kāi)發(fā)；鄭友琦負(fù)責(zé)論文整體設(shè)計(jì)，指導(dǎo)寫(xiě)作；王永平負(fù)責(zé)論文審閱；喬梁負(fù)責(zé)并行算法設(shè)計(jì)開(kāi)發(fā)；陶昱姍、劉嘯岳負(fù)責(zé)提供計(jì)算數(shù)據(jù)。