夏文勇 張滕飛 劉曉晶 熊進標 柴 翔

（上海交通大學(xué)核科學(xué)與工程學(xué)院 上海 200240）

核反應(yīng)堆物理計算中，堆芯由不同的組件構(gòu)成，若采用蒙特卡羅程序精細求解，計算時間開銷較大，實際求解中可采用確定論中的擴散或輸運方法。對于強吸收、高泄漏的中子學(xué)問題，往往需要采用較高精度的輸運計算［1］。常用的中子輸運方法包含有限元法［2-3］、離散縱標法［4］和節(jié)塊法［5-6］等。其中，變分節(jié)塊方法［7］為節(jié)塊法的一種，它起源于20世紀90年代，該方法從二階偶對稱形式的中子輸運方程出發(fā)，利用Galerkin變分和Lagrange乘子法在整個求解域建立一個包含節(jié)塊平衡關(guān)系和邊界條件的泛函，以空間上的正交多項式函數(shù)及角度上的球諧函數(shù)為基函數(shù)，采用Ritz離散方法展開泛函，獲得中子通量密度及節(jié)塊邊界處偏中子流密度的節(jié)塊響應(yīng)矩陣，并迭代求解。

中子輸運計算不可避免地存在計算量的問題，如何在保證輸運方法計算精度的情況下減少計算時間，這極大地影響著輸運方法程序的工程適用性。為提高變分節(jié)塊法在低階角度展開下的計算精度及計算效率，本文基于積分變分節(jié)塊方法［8-9］，以1/6堆芯計算替代全堆芯計算。當提高角度展開階數(shù)時，程序求解時間將迅速增加。通過采用分布式存儲方式MPI（Message Passing Interface）［10］使程序?qū)崿F(xiàn)并行計算功能，以進一步減小計算時間開銷。

1 理論模型

在各向同性散射假設(shè)下，三維中子輸運方程可以寫為：積分變分節(jié)塊方法中，偶中子角通量密度ψ+可離散為［8］：

式中：Σt、Σs分別表示中子總截面及中子散射截面，m-1；q 為中子源項，m-3·s-1；φ 為中子標通量密度，m-2·s-1；x、y、z為標準正交的空間多項式向量；ψ+為偶中子角通量密度；f為空間基函數(shù)；ψ(Ω)為偶中子角通量密度展開矩，此時角度依賴關(guān)系以隱式形式寫在展開矩中。

采用積分形式對輸運方程進行變分并結(jié)合邊界條件推導(dǎo)可得通量以及中子流密度方程（3）、（4）：

式中：V、C、B、R為響應(yīng)矩陣，與各個節(jié)塊內(nèi)部的材料以及空間、角度離散的基函數(shù)等有關(guān)［8］；φ為中子標通量密度的展開矩，φ=∫ψ(Ω)dΩ；q為中子源的展開矩；j±分別為出、入射中子流密度展開矩。由于離散形式的特殊性，V、C、B、R矩陣中包含形如∫dΩA(Ω)-1的矩陣運算過程，無法直接獲得解析解，因此在計算過程中采用角度求積組對其進行數(shù)值求積分。在物理上，式（3）主要表征了節(jié)塊內(nèi)部的中子平衡關(guān)系，式（4）主要表征了節(jié)塊表面的中子流連續(xù)性關(guān)系。

變分節(jié)塊法的求解過程主要包含響應(yīng)矩陣構(gòu)造以及矩陣方程求解兩個方面［11-12］。在響應(yīng)矩陣構(gòu)造過程中，不同典型節(jié)塊（截面或幾何區(qū)別于其他節(jié)塊）以及不同能群構(gòu)成了多套響應(yīng)矩陣的集合（V、C、B、R作為一套響應(yīng)矩陣集合）。不同的響應(yīng)矩陣集合的計算過程之間具有天然的脫耦性，可以由多計算核心并行執(zhí)行。

至于矩陣方程的求解，假設(shè)堆芯共有M層，每層徑向由若干小正六邊形組成。圖1給出全堆芯及1/6堆芯在軸向第一層的排列方式。首先對六角形幾何采用四色棋盤格式進行染色，即每個節(jié)塊都標記為綠、紅、橙、藍中的一種顏色，保證相鄰兩個節(jié)塊采用不同顏色。在所有節(jié)塊內(nèi)邊界，以相鄰節(jié)塊的出、入射流作為耦合條件，而在外表面，認為節(jié)塊自身是獨立的，即隨時可以根據(jù)邊界條件和出射流更新其入射流。

