空間核電推進球床反應堆熱工水力特性數值分析

張 冉，王成龍,2,*，秋穗正,2，田文喜,2，蘇光輝,2

(1.西安交通大學 核科學與技術學院，陜西 西安 710049；2.西安交通大學 動力工程多相流國家重點實驗室，陜西 西安 710049)

隨著我國空間安全和空間探索領域的進一步拓展，航天器將面臨更為復雜、嚴峻的空間環境[1]。傳統航天器采用的“太陽電池陣+蓄電池組”的電源系統配置已無法滿足任務需求，環境適應性強、長壽命、安全可靠的空間核動力系統成為必然甚至是唯一的選擇。

球床反應堆(PeBR, pellet bed reactor)空間核動力系統具有功率密度大、結構緊湊等優點[2-3]，可廣泛應用于空間核反應堆電源供電、核電推進(NEP，nuclear electric propulsion)和核熱推進(NTP, nuclear thermal propulsion)，尤其是MW級空間核動力的核電推進具有廣泛的應用前景。

針對用于核熱推進的PeBR[4-7]，新墨西哥大學的空間核動力研究所進行了詳細的堆芯物理計算分析，并開發了一維、二維穩態熱工水力分析程序和二維瞬態熱工水力分析程序。針對用于核電推進的PeBR，除提出相關概念設計外，并未進行深入研究。本文針對PeBR的結構特點，開發穩態熱工水力分析程序，對其全功率運行下的穩態工況進行分析，并研究冷、熱孔板對流量分配和溫度分布的影響以及入口壓力對壓損的影響，為堆芯進一步熱工水力優化設計提供參考依據。

1 核電推進球床反應堆

用于NEP的PeBR是以氦氣為冷卻劑的快中子堆[8-11]，圖1示出PeBR的軸向、徑向截面示意圖。如圖1b所示，堆芯是由內、外兩層孔板圍成的1個環形區域，球形燃料元件在區域內隨意堆積。孔板上開有大小不同的圓形小孔，孔徑通常較燃料元件直徑小10%～20%，即8～9 mm；孔板的平均孔隙率(孔板上所有小孔的面積占孔板表面積的比值)為10%～30%，內、外孔板也根據其溫度高低，稱為熱孔板和冷孔板。如圖1b所示，反應堆堆芯被隔板分為3個120°的扇形區域，每個區域均可獨立運行和被冷卻。每個區域與1套CBC系統相連，每個回路可產生5 MW電能。

從壓縮機出來的氦氣以約800 K的溫度進入反應堆，流經路徑如下：1) 氦氣從頂部反射層冷卻通道流過，同時冷卻頂部反射層；2) 氦氣向下流動進入冷卻劑入口通道，通過冷孔板流入堆芯,同時冷卻徑向反射層和控制鼓；3) 堆芯的氦氣與燃料球進行換熱后，通過熱孔板流出堆芯；4) 一小部分氦氣通過冷卻劑漏流管進入底部反射層冷卻通道，冷卻底部反射層后進入冷卻劑出口通道；5) 冷卻劑出口通道的氦氣流出反應堆去往渦輪機。系統及堆芯主要設計參數列于表1。

圖1 PeBR的軸向(a)和徑向(b)截面

圖2示出堆芯物理計算得出的徑向、軸向功率分布，歸一化功率密度定義為局部功率密度與平均功率的比值。將此功率分布應用于后面的熱工水力特性分析。

表1 系統及堆芯主要設計參數

圖2 堆芯的徑向(a)和軸向(b)功率分布

2 數學物理模型

圖3 PeBR的軸對稱物理模型

考慮堆芯結構的對稱性，以冷卻劑出口通道中心為對稱軸，建立PeBR的軸對稱物理模型，如圖3所示，z為堆芯軸向坐標，r為堆芯徑向坐標，H為堆芯的高度。模型中包含了反射層冷卻通道、冷卻劑入口通道、冷孔板、堆芯和熱孔板5個部分。針對氦氣在堆芯內的流動特性，作如下假設： 1) 氦氣在堆芯內的流動為沿徑向的一維流動；2) 不考慮堆芯的軸向導熱；3) 不考慮由黏性耗散引起的能量變化；4) 由于反應堆在太空失重環境中運行，不考慮重力。

