1/4柵元尺度兩級粗網(wǎng)有限差分加速方法

周曉宇， 王浩， 李佩軍， 郝琛,2

(1.哈爾濱工程大學 核安全與仿真技術重點學科實驗室，黑龍江 哈爾濱 150001；2.中國核電工程有限公司， 北京 100840)

粗網(wǎng)有限差分方法(coarse mesh finite difference，CMFD)[1]因其便于實施、加速效果好，得以在全堆芯輸運計算中廣泛研究和應用。但是，用于加速的低階系統(tǒng)求解往往引入額外的計算負擔，尤其瞬態(tài)模擬，CMFD求解時間與特征線方法(method of characteristics,MOC)中子輸運求解時間相當[2]，且存在數(shù)值不穩(wěn)定性的天然缺陷。為解決上述問題，學者們提出了多種更加穩(wěn)定有效的CMFD方法[3-7]。同時為有效降低CMFD自身的計算量，提出針對能量或空間的多級加速方法[8-9]。對于高保真中子輸運計算，無論采用何種類型的CMFD，僅采用基于柵元尺度的CMFD加速，都存在加速效果的瓶頸。以外迭代為MOC為例，無論CMFD怎樣收斂，外迭代次數(shù)會保持在一個特定的數(shù)值[10]。而對于高保真中子輸運計算來說，迭代次數(shù)是決定耗時的關鍵因素，如何打破這個瓶頸，在等價的前提下進一步減少迭代次數(shù)，成為限制高保真中子輸運方法發(fā)展的新挑戰(zhàn)。為研究有限差分法漸近收斂速度，傅里葉分析是一種廣泛應用的標準方法[11]，研究結果表明對于傳統(tǒng)CMFD方法，當適當減小光學厚度時可獲得更好的收斂性。

本文針對基于廣義平衡理論的粗網(wǎng)有限差分(gCMFD)方法[3]進行了傅里葉分析，依據(jù)傅里葉分析的結果基于全三維一步法精細化輸運求解程序HNET[12-13]重新構建以1/4柵元為單位的幾何建模模塊，并以此為基礎實現(xiàn)1/4柵元尺度gCMFD兩級加速，從而能夠從降低迭代次數(shù)和減少計算規(guī)模2個方面有效降低高保真輸運計算的計算負擔。采用C5G7 RA和VERA Problem #4三維基準題驗證方法的正確性和加速效果，并基于秦山一期壓水堆不同穩(wěn)態(tài)運行工況開展1/4模型穩(wěn)態(tài)臨界計算，綜合驗證程序的計算能力和加速效果。

1 傅里葉分析與兩級gCMFD加速

在HNET中，目前采用柵元尺度的兩級gCMFD方法加速2D/1D精細化中子輸運計算：二維采用MOC計算為一維NEM和三維多群gCMFD提供粗網(wǎng)均勻化宏觀截面、局部徑向節(jié)塊不連續(xù)因子(nodal discontinuity factor，NDF)和擴散系數(shù)修正因子(modified diffusion coefficient factor，MDF)，一維采用NEM計算為三維gCMFD提供局部軸向NDF和MDF，通過這些參數(shù)構建基于柵元均勻化網(wǎng)格的三維多群CMFD線性系統(tǒng)。再次求解全局三維多群gCMFD，為徑向二維MOC計算提供柵元的平均中子通量密度、軸向泄漏項、邊界偏中子流和特征值，同時，為軸向一維NEM計算提供節(jié)塊的平均通量和徑向泄漏項。然而，多群gCMFD線性系統(tǒng)自身條件數(shù)大，同樣面臨求解耗時嚴重的問題。基于同樣的思路，應用gCMFD方法建立與多群gCMFD完全等價的單群gCMFD線性系統(tǒng)，以加速多群gCMFD的迭代求解過程，能夠?qū)崿F(xiàn)在能量方面對多群gCMFD求解的有效加速，進而減少其計算負擔。相比于多群gCMFD，單群gCMFD線性系統(tǒng)自身條件數(shù)很低，計算代價也很低。因此，不會增加過多的額外計算負擔。

