劉鎮(zhèn)濤 肖 莉 和 琨 汪 壘 ,2)

* (中國地質(zhì)大學(xué)(武漢)數(shù)學(xué)與物理學(xué)院,武漢 430074)

? (中國地質(zhì)大學(xué)(武漢)數(shù)學(xué)科學(xué)中心,武漢 430074)

引言

電流體動力學(xué)(electro-hydrodynamics,EHD)屬于電磁學(xué)和流體動力學(xué)等多學(xué)科交叉領(lǐng)域,它描述了帶電粒子在庫侖力驅(qū)動下引起的電對流[1].作為EHD 的基本課題之一,電熱對流(electro-thermo convection,ETC)考慮電場、流場和溫度場等多個物理場耦合及帶電粒子與流體介質(zhì)之間的相互作用,也是EHD 的重要研究課題之一[2-3].近年來,隨著能源問題的日益突出,節(jié)能環(huán)保等相關(guān)問題越來越受到廣泛的關(guān)注,科研人員對于電熱對流強(qiáng)化傳熱應(yīng)用逐漸重視[4-5].ETC 與許多其他傳統(tǒng)強(qiáng)化傳熱方法相比,具有節(jié)約資源、提高傳熱效率、結(jié)構(gòu)簡單和安全等優(yōu)點,在強(qiáng)化傳熱領(lǐng)域受到許多研究人員的關(guān)注[6].

近年來眾多學(xué)者對ETC 問題開展了相關(guān)研究,早期對ETC 的研究主要是對簡化的理論模型和實驗數(shù)據(jù)進(jìn)行分析.Atten[7]的研究結(jié)果表明電熱對流可以顯著提高傳熱效率,傳熱效率提升了15 倍.Mccluskey 和Atten[8]的研究結(jié)果顯示,電熱對流對傳熱的影響比純自然對流要高一個數(shù)量級.Wu 等[9]研究了電熱對流中亞臨界或超臨界的分叉類型和有限幅度穩(wěn)定性標(biāo)準(zhǔn),并發(fā)現(xiàn)了超臨界點和亞臨界點的取值很大程度上取決于普朗特數(shù)Pr和無量綱離子遷移率M.隨后,Traore 等[10]提出一種求解ETC 控制方程的數(shù)值模擬方案,該方案能夠很好解決ETC 問題的多物理場耦合問題,基于該模型,作者還對二維腔體中的ETC 進(jìn)行了數(shù)值研究.Dantchi 等[11]采用有限體積法求解ETC 方程,模擬中考慮了二維腔體內(nèi)從下壁面電極和從側(cè)壁面電極注入電荷兩種情況,結(jié)合相關(guān)無量綱參數(shù)進(jìn)行分析,研究結(jié)果表明電場增強(qiáng)了傳熱,對于足夠大的T值,努塞特數(shù)與瑞利數(shù)無關(guān).Luo 等[12]設(shè)計了格子Boltzmann 模型來求解ETC 方程,使用4 個不同的分布函數(shù)分別計算流場、電勢、電荷密度分布和溫度場,還進(jìn)行了多尺度分析來將離散的格子Boltzmann 方程恢復(fù)到宏觀控制方程,并通過對幾種問題進(jìn)行驗證,最后得到的結(jié)果與解析解等其他數(shù)值結(jié)果一致.除此之外,Lu 等[13]采用了高分辨率逆風(fēng)方案(high resolution upwind scheme)探究了方腔內(nèi)介電液體中電熱對流的流動結(jié)構(gòu)和傳熱現(xiàn)象,通過研究相關(guān)無量綱參數(shù)如熱瑞利數(shù)、無量綱離子遷移率以及普朗特數(shù),最后發(fā)現(xiàn)了多種流動結(jié)構(gòu)以及傳熱規(guī)律.Hu 等[14]采用了格子Boltzamnn 方法研究了不同幾何結(jié)構(gòu)下的二維電熱對流問題,探討了瑞利數(shù)、電瑞利數(shù)、長寬比和普朗特數(shù)等無量綱參數(shù)對流場流線、等溫線、電荷密度分布和界面的平均努塞特數(shù)的影響,結(jié)果表明增加電瑞利數(shù)可以提升傳熱效率并改變了流場的結(jié)構(gòu).Wang 等[15]提出了一種譜元法(spectral element method)的新型數(shù)值方法,在驗證該數(shù)值方法能有效的解決電熱對流問題后,又模擬了正弦溫度邊界方腔內(nèi)的二維電熱對流問題,研究結(jié)果表明,與均勻溫度邊界相比,正弦溫度邊界能提高傳熱效率.還有一些關(guān)于電熱對流的拓展研究.Chand 等[16]研究了兩平行電極板之間的納米流體在電場影響下的熱不穩(wěn)定性.Srivastava[17]對多孔介質(zhì)中的二元流體電熱對流進(jìn)行了穩(wěn)定性分析.

