基于非均勻變分節(jié)塊法的pin-by-pin計(jì)算加速算法研究

張滕飛，吳宏春，曹良志，李云召，劉曉晶，熊進(jìn)標(biāo)，柴 翔

(1.上海交通大學(xué) 核科學(xué)與工程學(xué)院，上海 200240；2.西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

針對(duì)核反應(yīng)堆三維全堆芯計(jì)算，傳統(tǒng)的中子學(xué)模擬方法大多采用均勻化粗網(wǎng)的擴(kuò)散理論[1]。相比于均勻化粗網(wǎng)的擴(kuò)散理論，pin-by-pin計(jì)算既能更準(zhǔn)確地描述組件內(nèi)部的非均勻效應(yīng)，又能獲得柵元級(jí)別的高分辨中子學(xué)信息，有助于提高三維全堆芯模擬的計(jì)算精度和分辨率，是國(guó)內(nèi)外反應(yīng)堆物理研究的熱點(diǎn)之一。近年來(lái)，pin-by-pin計(jì)算的相關(guān)均勻化理論已有長(zhǎng)足發(fā)展，如SPH因子方法[2-3]、不連續(xù)因子方法[4]等，但另一方面，pin-by-pin計(jì)算的加速技術(shù)研究依然較少，大多情況下只能憑借并行計(jì)算手段提高計(jì)算速度[5-6]。實(shí)際的核反應(yīng)堆堆芯由數(shù)以萬(wàn)計(jì)的棒束構(gòu)成，采用pin-by-pin計(jì)算方法將產(chǎn)生百萬(wàn)量級(jí)的細(xì)網(wǎng)，對(duì)數(shù)值計(jì)算的收斂性和計(jì)算量提出了挑戰(zhàn)，也制約了pin-by-pin計(jì)算方法的應(yīng)用。研究表明，針對(duì)三維穩(wěn)態(tài)壓水堆堆芯問(wèn)題的單核pin-by-pin計(jì)算，需約24 h的計(jì)算時(shí)間才能完成[5]。

為從算法層面上提高pin-by-pin計(jì)算的效率，基于三維擴(kuò)散的非均勻變分節(jié)塊法，本文分別從內(nèi)迭代和外迭代兩方面入手，應(yīng)用粗網(wǎng)加速算法和粗網(wǎng)有限差分(CMFD)算法，研究加速pin-by-pin計(jì)算的方法。

1 理論模型

1.1 三維擴(kuò)散的非均勻變分節(jié)塊法

針對(duì)某一特定能群，在擴(kuò)散近似下，偶階中子擴(kuò)散方程可寫(xiě)為：

(1)

式中：φ為節(jié)塊內(nèi)部中子通量密度；Ω為方位角向量；Σt為中子總截面；Σr為中子移出截面；q為中子源項(xiàng)。

根據(jù)變分原理，在由若干節(jié)塊組成的整個(gè)求解域上，對(duì)應(yīng)擴(kuò)散方程的泛函可寫(xiě)作各節(jié)塊內(nèi)部及其表面上泛函的疊加貢獻(xiàn)。在擴(kuò)散近似下，各節(jié)塊泛函可寫(xiě)為：

(2)

式中：Γ為外部邊界；FV[φ，j]為各節(jié)塊的泛函；V為節(jié)塊體積；j為中子流；jn為每個(gè)節(jié)塊表面的中子流。

利用x-y方向的有限元形狀函數(shù)和z方向的正交多項(xiàng)式對(duì)節(jié)塊內(nèi)部中子通量密度、節(jié)塊表面中子流分別展開(kāi)，可得到節(jié)塊內(nèi)部的中子標(biāo)通量密度矩求解方程及代表節(jié)塊內(nèi)部和表面響應(yīng)關(guān)系的響應(yīng)矩陣方程：

φ=-C(I-Π)j+Hq

(3)

j=(I-RΠ)-1Bq

(4)

式中：φ為中子標(biāo)通量密度矩；q為中子源項(xiàng)展開(kāi)矩；j為對(duì)應(yīng)節(jié)塊表面的中子流密度矩；I為單位矩陣；Π為轉(zhuǎn)移矩陣，其中包含了問(wèn)題中相鄰節(jié)塊的表面出射流和表面入射流的映射關(guān)系；H、C、B、R為響應(yīng)矩陣，僅與節(jié)塊的幾何和材料分布相關(guān)。一般情況下，采用經(jīng)典的紅-黑迭代[7]格式求解式(4)以得到出、入射中子流密度矩，再代入式(3)得到中子標(biāo)通量密度矩。

由于非均勻變分節(jié)塊法在空間離散過(guò)程中采用了有限元形狀函數(shù)，因此具備構(gòu)造非均勻粗網(wǎng)的能力。圖1示出3種不同的粗網(wǎng)剖分方案。以圖1為例，1個(gè)17×17壓水堆UO2燃料組件在徑向上可有多種粗網(wǎng)剖分方案，即圖1中的3種形式。圖1中，灰色部分為控制棒導(dǎo)向管柵元，白色部分為UO2燃料柵元，黑色邊界代表實(shí)際的柵元(材料)邊界，紅色邊界為非均勻粗網(wǎng)邊界。在各粗網(wǎng)中，所有柵元統(tǒng)一采用4個(gè)直角三角形有限元進(jìn)一步剖分，并利用實(shí)際對(duì)應(yīng)的柵元均勻化截面進(jìn)行材料賦值。特別需要指出，17×17粗網(wǎng)剖分對(duì)應(yīng)了傳統(tǒng)pin-by-pin計(jì)算的均勻細(xì)網(wǎng)剖分方式，各粗網(wǎng)僅包含單一的柵元均勻化截面。

a——17×17；b——5×5；c——3×3圖1 3種粗網(wǎng)剖分方案

1.2 內(nèi)迭代的粗網(wǎng)加速方法

在非均勻變分節(jié)塊法中，利用紅-黑迭代格式求解式(4)的過(guò)程即為內(nèi)迭代過(guò)程。內(nèi)迭代處于整個(gè)求解過(guò)程的最內(nèi)層，它的計(jì)算時(shí)間消耗由中子流展開(kāi)向量j中所包含的未知量的個(gè)數(shù)決定，而未知量的個(gè)數(shù)與總的節(jié)塊表面數(shù)目和表面展開(kāi)階數(shù)有關(guān)。將多個(gè)均勻化柵元合并為1個(gè)粗網(wǎng)能有效降低節(jié)塊數(shù)目和節(jié)塊表面數(shù)目，從而顯著提高內(nèi)迭代的計(jì)算效率。

1.3 外迭代的加速算法

粗網(wǎng)加速方法的使用會(huì)增加式(3)中更新中子通量密度矩的計(jì)算負(fù)荷。然而，一方面，不同于內(nèi)迭代，每次外迭代僅進(jìn)行1次中子通量密度矩的更新，其計(jì)算代價(jià)相對(duì)較小；另一方面，外迭代的加速技術(shù)較為成熟，容易實(shí)現(xiàn)高的加速比。因此，配合相應(yīng)的外迭代加速算法，容易實(shí)現(xiàn)可觀的加速效果。

眾所周知，CMFD算法是最為常用且有效的外迭代加速算法之一，被廣泛應(yīng)用于非均勻計(jì)算方法如特征線方法(MOC)中[8-9]。CMFD算法利用有限差分假設(shè)下的粗網(wǎng)擴(kuò)散計(jì)算結(jié)果來(lái)修正細(xì)網(wǎng)輸運(yùn)計(jì)算的結(jié)果，從而以較低的計(jì)算代價(jià)獲得準(zhǔn)確的計(jì)算結(jié)果。其中為了保證細(xì)網(wǎng)輸運(yùn)計(jì)算及粗網(wǎng)擴(kuò)散計(jì)算之間的等價(jià)關(guān)系，利用細(xì)網(wǎng)輸運(yùn)的計(jì)算結(jié)果作為修正項(xiàng)對(duì)粗網(wǎng)擴(kuò)散計(jì)算進(jìn)行修正。本研究將CMFD算法[10]應(yīng)用于非均勻變分節(jié)塊法中，使每個(gè)非均勻粗網(wǎng)對(duì)應(yīng)1個(gè)CMFD均勻粗網(wǎng)，并建立七對(duì)角CMFD矩陣方程進(jìn)行求解。

矩陣分離(PM)算法[7]是變分節(jié)塊法系列方法的經(jīng)典加速方法之一，但傳統(tǒng)PM算法僅能應(yīng)用于內(nèi)迭代的加速計(jì)算。本研究中重新建立了PM算法，將其同時(shí)擴(kuò)展到內(nèi)、外迭代中，稱之為廣義矩陣分離(GPM)算法[11]。PM算法和GPM算法的迭代流程對(duì)比如圖2所示。其中，PM算法僅實(shí)現(xiàn)了對(duì)式(4)中響應(yīng)矩陣方程的矩陣分離，因此只能加速紅-黑迭代過(guò)程；而GPM算法則是將中子標(biāo)通量密度矩求解方程和響應(yīng)矩陣方程同時(shí)進(jìn)行矩陣分離，不僅對(duì)加速紅-黑迭代有效，且能加速裂變?cè)错?xiàng)的收斂速度，從而起到加速外迭代的效果。

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