2.1 堆芯模型

根據假設條件，堆芯內的控制方程如下。

氦氣質量守恒方程：

(1)

氦氣動量守恒方程：

(2)

氦氣能量守恒方程：

(3)

動量守恒方程中的球床阻力因子f采用von Krm?n關系式(式(4)～(6))計算[10]：

(4)

(5)

(6)

能量守恒方程中的asv為單位體積換熱面積：

(7)

式中：A為流通面積；W為氦氣質量流量；ρ為氦氣密度；p為氦氣壓力；dp為燃料球直徑；cp為氦氣比定壓熱容；u為氦氣流速；λv為氦氣熱導率；Tv為氦氣溫度；Ts為燃料球表面溫度；εcore為球床孔隙率，本文采用0.38，未考慮孔隙率在近壁面處的變化；φ為分段關系式；Rem為修正雷諾數；μ為氦氣動力黏度；hsv為氦氣與球床的對流換熱系數，可根據努塞爾數求得，本文采用式(8)計算：

(8)

(9)

式中，q?core為堆芯功率密度。

(10)

(11)

式中：λp為燃料球的有效熱導率[12]；λt為TRISO顆粒的等效熱導率；λc為石墨基質的熱導率；φ為燃料球的填充因子。

2.2 球床等效導熱模型

PeBR球床內充滿燃料球，球床的熱量除通過氦氣的對流換熱傳遞外，還有3種方式傳遞：燃料球通過球床空隙的輻射換熱、球床空隙內的氣體導熱和燃料球之間的接觸導熱。

球床能量守恒方程為：

(12)

式中，λeff為球床的有效導熱系數，本文使用以下3種有效導熱系數模型。

(13)

(14)

(15)

式中：λfuel為燃料球熱導率；B為變形因子；Λ為等效熱導率；εr為燃料球發射率；μp為橫向變形系數；Es為楊氏模量；NA為單位面積球個數；NL為單位長度球個數；p′為外部壓力；σ為黑體輻射常數；S和SF為結構常數。

2.3 冷、熱孔板

冷、熱孔板均當作單級多孔孔板處理，可采用下式[14]計算：

(16)

(17)

(18)

式中：ξ為流阻系數；S0為節流孔流通面積；S1為管道流通面積；τ為孔板外形調整系數，是孔板厚度與孔徑之比e/d0的函數；λ為孔內沿程阻力系數；v0為孔內流速；Δp為氦氣通過孔板的壓降；N為多孔孔板的孔數；d1為多孔孔板的孔徑。

3 數值方法

基于以上數學物理模型，根據PeBR的結構及運行特點開發了穩態熱工水力分析程序。程序使用有限容積法對方程進行求解，依照氦氣在堆芯內的流動特點，將堆芯先沿軸向均勻劃分為10個獨立氦氣通道，再將這10個通道沿徑向劃分成26個控制體。圖4示出PeBR的控制體劃分。

圖4 PeBR的控制體劃分

穩態下，依據并聯多通道模型對堆芯進行流量分配。程序計算流程如圖5所示：1) 輸入系統、結構等參數，根據堆芯參數初始化變量；2) 假設堆芯初始流量分配；3) 假設堆芯初始導熱換熱量分布；4) 計算各通道氦氣焓場、溫度場和壓力場；5) 計算球床的溫度分布；6) 利用有效導熱系數計算球床的導熱換熱量；7) 如果導熱換熱量與初始假設值大于誤差，則返回步驟3重新計算，直到滿足能量守恒為止；8) 如果各通道壓降不相等，則返回步驟2修改所假設的堆芯流量分配數據，重復上述計算過程，直到滿足壓降相等為止。

圖5 PeBR穩態熱工水力計算流程

4 結果與分析

4.1 孔隙率優化分析

PeBR的徑向和軸向功率分布不均，且堆芯裝載的是具有相同富集度的燃料球，無法通過調整燃料富集度來展平功率。冷、熱孔板不僅用來固定球床，還起到調節堆芯流量分配的作用，如果冷卻劑流量分配不合適，則堆芯內極有可能產生局部熱點，因此優化孔板孔隙率是一研究重點。

1) 孔隙率均勻分布