1.1 gCMFD傅里葉分析

傅里葉分析是研究有限差分法漸近收斂速度的一種廣泛應用的標準方法，已被廣泛應用于CMFD加速方法的分析。圖1為按照文獻[11]中的方法繪制的gCMFD譜半徑隨著光學厚度的變化曲線，并以C5G7基準題提供的宏觀截面和柵元尺寸數(shù)據(jù)，針對較難收斂的慢化劑材料標注平源區(qū)、1/4柵元以及柵元尺度對應的位置。其中，光學厚度ΣtΔ為分群總截面與柵元尺寸大小的乘積，θ為計算修正擴散系數(shù)時人為定義的系數(shù)，而使用修正擴散系數(shù)可以使 CMFD 的收更加平滑，譜半徑ρ的意義為衰減系數(shù)矩陣的特征值的絕對值的最大值：

(1)

(2)

若譜半徑小于1，則說明收斂；若大于1則說明會發(fā)散；且譜半徑越小說明加速效果越好。通過圖1可以看到，隨著粗網(wǎng)光學厚度減小，gCMFD譜半徑也在變小，意味著粗網(wǎng)格尺寸越小，收斂性越好，但當粗網(wǎng)格尺寸繼續(xù)減小時，收斂性變化趨于平緩，而實現(xiàn)難度和計算負擔會成倍增加。綜合考慮加速收益與實現(xiàn)復雜度，1/4柵元尺度的gCMFD為最值得嘗試的選擇。傅里葉分析的結果為通過建立基于1/4柵元的gCMFD加速系統(tǒng)來進一步改善收斂性能提供了重要的理論依據(jù)。

圖1 gCMFD在不同光學厚度下的譜半徑Fig.1 Spectral radius of gCMFD at different optical thicknesses

1.2 1/4柵元尺度的兩級gCMFD加速方法

1/4柵元尺度的多群gCMFD方法的推導從多群穩(wěn)態(tài)中子輸運方程的形式開始：

∑t,g(r)ψg(r,Ω)=Qg(r,Ω)

(3)

式中：ψ、∑t、Q分別是中子角通量密度、總截面和源項；下標g為能群編號。

將式(3)在空間和1/4 柵元內(nèi)積分即可得到1/4 柵元尺度的3D多群離散擴散方程：

(4)

交界面處從中心節(jié)塊c流向鄰近節(jié)塊n的中子流為：

(5)

(6)

(7)

MDF和NDF的求解可以通過使用來自徑向MOC和軸向NEM計算得到的1/4 柵元的表面中子流、通量和網(wǎng)格平均通量信息來表征，具體計算過程可以參考文獻[3]，在此不再贅述。

對于1/4柵元尺度的多群gCMFD方法中使用的均勻化參數(shù)，可以根據(jù)高階輸運計算的與區(qū)域相關的宏觀截面和通量信息來量化：

(8)

(9)

(10)

(11)

(12)

(13)

式中下標k是平源區(qū)的編號。如此，便可以建立1/4柵元尺度的多群gCMFD的線性系統(tǒng)：

(14)

由于多群gCMFD本身的條件數(shù)較大，如果利用Wielandt Shift加速逆冪迭代，則條件數(shù)會變得更大，導致一級多群gCMFD線性系統(tǒng)的收斂速度仍然較慢，而單群gCMFD線性系統(tǒng)的求解成本要低得多。因此基于同樣的思路，在1/4柵元尺度的多群gCMFD線性系統(tǒng)下通過在能群上歸并來自多群gCMFD的截面、通量和中子流信息來建立1/4柵元尺度的單群gCMFD線性系統(tǒng)，實現(xiàn)加速收斂：

(15)

1/4柵元尺度的多群gCMFD方法的建立流程與柵元尺度的多群gCMFD基本相同，只不過前者是基于1/4柵元為粗網(wǎng)大小進行求解，而不是基于完整的柵元，如圖2所示。

