TOPAZⅡ反應堆本體流固共軛傳熱數值模擬

鄒佳訊，郭春秋，孫 征

(中國原子能科學研究院 反應堆工程技術研究部，北京 102413)

TOPAZ-Ⅱ是由俄羅斯開發的熱離子型空間核動力反應堆的堆型[1,2]，TOPAZ-Ⅱ反應堆不同于常規的陸上反應堆，其堆型要求相對緊湊、小型或微型化。TOPAZ-Ⅱ堆內的冷卻劑為液態NaK合金，NaK共晶合金具有熔點低，熱導率高等優點。 國內外針對TOPAZ-Ⅱ開發了一系列專用系統分析程序，如俄羅斯開發了ENSY系統分析程序對啟動期間喪失冷卻劑事故進行了分析，美國開發了熱離子集成瞬態系統分析程序有TITAM，CENTRAR等[3～4]，國內開發的見諸報道的有西安交通大學開發的瞬態空間熱離子反應堆系統分析程序TASTIN等等[5～6]，這些系統分析程序可以分析穩態或者瞬態工況下系統的整體參數及動態響應情況，關于TOPAZ-Ⅱ詳細的堆本體三維流場溫場方面的信息較少。堆芯入口流量分配數據是反應堆熱工水力必不可少的計算輸入，從反應堆力學分析的角度來看，堆芯燃料和部件可能會因為自身的溫度梯度而產生熱應力，流動和傳熱計算的到的溫場及流場可為力學分析提供設計輸入，一般情況下，固體的熱應力分析和溫度計算分析可能是耦合的，但由于應力形變相對于幾何尺寸有時候是可以忽略的，所以在確定熱應力之前可單獨進行溫度分布計算。目前的商用CFD軟件無論是在國外還是國內，已越來越多應用于諸如水冷堆、液態金屬冷卻反應堆等反應堆熱工水力行為的數值計算和研究中[7～9]。本文利用計算流體力學技術，針對TOPAZ Ⅱ反應堆本體進行流固共軛傳熱模擬，獲得關鍵的熱工水力參數及堆本體溫度，為后續設計提供參考。

1 計算模型

1.1 TOPAZ-Ⅱ系統描述

該系統[10]由堆本體、電磁泵、輻射屏蔽(由鋼重屏蔽及氫化鋰輕屏蔽組成)、熱排放系統、穩壓系統以及重要的儀表控制、管道閥門等組成，如圖1所示。其堆本體活性區內有氫化鋯慢化劑，慢化劑上下為端部鈹反射層，慢化劑外層為容器筒體；筒體外側為側鈹反射層，側鈹反射層內有12個控制鼓其中3個為安全鼓，另外9個為調節轉鼓，鼓體材料為鈹，吸收體材料為碳化硼；活性區有 37個熱離子燃料元件可以產生共計約4.5～5.5 kW電源，37個熱離子元件分別坐落于均布在的慢化劑的37個豎直孔道內，孔道內由不銹鋼內外套管組成單獨的環形窄封流體通道，液態NaK合金經過37個環形窄封孔道豎直流動，不銹鋼內外套管之間的間隙為冷卻劑孔道，其中不銹鋼內套管外徑為24.5 mm,不銹鋼外套管內徑為25.9 mm，內外套管之間的間隙僅為0.7 mm。

37個流體通道兩端為上下集流腔室，其中下集流腔上柵板、下集流腔下柵板與反應堆筒體焊在一起后形成冷卻劑下集流腔；上集流腔上柵板、上集流腔下柵板焊在一起后形成冷卻劑上集流腔。在慢化劑孔道內的冷卻劑內外套管之間形成的流道與冷卻劑上下集流腔相通，其中充滿NaK冷卻劑。圖1(b)中1至37為燃料元件冷卻劑通道所在位置編號。系統主要設計參數[1-6]如表1所示。