上述研究在很大程度上加深了人們對電熱對流在傳熱過程中的內(nèi)在機(jī)理.但是應(yīng)當(dāng)指出的是,上述研究,模型大多數(shù)是水平放置的兩平行電極板[7,9,12,18]或者是局部單極注入電荷[19-20],這些研究中較少考慮溫度分布不均勻的情況,實際上在許多研究領(lǐng)域如對自然對流的研究,非均勻溫度邊界問題一直是熱門的研究課題.Poulikakos[21]研究了二維方腔左側(cè)壁面上半部分加熱,下半部分冷卻,其他壁面隔熱的情況,結(jié)果表明,當(dāng)加熱區(qū)域位于冷卻區(qū)域下方時,整體的傳熱增強(qiáng).Sathiyamoorthy等[22]研究了二維方腔兩側(cè)壁面的溫度是具有線性變化的情況,結(jié)果表明,當(dāng)壁面兩側(cè)的溫度呈線性變化時,方腔底部的努塞特數(shù)會變得十分不穩(wěn)定.此外,Bilgen 和Yedder[23]研究了二維矩形腔體垂直壁面的溫度呈正弦分布而其他壁面絕熱的情況.另外,考慮加熱器長度和位置的情況,Oztop 和Abu-Nada[24]研究了部分加熱的二維腔體中納米流體的自然對流,發(fā)現(xiàn)了加熱器的位置會影響流線結(jié)構(gòu)和溫度分布.最近,Wang 等[25]研究了局部加熱的二維方腔中納米流體的自然對流,結(jié)果表明傳熱效率隨著加熱器的長度增加而降低.從這些文獻(xiàn)中不難看出,不論是在電熱對流中考慮局部單極注入電荷還是在自然對流中考慮非均勻溫度邊界情況下,都會對傳熱效率產(chǎn)生影響,并產(chǎn)生復(fù)雜的流動結(jié)構(gòu).而電熱對流本質(zhì)是電對流和自然對流耦合的問題,因此對于局部單極注入電荷和局部加熱相結(jié)合的機(jī)制需要進(jìn)行更深入的研究,以便于理解非均勻溫度邊界的電熱對流的流體流動和傳熱特性.

本文采用了格子Boltzmann 方法(LBM)求解控制方程組,其中包括Navier?Stokes 方程,電荷密度守恒方程,電勢的泊松方程以及溫度方程,另外,為了減少計算時間以及提升計算效率,本文使用NVIDIA圖形處理平臺,編寫計算統(tǒng)一設(shè)備架構(gòu)程序.研究了在局部受熱條件下,方腔內(nèi)流體流動和系統(tǒng)傳熱效率的變化情況,并為研究其他非均勻溫度邊界的電熱對流問題提供參考數(shù)據(jù).

1 問題描述

1.1 物理模型和控制方程

如圖1 所示,介電液體被限制在長度為H的二維方腔中.左側(cè)壁面是加熱的電極板,溫度值固定為θ0,電勢保持為 φ0,長度為h,電極板的中心點到下壁面的距離用 δ 來表示.右壁面溫度固定為 θ1,電勢保持為 φ1,且 θ0>θ1,φ0>φ1,上下壁面是絕熱且絕緣的.為方便模擬,假設(shè)方腔中的流體是不可壓縮的牛頓流體,流體的密度變化是滿足標(biāo)準(zhǔn)的Boussinesq假設(shè)的,自由電荷的產(chǎn)生只發(fā)生在介電液體與電極界面之間的電化學(xué)反應(yīng),并且注入的電荷密度取恒定值為q0.另外,忽略了黏性熱耗散以及焦耳熱的影響[26].

圖1 物理模型示意圖Fig.1 Schematic diagram of the physical model

基于上述假設(shè),二維局部受熱腔體內(nèi)電熱對流問題可以通過連續(xù)方程、動量方程、電勢的泊松方程、電場的定義方程、電荷密度守恒方程和溫度方程共同描述[27]

在上述方程組中,u,E,g,q,θ,θref,φ 分別為流體速度、電場強(qiáng)度、重力加速度、電荷密度、溫度、參考溫度、電勢.p,υ,β,ε ,K,D,χ 分別代表壓力、運(yùn)動學(xué)黏性、熱膨脹系數(shù)、介電常數(shù)、離子遷移率、電荷擴(kuò)散系數(shù)以及熱擴(kuò)散系數(shù).為了更方便的研究電熱對流問題,引入以下無量綱參數(shù)[28]

其中,Ra為熱瑞利數(shù),表示浮升力與黏性力的比值.T為電瑞利數(shù),表示庫侖力與黏性力的比值.對于熱瑞利數(shù)Ra以及電瑞利數(shù)T而言,在實際實驗中可依據(jù)溫度差和電壓差靈活設(shè)置[29],為了系統(tǒng)探究其影響,本文將之分別約束在 0 ≤T≤1200,1.0×104≤Ra≤1.0×105區(qū)間.C為電荷注入強(qiáng)度,羅康[30]的理論研究和數(shù)值模擬研究中,一般取C=10.0作為其強(qiáng)注入的代表性值,為此本文中所涉及的電荷注入強(qiáng)度C在未加強(qiáng)調(diào)的情況下,其值均為10.0.無量綱離子遷移率M是所謂的水動力學(xué)遷移率與真實遷移率的比值.對于帶有注入粒子的介電液體而言,實驗顯示其值一般大于等于3[3],例如甲醇(H+,M=4.1)、氯苯(C l?,M=4.9)、乙醇(H+,M=8.8;C l?,M=26.5)、硝基苯(C l?,M=22)、碳酸丙烯(C l?,M=51),鑒于此,為了更加清晰地認(rèn)識無量綱離子遷移率M對于電熱對流的影響,本文所考察的無量綱離子遷移率M介于5.0 和50.0 之間.α 為無量綱電荷擴(kuò)散系數(shù),其值一般介于 1 0?3和1 0?4之間[28],在本文中參考了文獻(xiàn)[31],將其值固定為 1 .0×10?4.Pr為動量邊界層與熱邊界層之比,反映了流體物理性質(zhì)對熱傳遞的影響,本文重點關(guān)注一類液體形態(tài)的介電介質(zhì)[3],因此將其值設(shè)置在5.0 和50.0 之間.另外,還考慮了電極板的長度h和電極板的中心點到下壁面的距離δ,分別取值為 0 .5H≤h≤H,0 .25H≤δ ≤0.75H.對于溫度邊界、電勢邊界、電荷密度邊界以及速度邊界,它們的數(shù)學(xué)描述如下,左側(cè)壁面電極板邊界條件如下

而對于速度邊界而言,方腔的4 個壁面均為無滑移邊界

2 數(shù)值方法