圖2 1/4柵元尺度的兩級CMFD加速方法Fig.2 1/4 cell based multi-level CMFD acceleration method

即便是如此微小的處理，獲得的收益也是顯著的。在粗網(wǎng)尺寸變?yōu)?/4柵元后，MOC與gCMFD中傳遞的標通量信息增多，也更加精細，減小了其間的映射誤差，從而實現(xiàn)良好的加速性能，但gCMFD本身的計算負擔也同時增加。

2 基于1/4柵元的幾何建模及特征線布置

構造實體幾何(constructive solid geometry，CSG)方法是許多高級建模軟件(如CAD)的首選方法，利用不同結構化幾何間簡單的布爾運算就能很方便地得到不同的幾何體。因具有建模靈活，能夠降低幾何信息的存儲的特性，非常適合高度結構化且重復度很高的反應堆模型的建模。

2.1 1/4柵元為單位的幾何建模

HNET原有的幾何模塊是以柵元為單位的，只能構建完整的組件，無法適應當沿x或y方向的組件數(shù)為奇數(shù)，如17×17的布置時1/4堆芯的建模和計算。在這種情況下，在堆芯切割面處中存在1/2組件和1/4組件，而當組件內(nèi)柵元的布置也是奇數(shù)，如15×15時，切割面處便存在1/2柵元和1/4柵元。為實現(xiàn)以1/4柵元為單位來逐層組裝呈堆芯，利用CSG的思想，可通過定義柵元中慢化劑的邊界面、燃料組件的圓柱面即可構建一個柵元的模型，再將柵元的模型的中心創(chuàng)建2個互相垂直的面，將其劃分成為1/4柵元模型，并按照象限號進行區(qū)分和存儲，1/4柵元按照象限編號索引排列即可還原為完整的柵元，進而構成二維的組件，二維組件再通過軸向上的堆疊即可組成三維組件的模型，進而完成全堆芯或1/4堆芯的建模，結構層次如圖3所示。值得注意的是，當以1/4柵元為單位建模后，軸向的NEM計算也基于1/4柵元進行。

2.2 1/4柵元模塊化特征線布置

模塊化特征線可以大幅減少MOC的計算量內(nèi)存消耗，同時還可以精確地處理反射邊界條件。在柵元模塊化特征線模式中，只有外形尺寸和網(wǎng)格劃分完全相同的柵元，共用同一套特征線的幾何信息，如何保證特征線在相鄰的模塊間精準連接是模塊化特征線布置的重點和難點。針對對稱切割模型的特征線布置，合理地布置射線仍然可以使特征線精確連接，具體可參考文獻[14]。假設離散的方位角數(shù)量為N、初始射線間距為t、模塊寬度為w、高度為h。將整個方位角區(qū)域[0,2π]進行等分，若取n∈{1,2,…,N}，則方位角αn為：

(16)

圖3 1/4堆芯結構層次圖Fig.3 Hierarchical diagram of the 1/4 reactor core structure

得到方位角αn后，為保證整柵元和1/4柵元特征線布置方式相同，應保證在整柵元中沿x方向和y方向的特征線數(shù)目nx和ny為偶數(shù)，整柵元的nx和ny為：

(17)

(18)

對應的1/4柵元的nx和ny分別為：

(19)

(20)

得知nx和ny后，為了保證射線的均勻和對稱，需要對方位角αn進行修正，得到修正后的有效方位角αn,eff為：

(21)

同樣，修正后的射線間距teff為：

(22)

3 三維基準題數(shù)值驗證

本文首先基于C5G7-RA、VERA Problem #4三維基準題驗證建模和加速方法的正確性和效果。

3.1 C5G7-RA三維基準題驗證

C5G7三維基準題是OECD/NEA發(fā)布的用于驗證輸運求解器非均勻堆芯計算能力的的基準題。該基準題根據(jù)控制棒插入位置和深度，定義了3種堆芯布置，分別為提棒(UR)、插棒A(RA)和插棒B(RB)，3種堆芯布置的非均勻性和計算難度依次增大。本文選取RA來進行加速效果的驗證和分析。其1/4堆芯布置可參考文獻[15]。

