高径比对GaAs 熔体液桥热毛细对流失稳的影响1)

周 游 曾 忠, 刘 浩 张良奇

* (重庆大学航空航天学院,重庆 400044)

† (重庆交通大学西南水运工程科学研究所,重庆 400016)

引言

热毛细对流是沿流体-流体界面非均匀温度导致的表面张力梯度驱动的对流,其广泛存在于晶体生长[1-2]、液滴迁移[3-4]、3D 打印[5-6]等工业应用中,其中晶体生长是其重要的研究背景之一.浮区法是一种无坩埚高品质单晶生长技术,随着航空航天技术的发展,空间微重力环境为浮区法生长更大尺寸和更高品质的单晶提供了可能.在空间微重力环境中,热浮力对流极度衰减,热毛细对流成为浮区法晶体生长中熔体物质与热输运的主要方式.由两块不同温度同轴圆盘以及圆盘之间的液柱组成液桥模型,其可当成源于浮区法晶体生长技术的半浮区简化模型.近40 年来,液桥模型已被广泛运用于热毛细对流研究.

熔体普朗特数对液桥热毛细对流失稳机制和失稳模式有重要影响.Smith 和Davis[7]基于线性稳定性分析研究了不同普朗特数无限大平面液层中热毛细对流的不稳定性,发现了2 种对流失稳类型.Levenstam等[8-9]通过三维数值模拟和实验对低普朗特数熔体液桥热毛细对流研究时发现,热毛细对流第一次失稳是从二维轴对称定常流动转变为三维定常对流.Chen 等[10]采用微扰动法研究了低普朗特数液桥热毛细对流失稳,研究表明:当Pr≤0.1 时,液桥热毛细对流的第一次失稳为静态失稳.Ichiro 等[11]通过实验对高普朗特数熔体液桥热毛细对流进行了研究,研究表明随着温差的增大,热毛细对流第一次失稳是从二维轴对称定常流动转变为三维振荡对流.Wanschura 等[12]通过线性稳定性分析研究了不同普朗特数对液桥热毛细对流的影响,结果表明热毛细对流失稳存在两种失稳类型:当Pr≤0.05 时,液桥热毛细流动失稳从二维轴对称定常对流转变为三维定常对流,当0.5≤Pr≤4.8 时,液桥热毛细对流失稳从二维轴对称定常对流转变为三维振荡对流,两种失稳模式具有不同失稳机制.Chen 等[13]采用与Wanschura 等[12]同样的液桥模型,通过线性稳定性分析研究了10-10≤Pr≤8.0 的热毛细对流失稳,结果表明:当Pr≤0.06 时,热毛细对流失稳是由水动力学惯性机制引起,并且其失稳模式为静态失稳;当Pr≥0.1 时,失稳是由热毛细机制引起,并且其失稳模式为振荡失稳.2001 年,Levenstam 等[14]结合数值模拟和线性稳定性分析研究了0.001≤Pr≤7 范围内液桥热毛细流的失稳机制,他们的研究表明:当0.001≤Pr≤0.05 时,热毛细对流失稳是由水动力学惯性机制引起,其失稳模式为静态失稳,当0.057≤Pr≤0.068,热毛细对流失稳主要是由水动力学惯性机制引起,但其失稳模式转变为振荡失稳.当0.07≤Pr≤1.80 时,热毛细对流失稳是由水动力学惯性机制和热毛细机制共同作用引起,其失稳模式为振荡失稳,当0.183≤Pr≤7 时,热毛细对流失稳是主要是由热毛细机制引起,其失稳模式为振荡失稳.上述文献结果表明:液体普朗特数对液桥热毛细对流的失稳机制和失稳模式具有重要影响.对高普朗特数液体,热毛细对流第一次失稳是从二维轴对称定常对流转变为三维周期性振荡对流,其失稳机制是热毛细失稳机制;对低普朗特数液体,热毛细对流第一次失稳是从二维轴对称定常对流转变为三维定常对流,失稳机制为水动力学惯性失稳机制.

