康樂， 郝琛

(哈爾濱工程大學 核安全與仿真技術(shù)國防重點學科實驗室，黑龍江 哈爾濱 150001)

隨著高性能計算能力的快速提升，反應堆數(shù)值模擬技術(shù)在精細化方面的需求逐漸提高。然而，直接采用標準三維中子輸運方法開展全堆芯計算仍不現(xiàn)實，主要是因為典型反應堆堆芯的實際未知數(shù)總數(shù)過大，對于穩(wěn)態(tài)計算接近1015數(shù)量級，而對于瞬態(tài)計算來說則要更大。目前，一種實用且切實可行的方案是采用二維/一維方法[1-4]，即通過定義軸向泄漏項和徑向泄漏項，將全三維輸運計算轉(zhuǎn)為徑向二維輸運計算耦合軸向一維低階輸運計算。同時，精細化全三維瞬態(tài)輸運計算依賴于有效的全局加速方法的應用。其中，粗網(wǎng)有限差分方法(coarse mesh finite difference, CMFD)是一種流行的全局加速方法，特別針對二維/一維聚合算法，CMFD還為徑向二維輸運計算和軸向一維低階輸運計算提供了一種天然聚合計算框架。但是，傳統(tǒng)CMFD有收斂慢的缺陷，同時還存在收斂發(fā)散及耦合系數(shù)為負的潛在缺陷[5]。為了有效降低多群CMFD帶來的計算負擔，在研發(fā)HNET程序(high-fidelity neutron trasport program for 3D nuclear reactor core，HNET)過程中，提出了基于廣義等價理論的兩級瞬態(tài)CMFD加速方案，通過建立等價的單群CMFD來加速多群CMFD 瞬態(tài)固定源問題的收斂，并進一步加速精細化全三維瞬態(tài)輸運計算。基于CMFD并行計算框架，徑向二維輸運計算采用流行的特征線方法(method of characteristics, MOC)，軸向一維低階輸運計算采用節(jié)塊展開法(nodal expansion method, NEM)。

為了驗證HNET的時空中子動力學計算功能，本文采用C5G7-TD非均勻瞬態(tài)基準題對HNET的計算精度進行驗證，計算結(jié)果與高保真中子學程序MPACT進行對比。C5G7-TD系列基準題[6-7]包括一系列非均勻幾何時空動力學基準題，用于校驗動力學程序?qū)Ψ蔷鶆驇缀螁栴}的計算能力和計算精度。目前該基準問題已經(jīng)發(fā)布第一階段的系列問題，用于驗證時空中子學輸運算法與程序。

1 HNET瞬態(tài)計算方法

HNET采用二維MOC/一維NEM的聚合算法，并通過基于廣義等價理論的兩級瞬態(tài)CMFD加速全局三維瞬態(tài)輸運計算，同時還采用基于MPI的區(qū)域分解并行技術(shù)，進一步縮減運行時間，更多細節(jié)見文獻[5,8-11]。

1.1 三維瞬態(tài)輸運固定源方程

為了簡要說明三維瞬態(tài)固定源方程的推導過程，本文從瞬態(tài)多群中子輸運方程出發(fā)，在各向同性散射近似的條件下，輸運方程為：

Σt,g(r,t)φg(r,Ω,t)+Ss,g(r,Ω,t)+

(1)

緩發(fā)中子先驅(qū)核方程為：

(2)

式中：v為中子飛行速度；φ為角通量密度；Ck為第k緩發(fā)群的先驅(qū)核濃度；SF為經(jīng)過穩(wěn)態(tài)修正后的裂變源項；Ss,g為散射源項；χp和χd分別為瞬發(fā)和緩發(fā)中子能譜；βk為緩發(fā)中子份額；λk為緩發(fā)中子先驅(qū)核衰變常數(shù)；r為空間向量；Ω為角度向量；g為能群；Σt,g為總截面。

對瞬態(tài)角通量做指數(shù)變換：

(3)

其中，系數(shù)αn為：

(4)

將式(4)代入式(1)～(3)中，并且對時間導數(shù)項采用全隱式方法，則瞬態(tài)輸運方程可轉(zhuǎn)化為瞬態(tài)固定源方程為：

(5)

其中：

(6)

(7)

(8)

(9)

