陳捷超，李友榮，于佳佳

(重慶大學(xué)動力工程學(xué)院低品位能源利用技術(shù)及系統(tǒng)教育部重點實驗室，重慶 400044)

引 言

在混合溶液中，當自由表面上存在溫度和濃度不均勻時，會引起表面張力的不均勻，從而誘發(fā)耦合熱-溶質(zhì)毛細對流；同樣，當液層內(nèi)存在溫度和濃度不均勻時，會引起密度的不均勻，在重力場中會誘發(fā)浮力對流，或稱雙擴散對流。這兩種對流常常耦合共存于晶體生長、合金凝固、混合工質(zhì)相變傳熱等過程中，稱之為雙擴散毛細對流[1]。

迄今為止，對雙組分溶液的耦合熱-溶質(zhì)毛細對流和雙擴散浮力對流已進行了比較多的研究，取得了豐碩的成果。Ghorayeb等[2-3]和Xin等[4]通過分岔分析，討論了浮力比Rρ=-1時超臨界中心對稱的三渦胞流動結(jié)構(gòu)和次臨界順時針單胞流動結(jié)構(gòu)，分析了熱邊界層和溶質(zhì)邊界層的厚度與工質(zhì) Lewis數(shù)（Le）間的關(guān)系。Bardan等[5]和 Nishimura等[6]采用數(shù)值模擬、非線性穩(wěn)定性分析和能量分析等方法，闡述了不同浮力比下流動結(jié)構(gòu)的演變過程和流動失穩(wěn)的機理。Bergeon等[7]對大深寬比矩形腔內(nèi)的雙擴散對流進行了三維數(shù)值模擬和穩(wěn)定性分析，得到了豐富的流場結(jié)構(gòu)，詳細闡述了全局分岔研究中得到的混沌振蕩流動特征。Sezai等[8]對立方腔體內(nèi)的雙擴散對流進行了數(shù)值模擬，分析了浮力比在-2

對具有自由表面的液池在水平方向上施以大小相等、方向相反的熱毛細效應(yīng)和溶質(zhì)毛細效應(yīng)時，表面張力比Rσ=-1，此類邊界條件稱為滯止邊界條件[9]。Li等[10]通過對環(huán)形液池內(nèi)耦合熱-溶質(zhì)毛細對流的二維數(shù)值模擬，揭示了穩(wěn)態(tài)流動向周期振蕩流動轉(zhuǎn)變的物理機理。Chen等[9,11]對矩形液池內(nèi)耦合熱-溶質(zhì)毛細對流進行了三維數(shù)值模擬和線性穩(wěn)定性分析，討論了流動失穩(wěn)轉(zhuǎn)變的 Hopf分岔特性，揭示了周期性振蕩流動發(fā)生的物理機制。同時，Li等[12]分析了流動從周期性振蕩流動到混沌狀態(tài)的演變序列。

目前，對于環(huán)形液池內(nèi)純工質(zhì)的熱毛細對流和熱毛細-浮力對流已進行了比較多的研究[13-17]，但很少涉及到混合溶液的毛細對流和浮力對流的耦合。本文報道了一組表面張力比Rσ=-1時環(huán)形液池內(nèi)二元混合溶液雙擴散毛細對流的三維數(shù)值模擬結(jié)果，主要目的是了解環(huán)形液池內(nèi)雙擴散毛細對流的流動特征及失穩(wěn)機理，揭示液池深寬比對流型及其轉(zhuǎn)變過程的影響。

1 物理數(shù)學(xué)模型

物理模型如圖1所示，環(huán)形液池內(nèi)徑為ri，外徑為ro，深度為d，深寬比為A=d/(ro-ri)。液池底部為絕熱固壁，頂部為不變形絕熱自由界面，內(nèi)、外壁面分別維持恒定溫度Ti、To（Ti

圖1 物理模型Fig.1 Physical model