针对锡熔体(Pr=0.009)的热毛细对流,Li 等[15]数值模拟研究了液桥高径比对热毛细对流失稳的影响,结果显示:当高径比在0.6～2.2 之间时,热毛细对流失稳临界Marangoni 数(Mac)随着高径比的减小而单调增加,其热毛细对流失稳模式均为静态失稳.Rybicki 和Floryan[16]数值模拟研究了液桥的热毛细对流,结果表明高径比对液桥热毛细流动有重要影响.Velten 等[17]通过地面实验研究了高径比对液桥热毛细对流的影响,结果表明随着高径比的增加Mac逐渐减小.Chen 和Hu[18]通过线性稳定性分析研究了液桥的高径比对半浮区液桥热毛细对流的影响,随着高径比的增加Mac数和临界频率随之减小.Nishino 等[19]结合空间实验和线性稳定性分析方法,研究了高径比对高普朗特数液桥热毛细对流不稳定性的影响,其研究结果表明:当Pr=67 时,在高径比As=1.25 左右,临界频率和Mac出现突降现象,其失稳模式为振荡失稳.吴勇强等[20]通过地面实验研究了液桥高径比对液桥起振温差的影响,研究表明随着高径比的增大液桥内的起振温差逐渐呈减小趋势.王佳等[21]对大尺寸液桥的浮力-热毛细对流进行了地面实验,实验发现在高普朗特数情况下,在Mac附近,流场内会有行波现象出现,随着高径比的变化,其流动模式发生改变.Kang 等[22]在天宫二号上进行了高普朗特数液体桥热毛细对流的实验研究,以液桥高径比和体积比为特征,研究了几何参数对液桥热毛细对流失稳临界条件的影响,首次完整地得到了微重力条件下高径比-体积比参数空间内热毛细对流不稳定性的临界条件和振荡特性.Wang等[23]基于地面实验,研究了存在重力场作用时,几何参数(高径比和体积比)对大普朗特数大尺寸液桥热毛细对流失稳临界条件的影响,在不同的几何参数下出现了6 种不同的失稳状态.Liu 等[24]通过线性稳定性分析研究了高径比对环形液池热毛细对流(Pr=0.011 和Pr=1.4)的影响,研究发现:当高径比较小时,失稳模式为振荡失稳;当高径比较大时,失稳模式为静态失稳,这表示高径比对环形液池热毛细对流的失稳模式转变有影响,但是该现象尚未在液桥热毛细对流研究的文献中报道.

目前,虽然文献已有很多高径比对液桥热毛细流失稳影响的研究,但研究主要聚焦于高普朗特数熔体和低普朗特数熔体的热毛细对流失稳,而较少关注普朗特数在0.057≤Pr≤0.8[14]范围液桥的热毛细流失稳,特别缺乏这一范围内高径比对其热毛细失稳影响的研究.本文针对GaAs 熔体(Pr=0.068),采用谱元法对基态解进行求解,然后基于线性稳定性分析确定失稳临界参数Mac,最后结合能量分析方法研究对流失稳的物理机制,研究高径比对液桥热毛细对流失稳临界参数、失稳模式和失稳机制的影响.

1 数学模型和数值方法

1.1 物理模型与基本控制方程

采用图1 所示的液桥模型,液桥由中心轴线重合的圆盘和圆形液柱组成,圆柱形液桥高度为H,半径为R,高径比定义为As=H/R.液桥上端为低温固壁(φc),下端为高温固壁(φh),熔体当成不可压缩牛顿流体.假设表面张力为温度的线性函数,圆柱自由表面绝热并且忽略自由表面变形.

图1 半浮区液桥模型Fig.1 Floating half-zone model

