毛細力比為-1時環形液池內雙擴散毛細對流數值模擬

陳捷超，李友榮，于佳佳

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

引 言

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

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

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

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

1 物理數學模型

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

圖1 物理模型Fig.1 Physical model

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

取量綱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坐標。

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

式中，ν為動量擴散率。

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

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

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

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

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

2 計算結果與討論

2.1 穩態流動

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

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

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

2.2 臨界熱毛細Reynolds數

隨著水平溫差（濃度差）的增大，液池內的流動會由二維軸對稱穩態流動經 Hopf分岔轉變為二維軸對稱周期性振蕩流動，轉變點對應的熱毛細Reynolds數稱為臨界熱毛細Reynolds數ReTc。圖3給出了A=0.15時監測點P（R=1.5,θ=0,Z=0.145）處徑向速度振幅|VR,M|和量綱 1頻率F隨熱毛細Reynolds數ReT的變化，采用線性外推法可得臨界熱毛細 Reynolds數為ReTc=231.6，臨界頻率為Fc=10.6。

2.3 周期性振蕩流動

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

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

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

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

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

2.4 深寬比的影響

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

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

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

圖7 ReT=800時一個周期（τp）內自由表面溫度波動[(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

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

3 結 論

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

（1）由于熱浮力和溶質浮力不平衡，使得液池中只要有溫差就會出現二維軸對稱穩態流動。

（2）流動失穩后將轉變為二維軸對稱周期性振蕩流動，A=0.15時流動轉變的臨界熱毛細Reynolds數為ReTc=231.6，臨界頻率為Fc=10.6；隨著ReT繼續增大，流動會進一步轉變為三維周期性振蕩流動。

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

[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.