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

(1.東北大學(xué)材料電磁過程研究教育部重點實驗室，遼寧沈陽，110819；2.臨沂大學(xué)機械與車輛工程學(xué)院，山東臨沂，276000)

采用浮區(qū)法制備晶體過程中，在流體流動的驅(qū)動力為浮力、毛細力等。微重力環(huán)境下，浮力對流大大減弱，熱毛細對流成為主導(dǎo)對流形式。研究表明，振蕩熱毛細對流會導(dǎo)致晶體表面和內(nèi)部產(chǎn)生條紋，影響晶體的質(zhì)量。因此，在微重力環(huán)境下，半導(dǎo)體熔體內(nèi)熱毛細對流及流動結(jié)構(gòu)研究已成為空間材料生長的重要課題，對空間浮區(qū)法制備高質(zhì)量晶體具有重要意義。

液橋是浮區(qū)法制備晶體簡化模型。1979年，CHUN等[1]利用“光切割”技術(shù)，開展了半浮區(qū)液橋熱毛細對流實驗，發(fā)現(xiàn)了振蕩熱毛細對流。之后，針對液橋模型，國內(nèi)外開展了線性穩(wěn)定性分析[2-3]、三維數(shù)值模擬[4-6]和晶體生長實驗[7-9]。ZENG等[10]數(shù)值研究了三維振蕩熱毛細對流，結(jié)果表明熱毛細對流存在脈動振蕩和旋轉(zhuǎn)振蕩兩種不同的振蕩模式，溫度場隨時間周期性變化。YASUHIRO 等[11]通過數(shù)值模擬再現(xiàn)了實驗觀察到的不同頻率的超臨界熔體自由表面溫度振蕩，結(jié)果表明當溫差達到一定的臨界值時，表面溫度振蕩將向三維穩(wěn)定周期振蕩對流轉(zhuǎn)變。當溫差為更高的臨界值時，表面溫度振蕩將會發(fā)生向三維振蕩對流的第二次轉(zhuǎn)變。

另外，浮區(qū)法制備晶體過程中浮區(qū)熔體中存在雜質(zhì)，雜質(zhì)顆粒具有動態(tài)積累特性，影響制備晶體質(zhì)量。研究表明，在一定高徑比條件下，稀密度的小顆粒在液橋中聚集，并具有時間依賴性。這些顆粒繞熱毛細對流的軸旋轉(zhuǎn)呈變形的螺旋，形成顆粒動態(tài)積累結(jié)構(gòu)(PAS)[12]。SCHWABE 等[13]發(fā)現(xiàn)了PAS。UENO 等[14]也觀察到了類似的結(jié)果，同時得到閉合螺旋繞組數(shù)等于方位角波數(shù)。SCHWABE等[15]通過改變等密度下的顆粒直徑和等粒徑下顆粒密度與流體密度的比值，測量PAS 結(jié)構(gòu)的形成時間，解釋了PAS 的形成機制，即顆粒與流體的密度差影響PAS 的形成速度。MULDOON等[16]的研究表明，密度差和顆粒-自由表面相互作用共同影響PSA的形成。

目前，對PAS 的研究大多基于單向耦合，考慮三向耦合(同時考慮流體對顆粒的作用、顆粒對流體的反作用、顆粒之間的相互作用)的PAS 研究缺乏，因此，對PAS的形成機制仍不清楚。另外，國內(nèi)外關(guān)于熱毛細對流的研究主要集中在單一流體熱毛細對流研究，對夾雜顆粒的固液兩相流體熱毛細對流研究極其有限，對熱毛細對流的周期性變化與熱毛細對流中顆粒聚集規(guī)律之間相互關(guān)系的研究很少。采用浮區(qū)法制備單晶硅過程中熔區(qū)內(nèi)不可避免含有雜質(zhì)顆粒，針對這種夾雜顆粒的固液兩相流體熱毛細對流研究更具有理論意義與實際應(yīng)用價值。