采用柵元均勻化三維pin-by-pin問(wèn)題對(duì)本研究的計(jì)算精度和計(jì)算速度進(jìn)行驗(yàn)證。問(wèn)題的三維幾何與C5G7基準(zhǔn)題[12]的不插棒布置相同，軸向等距離地劃分為24層，截面為DRAGON產(chǎn)生的7群柵元均勻化截面[13]。在粗網(wǎng)加速計(jì)算中，采用圖1所示的3種粗網(wǎng)劃分方式進(jìn)行計(jì)算。值得注意的是，非均勻粗網(wǎng)節(jié)塊的表面中子流空間變化相比均勻細(xì)網(wǎng)計(jì)算更為劇烈，必須采用更高的徑向表面展開(kāi)階數(shù)。經(jīng)過(guò)對(duì)表面展開(kāi)階數(shù)的敏感性分析，本文針對(duì)17×17細(xì)網(wǎng)計(jì)算采用2階表面展開(kāi)，5×5粗網(wǎng)計(jì)算采用3階表面展開(kāi)，3×3粗網(wǎng)計(jì)算采用5階表面展開(kāi)，軸向統(tǒng)一采用2階展開(kāi)。以上離散階數(shù)能保證3種粗網(wǎng)劃分方式下的keff計(jì)算結(jié)果偏差低于10 pcm，柵元中子通量密度分布最大相對(duì)偏差不超過(guò)0.5%。此外，計(jì)算結(jié)果顯示，CMFD和GPM算法并不會(huì)對(duì)計(jì)算結(jié)果產(chǎn)生不可忽略的偏差，故在本文中僅關(guān)注計(jì)算效率的對(duì)比。

