湍流熱對流Prandtl數效應的數值研究?

包蕓 高振源 葉孟翔

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

1 引 言

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

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

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

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

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

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

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

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

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

3 Pr數對流動結構的影響

Prandtl數是流體本身的物理性質,表征流體黏性耗散和熱耗散的關系,在RB熱對流中,也被用于衡量黏性邊界層和溫度邊界層之比.換言之,Pr數越高,黏性耗散則越高,熱耗散越低,因此系統的流動結構是受Pr數控制的.

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

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

圖1 Ra=1010不同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.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]的 研究結果指出,Pr數較小時,系統的熱輸運主要由大尺度環流驅動,而Pr數較大時,羽流則起主導作用.本文的結果顯示,Pr＜2時,熱量主要通過不規則混亂的羽流形成的大尺度環流進行傳遞,而Pr≥2時,熱輸運的方式主要通過較穩定的、繞大尺度環流和角渦運動的細長狀羽流進行.

4 不同Pr數溫度邊界層特性

本文引用了Zhou等[5]關于溫度邊界層厚度的定義,由此計算得到RB系統的溫度邊界層厚度,討論Pr數對溫度特性的影響.圖2給出了Pr數為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數溫度分布在近底板的線性部分基本重合,在稍遠離后溫度分布出現差異,Pr數越大,Θ(z)越高.當z/H≥0.04時,不同Pr數的溫度分布趨于一致.

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

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

可見溫度邊界層厚度隨Pr數的變化不大,而且存在標度律關系.

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

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

5.1 不同Pr數的傳熱Nu數特性

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

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

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

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

5.2 湍流特性和Re數

Rayleigh-Bénard熱對流的湍流流動很復雜,速度的脈動特性反映了湍流強弱.本節討論Pr數對速度脈動的影響,圖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數在方腔x=0.5L,z=0.02H處的水平速度脈動情況,此處為大尺度環流的水平剪切區,流動主要為水平方向.從圖6中可以看出,對不同Pr數的速度脈動信號差異顯著.圖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數較大,圖6(c)中Pr=10時,速度脈動的振幅變得很小,速度平均值只有0.09,無論速度的平均值還是速度振幅都遠小于Pr=0.05時的情況.由此可見,Pr數越小,流場速度和速度脈動振幅越大,速度隨時間的變化程度更為劇烈,流場中的湍流度也越高.

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

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

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

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

6 結 論

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

1)隨著Pr數的變化,RB熱對流系統的流動結構發生改變;當Pr數較低時,羽流的形狀不規則,流場中存在較多的小尺寸羽流團,運動方式較為復雜;當Pr數較高時,羽流主要表現為規則細長狀,羽流運動使流場出現典型的大尺度環流和角渦;

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

3)在高Ra數時,當Pr＜1時傳熱Nu數隨著Pr數增加而增加,當Pr＞1時Nu數隨著Pr數變化不明顯;通過速度脈動信號和平均速度縱向分布最大值計算得到的Re〈u〉和ReUmax隨Pr數增加而減小,二者基本遵循同一標度律: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物理學報 20 204702]

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

[17]Bao Y,Ning H,Xu W 2014Acta Phys.Sin.63 154703(in Chinese)[包蕓,寧浩,徐煒2014物理學報 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力學學報 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