基于NTP PeBR的設計，冷孔板選擇孔隙率為0.1的軸向均勻分布，熱孔板選擇孔隙率為0.3的均勻分布[5]。圖6示出氦氣質量流量沿軸向的分布，軸向位置用相對高度(z/H)表示。可看出，氦氣在堆芯頂部、底部流量大，在堆芯中部流量小，與堆芯軸向功率密度分布恰恰相反。這是由于堆芯中部功率密度大，流過中部的氦氣溫度更高、比容更大，從而流速和壓降均增大，為了使各通道的壓降相等，堆芯流量重新分配時，就降低了堆芯中部的氦氣流量。

圖6 氦氣質量流量沿軸向的分布

圖7示出氦氣、燃料球表面和燃料球中心溫度在熱孔板處沿軸向的分布。由圖7可見：氦氣平均出口溫度為1 974 K，接近設計值2 000 K，軸向溫差648 K；燃料球中心最高溫度為2 500 K，雖低于燃料的熔點((3 123±50) K)，但軸向溫差高達732 K，堆芯中心出現熱點，且不均勻的溫度分布不利于反應堆熱工水力特性。圖8示出堆芯內氦氣溫度、壓力和燃料球中心溫度的分布，可明顯看出堆芯中心出現熱點。

圖7 熱孔板處溫度沿軸向的分布

a——堆芯內氦氣溫度；b——燃料球中心溫度；c——堆芯內氦氣壓力圖8 孔板孔隙率均勻分布的計算結果

2) 冷孔板孔隙率分析

為消除堆芯熱點，使堆芯溫度分布更均勻，希望孔板孔隙率沿軸向的分布與軸向功率密度分布相一致，使堆芯在低功率密度處的質量流量小，堆芯中心處有更大的質量流量。為研究冷孔板孔隙率對堆芯流量分配的影響，定義相對孔隙率為局部孔隙率與平均孔隙率的比值，孔隙率峰值比 (PPR, porosity peak ratio)為最大孔隙率與平均孔隙率的比值。

熱孔板采用孔隙率0.3的均勻分布。本文給出了5種冷孔板孔隙率分布，如圖9所示，PPR分別是1.5、1.9、2.0和2.2。PPR為1.0即孔隙率均勻分布。由圖9可見，PPR從1.5增大至2.2，冷孔板頂部、底部孔隙率明顯減小，中部孔隙率略微減小。

圖10示出冷孔板PPR對堆芯質量流量的影響，圖11、12分別示出冷孔板PPR對氦氣出口溫度和燃料球中心溫度的影響。由圖10～12可知，調節冷孔板的孔隙率，對堆芯的質量流量和溫度分布均沒有很大改善。這是因為氦氣在進氣通道內的溫度為870～950 K，經過冷孔板的氦氣沿軸向溫度變化不大，比容小，所以冷孔板孔隙率分布對堆芯質量流量分配影響不大。綜上，冷孔板孔隙率的設計參考值取為0.1。

圖9 不同PPR的冷孔板孔隙率沿軸向的分布

圖10 冷孔板PPR對質量流量的影響

圖11 冷孔板PPR對氦氣出口溫度的影響

3) 熱孔板孔隙率分析

為研究熱孔板孔隙率對堆芯流動和溫度的影響，同樣給出冷孔板采用孔隙率為0.1的均勻分布時，5種PPR的熱孔板孔隙率分布，如圖13所示。圖14示出熱孔板PPR對堆芯質量流量的影響，可看出，與孔隙率均勻分布相比，4種PPR的孔隙率均有效改善了堆芯的流量分布，隨PPR的增大，堆芯頂部、底部的流量越來越小，中部的流量越來越大。

圖12 冷孔板PPR對燃料球中心溫度的影響

圖13 不同PPR的熱孔板孔隙率沿軸向的分布

圖14 熱孔板PPR對質量流量的影響

圖15、16分別示出熱孔板PPR對氦氣出口溫度和燃料球中心溫度分布的影響。由圖16可知， 隨PPR從1.0增大到2.2，堆芯頂部、底部的溫度越來越高，中部的溫度越來越低。當PPR小于2.0時，堆芯中部溫度高于頂部和底部溫度，軸向溫差大；PPR為2.0時，軸向溫差最小、堆芯溫度分布最均勻；當PPR大于2.0時，出現堆芯中部溫度低、頂部和底部溫度高的情況，顯然這樣的溫度分布不合理。表2列出不同PPR的熱孔板處的溫度參數，經綜合比較，選定PPR= 2.0為設計參考值。

