側(cè)壁效應(yīng)對(duì)Rayleigh-Bénard對(duì)流影響的數(shù)值模擬研究

劉苗苗，王 啟，萬振華，孫德軍

(中國(guó)科學(xué)技術(shù)大學(xué) 近代力學(xué)系，合肥 230027)

0 引 言

由溫差導(dǎo)致的熱對(duì)流現(xiàn)象廣泛存在于自然界和各種工程應(yīng)用中，比如熱對(duì)流在海洋對(duì)流、大氣環(huán)流、地幔對(duì)流、太陽(yáng)對(duì)流等自然現(xiàn)象以及在電子設(shè)備的冷卻和凝固過程、核反應(yīng)堆設(shè)計(jì)、晶體生長(zhǎng)等工程應(yīng)用中扮演著重要角色。Rayleigh-Bénard（RB）對(duì)流模型是研究熱對(duì)流問題的經(jīng)典物理模型之一，RB對(duì)流指流體層中下壁面加熱、上壁面冷卻導(dǎo)致的流動(dòng)現(xiàn)象。該系統(tǒng)的控制參數(shù)主要是瑞利數(shù)Ra和普朗特?cái)?shù)Pr，響應(yīng)參數(shù)為努塞爾數(shù)Nu和雷諾數(shù)Re。當(dāng)Ra較小時(shí)，流動(dòng)由羽流主導(dǎo)，大尺度環(huán)流并不明顯[1]；隨著Ra的增加，系統(tǒng)逐漸進(jìn)入湍流狀態(tài)，此時(shí)上下壁面存在很薄的溫度邊界層，冷熱羽流在溫度邊界層內(nèi)生成并在浮力的作用下運(yùn)動(dòng)并自組織形成大尺度環(huán)流。隨著流動(dòng)復(fù)雜性的增加，對(duì)RB對(duì)流研究的重心也轉(zhuǎn)為對(duì)湍流狀態(tài)下系統(tǒng)各物理量統(tǒng)計(jì)特性的研究[2-4]。

在RB對(duì)流的研究中通常使用側(cè)壁絕熱的溫度邊界條件，在側(cè)壁絕熱的充分發(fā)展的RB湍流中，體區(qū)的溫度近似均勻，因此側(cè)壁的溫度在大部分區(qū)域近似為上下壁面溫度的平均值。實(shí)際的對(duì)流換熱裝置中，側(cè)壁的導(dǎo)熱性往往不能忽略[5-9]。極端條件下，當(dāng)側(cè)壁導(dǎo)熱系數(shù)相比流體很大時(shí)，可以認(rèn)為側(cè)壁是理想導(dǎo)熱壁面，此時(shí)側(cè)壁的溫度為線性分布[10-22]。

前 人 在 寬 度∶深 度∶高 度 之 比 為6∶2∶1[23]和6∶4∶1[24]的三維方腔內(nèi)研究了低Ra下不同側(cè)壁溫度對(duì)系統(tǒng)流動(dòng)結(jié)構(gòu)的影響。Corcione[25]在二維腔體研究了側(cè)壁溫度等于底壁溫度對(duì)流動(dòng)和傳熱的影響。Stevens等[6]研究了圓筒中不同側(cè)壁溫度對(duì)Nu的影響。Moon等[26]研究了側(cè)壁溫度等于上壁面溫度的立方腔體內(nèi)的流動(dòng)。Vasu等[27]研究了三維方腔中側(cè)壁等溫條件或等熱流條件對(duì)流動(dòng)模態(tài)的影響。也有研究討論了壁面厚度、導(dǎo)熱率對(duì)RB對(duì)流的影響，提出了側(cè)壁修正模型并進(jìn)行了分析[6-8,25,28-36]。工程應(yīng)用中熱對(duì)流的溫度邊界條件各種各樣，Pandey等[37]對(duì)方腔溫度邊界條件方面的研究進(jìn)行了詳細(xì)的總結(jié)，指出側(cè)壁溫度邊界條件主要分為絕熱、加熱、冷卻、溫度線性分布四種。目前，對(duì)側(cè)壁等溫且溫度為上下壁面溫度平均值的對(duì)流問題研究仍然不多。