圖1 TOPAZ-Ⅱ反應堆本體(停堆狀態)Fig.1 Diagram of TOPAZ-Ⅱ reactor complex (shutdown condition)

表1 TOPAZ-Ⅱ的主要設計參數Table 1 Main design parameters of TOPAZ-Ⅱ

1.2 網格化分

以圖2所示的TOPAZ-Ⅱ反應堆轉換器為對象，可以采用分區的方法進行網格化法，上下上腔室分別為一個區域，中間37個冷卻通道為第三個區域，中間的冷卻通道區域兩段分別用交界面與上腔室下柵格板、下腔室上柵格板相連。冷卻劑上下腔室區域生成7層附面層貼壁網格。冷卻劑周圍的主要固體區域之間在區域設置上相互獨立，以便靈活地的進行網格控制，固體區域在相貼的壁面上用交界面進行連接，用于數值模擬計算時熱傳導傳遞數據。計算結果與網格數目的敏感性方面，共劃分了四套網格(見表2)進行計算，選擇流體域NaK合金的最高溫度和固體域慢化劑塊的最高溫度進行比較，從表中可以看出，網格總數在300萬～400萬的三套網格計算得到NaK合金最高溫度、慢化劑最高溫度之間的差距不到1 ℃，滿足網格獨立性的要求，后續給出的分析均為基于第四套網格的計算結果，第四套網格中具體的網格數目組成包括：流體域NaK合金2 802 616單元數；二氧化碳修補氣體區域319 088單元數；慢化劑及端部鈹反射層區域361 767單元數；側鈹反射層區域301 203單元數；碳化硼吸收體區域54 096單元數；控制鼓體區域180 516單元數。

圖2 計算域網格示意圖Fig.2 Mesh display of calculation zone

表2 數值結果網格獨立性Table 2 Numerical results independence of numbers of mesh

1.3 材料物性

模擬計算中主要涉及的流體為NaK合金，固體有氫化鋯慢化劑、鈹反射層、不銹鋼等，穩態數值模擬計算中需要的關鍵熱物性熱中，三者的密度分別為5 615 kg/m3、1 830 kg/m3和7 900 kg/m3，熱導率分別為22.6 W/(m·K)、94.3 W/(m·K)、22.2W/(m·K)。流體換熱中需要用到NaK合金的密度、比熱容、熱導率、動力黏度等熱物性，固體的密度、比熱容和熱導率。冷卻劑熱物性參數(來自NaK工程手冊)包括密度、熱導率、比定壓熱容、動力黏度(見表3)，可以采取隨溫度的變化表進行插值或者用擬合函數載入求解。

表3 冷卻劑熱物性Table 3 Thermal and Transport Properties of Coolant NaK

堿性液態金屬的對流換熱數值模擬需要注意的是液態金屬非常低的普朗特數Pr。

(1)

從表2可以得到NaK的Pr隨著溫度的變化曲線及擬合函數如圖3所示，在高溫區750～950 K范圍內，Pr在0.006左右，非常低的Pr意味著分子熱傳導在整個傳熱過程中占據著較大的比重。

圖3 NaK合金普朗特數隨溫度的變化擬合曲線Fig.3 Fitting curve of Pr of liquidNaK with temperature

CFX前處理器提供了多種湍流模型，如常見的標準k-ε模型、k-ω、SSTk-ω等雷諾平均納維斯托克斯渦粘模型(RANS)，CFX前處理器中用自動近壁面處理方法(即壁面函數)對湍流流動中的近壁面的流動進行預測，而不需要精細化貼壁網格，在無量綱距離y+<300情況下，壁面函數均有效，且對y+沒有最小值要求。表4為額定流量下部分數值結果針對湍流模型的敏感性，從中可以看出，三種常見的湍流模型得到的溫度值相互之間的差距僅為0.2，但在反應流動和換熱的特性數據來看，SST總體上更為接近經驗關系式得到的值(見后文表5、表6)。

