胡 瑾，章盛祺，夏振華，*

(1. 浙江大學 流體動力與機電系統(tǒng)國家重點實驗室，杭州 310027；2. 浙江大學 航空航天學院，杭州 310027；3. 南方科技大學 廣東省湍流基礎研究與應用重點實驗室，深圳 518055；4. 北京大學 湍流與復雜系統(tǒng)國家重點實驗室，北京 100871)

0 引 言

熱對流是一種常見的物理現(xiàn)象，與人們的日常生活、生產(chǎn)息息相關[1]。其不僅存在于北極夏季融冰湖[2]、大氣及海洋環(huán)流[3-4]、地球內(nèi)部地幔形成[5]等自然現(xiàn)象中；還廣泛應用于核反應堆冷卻[6]、新型儲能材料[7]、芯片散熱[8]、磁流體傳熱[9]等工業(yè)工程問題中。而Rayleigh-Bénard對流（Rayleigh-Bénard convection，RBC）系統(tǒng)，則是從眾多復雜的自然現(xiàn)象及工程問題中抽象簡化出來的研究對流問題的經(jīng)典模型。同時，由于RBC系統(tǒng)中可以同時存在多種熱量傳遞方式，因此它也是傳熱領域非常經(jīng)典的課題，研究RBC系統(tǒng)對認識實際情況下各種各樣對流過程中的傳熱機制有著非常重要的意義。

狹義上的RBC通常指在一個封閉腔體內(nèi)，上表面冷卻，下表面加熱，形成恒定溫差，從而導致腔體內(nèi)流體產(chǎn)生特殊流動的系統(tǒng)，如圖1（來源van der Poel等[10]）所示。當系統(tǒng)充分發(fā)展后，上下壁面附近會形成溫度邊界層[11]，同時不斷生成冷熱羽流，并在腔體中形成穩(wěn)定的大尺度環(huán)流（large-scale circulation，LSC）結(jié)構。除了對經(jīng)典標度律[12-15]、擬序結(jié)構[16-17]及能量耗散[18-19]等的研究，科研工作者在近二十年對LSC的反轉(zhuǎn)問題也進行了深入的研究[20-22]。最近Jiang等[23]通過設置上下壁面的非對稱棘輪狀粗糙度和引入較小傾角對熱對流系統(tǒng)中的LSC和傳熱進行控制；Zhang等[24-25]通過增加側(cè)壁局部恒溫條件也實現(xiàn)了對LSC反轉(zhuǎn)及強弱的有效控制。

圖1 RBC系統(tǒng)示意圖[10]Fig. 1 Sketch of RBC system[10]

由于經(jīng)典的RBC系統(tǒng)假設上下邊界溫度恒定且分布均勻，這與一些實際系統(tǒng)的熱邊界條件并不相符。例如在地球物理研究中模擬地幔對流時，需要考慮熱傳導性差異較大的大陸板塊與海洋板塊，此時熱邊界條件就不再均勻。另一方面，工程應用中一些對流換熱優(yōu)化問題也需要考慮非均勻熱邊界條件。在對流換熱的優(yōu)化問題研究中，已有的增加對流傳熱效率的途徑主要分為兩種：一是改變換熱板表面的粗糙度，該方法已被證明可在有限Rayleigh數(shù)（Ra）范圍內(nèi)顯著增強對流熱流[26-29]；二是采用等溫和絕熱相結(jié)合的混合熱邊界條件，已有的數(shù)值和實驗研究結(jié)果表明，系統(tǒng)的傳熱效率會隨絕熱面積的增大而減小[30]，當絕熱面積超過50%時不同的絕熱位置設置會對熱傳遞產(chǎn)生顯著差異[31]。除此之外，也有科研人員發(fā)現(xiàn)，在非均勻熱邊界條件下，當單個絕熱塊特征尺度與溫度邊界層厚度相近時，其Nusselt數(shù)（Nu）會與經(jīng)典RBC系統(tǒng)相近[32-33]。近來，Vasiliev等[34]研究了Prandtl數(shù)Pr=6.46、Ra=1×107時僅下壁面局部加熱條件對方腔內(nèi)熱湍流系統(tǒng)的流動結(jié)構及傳熱規(guī)律的影響，其發(fā)現(xiàn)加熱面積相同時，LSC及形態(tài)與加熱塊數(shù)有較強的關聯(lián)，傳熱效率隨加熱塊數(shù)的增加而增大；Nandukumar等[35]研究了Pr= 1時，上下壁面均為一半加熱一半冷卻條件下熱對流系統(tǒng)與經(jīng)典RBC系統(tǒng)的關系，其發(fā)現(xiàn)此非均勻熱邊界條件在1×105＜Ra＜1×108時可明顯增加系統(tǒng)傳熱效率。由此可知，研究非均勻熱邊界條件下熱對流的流場及傳熱規(guī)律是更切合實際和富有價值的。

