基于模型的液態金屬三角形棒束數值傳熱研究

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

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

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

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

為改善SGDH方法的適用性，部分學者基于全局流動參數或局部湍流參數建立了用于計算液態金屬的非線性Prt關系式，如Cheng等[11]引入全局雷諾數Re和普朗特數Pr效應建立了全局Prt函數關系，Kays等[12]引入湍流運動黏度和普朗特數Pr建立了局部Prt函數關系。Duponcheel等[13]分析了現有的多個Prt關系式在簡單幾何中的計算精度，認為基于局部湍流參數建立的Prt關系式能較好地提高液體金屬的計算精度，但在復雜幾何中的適用性還需進一步分析。

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

1 湍流換熱模型

基于Boussinesq假設和SGDH的周期性充分發展RANS方程組如下：

(1)

(2)

uz4qw/ρcpubDh

(3)

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

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

表1 近壁級數展開Table 1 Series expansion near wall

將表1中的脈動量近壁關系式代入k、ε、kθ和εθ的定義式，可得：

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

本研究使用相應的阻尼函數如式(13)～(16)所示，模型系數為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]基于該模型計算了Pr=0.01～0.05的液態金屬在平板與圓管幾何中的湍流換熱，數值結果與DNS和實驗關系式符合良好。

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)

最后，為封閉計算湍流運動黏度與湍流熱擴散計算，采用引入Kolmogorov速度尺度的Abe湍流模型[20]計算νt，并采用適用于液態金屬近壁熱擴散行為的Manservisi傳熱模型[15-17]計算αt。

2 數值求解

3 數值結果

3.1 數值驗證

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

表2 物性及流動參數Table 2 Physical property and flow parameter

圖1 計算域示意圖Fig.1 Schematic diagram of computational domain

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

(17)

表3 三角形排列棒束努塞爾數經驗關系式Table 3 Correlation of Nusselt number for triangular bundle

3.2 數值結果

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

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

圖2 三角形棒束努塞爾數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的增加，間隙處的無量綱溫度減小。沿著間隙中心點B至通道中心點C的線上，無量綱溫度隨P/D的增加先減小后增加。從通道中心點C至壁面點D的線上，P/D越大，無量綱溫度越大。靠近通道中心的壁面點D附近(0°<φ<10°)的壁面無量綱溫度變化不明顯，而在靠近間隙處的壁面點A附近，無量綱溫度區分明顯，且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) 溫度波動及其耗散率分布

圖5 P/D=1.25的三角形棒束切面平均溫度波動與其耗散率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的三角形棒束線上無量綱溫度波動Fig.6 Dimensionless temperature fluctuation over lines for triangular bundles at Pe=2 000 with different P/D

3) 無量綱熱擴散系數分布

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

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

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

4) 湍流普朗特數分布

線段AB、BC與CD的湍流熱擴散系數Prt=vt/αt隨P/D的變化示于圖9。靠近壁面點A與點D附近形成湍流普朗特數局部較大值。經過壁面附近峰值區域后，Prt向間隙中心或通道中心逐漸衰減。在間隙AB處，Prt隨P/D的增加變化不明顯。而在遠離壁面的線段BC與CD區域，P/D越大，Prt整體分布越大。P/D為1.25時的棒束切面湍流普朗特數Prt隨Pe的變化示于圖10。從圖10可看出，隨著Pe的增加，湍流熱擴散作用增加，壁面附近的Prt較大值區域厚度向壁面衰減。

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

圖10 P/D=1.25的三角形棒束切面湍流普朗特數Fig.10 Turbulent Prandtl number on slice for triangular bundle with P/D=1.25

表4 不同P/D、Pe下的平均湍流普朗特數Table 4 Average turbulent Prandtl numbers under different P/D and Pe

3.3 模型適應性評價和輸入參數分析

4 結論