為簡化起見，假定：（1）流體為不可壓縮牛頓流體，流動為層流；（2）除表面張力和浮力項中的密度外，所有物性參數(shù)都為常數(shù)；（3）自由表面考慮熱毛細力和溶質(zhì)毛細力作用，其他固壁無滲透且滿足無滑移邊界條件；（4）浮力項中密度和表面張力都是溫度和濃度的線性函數(shù)，即

取量綱1參考長度、時間、速度和壓力分別為(ro-ri)、(ro-ri)2/ν、ν/(ro-ri)和ρν2/(ro-ri)2，則量綱 1控制方程為

邊界條件如下：

式中，Θ=(T-Ti)/(To-Ti)和Φ=(C-Ci)/(Co-Ci)分別為量綱1溫度和濃度，V為量綱1速度矢量，P為量綱1壓力，τ為量綱1時間，R和Z分別為量綱1坐標。

當液池自由表面上存在溫度或濃度的不均勻時，可用熱或溶質(zhì)毛細Reynolds數(shù)（Re）來表征熱或溶質(zhì)毛細力的大小；同樣地，當液池內(nèi)部存在溫度或濃度的不均勻分布時，可用熱或溶質(zhì) Grashof數(shù)（Gr）來表征熱或溶質(zhì)浮力的大小。熱毛細ReT和熱GrT分別定義為

式中，ν為動量擴散率。

溶質(zhì)毛細力和溶質(zhì)浮力的大小可由表面張力比Rσ和浮力比Rρ確定

液池內(nèi)外半徑比固定在ri/ro=0.5，工質(zhì)是質(zhì)量分數(shù)為0.2627/0.7373（對應(yīng)的摩爾分數(shù)為0.25/0.75）的甲苯/正己烷混合溶液，Prandtl數(shù)為Pr=ν/α=5.54，Lewis數(shù)為Le=α/D=25.78，這里，α為熱擴散率。其他物性參數(shù)如表1所示。

表1 甲苯/正己烷混合溶液物性參數(shù)Table 1 Properties of tolunen/n-hexane solution

控制方程采取有限容積法離散，對流項采用二階迎風(fēng)，擴散項采用二階中心差分，壓力-速度修正采用 SIMPLE方法，量綱 1時間步長為(0.12～1.1)×10-4，迭代收斂條件任意時間步長內(nèi)溫度、濃度和速度的最大相對誤差均小于10-4。

經(jīng)過反復(fù)驗證和比較，考慮到計算精度和計算時間，最終選定的計算網(wǎng)格為 100(R)×120(θ)×20(Z)。為了檢驗計算方法的正確性，首先對深寬比為0.5的矩形池內(nèi)雙擴散毛細對流進行了二維數(shù)值模擬，計算條件為Rσ=-1、Pr=5和Le=5，結(jié)果表明，計算獲得的等流函數(shù)線和等溫線分布與文獻[9]一致，監(jiān)測點處的振蕩頻率與文獻[9]給出值的最大偏差為2.4%。然后，在與Zhan等[11]相同條件下對矩形池內(nèi)雙擴散 Marangoni對流進行了計算，當Pr=5和Le=10時，計算得到的豎直壁面上平均Nusselt數(shù)（Nu）與文獻[11]給出的 Nusselt數(shù)最大偏差為2.8%。上述兩個算例表明，本文的計算方法是合理的，計算結(jié)果是可信的。

2 計算結(jié)果與討論

2.1 穩(wěn)態(tài)流動