但目前已有文獻中非均勻熱邊界條件主要考慮簡單對稱分布情況，所以進一步研究非對稱分布的復雜熱邊界條件對高Rayleigh數(shù)熱湍流系統(tǒng)流動結(jié)構及傳熱規(guī)律的影響是必要的。本文參考Zhang等[24-25]算例設置，基于Pr=2、Ra=1×108的經(jīng)典二維RBC系統(tǒng)，改變其下壁面加熱條件，初步研究了加熱長度為方腔邊長一半時，不同局部加熱位置對系統(tǒng)LSC結(jié)構和傳熱效率的影響。

1 控制方程及數(shù)值設置

1.1 控制方程

假設方腔高為H，滿足寬高比 Γ=1；u=(u,v)，其中u、v分別為水平方向及豎直方向速度分量；θ為溫度，θl和 θu分別為下壁面加熱區(qū)域及上壁面的恒定溫度； Δθ=θl?θu為 上下壁面溫差；給定運動黏度υ，熱擴散系數(shù)κ，導熱系數(shù)λ，重力加速度g，熱膨脹系數(shù)β；坐標原點定義在方腔左下角。引入特征速度（自由落體速度）特征時間T=H/U、特征長度H、特征溫度 Δθ，對控制方程進行無量綱化。假設系統(tǒng)滿足Boussinesq近似，則可得無量綱控制方程組：

其中Ra和Pr為RBC系統(tǒng)的無量綱控制參數(shù)，其定義如下：

后文中變量如無特別說明，均已無量綱化。圖2為本文計算的三種局部加熱系統(tǒng)及RBC系統(tǒng)示意圖。左右壁面為絕熱條件，即 ?θ/?x= 0；上壁面恒溫冷卻，即 θ=0；下壁面非加熱區(qū)域滿足絕熱條件，即?θ/?y= 0，加熱區(qū)域滿足 θ=1。

圖2 局部加熱系統(tǒng)及RBC系統(tǒng)計算示意圖Fig. 2 Sketches of local heating and RBC cases

1.2 數(shù)值方法

本文數(shù)值模擬均采用經(jīng)過部分修改的二階有限差分法代碼AFiD[10]進行計算，其中離散Poisson方程通過水平方向的離散余弦變換解耦并采用三對角求解器進行求解，時間推進采用二階顯式Adams-Bashforth格式實現(xiàn)。與de Vahl Davis等[36]的方腔自然對流結(jié)果對比，文獻在Pr= 0.71、Ra= 1×104、1×105、1×106條件下的壁面時均Nu數(shù)為2.243、4.519、8.800，而本文所用程序計算的結(jié)果為2.256、4.537、8.839；與Zhang等[25]的經(jīng)典RBC系統(tǒng)結(jié)果對比，文獻在Pr=2、Ra= 1×108條件下的壁面時均Nu數(shù)為24.90，而本文所用程序計算的結(jié)果為24.92?？梢哉J為與已有文獻結(jié)果基本一致，證明了現(xiàn)有代碼的準確性。

1.3 算例設置

本文在Pr= 2、Ra= 1×108條件下，以加熱左端點無量綱坐標X為位置參數(shù)，對全局加熱（經(jīng)典RBC）系統(tǒng)和局部加熱系統(tǒng)（加熱位置分別為X= 0、0.125、0.25，加熱長度l=0.5）進行統(tǒng)一的分析和比較。

所有算例中，水平方向為均勻網(wǎng)格，豎直方向為非均勻網(wǎng)格（靠近壁面處細化），保證邊界層內(nèi)網(wǎng)格大小Δ＜0.6min[ηK,ηB][37]；時間步長足夠小，足以確保Courant-Friedrichs-Lewy數(shù)小于0.2；達到統(tǒng)計穩(wěn)態(tài)后計算足夠長時間，RBC系統(tǒng)可觀察到反轉(zhuǎn)現(xiàn)象。文中 〈···〉表示對括號內(nèi)變量求平均，下標表示如何平均。

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