至此，式(6)可以使用任何標準的中子固定源問題求解器來解決。在HNET程序中，徑向二維輸運計算采用流行的MOC方法，軸向一維低階輸運計算采用節(jié)塊展開NEM方法，并通過基于廣義等價理論的兩級瞬態(tài)CMFD加速全局三維瞬態(tài)輸運計算。

1.2 基于廣義等價理論的兩級瞬態(tài)粗網(wǎng)有限差分加速方法

CMFD不僅為二維徑向MOC輸運計算和一維軸向NEM計算提供了天然的計算框架，同時還對全堆芯輸運計算起到加速作用，三維多群瞬態(tài)擴散方程為：

(10)

根據(jù)基于廣義等價理論的粗網(wǎng)有限差分方法，中心節(jié)塊c到相鄰節(jié)塊q的界面流為：

(11)

(12)

(13)

(14)

(15)

式中：W、E、N、S、T、B代表節(jié)塊在西、東、北、南、頂、底方向上的6個表面。

式(15)可以由任意固定源問題求解器來求解，且所用的特征值來自于穩(wěn)態(tài)輸運計算，其在整個瞬態(tài)過程中應保持不變。

由于多群CMFD線性系統(tǒng)的條件數(shù)很大，如若采用Wielandt shift方法加速逆冪迭代，其條件數(shù)會變得更大，因此，一級多群CMFD線性系統(tǒng)存在收斂慢的缺陷。一群CMFD線性系統(tǒng)條件數(shù)較低使得其容易求解，可利用一群CMFD獲得中子通量及裂變源加速多群CMFD?；诙嗳篊MFD獲得的中子通量、中子流及多群截面信息，可獲得一群CMFD線性系統(tǒng)為：

(16)

一群CMFD用以加速多群CMFD的特征值和裂變源項，求解多群CMFD線性系統(tǒng)，以計算得到新的多群中子通量分布和特征值，進而用于更新一群宏觀截面、中子通量和中子流。同時，計算得到的多群中子通量將進一步更新每一個節(jié)塊的界面流，而界面流則會被用于更新一群CMFD的節(jié)塊不連續(xù)因子和擴散系數(shù)修正因子，直到整個線性系統(tǒng)收斂。其中，當殘差范數(shù)下降到指定標準或最大迭代次數(shù)達到時，一群CMFD迭代或多群CMFD迭代相應退出，整體收斂通過檢查不連續(xù)因子的相對誤差是否達到收斂標準，具體細節(jié)可參考文獻[7]。

2 C5G7-TD基準題模型

C5G7-TD系列基準題二維算例包含TD0、TD1、TD2和TD3系列；三維算例包含TD4和TD5系列。每個系列都包含數(shù)個算例，這些算例包含控制棒移動和慢化劑密度變化等瞬態(tài)事件。控制棒移動是以組件為單位同時移動，即在同一組件內(nèi)所有控制棒均保持相同的高度。二維與三維的瞬態(tài)問題，更多細節(jié)見文獻[5]。

2.1 二維算例模型

1)TD0系列。

TD0系列算例描述了二維情況下控制棒瞬間插入與提出反應堆的過程。在初始時刻，控制棒突然插入到堆芯10%的位置，并停留1 s。隨后，控制棒提出到堆芯5%的位置，并再次停留1 s。最后，在2 s結(jié)束時，控制棒從堆芯中完全撤出。在二維模型中，這一過程被模擬為控制棒區(qū)1區(qū)截面的變化。

2)TD1與TD2系列。

TD1與TD2系列描述了二維情況下控制棒線性插入與提出反應堆的過程，通過線性地將導向管截面變化成控制棒截面來模擬控制棒移動過程。在初始時刻，所有的控制棒均在堆芯外部，在0～1 s過程中，控制棒從堆芯頂部勻速插入到堆芯的特定位置。其中，TD1為堆芯5%的位置，TD2為堆芯10%的位置。然后以相同的速度提出控制棒，并在2 s結(jié)束時完全提出堆芯區(qū)域。TD1和TD2基準題各有5道算例。

3)TD3系列。

TD3系列描述了二維情況下反應堆內(nèi)慢化劑密度變化過程。在起始時刻，燃料組件內(nèi)的慢化劑密度(不包括軸向反射層)為正常密度，在0～1 s過程中，慢化劑密度以一定速度降低。在1 s結(jié)束時，密度達到最小值，即初始密度的w值。然后慢化劑密度以同樣的速度增加，并在2 s結(jié)束時回到初始值。對于TD3中的不同算例，w的值有所區(qū)別，其中“TD3-1”的w為0.95；“TD3-2”的w為0.90；“TD3-3”的w為0.85；“TD3-4”的w為0.80。

