蘇興康，顧 龍,3,*，李顯文，張 璐，盛 鑫

(1.中國科學(xué)院 近代物理研究所，甘肅 蘭州 730000；2.中國科學(xué)院大學(xué) 核科學(xué)與技術(shù)學(xué)院，北京 100049；3.蘭州大學(xué) 核科學(xué)與技術(shù)學(xué)院，甘肅 蘭州 730000)

液態(tài)金屬，如鈉、鉛及鉛鉍共晶合金等，具有良好的流動(dòng)傳熱特性，被選為許多先進(jìn)能源系統(tǒng)的候選冷卻劑，如第四代堆[1]、次臨界堆[2]和太陽能發(fā)電站[3]等。然而大多數(shù)液態(tài)金屬屬于低普朗特?cái)?shù)(Pr)流體，如鉛鉍和鈉的Pr分別在0.03～0.01、0.01～0.006之間[4]。液態(tài)金屬的速度邊界層與溫度邊界層具有較大的差異性。在雷諾時(shí)均(RANS)方法中，通常使用Boussinesq假設(shè)將雷諾應(yīng)力轉(zhuǎn)化為湍流運(yùn)動(dòng)黏度的計(jì)算，而采用簡單梯度擴(kuò)散假設(shè)(SGDH)將雷諾熱流轉(zhuǎn)化為湍流熱擴(kuò)散系數(shù)的計(jì)算。對于傳統(tǒng)流體(Pr≈1)，雷諾比擬假設(shè)湍流運(yùn)動(dòng)黏度和湍流熱擴(kuò)散率的比值，即湍流普朗特?cái)?shù)(Prt)，可近似為常數(shù)(0.85～0.9)；而對于液態(tài)金屬(Pr?1)，其湍流普朗特?cái)?shù)分布與普朗特?cái)?shù)、流動(dòng)幾何和流動(dòng)狀態(tài)呈非線性相關(guān)，因此雷諾比擬假設(shè)定律不適用于其傳熱計(jì)算[5]。

為提高液態(tài)金屬傳熱的計(jì)算精度，在RANS方法中，可考慮不采用SGDH的方法來計(jì)算雷諾熱流，而是對雷諾熱流使用微分熱通量模型(DHFM)或代數(shù)熱通量模型(AHFM)來模擬。DHFM和AHFM對雷諾熱流矢量直接建立二階矩的多個(gè)微分輸運(yùn)方程或一階半的代數(shù)輸運(yùn)方程，來模擬液態(tài)金屬的速度場與溫度場之間的差異性[6]。許多學(xué)者對DHFM和AHFM的輸運(yùn)方程以及模型系數(shù)進(jìn)行了研究[7-10]。由于這些模型的建立方法多樣、輸運(yùn)方程數(shù)量較多和系數(shù)敏感性較強(qiáng)，尚未得到推廣。

為改善SGDH方法的適用性，部分學(xué)者基于全局流動(dòng)參數(shù)或局部湍流參數(shù)建立了用于計(jì)算液態(tài)金屬的非線性Prt關(guān)系式，如Cheng等[11]引入全局雷諾數(shù)Re和普朗特?cái)?shù)Pr效應(yīng)建立了全局Prt函數(shù)關(guān)系，Kays等[12]引入湍流運(yùn)動(dòng)黏度和普朗特?cái)?shù)Pr建立了局部Prt函數(shù)關(guān)系。Duponcheel等[13]分析了現(xiàn)有的多個(gè)Prt關(guān)系式在簡單幾何中的計(jì)算精度，認(rèn)為基于局部湍流參數(shù)建立的Prt關(guān)系式能較好地提高液體金屬的計(jì)算精度，但在復(fù)雜幾何中的適用性還需進(jìn)一步分析。