圖4 模塊間特征線連接示意Fig.4 Connection of characteristic rays among different modules

計算過程中，軸向燃料區(qū)域劃分為18層，反射層區(qū)域劃分為18層；燃料、控制棒及裂變室3種柵元，徑向共劃分為5環(huán)(其中燃料區(qū)域3環(huán))、輻向劃分為8個扇形分區(qū)；反射層區(qū)域內(nèi)采用矩形網(wǎng)格劃分；截面來源于基準題報告提供的7群宏觀截面；射線間距選擇0.01 cm，一個卦限內(nèi)選取18方位角和3極角，采用Tabuchi-Yamamoto極角求積組，收斂準則為keff和裂變源均小于1×10-5，核數(shù)布置為徑向1個核，軸向18個核，每層由1個核進行計算，總共18核。

為驗證1/4柵元尺度的多群gCMFD與柵元尺度的多群gCMFD2種方法的加速效果，同時探究1/4柵元尺度的單群gCMFD對1/4柵元尺度的多群gCMFD的加速效果，分別構建了柵元尺度的多群gCMFD(Cell-MG-gCMFD)、1/4柵元尺度的多群gCMFD(1/4Cell-MG-gCMFD)以及1/4柵元尺度的兩級gCMFD(1/4Cell-TL-gCMFD)3個算例進行計算，在每種計算條件下按恰好到達最小外迭代次數(shù)設定多群和單群的最大迭代次數(shù)，統(tǒng)計MOC和gCMFD的迭代次數(shù)與計算時間以及有效增殖因子，并以Cell-MG-gCMFD方案的總計算時間為基準，比較和統(tǒng)計了其他方案的總計算時間相較基準改進的倍數(shù)，定義為加速比，如表1中所示。

表1 C5G7 RA三維基準題計算結果Table 1 The computed result for the iterations and computational time for 3D C5G7 RA case

由表1可知，使用1/4Cell-MG-gCMFD方法，MOC的最少迭代次數(shù)從22次減為12次，成功打破了Cell-MG-gCMFD方法的加速瓶頸，加速比為1.434。然而，由于MG gCMFD的計算規(guī)模和迭代次數(shù)增加，其計算時間也變?yōu)樵瓉淼?2倍，但是這種耗時的增加又可以通過代價很小的1/4Cell-CG-gCMFD方法進一步限制下來，這樣兩級gCMFD加速方法構成的1/4Cell-TL-gCMFD加速系統(tǒng)的加速比進一步提升至1.518，且計算精度與單級相同。圖5展示了在柵元尺度和1/4柵元尺度下當令MG gCMFD完全收斂時MOC迭代次數(shù)變化，可以清晰觀察到1/4Cell-MG-gCMFD對MOC收斂顯著的加速效果。圖6展示了每一個外迭代內(nèi)最后一步MG gCMFD對應迭代次數(shù)與其收斂殘差的變化趨勢，進一步說明了1/4Cell-MG-gCMFD雖然可以加速外迭代的收斂，但自身收斂難度增加，而加入1/4柵元尺度的單群CMFD構成兩級CMFD加速方法后，其迭代次數(shù)和收斂效果都得到了明顯的改善。

圖5 不同加速方案MOC迭代收斂曲線Fig.5 The iteration history of MOC accelerated by different method

3.2 VERA Problem #4三維基準題驗證

為進一步探究1/4柵元尺度的單群gCMFD對1/4柵元尺度的多群gCMFD計算代價的縮減效果，同時驗證本文實現(xiàn)的1/4建模方法的正確性，本文采用基于HELIOS-47能群的VERA Problem #4基準題進行測試。此問題組件的軸向結構和徑向布置可參考文獻[16]。

圖6 MG-gCMFD迭代收斂歷史Fig.6 The convergence history of MG-CMFD

