TOPAZⅡ反應(yīng)堆本體流固共軛傳熱數(shù)值模擬

鄒佳訊，郭春秋，孫 征

(中國(guó)原子能科學(xué)研究院 反應(yīng)堆工程技術(shù)研究部，北京 102413)

TOPAZ-Ⅱ是由俄羅斯開(kāi)發(fā)的熱離子型空間核動(dòng)力反應(yīng)堆的堆型[1,2]，TOPAZ-Ⅱ反應(yīng)堆不同于常規(guī)的陸上反應(yīng)堆，其堆型要求相對(duì)緊湊、小型或微型化。TOPAZ-Ⅱ堆內(nèi)的冷卻劑為液態(tài)NaK合金，NaK共晶合金具有熔點(diǎn)低，熱導(dǎo)率高等優(yōu)點(diǎn)。 國(guó)內(nèi)外針對(duì)TOPAZ-Ⅱ開(kāi)發(fā)了一系列專(zhuān)用系統(tǒng)分析程序，如俄羅斯開(kāi)發(fā)了ENSY系統(tǒng)分析程序?qū)?dòng)期間喪失冷卻劑事故進(jìn)行了分析，美國(guó)開(kāi)發(fā)了熱離子集成瞬態(tài)系統(tǒng)分析程序有TITAM，CENTRAR等[3～4]，國(guó)內(nèi)開(kāi)發(fā)的見(jiàn)諸報(bào)道的有西安交通大學(xué)開(kāi)發(fā)的瞬態(tài)空間熱離子反應(yīng)堆系統(tǒng)分析程序TASTIN等等[5～6]，這些系統(tǒng)分析程序可以分析穩(wěn)態(tài)或者瞬態(tài)工況下系統(tǒng)的整體參數(shù)及動(dòng)態(tài)響應(yīng)情況，關(guān)于TOPAZ-Ⅱ詳細(xì)的堆本體三維流場(chǎng)溫場(chǎng)方面的信息較少。堆芯入口流量分配數(shù)據(jù)是反應(yīng)堆熱工水力必不可少的計(jì)算輸入，從反應(yīng)堆力學(xué)分析的角度來(lái)看，堆芯燃料和部件可能會(huì)因?yàn)樽陨淼臏囟忍荻榷a(chǎn)生熱應(yīng)力，流動(dòng)和傳熱計(jì)算的到的溫場(chǎng)及流場(chǎng)可為力學(xué)分析提供設(shè)計(jì)輸入，一般情況下，固體的熱應(yīng)力分析和溫度計(jì)算分析可能是耦合的，但由于應(yīng)力形變相對(duì)于幾何尺寸有時(shí)候是可以忽略的，所以在確定熱應(yīng)力之前可單獨(dú)進(jìn)行溫度分布計(jì)算。目前的商用CFD軟件無(wú)論是在國(guó)外還是國(guó)內(nèi)，已越來(lái)越多應(yīng)用于諸如水冷堆、液態(tài)金屬冷卻反應(yīng)堆等反應(yīng)堆熱工水力行為的數(shù)值計(jì)算和研究中[7～9]。本文利用計(jì)算流體力學(xué)技術(shù)，針對(duì)TOPAZ Ⅱ反應(yīng)堆本體進(jìn)行流固共軛傳熱模擬，獲得關(guān)鍵的熱工水力參數(shù)及堆本體溫度，為后續(xù)設(shè)計(jì)提供參考。

1 計(jì)算模型

1.1 TOPAZ-Ⅱ系統(tǒng)描述

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

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

圖1 TOPAZ-Ⅱ反應(yīng)堆本體(停堆狀態(tài))Fig.1 Diagram of TOPAZ-Ⅱ reactor complex (shutdown condition)

表1 TOPAZ-Ⅱ的主要設(shè)計(jì)參數(shù)Table 1 Main design parameters of TOPAZ-Ⅱ

1.2 網(wǎng)格化分

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

圖2 計(jì)算域網(wǎng)格示意圖Fig.2 Mesh display of calculation zone

表2 數(shù)值結(jié)果網(wǎng)格獨(dú)立性Table 2 Numerical results independence of numbers of mesh

1.3 材料物性

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

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

堿性液態(tài)金屬的對(duì)流換熱數(shù)值模擬需要注意的是液態(tài)金屬非常低的普朗特?cái)?shù)Pr。

(1)

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

圖3 NaK合金普朗特?cái)?shù)隨溫度的變化擬合曲線Fig.3 Fitting curve of Pr of liquidNaK with temperature

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

表4 部分?jǐn)?shù)值結(jié)果對(duì)湍流模型的敏感性Table 4 Sensitivity of part results with certain turbulent models

后文CFD分析使用的是基于K-Omega的 SST模型，模型中的湍流普朗特?cái)?shù)Prt(渦擴(kuò)散系數(shù)和渦熱傳導(dǎo)系數(shù)的比值)默認(rèn)值為0.85(FLUENT)/0.9(CFX)，適用于輕水或空氣的模擬，不完全適用于模擬液態(tài)金屬換熱，有若干文獻(xiàn)中提出或報(bào)道了計(jì)算Prt的經(jīng)驗(yàn)關(guān)系式或推薦值，如下式：