很多學(xué)者對零重力下表面張力比Rσ=-1時的毛細對流進行了研究，發(fā)現(xiàn)在溫差較小的液池內(nèi)，存在流動速度為零的純導(dǎo)熱狀態(tài)，并稱之為靜止平衡狀態(tài)[9,11]，只有當溫差足夠大時，液池內(nèi)才會出現(xiàn)流動。然而，在常重力條件下，如果熱浮力和溶質(zhì)浮力不平衡，則浮力的存在使得液池中只要有溫差，就會出現(xiàn)流動。由于液池中熱浮力與溶質(zhì)浮力同時存在，并且浮力比Rρ= -2.18≠-1，可見系統(tǒng)受到的體積力不平衡，溶質(zhì)浮力占主導(dǎo)地位，因此，液池中存在溫差（濃度差）時就會產(chǎn)生浮力對流。圖 2給出了ReT=10時R-Z截面上等流線、等溫線和等濃度線分布，此流動為二維軸對稱穩(wěn)態(tài)流動，液層上部近表面的流體從內(nèi)壁流向外壁后，從底部附近回流，在R-Z截面上形成一個順時針方向旋轉(zhuǎn)的渦胞。由于熱擴散速度大于溶質(zhì)擴散速度（Le>1），所以等溫線分布均勻，等濃度線分布受流體對流的影響較大。

圖2 A=0.15和ReT=10時R-Z截面等流函數(shù)線（上）、等溫線（中）和等濃度線（下）分布Fig.2 Streamlines (uppers), isotherms (middle) and isoconcentrations (bottom) at R-Z section at A=0.15 and ReT=10

圖3 監(jiān)測點P處徑向速度振幅與頻率的變化（符號☆為臨界點）Fig.3 Variations of amplitude and frequency of radial velocity at monitor point P (☆ denotes critical value)

2.2 臨界熱毛細Reynolds數(shù)

隨著水平溫差（濃度差）的增大，液池內(nèi)的流動會由二維軸對稱穩(wěn)態(tài)流動經(jīng) Hopf分岔轉(zhuǎn)變?yōu)槎S軸對稱周期性振蕩流動，轉(zhuǎn)變點對應(yīng)的熱毛細Reynolds數(shù)稱為臨界熱毛細Reynolds數(shù)ReTc。圖3給出了A=0.15時監(jiān)測點P（R=1.5,θ=0,Z=0.145）處徑向速度振幅|VR,M|和量綱 1頻率F隨熱毛細Reynolds數(shù)ReT的變化，采用線性外推法可得臨界熱毛細 Reynolds數(shù)為ReTc=231.6，臨界頻率為Fc=10.6。

2.3 周期性振蕩流動

圖4(a)給出了ReT=240和400時監(jiān)測點P處的徑向速度VR隨時間τ的變化，顯然，ReT=240時徑向速度會發(fā)生周期性振蕩，但總是大于零，即流體從內(nèi)壁流向外壁；而ReT=400時監(jiān)測點P處徑向速度隨時間周期性地改變流動方向，并且其振幅遠大于ReT=240時的振幅。ReT=400時內(nèi)部流場結(jié)構(gòu)如圖5所示，液層內(nèi)部出現(xiàn)了8個逆時針和順時針相間旋轉(zhuǎn)的小渦胞，這些渦胞在內(nèi)壁附近產(chǎn)生，沿徑向方向向外傳播，最后與外壁附近的一直保持順時針方向旋轉(zhuǎn)的渦胞合并。由于熱擴散速度較快，等溫線分布受對流影響較小，等濃度線則因質(zhì)量擴散速度較慢而在對流作用下變得傾斜，并且內(nèi)壁底部附近濃度梯度較大。流動失穩(wěn)的物理機制可這樣解釋，如果將圖2所示流動結(jié)構(gòu)視為基礎(chǔ)流場[14]，顯然，基礎(chǔ)流場中自由表面流體流速比液層內(nèi)部大，若液層內(nèi)部低速流體在某一擾動驅(qū)動下流至自由表面高速流動區(qū)，就會產(chǎn)生慣性滯后[18-19]，從而誘導(dǎo)順時針和逆時針兩個方向旋轉(zhuǎn)的小渦胞，同時，小渦胞在基礎(chǔ)流場上部流動驅(qū)動下往外壁傳播。

從圖4(b)可看出，ReT=400時監(jiān)測點P處切向速度幾乎為零，自由表面等溫線和等濃度線呈同心圓分布，如圖 6(a)所示，因此，液池內(nèi)流動為二維軸對稱周期性振蕩流動。