為進(jìn)一步提高SGDH方法的普適性，許多學(xué)者將計(jì)算雷諾熱流轉(zhuǎn)化為計(jì)算湍流熱擴(kuò)散系數(shù)后，類比使用兩方程k-ε湍流模型計(jì)算湍流運(yùn)動(dòng)黏度的方法，將湍流時(shí)間尺度引入湍流熱擴(kuò)散系數(shù)，從而建立兩方程kθ-εθ傳熱模型來計(jì)算湍流熱擴(kuò)散系數(shù)。用于封閉動(dòng)量方程的兩方程k-ε湍流模型和用于封閉能量方程的兩方程kθ-εθ傳熱模型統(tǒng)稱為四方程k-ε-kθ-εθ模型[14]。近年來，Manservisi等[15-17]基于Nagano等[18-19]和Abe等[20]建立的四方程模型，引入全局湍流普朗特?cái)?shù)0.9，發(fā)展了適用于液態(tài)鉛鉍湍流換熱的四方程模型，較好地預(yù)測了液態(tài)鉛鉍(Pr=0.025)在平行平板、圓柱形管道和棒束通道內(nèi)的充分發(fā)展換熱。何少鵬等[21]基于Manservisi等的四方程模型和開源計(jì)算流體力學(xué)程序OpenFOAM，計(jì)算了液態(tài)鉛鉍在帶繞絲棒束內(nèi)的傳熱特性。

1 湍流換熱模型

基于Boussinesq假設(shè)和SGDH的周期性充分發(fā)展RANS方程組如下：

(1)

(2)

uz4qw/ρcpubDh

(3)

式中：θ為周期性溫度；ub為流體平均速度；Dh為當(dāng)量直徑；qw為施加在棒表面的恒熱流密度；ρ、ν、α和cp分別為流體的物性參數(shù)密度、分子運(yùn)動(dòng)黏度、分子熱擴(kuò)散系數(shù)和比熱容；νt和αt分別為需要封閉計(jì)算的湍流運(yùn)動(dòng)黏度和湍流熱擴(kuò)散系數(shù)。相應(yīng)的周期性RANS方程分析可參考文獻(xiàn)[27-28]。

引入湍動(dòng)能k及其耗散率ε、溫度波動(dòng)kθ及其耗散率εθ來考慮νt和αt之間的差異性。采用泰勒級數(shù)展開法可獲得湍流變量的壁面漸進(jìn)行為，如表1所列，其中δ為離開壁面的距離。

表1 近壁級數(shù)展開Table 1 Series expansion near wall

將表1中的脈動(dòng)量近壁關(guān)系式代入k、ε、kθ和εθ的定義式，可得：

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

本研究使用相應(yīng)的阻尼函數(shù)如式(13)～(16)所示，模型系數(shù)為Cε1=1.5、Cε2=1.9、σk=1.4、σε=1.4、Cp1=0.925、Cp2=0.9、Cd1=1.1、Cd2=0.9、σkθ=1.4、σεθ=1.4。蘇興康等[29]基于該模型計(jì)算了Pr=0.01～0.05的液態(tài)金屬在平板與圓管幾何中的湍流換熱，數(shù)值結(jié)果與DNS和實(shí)驗(yàn)關(guān)系式符合良好。

fε={1-0.3exp[-(Rt/6.5)2]}

(13)

fw=(1-exp(-Rε/19))2

(14)

fd2=1/Cd2(Cε2fε-1)[1-exp(-Rε/5.7)]2

(15)

(16)

最后，為封閉計(jì)算湍流運(yùn)動(dòng)黏度與湍流熱擴(kuò)散計(jì)算，采用引入Kolmogorov速度尺度的Abe湍流模型[20]計(jì)算νt，并采用適用于液態(tài)金屬近壁熱擴(kuò)散行為的Manservisi傳熱模型[15-17]計(jì)算αt。

2 數(shù)值求解

3 數(shù)值結(jié)果

3.1 數(shù)值驗(yàn)證