2.1 局部加熱條件對流場變化的影響

首先考慮不同位置局部加熱條件對熱湍流達到統(tǒng)計穩(wěn)態(tài)后的流動影響，分析總動能的時間演化情況。本文計算總動能為：

其中下角標S表示對整個計算域面積求平均。通過計算可得不同加熱位置及經(jīng)典RBC的Ek(t)在達到統(tǒng)計穩(wěn)態(tài)后的變化情況，如圖3所示。

圖3 不同加熱位置及經(jīng)典RBC的Ek(t)變化情況Fig. 3 Time series of Ek(t) of local heating cases and traditional RBC case

其中RBC系統(tǒng)的總動能均值 〈Ek〉t=0.0200，最大值和最小值之間相差至少0.01，明顯比局部加熱系統(tǒng)的均值及振幅要大。但不同位置的局部加熱系統(tǒng)之間的均值差異并不明顯，因此我們進一步計算X=0、X=0.125、X=0.25三 個系統(tǒng)的 〈Ek〉t，可得其值分別為0.012 5、0.013 9、0.014 5，約 為RBC系統(tǒng)的62.5% ～72.5%。由于局部加熱系統(tǒng)三條曲線重合率較高，故由加熱位置從左至右分別給出系統(tǒng)Ek(t)的方差，結(jié)果為7.56×10?4、6.67×10?4、5.82×10?4，依次減小。綜上所述，我們可以認為，局部加熱系統(tǒng)中加熱位置越靠下壁面中心，系統(tǒng)總動能Ek(t)越大、其對應振幅越小。

根據(jù)前文所提，傳統(tǒng)RBC具有的流場結(jié)構特性之一便是當流動充分發(fā)展后，流場內(nèi)會形成較為穩(wěn)定的LSC結(jié)構且會發(fā)生反轉(zhuǎn)現(xiàn)象。為了分析局部加熱是否會影響LSC的反轉(zhuǎn)，本文采用角動量符號對此進行判斷。計算平均角動量定義如下[38]：

圖4給出局部加熱系統(tǒng)及RBC的L(t) 在1×104≤t≤4×104范圍內(nèi)的演化情況，其中L<0表示LSC呈順時針狀態(tài)，L>0為逆時針。從圖4中可以明顯地看到經(jīng)典RBC存在大渦反轉(zhuǎn)現(xiàn)象，而局部加熱系統(tǒng)均未出現(xiàn)反轉(zhuǎn)情況且LSC呈順時針狀態(tài)。放大局部加熱系統(tǒng)的L(t)演化情況，如圖5所示。結(jié)合圖5，進一步計算局部加熱系統(tǒng)達到統(tǒng)計穩(wěn)態(tài)后角動量均值，可得加熱位置從左壁面至下壁面中心 〈L〉t分別為?0.049 9、?0.052 6、?0.054 0，方差Lrms分別為0.002 4、0.001 8、0.001 5，由此可知局部加熱系統(tǒng)中加熱位置越靠近下壁面中心，其L(t)絕對值越大、振幅越小，系統(tǒng)主體LSC穩(wěn)定性越高，LSC反轉(zhuǎn)可能性越低。

圖4 不同加熱位置及經(jīng)典RBC的L(t)變化情況Fig. 4 Time series of L(t) of local heating cases and traditional RBC case

圖5 不同加熱位置的L(t)變化情況Fig. 5 Time series of L(t) of local heating cases