為了揭示振蕩流動時溫度和濃度的波動規(guī)律，定義波動值（δX）為

圖4 監(jiān)測點P處速度隨時間τ變化Fig.4 Variations of velocities with time τ at monitor point P

圖5 ReT =400時R-Z截面一個周期（τp）內(nèi)流線（上）、等溫線（中）和等濃度線（下）分布Fig.5 Streamlines (uppers), isotherms (middle) and isoconcentrations (bottom) at R-Z section in every quarter-period at ReT =400

圖6 某一時刻自由表面等溫線（實線）和等濃度線（虛線）分布Fig.6 Snapshots of distributions of surface isotherms (solid lines) and isoconcentrations (dotted lines)

式中，X為溫度或濃度。

圖7為ReT=800時自由表面的溫度和濃度波動圖，波動條紋如向日葵的花瓣呈相間分布，由中心向外分成4層，每一層有7個由里向外逐漸增大的“花瓣”。液層內(nèi)部的流動結(jié)構(gòu)與圖5所示的類似，區(qū)別在于ReT=800時R-Z截面上流胞個數(shù)變?yōu)?，并且流胞尺度比ReT=400時的大，流動也更強。由圖4(b)可知，ReT=800時監(jiān)測點P處的切向速度Vθ隨時間呈周期性振蕩，液池內(nèi)的流動為三維周期振蕩流動。由于溫度、濃度與速度相耦合，因而波動圖中的“花瓣”在內(nèi)壁附近產(chǎn)生，與液層內(nèi)部的渦胞運動方向相同，在徑向方向上向外壁傳播，同時，以量綱1角速度0.48沿周向順時針旋轉(zhuǎn)。由于熱擴散速度比質(zhì)擴散快，溫度擾動迅速地被耗散，因此，圖7中量綱1溫度擾動幅度只有量綱1濃度擾動幅度的一半左右。此外，ReT=800時液層內(nèi)部的溫度和濃度分布與ReT=400時（如圖5所示）的類似，等濃度線更容易受對流的影響；同理，圖6(b)中自由表面等溫線（實線）比等濃度線（虛線）分布更均勻，等濃度線的彎曲更明顯。

2.4 深寬比的影響

根據(jù)文獻[2,7]的定義，隨著溫差（濃度差）的增大，液池內(nèi)流體的流動狀態(tài)將由穩(wěn)態(tài)向非穩(wěn)態(tài)轉(zhuǎn)變，轉(zhuǎn)變過程中若溫差不大，此時系統(tǒng)的能量不足以維持一個有序的流動結(jié)構(gòu)，流動狀態(tài)則為振幅有限但頻率無規(guī)律的振蕩，此狀態(tài)稱為亞臨界非穩(wěn)態(tài)。圖8(a)～(c)所示為A=0.5和ReT=400時的亞臨界非穩(wěn)態(tài)流場結(jié)構(gòu)，靠近外壁面上部的流體豎直往下流動，至液池底部處開始流向液池中心及內(nèi)壁底部，然后沿內(nèi)壁面往上流，同時在對角線方向上回流至外壁面上部，最后形成順時針方向旋轉(zhuǎn)的扁長狀流胞。液池底部附近的流體在上述扁長狀流胞的誘導(dǎo)下形成一個逆時針方向旋轉(zhuǎn)的小渦胞；液層上部的流域有兩個逆時針方向旋轉(zhuǎn)的流胞，靠近內(nèi)壁的流胞比另一個靠近外壁的流胞稍大。等濃度線分布受對流影響較大，在外壁上方及內(nèi)壁底部形成了較大的局部溶質(zhì)梯度，等溫線因熱擴散較快而受對流影響相對較小，所以等溫線比等濃度線分布均勻；同時，從圖中還可看出，由于液層上部兩個逆時針方向旋轉(zhuǎn)渦胞較強，因此，上半?yún)^(qū)域等溫線的彎曲程度比下部區(qū)域大。