分别用H,υ/H,H2/υ 和 ρυ2/H2作为长度、速度、时间和压力的特征物理量,对流体的控制方程进行无量纲化,其中 υ 是运动学黏性系数,ρ 是密度,无量纲温度定义为 T=(φ-φc)/(φh-φc),其中φc为液桥低温圆盘固壁温度,φh为高温圆盘固壁高温,φ表示表示液桥内任意一点温度.热毛细对流失稳前是二维轴对称定常流,失稳后变为三维对流.采用柱坐标(r,θ,z),r 为径向坐标,θ 为周向坐标,z 为轴向坐标.无量纲控制方程为

无量纲边界条件可以表示为:

(1) 自由表面(r=1/As)

(2) 热壁(z=0)

(3) 冷壁(z=1)

其中u 表示径向速度,w 表示轴向速度,无量纲参数Pr 以及Ma 分别定义为

其中 κ 为热扩散系数,γT为流体表面张力系数,Δφ=φh-φc.

1.2 线性稳定性分析

为了得到高精度的失稳临界条件,采用基于时间分裂法的谱元法求解方程式(1)～式(4),得到定常轴对称的基态解,然后在所得二维基态解上施加一个三维小扰动[25]

其中 u0,p0和 T0表示速度、压力和温度的基态解,表示扰动速度、扰动压力和扰动温度.将式(7)代入到方程式(1)～式(4)中,并忽略高阶扰动项可导出扰动控制方程

扰动边界条件为:

(1) 自由表面(r=1/As)

(2) 热壁(z=0)

(3) 冷壁(z=1)

将扰动量表示为正则模形式

其中φ 表示u,v,w,T 和p,c.c.为复数虚部项的共轭.复数λ 可分为实部λr和虚部λi,其中实部λr表示为线性增长率,虚部λi表示为振荡频率ƒ,k 表示为失稳波数.其中λr＜0 表示流动稳定,λr＞0 表示流动不稳定,而λr=0 则代表流动失稳的临界状态,线性稳定性分析正是要捕捉λr=0 时的临界参数,在临界条件下,λi=0 表示失稳是静态失稳,λi≠0 表示失稳是振荡失稳.在将表达式(12)代入式(8)～式(11),采用谱元法对扰动方程进行离散,得到一个广义特征值问题,即

1.3 扰动能量分析

确定流动失稳临界参数后,可通过扰动能量分析的方法来研究流动失稳机制.对扰动方程(8)中的动量方程乘以û并对求整个求解域进行积分,可得到扰动动能变化率表达式[12,30-31]

其中 ∂ekin/∂t 表示扰动动能的分布,Dk表示黏性耗散率,Mz和Mφ分别表示由扰动引起的热毛细力沿液桥轴向和周向做功的功率,IV项表示从基态流场到扰动流场的动能传递率.则各个能量项的具体表达式为

其中

其中IV1,IV2和IV3表示基态径向速度向扰动场的能量传递,而IV4和IV5则表示基态轴向速度向扰动场的能量传递.V 表示液桥的体积,S 表示自由表面.

2 结果与讨论

2.1 网格依赖性和程序有效性验证

选取Pr=0.068,液桥高径比As=1,计算网格采用8r× 9z单元划分,单元内分别采用5r× 5z,6r× 6z,7r× 7z阶谱离散,对应网格节点为41 × 46,49 × 55,57 × 64.将不同网格的计算结果进行对比,如表1 所示,不同网格计算的波数吻合,Mac相对误差较小.综合考虑计算精度和计算效率,本文以下计算采用6r× 6z阶谱离散.

表1 基于不同网格分辨率得到的热毛细对流(Pr=0.068)失稳临界Marangoni 数(Mac)和波数kcTable 1 The critical Marangoni number (Mac) and wave number kc for the instability of thermocapillary flow with Pr=0.068 under different mesh resolution conditions

选取液桥高径比As=1 的热毛细对流进行程序有效性验证,计算网格采用8r× 9z单元划分,6r×6z阶谱离散,对应网格节点为49 × 55.将计算Mac转化为Rec(Rec=Mac/Pr)与Levestam 等[14]的结果进行对比,如表2 所示.