本文使用直接數(shù)值模擬（DNS）研究了Pr=5.3，寬高比 Γ=1的二維方腔和三維方腔內(nèi)的RB對(duì)流。主要研究側(cè)壁等溫且溫度為上下壁面平均溫度的系統(tǒng)（ISO算例），并與側(cè)壁絕熱的經(jīng)典RB流動(dòng)（RB算例）進(jìn)行對(duì)比。在二維方腔中，研究發(fā)現(xiàn)，Ra比較小時(shí)上下壁面Nu不相等，這與系統(tǒng)中出現(xiàn)的大尺度羽流在中心區(qū)域擺動(dòng)的流動(dòng)結(jié)構(gòu)有關(guān)；而在三維方腔中，也發(fā)現(xiàn)了上下壁面Nu不相等的現(xiàn)象。文中也討論了側(cè)壁等溫條件對(duì)系統(tǒng)的熱量輸運(yùn)（Nu）和動(dòng)量輸運(yùn)（Re）的影響。

1 控制方程和數(shù)值方法

腔體的高度為， 寬度為，寬高比固定為1，三維腔體深度為， 且寬高比同樣固定為1。上壁面溫度為、 下壁面溫度為。在Oberbeck-Boussinesq近似下，對(duì)流介質(zhì)的物性參數(shù)認(rèn)為是常數(shù)，僅在浮力項(xiàng)中考慮密度變化且其中為下的密度和熱膨脹系數(shù)。笛卡爾坐標(biāo)系中無量綱化的不可壓縮Navier-Stokes方程為：

其中，u為 無量綱速度矢量；t為 時(shí)間；p為 壓力，T為無量綱溫度；k為豎直方向單位矢量；Pr=為普朗特?cái)?shù)，表征了流體的動(dòng)量擴(kuò)散速度與熱擴(kuò)散速度間的比值；無量綱參數(shù)為瑞利數(shù)，刻畫了溫差驅(qū)動(dòng)湍流流動(dòng)的強(qiáng)度，其中為重力加速度，為上下壁面的溫度差值，為運(yùn)動(dòng)學(xué)黏性系數(shù)，為熱擴(kuò)散系數(shù)。無量綱化選取自由落體速度作為參考速度，腔體高度作為參考長(zhǎng)度，參考時(shí)間為，溫度通過進(jìn)行無量綱化。在各個(gè)固壁施加無滑移和無穿透的速度邊界條件。

數(shù)值求解使用課題組自主開發(fā)的lMn2d和lMn3d程序。該程序采用交錯(cuò)網(wǎng)格，將矢量放置在網(wǎng)格界面上，將標(biāo)量放在網(wǎng)格中心處，避免了速度和壓力解耦?？臻g離散使用二階精度的中心差分，對(duì)非線性項(xiàng)采用Adams-Bashforth格式，對(duì)黏性項(xiàng)和擴(kuò)散項(xiàng)采用Crank-Nicolson格式。壁面附近采用拉伸網(wǎng)格以更好地分辨邊界層。在二維方腔和三維方腔內(nèi)分別在1×106≤Ra≤1×1010、1×105≤Ra≤1×108范圍內(nèi)計(jì)算了34個(gè)、20個(gè)算例。表1和表2給出了二維和三維計(jì)算采用的網(wǎng)格數(shù)。

表1 不同Ra下二維數(shù)值模擬的網(wǎng)格數(shù)Table 1 Grid resolutions for two-dimensional simulations with different Ra

表2 不同Ra下三維數(shù)值模擬的網(wǎng)格數(shù)Table 2 Grid resolutions for three-dimensional simulations with different Ra

2 結(jié)果與討論

2.1 方腔中側(cè)壁等溫對(duì)Nu數(shù)的影響

熱量輸運(yùn)規(guī)律是熱對(duì)流系統(tǒng)中關(guān)注的重要問題。無量綱數(shù)Nu表征了對(duì)流傳熱的效率，其計(jì)算公式為：

其中，是單位面積的平板傳遞給對(duì)流介質(zhì)的實(shí)際熱量，是熱傳導(dǎo)狀態(tài)下的熱通量，為熱傳導(dǎo)系數(shù)。為 對(duì)流導(dǎo)致的熱通量和傳導(dǎo)導(dǎo)致的熱通量之和，即：

