液固两相流体热毛细对流中颗粒动态积累结构研究

黄鑫，梁儒全,2，范俊庚

(1.东北大学材料电磁过程研究教育部重点实验室，辽宁沈阳，110819；2.临沂大学机械与车辆工程学院，山东临沂，276000)

采用浮区法制备晶体过程中，在流体流动的驱动力为浮力、毛细力等。微重力环境下，浮力对流大大减弱，热毛细对流成为主导对流形式。研究表明，振荡热毛细对流会导致晶体表面和内部产生条纹，影响晶体的质量。因此，在微重力环境下，半导体熔体内热毛细对流及流动结构研究已成为空间材料生长的重要课题，对空间浮区法制备高质量晶体具有重要意义。

液桥是浮区法制备晶体简化模型。1979年，CHUN等[1]利用“光切割”技术，开展了半浮区液桥热毛细对流实验，发现了振荡热毛细对流。之后，针对液桥模型，国内外开展了线性稳定性分析[2-3]、三维数值模拟[4-6]和晶体生长实验[7-9]。ZENG等[10]数值研究了三维振荡热毛细对流，结果表明热毛细对流存在脉动振荡和旋转振荡两种不同的振荡模式，温度场随时间周期性变化。YASUHIRO 等[11]通过数值模拟再现了实验观察到的不同频率的超临界熔体自由表面温度振荡，结果表明当温差达到一定的临界值时，表面温度振荡将向三维稳定周期振荡对流转变。当温差为更高的临界值时，表面温度振荡将会发生向三维振荡对流的第二次转变。

另外，浮区法制备晶体过程中浮区熔体中存在杂质，杂质颗粒具有动态积累特性，影响制备晶体质量。研究表明，在一定高径比条件下，稀密度的小颗粒在液桥中聚集，并具有时间依赖性。这些颗粒绕热毛细对流的轴旋转呈变形的螺旋，形成颗粒动态积累结构(PAS)[12]。SCHWABE 等[13]发现了PAS。UENO 等[14]也观察到了类似的结果，同时得到闭合螺旋绕组数等于方位角波数。SCHWABE等[15]通过改变等密度下的颗粒直径和等粒径下颗粒密度与流体密度的比值，测量PAS 结构的形成时间，解释了PAS 的形成机制，即颗粒与流体的密度差影响PAS 的形成速度。MULDOON等[16]的研究表明，密度差和颗粒-自由表面相互作用共同影响PSA的形成。

目前，对PAS 的研究大多基于单向耦合，考虑三向耦合(同时考虑流体对颗粒的作用、颗粒对流体的反作用、颗粒之间的相互作用)的PAS 研究缺乏，因此，对PAS的形成机制仍不清楚。另外，国内外关于热毛细对流的研究主要集中在单一流体热毛细对流研究，对夹杂颗粒的固液两相流体热毛细对流研究极其有限，对热毛细对流的周期性变化与热毛细对流中颗粒聚集规律之间相互关系的研究很少。采用浮区法制备单晶硅过程中熔区内不可避免含有杂质颗粒，针对这种夹杂颗粒的固液两相流体热毛细对流研究更具有理论意义与实际应用价值。

本文作者基于一种三向耦合的CFD-DEM 方法，同时考虑流体对颗粒的作用、颗粒对流体的反作用、颗粒之间的相互作用，采用CFD 方法计算流体的流动行为，利用DEM方法处理颗粒运动和颗粒间的碰撞，对夹杂颗粒流体热毛细对流进行了数值模拟，探究不同高径比下的热毛细对流和PAS，分析其形成原因。

1 物理与数学模型

本研究采用的物理模型如图1所示。2 个半径相同的同心圆盘之间悬浮的是半浮区液桥，液桥的高度为L，半径为a，液桥高径比为L/a。上下圆盘之间存在温差，温度梯度的方向与重力的方向相反。

图1 液桥模型Fig.1 Liquid bridge model

考虑流体对颗粒的作用、颗粒对流体的反作用、颗粒之间的相互作用，三向耦合下液相的控制方程[17-23]如下。

连续性方程为

动量方程为

能量方程为

根据牛顿第二定律，单个颗粒的转动和平动方程分别为：

单个颗粒的能量守恒方程为

系统的边界条件为：下部圆盘(z=-L/2)处温度为T=Tc，上部圆盘(z=L/2)处温度为T=Th，上下圆盘绝热无滑移。自由面的径向温度梯度∂T/∂r=0。在自由面(R=a)上，考虑热毛细效应，采用下列方程：