在1/6堆芯計算問題中，圖1靠近全堆芯中心處節(jié)塊（如編號為1、2、3、……）的邊界具有相鄰節(jié)塊和節(jié)塊表面（節(jié)塊及表面編號滿足周期性對稱條件），故實際求解中需假想求解域周圍有一層節(jié)塊包覆。圖2展示了1/6堆芯中4色棋盤掃描及邊界處理，并標明節(jié)塊和節(jié)塊表面編號。由圖2（a）可知，節(jié)塊表面所標記的6種顏色對應(yīng)節(jié)塊的6個表面；圖2（b）中展示了求解域內(nèi)節(jié)塊1的徑向6個面入射流均為節(jié)塊2第一個面的出射流，這是由于在周期性對稱條件下，求解域內(nèi)節(jié)塊2第一個面繞節(jié)塊1的六個面逆時針依次旋轉(zhuǎn)60°。節(jié)塊（編號4、7、11、……）左邊相鄰的假想節(jié)塊（編號2、3、5、……）所對應(yīng)的表面編號為求解域內(nèi)節(jié)塊（編號2、3、5、……）表面編號逆時針旋轉(zhuǎn)60°，軸向各層情形類似，不再贅述。

圖1 全堆芯(a)及1/6堆芯(b)四色棋盤掃描Fig.1 Four-color scanning of full core(a)and 1/6 core(b)

圖2 1/6堆芯四色棋盤掃描及邊界處理(a)初始求解域布局，(b)考慮邊界條件后求解域布局Fig.2 Four-color scanning and boundary processing of the 1/6 core(a)Initial solution domain,(b)Considering the boundary conditions

程序的串行［13］及并行計算策略如圖3所示，串行計算將順序執(zhí)行，而并行計算則通過多進程進行單機或多機并行執(zhí)行。假設(shè)節(jié)塊數(shù)為N，能群數(shù)為G，則總的響應(yīng)矩陣集合數(shù)目為N×G。在并行計算過程中，無論是全堆芯還是1/6堆芯計算，對于響應(yīng)矩陣構(gòu)造部分，若將N×G套響應(yīng)矩陣均分給P個計算核心并行計算，則第p個核心所需計算的響應(yīng)矩陣集合為：

若N×G套響應(yīng)矩陣不能均分給P個核心并行計算（負載不均衡），令NG-·P=n則第p個核心所需計算的響應(yīng)矩陣集合為：

假設(shè)堆芯共包含28種典型節(jié)塊及4群能群，因此共有112套響應(yīng)矩陣集合。若分配給9個計算核心并行計算，則計算核心0～3各構(gòu)造13套響應(yīng)矩陣集合，計算核心4～8構(gòu)造12套響應(yīng)矩陣集合，引起負載不均衡，無法達到理想并行效率，此時理論上能達到的最高并行效率為（112/9）/13=95.7%。

對于矩陣方程求解部分，根據(jù)給定的計算核心數(shù)目，堆芯沿軸向自動劃分相應(yīng)數(shù)量的非重疊子區(qū)域。相同子區(qū)域內(nèi)節(jié)塊出射中子流作為相鄰節(jié)塊的入射中子流，不同子區(qū)域內(nèi)節(jié)塊出射中子流通過MPI通信傳遞給相鄰子區(qū)域的相鄰節(jié)塊作為入射流。

2 數(shù)值結(jié)果

為驗證加速求解策略正確性及其計算時間的優(yōu)勢，本文基于OECD/NEACRP發(fā)布的TAKEDA4基準題［5］進行數(shù)值準確性校核驗證，并比較快速計算策略帶來的時間效益。TAKEDA4是一個實驗快中子反應(yīng)堆（KNK-II），由于篇幅有限，本文僅對比控制棒半插情形，其堆芯徑向及軸向結(jié)構(gòu)如圖4所示。堆芯共包含12個材料區(qū)，四群宏觀截面，徑向每盒六角形組件邊長為7.5 cm，堆芯軸向總高度為190 cm。

由前述積分輸運理論推導(dǎo)、1/6堆芯邊界處理及并行區(qū)域分解，基于Fortran90編制了VITAS（Variational Integral Transport Analysis Solver）程序。VITAS程序可實現(xiàn)全堆芯及1/6堆芯求解，并可對此兩種情形進行并行計算。程序?qū)⒃诓⑿协h(huán)境下進行精度校核驗算及計算時間測試，并行效率定義為：

式中：E為并行效率；T1為單核計算時間；T2為多核計算時間；P為處理器數(shù)目。

針對算例TAKEDA4，節(jié)塊內(nèi)部展開階數(shù)為6，節(jié)塊表面展開階數(shù)為2，分別進行全堆芯（徑向剖分為169個節(jié)塊，軸向18層）、1/6堆芯（徑向剖分為29個節(jié)塊，軸向18層）兩種情形計算，并與由蒙特卡羅程序GMVP提供的有效增殖因子（0.983 9±0.000 4）對比。