本文作者基于一種三向耦合的CFD-DEM 方法，同時考慮流體對顆粒的作用、顆粒對流體的反作用、顆粒之間的相互作用，采用CFD 方法計算流體的流動行為，利用DEM方法處理顆粒運動和顆粒間的碰撞，對夾雜顆粒流體熱毛細對流進行了數(shù)值模擬，探究不同高徑比下的熱毛細對流和PAS，分析其形成原因。

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

本研究采用的物理模型如圖1所示。2 個半徑相同的同心圓盤之間懸浮的是半浮區(qū)液橋，液橋的高度為L，半徑為a，液橋高徑比為L/a。上下圓盤之間存在溫差，溫度梯度的方向與重力的方向相反。

圖1 液橋模型Fig.1 Liquid bridge model

考慮流體對顆粒的作用、顆粒對流體的反作用、顆粒之間的相互作用，三向耦合下液相的控制方程[17-23]如下。

連續(xù)性方程為

動量方程為

能量方程為

根據(jù)牛頓第二定律，單個顆粒的轉(zhuǎn)動和平動方程分別為：

單個顆粒的能量守恒方程為

系統(tǒng)的邊界條件為：下部圓盤(z=-L/2)處溫度為T=Tc，上部圓盤(z=L/2)處溫度為T=Th，上下圓盤絕熱無滑移。自由面的徑向溫度梯度?T/?r=0。在自由面(R=a)上，考慮熱毛細效應(yīng)，采用下列方程：

式中：r，θ和z為坐標；ur，uθ和uz分別為r，θ和z這3個方向上的速度分量；ρf，uf，p，Tf，cf和Г分別為流體的密度、速度、壓力、溫度、比熱容和熱擴散系數(shù)；ε，t，τ，g和σ分別為孔隙率、時間、氣相應(yīng)力張量、重力加速度和表面張力；mi，vi，ωi，Ii，Ti和cp,i分別為顆粒i的質(zhì)量、線速度、角速度、轉(zhuǎn)動慣量、溫度和比熱容；ffp為顆粒-流體作用力；Qf,i為顆粒與流體間的熱流量；ff,i為顆粒-流體作用力；fe,ij和fd,ij分別為彈性力矢量與黏性力矢量；Tn,ij，Tt,ij和Tr,ij分別為法向力矩、切向力矩和滾動摩擦力矩；Qi,j和Qi,f分別為顆粒間的換熱量和顆粒與流體的換熱量；MaT為熱Marangoni 數(shù)，MaT=rTTL/(μυ)；Pr為普朗特數(shù)，Pr=υ/α；ΔT為上下板間的溫差，ΔT=Th-Tc。

對3組不同高徑比下的含夾雜物的兩相流體熱毛細對流進行模擬，所用的流體介質(zhì)為0.65 號硅油，Pr=6.7。采用有限體積法對控制方程進行離散，通過PISO 算法求解控制方程。采用CFDDEM下的Eulerian-Lagrangian模型，模擬含夾雜物的兩相流。顆粒間的作用模型和顆粒與自由面的作用模型均為Ranz-Marshall 模型，模擬所用的顆粒為直徑70 μm的鋁球。在計算的起始時刻，液橋中均勻注入所有顆粒。計算所需的參數(shù)見表1。3種高徑比下均在液橋的(0.99a，0，0)處設(shè)置監(jiān)測點。

表1 流體物理性質(zhì)Table 1 Fluid physical property

為了驗證網(wǎng)格的相關(guān)性，在高徑比為1，ΔT=10 K 和MaT=12 457 時，分別采用34r×36z×60θ，34r×38z×80θ和36r×40z×100θ這3 種網(wǎng)格進行計算，監(jiān)測點(0.99a，0，0)處計算結(jié)果如表2所示。當網(wǎng)格數(shù)為34r×38z×80θ和36r×40z×100θ時，監(jiān)測點(0.99a，0，0)處溫度和速度的相對誤差均小于1%，因此，選用網(wǎng)格數(shù)更少的網(wǎng)格34r×38z×80θ作為數(shù)值計算網(wǎng)格。

表2 網(wǎng)格無關(guān)性檢驗Table 2 Grid dependence test

