精細(xì)化中子輸運(yùn)計(jì)算程序HNET動(dòng)力學(xué)計(jì)算初步驗(yàn)證

康樂(lè)， 郝琛

(哈爾濱工程大學(xué) 核安全與仿真技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，黑龍江 哈爾濱 150001)

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

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

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

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

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

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

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

(1)

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

(2)

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

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

(3)

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

(4)

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

(5)

其中：

(6)

(7)

(8)

(9)

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

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

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

(10)

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

(11)

(12)

(13)

(14)

(15)

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

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

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

(16)

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

2 C5G7-TD基準(zhǔn)題模型

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

2.1 二維算例模型

1)TD0系列。

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

2)TD1與TD2系列。

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

3)TD3系列。

TD3系列描述了二維情況下反應(yīng)堆內(nèi)慢化劑密度變化過(guò)程。在起始時(shí)刻，燃料組件內(nèi)的慢化劑密度(不包括軸向反射層)為正常密度，在0～1 s過(guò)程中，慢化劑密度以一定速度降低。在1 s結(jié)束時(shí)，密度達(dá)到最小值，即初始密度的w值。然后慢化劑密度以同樣的速度增加，并在2 s結(jié)束時(shí)回到初始值。對(duì)于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系列描述了三維情況下控制棒線性插入與提出反應(yīng)堆的過(guò)程。在初始時(shí)刻，所有控制棒均處于頂部反射層內(nèi)。隨后在0～8 s過(guò)程中，不同組件內(nèi)的控制棒以不同的速度勻速地插入與提出反應(yīng)堆。TD4系列共包含5個(gè)算例，分別模擬不同組件內(nèi)控制棒不同的移動(dòng)情況。

2)TD5系列。

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

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

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

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

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

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

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

圖3 TD1系列反應(yīng)性歷史Fig.3 TD1 reactivity history

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

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

圖6 TD2系列反應(yīng)性歷史Fig.6 TD2 reactivity history

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

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

圖9 TD3系列反應(yīng)性歷史Fig.9 TD3 reactivity history

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

表1 TD1～TD3系列算例運(yùn)行時(shí)間對(duì)比Table 1 Comparison of total run times for TD1～TD3 h

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

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

對(duì)于三維問(wèn)題TD4和TD5系列，堆芯相對(duì)功率變化和反應(yīng)性變化見(jiàn)圖10～15。其中，“TD4-4”和“TD4-5”描述了不同控制棒同時(shí)插入和提出的過(guò)程，這導(dǎo)致了相對(duì)復(fù)雜的裂變率變化和反應(yīng)性變化。表2給出了所有三維算例的運(yùn)行時(shí)間并與MPACT參考解對(duì)比。

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

表2 TD4、TD5算例的運(yùn)行時(shí)間Table 2 Run time summaries for TD4， TD5 h

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

圖12 TD4系列反應(yīng)性歷史Fig.12 TD4 reactivity history

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

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

圖15 TD5系列反應(yīng)性歷史Fig.15 TD5 reactivity history

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

4 結(jié)論

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

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

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

同時(shí)，用于瞬態(tài)計(jì)算的熱工反饋功能正在研發(fā)中，未來(lái)將進(jìn)一步測(cè)試HNET的時(shí)空動(dòng)力學(xué)計(jì)算能力。