(2)

Pe=Re×Pr

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

(3)

對(duì)于本文中的模擬對(duì)象TOPAZ-Ⅱ反應(yīng)堆液態(tài)NaK合金正常運(yùn)行工況其雷諾數(shù)Re>4 000，其流動(dòng)屬于環(huán)管內(nèi)的湍流流動(dòng)，且其Pe數(shù)約30。利用上述Reynolds式計(jì)算得到的Prt約為4.1和Xu Cheng關(guān)系式中的4.12基本一致，因此需要更改軟件的默認(rèn)值0.9為4.1。

1.4 計(jì)算條件

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

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

3 計(jì)算結(jié)果與分析

3.1 流量分配

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

(4)

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

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

3.2 摩擦因子

為了驗(yàn)證流動(dòng)數(shù)值計(jì)算的可靠性，從摩擦因子f的角度進(jìn)行分析，摩擦因子f反映壓降和流量的關(guān)系，其表達(dá)式如下：

(5)

其中ΔP/L為流動(dòng)方向長(zhǎng)度L上的壓降，De為等效水力直徑(通道特征尺寸)，ρ為流體密度，V為流體速度。對(duì)于光滑圓管湍流流動(dòng)，經(jīng)驗(yàn)關(guān)系式如下，其中Re為雷諾數(shù)：

(6)

德國(guó)學(xué)者V.Gnielinski[12]在光滑管湍流流動(dòng)摩擦因子關(guān)系式基礎(chǔ)針對(duì)環(huán)管湍流摩擦因子進(jìn)行了修訂，關(guān)系式如下：

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

(7)

(8)

(9)

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

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

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

3.3 努賽爾數(shù)

Nu是用于描述換熱的一個(gè)無(wú)量綱數(shù)[式(10)]，針對(duì)圓管/環(huán)管內(nèi)湍流流動(dòng)換熱，從二十世紀(jì)六七十年代至今，已有學(xué)者陸續(xù)提出計(jì)算液態(tài)金屬?gòu)?qiáng)迫對(duì)流換熱Nu若干經(jīng)驗(yàn)關(guān)系式， 其一般表達(dá)形式見(jiàn)式(11)或其修訂形式，在充分發(fā)展的流動(dòng)換熱流域中Nu一般隨著Pe數(shù)的增加而增加，在換熱穩(wěn)定段約為定值。

(10)

Nu=Nuo+a×Pe

(11)

式中，a,b分別為經(jīng)驗(yàn)關(guān)系式的擬合系數(shù)。

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

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

各流量臺(tái)階算例下的努賽爾數(shù)的對(duì)比見(jiàn)表5。從表中可以看出，三者之間的差別很小，表明了CFX數(shù)值模擬換熱分析的可靠性。

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

3.4 溫度場(chǎng)信息

環(huán)形通道內(nèi)冷卻劑溫度軸向分布、慢化劑、端部鈹反射層徑向沿某點(diǎn)沿著軸向上的溫度分布如圖8所，圖中縱坐標(biāo)為歸一化溫度(局部溫度/最高溫度)，橫坐標(biāo)為歸一化軸向位置(z/H)。冷卻劑在活性區(qū)高度375 mm距離上溫度成線性分布，在活性區(qū)外的兩端為溫度平臺(tái)，分別于出入口溫度接近，入口段溫度為743 K的設(shè)置值，出口溫度最高約846 K，僅與設(shè)計(jì)值843 K相差3 K。最高溫度出現(xiàn)在慢化劑區(qū)域，在中心燃料通道附近區(qū)域，最高溫度約為868 K，詳細(xì)的剖面溫度等值線見(jiàn)圖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 結(jié)論

本文利用數(shù)值模擬方法對(duì)TOPAZ-Ⅱ反應(yīng)堆進(jìn)行了流固共軛傳熱計(jì)算分析，通過(guò)計(jì)算分別對(duì)反應(yīng)流動(dòng)和壓降的摩擦因子以及反應(yīng)傳熱的努塞爾數(shù)進(jìn)行了比較，發(fā)現(xiàn)數(shù)值模擬計(jì)算得到的相應(yīng)結(jié)果與使用經(jīng)驗(yàn)關(guān)系式得到的計(jì)算結(jié)果非常接近，證明了數(shù)值模擬計(jì)算在反應(yīng)堆流動(dòng)與換熱分析上的可靠性與正確性，所得到的詳細(xì)流量分配數(shù)據(jù)可以為后續(xù)系統(tǒng)分析提供輸入，部件溫度三維場(chǎng)信息可以為力學(xué)計(jì)算提供接口，TOPAZ-Ⅱ的流固共軛傳熱數(shù)值模擬可以為后續(xù)相關(guān)堆型的優(yōu)化設(shè)計(jì)奠定基礎(chǔ)。