式中：r，θ和z为坐标；ur，uθ和uz分别为r，θ和z这3个方向上的速度分量；ρf，uf，p，Tf，cf和Г分别为流体的密度、速度、压力、温度、比热容和热扩散系数；ε，t，τ，g和σ分别为孔隙率、时间、气相应力张量、重力加速度和表面张力；mi，vi，ωi，Ii，Ti和cp,i分别为颗粒i的质量、线速度、角速度、转动惯量、温度和比热容；ffp为颗粒-流体作用力；Qf,i为颗粒与流体间的热流量；ff,i为颗粒-流体作用力；fe,ij和fd,ij分别为弹性力矢量与黏性力矢量；Tn,ij，Tt,ij和Tr,ij分别为法向力矩、切向力矩和滚动摩擦力矩；Qi,j和Qi,f分别为颗粒间的换热量和颗粒与流体的换热量；MaT为热Marangoni 数，MaT=rTTL/(μυ)；Pr为普朗特数，Pr=υ/α；ΔT为上下板间的温差，ΔT=Th-Tc。

对3组不同高径比下的含夹杂物的两相流体热毛细对流进行模拟，所用的流体介质为0.65 号硅油，Pr=6.7。采用有限体积法对控制方程进行离散，通过PISO 算法求解控制方程。采用CFDDEM下的Eulerian-Lagrangian模型，模拟含夹杂物的两相流。颗粒间的作用模型和颗粒与自由面的作用模型均为Ranz-Marshall 模型，模拟所用的颗粒为直径70 μm的铝球。在计算的起始时刻，液桥中均匀注入所有颗粒。计算所需的参数见表1。3种高径比下均在液桥的(0.99a，0，0)处设置监测点。

表1 流体物理性质Table 1 Fluid physical property

为了验证网格的相关性，在高径比为1，ΔT=10 K 和MaT=12 457 时，分别采用34r×36z×60θ，34r×38z×80θ和36r×40z×100θ这3 种网格进行计算，监测点(0.99a，0，0)处计算结果如表2所示。当网格数为34r×38z×80θ和36r×40z×100θ时，监测点(0.99a，0，0)处温度和速度的相对误差均小于1%，因此，选用网格数更少的网格34r×38z×80θ作为数值计算网格。

表2 网格无关性检验Table 2 Grid dependence test

2 结果与讨论

2.1 模拟结果验证

首先，将本研究的结果与SCHWABE等[12]的实验结果以及MELNIKOV 等[24]的数值模拟结果进行比较验证。当ΔT=8.5 K，Pr=13.5，L/a=1 时，在计算的起始时刻往液桥中均匀注入颗粒并施加温差，一段时间后液桥中出现稳定的温度和速度振荡，再经历一段时间后，颗粒呈现出动态积累，图2(a)所示为本研究模拟的液桥注入颗粒250 s 后的俯视图。颗粒分布与图2(b)中SCHWABE等[12]的实验结果和图2(c)中MELNIKOV 等[24]的模拟结果非常相似，也呈不对称的单环，PAS绕z轴逆时针旋转。图3所示为方位角波数模拟结果与SCHWABE 等[12]的实验结果对比。当ΔT=10 K，Pr=8 时，改变高径比，数值模拟测得对应的方位角波数。由图3可以看出：方位角波数数值模拟结果与实验结果相吻合，随着高径比增大，两者变化趋势相同，验证了本研究建立的模型与计算方法的有效性。

图2 颗粒动态积累结构对比Fig.2 Comparison of dynamic particle accumulation structure

图3 方位角波数模拟结果与实验结果对比Fig.3 Comparison of azimuthal wave number between simulation results and experimental results

2.2 热毛细对流中颗粒动态聚集特性

图4所示为液桥模型中水平面(z=0)的温度分布。从图4可见：水平面中存在着2，3和4倍的对称结构，对称结构与方位角波数m对应；当上下圆盘的温差ΔT=10 K，高径比为0.50，0.72 和1时，方位角波数m分别为4，3和2；随着高径比增大，方位角波数减少，水平面温度分布中低温冷区的数目减少，同时可以得到(L/a)·m≈2。图5所示为不同高径比液桥颗粒聚集结构的俯视图。从图5可见：中心出现颗粒较少区，形状分别为梭形(m=2)、三角形(m=3)、正方形(m=4)；在这个中心区域，颗粒速度较低。这个区域绕z轴以恒定的角速度旋转，旋转过程中形状保持不变。这种粒子较少的多边形区域的形状也与液桥的高径比有关。粒子链组成了多边形的边，多边形的尖端连接液桥的顶部和底部m次。中心区域的形状对应着温度场中的m倍对称结构。图6所示为不同高径比液桥中颗粒分布的鸟瞰图。从图6可见：颗粒聚集结构并不局限于1个平面内，而是分布在整个三维空间；当PAS 形成后，自由面和液桥体积的很大部分没有颗粒。颗粒沿闭合的环路积累，该环路在液桥中作三维螺旋运动，螺旋的形状不变；这个粒子链从热端开始，沿自由面向下，再随回流上升，直到再次接触到自由面；带着粒子的回流并没有穿过液桥中心，粒子链在自由面附近是陡峭的，在回流中则更加水平。粒子链不同的陡峭程度反映着颗粒不同的运动速度，即颗粒在自由表面附近速度更大，在回流中速度较小。