表4 部分數值結果對湍流模型的敏感性Table 4 Sensitivity of part results with certain turbulent models

后文CFD分析使用的是基于K-Omega的 SST模型，模型中的湍流普朗特數Prt(渦擴散系數和渦熱傳導系數的比值)默認值為0.85(FLUENT)/0.9(CFX)，適用于輕水或空氣的模擬，不完全適用于模擬液態金屬換熱，有若干文獻中提出或報道了計算Prt的經驗關系式或推薦值，如下式：

(2)

Pe=Re×Pr

Prt=4.12，Pe<1 000(Xu Cheng[11])

(3)

對于本文中的模擬對象TOPAZ-Ⅱ反應堆液態NaK合金正常運行工況其雷諾數Re>4 000，其流動屬于環管內的湍流流動，且其Pe數約30。利用上述Reynolds式計算得到的Prt約為4.1和Xu Cheng關系式中的4.12基本一致，因此需要更改軟件的默認值0.9為4.1。

1.4 計算條件

上腔室入口采取質量流量入口邊界，正常穩態額定工況總入口流量為1.3 kg/s，下腔室出口采取壓力出口邊界，設定為165000Pa，CFX中流體域中開啟K-Omega SST湍流模型。物性數據用CFX的表達語言載入相應模塊?；钚詤^375 mm高度上的環形管道內壁面使用均勻熱流密度，慢化劑釋熱率在徑向上分為5個區域，每個區域的軸向歸一化釋熱率分布如下圖4所示，其中縱坐標為歸一化釋熱率(針對慢化劑區域最大釋熱率進行歸一化，慢化劑最大體積釋熱率約1.0×106W/m3)，分別對各區域的軸向分布進行曲線擬合，得到與軸向位置有關的釋熱分布，然后在前處理中使用表達式語言載入，而端部鈹反射層(約0.1 kW)、側鈹反射層(約0.2 kW)、碳化硼(約0.4 kW)等其他固體區域發熱相對于慢化劑(約4 kW)來說比較小，因此使用各區域平均體積釋熱率作為輸入。側反射層外表面設置輻射換熱，發射系數為0.7，環境溫度設置為300 K。

圖4 慢化劑釋熱率分布Fig.4 Normalized heat generation rate distribution in moderator

3 計算結果與分析

3.1 流量分配

數值模擬計算完畢后，從CFD-POST提取37個燃料元件冷卻通道的流量數據，通道冷卻劑入口的流量分配系數由下式計算：

(4)

其中，εi、Mi分別為i通道的流量分配系數和流量。數值模擬得到流量分配結果如下圖5所示。圖中橫坐標上的數字1～37代表各通道編號，編號與圖1(b)中的號碼一一對應，縱坐標為流量分配系數，由式(4)計算得到。盡管冷卻劑從入口管進入上腔室這段流程流場較復雜，但由于37個燃料冷卻通道首先入口流通面積相等，各入口面基本處于等壓面，因此經分配后，各冷卻劑管道內的冷卻劑達到一相對均勻的狀態。從圖中可以看出，各通道流量分配因子在0.99～1.01；從堆芯中心至堆芯外圍區域入口流量分配系數逐漸變大，各圈流量分配系數分布較為均勻。εi>1.0的燃料元件主要分布在堆芯外圍區域，出口附近的幾個通道(編號為20、21、36、37；30、31、32、33；24、25、26、27)因子較大。

圖5 通道流量分配因子Fig.5 Mass-flow rate distribution factors

3.2 摩擦因子

為了驗證流動數值計算的可靠性，從摩擦因子f的角度進行分析，摩擦因子f反映壓降和流量的關系，其表達式如下：

(5)

其中ΔP/L為流動方向長度L上的壓降，De為等效水力直徑(通道特征尺寸)，ρ為流體密度，V為流體速度。對于光滑圓管湍流流動，經驗關系式如下，其中Re為雷諾數：

