高溫氣冷堆堆芯流場高效全局求解方法研究

周夏峰，李 富

(清華大學(xué) 核能與新能源技術(shù)研究院 先進(jìn)反應(yīng)堆工程與安全教育部重點(diǎn)實(shí)驗(yàn)室，北京 100084)

球床模塊式高溫氣冷堆核電站示范工程(HTR-PM)已列入國家科技重大專項(xiàng)，其中HTR-PM堆芯流場的計(jì)算模擬是專項(xiàng)的重要部分之一。在高溫氣冷堆(HTGR)中，由于冷卻劑氦氣具有溫度變化劇烈、密度變化大、流量分配上大流量和小流量并存，以及堆芯內(nèi)流動(dòng)區(qū)域相互耦合等特點(diǎn)，給HTGR的堆芯流場計(jì)算提出挑戰(zhàn)[1]。目前用于HTGR的安全分析程序是THERMIX軟件包，THERMIX系列程序中的KONVEK程序用于堆芯內(nèi)流場的計(jì)算。考慮到堆芯內(nèi)各部件的結(jié)構(gòu)形式和流動(dòng)特性，程序中設(shè)置了5種不同結(jié)構(gòu)的流動(dòng)區(qū)供用戶選擇，分別為不流動(dòng)區(qū)、球床區(qū)、豎管一維流動(dòng)區(qū)、換熱器區(qū)及空腔區(qū)，不同區(qū)之間的流場相互耦合[2]。該軟件包得到了德國GRS評(píng)審和我國核安全當(dāng)局的認(rèn)可，目前用于HTGR的安全分析[3]。

但由于程序編制年代較早，KONVEK程序采用了分塊迭代方法計(jì)算相互耦合的全堆芯流場，即上述5種不同結(jié)構(gòu)的流動(dòng)區(qū)分別獨(dú)立計(jì)算各自的流場，之后根據(jù)不同結(jié)構(gòu)流動(dòng)區(qū)之間的界面相互傳遞信息，不斷迭代更新，直到收斂。分塊迭代方法計(jì)算相互耦合的全堆芯流場時(shí)，部分工況收斂速度慢，耗時(shí)長，甚至?xí)霈F(xiàn)不收斂的問題。因此，本工作研究高效的全局求解方法計(jì)算HTGR堆芯流場分布。

1 數(shù)學(xué)物理模型

1.1 堆芯流場模型

HTGR一回路的冷卻劑是氦氣，由于氦氣密度小，流動(dòng)慣性力也非常小，因此它對(duì)流動(dòng)條件的響應(yīng)非常快，這意味著在每個(gè)短時(shí)間段都能調(diào)整到新平衡態(tài)，因此，在時(shí)間間距不太大時(shí)，可將堆芯流場按照準(zhǔn)穩(wěn)態(tài)處理，忽略時(shí)間相關(guān)項(xiàng)[4]。最終得到HTGR堆芯流場計(jì)算模型：

(1)

其中：Gr和Gz分別為質(zhì)量流量G在r和z方向上的分量，G=ρu，ρ為氣體密度，u為速度矢量；Φ為質(zhì)量源項(xiàng)；W為流動(dòng)阻力系數(shù)，也依賴于氣體流速，W=W(G)，由經(jīng)驗(yàn)關(guān)系式給出，堆芯內(nèi)不同流動(dòng)區(qū)的W不同；p為相對(duì)壓力；g為重力加速度。

由式(1)可得壓力Possion方程的計(jì)算關(guān)系式：

(2)

1.2 模型離散

圖1 壓力控制容積的網(wǎng)格離散示意圖

選取第(I，N)個(gè)網(wǎng)格，即壓力控制容積，對(duì)式(2)中的壓力場方程進(jìn)行積分得：

A1(I,N)(pI,N-pI,N-1)+A2(I,N)(pI,N-

pI-1,N)+A3(I,N)(pI,N-pI,N+1)+

A4(I,N)(pI,N-pI+1,N)=ψI,N+