图4 ΔT=10 K和MaT=12 457时不同高径比液桥水平面温度分布Fig.4 Temperature distribution in horizontal cross section at different aspect ratios when ΔT=10 K and MaT=12 457

图5 ΔT=10 K和MaT=12 457时不同高径比的液桥颗粒聚集结构的俯视图Fig.5 Top view of PAS at different aspect ratios when ΔT=10 K and MaT=12 457

图6 ΔT=10 K和MaT=12 457时不同高径比下颗粒分布的鸟瞰图Fig.6 Bird′s-eye view of PAS at different aspect ratios when ΔT=10 K and MaT=12 457

2.3 热毛细对流旋转振荡特性

图7所示为L/a=0.72(m=3)时，1个周期内水平面z=0处温度云图和点(0.99a，0，0)处的温度变化曲线。从图7可见：水平面z=0 处温度云图绕z轴旋转，旋转的周期Tmin=11.25 s。图7中的曲线为m=3时，监测点(0.99a，0，0)处温度变化曲线，振荡周期T′min=3.75 s。可以得到Tmin=3·T′min。当方位角波数m=2，3，4时，监测点处温度分别为T2，T3和T4，速度分别为u2，u3和u4。图8(a)所示为3种高径比液桥中监测点(0.99a，0，0)处温度变化曲线。可见3个监测点的温度出现了稳定的振荡，监测点温度T2

图7 L/a=0.72，ΔT=10 K和MaT=12 457时，1个周期内的温度云图和点(0.99a，0，0)处温度变化曲线Fig.7 Temperature field within one period and temperature-time curve at point(0.99a,0,0)when L/a=0.72,ΔT=10K and MaT=12 457

图8 ΔT=10 K和MaT=12 457时，不同高径比下点(0.99a，0，0)处温度-时间变化曲线和速度-时间变化曲线Fig.8 Temperature-time curves and velocity-time curves at points(0.99a,0,0)at different aspect ratios when ΔT=10 K and MaT=12 457

图9所示为高径比L/a=0.72(方位角波数m=3)时，1个周期内液桥中颗粒聚集结构的俯视图。在这个周期内，PAS绕z轴逆时针旋转。为清楚地展示PAS的旋转，三角形的1个顶点用圆点标记。旋转1周的时间为11.25 s，与液桥温度场的变化周期相同，所以，液桥的周期性旋转振荡与PAS 的转动相对应，PAS的旋转周期与液桥温度场的旋转周期相同，但两者的旋转方向相反。所有颗粒协同运动，造成PAS 随热毛细对流以相同的角速度旋转，这种运动趋势仅适用于整个PAS，与单个颗粒的运动轨迹有很大不同。

为了研究液桥中热毛细对流的速度振荡特性，图1所示的模型中设置了6个监测点，坐标见表3。图10(a)和图10(b)所示分别为L/a=1(m=2)时液桥上半部的3 个监测点的速度变化和功率谱密度图(PSD)。图10(a)表明3 个监测点的速度呈现稳定的振荡，监测点a，b和c的速度ua

图9 L/a=0.72，ΔT=10 K和MaT=12 457时，1个周期内液桥俯视图Fig.9 Top view of liquid bridge within one period when L/a=0.72,ΔT=10 K and MaT=12 457

图10 L/a=1，ΔT=10 K和MaT=12 457时液桥上半部3个监测点的速度变化和这一阶段的功率谱密度Fig.10 Evolution of velocities at three monitoring points at the upper part and power spectral density(PSD)at this stage L/a=1,ΔT=10 K and MaT=12 457

表3 固定监测点的坐标Table 3 Coordinates at fixed monitoring points

图11(a)和图11(b)所示分别为L/a=1(m=2)时液桥下半部分的3个监测点的速度变化和功率谱密度图。图11(a)所示下部3 个监测点的速度也呈现稳定振荡，但是，监测点d，e，f的速度uf