則公式（2）可寫為：

圖1給出了二維和三維方腔各個(gè)壁面Nu、|Nub/Nut?1|以 及各個(gè)側(cè)壁Nu與Nut的比值隨Ra的變化情況。從圖1（a、e）中發(fā)現(xiàn)Nub>NuRB而且Nut>NuRB，說明在二維和三維方腔中側(cè)壁等溫邊界條件增強(qiáng)了系統(tǒng)的傳熱效率。而且發(fā)現(xiàn)隨著Ra的增加，Nub、Nut和NuRB之間的差別逐漸減小。這是因?yàn)楫?dāng)側(cè)壁溫度為上下壁面平均溫度時(shí)，一部分熱量可以直接通過側(cè)壁來傳遞，從而繞開溫度邊界層，這個(gè)過程增強(qiáng)了系統(tǒng)的傳熱效率。而隨著Ra的增加，溫度邊界層的厚度急劇減小，導(dǎo)致側(cè)壁與溫度邊界層之間的傳熱減弱，這使得兩個(gè)系統(tǒng)的Nu逐漸相等。

當(dāng)側(cè)壁溫度固定為上下壁面溫度的平均值時(shí)，由于系統(tǒng)存在的上下對(duì)稱性，直覺告訴我們上下壁面的Nu應(yīng)該相等，然而有意思的是，從圖1（a、e）可以看出二維和三維方腔中上下壁面的Nu在某些比較低的Ra區(qū)間都出現(xiàn)了明顯的偏差，即Nub≠Nut，這說明系統(tǒng)的上下對(duì)稱性被打破，系統(tǒng)的流動(dòng)結(jié)構(gòu)可能發(fā)生了變化。

圖1（b、f）分別給出了二維方腔和三維方腔中Nub和Nut之間相對(duì)誤差隨Ra的變化。在二維方腔中發(fā) 現(xiàn) |Nub/Nut?1|在Ra=1×106時(shí) 最 大，并 且 隨著Ra的增加不斷減小，當(dāng)Ra≥1×107時(shí)|Nub/Nut?1|<0.5%，誤差可以忽略不計(jì)。在三維方腔中，上下壁面Nu不相等情形發(fā)生在有限的Ra區(qū)間內(nèi)。從圖1（c、d）可以看出，在二維方腔中，當(dāng)上下Nu發(fā)生較大偏差時(shí)，左右側(cè)壁Nu絕對(duì)值、側(cè)壁Nu與Nut比值的絕對(duì)值較大，而且側(cè)壁Nu正負(fù)相同。從圖1（g、h）可知在三維方腔中，當(dāng)Ra很小，比如Ra=1×105時(shí) ，側(cè)壁Nu絕對(duì)值、側(cè)壁Nu與Nut比值的絕對(duì)值較大，而且存在兩對(duì)側(cè)壁Nu互為相反數(shù)，導(dǎo)致側(cè)壁Nu之和為0，最終Nub=Nut；隨著Ra的增加，當(dāng)Ra=7×105和Ra=1×106時(shí)，側(cè)壁Nu絕對(duì)值仍然較大，但是不再互為相反數(shù)，側(cè)壁Nu之和不為0，而且側(cè)壁Nu與Nut比值的絕對(duì)值較大，從而導(dǎo)致Nub≠Nut；隨著Ra進(jìn)一步增加，側(cè)壁Nu絕對(duì)值迅速減小，與Nut的比值的絕對(duì)值也迅速減小，使得上下壁面Nu恢復(fù)相等。

圖1 二維和三維方腔中各Nu隨Ra的變化以及Nu之間關(guān)系隨Ra的變化Fig. 1 Nu numbers and the relationships between them as functions of Ra in two- and three-dimensional square cavities

2.2 方腔中側(cè)壁等溫邊界條件對(duì)時(shí)均流場(chǎng)的影響

時(shí)間平均的流場(chǎng)可以反映系統(tǒng)的整體流動(dòng)特性。為了探究圖1中上下壁面時(shí)均Nu不相等的原因，在二維方腔中選取Ra=1×106和Ra=2×107兩個(gè)典型算例討論了時(shí)均溫度場(chǎng)以及側(cè)壁Nu局部分布，在三維方腔中選取Ra=1×105、Ra=3×105和Ra=1×106三個(gè)典型算例研究了水平截面上的豎直速度分布。