表2 计算失稳临界参数(Rec,fc 和kc)与文献[14]结果对比Table 2 The critical values (Rec,fc,kc) from the present computations against Ref.[14]

针对Pr=0.009 熔体热毛细对流的失稳,基于我们自主开发谱元法程序计算结果与Li 等[15]的临界参数(失稳临界Marangoni 数(Mac)和波数kc)进行对比,如表3 所示.

表3 计算失稳临界参数(Mac 和kc)与文献[15]结果对比Table 3 The critical values (Mac and kc) from the present computations against Ref.[15]

以上与Levestam 等[14]和Li 等[15]的结果对比表明:本文计算热毛细对流失稳临界参数与他们结果吻合,我们程序的有效性得到了验证.

2.2 稳定性曲线与失稳模式

GaAs 熔体液桥热毛细对流的基态是二维轴对称定常对流.固定液桥的高度为1,液桥半径的变化范围为0.4≤R≤2.5,则对应液桥的高径比变化范围为2.5≥As≥0.4.离散采用6r× 6z阶谱离散,计算网格采用(8r-14r) × 9z单元划分,所以液桥轴向有54 个网格,径向根据高径比的不同有48～84 个网格.

如表4 所示,失稳波数随着高径比的增加而减小,失稳波数的变化符合Zeng 等[32]提出的高径比和失稳波数之间的关系1.57≤kcAs≤3.14.如图2 所示,临界Mac数随着高径比增加,整体呈下降趋势,Mac数在高径比为0.67 和高径比为1.18 处出现峰值点.表4 中的线性稳定性分析结果表明,在高径比0.40≤As≤2.50 区间内,存在两种失稳模式,当0.4≤As≤1.18 时,热毛细对流从二维轴对称定常对流转捩为三维振荡对流;当1.20≤As≤2.5 时,热毛细对流从二维轴对称定常对流转捩为三维定常对流.

图2 Mac 随高径比(As)的变化(Pr=0.068)Fig.2 The variation of Mac versus the aspect ratio As for Pr=0.068

表4 Pr=0.068 时不同高径比下的失稳临界参数Table 4 The critical values at different aspect ratios for Pr=0.068

在图3 中 IV项表示IV1项到IV5项的总和,M 项表示Mz和Mφ的总和.由图3 可知,失稳主要是由水动力学惯性机制主导.结合表5 扰动能量平衡分析的结果可知,IV4和IV5是造成流动失稳的主要原因,且Mz和Mφ也占据一定比例,表明失稳是由水动力学惯性机制占主导,同时热毛细机制对失稳有所贡献.在高径比As=1.82 处,由水动力学惯性机制主导的IV4+IV5的占比最大,而热毛细机制引起的Mz+Mφ的占比最小,Levestam 等[14]针对高径比As=1熔体液桥热毛细对流失稳研究表明:当Pr=0.068时,液桥热毛细对流是振荡失稳,其失稳机制是水动力学的惯性失稳机制,这与本文的分析结果一致.下面具体介绍随高径比变化的两种失稳模式.

图3 扰动能量随高径比As 变化Fig.3 The variation of the kinetic energy budgets with respect to the aspect ratio

图3 扰动能量随高径比As 变化(续)Fig.3 The variation of the kinetic energy budgets with respect to the aspect ratio (continued)

表5 Pr=0.068 时不同高径比能量分析结果Table 5 The energy analysis results at different aspect ratios for Pr=0.068

2.3 失稳机制

2.3.1 失稳模式I (振荡失稳)