Sugiyama等[38]曾描述了角渦在LSC反轉(zhuǎn)過程中的作用，即由于角落邊界層分離羽流的能量補給，小的角渦會逐漸增大，直至達到腔體高度的一半，然后破壞主體LSC結(jié)構，并在相反方向上建立另一個新的LSC。由此我們知道反轉(zhuǎn)現(xiàn)象的產(chǎn)生與角渦的生長密不可分，為了觀察角渦及中心大渦的變化情況，繪制了流場云圖及動畫（在網(wǎng)刊資源附件中提供）。從動畫及圖6中可看出無論是否局部加熱，其基本流場結(jié)構均主要由中心大渦及角渦組成。其中局部加熱系統(tǒng)中下壁面角渦不穩(wěn)定，且大小受限。雖然RBC系統(tǒng)中的角渦關于方腔中心的對稱性更強，θ、u、v的數(shù)值分布均具有更高的對稱性，但其中心大渦也被更劇烈地擠壓，更容易誘導渦的破碎和反轉(zhuǎn)。而在局部加熱系統(tǒng)演化過程中，系統(tǒng)熱羽流只能沿某一（距離加熱區(qū)域更近的）側(cè)壁上升，遇到下降的冷羽流時發(fā)生分離，在上壁面附近形成與RBC系統(tǒng)中相似的正常形態(tài)角渦，而在另一側(cè)由于缺少上升的熱羽流，只能在下壁面附近產(chǎn)生僅由冷羽流誘導的角渦。此時系統(tǒng)中下壁面角渦由于缺乏能量補給條件，顯然無法正常生長，使得系統(tǒng)LSC無法發(fā)生反轉(zhuǎn)。我們猜測，本文研究的局部加熱系統(tǒng)里由于下壁面的角渦缺乏持續(xù)的冷熱羽流相互競爭機制，導致該角渦大小無法突破LSC反轉(zhuǎn)的閾值，即達到腔體高度的一半，從而無法發(fā)生LSC反轉(zhuǎn)，這與Sugiyama等[38]描述的機制一致。此外，局部加熱系統(tǒng)中角渦的大小雖然受限制，但仍然會隨著時間演化不斷波動。

圖6 不同加熱位置及經(jīng)典RBC瞬時溫度云圖Fig. 6 Temperature contours of local heating cases and traditional RBC case

2.2 不同位置的局部加熱條件對傳熱效率的影響

由于左右壁面均為絕熱條件，故只考慮豎直方向傳熱情況，通過Nu數(shù)對其進行量化，本文定義有效Nu數(shù)為：

其中下邊界絕熱壁面 ?θ/?y=0，由于進出熱流的平衡，滿足下邊界加熱區(qū)域平均熱流密度為上邊界的兩倍，故根據(jù)上述公式仍能得到可靠的系統(tǒng)平均Nu數(shù)。圖7給出了系統(tǒng)達到統(tǒng)計穩(wěn)態(tài)后，不同加熱位置及經(jīng)典RBC系統(tǒng)的Nu(t)隨時間演化情況。為了更直觀地展示和對比，我們將RBC系統(tǒng)的Nu(t)上移4，而將X=0.125和X=0.25的Nu(t)下移了8和16。

圖7 不同加熱位置及經(jīng)典RBC的上下壁面Nu(t)變化情況Fig. 7 Time series of Nu(t) at the bottom wall (a) and top wall (b)from different cases

其 中RBC系 統(tǒng) 〈Nu〉t=24.92，如 上 文 所 提，與Zhang等[25]文獻給出的24.90基本一致。從圖中可以明顯看出的是，RBC系統(tǒng)的 〈Nu〉t比局部加熱系統(tǒng)振蕩更為劇烈，而局部加熱系統(tǒng)內(nèi)部則是隨著加熱位置越靠近中心而越小。為了得到定量的數(shù)據(jù)，我們計算X=0、X=0.125、X=0.25三個局部加熱系統(tǒng)在統(tǒng)計穩(wěn)態(tài)的 〈Nu〉t，結(jié)果分別為16.97、17.46、18.23。可以發(fā)現(xiàn)，雖然局部加熱系統(tǒng) 〈Nu〉t明顯小于RBC系統(tǒng)，但隨著加熱位置越靠近下壁面中心， 〈Nu〉t的值也會隨之越來越大，傳熱效率依次升高，這與前文總動能 〈Ek〉t、角動量絕對值 |〈L〉t|的相對大小保持一致。其中雖然局部加熱系統(tǒng)的加熱長度l=0.5，僅為RBC系統(tǒng)的一半，但其 〈Nu〉t可以達到傳統(tǒng)RBC的68.1% ～ 73.2%，遠高于設想的50%。而且通過比較分析，我們發(fā)現(xiàn)X=0.25系 統(tǒng) 〈Nu〉t比X=0.125系 統(tǒng) 增 長 了4.4%，比X=0系統(tǒng)增長了7.4%。也就是說，在局部加熱區(qū)域長度確定的情況下，除了Bakhuis等[33]通過將導熱及絕熱區(qū)域分布為數(shù)量盡量多的條紋圖案以提升系統(tǒng)傳熱效率外，我們也可以通過控制整體加熱位置，來使系統(tǒng)傳熱效率最大化。