圖 9(a)、(b)分別給出了A=0.5和ReT=400時的自由表面溫度和濃度波動。顯然，A=0.5和ReT=400時的流動狀態(tài)為亞臨界非穩(wěn)態(tài)[2,7]，自由表面波動表現(xiàn)為 12片大小不一的沿圓周方向分布的“花瓣”，外部為一條“圓環(huán)帶”，分別對應(yīng)著自由表面下兩個逆時針方向旋轉(zhuǎn)渦胞，從波動圖顏色深淺程度可看出，“花瓣”的擾動幅值比“圓環(huán)帶”的大。

當A=1和ReT=400時，液池內(nèi)的流動為三維穩(wěn)態(tài)流動，如圖8(d)～(f)所示，兩個逆時針方向旋轉(zhuǎn)的大流胞分別占據(jù)液層上部和下部流域，兩者間交界區(qū)域出現(xiàn)了一個沿對角線方向傾斜的扁長狀流胞；流場中等溫線分布受流動影響呈“Z”形，內(nèi)壁下半部及外壁上半部附近的等濃度線較密，出現(xiàn)厚度約為A/20的濃度邊界層。圖9(c)、(d)所示為自由表面上的波動圖，顯然，其為典型徑向波，波數(shù)為12。

圖7 ReT=800時一個周期（τp）內(nèi)自由表面溫度波動[(a)～(d)]和濃度波動[(e)～(h)]Fig.7 Snapshots of surface temperature fluctuation [(a)—(d)] and concentration fluctuation [(e)—(h)] in every quarter-period at ReT =800

圖8 ReT =400時R-Z截面流線（左）、等溫線（中）和等濃度線（右）分布Fig.8 Streamlines (left), isotherms (middle) and isoconcentrations (right) at R-Z section at ReT =400

圖9 ReT =400時自由表面溫度和濃度波動Fig.9 Snapshots of surface temperature fluctuation and surface concentration fluctuation

在立方體液池內(nèi)，當浮力比為-2時穩(wěn)態(tài)雙擴散自然對流中沿對角線方向傾斜的流胞占據(jù)液池大部分流域，在其誘導(dǎo)下液層上部及下部出現(xiàn)了兩個小流胞[8]。相比之下，在具有自由表面的環(huán)形液池中，毛細力效應(yīng)使得液層上部的流胞流動增強，如圖8(d)所示，液層上部流胞幾乎占據(jù)液池的1/3區(qū)域。當熱毛細Reynolds數(shù)為ReT=400時，A=0.15液層內(nèi)有8個流胞（如圖5所示），A=0.5和A=1液池中的流胞個數(shù)分別為2和1，可以推斷，淺液池中流動受到空間限制而流胞數(shù)量多，相比之下，深液池內(nèi)的流動比淺液池中的流動更穩(wěn)定。

3 結(jié) 論

通過對常重力條件下、表面張力比Rσ= -1的環(huán)形液池內(nèi)二元混合溶液雙擴散毛細對流的三維數(shù)值模擬，可以得到以下結(jié)論。

（1）由于熱浮力和溶質(zhì)浮力不平衡，使得液池中只要有溫差就會出現(xiàn)二維軸對稱穩(wěn)態(tài)流動。

（2）流動失穩(wěn)后將轉(zhuǎn)變?yōu)槎S軸對稱周期性振蕩流動，A=0.15時流動轉(zhuǎn)變的臨界熱毛細Reynolds數(shù)為ReTc=231.6，臨界頻率為Fc=10.6；隨著ReT繼續(xù)增大，流動會進一步轉(zhuǎn)變?yōu)槿S周期性振蕩流動。

（3）當保持ReT不變時，液池內(nèi)流體的流動隨深寬比的改變將分別出現(xiàn)周期振蕩、亞臨界非穩(wěn)態(tài)以及三維穩(wěn)態(tài)等流動狀態(tài)。

