陽建林 張義招 周 全

(1.上海大學力學與工程科學學院, 上海 200444;2.上海大學上海市應用數學和力學研究所, 上海 200072)

湍流作為一種流體流動現象存在于自然界中, 我們生產生活中見到的流體流動絕大部分都處于湍流狀態(tài).Rayleigh-B′enard(RB)湍流熱對流是抽象出來的研究對流問題的經典流體力學模型, 研究RB 湍流熱對流可以理解湍流運動的物理本質.RB湍流熱對流的研究最早可以追溯到20 世紀初, 1900 年B′enard[1]進行了實驗研究, 1916 年Rayleigh 等[2]在理論上對浮力驅動的分層流動進行了描述.RB系統(tǒng)由一個封閉的對流槽構成, 對流槽頂部和底部放置有可以傳導熱量的導板, 對流槽內部充滿流體介質.對底部導板加熱并對頂部導板進行冷卻, 底部流體介質受熱膨脹, 頂部流體介質遇冷收縮, 就會產生密度差.在浮力和重力的作用下, 底部熱流體上升, 頂部冷流體下降, 開始流動.根據上下導板的溫度差?T不同, 對流槽內流體會出現不同流動狀態(tài).隨著?T增大, 對流槽內流體會由靜止發(fā)展為規(guī)則的對流狀態(tài);當?T繼續(xù)增大到一定程度時, 對流槽內流體的運動轉變?yōu)橥牧鳡顟B(tài).在充分發(fā)展的湍流熱對流系統(tǒng)中, 冷熱羽流分別從上下導板中生成、分離, 并在運動過程中通過自組織運動形成大尺度環(huán)流[3].

1 RB 湍流熱對流

1.1 控制方程、控制參數和響應參數

RB 湍流熱對流是一個復雜的強非線性系統(tǒng), 流體的溫度與速度耦合, 溫度變化會引起密度的變化.在流體溫度梯度變化不大的情況下, 通過引入Oberbeck-Boussinesq 近似[4-5]簡化模型, RB 系統(tǒng)的控制方程組使用特征長度H、特征速度κ/H、特征溫度T、特征時間H2/κ無量綱后表示如下:

式中:u是速度矢量,θ是溫度,p是壓力,z為豎直方向的單位矢量.式(1)～(3)為含有熱浮力項的Navier-Stokes(NS)方程連續(xù)性方程和熱運輸方程.由上述無量綱化的方程組可以得出RB 系統(tǒng)的兩個重要的控制參數Rayleigh 數(Ra)和Prandtl 數(Pr), 分別為

式中:α為膨脹系數;g為重力加速度;ν為流體黏性系數;κ為熱擴散系數.Ra數表示驅動流體流動的熱浮力與抑制流體運動的流體耗散力的比值, 與對流槽的上下板的溫差?T和對流槽高度H的三次方成正比.Pr數代表流體本身的物理性質, 是動量擴散率與熱擴散率之比,也反映黏性邊界層厚度δu和溫度邊界層厚度δth之間的關系.流槽的幾何形狀采用寬高比Γ描述,

式中:D為對流槽的水平方向長度;H為豎直方向長度.

RB 對流系統(tǒng)還有兩個重要的響應參數, 分別為Nusselt 數(Nu)和Reynolds 數(Re),

式中:Nu用于描述系統(tǒng)的傳熱效率;λ是流體的熱傳導系數;Q為通過對流槽的熱通量;Re表示流體所受的慣性力和黏性力之比, 是系統(tǒng)的無量綱速度;U可以是大尺度環(huán)流的特征速度,也可以是流場中速度均方根的值.

1.2 研究現狀

