王苗苗，張傳成，任 浩，唐緒兵，丁守軍，鄒 勇，黃護林

(1.安徽工業(yè)大學微電子與數(shù)據(jù)科學學院，馬鞍山 243032；2.南京航空航天大學航天學院，南京 210016)

0 引 言

流體在表面張力梯度的驅(qū)動下會發(fā)生流動，這就是熱毛細對流。在浮區(qū)法生長晶體過程中，微重力環(huán)境下熱毛細對流是主要的對流方式。振蕩熱毛細對流是影響晶體生長質(zhì)量的主要原因之一[1]。在過去的十幾年中，研究人員在空間站和航天飛機內(nèi)進行了大量的微重力條件下液橋中熱毛細對流的研究。研究發(fā)現(xiàn)液橋中熱毛細對流結(jié)構(gòu)隨著Marangoni(Ma)數(shù)的變化而不斷演化[2-4]。Jin等[5]研究發(fā)現(xiàn)流速場沿自由面呈現(xiàn)兩個或四個主渦，這些對流渦在較高的Ma數(shù)下變得更強。當流體的Prandtl(Pr)數(shù)較高時，由表面張力梯度引起的熱毛細對流將從二維定常軸對稱對流過渡到振蕩對流，且振蕩對流可通過周向的三維波運動進行描述。Jayakrishnan等[6]對微重力條件下的高普朗特數(shù)5 cSt硅油(Pr=67)液橋中熱毛細對流的振蕩模態(tài)進行分解，分析結(jié)果表明，液橋自由表面的三維溫差所引起的擾動可以分解成各種動態(tài)模式，有助于揭示振蕩熱毛細流動背后的物理機制。Varas等[7]研究了相變材料在微重力下熔化過程中液橋的熱毛細效應(yīng)，結(jié)果表明隨著Ma數(shù)的增大，熱毛細對流從穩(wěn)定對流過渡到振蕩對流，并且在振蕩對流區(qū)，隨著Ma數(shù)增大先觀察到的是熱液行波(HTW)模式，然后觀察到振蕩駐波(OSW)模式。而當流體的Pr數(shù)較低時，其表面張力流隨著Ma數(shù)的增大將會先后出現(xiàn)兩次失穩(wěn)：第一次對流失穩(wěn)是二維定常軸對稱對流過渡到三維定常非軸對稱對流；另一次是向三維振蕩對流轉(zhuǎn)變。在天宮二號空間站，我國研究人員完成了液橋熱毛細對流不穩(wěn)定波模的空間實驗，研究了高徑比(Ar)和體積比(Vr)對大Pr數(shù)液橋的影響[8]。

對于固定幾何比率的圓柱形液橋，三維定常對流的臨界Reynolds(Re)數(shù)會隨著Pr數(shù)增大而增大。液橋表面張力流失穩(wěn)的臨界Re數(shù)受液橋幾何比率、液橋體積、流體物理性質(zhì)等多種因素的影響。Le等[9]通過采用三維瞬態(tài)數(shù)值模擬方法，研究了微重力下具有非等圓盤的低Pr數(shù)流體(Pr=0.01)半浮區(qū)液橋熱毛細對流不穩(wěn)定性的發(fā)生和一般特征，介紹了三維穩(wěn)態(tài)熱毛細對流和三維振蕩熱毛細對流的特點。結(jié)果表明,在不同圓盤半徑比的液橋中，方位流振幅之間的偏差很小，且振動頻率隨圓盤半徑比的增大而減小。