Ra=1×106

圖2給出了二維方腔中 時(shí)ISO系統(tǒng)和RB系統(tǒng)的時(shí)均溫度場(chǎng)和速度場(chǎng)，以及ISO系統(tǒng)時(shí)均流場(chǎng)左右側(cè)壁的局部Nu分布。通過圖2可以觀察到，ISO系統(tǒng)和RB系統(tǒng)的時(shí)均場(chǎng)的共同特點(diǎn)是，它們均呈現(xiàn)出四渦模態(tài)。然而ISO系統(tǒng)流場(chǎng)上下出現(xiàn)不對(duì)稱，而經(jīng)典RB系統(tǒng)的時(shí)均流場(chǎng)則呈現(xiàn)出完美的對(duì)稱四渦結(jié)構(gòu)。由圖2（c）可知，ISO系統(tǒng)時(shí)均流場(chǎng)左側(cè)壁局部Nu(NuLl) 和右側(cè)壁局部Nu(NuLr)分布完全重合，但是NuLl和NuLr上下不對(duì)稱，進(jìn)行積分后在左右側(cè)壁處得到Nul<0且Nur<0。系統(tǒng)存在能量的守恒，即Nub?Nut+Nul+Nur=0 ，最終導(dǎo)致Nub≠Nut。

圖2 二維數(shù)值模擬中R a=1×106時(shí)的時(shí)均溫度場(chǎng)分析（圖中黑色箭頭代表速度矢量）Fig. 2 Analysis of time-averaged temperature field when R a=1×106 in two-dimensional simulation(Black arrows represent velocity vectors)

圖3展示了Ra=2×107時(shí)ISO系統(tǒng)和經(jīng)典RB系統(tǒng)的時(shí)均流場(chǎng)以及ISO系統(tǒng)時(shí)均流場(chǎng)的左右側(cè)壁Nu分布。在計(jì)算過程中，兩類系統(tǒng)都沒有發(fā)生大尺度結(jié)構(gòu)的流動(dòng)反轉(zhuǎn)。通過圖3（a、b）可以發(fā)現(xiàn)兩個(gè)系統(tǒng)均存在一個(gè)大的環(huán)流結(jié)構(gòu)，且在對(duì)角線上存在一對(duì)角渦。圖3（c）中左右側(cè)壁局部Nu呈現(xiàn)出反對(duì)稱分布，即NuLl(z)=?NuLr(1?z)， 最終積分得到Nul=?Nur，側(cè)壁整體傳熱為0。由于系統(tǒng)能量平衡，所以上、下壁面的Nu必然相等。

圖3 二維數(shù)值模擬中R a=2×107 時(shí)的時(shí)均溫度場(chǎng)分析Fig. 3 Analysis of time-averaged temperature field when R a=2×107 in two-dimensional simulation

圖4展示了三維RB算例和ISO算例在計(jì)算的Ra區(qū)間內(nèi)定常、非定常的流動(dòng)狀態(tài)相圖?？梢园l(fā)現(xiàn)上下壁面Nu的不相等情形主要發(fā)生在側(cè)壁等溫系統(tǒng)剛剛由定常轉(zhuǎn)為非定常流動(dòng)的區(qū)間內(nèi)。

圖4 三維方腔內(nèi)RB對(duì)流的流動(dòng)狀態(tài)相圖Fig. 4 Phase diagram of flow states in three-dimensional RB convection

Pallares等[14]指出可以根據(jù)水平截面上豎直速度分布判斷系統(tǒng)內(nèi)的流動(dòng)結(jié)構(gòu)。圖5是Ra=1×105、Ra=3×105、Ra=1×106時(shí)三維ISO算例的時(shí)均流場(chǎng)在z= 0.15、0.5、0.85處的豎直速度分布，可看出系統(tǒng)始終存在單渦模態(tài)。觀察圖5（a）中z= 0.5平面內(nèi)的豎直速度分布，可以發(fā)現(xiàn)系統(tǒng)的定常上升流動(dòng)和下降流動(dòng)關(guān)于對(duì)角面對(duì)稱。隨著Ra增加，系統(tǒng)進(jìn)入非定常流動(dòng)狀態(tài)，圖5（b、c）中z= 0.5面上的豎直速度分布表明上升流和下降流之間失去了對(duì)稱性，該非對(duì)稱結(jié)構(gòu)導(dǎo)致側(cè)壁整體的Nu之和不為0，從而上下壁面的Nu不再相等，這與前面提到的二維情形比較類似。隨著Ra進(jìn)一步增加，側(cè)壁Nu迅速減小，相對(duì)于上下Nu可以忽略不計(jì)，此時(shí)上下壁面的Nu再次相等。