在快堆動(dòng)力系統(tǒng)中，液態(tài)鈉的常用Pr范圍為0.006～0.01，而液態(tài)鉛鉍為0.01～0.03[4]。為進(jìn)一步分析各向同性四方程模型在快堆常見燃料組件設(shè)計(jì)——三角形排列棒束內(nèi)的適用性，取Pr=0.01的液態(tài)金屬在其中的充分發(fā)展流動(dòng)傳熱進(jìn)行數(shù)值計(jì)算。圖1為三角形排列棒束的示意圖，線段AB、BC、CD和弧線DA是本研究分析的局部位置。Pr=0.01的典型物性取值[4]及流動(dòng)參數(shù)列于表2。與流體接觸的棒表面由均勻壁通量qw加熱。在其他面上應(yīng)用對稱條件來模擬具有三角形排列的無限棒束區(qū)域。進(jìn)出口采用周期性邊界條件cyclic，熱充分發(fā)展的流動(dòng)長度Lz設(shè)置為10Dh[27]。離開壁面的第1個(gè)網(wǎng)格距離設(shè)置為10-3mm，以滿足各向同性四方程模型的近壁湍流行為要求，即使在雷諾數(shù)達(dá)40萬的情況下，無量綱距離y+也滿足小于1的要求。經(jīng)網(wǎng)格無關(guān)性分析，最終P/D分別為1.25、1.34和1.46的計(jì)算域的六面體網(wǎng)格參數(shù)(徑向×周向×軸向)設(shè)置為60×80×180、70×80×240和80×80×300。

表2 物性及流動(dòng)參數(shù)Table 2 Physical property and flow parameter

圖1 計(jì)算域示意圖Fig.1 Schematic diagram of computational domain

取工程上重要的傳熱評估量——無量綱努塞爾數(shù)Nu進(jìn)行分析，Nu按以下定義進(jìn)行計(jì)算：

(17)

表3 三角形排列棒束努塞爾數(shù)經(jīng)驗(yàn)關(guān)系式Table 3 Correlation of Nusselt number for triangular bundle

3.2 數(shù)值結(jié)果

1) 冷卻劑與壁面溫度分布

P/D=1.25時(shí)三角形棒束的流向切面無量綱溫度分布云圖示于圖3，其中無量綱溫度θ+定義為λ(T-Tb)/qwDh。由圖3可見，間隙中心處的流體溫度大于通道中心處，且通道中心的流體溫度低于平均流體溫度Tb。壁面溫度沿周向不均勻分布，且靠近間隙的壁面溫度高于靠近通道中心的壁面溫度。如圖2a所示，在相同加熱功率和流動(dòng)面積的情況下，努塞爾數(shù)隨Pe的增加而增加，即通道內(nèi)的流體隨流速的增加得到充分加熱而壁面得到充分冷卻，使得流體平均溫度Tb增加而壁面溫度減小，因此在圖3中可觀察到Pe由1 000增加至4 000時(shí)，壁面無量綱溫度減小。

圖2 三角形棒束努塞爾數(shù)Fig.2 Nusselt number for triangular bundle

圖3 P/D=1.25的三角形棒束切面無量綱溫度Fig.3 Dimensionless temperature on slice for triangular bundle with P/D=1.25

線段AB、BC與CD和弧線DA的無量綱溫度隨P/D的變化示于圖4。隨著P/D的增加，間隙處的無量綱溫度減小。沿著間隙中心點(diǎn)B至通道中心點(diǎn)C的線上，無量綱溫度隨P/D的增加先減小后增加。從通道中心點(diǎn)C至壁面點(diǎn)D的線上，P/D越大，無量綱溫度越大?？拷ǖ乐行牡谋诿纥c(diǎn)D附近(0°<φ<10°)的壁面無量綱溫度變化不明顯，而在靠近間隙處的壁面點(diǎn)A附近，無量綱溫度區(qū)分明顯，且P/D越小，該處無量綱溫度越大。