由于不同圓柱體液橋的高徑比存在差異，其內(nèi)部的熱毛細對流結(jié)構(gòu)也受到很大影響，同時，在晶體生長過程中高徑比明顯會限制晶體生長尺寸。Rybicki和Floryan[10]對液橋熱毛細對流進行了數(shù)值模擬，發(fā)現(xiàn)高徑比對液橋熱毛細流動有重要影響。Velten等[11]通過地面實驗研究了高徑比對液橋熱毛細對流的影響，結(jié)果表明當高徑比增加時，臨界Ma數(shù)逐漸減小。Chen和Hu[12]利用線性不穩(wěn)定性分析方法，研究了微重力環(huán)境下半浮區(qū)對流的毛細不穩(wěn)定性，分析了液橋體積和高徑比對臨界Ma數(shù)的影響。研究發(fā)現(xiàn)，在大Pr數(shù)和小高徑比的情況下，臨界Ma數(shù)與液橋體積的關(guān)系曲線中存在一個間隙區(qū)域，在該區(qū)域內(nèi)不會觀察到振蕩對流。Fan等[13]數(shù)值模擬了不同高徑比對熱毛細對流形成過程中的流動結(jié)構(gòu)影響，其中定義了強對流區(qū)(靠近界面)與弱對流區(qū)(位于液橋中心)的概念，得出強對流區(qū)的范圍在高徑比增大到1以后趨于穩(wěn)定。Kuhlmann和Nienhüser[14-15]報道了液橋自由界面的動態(tài)變形是由流體動壓、黏性應(yīng)力、液橋體積、表面張力和流體靜壓等因素共同決定的，而熱毛細對流的不穩(wěn)定性則取決于液橋的高徑比和自由表面的形狀。

熔體中產(chǎn)生的振蕩熱毛細對流對生長晶體的質(zhì)量有極大負面影響。對于具有高導(dǎo)電性的硅、鍺熔體，這類晶體在生長過程中大多選擇施加磁場來抑制有害流動。Dold等[16]實驗研究了旋轉(zhuǎn)磁場強度對摻雜硅浮區(qū)生長的影響，同時對流體流動進行了數(shù)值模擬。研究結(jié)果表明：在旋轉(zhuǎn)磁場作用下，時變?nèi)S流動轉(zhuǎn)變?yōu)闇瘦S對稱運動。Walker等[17]研究分析了旋轉(zhuǎn)磁場下浮動區(qū)熔體流動的線性穩(wěn)定性。結(jié)果表明，在沒有旋轉(zhuǎn)磁場的情況下，熔體流動會存在一系列不穩(wěn)定性，使晶體中出現(xiàn)生長條紋。由旋轉(zhuǎn)磁場的洛倫茲力引起的方位角速度可以降低這些不穩(wěn)定性，他們的結(jié)論與Dold等[16]的實驗結(jié)果一致。

綜上所述，影響液橋熱毛細對流不穩(wěn)定性的因素有很多，如高徑比、溫差、磁場構(gòu)形等。目前大多數(shù)文獻報道都是采用半浮區(qū)液橋模型進行的研究。半浮區(qū)液橋模型結(jié)構(gòu)簡單、邊界條件容易確定。研究者們通過固定半浮區(qū)液橋半徑，改變圓柱形液橋高度的方法對不同高徑比下熱毛細對流的穩(wěn)定性作了較為細致的描述，但關(guān)于振蕩熱毛細對流的產(chǎn)生及演化規(guī)律仍需進一步研究。本文通過固定半浮區(qū)液橋的高度，改變圓柱形液橋的半徑，探索熱毛細對流的穩(wěn)定性與高徑比之間的關(guān)系。

1 物理數(shù)學模型

計算所使用的物理模型為柱坐標系(r,θ,z)下的半浮區(qū)液橋模型，如圖1所示。不考慮零重力條件下自由面的形變，液橋自由面可近似為圓柱面。液橋懸浮在兩個半徑均為R的同軸圓盤之間，其高度為L。在位于z=-L/2和z=L/2處，視固液界面處液體速度為零，并使溫差ΔT(ΔT=Th-Tc)保持不變。

圖1 物理模型

σ=σ0+γTT

假設(shè)熔體是一種不可壓縮性的牛頓型黏性流體，作用于自由表面上的表面張力的大小可表示為以σ0為初值、以γT為比例系數(shù)的線性函數(shù)：

(1)

式中：γT=?σ/?T，熔體的其他物理特性不變。

在上述假設(shè)條件下，熔區(qū)內(nèi)N-S方程組可表示為：

(2)

(3)

(4)

邊界條件和初始條件設(shè)定如下：

當z=-L/2時，T=Tc

(5)

當z=L/2時，T=Th

(6)

(7)

(8)

式中：n是熔體表面的法向單位矢量，方向指向面外。硅熔體物性參數(shù)[18-19]及計算需要的其他參數(shù)如表1所示。

表1 硅熔體物性參數(shù)[18-19]和其他計算所需參數(shù)

2 計算方法