由表4 可知,当0.40≤As≤1.18 时,失稳临界频率不等于零,即热毛细对流从定常二维轴对称流动失稳变为三维周期性振荡对流.随着高径比的增加,失稳的临界波数从kc=6 逐渐减少为kc=2.由表5和图3 可知,在失稳模式I 中,失稳主要是由水动力学惯性机制作用,其中IV4和IV5是造成水动力学惯性失稳的主要原因;在这一失稳模式中,热毛细机制对失稳也有一定贡献.通过具体分析,失稳模式I 在不同高径比范围内,IV4项和IV5项对失稳的贡献略微有所不同.当0.40≤As≤0.59 时,液桥流动失稳的能量主要是由IV4和IV5提供(约占100%),且二者比例相当;IV2和IV3为负值,代表其对对流具有稳定作用;同时扰动引起的热毛细力做功(M 项)的贡献占比小于8%(表5).当0.67≤As≤1.00 时,IV5项对液桥流动失稳的贡献在60%左右,是造成流动失稳的主要原因;IV2已经为正值,IV3仍然为负值,热毛细力做功的贡献占比为6.52%～9.86%.当1.05≤As≤1.18 时,液桥流动失稳主要是由IV4项引起,其贡献在65%左右,在这一高径比区间Mz和Mφ的对失稳的贡献有所增加,热毛细力做功的贡献占比为15.31%～27.35%,即热毛细机制有所增强.下面将详细分析上述3 个高径比区间的失稳机制.

当0.40≤As≤0.59 时,由表5 可知,在这一高径比区间中,Mz和Mφ对液桥的流动失稳贡献较小,IV项(IV项表示IV1项到IV5项的总和)是造成流动失稳的主要原因.由图4(a)扰动能量平衡图可以直观的看出,失稳主要由IV4和IV5造成,即基态的轴向速度径向和轴向梯度为引起失稳的扰动动能传递了主要的能量.由局部扰动动能极值图(图4(b))可知,液桥热毛细流动最危险的失稳区域出现在距液桥自由液面R/3 处并靠近热壁附近.在液桥自由液面的扰动温度图(图5(b))上分布了5 个热极和5 个冷极,其中红色表示热极,蓝色表示冷极,自由液面的扰动速度从热极指向冷极,与扰动热毛细力的方向相同.从图5(a) 中可以看出,在靠近自由液面约R/3 处出现扰动温度强斑,截面处的扰动速度从热极指向冷极,并在冷极处形成扰动速度流胞.

图4 高径比As=0.48,Pr=0.068 的临界状态(Mac=1936.3)下:(a)扰动能量平衡条状图,(b)局部扰动动能极值图Fig.4 (a) The disturbance energy balance and (b) distribution of the local disturbance kinetic energy extreme value under critical condition(Mac=1936.3) at As=0.48,Pr=0.068

图5 高径比As=0.48,Pr=0.068 的临界状态 (Mac=1936.3) 下:(a) z=0.51 截面扰动温度 (云图) 和扰动速度矢量 (箭头),(b)自由液面扰动温度 (云图) 和扰动速度矢量 (箭头)Fig.5 The disturbance temperature (colored contour) and disturbance velocity vector (arrow):(a) the z=0.51 cross section,(b) the free surface under the critical condition (Mac=1936.3) at As=0.48,Pr=0.068

在0.67≤As≤1 区域内,由图6 的扰动能量平衡图可知,当高径比在这一范围内,失稳依然是由水动力学惯性机制占主导,但与0.40≤As≤0.59 这一高径比范围内失稳的具体扰动动能项存在一些差异,该高径比范围内液桥的流动失稳主要是由IV5项主导,即基态轴向速度的轴向梯度为流动失稳提供了主要能量.流动最危险的失稳位置出现在靠近热壁处距离自由液面R/2 附近(图6(b)).如图7(a)中扰动温度所示,该高径比范围内扰动温度的特征与0.40≤As≤0.59 时不同,在这一高径比范围内,扰动温度的特征是:在靠近对称轴R/3 处和靠近自由液面附近出现强斑.在自由液面处扰动速度从热极指向冷极,与扰动热毛细力的方向相同(图7(b)).

图6 高径比As=0.77,Pr=0.068 的临界状态 (Mac=1083.9) 下:(a)扰动能量平衡条状图,(b)局部扰动动能极值图Fig.6 (a) The disturbance energy balance and (b) distribution of the local disturbance kinetic energy extreme value under critical condition(Mac=1083.9) at As=0.77,Pr=0.068