圖5 側(cè)壁等溫邊界條件的三維RB對(duì)流時(shí)均豎直方向速度場(chǎng)在水平截面高度為z = 0.15，z = 0.5和z = 0.85處的分布Fig. 5 Time-avraged vertical velocity for the three-dimensional RB convection with isothermal sidewalls at z = 0.15, 0.5, 0.85

2.3 腔側(cè)壁等溫對(duì)Re的影響

雷諾數(shù)是熱對(duì)流系統(tǒng)另一個(gè)重要的響應(yīng)參數(shù)，可以用于刻畫湍流的強(qiáng)度，其計(jì)算公式為：

我們分別使用ReISO和ReRB分別表示ISO系統(tǒng)和經(jīng)典RB系統(tǒng)的Re，探究側(cè)壁等溫邊界條件對(duì)Re的影響。圖6給出了二維和三維方腔流動(dòng)中ReISO、ReRB和ReISO/ReRB隨Ra的變化。通過圖6（b、d）可以發(fā)現(xiàn)ReISO/ReRB<1幾乎在所有計(jì)算Ra值區(qū)間下成立，只在二維方腔中Ra=3×107時(shí)大于1，而且與1僅相差0.90%，說明側(cè)壁等溫一般會(huì)減小系統(tǒng)Re。

圖6 Re和R eISO/ReRB隨著Ra的變化Fig. 6 Re and R eISO/ReRB as functions of Ra

2.4 二維方腔中側(cè)壁等溫對(duì)瞬時(shí)流場(chǎng)結(jié)構(gòu)的影響

本小節(jié)進(jìn)一步討論上下壁面Nu出現(xiàn)不相等時(shí)對(duì)應(yīng)的瞬時(shí)流動(dòng)狀態(tài)。角動(dòng)量的符號(hào)可以表征大尺度環(huán)流的瞬時(shí)流動(dòng)方向，其計(jì)算公式為：

其中 (v,w)是 速度矢量的分量，表示對(duì)整個(gè)計(jì)算區(qū)域進(jìn)行平均。

Ra=1×106

圖7展示了 時(shí)二維ISO系統(tǒng)和RB系統(tǒng)的角動(dòng)量變化以及典型的瞬時(shí)溫度場(chǎng)。如圖7（a）所示，可以看到經(jīng)典RB對(duì)流系統(tǒng)的角動(dòng)量隨時(shí)間發(fā)生頻繁的正負(fù)變化，說明系統(tǒng)在連續(xù)不斷的發(fā)生流動(dòng)反轉(zhuǎn)。圖7（b、c、d）展示了經(jīng)典RB對(duì)流系統(tǒng)發(fā)生反轉(zhuǎn)的連續(xù)過程，可以觀察到此時(shí)流動(dòng)由羽流主導(dǎo)，瞬時(shí)流場(chǎng)中存在大渦結(jié)構(gòu)，隨著角渦的增長(zhǎng)并發(fā)生融合，大尺度環(huán)流發(fā)生反轉(zhuǎn)。圖7（e）給出了這個(gè)反轉(zhuǎn)過程中角動(dòng)量的變化。圖7（f）展示了ISO系統(tǒng)的角動(dòng)量變化，發(fā)現(xiàn)角動(dòng)量同樣頻繁發(fā)生正負(fù)變化，但是在瞬時(shí)溫度場(chǎng)中沒有發(fā)現(xiàn)大尺度環(huán)流結(jié)構(gòu)以及流動(dòng)反轉(zhuǎn)，此時(shí)流動(dòng)也由羽流主導(dǎo)，沒有形成明顯的大尺度環(huán)流，而是出現(xiàn)如圖7（g～i）所示的三種典型的瞬時(shí)溫度場(chǎng)，分別為大尺度熱羽流占據(jù)中心區(qū)域并進(jìn)行擺動(dòng)的流場(chǎng)、大尺度冷羽流占據(jù)腔體中心區(qū)域并進(jìn)行擺動(dòng)的流場(chǎng)和冷熱羽流無序發(fā)射的流場(chǎng)。