2.4 热毛细对流中颗粒振荡机理分析

图12所示为近似的m倍对称结构，速度矢量可以用m对扰动涡描述，自由面的流动沿方位角从热点流向冷点。当马赫数Ma足够大时(Ma>Mac，Mac为不稳定起点的马赫数)，与热传导相比对流效应占主导地位，将产生“过冷或过热”并导致不稳定。由于自由表面存在温度梯度，局部冷点生成轴向扰动速度w′和方位角扰动速度u′。流动从热点流向冷点。根据连续性原理，方位角扰动速度u′形成了1对扰动涡。轴向扰动速度w′则从上盘带来热的流体来加热这个冷点。温度和速度扰动之间的相位差产生过热的热流，热流使局部冷点变为局部热点。新的局部热点也生成轴向扰动速度w′和方位角扰动速度u′，与冷点产生的扰动速度方向相反。局部热点产生的方位角扰动速度u′也形成了1对扰动涡，内部的冷流通过这个涡被带到表面来冷却这个热点。温度和速度扰动之间的相位差产生过冷的冷流，冷流使局部热点变为局部冷点。当Ma>Mac时，这个过程交替进行，冷点和热点之间的振荡一直保持甚至加剧。这种持续的振荡是3种高径比液桥中监测点(0.99a，0，0)处温度与速度曲线保持振荡的原因。

图11 L/a=1，ΔT=10 K和MaT=12 457时液桥下半部3个监测点的速度变化和这一阶段的功率谱密度Fig.11 Evolution of velocities at three monitoring points at the lower part and power spectral density(PSD)at this stage when L/a=1,ΔT=10 K and MaT=12 457

图12 扰动涡对模型Fig.12 Models of disturbance vortex pairs

若冷点(热点)继续保持或放大，则在方位角方向上，这个冷点(热点)发展成为m个冷热扰动，m取决于扰动涡的大小。以热点为例，3种高径比液桥中u4

图13 覆盖整个区域的扰动涡对Fig.13 Vortex pairs formed in the whole circumference region

图14(a)所示为单个颗粒的运动轨迹图，单个粒子的轨迹不是圆形的，也不是闭合的，而为花朵状的连续曲线。独立运动的单个颗粒在自由面附近向下运动，在冷盘附近改变运动方向，然后随回流向上运动，在热的圆盘与自由面交界处再次向下运动。颗粒的运动轨迹可以看作在r-z平面的二维翻转和在方位角方向运动的叠加。粒子运动轨迹不是闭合的回路，颗粒的运动轨迹和PAS形状结构有很大的不同。在自由面附近B点，颗粒速度最大，轴向温度梯度使颗粒在自由面附近获得较大的速度；在热毛细对流的中心A点和C点，颗粒速度最小。所以，粒子链在自由面附近是陡峭的，回流中更加水平。图14(b)所示为颗粒(该颗粒与图14(a)中的颗粒相同)速度变化曲线。颗粒速度振荡变化，在自由面附近速度变化更为剧烈，在回流中速度变化较小。颗粒速度曲线中微小变化的部分对应颗粒和其他颗粒碰撞，速度曲线产生不规则变化。许多粒子沿与图示类似的轨迹运动，这是PSA形成的直接原因。

图14 L/a=0.72，ΔT=10 K和MaT=12 457时PAS组成颗粒的轨迹俯视图与侧视图和颗粒的速度变化曲线Fig.14 Top view and side view for path lines of PASforming particle and velocity-time curve of the particle when L/a=0.72,ΔT=10 K and MaT=12 457

3 结论

1)热毛细对流的方位角波数m与液桥高径比L/a之间的关系为(L/a)·m≈2。三向耦合下的颗粒分布也显示出PAS。颗粒沿闭合环路积累，粒子链在空间做三维螺旋运动。

2)内部流场呈周期性变化，温度场绕z轴顺时针旋转，PAS 以相同的速度绕z轴旋转。水平面z=0 处温度场旋转1 周，液桥固定点处温度周期性变化m次。整个液桥内温度场与速度场之间存在强耦合，两者的时空特性相似。沿着液桥自由面，流体的速度从热端到冷端先增加后减小，速度的振荡强度则逐渐增强。

3)冷点与热点的交替转换，使振荡保持甚至加剧，最后产生了覆盖整个液桥的m对扰动涡，这是温度场呈m倍对称分布的原因。所有颗粒按照近似的轨迹协同运动，使PAS 随热毛细对流以相同的角速度旋转。