粗糙邊界作為一個增強傳熱的方法廣泛應用于湍流熱對流中, Shen 等[6]的實驗以水為對流介質, 在對流槽上下導板表面均勻地排布了金字塔狀的粗糙元, 發(fā)現當Ra數超過某一臨界值時, 系統(tǒng)的整體熱輸運效率可以增強約20%.Du 等[7-8]的工作更進一步地表明, 大尺度環(huán)流引發(fā)的強水平剪切流動與粗糙元相互作用會在粗糙元之間形成二次渦, 有利于羽流的形成, 從而可以極大地提高系統(tǒng)的整體熱輸運效率, 實驗測得Nu數增量超過76%.Shishkina 等[9]對底部粗糙元剖面為長方形的方形對流槽進行數值模擬, 并基于二維Prandtl-Blasius 邊界層方程提出了一種計算粗糙元引起的Nu數偏差的分析模型.Liot 等[10]對頂板保持光滑, 下底板布置有方形粗糙元的方形對象槽進行實驗, 值得一提的是此工作中的實驗裝置非常大, 高度和寬度達到2.5 m, 厚度為0.5 m, 能夠測得非常精細的速度場, 研究發(fā)現粗糙元之間流體存在外部對流和流體交換時能顯著增強傳熱.Xie 等[11]對底部粗糙元為金字塔型的圓柱形對流槽進行實驗, 通過定義粗糙元高寬比λ, 7.5×107≤Ra≤1.31×109,Pr數從3.57 增加至23.34, 發(fā)現隨著λ從0.5 增加至4,Nu和Ra的標度律關系大致可以分為3 個區(qū)域(區(qū)域一標度律和光滑相比并無明顯變化, 區(qū)域二標度律從0.36 上升至0.59 是由熱邊界層主導的, 區(qū)域三從0.30 上升至0.50 是由黏性邊界層主導的), 并給出Re和Ra的標度律關系從0.471 增加至0.551 是由大尺度環(huán)流的動力學發(fā)生變化引起的解釋.Stringano 等[12]研究了粗糙元尺寸h對系統(tǒng)傳熱的影響, 結果表明當系統(tǒng)的平均溫度邊界層厚度小于粗糙元尺寸時可以觀測到傳熱增強.上述研究結果都表明粗糙元能夠增強系統(tǒng)的傳熱.

粗糙元被認為可能有助于觀測到“終極區(qū)間”, 因為粗糙元的存在會破壞邊界層結構.Toppaladoddi 等[13]使用直接數值模擬計算帶有正弦形狀粗糙元的RB 系統(tǒng), 通過改變粗糙元波長, 發(fā)現存在一個最優(yōu)波長使得Nu數與Ra數接近滿足1/2 的標度律.根據“終極區(qū)間”理論, 只有在Ra數很大時才能觀察到1/2 標度律.而Zhu 等[14]采用相同的粗糙元將Ra數擴展到1012, 發(fā)現在Ra數相對較低的區(qū)間內能夠觀察到1/2標度律, 但是隨著Ra數增加,標度律指數會重新回到1/3,Nu數與Ra數1/2 標度律只能存在小范圍Ra數內, 以此證明Toppaladoddi 等[13]觀測到的并不是“終極區(qū)間”.最近, Zhu 等[15]采用多尺度粗糙元, 成功實現了在廣范圍Ra數內保持1/2 的標度律, Xia[16]則認為這是一種通過調整RB 系統(tǒng)邊界形狀來控制傳熱的新方法.

Zhang 等[17]通過改變三角形粗糙元高度h, 模擬二維方形對流槽, 結果表明: 對流槽加入粗糙元后并不是總意味著傳熱增加, 在粗糙元高度較小時反而會使傳熱減小;當粗糙元的高度h較小時, 相鄰粗糙元間腔體內部的流動主要受流體的黏性所控制, 這時熱流體被限制在這些腔體內部無法充分混合, 導致腔體內部溫度邊界層厚度增加, 從而阻礙了系統(tǒng)的整體熱輸運效率;相反, 當h足夠大時腔體內部逐漸增強的流體流動可以有效地混合冷熱流體, 導致更多更強的羽流結構生成, 使Nu數急劇增加.

上述工作主要側重于研究粗糙壁面對傳熱的影響, 而對Re數的關注度有所欠缺, 因此本工作主要在前人工作基礎上使用數值模擬的方法研究粗糙元高度對Re數的影響.

2 計算模型及參數

圖1 為RB 對流系統(tǒng)計算模型.對流槽寬度為D、高為H, 底板加熱溫度為無量綱溫度0.5, 頂板冷卻溫度為無量綱溫度?0.5;對流槽的側壁保持絕熱, 所有邊界均為無滑移邊界條件.對流槽上下板分布有等腰直角三角形形狀的粗糙元, 粗糙元的頂角保持直角不變, 通過改變粗糙元的數量來改變粗糙元的高度h, 且保持粗糙元與對流介質的接觸面積不變, 即達到改變邊界粗糙程度又不影響加熱接觸面積的效果.易知加入粗糙元后加熱面積是沒有粗糙元(光滑)情況下的1.41 倍.

圖1 RB 對流系統(tǒng)計算模型示意圖Fig.1 Schematic diagram of calculation model on RB

