包蕓 高振源 葉孟翔

(中山大學(xué)力學(xué)系,廣州 510275)

1 引 言

熱對流現(xiàn)象廣泛存在于自然界和工業(yè)設(shè)計中,熱對流特性的研究有著重要的意義.Rayleigh-Bénard(RB)熱對流是從眾多熱對流過程中抽象出來的典型物理模型之一,是在一個封閉的空間內(nèi)下底板加熱上底板冷卻產(chǎn)生熱對流運動和熱輸運的系統(tǒng)[1].RB熱對流系統(tǒng)存在豐富而復(fù)雜的流動和熱輸運現(xiàn)象,一直受到國內(nèi)外學(xué)者的關(guān)注和研究.

國內(nèi)外關(guān)于RB熱對流的研究有一百多年的歷史,成果非常豐富.在理論研究方面,Grossman和Lohse[2,3]提出的GL理論最為成功,可以準確描述大范圍內(nèi)傳熱Nu數(shù)與Rayleigh(Ra)數(shù)和Pranntl(Pr)數(shù)的關(guān)系.Shishkina等[4]考慮了湍流脈動的影響,得到了湍流RB熱對流溫度邊界層方程.Wang等[5]研究了溫度脈動邊界層方程,并分析了溫度脈動溫度剖面的分布特性.RB對流的實驗[6?8]和數(shù)值模擬[9?11]研究也有許多進展.最新的實驗研究中,He等[12]進行了Pr=0.8極高Ra數(shù)的系列實驗,1012≤Ra≤1015,研究結(jié)果發(fā)現(xiàn),當Ra≤1013時,Nu～Ra0.31,Re～Ra0.43,當Ra數(shù)接近1015時,系統(tǒng)參數(shù)的標度率發(fā)生變化,Nu～Ra0.38,Re～Ra0.50,與GL理論預(yù)測的極高Ra數(shù)的結(jié)果一致.Zhou和Xia[13]對RB系統(tǒng)的黏性邊界層進行了研究,發(fā)現(xiàn)溫度邊界層分布與Prandtl-Blasius理論一致.Stevens等[14]通過直接數(shù)值模擬(DNS)方法對三維圓柱RB熱對流系統(tǒng)進行了模擬,Ra數(shù)達到2×1012,Pr數(shù)的范圍是0.5＜Pr＜10.黃茂靜和包蕓[15]計算了二維和三維的情況,發(fā)現(xiàn)二維熱對流和三維展向平均熱對流的時間平均場中都存在大尺度環(huán)流和角渦,但角渦的個數(shù)不同,二維和三維溫度邊界層厚度關(guān)于Ra數(shù)的變化都存在基本一致的標度率關(guān)系.鄒鴻岳等[16]討論了熱對流的軟硬湍流特性.包蕓等[17]討論了湍流熱對流大尺度環(huán)流反轉(zhuǎn)過程中角渦的變化特性.van der Poel等[18]通過DNS模擬計算,對比了Ra=106,108時的系列Pr數(shù)對傳熱的影響.Bao等[19]和Chen等[20]研究了加入豎直隔板并在頂端留有很小縫隙后整體傳熱效率的倍增現(xiàn)象.

本文對高Ra=1010的系列Pr數(shù)的二維湍流熱對流進行DNS計算研究.根據(jù)數(shù)值結(jié)果,探討不同Pr數(shù)情況下瞬時溫度場的羽流特性,以及Pr數(shù)對溫度邊界層和速度脈動的影響,并討論二維湍流熱對流在高Ra數(shù)下傳熱Nu數(shù)和湍流特性Re數(shù)隨Pr數(shù)的變化關(guān)系.

2 湍流熱對流DNS的并行直接求解

在Oberbeck-Boussinesq近似下,無量綱化的熱對流方程為

其中,V為無量綱速度矢量,θ為無量綱溫度,k為單位垂向矢量,p為壓力.無量綱參數(shù)Ra=(βg?θH3)/(κν)為瑞利數(shù),Pr=ν/κ為普朗特數(shù),Γ=W/H反映了對流系統(tǒng)的幾何尺寸.β為熱膨脹系數(shù),g為重力加速度,?θ為上下壁面溫差,H為系統(tǒng)裝置的高度,W為寬度,ν為運動黏性系數(shù),κ為熱擴散率.

由于湍流流動和熱羽流運動的復(fù)雜性,高Ra數(shù)湍流熱對流的DNS模擬計算工作量規(guī)模巨大,一直是DNS數(shù)值模擬研究湍流熱對流的困難所在.利用超級計算機計算資源,建立高效并行計算方法是實現(xiàn)高Ra數(shù)湍流熱對流DNS模擬研究的有效途徑.

