二維湍流熱對流羽流運動路徑對傳熱特性的影響*

包蕓 何建超 高振源

(中山大學航空航天學院,廣州 510275)

1 引 言

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

國內外關于RB熱對流的研究成果非常豐富,主要集中在傳熱特性變化規律以及溫度邊界層特性等問題.在理論研究方面,Grossman和Lohse[2,3]提出的GL理論最為成功,可以準確預測和描述大范圍內傳熱Nusselt(Nu)數與Rayleigh(Ra)數和Prandtl(Pr)數的關系.實驗研究中,He等[4]進行了Pr=0.8極高Ra數的系列實驗,發現當Ra數接近1015時,反映系統參數變化關系的標度律有所增大,與GL理論預測的極高Ra數的結果一致.Zhou和Xia[5]對RB系統的溫度邊界層進行研究,發現溫度邊界層分布與Prandtl-Blasius理論一致.Stevens等[6]通過直接數值模擬(DNS)方法對三維圓柱RB熱對流系統進行了模擬,Ra數達到2×1012.Zhu等[7]將二維湍流熱對流的DNS模擬Ra數提高到Ra=1014,并發現羽流發射特性在極高Ra數發生變化,引起傳熱增強.Stevens等[8]研究了湍流熱對流中熱分布的超結構.黃茂靜和包蕓[9]計算了二維和三維的情況,發現二維熱對流和三維展向平均熱對流的時間平均場中都存在大尺度環流和角渦,兩者溫度邊界層厚度關于Ra數的變化存在基本一致標度率關系.包蕓等[10]討論了Pr數對二維湍流熱對流流動和傳熱特性的影響.Bao等[11]和Chen等[12]研究了隔板對流系統,僅在對流槽中加入一定數量頂端留有狹縫的豎直隔板,結果發現可以成倍增強傳熱效率.林澤鵬和包蕓[13,14]對隔板對流系統的幾何參數、壓力特性等展開了系統的研究.Shishkina 等[15]考慮湍流脈動的影響,得到了湍流RB熱對流溫度邊界層方程.何鵬等[16]利用湍流溫度邊界層方程解,與不同Ra數和Pr數下二維DNS數值解溫度邊界層擬合,得到了很好的效果,并對擬合參數特性進行了討論.Gao等[17]在研究湍動能沿縱向分布時發現,二維和三維湍動能的差異僅發生在溫度邊界層內部.

Van der Poel等[18]對比了二維和三維熱對流DNS計算結果,二維計算 R a≤1010,發現三維結果與GL理論的預測吻合良好.由于二維的Nu數值均小于三維的結果,二維Nu數在Ra < 1010的范圍可以與向下平移的GL理論預測倍數線相符合,表明二維Nu數的變化規律與GL理論預測有一致性.但進一步增高Ra數時二維的傳熱結果偏離了GL理論預測倍數線.本文大范圍Ra數二維湍流熱對流的DNS計算結果研究發現,當Ra數進一步提高,傳熱Nu數隨Ra數的變化規律出現新的轉折,又回到了另一條GL理論預測倍數線上,并且這種傳熱特性的轉折變化規律與大尺度環流羽流運動路徑周長的變化特性相關聯.

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

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

其中V為無量綱速度矢量,θ為無量綱溫度,p為壓力,k為垂向單位矢量.無量綱參數Rayleigh(Ra)數為瑞利數,表征浮力驅動力與阻礙運動的力兩者相對大小.Prandtl(Pr)數為普朗特數,表征流體黏性耗散和熱耗散的關系.數值計算中邊界條件為壁面速度均采用無滑移條件,溫度為側壁采用絕熱條件,上下底板采用恒溫條件,上底板冷卻θ=—0.5,下底板加熱θ=0.5.