本工作采用交錯網格四階差分格式求解控制方程, 采用浸入邊界法處理粗糙元的邊界問題, 采用非等距網格跟蹤粗糙面附近的細節(jié), 并在求解過程中確保熱邊界層內有16 個以上的網格點, 以滿足網格分辨率[18].本工作計算所使用代碼在Zhang 等[17]、Bao 等[18]和Chen等[20]的工作中已得到了驗證, 計算過程主要在“天河二號”超算中心進行, 所有算例保持寬高比Γ=1 不變(D=1,H=1).Ra數的計算區(qū)間為107≤Ra≤109,Pr=0.7, 粗糙元數量由2個增加到70 個, 對應粗糙元高度h由0.125 減小至0.014.

3 結果分析和討論

Re數和系統(tǒng)內流體的運動相關, 可以表征系統(tǒng)的無量綱流動速度.本工作中Re數使用速度均方根Urms進行計算, 結合式(6), 有

式中:表示時間和空間上取平均.

圖2(a)為Ra=108,Pr=0.7 時Re數隨粗糙元高度h的變化規(guī)律, 可以發(fā)現: 隨著粗糙元高度h的增加,Re數首先有一個減小的過程, 達到極小值后Re數開始增加, 拐點位于粗糙元高度h=0.05 處, 對應的粗糙元數目為10;當Re數增加到極大值后, 隨著高度h繼續(xù)增加,Re數又有一個減小的過程, 此時拐點位于粗糙元高度h=0.1 處, 對應粗糙元數目為5.

光滑情況(粗糙元高度h= 0) 下的Re數記為Re(0), 對粗糙情況下Re數進行標準化, 結果如圖2(b)所示, 其中虛線表示Re(h)/Re(0)= 1, 即光滑情況Re數與粗糙情況Re數相等.可以發(fā)現: 粗糙元對Re數的影響很大, 隨著粗糙元高度不同, 對Re數既有增強作用也有抑制作用.只有當粗糙元高度h足夠大時, 粗糙元才能對Re數有增強作用, 但是當粗糙元高度h極大時,Re數又會受抑制.

粗糙元對Re數的作用非常復雜, 根據Re數隨粗糙元高度h變化的2 個拐點可以大致分為3 個區(qū)間(見圖2(b)).那么Re數隨粗糙元高度h如此變化的物理本質是什么呢？Re數由速度均方根Urms、對流槽高H和流體黏性系數ν共同決定, 算例中固定對流槽寬高比Γ和Pr數不變, 只改變粗糙元高度h, 那么一定是h影響Urms, 從而影響Re數.

圖2 Re 數隨粗糙元高度h 的分布(Ra=108, Pr =0.7)Fig.2 Re as a function of the roughness height h (Ra=108, Pr =0.7)

圖3 是速度時間平均場的分布, 其中箭頭代表速度的方向.區(qū)域一(見圖2(b))內Re數減小對應圖3(a)變化到圖3(b), 從圖中可以發(fā)現都存在一個明顯的大尺度環(huán)流, 但是在大尺度環(huán)流作用下, 圖3(a)中速度均方根Urms最大值為0.502, 而(b)中最大值為0.423.從Urms分布上來看, 圖3(a)中Urms值大的點(紅色區(qū)域)明顯多于圖3(b), 主要沿大尺度環(huán)流外側分布, 分別呈橢圓狀和呈圓形狀.隨著粗糙元高度h增加, 大尺度環(huán)流發(fā)現空間受限, 大尺度環(huán)流形狀逐漸由橢圓轉變?yōu)閳A形,Urms隨著減小導致Re數隨h的減小, 這也印證了Re數與大尺度環(huán)流的動力學密切相關[12].此外二次渦的作用也不可忽略, 二次渦是在大尺度環(huán)流剪切作用下粗糙元之間流體形成類似于“頂蓋流”[21]的流動.如圖4(a)所示, 粗糙元高度h較小時二次渦發(fā)展受限, 混合作用弱, 流動速度小[17], 這導致大量流體“滯留”在粗糙元之間, 對Re數貢獻非常小.

區(qū)域二(見圖2(b))內Re數增加對應圖3(b)變化到圖3(c), 在此區(qū)域內Re數隨著粗糙元高度h增加而增加, 此時圖2(c)中粗糙元高度h比圖3(b)中更大.隨著h的增加, 粗糙元對大尺度環(huán)流的抑制作用應該更強, 為何會如此反常呢.我們認為主要有兩個方面原因: 首先是羽流的作用, 羽流之間互相作用形成更大尺度的渦流進而最終形成大尺度環(huán)流[3];其次是二次渦的作用.Du 等[7-8]的研究結果表明粗糙元有利于羽流的生成, 隨著粗糙元的高度增加, 粗糙元對羽流生成作用越有利, 隨著羽流大量生成, 又為大尺度環(huán)流注入了動力, 這就導致了Urms增加從而Re數增加.圖4 所示的是粗糙元附近的Urms, 其中圖4(a)對應圖3(b)下底板正中間局部放大部分, 圖4(b) 對應圖3(c)下底板正中間局部放大部分.與區(qū)域一不同, 隨著粗糙元高度h增加, 粗糙元之間腔體越大, 更有利于二次渦發(fā)展, 二次渦越強[21].對比圖4(a)和(b)明顯可以發(fā)現, 由于二次渦變強圖, 4(b)中的Urms無論從數值大小還是分布都大于圖4(a).