二維熱對流DNS的求解過程采用投影法.動量方程和溫度對流擴散方程的計算,時間采用一階格式,空間采用二階格式,容易實現(xiàn)并行計算.由連續(xù)方程推導(dǎo)出的壓力泊松方程需要全流場聯(lián)立求解.在以往小規(guī)模的二維熱對流DNS模擬計算中,利用快速傅里葉變換解耦泊松方程,用追趕法求解三對角方程組的泊松方程直接解法[21],在單線程的計算上比用迭代求解方法有效很多.但三對角方程組的追趕法不易實現(xiàn)高效規(guī)模并行.Sun和Zhang[22]針對強對角占優(yōu)的三對角方程組提出了可并行計算的PDD算法.結(jié)合PDD算法建立實現(xiàn)規(guī)模并行的泊松方程直接求解方法.利用以上高效的壓力泊松方程并行直接求解方法,聯(lián)合其他易并行的動量方程等計算,創(chuàng)建了高效二維熱對流的并行直接求解方法(parallel direct method of DNS,PDM-DNS),具有很好的并行效率[23],并將此方法推廣到三維熱對流DNS模擬中.

本文采用PDM-DNS方法數(shù)值模擬二維高Ra數(shù)湍流熱對流,計算Ra=1010,0.05≤Pr≤20的十個算例.計算采用交錯網(wǎng)格,在上下邊界為非等距加密網(wǎng)格,當Pr≥0.4時,網(wǎng)格數(shù)為1024×1152(nx×ny),當Pr＜0.4時,網(wǎng)格數(shù)為1024×1408.速度邊界設(shè)為無滑移條件,左右兩側(cè)絕熱,上邊界低溫,給定為?0.5,下邊界高溫,給定為0.5.計算時間迭代步數(shù)均大于1000萬步.全部算例的計算時長都超過400個無量綱時間單位.采用MPI+OpenMP并行模式,在廣州超算中心“天河二號”超級計算機上完成所有的計算工作.同時還進行了三維情況下兩個Pr數(shù)的DNS計算,網(wǎng)格數(shù)為1024×128×1152.

3 Pr數(shù)對流動結(jié)構(gòu)的影響

Prandtl數(shù)是流體本身的物理性質(zhì),表征流體黏性耗散和熱耗散的關(guān)系,在RB熱對流中,也被用于衡量黏性邊界層和溫度邊界層之比.換言之,Pr數(shù)越高,黏性耗散則越高,熱耗散越低,因此系統(tǒng)的流動結(jié)構(gòu)是受Pr數(shù)控制的.

為了探討Pr數(shù)對高Ra數(shù)RB熱對流系統(tǒng)的影響,首先研究不同Pr數(shù)的流動結(jié)構(gòu)之間的差異.通過可視化技術(shù),給出不同Pr數(shù)情況下瞬時溫度場的分布,能夠反映Pr數(shù)對流動結(jié)構(gòu)的形態(tài)的影響.文中討論的溫度均為無量綱溫度.

圖1給出了八個典型Pr數(shù)的二維熱對流瞬時溫度場.為了更清晰地反映羽流運動情況,圖中溫度色標選用?0.1—0.1.從圖1中可以看出,不同Pr數(shù)的流動結(jié)構(gòu)存在明顯的差異.從上排Pr數(shù)較小的瞬時溫度場中看到,雜亂無序的冷羽流沿著左側(cè)邊壁下降,熱羽流則沿著右側(cè)邊壁上升,形成大尺度環(huán)流,在上下底板附近羽流的分布無明顯規(guī)律.此時羽流形狀多數(shù)呈不規(guī)則,分布較為混亂,而且隨著Pr數(shù)的增大,羽流的混亂情況在減弱.下排Pr數(shù)較大的瞬時溫度場顯示出隨Pr數(shù)的增大流動狀態(tài)的明顯變化.Pr=2.0時,方腔中出現(xiàn)細長的條狀羽流.當Pr=4.3和Pr=10.0時,方腔的兩側(cè)邊壁羽流都呈細長狀,具有較好的規(guī)律性.在大Pr數(shù)的瞬時溫度場中,冷羽流沿著左側(cè)邊壁下降,遇到上升的熱羽流后,流向下底板中間區(qū)域,熱羽流運動則相反,形成較穩(wěn)定的大尺度環(huán)流和角渦流動結(jié)構(gòu).

