毛细力比为-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.