圖15 熱孔板PPR對氦氣出口溫度的影響

圖16 熱孔板PPR對燃料球中心溫度的影響

表2 不同PPR的熱孔板處溫度參數

4.2 最優值分析

圖17 氦氣溫度分布

圖18 燃料球中心溫度分布

冷孔板孔隙率為0.1、熱孔板PPR為2.0時，圖17、18分別示出氦氣和燃料球中心的溫度分布。氦氣的平均出口溫度為2 003 K，軸向溫差為10 K。燃料球中心最高溫度為2 250 K，軸向溫差為81 K。與孔隙率沿軸向均勻分布相比，堆芯溫度沿軸向分布更加均勻。燃料球中心最高溫度距ZrC包殼的失效溫度3 200 K還有很大裕量(950 K)，ZrC包殼作為TRISO顆粒球的最后一道安全屏障，可完全包容所有的裂變產物。堆芯球床的徑向溫度梯度最大僅為3.5 K/mm，燃料球的局部熱點和熱應力可忽略不計。

圖19示出氦氣的壓力分布。圖20示出不同軸向位置氦氣壓力沿徑向的分布。可看出，越靠近熱孔板，氦氣壓力降低得越快，這是由于隨氦氣向熱孔板的流動，流通截面減小，使得氦氣流速增加，從而壓損增大。同時，流向熱孔板的過程中，氦氣的溫度升高，使得氦氣的比容增大，導致更大的流速和壓損。氦氣流過熱孔板時，壓力有顯著的變化，體現了熱孔板孔隙率對堆芯流量分配的有效調節作用。

圖19 氦氣壓力分布

圖20 不同軸向位置的徑向壓力分布

孔隙率均勻分布時，氦氣流過反應堆的總壓損(包括反射層冷卻通道、進氣通道、冷孔板、堆芯和熱孔板)為0.132 MPa，堆芯壓損(包括冷、熱孔板和堆芯)為0.106 MPa。孔隙率優化后，總壓損為0.172 MPa，增加了30%，堆芯壓損為0.144 MPa，增加了36%。雖然壓損有所增大，但由于氦氣在堆芯內為徑向流動，反應堆總壓損仍很小，僅為入口壓力的17%。

4.3 入口壓力對氦氣流動的影響

PeBR中氦氣沿徑向流動，大幅減小了堆芯內壓損，因此反應堆可在較低的系統壓力下運行，一般在1～5 MPa之間。表3列出不同氦氣入口壓力下的壓損。

表3 入口壓力對壓損的影響

由表3可知，隨入口壓力的減小，總壓損和堆芯壓損均越來越大。但系統壓力降低不僅有利于減小壓氣機的耗功量，提高能量轉化效率，且可減輕管道、壓力容器等承壓設備質量，從而降低系統的尺寸和質量。

5 結論

針對應用于NEP的PeBR開發了穩態熱工水力程序，并對PeBR進行了熱工安全分析和優化設計，所得的主要結論如下。

1) 反應堆全功率穩態運行工況下，燃料球中心最高溫度遠低于燃料的熔點，具有873 K的安全裕量。球床的徑向溫度梯度為3.5 K/mm，燃料球的局部熱點和熱應力可忽略不計。

2) 冷孔板孔隙率分布對調節堆芯流量分布作用不大，設計中可選擇孔隙率為0.1的均勻分布。熱孔板孔隙率分布對堆芯流量分配有很大影響，PPR為2.0的熱孔板不僅可避免堆芯熱點，且可使氦氣溫度、堆芯溫度沿軸向分布更均勻，提高了反應堆的熱工水力特性。

3) 由于熱孔板處流量限制，與孔隙率均勻分布時相比，優化后的熱孔板將使壓損增大，但由于氦氣在堆芯內為徑向流動，反應堆總壓損仍很小。隨氦氣入口壓力的增大，反應堆總壓損和堆芯壓損均減小，但較大的系統壓力會相應增加系統的尺寸和質量。