在數(shù)值方法選擇上,本文采用LBM 進(jìn)行模擬,與有限體積法和有限元方法等傳統(tǒng)方法不同,LBM 在處理復(fù)雜邊界問題時具有高效和簡單的優(yōu)點[32],又因為LBM 具有天然并行性[33],所以使得并行計算成為可能[34-35].經(jīng)過數(shù)十年的發(fā)展,LBM 的應(yīng)用領(lǐng)域從單相問題延伸到多相流[36]、多孔介質(zhì)[37-38]、在電流體動力學(xué)[9,12,30]等諸多領(lǐng)域,并受到了許多研究人員的廣泛關(guān)注.在本文中,擬并采用4 個演化方程分別求解Navier?Stokes 方程,電勢的泊松方程,電荷密度守恒方程和溫度方程.接下來將依次介紹.

2.1 流場格子Boltzmann 模型

在假設(shè)流體為不可壓牛頓流體的前提下,采用的是Guo 等[39]提出的格子Boltzmann 模型進(jìn)行求解,該模型的演化方程可以表示如下

其中外力F=qE+gβ(θ ?θref).τf為流場的無量綱松弛時間,它與 υ 的關(guān)系如下

最后計算宏觀速度u以及壓力項p,由以下公式計算得到

2.2 電勢場格子Boltzmann 模型

本質(zhì)上而言,電勢控制方程(3)為一類典型的泊松方程.在這里采用下述演化方程[41]

其中,ζ 為擴(kuò)散系數(shù),取 ζ 等于0.5 以保持?jǐn)?shù)值的穩(wěn)定性.此外,源項R=q/ε.平衡態(tài)分布函數(shù)定義如下

因為電勢的平衡態(tài)分布函數(shù)為線性的平衡態(tài)形式,為加快運(yùn)算效率,采用D2Q5 的格子模型[41],該模型的權(quán)系數(shù)取值如下

而電場E則可以通過式(4)和電勢的一階分布函數(shù)計算得到[12],具體計算表達(dá)式如下

2.3 電荷密度的格子Boltzmann 模型

從文獻(xiàn)[30]中可知電荷的傳輸機(jī)制有三種,第一種是在電場作用下,電荷會發(fā)生漂移;第二種是電荷在流體速度的影響下運(yùn)動;第三種則是電荷的擴(kuò)散.而電荷密度守恒方程屬于強(qiáng)對流方程,對于該方程的求解在算法中起著重要作用,因此為了能很好的求解該方程,使用羅康[30]提出的格子Boltzmann模型來求解,其演化方程表示如下

上 述等式中的權(quán)重 ωi,離散速度ci都與流場相同.

2.4 溫度場格子Boltzmann 模型

本文中的溫度方程是忽略了黏性熱耗散和焦耳熱后的簡化方程,從實現(xiàn)的角度和精度的滿足兩個方面,簡化的溫度方程能很好地應(yīng)用于各種傳熱問題中[42].因此采用文獻(xiàn)[12]中的方法,溫度的演化方程如下

3 邊界處理和程序流程

在格子Boltzmann 方法中,邊界條件的處理一般采用Guo 等[44]提出的非平衡態(tài)外推格式,所得的分布函數(shù)精度為二階,具有很好的數(shù)值穩(wěn)定性.另外,程序的具體步驟為:(1)初始化為系統(tǒng)平衡態(tài),此時電荷運(yùn)輸和溫度傳遞處于平衡狀態(tài),流場處于靜止?fàn)顟B(tài);(2)進(jìn)行碰撞流動,并通過式(19)計算速度u;(3)計算電勢的演化方程(21),結(jié)合式(25)和式(26)得到電勢 φ 和電場強(qiáng)度E;(4)計算電荷密度分布的演化方程(27),通過式(30)得到電荷密度q;(5)計算溫度的演化方程(31),并計算式(34) 得到溫度 θ ;(6) 計算式(17) 得到結(jié)果并代入流場的演化方程(13)中,依次重復(fù)步驟(2)～ 步驟(6)直到系統(tǒng)穩(wěn)定并收斂為止.最后,定義二維問題中的局部努塞特數(shù)Nuy和平均努塞特數(shù)Nuav[34]以及最大的流體速度Vmax分別用來表示傳熱效率和流動強(qiáng)度,通過如下公式計算得到

4 數(shù)值驗證

如圖2 所示,對系統(tǒng)處于平衡態(tài)時的電荷密度q和x方向的電場強(qiáng)度Ex的數(shù)值解與解析解進(jìn)行了比較.其中,平衡態(tài)電荷密度q和電場強(qiáng)度Ex的解析解表示如下

圖2 電荷密度 q,水平方向電場強(qiáng)度 E x 的數(shù)值解與解析解對比Fig.2 Comparisons of the charge density q and horizontal electric field strength E x between the analytical solutions and numerical solutions

上式中的參數(shù)a和b的取值與電荷注入強(qiáng)度C有關(guān)[12],在驗證中取C=10.

結(jié)果表明,通過當(dāng)前LBM 代碼計算得到的數(shù)值解 與解析解吻合較好.

5 數(shù)值模擬結(jié)果與討論

5.1 庫侖力和浮升力作用下的電熱對流問題探究