(6)

德國學者V.Gnielinski[12]在光滑管湍流流動摩擦因子關系式基礎針對環管湍流摩擦因子進行了修訂，關系式如下：

fg=(1.8*lg10(Re*)-1.5)-2

(7)

(8)

(9)

式(8)即Gnienlinski增加的修訂部分，α為環管內徑di和外徑do的比值，

為了對數值模擬得到的環管摩擦因子進行比較，在CFD-POST中的利用式(5)求得不同流量算例下環管的摩擦因子f_CFD，分別和利用經驗關系式(6)和式(8)得到曲線進行比較，如圖6所示，圖中橫坐標為額定流量65%～150%范圍內雷諾數Re，縱坐標表對應的摩擦系數。從圖中可以看出，數值模擬得到的摩擦系數與關系式(8)得到的摩擦系數fg之間的整體偏差不到5%，驗證了流動計算的可靠性。

圖6 不同雷諾數下摩擦系數比較Fig.6 Comparision of frcition factors under different Re number

3.3 努賽爾數

Nu是用于描述換熱的一個無量綱數[式(10)]，針對圓管/環管內湍流流動換熱，從二十世紀六七十年代至今，已有學者陸續提出計算液態金屬強迫對流換熱Nu若干經驗關系式， 其一般表達形式見式(11)或其修訂形式，在充分發展的流動換熱流域中Nu一般隨著Pe數的增加而增加，在換熱穩定段約為定值。

(10)

Nu=Nuo+a×Pe

(11)

式中，a,b分別為經驗關系式的擬合系數。

在CFD-POST中利用式(10)求得發熱段入口一段距離上Nu的分布如圖7所示，從圖中看出其換熱穩定段的Nu約為6.1。

圖7 額定工況下Nu沿受熱段軸向的分布Fig.7 Axial distribution of Nu number in inlet regions of NAK channel under nominal condition

各流量臺階算例下的努賽爾數的對比見表5。從表中可以看出，三者之間的差別很小，表明了CFX數值模擬換熱分析的可靠性。

表5 各流量下努賽爾數Table 5 Nu under different mass flow rate conditions

3.4 溫度場信息

環形通道內冷卻劑溫度軸向分布、慢化劑、端部鈹反射層徑向沿某點沿著軸向上的溫度分布如圖8所，圖中縱坐標為歸一化溫度(局部溫度/最高溫度)，橫坐標為歸一化軸向位置(z/H)。冷卻劑在活性區高度375 mm距離上溫度成線性分布，在活性區外的兩端為溫度平臺，分別于出入口溫度接近，入口段溫度為743 K的設置值，出口溫度最高約846 K，僅與設計值843 K相差3 K。最高溫度出現在慢化劑區域，在中心燃料通道附近區域，最高溫度約為868 K，詳細的剖面溫度等值線見圖9，10。

圖8 慢化劑、端部鈹反射層以及冷卻劑軸向溫度分布Fig.8 Axial temperature distribution of moderator、end reflectors and coolant

圖9 橫截面溫度示意圖Fig.9 Temperature distribution at axial cross section

圖10 縱剖面溫度示意圖Fig.10 Temperature contour in the longitudinal cross section

4 結論

本文利用數值模擬方法對TOPAZ-Ⅱ反應堆進行了流固共軛傳熱計算分析，通過計算分別對反應流動和壓降的摩擦因子以及反應傳熱的努塞爾數進行了比較，發現數值模擬計算得到的相應結果與使用經驗關系式得到的計算結果非常接近，證明了數值模擬計算在反應堆流動與換熱分析上的可靠性與正確性，所得到的詳細流量分配數據可以為后續系統分析提供輸入，部件溫度三維場信息可以為力學計算提供接口，TOPAZ-Ⅱ的流固共軛傳熱數值模擬可以為后續相關堆型的優化設計奠定基礎。