图7 高径比As=0.77,Pr=0.068 的临界状态 (Mac=1083.9) 下:(a) z=0.51 截面扰动温度 (云图) 和扰动速度矢量 (箭头),(b)自由液面扰动温度 (云图) 和扰动速度矢量 (箭头)Fig.7 The disturbance temperature (colored contour) and disturbance velocity vector (arrow):(a) the z=0.51 cross section,(b) the free surface under the critical condition (Mac=1083.9) at As=0.77,Pr=0.068

当1.05≤As≤1.18 时,由扰动能量平衡图(图8(a))可知,在这一高径比范围内,热毛细机制对失稳的贡献增大,但液桥热毛细流动失稳仍是由水动力学的惯性机制主导.与上述两种情况不同,在这一高径比范围内IV4项是造成失稳的主要原因,但 IV2项对流动失稳有较为明显的抑制作用,即基态轴向速度的径向梯度为流动失稳提供了主要能量,而基态径向速度的轴向梯度有抑制流动失稳的作用.根据图8(b)可知,最危险的失稳位置出现在热壁附近,并距离对称轴R/3 处.从自由液面展开图(图9(a))显示:在靠近冷壁附近扰动速度从热极指向冷极,而在靠近热壁附近扰动速度从冷极指向热极.根据z=0.20,z=0.51 和z=0.74 截面扰度温度和扰动速度矢量图(图9(b)～图9(d))可知,在z=0.20 截面处扰动速度从冷极指向热极附近,扰动温度模式的特征为轮辐状强斑,并且布满整个液桥横截面;在z=0.51 截面处,扰动速度在靠近对称轴附近是由热热极指向冷极,在靠近自由液面附近是由冷极指向热极,从而扰动速度在热极和冷极之间形成一个扰动流动循环,扰动温度模式特征为:在靠近自由液面附近形成椭圆形强斑;在z=0.74 截面中扰动速度从热极指向冷极,扰动温度特征为:在自由液面附近和靠近自由液面R/3 处形成强斑.

图8 高径比As=1.11,Pr=0.068 的临界状态 (Mac=1273.6) 下:(a)扰动能量平衡条状图,(b)局部扰动动能极值图Fig.8 (a) The disturbance energy balance and (b) distribution of the local disturbance kinetic energy extreme value under critical condition(Mac=1273.6) at As=1.11,Pr=0.068

图9 高径比 As=1.11，Pr=0.068 的临界状态(Mac=1273.6)下扰动温度(云图)和扰动速度矢量(箭头)Fig.9 The disturbance temperature (colored contour) and disturbance velocity vector (arrow) under the critical condition (Mac=1273.6) for the case of As=1.11 and Pr=0.068

2.3.2 失稳模式II (静态失稳)

根据表4 可知,在1.20≤As≤2.50 区域内,失稳临界频率等于零,表示在这一高径比范围内,热毛细对流从二维轴对称定常对流失稳变为三维定常对流.随着高径比的增加,失稳的临界波数kc=2 转变为kc=1.由表5 和图3 可知,液桥的流动失稳依然是由IV项主导,这表明失稳模式II 仍然是水动力学惯性失稳机制主导,其中M 项在失稳中的占比随高径比的增加先增大,再减小最后再增大(图3(a)).在失稳模式II 中IV4和IV5是引起液桥流动失稳的主要原因,但在不同的高径比下 IV4项和IV5项的占比存在差异.在1.20≤As≤1.54 范围内,液桥的流动失稳是由IV4项主导(约占43%),但IV5项对失稳的贡献同样较大(约占28%),M 项的占比在13.39%～22.18%;在1.67≤As≤2.50,IV4项对失稳的贡献增大,而IV5项和M 项的在失稳中的占比减小,同时IV2和IV3变为负值,所以IV4项是引起流动失稳的主要原因.下面对失稳模式II 进行具体分析.