本節(jié)中,對二維局部受熱腔體結(jié)構(gòu)中的電熱對流問題進(jìn)行了模擬,并初步對該結(jié)構(gòu)下系統(tǒng)的傳熱效率、流體流動以及電熱對流的分岔結(jié)構(gòu)進(jìn)行了分析.實驗中將采用 3 00×300的網(wǎng)格分辨率,其可滿足網(wǎng)格無關(guān)性驗證.如圖3(a)所示,首先觀察了在熱對流很弱時(熱瑞利數(shù)Ra=400),電對流主導(dǎo)流體流動的情況,其他參數(shù)分別為C=10,M=10,Pr=10,h=0.5H,δ=0.5H.圖3(a)中描述了隨著電瑞利數(shù)T的增加和降低,相應(yīng)傳熱效率Nuav的變化情況,圖3(a) 中還給出了相關(guān)分岔點處的電荷密度分布圖.對于電瑞利數(shù)T的值逐漸增加的情況,當(dāng)電瑞利數(shù)T低于臨界值[45](Tc≈174,Pr=10,M=10) 時,Nuav的變化不是很大,在前半段幾乎保持水平狀態(tài),并且電荷密度分布比較均勻,這都表明流體沒有形成強(qiáng)烈的對流,庫侖力無法有效的激發(fā)流體流動,因此傳熱效率變化速度緩慢.然而當(dāng)電瑞利數(shù)T超過臨界值Tc以后,Nuav的值會突然向上跳躍式增加,并且電荷分布圖中會出現(xiàn)空電荷區(qū),空電荷區(qū)的形成與電荷運(yùn)輸過程中的漂移和對流機(jī)制之間的競爭有關(guān)[7,9,12,46].從電荷密度守恒方程(5)中的對流速度(KE+u)可以知道,一旦該區(qū)域的最大速度超過了電荷的漂移速度KE,強(qiáng)對流結(jié)構(gòu)就會阻礙擴(kuò)散到這些區(qū)域的電荷,結(jié)果就出現(xiàn)了空電荷區(qū),并且會發(fā)生傳熱增強(qiáng)的現(xiàn)象.對于逐漸減少電瑞利數(shù)T的情況,當(dāng)從電瑞利數(shù)T=300 開始逐漸減少時,Nuav并不能沿著增長時的曲線完全返回,這意味著強(qiáng)大的對流可以維持很大范圍內(nèi)的電瑞利數(shù)T值.但是當(dāng)電瑞利數(shù)T接近148 附近時,Nuav會突然下降,可以觀察到電瑞利數(shù)T=148 和T=174之間會形成所謂的回滯區(qū),并且如預(yù)期中那樣,空電荷區(qū)消失,Nuav會沿著曲線原路返回.接著觀察在電對流很弱時(電瑞利數(shù)T=80),熱對流主導(dǎo)流體流動的情況,圖3(b)中隨著Ra的增加和降低,往返的Nuav共用相同的演化曲線,表現(xiàn)出超臨界特征(Pr=10,M=10).另外,Nuav變化過程中并沒有出現(xiàn)回滯區(qū),這說明流體流動表現(xiàn)出規(guī)律且單一的機(jī)制.在圖3(b)中,還插入了熱瑞利數(shù)Ra=800時的電荷密度分布和溫度分布圖,通過觀察發(fā)現(xiàn)它們并沒有形成空電荷區(qū),電荷密度分布與溫度分布相似,說明了流體運(yùn)動十分弱,速度無法超過電荷的漂移速度.這也說明當(dāng)電瑞利數(shù)T的取值很低時,只憑借浮升力的作用無法產(chǎn)生電熱對流中分岔結(jié)構(gòu).通過對兩種不同情況下的電熱對流的穩(wěn)定性分析,能夠?qū)﹄姛釋α髦性S多現(xiàn)象的產(chǎn)生有一定的初步認(rèn)識,接下來會探討在庫侖力和浮升力共同作用下的電熱對流問題.

圖3 (a) 熱瑞利數(shù) Ra=400,變化電瑞利數(shù) T 的分岔結(jié)構(gòu)和(b) 電瑞利數(shù)T=80,變化熱瑞利數(shù)Ra 的分岔結(jié)構(gòu) (C=10,M=10,Pr=10,h=0.5H,δ=0.5H)Fig.3 (a) The distributions of the maximum speed Vmax vs.electric Rayleigh number T .(b) The distributions of the maximum speed Vmax vs.Rayleigh number R a (C=10,M=10,P r=10,h=0.5H,δ=0.5H)

固定參數(shù)Ra=1×104,C=10,M=10,Pr=10,h=0.5H,δ=0.5H,圖4 描繪了電瑞利數(shù)T在100～1200 范圍內(nèi),不同強(qiáng)度的庫侖力作用下的系統(tǒng)傳熱效率和流體流動的變化情況.可以看出,由于庫侖力和浮升力之間的相互競爭,每條曲線都對應(yīng)著不同的流動狀態(tài).隨著電瑞利數(shù)T的增加,庫侖力變大,Nuav和Vmax顯著增加,系統(tǒng)的傳熱效率和流動強(qiáng)度都會變大.值得注意的是,當(dāng)電瑞利數(shù)T較小(例如電瑞利數(shù)T=100)時,系統(tǒng)始終處于穩(wěn)定的對流狀態(tài),此時的傳熱效率和流動強(qiáng)度的增強(qiáng)是有限的.接著當(dāng)T增加時,則會在演化的最后出現(xiàn)一種不穩(wěn)定的對流狀態(tài)(T=200),并且Nuav會圍繞某一值波動.但是當(dāng)電瑞利數(shù)T的值增大時(5 00 ≤T≤800),系統(tǒng)再次出現(xiàn)穩(wěn)定的對流狀態(tài)(電瑞利數(shù)T=500)以及周期性流動狀態(tài)(電瑞利數(shù)T=600,T=800時),主要原因是庫侖力逐漸增大后,會在對流中占據(jù)主導(dǎo)地位.隨著電瑞利數(shù)T的增加(當(dāng)電瑞利數(shù)T=1200時),流動變得十分不穩(wěn)定,此時的庫侖力急劇增加并由穩(wěn)定或周期性的對流狀態(tài)變成混沌狀.