二維熱對流數值求解過程采用投影法.時間和空間均采用二階格式,顯式計算動量方程和溫度對流擴散方程,壓力泊松方程則需要全流場聯立求解.利用FFT變換解耦泊松方程,而后求解三對角方程組所建立的直接解法,在大規模湍流熱對流計算中必須通過并行計算進行.利用壓力泊松方程的高效并行直接求解PDD算法,聯合其他顯示容易并行計算的動量方程等,建立了熱對流DNS的并行直接求解方法(PDM-DNS,parallel direct method of DNS),并具有很好的并行效率[19].

3 高Ra數2D湍流熱對流傳熱特性

熱對流系統研究的核心問題是由熱浮力引起的湍流流動傳輸熱量的能力,用系統整體傳熱Nusselt數(Nu)表示.湍流熱對流傳熱特性最有成效的理論研究成果是GL理論[2,3],GL理論很好地預測了傳熱Nu數隨Ra數和Pr數的變化特性.Stevens等[20]通過近年大量新的實驗和計算數據,重新對GL理論結果中的參數進行了修正,使得在很大的Ra數范圍內Nu數對Ra數和Pr數的變化關系都能很好地滿足理論預測結果.在Ra數的大范圍中,傳熱Nu數隨Ra數的變化標度律并不是常值,是隨Ra數的變化而變化的.

當Ra≥1012以后,無論是實驗還是DNS數值模擬研究都碰到過極大的困難,困擾研究的發展.因此一般稱當Ra≥1012時為極高Ra情況.實驗研究比DNS模擬較早突破并進行極高Ra數湍流熱對流研究[4].Zhu等[7]在改進計算技術后將并行計算效率提高數十倍,從而在2018年實現了Ra=1014的二維湍流熱對流DNS模擬,所花計算資源超過500萬核時.

本文采用新的計算方法PDM-DNS,在“天河二號”超級計算機上進行了系列Ra數的二維湍流熱對流DNS模擬,對應兩組Pr=0.7和4.3.計算從Ra=107開始,最高Ra數達到Ra=1013.兩組計算的Ra數如表1所列.

表1 計算Ra數Table 1. The Ra numbers.

大量的計算結果都發現,二維熱對流的傳熱Nu數在同樣的Ra數要小于三維的結果,但具有與三維結果隨Ra數變化規律一致的特性.因此在探討二維湍流熱對流的傳熱特性時,用GL理論預測值乘以小數使其向下平移的倍數線,來反映二維傳熱特性的變化規律[18].

圖1中給出了本文計算的兩組Pr數情況下的二維湍流熱對流的傳熱Nu隨Ra數的變化情況,同時給出了四個三維計算結果.圖中顯示的是Nu/Ra0.3隨Ra數的變化,其中黑線是GL理論預測線,虛線是GL理論預測倍數線.紅色圓點是Pr=0.7的結果,藍色三角是Pr=4.3的計算結果.圖中還對應給出了Van der Poel等[18]和Zhang等[21]的二維計算結果以及Stevens等[6]的三維計算結果.

從圖1中可以看到,對Pr=0.7的二維結果,原本隨Ra數變化而增大的 Nu/Ra0.3在 Ra ＞109時開始減小,到Ra≈1010時達到最小,而后隨著Ra數的增加,Nu/Ra0.3隨Ra數的變化又開始相應地變大,回到另一條GL理論預測倍數線上,變化過程中出現兩個轉折.Pr=4.3的計算結果也出現同樣的變化特征,只是比Pr=0.7時對應的Ra數后移一個量級.在Ra ≤ 1010的范圍內Pr=4.3的Nu數與Ra數之間有較好的0.3標度律關系,當 Ra ＞ 1010后與 Van der Poel等[18]的二維結果一樣,Nu/Ra0.3隨Ra數的變化向下減小,但減小的幅度并沒有隨Ra數的增加而一直增大,當 Ra ≥ 2×1011時 Nu/Ra0.3的值又開始增加,并與GL理論預測倍數線變化走向一致.