A2(I,N)ρI-1gΔzI-A4(I,N)ρIgΔzI+1

(3)

其中：

A3(I,N)=A1(I,N+1)

A4(I,N)=A2(I+1,N)

1.3 KONVEK程序中堆芯流場計(jì)算流程

式(3)中A1(I,N)、A2(I，N)、A3(I，N)含有阻力系數(shù)W，W為速度場的函數(shù)，因此堆芯流場的計(jì)算必然帶有迭代性質(zhì)。HTGR堆芯流場分布和耦合情況示于圖2。由圖2可知，HTGR堆芯內(nèi)具有不同結(jié)構(gòu)形式的流動(dòng)區(qū)，各區(qū)的流場之間相互耦合，這使得在計(jì)算時(shí)需考慮耦合問題的求解方法。KONVEK程序采用分塊迭代方法計(jì)算相互耦合的全堆芯流場，即不同結(jié)構(gòu)流動(dòng)區(qū)分別獨(dú)立計(jì)算各自流場，之后根據(jù)不同結(jié)構(gòu)流動(dòng)區(qū)之間界面相互傳遞信息，不斷迭代更新，直至收斂。流場計(jì)算流程如圖3所示。

圖2 HTGR堆芯流場分布和耦合情況

2 目前計(jì)算方法存在的問題

對(duì)于如圖2所示的堆芯流場分布與耦合情況，堆芯中具有不流動(dòng)區(qū)、球床區(qū)、豎直一維流動(dòng)區(qū)和空腔區(qū)，各流動(dòng)區(qū)之間通過邊界相互耦合。依據(jù)堆芯不同結(jié)構(gòu)的流動(dòng)區(qū)，離散后的方程組寫成以下形式：

(4)

其中：A為方程離散后形成的矩陣；X1為所有空腔區(qū)內(nèi)各網(wǎng)格點(diǎn)壓力所形成的向量；X2為所有球床區(qū)內(nèi)各網(wǎng)格點(diǎn)壓力所形成的向量；X3為所有豎直一維流動(dòng)區(qū)內(nèi)各網(wǎng)格點(diǎn)壓力所形成的向量；向量b1、b2、b3分別為X1、X2、X3對(duì)應(yīng)的已知源項(xiàng)。

圖3 KONVEK程序中堆芯流場計(jì)算流程

2.1 分塊迭代方法存在的問題分析

(5)

其中，n為迭代步。

對(duì)式(5)進(jìn)行求解可得：

設(shè)G為：

為保持分塊迭代的收斂性，需要使G的譜半徑ρr(G)<1，且ρr(G)越小，收斂越快[6]，這與傳統(tǒng)意義上的Gauss-Seidel迭代方法的原理相同，存在一定局限性。ρr(G)接近1或ρr(G)≥1時(shí)，會(huì)出現(xiàn)收斂慢，甚至不收斂的現(xiàn)象，計(jì)算時(shí)間會(huì)很長，從而影響計(jì)算精度。

2.2 壓力場和流場計(jì)算方法存在的問題分析

由圖3可知，在特定溫度場下，首先根據(jù)上一步計(jì)算出的速度場得到壓力方程的相關(guān)系數(shù)，之后利用分塊迭代方法分別求解堆芯內(nèi)不同結(jié)構(gòu)流動(dòng)區(qū)的壓力場，之后得到各網(wǎng)格邊界的壓力梯度，再根據(jù)壓力梯度計(jì)算網(wǎng)格邊界的流量，更新速度場，反復(fù)進(jìn)行上述過程，直到收斂為止。也就是說只要計(jì)算出壓力場，速度場即可得到，由此可見，速度場的計(jì)算精度嚴(yán)重依賴于壓力場的求解精度。

此外，壓力場輕微變化將會(huì)引起堆芯內(nèi)速度場的較大變化，所以壓力場的收斂精度不夠或不收斂時(shí)，會(huì)引起速度場始終有相當(dāng)大的波動(dòng)，最終導(dǎo)致速度場收斂較慢，甚至不收斂。