圖4 在不同電瑞利數(shù)T 下,(a) 平均努塞特數(shù) N uav 和(b) 最大速度 Vmax 隨時間變化的曲線 (R a=1×104,C=10,M=10,P r=10,h=0.5H,δ=0.5H)Fig.4 The distributions of (a) the average Nusselt number N uav vs.time and (b) the maximum speed Vmax vs.t ime under different electric Rayleigh number T (R a=1×104,C=10 ,M=10,P r=10,h=0.5H,δ=0.5H)

為了直觀地理解電瑞利數(shù)T對電熱對流的影響,如圖5 所示,給出了穩(wěn)定和非穩(wěn)定情況下的電荷密度分布、溫度分布和流線分布圖.圖5(a)～ 圖5(e)分別為電瑞利數(shù)T=0,2 00,5 00,8 00,1 200時,其余無量綱參數(shù)為Ra=1×104,C=10,M=1 0,Pr=10,h=0.5H,δ=0.5H.可以看到,當(dāng)電瑞利數(shù)T=0時,表示為純自然對流的情況,方腔內(nèi)只有一個沿順時針方向旋轉(zhuǎn)的漩渦,且上下對稱.當(dāng)電瑞利數(shù)T=200 時,流線分布圖與電瑞利數(shù)T=0時相似,都出現(xiàn)了一個漩渦,雖然在方腔頂部出現(xiàn)了兩個次級渦,但是庫侖力還比較弱,對流動和傳熱產(chǎn)生的貢獻(xiàn)很小,電荷只能在浮升力的影響下朝方腔頂部移動,下方的電荷密度則會變小.當(dāng)電瑞利數(shù)T=500時,流線變得十分復(fù)雜,不僅由一個漩渦分裂成兩個漩渦,而且方腔的4 個角落均出現(xiàn)了4 個小渦,此時的庫侖力急劇增加,對流動和傳熱效率的影響大于浮升力,因此引導(dǎo)電荷逐漸擺脫浮升力的束縛,進(jìn)而形成了上下對稱的空電荷區(qū),此時溫度的分布受到了庫侖力的影響,溫度的熱羽流與空電荷區(qū)形狀相似.當(dāng)電瑞利數(shù)T繼續(xù)增加到800 和1200 時,流線中的上下兩個漩渦都分別超過了分界線,流體流動的平衡再次被打破,由于強(qiáng)大的庫倫力導(dǎo)致了對流變得十分不穩(wěn)定,空電荷區(qū)變大,傳熱效率顯而易見的增加.

圖5 電瑞利數(shù)(a) T=0,(b) T=200,(c) T=500,(d) T=800,(e) T=1200 時的電荷密度分布(上)、溫度分布(中)和流線分布(下)圖(R a=1×104,C=10 ,M=10,P r=10,h=0.5H,δ=0.5H)Fig.5 The distributions of the charge density (top),temperature (middle) and streamlines (bottom) under different electrical Rayleigh number T.(a) T=0,(b) T=200,(c) T=500,(d) T=800,(e) T=1200 (R a=1×104,C=10,M=10,P r=10,h=0.5H,δ=0.5H)

最后,在圖6(a) 中,給出了在熱瑞利數(shù)Ra=1×104,5 ×104,1 ×105情況下的平均努塞特數(shù)Nuav隨電瑞利數(shù)T變化的曲線,并固定其他無量綱參數(shù)C=10,M=10,Pr=10,h=0.5H,δ=0.5H.通過上述分析,知道電瑞利數(shù)T是存在一個臨界閾值Tc,當(dāng)電瑞利數(shù)T超過該值時,Nuav會急劇增加.在圖6(a)中,可以發(fā)現(xiàn)類似的現(xiàn)象,熱瑞利數(shù)Ra=1×105時,此時電瑞利數(shù)T的臨界閾值在400～ 500 之間,電瑞利數(shù)T的取值超過500 以后,Nuav會發(fā)生明顯的變化.

圖6 (a) 探究不同電瑞利數(shù)T 下的平均努塞特數(shù) N uav 的變化過程和(b) 探究不同熱瑞利數(shù) R a 下的平均努塞特數(shù) N uav 的變化過程(C=10,M=10,P r=10,h=0.5H,δ=0.5H)Fig.6 (a) The distributions of the average Nusselt number N uav vs.electrical Rayleigh number T under different Rayleigh number R a .(b) The distributions of the average Nusselt number N uav vs.Rayleigh number R a under different electrical Rayleigh number T (C=10,M=10,P r=10,h=0.5H,δ=0.5H)

需要指出的是,熱瑞利數(shù)Ra=5×104時,電瑞利數(shù)T在 3 00～400之間出現(xiàn)的凹點是因為此區(qū)間內(nèi)的庫侖力和浮升力對腔體下半?yún)^(qū)域的對流影響相當(dāng),加之二者所致的渦流方向在此區(qū)域又相反[27],如此二者相互競爭尤為激烈,腔體內(nèi)的總體能量減少(主要體現(xiàn)在動能上),進(jìn)而導(dǎo)致傳熱效率相對電瑞利數(shù)T≤300時略有降低,從而出現(xiàn)凹點.而當(dāng)熱瑞利數(shù)Ra=104,電瑞利數(shù)T=300～400時,因此時浮升力相對庫侖力的影響較小,系統(tǒng)內(nèi)庫侖力影響絕對占優(yōu),因此方腔內(nèi)的傳熱效率隨著電瑞利數(shù)T的增加而增加,從而也不會出現(xiàn)上述凹點.當(dāng)電瑞利數(shù)T小于臨界閾值Tc時,浮升力的作用會增強(qiáng)系統(tǒng)的傳熱效率,但是當(dāng)電瑞利數(shù)T大于臨界閾值Tc以后,Nuav的上升是由于庫侖力增強(qiáng)導(dǎo)致的.當(dāng)電瑞利數(shù)T足夠大時,因為庫侖力在流動中發(fā)揮主導(dǎo)作用,反觀浮升力的影響可以忽略,所以傳熱增強(qiáng)不再取決于熱瑞利數(shù)Ra,可以看到3 條曲線最后收斂為一個點.在圖6(b)中,探究了不同電瑞利數(shù)T值情況下,Nuav隨著熱瑞利數(shù)Ra的演化曲線,當(dāng)電瑞利數(shù)T的取值大于500 以后,曲線的斜率是逐漸下降的,這是庫侖力在主導(dǎo)熱傳遞引起的,浮升力的影響在逐漸減弱,特別是當(dāng)電瑞利數(shù)T=800時,曲線幾乎是水平的,Nuav的變化與熱瑞利數(shù)Ra的取值無關(guān).可以注意到,與純自然對流問題(電瑞利數(shù)T=0)情況相比,Nuav顯著增加,當(dāng)熱瑞利數(shù)Ra=1×104時,在電瑞利數(shù)T的取值為500 的情況下,較于純自然對流傳熱大約增強(qiáng)了 8 4.31%,說明加入電場后會強(qiáng)化傳熱.在本節(jié)中,探究了電熱對流中的基本特征,并加深了對電熱對流中庫侖力和浮升力共同作用下流體流動結(jié)構(gòu)、溫度分布以及電荷的運(yùn)輸機(jī)制的理解.接下來,會對其 他無量綱參數(shù)進(jìn)行研究.