表1給出了堆芯在控制棒半插棒情形時，各區(qū)域平均中子通量密度的數(shù)值與參考值間相對誤差。通過對比我們可以知道：兩種不同尺寸的堆芯幾何，在不同PN角度展開階數(shù)下各區(qū)域的中子通量密度的誤差基本一致，證明堆芯簡化計算的正確性。當角度展開階數(shù)為P1時，控制棒區(qū)G1能群處中子通量密度誤差的絕對值最大，最大值為3.3%。若提高角度展開階數(shù)至P3，誤差降至1%。當角度展開階數(shù)為P3時，平均中子通量密度的相對誤差數(shù)值均控制在1.2%以內(nèi)，表明積分方法具有很高的計算精度。

圖3 串行(a)和并行(b)計算策略對比Fig.3 Comparison of serial(a)and parallel(b)computing strategies

圖4 KNK-II實驗快堆基準題堆芯布置圖 (a)徑向分布，(b)軸向分布Fig.4 KNK-II experimental fast reactor benchmark core layout (a)Radial distribution,(b)Axial distribution

圖5 keff與PN展開階關(guān)系Fig.5 The relationship of keffand PNexpansion order

圖5給出了控制棒半插情形下堆芯有效增殖因子隨PN開階數(shù)的變化，從圖5中可知，在各PN角度展開下，兩種幾何尺寸堆芯的有效增殖因子數(shù)值基本一致。PN展開階數(shù)由P1增至P3，有效增殖因子與參考值誤差迅速由441 pcm降至21 pcm，獲得很好的計算精度。

表2是TAKEDA4在控制棒半插情形下，1/6堆芯及全堆芯中積分變分節(jié)塊方法并行計算結(jié)果。其中，響應(yīng)矩陣方程的求解部分采用了分塊響應(yīng)矩陣加速（PM［8］加速）。從可以看到，兩種堆芯計算方式中響應(yīng)矩陣構(gòu)造時間基本一致，這是由于1/6堆芯改變的是堆芯幾何尺寸，不影響典型節(jié)塊數(shù)及能群數(shù)。響應(yīng)矩陣構(gòu)造部分的并行效率低下，一方面是由于并行效率不能達到理想值（3個計算核心時理論并行效率最高為98.25%，9個計算核心時理論并行效率最高為95.73%）；另一方面是由于大維度數(shù)組的通信時間開銷制約了并行計算帶來的時間效益。在響應(yīng)矩陣方程的求解中，由于堆芯尺寸僅為全堆芯的1/6，求解節(jié)塊數(shù)目降低近5/6（522/3 042），故1/6堆芯求解時間開銷較小。隨著并行計算核心數(shù)目增加，總程序求解時間迅速下降，并行效率也隨之降低。這主要是由于計算核心數(shù)目的增加導(dǎo)致通信開銷增加，使得通信時間占比逐漸增大。1/6堆芯幾何中，當并行核心數(shù)目為9時，隨著角度展開階數(shù)增加，并行效率并不單調(diào)：一方面是由于矩陣維度增加，通信時間開銷增加；另一方面是由于分配給各計算核心的計算量增加，計算時間增加。所以并行效率是由計算時間與通信時間共同決定的。同時，并行效率還與并行機此時的數(shù)據(jù)傳輸速度（帶寬）、磁盤I/O等有關(guān)。對于全堆芯幾何與1/6堆芯幾何計算，當并行計算核心數(shù)目為9、角度展開階數(shù)為P5時，程序的總并行效率達到70%左右。

表1 KNK-II基準題半插棒情形時平均中子通量密度誤差（%）Table 1 Error of average neutron flux density of the KNK-II benchmark rod-in-half

表2 TAKEDA4基準題兩種堆芯尺寸下并行計算結(jié)果Table 2 Parallel calculation results of the TAKEDA4 benchmark under two core sizes

3 結(jié)語

本文從變分積分節(jié)塊方法出發(fā)，采用60°周期性對稱邊界條件將全堆芯計算問題簡化為1/6堆芯計算問題。通過比較兩種幾何尺寸下有效增殖因子及各區(qū)的平均中子通量密度，校核了1/6堆芯計算與全堆芯計算結(jié)果的一致性，并通過與GMVP蒙特卡羅程序結(jié)果做對比，驗證了積分變分節(jié)塊方法具有很好的計算精度。同時，初步實現(xiàn)對兩種幾何尺寸堆芯的并行計算。在9個計算核心并行下，相比單核全堆芯計算，1/6堆芯總計算時間開銷在P1、P3、P5下分別約降低1/26、1/25、1/22，初步實現(xiàn)六角形節(jié)塊的快速求解。以上計算并行效率總體偏低，如何提高程序的并行效率還需對并行算法及程序進一步優(yōu)化。