陽祥,陶文銓

(1.國家核電技術(shù)有限公司北京研發(fā)中心, 100190, 北京; 2.西安交通大學能源與動力工程學院, 710049, 西安)

高瑞利數(shù)下封閉腔內(nèi)自然對流的數(shù)值模擬

陽祥1,陶文銓2

(1.國家核電技術(shù)有限公司北京研發(fā)中心, 100190, 北京; 2.西安交通大學能源與動力工程學院, 710049, 西安)

為了推廣應(yīng)用高瑞利數(shù)下的自然對流換熱技術(shù),有必要對自然對流流動與換熱特性進行深入研究。采用不引入人工擾動的直接數(shù)值模擬方法,對發(fā)生在高寬比為4的封閉腔內(nèi)的自然對流流動與換熱進行了研究,分析了平均溫度、平均主流速度、渦量和局部努塞爾數(shù)的分布特性。研究結(jié)果表明:從靜止等溫流體初始條件出發(fā),不引入任何人工擾動自然對流可以順利發(fā)展到湍流,節(jié)約了計算資源;即便瑞利數(shù)等于1010,自然對流的平均溫度、平均主流速度、渦量和局部努賽爾數(shù)分布都具有邊界層型流動和換熱的特征;在普朗特數(shù)為0.71～500的范圍,當封閉腔內(nèi)自然對流換熱出現(xiàn)湍流換熱特征時,局部瑞利數(shù)處于107～108量級。

封閉腔;湍流;自然對流;直接數(shù)值模擬

豎直封閉腔內(nèi)的自然對流是傳熱與流動學科中的一個經(jīng)典問題,有很多工程應(yīng)用背景,如核反應(yīng)堆中的換熱、放射性廢料的冷卻、房間通風、太陽能集熱器和電子器件冷卻等。豎直壁面存在溫差的封閉腔內(nèi)沿高溫、低溫豎直壁面會分別形成上升流動和下降流動,并且在腔內(nèi)形成回流,高瑞利數(shù)Ra流動會形成層流、過渡和湍流區(qū)域。Tian和Karayiannis對弱湍流自然對流進行了實驗研究,首次發(fā)表了加熱封閉腔內(nèi)溫度等值線分布和速度矢量分布的實驗結(jié)果[1-2]。Giel等對幾種高寬比的封閉腔內(nèi)自然對流傳熱(高Ra)進行了試驗研究,測量了速度和溫度的平均值和脈動值分布,并進一步獲得了脈動速度和溫度的頻譜圖;在分析試驗結(jié)果的基礎(chǔ)上他們認為高寬比大于等于4的封閉腔內(nèi)發(fā)生湍流自然對流時展向的尺寸大小對流動和換熱影響很小[3]。在數(shù)值模擬方面的研究要少得多,尤其是直接模擬和大渦模擬的研究更少,豎直壁面自然對流通常在高Ra下才會轉(zhuǎn)變?yōu)橥牧?并且從層流過渡到湍流很難處理,對計算技術(shù)和資源要求很高。Paolucci首次對豎直壁面存在溫差的方腔自然對流問題進行直接模擬,他采用有限差分法對二維封閉方腔內(nèi)湍流自然對流進行了模擬,證實了自然對流從層流到湍流直接模擬的可行性[4]。Trias等采用直接數(shù)值模擬方法研究了高Ra下二維封閉腔和三維封閉腔的自然對流,結(jié)果表明二維和三維計算在模擬封閉腔內(nèi)自然對流和換熱的一般特征時基本沒有差別,二維計算更為經(jīng)濟實用,三維計算具有更豐富的信息,特別是在湍流高階統(tǒng)計量方面[5]。

本文對幾種物質(zhì)(不同普朗特數(shù)Pr)在高寬比為4的封閉腔內(nèi)發(fā)生的自然對流采用經(jīng)濟實用的二維直接數(shù)值模擬,分析了包括速度場、溫度場、渦量場和局部努塞爾數(shù)Nuy分布的特點,揭示了自然對流從層流到湍流的流動與換熱特性。

1 物理模型和數(shù)學模型

1.1 物理模型

圖1 物理模型示意圖