4. B短語辨析。A.為……擔(dān)心；B.為……感到自豪；C.害怕；D.為……感到遺憾。聯(lián)系前一句描述，可知此處指的是我為你感到自豪。故選B。

5.2 電極板位置影響探究

本節(jié)中,研究了電極板位置的影響.圖7 描繪了電極板位置對在不同電瑞利數(shù)T下的傳熱效率的影響,其余各無量綱參數(shù)為Ra=1×104,C=10,Pr=10,M=10,h=0.5H.其中,當(dāng)電瑞利數(shù)T=800,熱瑞利數(shù)Ra=1×104時,相比于浮升力,庫侖力影響絕對占優(yōu);將加熱壁面從下半?yún)^(qū)域 (δ=0.25H)更改為上半?yún)^(qū)域 (δ=0.75H)后,浮升力變化帶來的影響基本可以忽略,此時系統(tǒng)傳熱效率主要由庫侖力決定,又因兩種情況下流動結(jié)構(gòu)以及庫侖力分布等基本鏡面對稱,從而T=800時,上下兩種位置對應(yīng)的傳熱效率基本重合.另外,當(dāng)加熱壁面處于中間位置(δ=0.5H) 時,上下兩部分的速度分布相對更加均勻,流體流動速度比另外兩個位置下的要大,對流強(qiáng)度相對更強(qiáng),溫度邊界層也相對更薄,因而此位置下的傳熱效率一般優(yōu)于另外兩種情況.另外,在 0 ≤T≤225 區(qū)間,電極板在上半?yún)^(qū)域時(δ=0.75H)的Nuav小于在下半?yún)^(qū)域時(δ=0.25H)的Nuav.在圖7中畫出了電瑞利數(shù)T=200 時,A和B兩點處的平均努塞特數(shù)變化曲線圖,可以看出B點處的對流十分不穩(wěn)定,這也說明庫侖力和浮升力競爭相較于A點處更加激烈,并在極大程度上影響了系統(tǒng)的傳熱效率,因此 δ=0.75H時的傳熱效率低于 δ=0.25H時的傳熱效率,但是隨著電瑞利數(shù)T的取值繼續(xù)增加,庫侖力在對流中的影響增強(qiáng),當(dāng)電瑞利數(shù)T=800時,庫侖力占據(jù)絕對優(yōu)勢,此時系統(tǒng)的傳熱效率與熱瑞利數(shù)Ra無關(guān),所以兩個位置重合.最后,在電瑞利數(shù)T=200 到T=300之間,3 條曲線分別出現(xiàn)了一個臨界閾值,這與圖6(a)中熱瑞利數(shù)Ra=1×104時出現(xiàn)的臨界閾值范圍相同.

圖7 電極板在不同位置下,隨著電瑞利數(shù) T 增加,平均努塞特數(shù)Nuav 的變化曲線 (R a=1×104,C=10,P r=10,M=10,h=0.5H)Fig.7 The distributions of the average Nusselt number N uav vs.electrical Rayleigh number T under different heating locations δ (R a=1×104,C=10,P r=10,M=10,h=0.5H)

如圖8 所示,畫出了當(dāng)電瑞利數(shù)T=250時,δ=0.25H和 δ=0.75H時的電荷密度分布、溫度分布和流線分布圖.在這兩種情況下,電荷密度分布圖中都出現(xiàn)了空電荷區(qū),根據(jù)在5.1 節(jié)中的分析,可以知道,此時的庫侖力增大并與浮升力相互競爭,通過觀察和分析兩種情況下的流線分布(圖8(c1)和圖8(c2)),當(dāng)電極板在下半?yún)^(qū)域時,由于庫侖力增強(qiáng),導(dǎo)致了下部分對流減弱,并且圖8(c1)下部分的漩渦比圖8(c2)的漩渦小,因此系統(tǒng)的傳熱效率也降低了.所以從圖7 中可以發(fā)現(xiàn)電極板在上半?yún)^(qū)域的Nuav高于在下半?yún)^(qū)域時的Nuav,并且在之后的發(fā)展中一直保持(即使差異會越來越小).

圖8 電極板在上半?yún)^(qū)域和下半?yún)^(qū)域時的電荷密度分布圖(上)、溫度分布圖(中)和流線分布圖(下) (T=250,R a=1×104,C=10,Pr=10,M=10,h=0.5H)Fig.8 The distributions of the charge density (top),temperature(middle) and streamlines (bottom) as the electrode plate is at the top and bottom positions (T=250,R a=1×104,C=10 ,P r=10,M=10,h=0.5H)