2.2 三維算例模型

1)TD4系列。

TD4系列描述了三維情況下控制棒線性插入與提出反應堆的過程。在初始時刻，所有控制棒均處于頂部反射層內(nèi)。隨后在0～8 s過程中，不同組件內(nèi)的控制棒以不同的速度勻速地插入與提出反應堆。TD4系列共包含5個算例，分別模擬不同組件內(nèi)控制棒不同的移動情況。

2)TD5系列。

TD5系列描述了三維情況下反應堆內(nèi)慢化劑密度變化過程。在整個瞬態(tài)過程中，所有控制棒均保持在燃料組件頂部的反射層中。隨后在0～8 s過程中，不同組件內(nèi)的慢化劑密度以組件為單位同步地進行變化。TD5系列共包含4個算例，分別模擬不同組件內(nèi)慢化劑變化情況。

3 HNET數(shù)值結(jié)果與分析

對于二維算例，本文展示TD1、TD2與TD3的計算結(jié)果；對于三維算例，本文展示TD4與TD5的計算結(jié)果。對于二維算例，瞬態(tài)輸運計算的步長為10 ms，對于三維算例，時間步長為25 ms。對于所有算例，徑向網(wǎng)格的劃分方式是一致的。對于燃料組件內(nèi)的非均勻柵元，在燃料區(qū)域徑向劃分3圈且周向劃分8區(qū)，在慢化劑區(qū)域徑向劃分1圈且周向劃分8區(qū)；對于反射層內(nèi)的慢化劑柵元，靠近燃料組件的部分柵元劃分為6×6的平源區(qū)網(wǎng)格，遠離燃料組件的部分柵元劃分為1×1的平源區(qū)網(wǎng)格。對于所有三維算例，軸向網(wǎng)格尺寸為 5.355 cm，特征線間距為0.03 cm，軸向求積組采用Yamamoto求積組[12]，半空間內(nèi)采用64個方位角與3個極角的組合。特征值的收斂準則為1.0×10-6，中子標通量采用無窮范數(shù)誤差，其收斂準則為1.0×10-4。算例所使用的CPU型號均為2.60 GHz Intel Xeon E5-2690 v4，二維算例均采用徑向9核的區(qū)域分解并行計算方案，三維維算例均采用軸向32核的區(qū)域分解并行計算方案。為了保證數(shù)值結(jié)果對比的有效性，上述參數(shù)除求積組和CPU型號外均與MPACT的參考解保持一致。

3.1 二維算例計算結(jié)果

對于C5G7-TD基準題的穩(wěn)態(tài)特征，MPACT程序的參考解為1.186 673，HNET程序的計算結(jié)果為1.186 867，HNET結(jié)果與參考解的相對誤差約為0.016%，表明HNET計算得到的初始時刻的堆芯狀態(tài)是正確的。二維算例的堆芯相對功率歷史、與MPACT參考解的堆芯功率誤差和反應性歷史如圖1～9所示。

圖1 TD1系列堆芯相對功率歷史Fig.1 TD1 fractional core power history

圖2 對比MPACT程序的TD1系列堆芯功率相對誤差Fig.2 Relative errors of TD1 fractional core power compared to MPACT results

圖3 TD1系列反應性歷史Fig.3 TD1 reactivity history

圖4 TD2系列堆芯相對功率歷史Fig.4 TD2 fractional core power history

圖5 對比MPACT程序的TD2系列堆芯功率相對誤差Fig.5 Relative error of TD2 fractional core power compared to MPACT results

圖6 TD2系列反應性歷史Fig.6 TD2 reactivity history

圖7 TD3系列堆芯相對功率歷史Fig.7 TD3 fractional core power history

圖8 對比MPACT程序的TD3系列堆芯功率相對誤差Fig.8 Relative error of TD3 fractional core power compared to MPACT results

圖9 TD3系列反應性歷史Fig.9 TD3 reactivity history