本文采用如圖1所示的物理模型,x坐標軸沿水平方向,y坐標軸沿豎直方向,坐標原點在二維封閉腔的中心,方腔高寬比為4,水平壁面絕熱,兩豎直壁面溫度分別保持恒定,并且左壁面溫度高右壁面溫度低。

1.2 數(shù)學模型

采用Boussinesq假設(shè),歸一化的二維質(zhì)量守恒方程、動量守恒方程和能量守恒方程如下

(1)

(2)

RaPrΘ

(3)

(4)

1.3 控制方程的離散及數(shù)值方法

采用分步算法求解非穩(wěn)態(tài)N-S方程,采用有限差分法離散控制方程,對流項和擴散項都采用二階相容中心差分格式,在時間離散上采用二階半隱半顯的混合Adams-Bashforth/Crank-Nicholson格式。

湍流自然對流的直接數(shù)值模擬需要滿足空間和時間尺度要求。根據(jù)George等給出的經(jīng)驗公式[6]來估算自然對流中的導熱次層厚度,如下式所示

λ≈1.7(PrRam)-1/3

(5)

式中:Ram是基于ΔTm=Th-T0定義的,其中T0=0.5(Th+Tc)。時間尺度也應(yīng)滿足一定的要求,根據(jù)Paolucci的經(jīng)驗[3],時間尺度滿足下式要求

Δt≤8π(Ra3Pr)-1/4

(6)

按照上述要求,本文中選取的網(wǎng)格數(shù)沿x和y方向分別為160×320,x方向的最大網(wǎng)格距離為9.8×10-3,y方向最大網(wǎng)格距離為1.96×10-2,最小網(wǎng)格距離在貼近壁面處,沿豎直壁面和水平壁面的最小網(wǎng)格距離均為1×10-5,選取的時間步長為Δt=1×10-8。

速度采用無滑移邊界條件,高溫壁面的溫度邊界條件為Θ=0.5,低溫壁面為Θ=-0.5,水平壁面上絕熱。計算初場為0,為了得到統(tǒng)計解,在完全進入湍流狀態(tài)后,再進行至少200萬步左右的時間推進,每隔1萬步進行數(shù)據(jù)采集,對采集到的瞬時數(shù)據(jù)進行時間平均,得到統(tǒng)計解。在數(shù)據(jù)處理中首先完成了數(shù)值結(jié)果對網(wǎng)格數(shù)量的依賴性,文中展示的計算數(shù)據(jù)都是完成網(wǎng)格無關(guān)性考核后的結(jié)果。

1.4 程序驗證

圖2 邊界層外區(qū)域歸一化主流速度分布

2 計算結(jié)果及分析

2.1 流場和溫度場

圖3給出了不同Pr下的等溫線分布,由圖可知等溫線的拓撲結(jié)構(gòu)關(guān)于x=0反對稱。沿著豎直高溫壁面流體朝上呈邊界層型流動,在底部區(qū)域溫度在很薄的流體層內(nèi)變化,超過封閉腔一半高度后等溫線明顯朝上凸出,溫度邊界層明顯增厚,自然對流失穩(wěn)并進入湍流狀態(tài)。低溫壁面?zhèn)攘黧w的溫度變化與高溫側(cè)相似,但是方向發(fā)生了變化。在封閉腔的中心區(qū)域內(nèi)等溫線沿x方向的變化較小,隨著Pr增加沿水平方向溫度變化滿足?Θ/?x≈0,流體呈明顯的分層狀態(tài)。這是因為隨著Pr增加,流體黏滯作用相對增強,溫度邊界層減薄的緣故。

圖4給出了平均主流速度V和平均溫度Θ沿y=0的分布,V和Θ關(guān)于x=0反對稱,文中只給出高溫側(cè)的溫度和速度分布。在y<0的區(qū)域內(nèi)自然對流已經(jīng)完成從層流到湍流的過渡,因此圖中的速度和溫度分布都是湍流狀態(tài)下的。從圖4a可看出,速度分布定性上與豎直平板自然對流的速度分布相似:當-0.5-0.35時主流速度幾乎為0。圖4b中是平均溫度的分布,溫度在-0.5-0.46時溫度已降至0。比較圖4b與圖4a可知,溫度邊界層厚度明顯小于速度邊界層厚度。