在1.20≤As≤1.54 范围,失稳临界波数kc=2.根据图10(a)所示,IV4和IV5是造成流动失稳的主原因,即基态轴向速度的径向梯度和轴向梯度为引起失稳的扰动动能传递了主要的能量,其中IV4项比IV5项的贡献更大,Mz和Mφ项在流动失稳中的占比也相对较大.流动失稳最危险的位置出现在液桥中部偏上附近并靠近自由液面处(图10(b)).从自由液面展开图(图11(a))可知,扰动速度在靠近冷壁处是由热极指向冷极,在热壁附近则是由冷极指向热极,这与图12 所示的低普朗特数(Pr=0.009)液桥在自由液面处扰动速度的指向不同,在低普朗特数液桥中自由有液面的扰动速度都是从冷极指向热极,且失稳都是水动力学惯性机制作用下的失稳,热毛细机制对失稳没有贡献.由图11(b)～图11(d)可知,z=0.20 截面处扰动温度的特征为:靠近自由液面R/3 形成强斑,截面处的扰动速度由冷极流向热极;在z=0.51 截面处,扰动温度特征为:在自由液面R/3 形成椭圆形强斑,并衍生出来的弱斑向自由液面靠近,截面处的扰动速度在靠近对称轴处是由热级指向冷极,在靠近自由液面附近是由冷极指向热极;在z=0.74 截面处,扰动温度特征为:扰动温度强斑向自由液面靠近,截面处的扰动速度是由热极流向冷极.

图10 高径比As=1.43,Pr=0.068 的临界状态 (Mac=531.9) 下:(a)扰动能量平衡条状图,(b) 局部扰动动能极值图Fig.10 (a) The disturbance energy balance and (b) distribution of the local disturbance kinetic energy extreme value under critical condition(Mac=531.9) at As=1.43,Pr=0.068

图11 高径比As=1.43,Pr=0.068 的临界状态(Mac=531.9)下扰动温度(云图)和扰动速度矢量(箭头)Fig.11 The disturbance temperature (colored contour) and disturbance velocity vector (arrow) under the critical condition (Mac=531.9) for the case of As=1.43 and Pr=0.068

图12 高径比As=1.43,Pr=0.009 临界状态(Mac=531.9)下自由液面处扰动温度(云图)和扰动速度矢量(箭头)Fig.12 The disturbance temperature (colored contour) and disturbance velocity vector (arrow) at the free surface under the critical condition(Mac=531.9) for the case of As=1.43 and Pr=0.009

在1.67≤As≤2.50 区间内,失稳临界波数kc=1.由图13(a)可知,IV4是造成流动失稳的主要原因,即基态轴向速度的径向梯度为流动失稳提供了最主要能量,而Mz也对液桥的流动失稳有所影响,这表明失稳仍然是由水动力学惯性机制主导,同时热毛细机制对失稳有所影响.与上述所有失稳机制不同的是,在该高径比范围内Mφ项起抑制流动失稳的作用.由图13(b)可知,最危险的失稳位置出现在液桥中部并靠近自由液面.在图14(a)中,扰动速度在靠近液桥冷壁处是由热极指向冷极,在靠近液桥热壁处扰动速度从冷极指向热极,在1.05≤As≤1.18 和1.20≤As≤1.54 这两个高径比范围内都曾出现这一现象,但在该高径比范围内这一现象更加明显,扰动速度更多地从冷极指向热极.由图14(a)和图15 对比分析可知,低普朗特数(Pr=0.009)液桥自由液面的扰动速度是从冷极指向热极,这与GaAs 熔体液桥自由液面的扰动速度指向不同,由图14(b)～图14(d)可知,在z=0.20,z=0.51,z=0.74 的横截面的边缘处,扰动速度都是从冷极指向热极,同样验证了自由液面处扰动速度更多的是从冷极指向热极;这三处横截面的扰动温度特征为:随着z 坐标的增大,扰动温度强斑逐渐向自由液面靠近.