對于TD1、TD2和TD3系列，在0～1 s和1～2 s的過程中，宏觀截面以相同的速度增減變化，使得裂變率和反應性的增減也出現(xiàn)了對稱的趨勢。同時，TD2系列的控制棒插入深度大于TD1系列，導致TD2系列引入了數(shù)值更大的負反應性。圖1～9表明相對功率和反應性的變化規(guī)律與MPACT的結(jié)果相符。同時，表1列出了每個二維算例的10 s瞬態(tài)過程的運行時間并與MPACT參考解對比。

表1 TD1～TD3系列算例運行時間對比Table 1 Comparison of total run times for TD1～TD3 h

對于TD1系列，功率峰值的最大偏差為0.48%；功率的最大偏差為0.48%，發(fā)生在“TD1-5”的1 s時刻；反應性峰值最大偏差為0.26%，發(fā)生在“TD1-1”。對于TD2系列，功率峰值的最大偏差為0.21%；功率的最大偏差為0.14%，發(fā)生在“TD2-5”的10 s時刻；反應性峰值的最大偏差為1.26%，發(fā)生在“TD2-3”。對于TD3系列，功率峰值的最大偏差為-0.72%；功率的最大偏差為-0.74%，發(fā)生在“TD3-1”的1.05 s時刻；反應性峰值的最大偏差為1.14%，發(fā)生在“TD3-1”。

3.2 三維算例計算結(jié)果

對于三維問題TD4和TD5系列，堆芯相對功率變化和反應性變化見圖10～15。其中，“TD4-4”和“TD4-5”描述了不同控制棒同時插入和提出的過程，這導致了相對復雜的裂變率變化和反應性變化。表2給出了所有三維算例的運行時間并與MPACT參考解對比。

圖10 TD4系列堆芯相對功率歷史Fig.10 TD4 fractional core power history

表2 TD4、TD5算例的運行時間Table 2 Run time summaries for TD4， TD5 h

圖11 對比MPACT程序的TD4系列堆芯功率相對誤差Fig.11 Relative error of TD4 fractional core power compared to MPACT results

圖12 TD4系列反應性歷史Fig.12 TD4 reactivity history

圖13 TD5系列堆芯相對功率歷史Fig.13 TD5 fractional core power history

圖14 對比MPACT程序的TD5系列堆芯功率相對誤差Fig.14 Relative error of TD5 fractional core power compared to MPACT results

圖15 TD5系列反應性歷史Fig.15 TD5 reactivity history

對于TD4系列，功率峰值的最大偏差為-0.21%，出現(xiàn)在“TD4-1”；功率的最大偏差為-0.55%，發(fā)生在“TD4-5”的3.375 s時刻；穩(wěn)定之后相對偏差為0.10%。反應性峰值最大相對偏差為2.21%，發(fā)生在“TD4-4”。對于TD5系列，功率峰值的最大偏差為-0.21%，出現(xiàn)在“TD5-1”；功率最大偏差為-0.23%，發(fā)生在“TD5-1”的2.3 s時刻；穩(wěn)定之后相對偏差為-0.05%。反應性峰值最大相對偏差為0.66%，發(fā)生在“TD5-4”。

4 結(jié)論

1) 對于二維問題，在瞬態(tài)事件期間堆芯功率的相對誤差有所波動，瞬態(tài)事件結(jié)束后相對誤差逐漸趨于平穩(wěn)。雖然計算模型的細微差別可能引入了一些偏差，如空間幾何模型、離散策略、動力學參數(shù)的處理和加速方法等，但是總體而言，HNET的數(shù)值結(jié)果與MPACT基準解吻合得很好，堆芯功率的相對誤差最大不超過0.8%，在保持足夠精度的同時HNET程序算例擁有明顯的效率優(yōu)勢。

2) 對于三維問題，在TD4問題上HNET的相對誤差不超過0.55%，在TD5問題上相對誤差小于0.23%，在輸運時間步長相同的情況下，HNET程序的計算總耗時更少，與MPACT程序相比HNET程序只需要50%～85%的計算時間便可實現(xiàn)相同精度的瞬態(tài)輸運計算。

3) 總體而言，HNET瞬態(tài)過程中堆芯總功率變化與反應性變化均與MPACT符合良好，已完全具備三維瞬態(tài)精細化計算的能力。

同時，用于瞬態(tài)計算的熱工反饋功能正在研發(fā)中，未來將進一步測試HNET的時空動力學計算能力。