荷里路德宮的一些建筑是三層的，一層是長長的連廊，二層和三層是房間。我驚奇地發現，宮殿的墻壁上居然有三種不同的裝飾柱式，爸爸說，一層為多立克柱式，二層為愛奧尼柱式，三層則為科林斯柱式。在同一棟建筑上同時呈現出三種柱式，這也算是對古希臘建筑藝術風格的致敬吧！

圖1 熱對流傳熱Nu數隨Ra數的變化情況Fig.1.The variation of the Nu number with Ra number for the turbulent convection.

圖1中二維湍流熱對流傳熱特性的這種在一定的Ra數范圍,Nu/Ra0.3隨Ra數變化的標度律相應減小而后增加的特殊情況,可能與二維湍流熱對流的流動特性有關.為此,我們從分析二維湍流熱對流的流場特性出發,探討二維湍流熱對流傳熱特性的特殊變化與二維湍流熱對流流動特征之間的相關關系.

4 高Ra數2D湍流熱對流流動特性

先以Pr=0.7為例,給出四個典型Ra數的二維湍流熱對流的瞬時溫度場以及溫度和速度的平均場,討論溫度羽流在不同Ra數的變化情況.

瞬時溫度場中的溫度分布可以反映冷熱羽流的狀態.圖2給出了Pr=0.7時四個典型Ra數瞬時溫度場,可以看到在不同Ra數時羽流的形態變化.給出的四個典型Ra數都已不存在較穩定的大尺度環流和角渦.在Ra=109,角渦脫離角區隨大尺度環流繞行,把羽流拉成條狀并隨機運動,并將部分羽流擠到方腔的中部,使得流動的最大速度不再沿側壁而向方腔中部偏移.Ra=1010時仍以條狀羽流為主,出現個別的渦狀羽流.到Ra=1011和極高Ra數Ra=1012,羽流基本上全部變為小尺寸的渦狀羽流,并相互纏繞,隨大尺度環流做宏觀的繞行.Ra數越高,渦狀羽流尺寸越小,數量越多,在流動中繞行的時間越長.

平均場特性將濾掉羽流運動的瞬時影響,給出總體傳熱特性和運動特征.圖3是上面四個典型Ra數的平均溫度場和流線圖.雖然四個不同Ra數的瞬時羽流形態和運動不一樣,但它們的平均場分布較為相似.近底板的溫度分布隨Ra數的增高而越來越接近底板,表明溫度邊界層越來越薄.流線圖形狀基本一致,平均場流動都接近圓形,有四個小角渦存在.但平均場流線圖并沒有反映出平均速度場的大小分布和變化情況.

圖4給出了對應的平均場的速度值分布.四個典型Ra數的平均場的速度值分布都基本為圓形,但速度值的最大值位置在不同Ra數時是不同的.Ra=1010速度值較大區域分布明顯向方腔中部聚攏,速度的最大值位置隨Ra數的變化是先向中部聚集再向外部擴展的過程.二維高Ra數的平均場速度分布的變化過程與瞬時羽流的運動形態相關聯.高Ra數時的角渦脫離角區,將羽流擠入方腔的中部,使得由溫度加速的流動最大速度分布區域也向方腔中部聚集.當Ra更高時羽流形態變化為小尺寸旋渦團狀羽流,團狀羽流隨大尺度環流繞行不再擠壓羽流運動,最大速度分布再次向方腔四周擴展.

圖2 Pr=0.7典型Ra數的二維熱對流瞬時溫度場Fig.2.The instantaneous temperature fields with typical Ra number at Pr=0.7.

圖3 Pr=0.7典型Ra數的二維熱對流平均溫度場和流線圖Fig.3.The average temperature fields and stream lines with typical Ra number at Pr=0.7.

圖4 Pr=0.7典型Ra數的平均速度場分布Fig.4.The average velocity field distributions with typical Ra number at Pr=0.7.