采用結(jié)構(gòu)化網(wǎng)格將計算區(qū)域進行劃分，為提高計算精度，在上下底面及熔體表面處進行局部加密，如圖2所示。當R=10 mm時，驗證網(wǎng)格獨立性，其獨立性驗證曲線如圖3所示，綜合考慮計算精度和計算時間，挑選40r×80θ×80z(對半徑R=5 mm和R=4 mm的液橋選取30r×80θ×80z)的網(wǎng)格開展數(shù)值計算。在計算中，設(shè)置時間步長為10-3s。采用有限體積方法對控制方程進行離散，非穩(wěn)態(tài)項采取二階隱式推進法。采用PISO算法的壓力速度耦合，對流項采用QUICK格式，擴散項采用二階中心差分格式。將該算法的結(jié)果同Levenstam等[20]的結(jié)果進行了比較，表2為在Pr=0.01，Re=3 500時，兩者的對比結(jié)果。從中可以看出吻合度很高。

圖2 網(wǎng)格劃分

圖3 網(wǎng)格獨立性驗證曲線

表2 本文和 Levenstam[20]的計算結(jié)果對比

3 結(jié)果與討論

在零重力環(huán)境下，由重力引起的自然對流消失，此時，由表面張力梯度引起的表面張力流起主導(dǎo)作用。熔體的表面張力與溫差密切相關(guān)。當熔體的上壁面和下壁面之間存在溫差時，在流體的自由表面上將會產(chǎn)生表面張力梯度，因此，驅(qū)動著自由表面附近的流體從熱端流向冷端，根據(jù)質(zhì)量守恒原理，從冷端流向熱端的回流在熔體內(nèi)部產(chǎn)生，從而形成一個對流循環(huán)。當高徑比Ar=0.5時，熔體對流表現(xiàn)出明顯的三維對流特征，此時z=0切面上的速度分布如圖4(a)～(c)所示。因為表面張力驅(qū)動的軸向速度并不相等，導(dǎo)致熔體回流速度在自由表面附近也具有方向角性，在熔體中心軸周圍，軸向速度又一次變?yōu)樨撝?。徑向速度在切面上形成兩對中心對稱分布的徑向波，切面上的周向速度形成的渦胞也成對出現(xiàn)，在熔體中心軸附近形成兩對中心對稱分布的渦胞，在外圍區(qū)域則形成與之對應(yīng)的兩對中心對稱分布的渦胞，形成完整的對流胞單元。數(shù)值結(jié)果表明，在切面中心區(qū)域，熔體沿周向從高溫區(qū)流向低溫區(qū)，從而形成了逆時針和順時針交替進行的四個對流渦胞單元；而在熔體表面附近區(qū)域，對流則從低溫區(qū)流向高溫區(qū)，這同表面張力流方向恰好相反，說明周向渦胞的形成與周向上的溫度差無關(guān)。

由于對流的影響，如圖4(d)所示，熔體z=0切面上的溫度分布表現(xiàn)為明顯的非軸對稱性：在切面外圍邊沿區(qū)域形成了一對低溫區(qū)和一對高溫區(qū)，由于高溫區(qū)附近的徑向溫度梯度遠大于低溫區(qū)域的徑向溫度梯度，這樣在熔體內(nèi)部也形成了一對與表面處的高溫區(qū)相對應(yīng)的低溫區(qū)。對比圖4(a)和圖4(d)不難發(fā)現(xiàn)，溫度等值線結(jié)構(gòu)和軸向速度等值線結(jié)構(gòu)具有很高的相似性，說明切面上溫度分布受軸向速度影響很大，在軸向速度出現(xiàn)極值的位置，溫度也相應(yīng)地出現(xiàn)極值。

在熔體表面附近和內(nèi)部區(qū)域各選取一個監(jiān)測點，監(jiān)測點的溫度和軸向速度隨時間的變化過程如圖4(e)、(f)所示。從圖4 (e)、(f)可以看出，流動時間t>125 s后，監(jiān)測點溫度和速度都不隨時間的變化而發(fā)生變化，熔體處于穩(wěn)態(tài)對流模式。

圖4 Ar=0.5時z=0切面上的速度、溫度分布和監(jiān)測點溫度、軸向速度隨時間的變化