2 結(jié)果與討論

2.1 模擬結(jié)果驗證

首先，將本研究的結(jié)果與SCHWABE等[12]的實驗結(jié)果以及MELNIKOV 等[24]的數(shù)值模擬結(jié)果進行比較驗證。當ΔT=8.5 K，Pr=13.5，L/a=1 時，在計算的起始時刻往液橋中均勻注入顆粒并施加溫差，一段時間后液橋中出現(xiàn)穩(wěn)定的溫度和速度振蕩，再經(jīng)歷一段時間后，顆粒呈現(xiàn)出動態(tài)積累，圖2(a)所示為本研究模擬的液橋注入顆粒250 s 后的俯視圖。顆粒分布與圖2(b)中SCHWABE等[12]的實驗結(jié)果和圖2(c)中MELNIKOV 等[24]的模擬結(jié)果非常相似，也呈不對稱的單環(huán)，PAS繞z軸逆時針旋轉(zhuǎn)。圖3所示為方位角波數(shù)模擬結(jié)果與SCHWABE 等[12]的實驗結(jié)果對比。當ΔT=10 K，Pr=8 時，改變高徑比，數(shù)值模擬測得對應(yīng)的方位角波數(shù)。由圖3可以看出：方位角波數(shù)數(shù)值模擬結(jié)果與實驗結(jié)果相吻合，隨著高徑比增大，兩者變化趨勢相同，驗證了本研究建立的模型與計算方法的有效性。

圖2 顆粒動態(tài)積累結(jié)構(gòu)對比Fig.2 Comparison of dynamic particle accumulation structure

圖3 方位角波數(shù)模擬結(jié)果與實驗結(jié)果對比Fig.3 Comparison of azimuthal wave number between simulation results and experimental results

2.2 熱毛細對流中顆粒動態(tài)聚集特性

圖4所示為液橋模型中水平面(z=0)的溫度分布。從圖4可見：水平面中存在著2，3和4倍的對稱結(jié)構(gòu)，對稱結(jié)構(gòu)與方位角波數(shù)m對應(yīng)；當上下圓盤的溫差ΔT=10 K，高徑比為0.50，0.72 和1時，方位角波數(shù)m分別為4，3和2；隨著高徑比增大，方位角波數(shù)減少，水平面溫度分布中低溫冷區(qū)的數(shù)目減少，同時可以得到(L/a)·m≈2。圖5所示為不同高徑比液橋顆粒聚集結(jié)構(gòu)的俯視圖。從圖5可見：中心出現(xiàn)顆粒較少區(qū)，形狀分別為梭形(m=2)、三角形(m=3)、正方形(m=4)；在這個中心區(qū)域，顆粒速度較低。這個區(qū)域繞z軸以恒定的角速度旋轉(zhuǎn)，旋轉(zhuǎn)過程中形狀保持不變。這種粒子較少的多邊形區(qū)域的形狀也與液橋的高徑比有關(guān)。粒子鏈組成了多邊形的邊，多邊形的尖端連接液橋的頂部和底部m次。中心區(qū)域的形狀對應(yīng)著溫度場中的m倍對稱結(jié)構(gòu)。圖6所示為不同高徑比液橋中顆粒分布的鳥瞰圖。從圖6可見：顆粒聚集結(jié)構(gòu)并不局限于1個平面內(nèi)，而是分布在整個三維空間；當PAS 形成后，自由面和液橋體積的很大部分沒有顆粒。顆粒沿閉合的環(huán)路積累，該環(huán)路在液橋中作三維螺旋運動，螺旋的形狀不變；這個粒子鏈從熱端開始，沿自由面向下，再隨回流上升，直到再次接觸到自由面；帶著粒子的回流并沒有穿過液橋中心，粒子鏈在自由面附近是陡峭的，在回流中則更加水平。粒子鏈不同的陡峭程度反映著顆粒不同的運動速度，即顆粒在自由表面附近速度更大，在回流中速度較小。