圖5 Ra=1012時不同Pr數的瞬時溫度場及羽流結構Fig.5.The instantaneous temperature field and plume structure with different Pr number at Ra=1012.

5 高Ra數2D湍流熱對流大尺度環流路徑與傳熱Nu數之間的相關性

定義大尺度環流路徑,用于描述2D湍流熱對流的平均場速度羽流運動路徑特性.在平均場速度分布中找到速度值最大的那點,畫出過這點的流線圖,并以此流線反映大尺度環流路徑[22].

圖6給出了兩個典型熱對流的大尺度環流路徑.可以看到,熱對流大尺度環流路徑在Ra=108是橢圓封閉曲線,路徑較長,而在Ra=1011是一個近似圓形的封閉流線,路徑較短.隨不同Ra數大尺度環流的路徑周長CLSC發生變化.

計算Pr=0.7時不同Ra數的大尺度環流路徑周長CLSC,并給出它們隨Ra數的變化,如圖7中上半部分的紅色圓點所示.圖中在Ra數較低時橢圓形大尺度環流路徑周長CLSC值較大.在Ra ≈109時有一個突然的減小,此時對應發生的是流態中角渦脫落羽流被擠向方腔的中心,大尺度環流變為圓形.在Ra ≈ 109至Ra ≈ 1010的區域大尺度環流路徑周長CLSC值較小.隨著Ra數的增加,羽流形態發生變化形成團狀羽流,大尺度環流路徑周長CLSC逐漸增大,到Ra=1013.

圖6 Pr=0.7時大尺度環流路徑圖Fig.6.The large scale circulation path at Pr=0.7.

圖7 大尺度環流路徑周長與傳熱Nu數隨Ra變化及其相關性,Pr=0.7Fig.7.The variation of large scale circulation path length and Nu number with Ra and its correlation at Pr=0.7.

討論Nu數隨Ra數的變化與大尺度環流路徑變化的關聯,將圖2中Pr=0.7時二維湍流熱對流的傳熱Nu/Ra0.3數隨Ra數的變化局部放大并入圖7中的下部.圖中可以看到兩者具有很好的相關性.Nu/Ra0.3數和大尺度環流路徑周長CLSC隨Ra數的變化都存在兩個明顯的轉折點.Nu/Ra0.3數隨Ra數的變化有個先增加再下降而后又增加的過程.當大尺度環流路徑周長CLSC在Ra ≈109突然變小時,Nu/Ra0.3數的變化出現減小的轉折.隨著Ra數的增加到Ra ≈ 1010,大尺度環流路徑周長CLSC開始增大,Nu/Ra0.3數的變化出現第二個轉折并開始增加,而且與GL理論預測0.72的倍數線符合良好.這一結果表明當Ra數很高時,二維湍流熱對流的傳熱特性變化規律仍可以與GL理論預測變化趨勢一致.

圖8中給出了Pr=4.3的二維湍流熱對流大尺度環流路徑周長CLSC隨Ra數的變化,同樣可以看到類似于圖7中Pr=0.7的結果,存在明顯的兩個轉折.

圖8 大尺度環流路徑周長與傳熱Nu數隨Ra變化及其相關性,Pr=4.3Fig.8.The variation of large scale circulation path and Nu number with Ra and its correlation at Pr=4.3.

圖8中同時給出了Van der Poel等[18]的二維計算結果,計算Pr=4.38,最大Ra=1011.Nu/Ra0.3數隨Ra數的變化局部放大后可以明顯的看到,當Ra數增高至Ra=1010時出現第一個轉折,二維的結果偏離了GL理論預測倍數線,使得高Ra二維的計算結果不再與GL理論預測變化趨勢一致,而且偏離情況隨Ra數越來越大,因此他們認為二維的計算結果在高Ra數時可能沒有意義.但他們沒有繼續提高Ra數的計算,錯過第二個轉折點的發現.