由2.1節(jié)分析可知：采用分塊迭代方法處理相互耦合的堆芯流場，當(dāng)分塊迭代中迭代矩陣G的譜半徑ρr(G)接近1或ρr(G)≥1時(shí)，會(huì)出現(xiàn)收斂慢，甚至不收斂的現(xiàn)象，這樣就很難保證壓力場的求解精度；由于KONVEK程序編制年代較早，在分別求解每一堆芯流動(dòng)區(qū)的壓力場時(shí)，使用了Gauss-Seidel迭代法來進(jìn)行矩陣求解，部分工況計(jì)算時(shí)間較長，同樣也很難保證計(jì)算的收斂性和計(jì)算精度。

綜上所述，計(jì)算HTGR堆芯流場所使用的模型對(duì)壓力場的求解精度需求很高，每次計(jì)算壓力場時(shí)，需保證壓力場計(jì)算方法的收斂性，并盡可能準(zhǔn)確地求解壓力場；分塊迭代方法求解相互耦合的問題可能會(huì)出現(xiàn)收斂速度慢，甚至不收斂的問題，同時(shí)，選用Gauss-Seidel迭代法作為壓力場的矩陣求解器，同樣會(huì)出現(xiàn)上述問題。基于此，選取高效的全局求解方法處理堆芯相互耦合的問題，使用高效的直接求解方法求解壓力場方程。

3 全局求解策略

分塊迭代方法的本質(zhì)是采用Jacobi或Gauss-Seidel方法求解式(4)，這樣可能出現(xiàn)不收斂或收斂速度慢的問題，從而可能導(dǎo)致很難得到其收斂解。全局求解方法的本質(zhì)則是將式(4)看作整體，聯(lián)立求解。

由于直接消去法求解線性方程組時(shí)對(duì)譜半徑和矩陣的條件數(shù)不敏感[7]，可得到高精度的數(shù)值解，即使對(duì)于較為病態(tài)的線性方程組，也能得到相對(duì)合理的解，同時(shí)直接消去法可實(shí)現(xiàn)壓力場的全局求解。因此，為保證壓力場計(jì)算的穩(wěn)定性和準(zhǔn)確性，避免譜半徑對(duì)收斂性的影響，本文選取了選主元的高斯消去法作為堆芯流場計(jì)算的全局求解方法。

4 結(jié)果與分析

4.1 計(jì)算結(jié)果

為驗(yàn)證全局求解方法和程序編寫的正確性，選取HTR-PM滿功率穩(wěn)態(tài)運(yùn)行狀態(tài)為計(jì)算對(duì)象，基本運(yùn)行參數(shù)如下：滿功率為250 MW；堆芯系統(tǒng)壓力為7×106Pa；堆芯氦氣入口溫度為250 ℃；堆芯氦氣入口流量為96 kg/s。分別用KONVEK程序和新開發(fā)的全局求解程序?qū)x取的問題進(jìn)行計(jì)算，圖4示出計(jì)算的堆芯壓力場分布情況。由圖4可見：新開發(fā)的全局求解程序計(jì)算出的壓力場分布與KONVEK程序計(jì)算結(jié)果吻合很好，只有極少數(shù)的網(wǎng)格點(diǎn)出現(xiàn)0.1%～0.5%的相對(duì)誤差，這是由于這些地方的相對(duì)壓力很小，幾乎接近零的緣故。

a——KONVEK程序；b——新開發(fā)的程序；c——計(jì)算結(jié)果的相對(duì)誤差

4.2 計(jì)算效率

圖5 收斂特性對(duì)比