圖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時不同高徑比的液橋顆粒聚集結(jié)構(gòu)的俯視圖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 熱毛細對流旋轉(zhuǎn)振蕩特性

圖7所示為L/a=0.72(m=3)時，1個周期內(nèi)水平面z=0處溫度云圖和點(0.99a，0，0)處的溫度變化曲線。從圖7可見：水平面z=0 處溫度云圖繞z軸旋轉(zhuǎn)，旋轉(zhuǎn)的周期Tmin=11.25 s。圖7中的曲線為m=3時，監(jiān)測點(0.99a，0，0)處溫度變化曲線，振蕩周期T′min=3.75 s。可以得到Tmin=3·T′min。當方位角波數(shù)m=2，3，4時，監(jiān)測點處溫度分別為T2，T3和T4，速度分別為u2，u3和u4。圖8(a)所示為3種高徑比液橋中監(jiān)測點(0.99a，0，0)處溫度變化曲線。可見3個監(jiān)測點的溫度出現(xiàn)了穩(wěn)定的振蕩，監(jiān)測點溫度T2

圖7 L/a=0.72，ΔT=10 K和MaT=12 457時，1個周期內(nèi)的溫度云圖和點(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(方位角波數(shù)m=3)時，1個周期內(nèi)液橋中顆粒聚集結(jié)構(gòu)的俯視圖。在這個周期內(nèi)，PAS繞z軸逆時針旋轉(zhuǎn)。為清楚地展示PAS的旋轉(zhuǎn)，三角形的1個頂點用圓點標記。旋轉(zhuǎn)1周的時間為11.25 s，與液橋溫度場的變化周期相同，所以，液橋的周期性旋轉(zhuǎn)振蕩與PAS 的轉(zhuǎn)動相對應(yīng)，PAS的旋轉(zhuǎn)周期與液橋溫度場的旋轉(zhuǎn)周期相同，但兩者的旋轉(zhuǎn)方向相反。所有顆粒協(xié)同運動，造成PAS 隨熱毛細對流以相同的角速度旋轉(zhuǎn)，這種運動趨勢僅適用于整個PAS，與單個顆粒的運動軌跡有很大不同。

為了研究液橋中熱毛細對流的速度振蕩特性，圖1所示的模型中設(shè)置了6個監(jiān)測點，坐標見表3。圖10(a)和圖10(b)所示分別為L/a=1(m=2)時液橋上半部的3 個監(jiān)測點的速度變化和功率譜密度圖(PSD)。圖10(a)表明3 個監(jiān)測點的速度呈現(xiàn)穩(wěn)定的振蕩，監(jiān)測點a，b和c的速度ua

圖9 L/a=0.72，ΔT=10 K和MaT=12 457時，1個周期內(nèi)液橋俯視圖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個監(jiān)測點的速度變化和這一階段的功率譜密度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 固定監(jiān)測點的坐標Table 3 Coordinates at fixed monitoring points

圖11(a)和圖11(b)所示分別為L/a=1(m=2)時液橋下半部分的3個監(jiān)測點的速度變化和功率譜密度圖。圖11(a)所示下部3 個監(jiān)測點的速度也呈現(xiàn)穩(wěn)定振蕩，但是，監(jiān)測點d，e，f的速度uf

2.4 熱毛細對流中顆粒振蕩機理分析

圖12所示為近似的m倍對稱結(jié)構(gòu)，速度矢量可以用m對擾動渦描述，自由面的流動沿方位角從熱點流向冷點。當馬赫數(shù)Ma足夠大時(Ma>Mac，Mac為不穩(wěn)定起點的馬赫數(shù))，與熱傳導(dǎo)相比對流效應(yīng)占主導(dǎo)地位，將產(chǎn)生“過冷或過熱”并導(dǎo)致不穩(wěn)定。由于自由表面存在溫度梯度，局部冷點生成軸向擾動速度w′和方位角擾動速度u′。流動從熱點流向冷點。根據(jù)連續(xù)性原理，方位角擾動速度u′形成了1對擾動渦。軸向擾動速度w′則從上盤帶來熱的流體來加熱這個冷點。溫度和速度擾動之間的相位差產(chǎn)生過熱的熱流，熱流使局部冷點變?yōu)榫植繜狳c。新的局部熱點也生成軸向擾動速度w′和方位角擾動速度u′，與冷點產(chǎn)生的擾動速度方向相反。局部熱點產(chǎn)生的方位角擾動速度u′也形成了1對擾動渦，內(nèi)部的冷流通過這個渦被帶到表面來冷卻這個熱點。溫度和速度擾動之間的相位差產(chǎn)生過冷的冷流，冷流使局部熱點變?yōu)榫植坷潼c。當Ma>Mac時，這個過程交替進行，冷點和熱點之間的振蕩一直保持甚至加劇。這種持續(xù)的振蕩是3種高徑比液橋中監(jiān)測點(0.99a，0，0)處溫度與速度曲線保持振蕩的原因。