本文Pr=4.3的二維湍流熱對流DNS結果表明,在一定的高Ra數范圍內2D計算結果Nu/Ra0.3數隨Ra數的變化會出現向下偏離GL理論預測倍數線的現象.但同時可以看到,這個范圍的大尺度環流路徑周長也是減小的.繼續提高計算Ra數,當Ra ＞ 2X1011后大尺度環流路徑周長開始增大,Nu/Ra0.3數的變化出現第二個轉折,并且與GL理論預測0.7倍數線符合良好,又恢復到與GL理論預測變化基本一致.

Pr=4.3的兩個轉折點對應的Ra數與Pr=0.7時明顯增大,分別為 Ra ≈ 5X109和Ra ≈2X1011.這個變化與流動的湍流特性有關,Pr數越小流動的湍流程度越強[9],因此對于較小Pr數,熱對流流動更加容易發生角渦脫落和出現團狀羽流.不同Pr數二維的結果在第二個轉折點后相符合的GL理論預測倍數線略有不同,Pr=0.7時為0.72×GL理論預測值,Pr=4.3時為 0.70×GL理論預測值.Pr數較大GL理論預測倍數線略為下移.

以上研究結果可見,二維湍流熱對流DNS計算結果Nu數隨Ra數的變化,與三維的結果有些不同.當一定的Ra數范圍由于羽流的不規則運動而向方腔中部聚攏,使得大尺度環流路徑周長變小,傳熱Nu數隨Ra數的變化會向下偏離GL理論預測倍數線.其原因可能是與GL理論的四壁邊界層和中心區的假設不完全符合造成的.二維湍流熱對流傳熱特性隨Ra數變化出現兩個轉折現象的物理因素,還需要更多深入的研究.

綜上所述,二維湍流熱對流的計算結果在極高Ra數的情況下,傳熱Nu數隨Ra數的變化依然可以與GL理論預測倍數線符合良好,即保持與GL理論預測相同的變化規律.這對進一步開展極高Ra數的二維湍流熱對流DNS計算以及湍流熱對流“終極態”的研究,有著重要的意義.

6 結 論

本文采用新的計算方法完成了系列二維高和極高Ra數的湍流熱對流的DNS模擬.通過研究二維高和極高Ra數的傳熱變化特性,發現在一定的Ra數范圍Nu數隨Ra數的變化標度律減小,與GL理論預測倍數線偏離,但隨著Ra數進一步提高偏離現象消失,而這一傳熱特性隨Ra數的變化規律與大尺度環流路徑的大小值相關.通過研究得到以下結論.

1)創建了高效二維熱對流的并行直接求解方法(PDM-DNS),完成了系列不同Ra數的二維湍流熱對流DNS模擬,包括Ra=1013的DNS模擬.

2)二維湍流熱對流中,傳熱Nu/Ra0.3數隨Ra數的變化與反映羽流運動的大尺度環流路徑的變化具有很好的相關性,具有兩個轉折點.Pr=0.7時,Ra ≈ 109大尺度環流變為圓形,大尺度環流周長CLSC隨Ra數變化突然減小,出現第一轉折點,在Ra ≈ 1010時CLSC最小,出現第二轉折點,而后隨Ra增加變大.同樣傳熱特性隨Ra數的變化存在對應的兩個轉折.當Ra ＞ 109時Nu/Ra0.3隨Ra變化的標度律減小,出現偏離GL理論預測倍數線的現象.在Ra數大于第二轉折點Ra ≈1010,Nu/Ra0.3隨Ra變化的標度律變為增加并再次與GL理論預測倍數線相符合,直到Ra=1013.Pr=4.3時的結果與上述變化規律一致,只是對應轉折點的Ra數較高.

3)當Ra數大于第二轉折點時,2D湍流熱對流的傳熱Nu/Ra0.3隨Ra數的變化與GL理論預測倍數線符合良好,表明極高Ra數情況下2D熱對流的傳熱特性與GL理論預測變化趨勢基本一致.