觀察圖7（b）并與圖7（a）作對比，可以看出在達到統(tǒng)計穩(wěn)態(tài)后，RBC系統(tǒng)上下壁面Nu(t)特征相似；在局部加熱系統(tǒng)中，雖然上壁面Nu(t)振幅區(qū)別不大，但加熱位置越靠近下壁面中心，下壁面Nu(t)幅值越小導致上下壁面之間差距越明顯。這是因為上壁面附近以側(cè)壁分離出的冷熱羽流共同作用為主，對其影響相近。而另一側(cè)由于沒有熱羽流的阻礙，冷羽流下降后會先到達下壁面而后橫向運動，再與相遇的熱羽流一同上升。結(jié)合云圖及動畫可知，加熱位置越靠近中心，即加熱區(qū)域離冷羽流到達下邊界位置越近，受其沖擊作用和橫向運動影響的范圍就越大，從而導致自由生長的熱羽流數(shù)量越少，因此加熱區(qū)域溫度脈動就越小，Nu(t)振幅也越小。

圖8展示了統(tǒng)計穩(wěn)態(tài)時，下壁面加熱區(qū)域的時均溫度梯度分布，其證明加熱位置越靠近冷羽流沖擊點，其受冷羽流影響端溫度梯度就越大，即冷羽流作用程度越強。而中心加熱系統(tǒng)中，加熱位置最靠近冷羽流發(fā)生變向運動區(qū)域，其受影響最大，熱羽流生成數(shù)量始終維持在1～2個，因此熱羽流可以以近似一簇的狀態(tài)快速從左側(cè)壁上升進而與冷羽流進行相互作用，又因為其對于方腔中心而言具有更穩(wěn)定的LSC結(jié)構，因此具有更高效率的對流換熱，故整體傳熱效率更高。

圖8 不同加熱位置系統(tǒng)的下壁面時均溫度梯度Fig. 8 Temperature gradient of the bottom wall of different heating position cases

3 結(jié) 論

本文以二維方腔RBC系統(tǒng)為基礎，通過直接數(shù)值模擬計算了Pr=2、Ra=1×108條件下，加熱區(qū)域長度l=0.5，加熱位置分別為X= 0、0.125、 0.25的三種局部加熱系統(tǒng)。分別研究了不同位置的局部加熱條件對系統(tǒng)總動能、角動量、LSC反轉(zhuǎn)及傳熱效率的影響。對結(jié)果進行分析，得出以下結(jié)論：

1）局部加熱系統(tǒng)中，加熱位置越靠近下壁面中心，其總動能和角動量絕對值越大、振幅越小，系統(tǒng)主體LSC穩(wěn)定性越高；

2）局部加熱條件會對系統(tǒng)下壁面角渦的生長產(chǎn)生限制，進而抑制LSC反轉(zhuǎn)；

3）雖然加熱長度為RBC系統(tǒng)的一半，但局部加熱系統(tǒng)總動能〈Ek〉t可達到后者的62.5%～72.5%，傳熱效率可達到68.1%～73.2%；

4）在一定的加熱長度下，可以通過調(diào)整加熱位置令傳熱效率最大化，其中X=0.25系 統(tǒng)就比X=0系統(tǒng)增長了7.4%；

5）在沒有熱羽流上升側(cè)，冷羽流下降后到達壁面處改為橫向運動，該位置與加熱區(qū)域越近則對熱羽流數(shù)量和對應溫度脈動影響越明顯，其對應端點溫度梯度也越大。

由于本文計算選取局部加熱位置數(shù)量有限，因此想要得到更為確切的定量關系，還需大量不同加熱位置及加熱長度的計算結(jié)果。本文初步證明了加熱長度確定的情況下，可以通過加熱位置來控制傳熱效率。關于局部加熱條件對LSC反轉(zhuǎn)的影響，后續(xù)還會深入研究。除此之外，本文作為對非對稱局部加熱問題的初步研究，目前僅開展了二維直接數(shù)值模擬研究，后續(xù)我們將進一步開展三維問題的直接數(shù)值模擬與分析，以及相應的實驗研究。