圖2 PM算法(a)和GPM算法(b)的迭代流程對(duì)比

2.1 整體計(jì)算效率對(duì)比

pin-by-pin問(wèn)題計(jì)算效率的對(duì)比列于表1。其中，CMFD/GPM計(jì)算時(shí)間代表求解CMFD或GPM算法所需的計(jì)算時(shí)間，整體求解時(shí)間代表求解整個(gè)三維問(wèn)題所耗費(fèi)的計(jì)算時(shí)間。由表1可見(jiàn)：不采用任何外迭代加速算法時(shí)，采用5×5和3×3粗網(wǎng)加速分別可將整體求解時(shí)間降低至細(xì)網(wǎng)計(jì)算的27%和34%；采用CMFD和GPM算法皆可有效地降低外迭代的次數(shù)。以17×17為例，采用CMFD和GPM算法分別將外迭代次數(shù)降低了88%和86%。然而，由于CMFD算法中需要直接求解耦合了三維粗網(wǎng)節(jié)塊通量的七對(duì)角矩陣，在問(wèn)題中節(jié)塊數(shù)目過(guò)多時(shí)會(huì)嚴(yán)重影響計(jì)算效率：在17×17(pin-by-pin細(xì)網(wǎng))的劃分下，求解CMFD算法的計(jì)算時(shí)間甚至高達(dá)直接細(xì)網(wǎng)計(jì)算的12.70倍。此外，一方面，網(wǎng)格尺寸增大時(shí)，網(wǎng)格總數(shù)量降低，因此求解紅-黑迭代(內(nèi)迭代)的時(shí)間會(huì)有所減少；另一方面，網(wǎng)格尺寸增大時(shí)，外迭代收斂效果變差，外迭代次數(shù)增加。二者的平衡導(dǎo)致出現(xiàn)了5×5粗網(wǎng)加速效果最優(yōu)的現(xiàn)象。綜上所述，雖然GPM算法的外迭代加速效果稍差于CMFD算法，但其更適用于大規(guī)模問(wèn)題的加速計(jì)算，且研究表明GPM算法的穩(wěn)定性更優(yōu)。當(dāng)采用粗網(wǎng)加速時(shí)，CMFD和GPM算法皆可獲得理想的加速比，最高可達(dá)10倍以上。