圖1 Ra=1010不同Pr數(shù)的瞬時溫度場 (a)Pr=0.05;(b)Pr=0.1;(c)Pr=0.4;(d)Pr=0.7;(e)Pr=1.2;(f)Pr=2.0;(g)Pr=4.3;(h)Pr=10.0Fig.1.Snapshots of the instantaneous temperature fi elds for Ra=1010with different Pr:(a)Pr=0.05;(b)Pr=0.1;(c)Pr=0.4;(d)Pr=0.7;(e)Pr=1.2;(f)Pr=2.0;(g)Pr=4.3;(h)Pr=10.0.

Breuer等[24]及Verizicco和Camussi[25]的 研究結(jié)果指出,Pr數(shù)較小時,系統(tǒng)的熱輸運主要由大尺度環(huán)流驅(qū)動,而Pr數(shù)較大時,羽流則起主導(dǎo)作用.本文的結(jié)果顯示,Pr＜2時,熱量主要通過不規(guī)則混亂的羽流形成的大尺度環(huán)流進行傳遞,而Pr≥2時,熱輸運的方式主要通過較穩(wěn)定的、繞大尺度環(huán)流和角渦運動的細長狀羽流進行.

4 不同Pr數(shù)溫度邊界層特性

本文引用了Zhou等[5]關(guān)于溫度邊界層厚度的定義,由此計算得到RB系統(tǒng)的溫度邊界層厚度,討論Pr數(shù)對溫度特性的影響.圖2給出了Pr數(shù)為2時,溫度的縱向分布以及靠近底板附近的局部放大圖.

圖2 Pr=2平均場溫度縱向分布,橫坐標為無量綱高度Fig.2.Averaged vertical pro fi le of temperature fi elds for Pr=2.The abscissa stands for non-dimensional height.

圖3給出了0.05≤Pr≤20,Ra=1010的溫度縱向分布,綠線為Pr=20,位于最上方,黑線為Pr=0.05,位于最下方.從圖中可以看出,不同Pr數(shù)溫度分布在近底板的線性部分基本重合,在稍遠離后溫度分布出現(xiàn)差異,Pr數(shù)越大,Θ(z)越高.當z/H≥0.04時,不同Pr數(shù)的溫度分布趨于一致.

圖3 近底板區(qū)溫度分布,不同顏色表示不同的Pr數(shù)Fig.3.Temperature distribution close to the bottom plate.The different colors of curves represent different Pr.

圖4給出了溫度邊界層厚度δ與Pr數(shù)的關(guān)系,可以看出,溫度邊界層厚度δ隨Pr數(shù)的變化緩慢,隨著Pr數(shù)的增加,溫度邊界層厚度逐漸減小.溫度邊界層厚度隨Pr數(shù)的變化滿足標度律關(guān)系,

可見溫度邊界層厚度隨Pr數(shù)的變化不大,而且存在標度律關(guān)系.

圖4 不同Pr數(shù)的溫度邊界層厚度Fig.4.Thickness of thermal boundary layer for different Pr.

5 Pr數(shù)對傳熱特性Nu數(shù)和湍流特性Re數(shù)的影響

5.1 不同Pr數(shù)的傳熱Nu數(shù)特性

傳熱效率Nu數(shù)和Re數(shù)是湍流熱對流系統(tǒng)的兩個重要的響應(yīng)參數(shù),GL理論給出了當前最有效的理論預(yù)測.Stevens等[26]結(jié)合大量的實驗數(shù)據(jù)重新修正了GL理論給出的Nu(Ra,Pr)關(guān)系中的系數(shù).GL理論關(guān)于Nu數(shù)的表述為:

根據(jù)Stevens給出的參數(shù),數(shù)值求解上述方程組,得到Ra=1010時的Nu(Pr)關(guān)系.

圖5 Nu數(shù)與Pr數(shù)的關(guān)系,空心菱形為二維結(jié)果,實心菱形為三維結(jié)果Fig.5.Nu as a function of Pr for 2D(hollow diamond)and 3D(solid diamond).