圖9 給出了電瑞利數(shù)T=800時,電極板在下半?yún)^(qū)域圖9(a),中間圖9(b),上半?yún)^(qū)域圖9(c)位置時的電荷密度分布、溫度分布和流線分布圖.可以看到在電荷分布圖中都形成了一個很明顯的空電荷區(qū),此時的庫侖力在流動過程中占據(jù)主導(dǎo)地位.當(dāng)電極板在上半?yún)^(qū)域和下半?yún)^(qū)域時,庫侖力分布不均勻,從而導(dǎo)致速度大小分布不均勻.電極板在下半?yún)^(qū)域時,下半部分的速度大小會比上半部分的小,因此方腔下半部分的熱量轉(zhuǎn)移相比上半部分的熱量轉(zhuǎn)移要少,從溫度分布圖中的熱羽的形狀也可以發(fā)現(xiàn)靠近壁面的溫度邊界層變厚.無論是電極板在上半?yún)^(qū)域還是在下半?yún)^(qū)域,與電極板在中間時相比,都會使得靠近壁面一側(cè)的溫度梯度降低,從而傳熱效率變低.通過以上分析說明了當(dāng)電極板放在中間位置時,傳熱效率最佳.另外,電極板在上半?yún)^(qū)域和下半?yún)^(qū)域時的漩渦結(jié)構(gòu)類似,這也說明了為什么傳熱效率會隨著電瑞利數(shù)T的取值增加而最終接近.通過對比圖8 和圖9 中的電荷密度分布,可以發(fā)現(xiàn)電荷分布會受到庫侖力和浮升力的影響,當(dāng)電瑞利數(shù)T的取值很小時,庫侖力較弱,左壁面的電荷會由于浮升力作用向上壁面運(yùn)動,而在右壁面一定范圍內(nèi)的電荷向下運(yùn)動后,接著會沿著下壁面向左壁面流動,從而形成一個順時針的對流.雖然庫侖力始終向右并會對下部分一定區(qū)域的對流造成阻礙,但是電瑞利數(shù)T的取值很小,浮升力大于庫侖力,因此流動的方向仍然和在浮升力影響下的方向一致.但是隨著電瑞利數(shù)T的取值增加,庫侖力逐漸增加,流動的方向逐漸受到庫侖力的影響,最終的對流方向會被庫侖力主導(dǎo),這也可以從圖9 中的溫度分布圖中的熱羽形狀和電荷分布形狀相似看出.

圖9 電極板在(a)下半?yún)^(qū)域,(b)中間,(c)上半?yún)^(qū)域,電荷密度分布(左)、溫度分布(中)和流線分布(右)圖 (T=800,R a=1×104,C=10,Pr=10,M=10,h=0.5H)Fig.9 The distributions of the charge density (left),temperature (middle) and streamlines (right) as the electrode plate is at the top,middle and bottom positions (T=800,R a=1×104,C=10,P r=10,M=10,h=0.5H)

5.3 電極板長度影響探究

本節(jié)中,考慮了電極板長度h對傳熱的影響,相應(yīng)的數(shù)據(jù)結(jié)果如圖10 所示,其余各無量綱參數(shù)為Ra=1×104,C=10,Pr=10,M=10,δ=0.5H.從圖10中可以看出,對于給定的電瑞利數(shù)T,傳熱效率會隨著電極板的長度增加而降低,并且3 種情況下電瑞利數(shù)T的臨界閾值都在 2 00～300之間,在大于臨界閾值以后,它們的變化趨勢十分相似,再結(jié)合上一節(jié)中的臨界閾值,發(fā)現(xiàn)電極板長度h和位置 δ 對臨界閾值的范圍影響不大.

圖10 電極板長度不同時,隨著電瑞利數(shù) T 增加,平均努塞特數(shù)Nuav 的變化曲線 (R a=1×104,C=10,P r=10,M=10,δ=0.5H)Fig.10 The distributions of the average Nusselt number N uav vs.electrical Rayleigh number T under different length of electrode plate h(R a=1×104,C=10 ,P r=10,M=10,δ=0.5H)

分別畫出了電瑞利數(shù)T=250(圖11)和電瑞利數(shù)T=300(圖12)時的電荷密度分布、溫度分布和流線分布圖.在圖11 中,隨著電極板的長度增加,溫度分布圖中的熱羽流變小至一半并類似于純自然對流(電瑞利數(shù)T=0),空電荷區(qū)向下移動并逐漸遠(yuǎn)離電極板,從流線分布圖中可以看出下方的一級漩渦最終退化為次級漩渦,下方的對流由有序逐漸變?yōu)闊o序,從而整體的傳熱效率逐漸下降.又如圖12 所示,當(dāng)電瑞利數(shù)T=300時,此時的庫侖力在與浮升力的競爭中占據(jù)了優(yōu)勢,并且在方腔的中心區(qū)域形成了空電荷區(qū).但是隨著電極板兩端的增加,電場力減弱,熱壁面的速度下降,無法將更多的熱量轉(zhuǎn)移到介電液體中,同時熱壁面的溫度邊界層變厚,溫度梯度降低,因此傳熱效率會隨著電極板的長度增加而降低.

圖11 不同長度下,(a) h=0.5H,(b) h=0.75H,(c) h=H,電荷密度分布圖(左)、溫度分布圖(中)和流線分布圖(右) (T=250,R a=1×104,C=10 ,P r=10,M=10,δ=0.5H)Fig.11 The distributions of the charge density (left),temperature (middle) and streamlines (right) under different length of electrode plate.(a) h=0.5H,(b) h=0.75H,(c) h=H (T=250,R a=1×104,C=10,P r=10,M=10,δ=0.5H)