圖4 不同P/D、Pe=2 000的三角形棒束線上溫度Fig.4 Dimensionless temperature over lines for triangular bundles at Pe=2 000 with different P/D

2) 溫度波動(dòng)及其耗散率分布

圖5 P/D=1.25的三角形棒束切面平均溫度波動(dòng)與其耗散率Fig.5 Average temperature fluctuaion and its isotropic dissipation on slice for triangular bundle with P/D=1.25

圖6 不同P/D、Pe=2 000的三角形棒束線上無量綱溫度波動(dòng)Fig.6 Dimensionless temperature fluctuation over lines for triangular bundles at Pe=2 000 with different P/D

3) 無量綱熱擴(kuò)散系數(shù)分布

線段AB、BC與CD的無量綱熱擴(kuò)散系數(shù)α+=αt/α隨P/D的變化示于圖7。其中，α+為由湍流熱擴(kuò)散引起的湍流熱通量與由分子熱傳導(dǎo)引起的分子熱通量之比。α+值越大，表示該處流場的湍流熱擴(kuò)散作用越強(qiáng)而分子熱傳導(dǎo)作用影響越小。在靠近壁面點(diǎn)A與點(diǎn)D附近，由于液態(tài)金屬的分子熱邊界較厚，α+?1，即該處的熱通量主要由分子熱傳導(dǎo)產(chǎn)生。在間隙AB處，α+隨P/D的增加變化不明顯。而在線段BC與CD上，P/D越大，α+整體越小。P/D為1.25時(shí)的棒束切面無量綱熱擴(kuò)散系數(shù)α+隨Pe的變化示于圖8。從圖8可看出，在低Pe下，如Pe=250和1 000，全流場α+均小于1；隨著Pe的增加，湍流熱擴(kuò)散作用增加，α+的峰值逐漸增加，開始出現(xiàn)湍流熱擴(kuò)散作用大于分子熱傳導(dǎo)作用的區(qū)域，即α+≥1的區(qū)域。

圖7 不同P/D、Pe=2 000的三角形棒束線上無量綱熱擴(kuò)散系數(shù)Fig.7 Dimensionless thermal diffusivity over lines for triangular bundles at Pe=2 000 with different P/D

圖8 P/D=1.25的三角形棒束切面無量綱熱擴(kuò)散系數(shù)Fig.8 Dimensionless thermal diffusivity on slice for triangular bundle with P/D=1.25

4) 湍流普朗特?cái)?shù)分布

線段AB、BC與CD的湍流熱擴(kuò)散系數(shù)Prt=vt/αt隨P/D的變化示于圖9?？拷诿纥c(diǎn)A與點(diǎn)D附近形成湍流普朗特?cái)?shù)局部較大值。經(jīng)過壁面附近峰值區(qū)域后，Prt向間隙中心或通道中心逐漸衰減。在間隙AB處，Prt隨P/D的增加變化不明顯。而在遠(yuǎn)離壁面的線段BC與CD區(qū)域，P/D越大，Prt整體分布越大。P/D為1.25時(shí)的棒束切面湍流普朗特?cái)?shù)Prt隨Pe的變化示于圖10。從圖10可看出，隨著Pe的增加，湍流熱擴(kuò)散作用增加，壁面附近的Prt較大值區(qū)域厚度向壁面衰減。

圖9 不同P/D、Pe=2 000的三角形棒束線上湍流普朗特?cái)?shù)Fig.9 Turbulent Prandtl number over lines for triangular bundles at Pe=2 000 with different P/D

圖10 P/D=1.25的三角形棒束切面湍流普朗特?cái)?shù)Fig.10 Turbulent Prandtl number on slice for triangular bundle with P/D=1.25

表4 不同P/D、Pe下的平均湍流普朗特?cái)?shù)Table 4 Average turbulent Prandtl numbers under different P/D and Pe

3.3 模型適應(yīng)性評價(jià)和輸入?yún)?shù)分析

4 結(jié)論