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

1（西安交通大學 西安 710049）

2（中國核動力研究設計院 成都 610041）

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

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

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

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

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

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

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

1 計算方法

1.1 任意三棱柱網格的離散縱標節(jié)塊法

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

多群中子輸運方程的形式如下所示［20］：

其中：

對于式（1）考慮各向同性散射，角度使用SN離散，在一三棱柱內的多群中子輸運方程形式可寫為：

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

中子源項Q?g(x，y，z)包括裂變源和散射源，cm-3·s-1，具體表達式可寫為：

令：

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

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

圖2 直角坐標系下等邊三角形(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。

對于式（4），略去x′和y′的上標，可寫為：

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

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

式中：

對于式（6）和式（7）進行微分方程求解即可得到相應方向的中子通量密度。

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

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

1.2 幾何建模

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

圖3展示了一些非結構堆芯的徑向建模及網格剖分示意圖。

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

1.3 并行算法

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

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

2 數值驗證與應用

2.1 轉鼓反應堆計算驗證

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

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

對于ACMIR的建模方式，燃料組件采用六角形陣列模式，其余部分采用CSG 進行真實建模，SARAX的徑向建模及網格剖分如圖5所示。

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

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

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

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

表 1 空間堆特征值計算結果Table 1 The keff results of ACRIM

2.2 熱管堆計算驗證

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

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

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

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

表 2 熱管堆特征值計算結果Table 2 The keff results

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

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

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

2.3 并行效率測試

在天河二號服務器上對SARAX 的并行算法的強并行效率進行了測試，即問題規(guī)模不變，不斷增加并行核數。對于測試的問題，選擇了網格數為5 617×14 的模型A 和網格數為7 132×9 的模型B 進行測試。計算能群數均為33群，S4P1。模型A和模型B的并行計算結果如表3所示。

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

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

3 結語

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

作者貢獻聲明肖博文負責論文編寫，程序設計開發(fā)；鄭友琦負責論文整體設計，指導寫作；王永平負責論文審閱；喬梁負責并行算法設計開發(fā)；陶昱姍、劉嘯岳負責提供計算數據。