圖7 二維計(jì)算中R a=1×106時(shí)RB算例和ISO算例的角動(dòng)量隨時(shí)間的演化以及典型瞬時(shí)溫度場(chǎng)，其中（b～d）為RB算例發(fā)生反轉(zhuǎn)時(shí)不同時(shí)刻的瞬時(shí)流場(chǎng)，（g～i）為ISO算例三個(gè)不同時(shí)刻的典型瞬時(shí)溫度場(chǎng)Fig. 7 Evolution of the angular momentum and typical instantaneous temperature fields for the adiabatic-sidewall case and the isothermal-sidewall case with R a=1×106 in 2D simulations. Figs. (b～d) instantaneous temperature fields during a flow reversal process for the adiabatic-sidewall case. Figs. (g～i) three typical instantaneous temperature fields for the isothermal-sidewall case

當(dāng)系統(tǒng)處于大尺度熱羽流占據(jù)中心區(qū)域并進(jìn)行擺動(dòng)的流動(dòng)結(jié)構(gòu)時(shí)，下壁面Nu小于上壁面Nu；當(dāng)系統(tǒng)處于大尺度冷羽流占據(jù)中心區(qū)域并進(jìn)行擺動(dòng)的流動(dòng)結(jié)構(gòu)時(shí)，上壁面Nu小于下壁面Nu；當(dāng)系統(tǒng)處于冷熱羽流無序發(fā)射的流場(chǎng)結(jié)構(gòu)時(shí)，上下壁面Nu幾乎相等。由于三種流動(dòng)結(jié)構(gòu)持續(xù)的時(shí)間不同，導(dǎo)致了時(shí)均Nu和時(shí)均流場(chǎng)的不對(duì)稱。而隨著Ra的逐漸增加，系統(tǒng)流動(dòng)加強(qiáng)，瞬時(shí)流場(chǎng)呈現(xiàn)出穩(wěn)定的大渦結(jié)構(gòu)，上述三種結(jié)構(gòu)消失，此時(shí)上下壁面的Nu相等。

3 結(jié) 論

本文結(jié)合二維和三維直接數(shù)值模擬方法研究了側(cè)壁恒溫且溫度為上下壁面平均值的邊界條件對(duì)RB對(duì)流的影響，并與傳統(tǒng)的側(cè)壁絕熱條件下的RB對(duì)流進(jìn)行了比較。寬高比 Γ固定為1，普朗特?cái)?shù)固定為Pr= 5.3，二維和三維研究的Ra區(qū)間分別為1×106≤Ra≤1×1010、1 ×105≤Ra≤1×108。主要結(jié)論歸結(jié)如下：

1）Ra較小時(shí)，在特定Ra區(qū)間內(nèi)存在上下壁面Nu不相等的現(xiàn)象。該現(xiàn)象出現(xiàn)時(shí)，流場(chǎng)結(jié)構(gòu)存在明顯的上下不對(duì)稱，導(dǎo)致側(cè)壁的Nu之和并不為0，從而上下壁面的Nu不相等。隨著Ra的增加，側(cè)壁的Nu減小且可以忽略，此時(shí)上下壁面的Nu恢復(fù)成相等。

2）側(cè)壁恒溫時(shí)系統(tǒng)的傳熱效率要高于側(cè)壁絕熱情況，這是由于側(cè)壁恒溫條件下，一部分熱量可以直接通過側(cè)壁來傳遞，從而繞開溫度邊界層。

3）當(dāng)采用恒溫邊界條件時(shí)，系統(tǒng)的Re要小于側(cè)壁絕熱條件下對(duì)應(yīng)的Re，從而表明側(cè)壁恒溫邊界條件會(huì)降低系統(tǒng)的流動(dòng)強(qiáng)度。