圖12 不同長度下,(a) h=0.5H,(b) h=0.75H,(c) h=H,電荷密度分布圖(左)、溫度分布圖(中)和流線分布圖(右) (T=300,R a=1×104,C=10 ,P r=10,M=10,δ=0.5H)Fig.12 The distributions of the charge density (left),temperature (middle) and streamlines (right) under different length of electrode plate.(a) h=0.5H,(b) h=0.75H,(c) h=H (T=300,R a =1×104,C=10,P r=10,M=10,δ=0.5H)

5.4 無量綱離子遷移率和普朗特數(shù)的聯(lián)合影響

如圖13 所示,還研究了無量綱離子遷移率M和普朗特數(shù)Pr的聯(lián)合影響,其余無量綱參數(shù)為Ra=1×104,T=300,C=10,δ=0.5H,h=0.5H.從圖13 可以看出在給定Pr的情況下,隨著M的增加,平均努塞特數(shù)Nuav會從一開始的急劇下降轉(zhuǎn)變?yōu)橼厔萜骄彽那闆r.M越小,Nuav會隨著Pr增加而增加,當(dāng)M很大時(M=50),Nuav幾乎保持不變.事實上,如文獻(xiàn)[9]指出,M的大小決定了靜電方程組與Navier-Stokes 方程組的耦合緊密程度,M越大,兩個方程組的耦合程度越弱,進(jìn)而庫侖力對流場的作

用減弱.另外,不同的M和Pr數(shù)的流體流動狀態(tài)也存在一定差異,主要分為穩(wěn)態(tài)和非穩(wěn)態(tài),具體可見圖13中的標(biāo)記.

圖13 固定普朗特數(shù)Pr,在不同無量綱離子遷移率M 下,平均努塞特數(shù) N uav 的變化曲線 (R a=1×104,T=300,C=10,δ=0.5H,h=0.5H)Fig.13 The distributions of the average Nusselt number N uav vs.dimensionless ion mobility M under different Prandtl number P r (R a=1×104,T=300,C=10,δ=0.5H,h=0.5H)

最后,取圖13 中Pr=50的曲線肘部處和端點處的值即M=25,M=30和M=50,畫出這3 點處的電荷密度分布、溫度分布和流線分布圖,如圖14所示,當(dāng)M=25和M=30時,空電荷區(qū)會變得不明顯,還會偏移中心向下且遠(yuǎn)離電極板,這都說明庫侖力的影響在減弱,而且流線分布圖中的漩渦結(jié)構(gòu)的差異更加明顯,由兩個反向旋轉(zhuǎn)的漩渦變成了一個漩渦結(jié)構(gòu),這不僅減弱了對流而且破壞了有序的對流結(jié)構(gòu),因此降低了系統(tǒng)的傳熱效率.當(dāng)M=50時,差異會變得十分明顯,雖然離子遷移率增加導(dǎo)致大量電荷注入到方腔中,但是兩方程組的耦合程度減弱,庫侖力減弱導(dǎo)致系統(tǒng)的對流減弱.所以當(dāng)M越大時,系統(tǒng)的傳熱效率反而會減小.

圖14 不同M 值下的電荷密度分布(左)、溫度分布(中)和流線分布(右)圖.(a) M=25,(b) M=30,(c) M=50 (P r=50,R a=1×104,T=300,C=10,δ=0.5H,h=0.5H)Fig.14 The distributions of the charge density (left),temperature (middle) and streamlines (right) under different dimensionless ion mobility M .(a) M=25,(b) M=30,(c) M=50 (P r=50,R a=1×104,T=300,C=10,δ=0.5H,h=0.5H)

6 結(jié)論

本文采用了格子Boltzmann 方法對二維局部受熱腔體結(jié)構(gòu)中電熱對流問題進(jìn)行模擬,首先對系統(tǒng)的傳熱效率、流體流動以及電熱對流的分岔結(jié)構(gòu)進(jìn)行了分析,對在庫侖力和浮升力共同作用下流體的流動、溫度的分布以及電荷的運(yùn)輸機(jī)制有了初步理解,然后探究了電極板的位置和長度的影響,最后探討了無量綱離子遷移率M和普朗特數(shù)Pr的聯(lián)合影響.根據(jù)以上研究中呈現(xiàn)的數(shù)值結(jié)果,主要結(jié)論如下:

(1) 固定熱瑞利數(shù)Ra,平均努塞特數(shù)Nuav會隨著電瑞利數(shù)T的增加而增加,這主要是由于庫侖力增強(qiáng)導(dǎo)致的.當(dāng)電瑞利數(shù)T超過臨界閾值(Tc,Pr=10,M=10)后,Nuav會急劇增加,并且熱瑞利數(shù)Ra越大,臨界閾值也會增加.通過穩(wěn)定性分析可以發(fā)現(xiàn),電瑞利數(shù)T會存在一個亞臨界點Tf,然而熱瑞利數(shù)Ra不會出現(xiàn)分岔結(jié)構(gòu),會表現(xiàn)為超臨界特征.另外,當(dāng)電瑞利數(shù)T足夠大時,庫侖力占據(jù)絕對優(yōu)勢,浮升力對傳熱的影響減弱,Nuav的值不再依賴于熱瑞利數(shù)Ra.

(2) 傳熱效率很大程度上受到電極板位置的影響,并且在中間位置時的傳熱效率最佳,這主要是由于電極板在中間位置時的對流最強(qiáng)導(dǎo)致的.

(3) 傳熱效率會隨著電極板的長度增加而降低,這主要是由于電極板長度增加以后,方腔底部的對流減弱,從而導(dǎo)致傳熱效率降低.

(4) 當(dāng)M越大時,系統(tǒng)的傳熱效率越低,這是因為隨著M的增加,庫侖力會對流場的影響減弱,從而導(dǎo)致對流減弱和傳熱效率降低.