圖11 L/a=1，ΔT=10 K和MaT=12 457時液橋下半部3個監(jiān)測點的速度變化和這一階段的功率譜密度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

若冷點(熱點)繼續(xù)保持或放大，則在方位角方向上，這個冷點(熱點)發(fā)展成為m個冷熱擾動，m取決于擾動渦的大小。以熱點為例，3種高徑比液橋中u4

圖13 覆蓋整個區(qū)域的擾動渦對Fig.13 Vortex pairs formed in the whole circumference region

圖14(a)所示為單個顆粒的運動軌跡圖，單個粒子的軌跡不是圓形的，也不是閉合的，而為花朵狀的連續(xù)曲線。獨立運動的單個顆粒在自由面附近向下運動，在冷盤附近改變運動方向，然后隨回流向上運動，在熱的圓盤與自由面交界處再次向下運動。顆粒的運動軌跡可以看作在r-z平面的二維翻轉(zhuǎn)和在方位角方向運動的疊加。粒子運動軌跡不是閉合的回路，顆粒的運動軌跡和PAS形狀結(jié)構(gòu)有很大的不同。在自由面附近B點，顆粒速度最大，軸向溫度梯度使顆粒在自由面附近獲得較大的速度；在熱毛細對流的中心A點和C點，顆粒速度最小。所以，粒子鏈在自由面附近是陡峭的，回流中更加水平。圖14(b)所示為顆粒(該顆粒與圖14(a)中的顆粒相同)速度變化曲線。顆粒速度振蕩變化，在自由面附近速度變化更為劇烈，在回流中速度變化較小。顆粒速度曲線中微小變化的部分對應(yīng)顆粒和其他顆粒碰撞，速度曲線產(chǎn)生不規(guī)則變化。許多粒子沿與圖示類似的軌跡運動，這是PSA形成的直接原因。

圖14 L/a=0.72，ΔT=10 K和MaT=12 457時PAS組成顆粒的軌跡俯視圖與側(cè)視圖和顆粒的速度變化曲線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 結(jié)論

1)熱毛細對流的方位角波數(shù)m與液橋高徑比L/a之間的關(guān)系為(L/a)·m≈2。三向耦合下的顆粒分布也顯示出PAS。顆粒沿閉合環(huán)路積累，粒子鏈在空間做三維螺旋運動。

2)內(nèi)部流場呈周期性變化，溫度場繞z軸順時針旋轉(zhuǎn)，PAS 以相同的速度繞z軸旋轉(zhuǎn)。水平面z=0 處溫度場旋轉(zhuǎn)1 周，液橋固定點處溫度周期性變化m次。整個液橋內(nèi)溫度場與速度場之間存在強耦合，兩者的時空特性相似。沿著液橋自由面，流體的速度從熱端到冷端先增加后減小，速度的振蕩強度則逐漸增強。

3)冷點與熱點的交替轉(zhuǎn)換，使振蕩保持甚至加劇，最后產(chǎn)生了覆蓋整個液橋的m對擾動渦，這是溫度場呈m倍對稱分布的原因。所有顆粒按照近似的軌跡協(xié)同運動，使PAS 隨熱毛細對流以相同的角速度旋轉(zhuǎn)。