計算時，特征線間距為0.03 cm，每卦限內(nèi)采用16個方位角和3個極角，特征值的收斂標準為1.0×10-5，裂變源采用無窮范數(shù)誤差，其收斂標準為1.0×10-5。MG gCMFD和CG gCMFD的特征值收斂標準為1.0×10-8。徑向為反射邊界，軸向為真空邊界。每層分配一個核，共使用58核進行計算。

針對此基準題，本文基于HNET程序構建全模型下基于柵元的兩級gCMFD(Cell-TL-gCMFD-全模型)、全模型下基于1/4柵元的兩級gCMFD加速(1/4Cell-TL-gCMFD-全模型)以及1/4模型下基于1/4柵元的兩級gCMFD加速(1/4Cell-TL-gCMFD-1/4模型)3個算例進行計算，其中全模型指計算完整的3×3組件模型，而1/4模型為3×3模組件型的1/4模型。在每種算例下按單群CMFD充分收斂(最大迭代次數(shù)設為1 000)的情況下恰好到達最小外迭代次數(shù)來設定多群CMFD最大迭代次數(shù)，統(tǒng)計MOC和CMFD的迭代次數(shù)與計算時間以及有效增殖因子，并計算了加速比，如表2中所示。Cell-TL-gCMFD-全模型與1/4Cell-TL-gCMFD-全模型2個算例的有效增殖因子與參考值的偏差均為15×10-5，表明本文構建的1/4模型與原有幾何構建方法計算精度相同。同時，由于應用基于1/4柵元的兩級CMFD加速，MOC的迭代次數(shù)從11次減少到8次，但由于CMFD的計算量大幅增加，兩級CMFD的計算時間變?yōu)樵瓉淼?倍，然而在增加了特征線交點信息與CMFD的計算量的情況下加速比仍達到了1.651，較7群有了明顯的提高。而在1/4模型下，由于MOC的計算量變?yōu)槿Ｐ偷?/4，相較基于柵元的兩級CMFD加速計算，平均每個迭代步下MOC的計算時間從81.83 s減少為19.875 s，減少了3/4，同時總迭代次數(shù)減少了3次，CMFD的計算量相差無幾，但在1/4柵元尺度上收斂性能提升，整體計算加速比達到4.294。

表2 VERA Problem #4三維基準題計算結果Table 2 The computed result for the iterations and computational time for VERA Problem #4 case

在C5G7基準題中，Cell-MG-gCMFD方案平均每個迭代步下MOC計算時間為19.93 s，1/4Cell-MG-gCMFD方案為23.51 s，時間稍有所提升，而在VERA基準題中，時間分別為81.83 s和68.16 s，時間明顯下降。實際上，在以1/4柵元為單位構建幾何模型后，模塊化特征線的單位也變?yōu)?/4柵元，MOC的總計算量和信息存儲量沒有發(fā)生改變，但是索引的序列變?yōu)樵瓉淼?倍，導致在特征線連接時時間稍有所增加，但不會增加很多。而在VERA基準題中，單次外迭代的計算更為困難，但基于1/4柵元的多級CMFD加速可以為其提供更好的迭代初值，能夠顯著的加速其收斂計算，總體的收益遠大于由于索引序列增加帶來的計算時間增加。

4 實堆驗證

為進一步表明基于1/4柵元的多級CMFD加速方法計算1/4模型的優(yōu)越性，本文針對秦山一期核電廠首爐堆芯在壽期初、熱態(tài)零功率工況下，展開了1/4堆芯的穩(wěn)態(tài)計算。

4.1 秦山1/4堆芯模型

秦山一期核電廠由313個組件組成，其中每個組件均由15×15的柵元組成，柵元間距為1.33 cm，活性區(qū)軸向高度為290 cm。本文構建的秦山1/4堆芯模型，完整包含了燃料棒、毒物棒、控制棒、導向管、定位格架、圍板和反射層等結構。同時考慮到上下反射層區(qū)域等，在軸向上劃分為62層。堆芯內(nèi)燃料組件按秦山一期核電廠首爐堆芯燃料裝載方案布置，模型的徑向幾何示意圖如圖7所示。

