陽(yáng)祥,陶文銓

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

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

陽(yáng)祥1,陶文銓2

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

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

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

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

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

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

1.1 物理模型

圖1 物理模型示意圖

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

1.2 數(shù)學(xué)模型

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

(1)

(2)

RaPrΘ

(3)

(4)

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

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

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

λ≈1.7(PrRam)-1/3

(5)

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

Δ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,選取的時(shí)間步長(zhǎng)為Δt=1×10-8。

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

1.4 程序驗(yàn)證

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

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

2.1 流場(chǎng)和溫度場(chǎng)

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

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

圖5是某一時(shí)刻封閉腔內(nèi)自然對(duì)流的渦量分布,渦量的定義為ω=?v/?x-?u/?y,圖中深色表示渦量大,淺色表示渦量小。渦量的分布展現(xiàn)了在浮升力驅(qū)動(dòng)下流體運(yùn)動(dòng)沿豎直壁面的發(fā)展軌跡,下面以高溫壁面(每幅圖的左側(cè)壁面)流動(dòng)為例進(jìn)行說(shuō)明。在高溫表面下部,渦量分布很規(guī)律,流動(dòng)處于層流狀態(tài);靠近中間高度處,渦量分布變化較劇烈,流體呈拋射狀向上運(yùn)動(dòng),流動(dòng)已經(jīng)失穩(wěn)伴有湍流的特征;超過中間高度后有明顯的間歇性渦結(jié)構(gòu)出現(xiàn),渦量變化劇烈,表現(xiàn)出明顯的湍流特征,流動(dòng)已發(fā)展成旺盛的湍流;在靠近頂部的區(qū)域,由于上頂面的阻擋作用,渦量變化較劇烈,流體的脈動(dòng)性減弱。在低溫壁面,邊界層型流動(dòng)向下發(fā)展,壁面附近流體的渦量分布特性與高溫壁面?zhèn)认嗨啤ｋS著Pr的變化,渦尺度發(fā)生了較明顯的變化:Pr=0.71時(shí)湍流邊界層內(nèi)的渦尺度較大,可以清楚地觀察到流體邊界層在空間的演變,在左右兩個(gè)流動(dòng)邊界層交匯的地方流體邊界層之間發(fā)生了輕微的干涉;當(dāng)Pr≥10時(shí),渦尺度明顯變小且呈長(zhǎng)條狀,數(shù)量增多,但流動(dòng)邊界層從層流到湍流的發(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個(gè)區(qū)間組成。圖6a給出了Pr=0.71時(shí)Nuy隨Ray的變化情況,從高溫壁面的底部起的一段區(qū)域內(nèi),流動(dòng)為層流狀態(tài),Nuy呈線性下降;流動(dòng)發(fā)展到過渡狀態(tài)后,Nuy曲線發(fā)生振蕩但頻率較低、幅度較小;隨著流動(dòng)發(fā)展到旺盛湍流,Nuy曲線出現(xiàn)明顯的高頻次振蕩;隨著高度的進(jìn)一步增加,壁面阻擋作用開始顯現(xiàn),流體脈動(dòng)減弱,Nuy曲線又重新回歸到接近直線分布,并且Nuy值急劇下降。隨著Pr的增加,Nuy曲線振蕩的幅值和頻率都有所下降,這是因?yàn)轲饔孟鄬?duì)增加了;在接近上壁面區(qū)域內(nèi),Nuy值下降得更快,并有層流化趨勢(shì)。

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

(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é) 論

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

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

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

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

[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。 作者簡(jiǎn)介: 陽(yáng)祥(1980—),男,工程師。 基金項(xiàng)目: 國(guó)家自然科學(xué)基金重點(diǎn)資助項(xiàng)目(51136004)。

時(shí)間: 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