當Ar=1時，熔體中的對流模式從三維穩(wěn)定流動轉(zhuǎn)變?yōu)槎囝l振蕩流動。首先來分析z=0切面上的速度分布。從圖5(a)可以看出，在某時刻，熔體的軸向速度在z=0截面上的分布具有方向性，并且在熔體表面，表面張力流動的速度高于Ar=0.5時的速度；與Ar=0.5時中心區(qū)域軸向速度為負值的情況不同，Ar=1時，切面上靠近中心軸的區(qū)域完全被回流占據(jù)。圖5(b)和圖5(c)顯示了z=0切面上相應(yīng)的徑向和周向速度分布。與Ar=0.5時的徑向速度分布相比，Ar=1時的徑向波結(jié)構(gòu)也是中心對稱分布，但軸向波的方位不同；隨著不穩(wěn)定性的增加，自由界面處的周向波被熔體中心形成的周向波擠壓，其占據(jù)面積明顯小于Ar=0.5時的面積。溫度分布如圖5(d)所示，熔體表面附近存在一對熱區(qū)和冷區(qū)，這與Ar=0.5時相同。不同之處在于，Ar=1時，因為回流占據(jù)了熔體的整個中心區(qū)域，所以最低溫度點出現(xiàn)在切面中央位置附近。

圖5 Ar=1某時刻z=0切面上速度、溫度分布和監(jiān)測點軸向速度、溫度隨時間變化與頻譜圖

同樣在液橋自由表面附近和內(nèi)部區(qū)域分別選取兩個監(jiān)測點，監(jiān)測點軸向速度和溫度隨時間變化情況如圖5(e)、(f)所示，其對應(yīng)的頻譜如圖5(g)、(h)所示。軸向速度和溫度隨時間呈周期性變化且振蕩頻率相同，這個結(jié)果表明，溫度振蕩是由熱毛細對流振蕩引起的。自由表面附近監(jiān)測點(4.9 mm,0,0)的振蕩主頻為f1=0.174 Hz，二次諧波頻率為f2=0.347 Hz，熔體內(nèi)部監(jiān)測點(2.5 mm,0,0)處振蕩主頻為f2，二次諧波頻率為f1，表明在熔體內(nèi)部熱毛細對流的不穩(wěn)定性被放大，頻率增加了一倍。此外，還發(fā)現(xiàn)了多個諧波頻率，諧波頻率之間均滿足倍頻關(guān)系fn=nf1。

溫差ΔT=5 K保持不變，半徑繼續(xù)減小至R=4 mm (Ar=1.25)。圖6顯示了不同時期監(jiān)測點軸向速度和溫度的變化以及相應(yīng)的頻譜圖。從圖6(a)可以看出，當流動時間t=25 s時，監(jiān)測點速度和溫度已經(jīng)呈現(xiàn)周期性振蕩，隨著時間的增加，速度平均值也在緩慢上升。當t>100 s時，速度和溫度的平均值均呈顯著上升趨勢，監(jiān)測點軸向速度分別在t=275 s和t=420 s時達到峰值，然后略有下降。當t>500 s時，軸向速度的平均值及振幅均不再隨時間變化。如圖6(b)顯示各監(jiān)測時段溫度振蕩的平均值和振幅均在緩慢下降；t>450 s時，溫度振蕩完全消失。這是因為軸向速度在等幅振蕩模式的振幅很小，不足以使溫度發(fā)生振蕩。

圖6 Ar=1.25時各個時間段內(nèi)監(jiān)測點軸向速度、溫度變化情況和頻譜圖

如上所述，在流動的初始階段(t<400 s)，軸向速度和溫度不僅隨時間振蕩，而且它們的平均值也發(fā)生變化。但根據(jù)頻譜分析發(fā)現(xiàn)，溫度和速度的振蕩頻率是相同的，不同監(jiān)測點的振蕩頻率也是相同的，都屬于單頻振蕩模式。圖6(c)顯示了t=25～100 s不同監(jiān)測點的頻譜圖，從圖中可以看出，自由表面和熔體內(nèi)部的溫度以及軸向速度的振蕩頻率f=0.398 Hz。當400 s500 s)的振蕩主頻率完全相同。