图13 高径比 As=2.00，Pr=0.068 的临界状态 (Mac=234.6) 下:(a)扰动能量平衡条状图,(b)局部扰动动能极值图Fig.13 (a) The disturbance energy balance and (b) distribution of the local disturbance kinetic energy extreme value under critical condition(Mac=234.6) at As=2.00,Pr=0.068

图14 高径比As=2.00,Pr=0.068 的临界状态(Mac=234.6)下扰动温度(云图)和扰动速度矢量(箭头)Fig.14 The disturbance temperature (colored contour) and disturbance velocity vector (arrow) under the critical condition (Mac=234.6) for the case of As=2.00 and Pr=0.068

图15 As=2.00 Pr=0.009 时临界状态(Mac=234.6)下自由液面处扰动温度(云图)和扰动速度矢量(箭头)Fig.15 The disturbance temperature (colored contour) and disturbance velocity vector (arrow) at the free surface under the critical condition(Mac=234.6) for the case of As=2.00 and Pr=0.009

2.3.3 基态流函数和温度分布

由图16 和图17 可知,在失稳模式I 和失稳模式II 中,随着高径比的增大基态涡胞的位置始终靠近自由液面附近.从图16 和图17 可以看出,当高径比较小时,温度场在靠近自由液面R/3 处向热壁处挤压,同时温度场也在自由液面处向冷壁处挤压,从而形成较大温度梯度;但随着高径比增加,这种温度场的挤压相对减缓,从而随着高径比的增加温度场沿轴向分布相对均匀.

图16 失稳模式I 在临界状态下的基态流函数(云图)和基态温度场(等值线)Fig.16 The basic flow (colored contour) and temperature (isotherm) of instability mode I under critical condition

图17 失稳模式II 在临界状态下的基态流函数(云图)和基态温度场(等值线)Fig.17 The basic flow (colored contour) and temperature (isotherm) of instability mode II under critical condition

图17 失稳模式II 在临界状态下的基态流函数(云图)和基态温度场(等值线) (续)Fig.17 The basic flow (colored contour) and temperature (isotherm) of instability mode II under critical condition (continued)

3 结论

基于已有文献可知,液桥模型中熔体普朗特数对热毛细对流失稳机制和失稳模式具有重要影响:对典型高普朗特数熔体(例如Pr＞1),热毛细对流第一次失稳是从二维轴对称定常对流转变三维周期性振荡对流,其失稳机制是热毛细失稳机制;对典型低普朗特数熔体(例如Pr=0.011),热毛细对流第一次失稳是从二维轴对称定常对流转变为三维定常对流,失稳机制是水动力学惯性失稳.同时文献工作显示液桥高径比对热毛细流对流的失稳临界参数具有重要影响,但是未发现高径比会影响热毛细对流失稳机制和失稳模式(三维静态失稳和三维振荡失稳).

本文在基于谱元法线性稳定性分析对GaAs 熔体(Pr=0.068)液桥热毛细对流研究发现,在同一普朗特数不同高径比的情况下,随着高径比的改变液桥出现两种失稳模式:当0.4≤As≤1.18 时,热毛细对流是由二维轴对称定常对流直接转变为三维周期性振荡对流(振荡失稳);当1.20≤As≤2.5 时,热毛细对流是由二维轴对称定常对流转变为三维定常流动(静态失稳).基于能量分析我们发现,GaAs 熔体液桥的热毛细对流失稳是由水动力学惯性机制占主导,两种机制(水动力学惯性机制和热毛细机制)共同作用造成.两种机制对失稳的贡献随着高径比的变化而变化,其失稳现象更加复杂,从而随着高径比的改变出现了两种不同的失稳模式.由IV项各项的占比可知,IV4和IV5是造成流动失稳的主要原因,即基态轴向速度径向梯度和轴向梯度为扰动场的失稳转递了主要的能量.在失稳模式I(振荡失稳)中,最危险的失稳位置出现在热壁附近,并随着高径比的增加逐渐向对称轴靠近;在失稳模式II(静态失稳)中,最危险的失稳位置液桥中部横截面并靠近自由液面附近.