圖3 速度時間平均場對應不同粗糙元高度h 的分布(Ra=108, Pr =0.7)Fig.3 Time aver velocity fields for differnet roughness height h (Ra=108, Pr =0.7)

圖4 粗糙元附近速度時間平均場(Ra=108, Pr =0.7)Fig.4 Time aver velocity fields near rough element (Ra=108, Pr =0.7)

由圖2(b)中的區(qū)域三可知,Re數隨著粗糙元高度h增加而減小.對比圖3(c)和(d)可以發(fā)現, 隨著粗糙元的進一步增加, 大尺度環(huán)流的發(fā)展空間變得很小, 這導致大尺度環(huán)流變弱, 此時Urms最大值僅為0.399.并且大尺度環(huán)流的影響范圍也變得十分有限, 這導致Urms在空間分布上不如圖3(c).

總體來說,Re數隨粗糙元高度h的分布主要是粗糙元對大尺度環(huán)流抑制作用以及與粗糙元對羽流的生成和二次渦的促進作用競爭的結果.當粗糙元對羽流的生成和二次渦的促進作用大于粗糙元存在對大尺度環(huán)流的抑制作用時,Re數增加, 反之則減小.

圖5(a)給出了不同Ra數時Re數隨著粗糙元高度h的變化規(guī)律.可以發(fā)現: 對于所有Ra數依舊能夠符合隨著粗糙元高度h增加Re數有一個先減小后增加再減小的變化規(guī)律;但是在不同Ra數下,Re數分布的拐點有所不同,Ra=107時拐點對應粗糙元高度最大,Ra=109時拐點對應粗糙元高度最小.這是因為隨著Ra數增加, 大尺度環(huán)流剪切能力越強, 越易于增加二次渦的強度, 也有利于羽流的產生, 所以當Ra數很大時, 粗糙元高度h很小就能對Re數增加起到促進作用.

圖5 不同Ra 數下Re 數隨粗糙元高度h 的分布(Pr =0.7)Fig.5 Re as a function of the roughness height h for different Ra (Pr =0.7)

對于光滑情況,Ra數與Re數大致滿足Re～Ra0.6的標度律關系.Zhang 等[22]的研究結果表明, 當Pr=0.7 時Re～Ra0.59, 本工作中Re～Ra0.57, 二者結果相近.粗糙元對標度律的影響如圖5(b)所示, 其中虛線為光滑時Ra數與Re數的標度律0.57.在粗糙元高度h= 0.05 時,粗糙元能提高Ra數與Re數的標度指數至0.62.固定h= 0.05,Ra= 109時粗糙元已經能促進Re數增加, 而Ra= 108是抑制Re增加,Ra= 107則抑制更多.這意味著粗糙元高度h= 0.05 時,Ra= 107的Re數是小于光滑情況的, 而Ra= 109的Re數是大于光滑情況的,直接導致了標度律指數相比于光滑情況有所增加.同理, 如果粗糙元對Ra= 109時Re數的抑制作用比Ra= 107時更大, 則表現為標度律指數減小, 如圖5(b)前半段.標度律指數變化的深層原因還是因為隨著Ra數增加, 二次渦在很小的粗糙元高度h就能對Re數起促進作用.

4 結 論

本工作通過數值模擬的方法研究二維粗糙邊界RB熱對流系統(tǒng)中粗糙元高度h對Re數的影響.

(1) 同一Ra數下, 隨著h的增加,Re數有一個先減小后增加再減小的過程, 這是粗糙元對大尺度環(huán)流的抑制作用與粗糙元對羽流的生成和二次渦的促進作用相互競爭所導致的.當粗糙元對羽流的生成和二次渦的促進作用大于粗糙元存在對大尺度環(huán)流的抑制作用時,Re數會增加.

(2) 不同Ra數下, 粗糙元對羽流的生成和二次渦的促進作用所需粗糙元高度h不同, 導致了粗糙情況下Ra數與Re數的標度律關系異常.