隨著Ar的變化，熔體內(nèi)的對流會呈現(xiàn)出多種形態(tài)。圖7(a)～(c)為t=500 s時，不同高徑比下子午面θ=0上的流線分布和溫度等值線。當Ar=0.5時，子午面θ=0上的流線左右對稱，兩邊的渦心在同一高度上。在熔體內(nèi)部，熔體水平流向中心區(qū)域，在中軸線附近交匯后一分為二，大部分熔體流向低溫面，在低溫面附近形成不完全的回流渦，剩余小部分熔體在交匯后向上底面流去。在對流影響下，溫度等值線在熔體中心軸和自由表面附近向下彎曲，而在熔體向上流動的區(qū)域，溫度等值線則向上凸出，因此等溫線的輪廓呈現(xiàn)字母“M”的形狀。當Ar=1時，左右兩邊對流渦胞的渦中心沿軸向稍微向下移動，渦中心高度相同，仍位于z=0橫切面附近。此時中心區(qū)域有熔體流過，左右對稱的對流結(jié)構(gòu)遭到破壞。在靠近低溫面的中心帶，有一對對流較弱的回流渦胞發(fā)育形成。子午面上的等溫線在自由表面附近向下彎曲，在中心附近區(qū)域則向上隆起。當Ar=1.25時，對應(yīng)R

圖7 不同高徑比下子午面(θ=0、π/2)上流線和溫度等值線

因為ΔT=5 K時的流動均表現(xiàn)為三維結(jié)構(gòu)，子午面θ=π/2上的對流模式和溫度分布相較于θ=0面上的有所不同，結(jié)果如圖7(d)～(f)所示。當Ar=0.5時，流線結(jié)構(gòu)圍繞中心軸對稱分布。位于中軸線上的高溫區(qū)域的熔體由左右分別流向自由表面，對流渦心在子午面兩底角位置。由于對流渦胞的核心區(qū)域只占子午面極小部分，所以溫度分布主要由熱擴散決定，在對流渦胞影響較大的兩底角附近，溫度等值線向下彎曲，其余部分則保持水平狀態(tài)。當Ar=1時，流線結(jié)構(gòu)與θ=0面上的流線結(jié)構(gòu)相似，但中心區(qū)域的流體的流向同θ=0面上的正好相反。溫度分布與θ=0面上的分布大致相同。當Ar=1.25時，整個熔體空間被左右兩個熱毛細對流渦胞占據(jù)，在左半部分表現(xiàn)逆時針流動而在右半部分表現(xiàn)順時針流動。此時溫度等值線在靠近中心軸處向上凸出的程度增加，在靠近自由表面處則向下彎曲程度更高。

4 結(jié) 論

對于熔體高度為5 mm的半浮區(qū)模型，采用數(shù)值模擬的方法，研究了零重力條件下溫差為5 K的熔體中熱毛細對流的動態(tài)特性，分析了高徑比變化對熱毛細對流的影響。闡明了熱毛細流動的轉(zhuǎn)捩過程及振蕩模式，加深了對浮區(qū)內(nèi)熱毛細流動特性的認識，模擬結(jié)果可為晶體生長工藝優(yōu)化提供指導(dǎo)和借鑒，具有重要科學意義和應(yīng)用價值。主要結(jié)論如下：

(1)當Ar=0.5時，熔體對流表現(xiàn)為三維穩(wěn)態(tài)結(jié)構(gòu)；當Ar=1時，熔體對流由穩(wěn)態(tài)對流變化為多頻振蕩對流，不同監(jiān)測點的物理量振蕩的主頻率各不相同，但同一監(jiān)測點的速度和溫度具有相同的頻譜結(jié)構(gòu)，主頻和各諧頻之間滿足倍頻關(guān)系fn=nf1；當Ar=1.25時，熔體開始時表現(xiàn)為單頻振蕩模式，其振蕩幅度隨時間逐漸衰減，經(jīng)過一段時間，溫度振蕩消失，最終，監(jiān)測點軸向速度演化為小幅振蕩模式。

(2)熔體內(nèi)溫度分布受對流影響很大，溫度等值線隨熔體流向發(fā)生彎曲。不同高徑比下，在接近自由表面的位置均形成兩對冷熱交替的區(qū)域，無論是穩(wěn)態(tài)流動還是振蕩對流，冷熱區(qū)域的方位角位置都幾乎不隨時間變化而改變。