圖5是某一時刻封閉腔內(nèi)自然對流的渦量分布,渦量的定義為ω=?v/?x-?u/?y,圖中深色表示渦量大,淺色表示渦量小。渦量的分布展現(xiàn)了在浮升力驅(qū)動下流體運動沿豎直壁面的發(fā)展軌跡,下面以高溫壁面(每幅圖的左側(cè)壁面)流動為例進行說明。在高溫表面下部,渦量分布很規(guī)律,流動處于層流狀態(tài);靠近中間高度處,渦量分布變化較劇烈,流體呈拋射狀向上運動,流動已經(jīng)失穩(wěn)伴有湍流的特征;超過中間高度后有明顯的間歇性渦結(jié)構(gòu)出現(xiàn),渦量變化劇烈,表現(xiàn)出明顯的湍流特征,流動已發(fā)展成旺盛的湍流;在靠近頂部的區(qū)域,由于上頂面的阻擋作用,渦量變化較劇烈,流體的脈動性減弱。在低溫壁面,邊界層型流動向下發(fā)展,壁面附近流體的渦量分布特性與高溫壁面?zhèn)认嗨?。隨著Pr的變化,渦尺度發(fā)生了較明顯的變化:Pr=0.71時湍流邊界層內(nèi)的渦尺度較大,可以清楚地觀察到流體邊界層在空間的演變,在左右兩個流動邊界層交匯的地方流體邊界層之間發(fā)生了輕微的干涉;當Pr≥10時,渦尺度明顯變小且呈長條狀,數(shù)量增多,但流動邊界層從層流到湍流的發(fā)展過程仍清晰可辨,渦量在中心區(qū)域內(nèi)較弱,流體處于分層狀態(tài)。

(a)Pr=0.71(b)Pr=10(c)Pr=50(d)Pr=100 (e)Pr=500

(a)V沿y=0水平線的分布

(b)Θ沿y=0水平線的分布

2.2 傳熱特性

圖6展示了沿高溫壁面局部努塞爾數(shù)Nuy隨局部瑞利數(shù)Ray的變化曲線,Nuy定義為Nuy=hy/k,Ray定義為Ray=gβΔTy3/αν。由圖可知,每條Nuy-Ray曲線都由層流區(qū)、過渡區(qū)、旺盛湍流區(qū)、頂部區(qū)4個區(qū)間組成。圖6a給出了Pr=0.71時Nuy隨Ray的變化情況,從高溫壁面的底部起的一段區(qū)域內(nèi),流動為層流狀態(tài),Nuy呈線性下降;流動發(fā)展到過渡狀態(tài)后,Nuy曲線發(fā)生振蕩但頻率較低、幅度較小;隨著流動發(fā)展到旺盛湍流,Nuy曲線出現(xiàn)明顯的高頻次振蕩;隨著高度的進一步增加,壁面阻擋作用開始顯現(xiàn),流體脈動減弱,Nuy曲線又重新回歸到接近直線分布,并且Nuy值急劇下降。隨著Pr的增加,Nuy曲線振蕩的幅值和頻率都有所下降,這是因為黏滯作用相對增加了;在接近上壁面區(qū)域內(nèi),Nuy值下降得更快,并有層流化趨勢。

根據(jù)Nuy-Ray曲線的分布特點,發(fā)現(xiàn)了一個非常有意思的現(xiàn)象:在Pr=0.71～500的范圍內(nèi),Nuy都在Ray=107～108這個范圍內(nèi)開始發(fā)生振蕩,此Ray范圍對應(yīng)的高度與溫度、速度和渦量分布圖中出現(xiàn)湍流特征的高度是一致的。因此,通過觀察Nuy-Ray曲線,可以粗略地判斷在Ray=107～108范圍內(nèi),自然對流換熱具有湍流換熱的特點。

(a)Pr=0.71 (b)Pr=10(c)Pr=50(d)Pr=100 (e)Pr=500

(a)Pr=0.71 (b)Pr=10 (c)Pr=50

(d)Pr=100 (e)Pr=500

3 結(jié) 論

本文對高寬比為4的封閉腔內(nèi)自然對流進行了直接數(shù)值模擬,研究了自然對流從層流到湍流的流動與傳熱特性,得出以下結(jié)論:

