夏文勇 張滕飛 劉曉晶 熊進(jìn)標(biāo) 柴 翔

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

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

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

1 理論模型

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

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

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

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

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

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

在1/6堆芯計(jì)算問題中，圖1靠近全堆芯中心處節(jié)塊（如編號(hào)為1、2、3、……）的邊界具有相鄰節(jié)塊和節(jié)塊表面（節(jié)塊及表面編號(hào)滿足周期性對(duì)稱條件），故實(shí)際求解中需假想求解域周圍有一層節(jié)塊包覆。圖2展示了1/6堆芯中4色棋盤掃描及邊界處理，并標(biāo)明節(jié)塊和節(jié)塊表面編號(hào)。由圖2（a）可知，節(jié)塊表面所標(biāo)記的6種顏色對(duì)應(yīng)節(jié)塊的6個(gè)表面；圖2（b）中展示了求解域內(nèi)節(jié)塊1的徑向6個(gè)面入射流均為節(jié)塊2第一個(gè)面的出射流，這是由于在周期性對(duì)稱條件下，求解域內(nèi)節(jié)塊2第一個(gè)面繞節(jié)塊1的六個(gè)面逆時(shí)針依次旋轉(zhuǎn)60°。節(jié)塊（編號(hào)4、7、11、……）左邊相鄰的假想節(jié)塊（編號(hào)2、3、5、……）所對(duì)應(yīng)的表面編號(hào)為求解域內(nèi)節(jié)塊（編號(hào)2、3、5、……）表面編號(hào)逆時(shí)針旋轉(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］及并行計(jì)算策略如圖3所示，串行計(jì)算將順序執(zhí)行，而并行計(jì)算則通過多進(jìn)程進(jìn)行單機(jī)或多機(jī)并行執(zhí)行。假設(shè)節(jié)塊數(shù)為N，能群數(shù)為G，則總的響應(yīng)矩陣集合數(shù)目為N×G。在并行計(jì)算過程中，無(wú)論是全堆芯還是1/6堆芯計(jì)算，對(duì)于響應(yīng)矩陣構(gòu)造部分，若將N×G套響應(yīng)矩陣均分給P個(gè)計(jì)算核心并行計(jì)算，則第p個(gè)核心所需計(jì)算的響應(yīng)矩陣集合為：

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

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

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

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

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

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

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

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

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

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

圖4 KNK-II實(shí)驗(yàn)快堆基準(zhǔn)題堆芯布置圖 (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，獲得很好的計(jì)算精度。

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

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

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

3 結(jié)語(yǔ)

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