注：1) 計(jì)算時(shí)間按照17×17細(xì)網(wǎng)計(jì)算時(shí)間為1.00歸一

2.2 外迭代的收斂趨勢(shì)對(duì)比

不同算法外迭代收斂趨勢(shì)的對(duì)比如圖3所示。每次內(nèi)迭代的時(shí)間體現(xiàn)為兩次外迭代之間的時(shí)間差。由圖3可見(jiàn)：粗網(wǎng)加速算法并不會(huì)對(duì)外迭代的總次數(shù)產(chǎn)生明顯變化，但能有效減少每次內(nèi)迭代所耗費(fèi)的計(jì)算時(shí)間，從而加速整體的計(jì)算過(guò)程；兩種外迭代加速算法均對(duì)減少外迭代總次數(shù)起到良好的效果，但顯而易見(jiàn)，CMFD和GPM算法的求解過(guò)程將增加單次外迭代的時(shí)間。整體上，外迭代次數(shù)的減少仍占據(jù)主要影響，且可預(yù)見(jiàn)，當(dāng)采用輸運(yùn)計(jì)算時(shí)，CMFD和GPM算法的求解時(shí)間占總計(jì)算時(shí)間的份額將進(jìn)一步減小，從而實(shí)現(xiàn)更大的加速比。當(dāng)采用GPM算法的5×5粗網(wǎng)和CMFD算法的3×3粗網(wǎng)計(jì)算時(shí)，加速效果最優(yōu)。

3 結(jié)論與展望

基于三維中子擴(kuò)散的非均勻變分節(jié)塊法，本文提出一種加速pin-by-pin計(jì)算的思路，分別從內(nèi)迭代和外迭代兩個(gè)層面加速計(jì)算。利用非均勻幾何的描述能力建立非均勻粗網(wǎng)剖分以加速內(nèi)迭代，并應(yīng)用了CMFD和GPM兩種算法加速外迭代。數(shù)值結(jié)果表明，該方法能在保證精度的前提下，將pin-by-pin計(jì)算所耗費(fèi)的時(shí)間降低10倍以上。后續(xù)的工作包括低階輸運(yùn)計(jì)算功能及并行計(jì)算功能的實(shí)現(xiàn)。