(1)實現(xiàn)了從靜止等溫流體初始條件出發(fā),在不引入任何人工擾動的情況下完成了封閉腔內(nèi)自然對流的直接數(shù)值模擬,大大縮短了計算時間,節(jié)約了計算資源;

(2)封閉腔內(nèi)自然對流的主要物理量沿等溫壁面的分布,包括平均溫度、平均主流速度、渦量和局部努賽爾數(shù),具有邊界層型流動和換熱的特點,這些物理量的分布都展示了封閉腔內(nèi)自然對流從層流到湍流的發(fā)展演變過程;

(3)不同Pr下封閉腔內(nèi)自然對流換熱有一定的共性,在Pr=0.71～500的范圍內(nèi),當Ray=107～108量級時,自然對流換熱具有明顯的湍流換熱特征。

[1] TIAN Y S, KARAYIANNIS T G. Low turbulence natural convection in an air filled square cavity: Part I The thermal and fluid flow fields [J]. International Journal of Heat and Mass Transfer, 2000, 43(6): 849-866.

[2] TIAN Y S, KARAYIANNIS T G. Low turbulence natural convection in an air filled square cavity: Part II The turbulence quantities [J]. International Journal of Heat and Mass Transfer, 2000, 43(6): 867-884.

[3] GIEL P, SCHMIDT F. An experimental study of high Rayleigh natural convection in an enclosure [C]∥Proceedings of the 8th International Heat Transfer Conference. San Francisco, USA: Hemisphere Publishing Corporation, 1986: 1459-1464.

[4] PAOLUCCI S. Direct numerical simulation of two-dimensional turbulent natural convection in an enclosed cavity [J]. Journal of Fluid Mechanics, 1990, 215: 229-262.

[5] TRIAS F X, SORIA M, OLIVA A. Direct numerical simulations of two-and three-dimensional turbulent natural convection flows in a differentially heated cavity of aspect ratio 4 [J]. J Fluid Mech, 2007, 586(3): 259-293.

[6] GEORGE W K. A theory for natural convection turbulent boundary layers next to heated vertical surfaces [J]. International Journal of Heat and Mass Transfer, 1979, 22(6): 813-826.

[7] FARHANGNIA M, BIRINGEN S, PELTIER L J. Numerical simulation of two-dimensional buoyancy-driven turbulence in a tall rectangular cavity [J]. International Journal for Numerical Methods in Fluids, 1998, 23(12): 1311-1326.

(編輯 荊樹蓉)

NumericalSimulationsforNaturalConvectionwithHighRayleighNumberinaTallRectangularCavity

YANG Xiang1,TAO Wenquan2

(1. Beijing R&D Center of State Nuclear Power Technology Corporation, Beijing 100190, China;2. School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China)

To generalize heat transfer technology of natural convection with high Rayleigh number, it is necessary to extensively investigate the characteristics of fluid flow and heat transfer by natural convection. For a tall rectangular cavity of an aspect ratio 4, direct numerical simulations are carried out without any artificial perturbations on the flow field, and the distributions of averaged physical quantities such as temperature, mainstream velocity, vorticity and local Nusselt number are analyzed. The results show that: 1) Starting under quiescent and isothermal flow conditions, the flow can be driven to turbulence without any artificial perturbations, which saves computational resources; 2) As Rayleigh numbers up to 1010, the distributions of temperature, mainstream velocity, vorticity and local Nusselt numbers for natural convection in the tall cavity are endowed with the similar features of those in boundary-layer flow on a heated vertical plate; and 3) With Prandtl number in range of 0.71-500 and local Rayleigh number of magnitudes 107-108, heat transfer on the isothermal surfaces of tall cavity already presents the characteristics of turbulent flow.

enclosure; turbulence; natural convection; direct numerical simulation

10.7652/xjtuxb201405005

2013-08-21。 作者簡介: 陽祥(1980—),男,工程師。 基金項目: 國家自然科學基金重點資助項目(51136004)。

時間: 2014-02-26 網(wǎng)絡(luò)出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20140226.1203.031.html

TK124

:A

:0253-987X(2014)05-0027-05