圖5給出了Ra=1010時的傳熱Nu數(shù)隨Pr數(shù)的變化情況,其中黑線為GL理論的預(yù)測結(jié)果,三維結(jié)果為Pr=0.7,4.3的兩個計算結(jié)果,二維為Pr=0.05,0.1,0.2,0.5,0.7,1.2,2.0,4.3,10.0,20.0的十個計算結(jié)果.從圖中可以看出,Ra=1010時,二維與三維的傳熱Nu數(shù)有較明顯的差異.三維結(jié)果與GL理論預(yù)測結(jié)果一致,而二維結(jié)果普遍低于GL理論的預(yù)測.二維結(jié)果的Nu數(shù)隨Pr數(shù)的變化過程,當Pr＜0.7時,Nu數(shù)隨Pr增加而增加,并在低Pr數(shù)時逐漸趨近GL理論的預(yù)測曲線;當0.7≤Pr≤4.3時,Nu數(shù)與Pr數(shù)的變化很小;當Pr＞10時,Nu數(shù)則有略微增加,并逐步接近GL理論預(yù)測值.在Stevens等[26]得到GL理論預(yù)測Nu(Ra,Pr)關(guān)系表明,當Pr數(shù)較大時,Nu數(shù)對Pr數(shù)不敏感,Pr數(shù)較小時,Nu數(shù)隨Pr數(shù)增加而增加.本文計算的Ra=1010二維結(jié)果與van der Poel等[18]DNS模擬計算較低Ra=106,108的系列Pr數(shù)對傳熱的影響,二維和三維出現(xiàn)的Nu數(shù)隨Pr數(shù)不同變化特性基本一致.

5.2 湍流特性和Re數(shù)

Rayleigh-Bénard熱對流的湍流流動很復(fù)雜,速度的脈動特性反映了湍流強弱.本節(jié)討論Pr數(shù)對速度脈動的影響,圖6給出了Pr=0.05,0.7,10的水平速度脈動信號,信號源的位置為靠近下底板的中心處,時長為100個無量綱單位時間.

圖6 x=0.5,z=0.02處水平速度脈動,橫坐標為無量綱時間,縱坐標為瞬時速度 (a)Pr=0.05;(b)Pr=0.7;(c)Pr=10Fig.6.Fluctuation of horizontal velocity at x=0.5 and z=0.02:(a)Pr=0.05;(b)Pr=0.7;(c)Pr=10.

圖6為不同Pr數(shù)在方腔x=0.5L,z=0.02H處的水平速度脈動情況,此處為大尺度環(huán)流的水平剪切區(qū),流動主要為水平方向.從圖6中可以看出,對不同Pr數(shù)的速度脈動信號差異顯著.圖6(a)中,Pr=0.05時,速度的變化范圍為?0.96≤u≤1.78,速度平均值約為0.56.當Pr=0.7時,圖6(b)中,速度的變化范圍?0.47≤u≤0.52,平均值約為0.20.當Pr數(shù)較大,圖6(c)中Pr=10時,速度脈動的振幅變得很小,速度平均值只有0.09,無論速度的平均值還是速度振幅都遠小于Pr=0.05時的情況.由此可見,Pr數(shù)越小,流場速度和速度脈動振幅越大,速度隨時間的變化程度更為劇烈,流場中的湍流度也越高.

流場中的Re數(shù)是反映速度場湍流特性的重要參數(shù).選取不同的特征速度,可計算出不同的Re數(shù).本文采用兩種計算Re數(shù)的方法,一種是通過給定空間點處速度脈動的時均值〈u〉計算得到的Re〈u〉,另一種是通過平均場的速度Umax計算得到的ReUmax.

圖7 Re數(shù)與Pr數(shù)的關(guān)系,藍色三角為Re〈u〉,紅色圓為ReUmaxFig.7.Re〈u〉 (blue triangle)and ReUmax(red circle)as a function of Pr.

圖7給出了高Ra=1010時兩種Re數(shù)隨Pr數(shù)的變化,藍色表示由脈動速度〈u〉計算得到的Re〈u〉,紅色表示由平均場最大速度Umax計算得到的ReUmax.圖中可見,ReUmax的值大于Re〈u〉的值.隨著Pr數(shù)的增加,兩種Re數(shù)均減少:當Pr=0.05時,Re〈u〉=2.37×105,ReUmax=4.04×105;當Pr=20時,Re〈u〉=1.84×103,ReUmax=3.53×103.并且Re數(shù)隨Pr數(shù)變化存在標度律關(guān)系,計算得到

二者的標度率指數(shù)基本相同.由此可見,在高Ra系列Pr數(shù)的計算中,兩種Re數(shù)的大小值存在明顯的差異,但隨Pr數(shù)的變化滿足一致的標度律關(guān)系.

6 結(jié) 論

Rayleigh-Bénard熱對流是一個經(jīng)典的湍流熱對流問題,直接數(shù)值模擬是重要的研究手段.通過數(shù)值計算可以模擬任意Pr數(shù)和Ra數(shù)的情況,分析實驗中不易測量的物理量,對熱對流領(lǐng)域的發(fā)展具有重要的研究意義.本文使用高效并行直接數(shù)值求解方法模擬了高Ra=1010時系列Pr數(shù)(0.05≤Pr≤20)湍流熱對流,對比分析了不同Pr數(shù)下流動結(jié)構(gòu)、溫度邊界層特性和速度脈動特性,討論了傳熱效率Nu數(shù)和Re數(shù)與Pr數(shù)的關(guān)系.得到以下結(jié)論:

1)隨著Pr數(shù)的變化,RB熱對流系統(tǒng)的流動結(jié)構(gòu)發(fā)生改變;當Pr數(shù)較低時,羽流的形狀不規(guī)則,流場中存在較多的小尺寸羽流團,運動方式較為復(fù)雜;當Pr數(shù)較高時,羽流主要表現(xiàn)為規(guī)則細長狀,羽流運動使流場出現(xiàn)典型的大尺度環(huán)流和角渦;

2)在高Ra數(shù)時,不同Pr數(shù)的溫度Θ(z)縱向分布差異較小;對0.05≤Pr≤20,線性區(qū)溫度曲線基本重合,離開線性區(qū)后,Pr數(shù)越高,溫度曲線分布越高;溫度邊界層厚度隨Pr數(shù)的變化緩慢減小,遵循標度率δ～Pr?0.035;速度脈動信號在Pr數(shù)較小時振幅較大,波動更為強烈,湍流強度更高;

3)在高Ra數(shù)時,當Pr＜1時傳熱Nu數(shù)隨著Pr數(shù)增加而增加,當Pr＞1時Nu數(shù)隨著Pr數(shù)變化不明顯;通過速度脈動信號和平均速度縱向分布最大值計算得到的Re〈u〉和ReUmax隨Pr數(shù)增加而減小,二者基本遵循同一標度律:Re～Pr?0.81.

[1]Ahlers G,Grossmann S,Lohse D 2009Rev.Mod.Phys.81 503

[2]Grossmann S,Lohse D 2000J.Fluid Mech.407 27

[3]Grossmann S,Lohse D 2001Phys.Rev.Lett.86 3316

[4]Shishkina O,Horn S,Wagner S,Ching E S C 2015Phys.Rev.Lett.114 114302

[5]Wang Y,He X Z,Tong P 2016Phys.Rev.Fluids1 08230

[6]Salort J,Liot O,Rusaouen E 2014Phys.Fluids26 015112

[7]Wei P,Chan T S,Ni R,Zhao X Z 2014J.Fluid Mech.740 28

[8]Xie Y C,Xia K Q 2017J.Fluid Mech.825 573

[9]Wagner S,Shishkina O 2015J.Fluid Mech.763 109

[10]Toppaladoddi S,Succi S,Wettlaufer J S 2017Phys.Rev.Lett.118 074503

[11]Zhu X J,Stevens R J A M,Verzicco R,Lohse D 2017Phys.Rev.Lett.119 154501

[12]He X Z,Funfschilling D,Nobach H,Bodenschatz E,Ahlers G 2012Phys.Rev.Lett.108 024502

[13]Zhou Q,Xia K Q 2010Phys.Rev.Lett.104 104301

[14]Stevens R J A M,Lohse D,Verzicoo R 2011J.Fluid Mech.688 31

[15]Huang M J,Bao Y 2016Acta Phys.Sin.20 204702(in Chinese)[黃茂靜,包蕓 2016物理學(xué)報 20 204702]

[16]Zhou H Y,Xu W,Bao Y 2014J.Hydrodyn.29 34(in Chinese)[鄒鴻岳,徐煒,包蕓2014水動力學(xué)研究與進展A輯29 34]

[17]Bao Y,Ning H,Xu W 2014Acta Phys.Sin.63 154703(in Chinese)[包蕓,寧浩,徐煒2014物理學(xué)報 63 154703]

[18]van der Poel E P,Stevens R J A M,Lohse D 2013J.Fluid Mech.736 177

[19]Bao Y,Chen J,Liu B F 2015J.Fluid Mech.784 R5

[20]Chen J,Bao Y,Yin Z X 2017Int.J.Heat Mass Tran.115 556

[21]Xu W,Bao Y 2013Acta Mech.Sin.45 1(in Chinese)[徐煒,包蕓2013力學(xué)學(xué)報 45 1]

[22]Sun X H,Zhang W 2004IEEE Trans.Parall.Distr.15 97

)

[24]Breuer M,Wessling S,Schmalzl J,Hansen U 2004Phys.Rev.E69 026302

[25]Verzicco R,Camussi R 1999J.Fluid Mech.383 55

[26]Stevens R J A M,van der Poel E P,Grossmann S,Lohse D 2013J.Fluid Mech.730 295