圖7 秦山1/4堆芯徑向幾何示意Fig.7 Radial geometry configurtion for Qinshan 1/4 core

計算時，特征線間距為0.05 cm，每卦限內(nèi)采用10個方位角和3個極角的組合，特征值的收斂標準為1.0×10-6，裂變源采用無窮范數(shù)誤差，其收斂標準為1.0×10-4。MG gCMFD的特征值收斂標準為1.0×10-10，最大迭代次數(shù)為4；CG gCMFD的特征值收斂標準為1.0×10-10，最大迭代次數(shù)為50。徑向?qū)ΨQ軸處為反射邊界，其余均為真空邊界。在山河超算平臺使用62個節(jié)點進行計算，每個節(jié)點使用16核計算一層堆芯模型，共使用992核進行計算。

4.2 數(shù)值結果

本文分別針對ARO工況和4種控制棒束插入工況進行了穩(wěn)態(tài)臨界計算，其中4根控制棒分別命名為T1、T2、T3、T4，初始狀態(tài)為未插入，插入狀態(tài)為in。在令有效增殖因子為1(誤差為10-7×10-5)時，得到了各工況下的末端硼濃度。對比《數(shù)字化反應堆30萬千瓦核電廠基準問題模擬計算報告》中的參考值，各工況下的末端硼濃度如表3所示，最大偏差在-21×10-6，低于通用的設計準則50×10-6及設計院針對程序驗證準則30×10-6，滿足精度要求。

表3 秦山1/4三維堆芯模型末端硼濃度計算結果

同時針對ARO工況采用報告中給出的硼濃度數(shù)據(jù)進行高保真穩(wěn)態(tài)臨界計算，得到特征值和歸一化柵元功率分布。為了表明本文提出并實施的方法在實堆應用中的效果，分別統(tǒng)計了原有的柵元尺度下全模型以及1/4柵元尺度下1/4模型的總計算時間。臨界偏差為-191×10-5，特征值為0.998 09。如表4所示，2種計算方案得到的特征值均為0.998 09，與臨界值1的偏差為-191×10-5,低于通用的設計準則500×10-5，滿足精度要求。同時，基于1/4柵元多級CMFD加速方法的1/4模型高保真全堆芯輸運計算外迭代次數(shù)從12次減少為9次，總時間為88.3核時，加速比為2.637。

表4 秦山核電站三維堆芯模型ARO工況臨界計算結果

計算得到的ARO工況下的歸一化柵元功率分布如圖8所示，功率在軸向呈現(xiàn)出明顯的中心高，兩側低的趨勢；而在徑向上，由于燃料分區(qū)布置和毒物棒對功率分布的影響，呈現(xiàn)出功率峰的交替分布。

圖8 不同工況下歸一化柵元功率的三維分布Fig.8 3D distribution of normalized pin power under different conditions

5 結論

1)針對基于1/4柵元的多群CMFD方法計算成本升高，進一步提出并實施基于1/4柵元的單群CMFD方法加速多群CMFD的收斂構成兩級CMFD加速，以小代價獲取大的時間收益。

2)采用三維基準題驗證方法的正確性和加速效果，計算結果表明，本文提出的1/4柵元尺度多級CMFD加速方法成功突破原有加速瓶頸，加速效果顯著。同時基于秦山一期壓水堆模型不同穩(wěn)態(tài)運行工況開展1/4模型高保真穩(wěn)態(tài)臨界計算，驗證了HNET的實堆計算能力，具有較好的計算精度和計算效率。然而1/4柵元尺度的單群CMFD計算代價仍然巨大，未來可進一步基于1/4柵元幾何建模功能開發(fā)基于整柵元的單群CMFD加速1/4柵元多群CMFD實現(xiàn)空間-能量混合加速，進一步提升加速效果。