[1] Abbasoglu S, Sezai I. Three-dimensional modeling of melt flow and segregation during Czochralski growth of GexSi1-xsingle crystals [J].International Journal of Thermal Sciences, 2007, 46:561-572.

[2] Ghorayeb K, Mojtabi A. Double diffusive convection in a vertical rectangular cavity [J].Physics of Fluids, 1997, 9(8): 2339-2348.

[3] Ghorayeb K, Khallouf H, Mojtabi A. Onset of oscillatory flows in double-diffusive convection [J].International Journal of Heat and Mass Transfer, 1999, 42(4): 629-643.

[4] Xin S, Le Quéré P, Tuckerman L S. Bifurcation analysis of double-diffusive convection with opposing horizontal thermal and solutal gradients [J].Physics of Fluids, 1998, 10(4): 850-858.

[5] Bardan G, Bergeon A, Knobloch E, Mojtabi A. Nonlinear doubly diffusive convection in vertical enclosures [J].Physica D:Nonlinear Phenomena, 2000, 138(1): 91-113.

[6] Nishimura T, Wakamatsu M, Morega A M. Oscillatory double-diffusive convection in a rectangular enclosure with combined horizontal temperature and concentration gradients [J].International Journal of Heat and Mass Transfer, 1998, 41(11): 1601-1611.

[7] Bergeon A, Knobloch E. Natural doubly diffusive convection in three-dimensional enclosures [J].Physics of Fluids, 2002, 14(9):3233-3250.

[8] Sezai I, Mohamad A A. Double diffusive convection in a cubic enclosure with opposing temperature and concentration gradients [J].Physics of Fluids, 2000, 12(9): 2210-2223.

[9] Chen Z W, Li Y, Zhan J M. Double-diffusive Marangoni convection in a rectangular cavity: onset of convection [J].Physics of Fluids,2010, 22(3): 034106.

[10] Li Y R, Zhou Y L, Tang J W, Gong Z X. Two-dimensional numerical simulation for flow pattern transition of thermal-solutal capillary convection in an annular pool [J].Microgravity Science and Technology, 2013, 25(4): 225-230.

[11] Zhan J M, Chen Z W, Li Y S, Nie Y H. Three-dimensional double-diffusive Marangoni convection in a cubic cavity with horizontal temperature and concentration gradients [J].Physical Review E, 2010, 82(6): 066305.

[12] Li Y S, Chen Z W, Zhan J M. Double-diffusive Marangoni convection in a rectangular cavity: transition to chaos [J].International Journal of Heat and Mass Transfer, 2010, 53(23):5223-5231.

[13] Lappa M. Thermal convection and related instabilities in models of crystal growth from the melt on earth and in microgravity: past history and current status [J].Crystal Research and Technology, 2005,40(6): 531-549.

[14] Li Y R, Peng L, Akiyama Y, Imaishi N. Three-dimensional numerical simulation of thermocapillary flow of moderate Prandtl number fluid in an annular pool [J].Journal of Crystal Growth, 2003, 259(4):374-387.

[15] Yu J J, Ruan D F, Li Y R, Chen J C. Experimental study on thermocapillary convection of binary mixture in a shallow annular pool with radial temperature gradient [J].Experimental Thermal and Fluid Science, 2015, 61: 79-86.

[16] Li Y R, Zhou Y L, Tang J W, Gong Z X. Two-dimensional numerical simulation for flow pattern transition of thermal-solutal capillary convection in an annular pool [J].Microgravity Sci. Technol., 2013,25: 225-230.

[17] Li Y R, Peng L, Shi W Y, Imaishi N. Convective instability in anular pools [J].FDMP, 2006, 2(3): 153-165.

[18] Smith M K, Davis S H. Instabilities of dynamic thermocapillary liquid layers(Part 1): Convective instabilities [J].Journal of Fluid Mechanics, 1983, 132: 119-144.

[19] Smith M K. Instability mechanisms in dynamic thermocapillary liquid layers [J].Physics of Fluids, 1986, 29(10): 3182.