為能夠充分體現(xiàn)全局求解方法的計(jì)算效率，選取HTR-PM瞬態(tài)運(yùn)行狀態(tài)作為計(jì)算對(duì)象，基本運(yùn)行參數(shù)如下：堆芯系統(tǒng)壓力為7×106Pa；堆芯氦氣入口溫度為250 ℃(恒定)；堆芯氦氣入口流量在0～1 h時(shí)為96 kg/s，1～1.5 h時(shí)從96 kg/s線性降至0.01 kg/s，1.5～3 h時(shí)保持0.01 kg/s不變。通過對(duì)比KONVEK程序和新開發(fā)的全局求解程序的收斂特性和計(jì)算時(shí)間，來分析新開發(fā)程序的計(jì)算效率。圖5示出HTR-PM瞬態(tài)運(yùn)行過程中1.407 h工況時(shí)間段內(nèi)流場的收斂特性曲線。在該工況時(shí)間段內(nèi)，溫度場與流場之間迭代了2次，每次流場計(jì)算的收斂特性如圖5所示。由圖5可知：當(dāng)速度場的收斂精度達(dá)到0.01附近后，速度場收斂很緩慢，幾乎不收斂，很難達(dá)到更高的精度；而新開發(fā)的程序僅需要2步迭代，速度場的收斂精度就可達(dá)到5×10-5。新開發(fā)的程序在收斂速率上有明顯優(yōu)勢(shì)。

將新開發(fā)的全局求解模塊應(yīng)用于THERMIX程序包，分析全局求解方法對(duì)堆芯流場和溫度場總計(jì)算效率的影響，表1列出堆芯流場和溫度場的最終收斂情況。由表1可知：原THERMIX程序包在計(jì)算流場時(shí)收斂速度慢，壓力場與速度場之間迭代900次仍不收斂，從而導(dǎo)致溫度場和流場之間迭代100次也不收斂，在表1所列每個(gè)工況時(shí)間段的計(jì)算效率均相對(duì)較低，耗時(shí)長；新開發(fā)的全局求解程序大幅提高了計(jì)算效率，流場和溫度場均在最大迭代次數(shù)之前很快達(dá)到收斂標(biāo)準(zhǔn)，計(jì)算效率提高了幾十倍，甚至上百倍，說明了該全局求解方法的高效性。

5 總結(jié)

對(duì)分塊迭代方法的收斂性分析研究表明：當(dāng)分塊迭代中迭代矩陣G的譜半徑ρr(G)接近1或ρr(G)≥1時(shí)，就會(huì)出現(xiàn)收斂慢，甚至不收斂的現(xiàn)象，并難以保證壓力場的求解精度。本文開發(fā)了適用于計(jì)算HTGR堆芯流場的高效全局求解方法，并對(duì)HTR-PM穩(wěn)態(tài)和瞬態(tài)運(yùn)行狀態(tài)進(jìn)行了驗(yàn)證計(jì)算，與分塊迭代方法相比，高效全局求解方法的收斂性和計(jì)算效率均有很大提高。

表1 計(jì)算效率對(duì)比

參考文獻(xiàn)：

[1] 吳宗鑫，張作義. 先進(jìn)核能系統(tǒng)和高溫氣冷堆[M]. 北京：清華大學(xué)出版社，2004：238-245.

[2] CLEVELAND J C, GREENE S R. Application of the THERMIX-KONVEK code to accident analyses of modular pebble bed high temperature reactors[R]. USA: Oak Ridge National Laboratory, 1986.

[3] 清華大學(xué)核能與新能源技術(shù)研究院. 華能山東石島灣核電廠高溫氣冷堆核電站示范工程初步安全分析報(bào)告[R]. 北京：清華大學(xué)核能與新能源技術(shù)研究院，2008.

[4] 王登營. 高溫氣冷堆多回路系統(tǒng)的模擬[D]. 北京：清華大學(xué)核能與新能源技術(shù)研究院，2011.

[5] 陶文銓. 數(shù)值傳熱學(xué)[M]. 2版. 西安：西安交通大學(xué)出版社，2001：198-203.

[6] HOOPER R, HOPKINS M, PAWLOWSKI R, et al. Final report on LDRD project: Coupling strategies for multi-physics applications, SAND2007-7146[R]. USA: Sandia National Laboratories, 2007.

[7] 關(guān)治，陸金甫. 數(shù)值分析基礎(chǔ)[M]. 2版. 北京：高等教育出版社，2010：35-66.