微重力雙向溫差作用下Czochralski法硅熔體中的熱毛細對流＊

魏程星云，彭 嵐，張全壯

(重慶大學 動力工程學院，低品位能源利用技術及系統教育部重點實驗室，重慶 400044)

微重力雙向溫差作用下Czochralski法硅熔體中的熱毛細對流＊

魏程星云，彭 嵐，張全壯

(重慶大學 動力工程學院，低品位能源利用技術及系統教育部重點實驗室，重慶 400044)

采用數值模擬的方法研究了微重力條件下Czochralski法生長硅晶體過程中熔體熱毛細對流的基本特征，探討了水平和垂直溫度梯度的耦合對熔體流動的影響。熔體自由表面與外界輻射換熱，水平溫度梯度Marangoni(Ma)數選取(0～3 000)，底部熱流Q選取(1.39×10-2～1.76×10-2)。結果表明，當Q和Ma數均較小時，流動為穩態，液池內產生3個流胞，熔體流動由Q主導，減小Q或增大Ma數可使流動更穩定。當Ma數增大到一定值時，流動從穩態轉變為非穩態，流動的臨界Mac數隨Q的增大而顯著減小。流動失穩后，出現了新的流動轉變方式，Ma數為影響表面波動形式的關鍵因素，Q會改變熱流體波數，是晶體附近的熱流體波產生的決定因素。隨著Ma數和Q的不斷增強，自由表面最終形成彎曲條幅狀熱流體波。

熱毛細對流；Cz結構；雙向溫度梯度；數值模擬;表面張力

0 引 言

硅晶體因其優越的熱學、電學、光學等特性被廣泛應用于光伏工程、電子技術、航空航天技術等領域。目前，Czochralski生長法(Cz法)是制備硅晶體的最主要的方式之一。在Cz法制備硅單晶的過程中，有兩種不同形式的表面張力驅動硅熔體進行流動，一種是水平溫度梯度形成的表面張力所驅動的熱毛細對流；另一種是垂直溫度梯度所驅動的Marangoni對流。這兩種流動相互作用，對Cz法生成晶體的質量有著很大的影響。近些年來，許多學者對液池內的熱毛細對流進行了大量的研究[1-5]。Kyung-Woo Yi等[6]采用數值模擬的方法，發現熱毛細對流會影響Czochralski結構中的流體溫度和流動結構的不對稱性，并且認為Marangoni-Bénard不穩定性和Rayleigh-Bénard不穩定性是在不同表面張力系數下引起流動失穩的主要原因。Nepomnyashchy[7]最早使用數值模擬方法及線性穩定性理論對微重力條件下的矩形液池內的Bénard-Marangoni對流進行了研究，闡述了垂直和水平的溫度梯度的相對大小對液層內的流胞運動方向的影響。Jing等[8]發現Marangoni效應對大尺寸Cz結構內熔體自由表面的波動形式有著決定性的作用。Yao等[9-10]用數值模擬的方法對液橋內的熱毛細對流進行研究，結果表明磁場可以很好的抑制熱毛細對流，且普朗特數的大小對熱毛細對流起著重要的作用。Peng等[11]對方形池內的硅油進行了不同溫差下的實驗研究，實驗結果表明浮力對流可以抑制熱毛細對流。Y.Takagi等[12]發現坩堝旋轉和磁場相結合對熱毛細對流產生的的熱流體波有很好的抑制作用。T.Yamamoto等[13]發現通過改變液池的寬深比可以控制熔體內熱毛細對流的流動方向。Wang等[14]利用三維數值模擬研究了環形淺液池內垂直溫度梯度對硅油熱毛細對流的影響，結果表明增加底部熱流可以改變自由表面振蕩波波形，流動最不穩定的區域在環形池內壁附近。

以往國內外的大部分研究主要集中于Marangoni對流或者熱毛細對流單獨存在下的情況，很少有研究報道流體中兩種流動同時存在的情況，然而在晶體的實際生長過程中由于存在雙向溫差使得熔體中同時存在著兩種流動，稱為Marangoni-熱毛細對流；并且大部分學者研究的液池均為淺液池，而Cz法制備晶體的坩堝深度一般較深。因此，對雙向溫差下深液池中的熱毛細對流的研究顯得尤為重要。為了了解實際過程中雙向溫差耦合作用對熔體熱毛細對流的影響，本文采用有限差分法對Cz法生長硅晶體過程中的熱毛細對流進行了三維數值模擬，得出了不同溫差和底部熱流密度耦合的情況下的硅熔體流動的基本特征，重點分析了水平和豎直方向的溫度梯度耦合作用對熔體內部流動的影響。

1 數學物理模型

物理模型如圖1所示，坩堝內半徑rc=50 mm，晶體半徑rs=15 mm，液池深度h=8 mm，坩堝底部以恒定均勻熱流q加熱，自由邊界為不變形的平面且與環境進行輻射換熱，坩堝壁維持恒溫Th，結晶界面溫度為Tc=Tm，水平溫差ΔT=Th-Tc。

圖1 物理模型

為簡化計算，作出如下假定：(1) 熔體為不可壓縮牛頓流體；(2) 流速較低，流動為層流；(3) 所有固-液界面滿足無滑移條件；(4) 液體自由表面平整無變形；(5) 液體自由表面考慮熱毛細力的作用并且與外界進行輻射換熱。

在上述假定條件下，引入無量綱參數溫度

和無量綱熱流密度

·V=0

(1)

(2)

(3)

邊界條件為：

自由表面

(4)

(5)

熔體-晶體界面

Z=H,0≤R≤Rs,0≤θ≤2π

Vθ=VR=VZ=0,Θ=1

(6)

坩堝底部

Z=0,0≤R<1,0≤θ≤2π

(7)

坩堝側壁

0

Vθ=VR=VZ=0,Θ=Θh

(8)

初始條件(τ=0)

VR=VZ=Vθ=0,

(9)

其中

和

分別為普朗特數，毛細雷諾數和熱輻射數。定義Marangoni數為

式中

為表面張力溫度系數，v為動量擴散率，α為熱擴散率，σ為斯忒藩-玻耳茲曼常數，ε為發射率，λ為導熱系數，ρ為熔體密度。熔點溫度Tm=1 683 K下硅熔體的物性參數如表1所示。

表1Tm=1 683 K時硅熔體物性參數

Table 1 Physical properties of the silicon melt atTm=1 683 K

參數,符號單位數值表面張力溫度系數,γTN/(m·K)-7.0×10-5密度,ρkg/m32530動力粘度,μkg/(m·s)7.0×10-4動量擴散率,vm2/s2.767×10-7比熱,CpKJ/(kg·K)1000導熱系數,λW/(m·K)64熱擴散率,αm2/s2.53×10-5輻射率,ε0.318

利用有限容積法對基本方程進行離散，擴散項采用中心差分，對流項采用QUICK格式，壓力-速度修正采用SIMPLE方法，計算采用80r×40z×60θ的非均勻網格，時間步長取為1.19×10-3s，在每個時間步長內，如果控制容積上的最大殘差<10-8，則認為計算收斂。計算程序的正確性已被大量的計算結果所證實[15]，網格無關性驗證如表2所示，因此計算中選擇的網格為80r×40z×60θ。

表2 網格收斂性的檢驗(Ma=1 581，Q=1.39×10-2)

Table 2 Validation of grid convergence (Ma=1 581，Q=1.39×10-2)

網格(r×z×θ)振蕩頻率波數80×32×400.084480×40×600.0833480×60×800.08354

2 結果與討論

2.1 基本流動

當底部熱流密度和溫差較小時，液池內硅熔體的流動較弱為穩態軸對稱流動，豎直平面內的流動流型可用等流函數線表示(圖中正、負號分別代表逆時針和順時針方向)。無因次流函數定義如下

(10)

2.1.1Ma數對基本流的影響

當水平溫差較小時，流動為穩態。自由表面的熔體在熱毛細力的驅動下從坩堝壁向晶體流動，在熔體中產生逆時針流胞。圖2給出了Q=1.53×10-2時，不同Ma數下液池內R-Z平面的等流函數線和等溫線分布圖。當Ma數較小時，垂直溫差主導熔體內部流動，產生了占據主流區的順時針大流胞。晶體附近的熔體在水平溫差作用下形成一個逆時針小流胞，熔體內形成雙胞流。隨著Ma數增大，左側小流胞流動增強，右側大流胞流動被抑制，液池右上角出現微小的逆時針小流胞。熔體內流函數最大值隨Ma數增加而減小。由等溫線圖可知，隨Ma數的增加，熔體最高溫度和自由表面最高溫度增加。晶體下方等溫線保持水平，熔體幾乎無流動。

圖2 R-Z平面等流函數線和等溫線圖

圖3給出了Q=1.53×10-2時，自由表面徑向速度分布圖和R=0.4處的徑向速度在Z方向的分布圖。

圖3 自由表面徑向速度分布和R=0.4處的徑向速度沿Z向分布

Fig 3 Distribution of radial velocity along the free surface andZ

由圖3可見，自由表面上的徑向速度出現了3個極值?？拷w和坩堝壁的逆時針流胞使得自由表面形成了正極大值，而右側順時針流胞使得自由表面形成了負極大值。隨著Ma數的增加，自由表面正向流動增強而負向流動減弱，說明熔體中逆時針流胞流動增強而順時針流胞流動被抑制。從R=0.4處的徑向速度分布圖可見，隨著Ma數的增大徑向速度也增大。自由表面張力誘導的流動主要發生在0.1

2.1.2 熱流密度Q對基本流的影響

底部熱流加熱和自由表面的輻射散熱形成了貫穿液池的垂直方向的溫度梯度，誘發了自下而上順時針流動的右側大流胞。圖4給出了Ma=198時，液池內R-Z平面的等流函數線和等溫線分布圖。隨著Q增大，熔體流動強度增強，右上角小流胞消失，坩堝內剩下兩個反向旋轉的流胞。底部和自由表面的最高溫度逐漸坩堝壁移動。熔體中心等溫線發生劇烈變化而晶體下部等溫線保持水平。

當Ma=198時，垂直和水平溫差在靠近晶體側的驅動力方向一致，此處的徑向速度達到了最大值，如圖5所示。隨著Q增加，左側正向速度和右側負向速度增加，右側正向速度減小。這說明增大Q促進了左側流胞和右側大流胞的流動，抑制了右側小流胞的流動。R=0.4處熔體流動為逆時針流動，自由表面溫度梯度驅動的徑向內流集中在0.1

2.2 流動轉變

當底部熱流密度和自由表面溫度梯度不斷增大直到超過一定值時，流場中任何一個微小的擾動都會不斷地放大，使得流動狀態從穩態變為非穩態。采用逼近法確定了流動轉變的臨界Mac數。如圖6所示，增大Q可以使得流動轉變的臨界值減小，且較小的Q的改變可以造成較大臨界Mac數變化。

圖4 R-Z平面等流函數線和等溫線圖

圖5 自由表面徑向速度分布和R=0.4處的徑向速度沿Z向分布

Fig 5 Distribution of radial velocity along the free surface andZ

圖6 臨界Mac數隨Q的變化

2.3 三維振蕩流動

當Ma數超過臨界值后，熔體的流動轉變為三維非穩態流動。此時硅熔體的流動較為復雜，為了更直觀的了解熔體的周向速度的時相關特性，研究水平和垂直溫度梯度對流動的影響，定義自由表面溫度波動δT和自由表面速度波動δV為

(11)

(12)

2.3.1Ma數對振蕩流的影響

由圖7可見，隨著Ma數的增大，自由表面的最大溫度波動值不斷增加。圖8為Q=1.53×10-2，Ma=692時，一個周期內流動演變的過程。由圖8可見，自由表面上速度波動呈波數為5的單或雙層花瓣狀振蕩波，速度波動集中于晶體附近。振蕩由外向內傳遞，為單、雙層交替產生的結構。圖9為R=0.35處的自由表面速度波動隨時間變化的STD圖，圖中可見一簇明暗相間的直線且每條線從下往上由亮變暗，說明該點振蕩在徑向上有周期性的強弱變化，周向上沒有運動。

圖7 不同Q下的自由表面溫度波動最大值隨Ma數的變化

Fig 7 The max surface azimuthal temperature changes along the Ma at differentQ

圖8 Ma=642, Q=1.53×10-2時，一個周期內表面速度波動隨時間的變化

圖9Ma=642時R=0.35處表面速度波動的STD圖

Fig 9 The space-time diagram (STD) of surface azimuthal velocity along the circumference atR=0.35 forMa=642

圖10給出了Q=1.53×10-2時不同Ma數下自由表面的速度波動圖和STD圖。當Ma=791時，表面波動轉變為晶體附近振蕩的波數為5的花瓣狀熱流體波。隨著Ma數的進一步增大，振蕩從晶體附近擴散到了整個熔體表面。當Ma=988時，晶體附近熱流體波的流動區域被壓縮，在熱流體波外層出現了波數為7的振蕩波，由STD圖可見，外部振蕩為三維穩態結構。隨著Ma數增加到2 174，表面波動轉變為波數為5的熱流體波。從STD圖可以看出，此時熱流體波振幅呈周期性的變化。由此可見，自由表面振蕩隨Ma數的增加從晶體附近擴散到整個表面，隨著Ma數增大，流動越來越不穩定。

2.3.2 熱流密度Q對振蕩流的影響

在液池底部所加的熱流密度和表面散熱產生了垂直方向的溫度梯度，由圖7可見，自由表面的最大溫度波動值隨Q的增加不斷增加。

圖11給出了Q=1.62×10-2和Ma=692時，自由表面速度波動和STD圖?？梢钥闯?，隨著Q的增大，熱流體波波數和圖10(a)相比減少到4個，波動集中在晶體附近，STD圖由四對明暗相間的傾斜條紋組成。

圖12給出了Q=1.62×10-2,Ma=1 383時的表面速度波動和STD圖，在晶體附近存在波數為4的熱流體波。其流動區域較圖11明顯減小，在熱流體波外側出現了7條直幅狀振蕩波。

圖11Q=1.62×10-2,Ma=988,R=0.35表面速度波動圖和STD圖

Fig 11 Snapshots of surface azimuthal velocity and STD

由R=0.35處的STD圖可以看出，晶體附近波動為熱流體波，從R=0.75處的STD圖可知，外圍振蕩波為三維穩態結構。與圖10(b)對比可以發現，當流動轉變為外部三維穩態而內部為熱流體波的雙層結構時，隨著Q的增大，晶體附近的熱流體波流動區域增大，流動強度增強。由圖13可見，當Q較小而Ma數較大時，自由表面波動主要為三維穩態，而熱流體波則被壓縮到晶體附近區域幾乎消失，當Ma數更大時，熱流體波完全消失。

圖12 Q=1.62×10-2, Ma=1 383表面速度波動圖和STD圖

圖13Q=1.46×10-2,Ma=691時的表面速度波動圖和STD圖

Fig 13 Snapshots of surface azimuthal velocity fluctuation and STD

3 結 論

采用數值模擬的方法對雙向溫度梯度下Czochralski法生長硅晶體中的熔體熱毛細對流進行了研究。分析了垂直和水平溫差對熔體流動的影響，確定了流動轉變的臨界值。結果表明：

(1) 當Ma數和Q較小時，流動為穩態。隨Ma數的增加和Q的減小，熔體右上角出現一個微小的流胞，晶體下部熔體幾乎沒有流動。此時的熔體流動Q起主導作用，右上角小流胞和右側大流胞相互抑制，增大Ma數和減小Q可使得熔體整體流動減弱，流動更穩定；

(2) 隨著Ma數的增大，流動轉變為非穩態。臨界Mac數隨Q的增大而減小，且Q微小的改變能造成臨界Mac數的顯著改變。硅晶體生長過程中使用較小的Q可以延緩流動失穩現象。

(3) 當流動失穩轉變為非穩態后，Ma數為改變自由表面波形的關鍵因素。隨著Ma數的增加，自由表面波形由單、雙層花瓣互相轉換的振蕩波轉變為花瓣狀熱流體波，接著轉變為外層直幅波和內層熱流體波同時存在的結構并最終轉變為彎曲的條幅狀熱流體波，隨著Ma數的增大熔體自由表面波動不斷增強；

(4) 流動失穩后，隨著Q的增大，熱流體波波數減小，熔體自由表面波動不斷增強。當Q相對于Ma數較小時，靠近晶體側熱流體波消失，晶體附近流動較穩定；而增大Q，晶體附近速度波動增加，選用較大的Ma數和較小的Q有利于增加晶體附近流體流動的穩定性。

綜上所述，底部熱流Q對熔體內Marangoni-熱毛細對流的穩定性有著極為重要的影響，Q過大容易導致流動失穩并且會造成結晶界面附近劇烈振蕩，影響晶體質量。而選擇較大的Ma數可以抑制結晶界面附近的振蕩。因此在晶體生長過程中可以選用較小的Q搭配較大的Ma數來減弱結晶界面周圍的流體振蕩，提高晶體質量。

[1] Kumar V, Biswas G, Brenner G, et al. Effect of thermocapillary convection in an industrial Czochralski crucible: numerical simulation [J]. International Journal of Heat and Mass Transfer, 2003, 46(9):1641-1652.

[2] Anilkumar A V, Grugel R N, Bhowmick J, et al. Suppression of thermocapillary oscillations in sodium nitrate floating half-zones by high-frequency end-wall vibrations[J]. Journal of Crystal Growth, 2004, 276(s1-2):194-203.

[3] Teitel M, Schwabe D, Gelfgat A Y. Experimental and computational study of flow instabilities in a model of Czochralski growth [J]. Journal of Crystal Growth, 2008, 310(7-9):1343-1348.

[4] Hintz P, Schwabe D, Wilke H. Convection in a Czochralski crucible-part1: non-rotating crystal [J]. Journal of Crystal Growth, 2001, 222(1-2):343-355.

[5] Schwabe D. Buoyant-thermocapillary and pure thermocapillary convective instabilities in Czochralski systems [J]. Journal of Crystal Growth, 2002, 237-239:1849-1853.

[6] Yi K W, Kakimoto K, Eguchi M, et al. Spoke patterns on molten silicon in Czochralski system[J]. Journal of Crystal Growth, 1994, 144(1): 20-28.

[7] Nepomnyashchy A A, Simanovskii I B, Braverman L M. Stability of thermocapillary flows with inclined temperature gradient[J]. Journal of Fluid Mechanics, 2001, 442: 141-155.

[8] Jing C J, Tsukada T, Hozawa M, et al. Numerical studies of wave pattern in an oxide melt in the Czochralski crystal growth[J]. Journal of Crystal Growth, 2004, 265(265):505-517.

[9] Yao L, Zeng Z, Li X, et al. Effects of rotating magnetic fields on thermocapillary flow in a floating half-zone[J]. Journal of Crystal Growth, 2011, 316(1):177-184.

[10] Yao L, Zeng Z, Chen J, et al. Investigation of convection control under the non-uniform RMF in a liquid bridge[J]. Procedia Engineering, 2012, 31(4):659-664.

[11] Peng Z, Li D, Qi K. Transition to chaos in thermocapillary convection[J]. International Journal of Heat & Mass Transfer, 2013, 57(2):457-464.

[12] Takagi Y, Okano Y, Minakuchi H, et al. Combined effect of crucible rotation and magnetic field on hydrothermal wave[J]. Journal of Crystal Growth, 2014, 385(1):72-76.

[13] Yamamoto T, Takagi Y, Okano Y. Numerical investigation of the effect of free surface shape on the direction of thermocapillary flow in a thin circular pool[J]. Journal of Chemical Engineering of Japan, 2015, 48(6):407-417.

[14] Wang F, Peng L, Zhang Q Z, et al. Effect of the vertical heat transfer on the thermocapillary convection in a shallow annular pool[J]. Microgravity Science & Technology, 2015, 27(2):107-114.

[15] Li Y R, Imaishi N, Azami T, et al. Three-dimensional oscillatory flow in a thin annular pool of silicon melt[J]. Journal of Crystal Growth, 2004, 260(1):28-42.

Thermocapillary convection in Czochralski silicon growth by bidirectional temperature gradients under microgravity

WEI Chengxingyun，PENG Lan, ZHANG Quanzhuang

(Key Laboratory of Low-Energy Utilization Technologies and Systems of Ministry of Education, College of Power Engineering, Chongqing University, Chongqing 400044,China)

The bidirectional temperature gradients play an important part in the convection of crystal growth. However, most of the researches before only focused on unidirectional temperature. In this paper, under microgravity, a serious of numerical simulations was conducted to study the thermocapillary convection in Czochralski growth of silicon crystal. The coupling effects of vertical and horizontal temperature gradients were discussed. The radiation on the melt surface was considered, the range of horizontal and vertical temperature gradient wereMa=(0-3000) andQ=(1.39×10-2-1.76×10-2) respectively. Results show that flow is steady for smallQandMa. In this case, there are three flow cells in the melt, the melt flow is dominated byQ; the flow is more stable due to the decrease of theQand the increase of theMa. When the Ma exceeds a threshold value, an unsteady three dimensional flow is developed. With the increase of theQ, the critical Ma significantly decreased. When the flow is unsteady, a new way of the flow pattern transtion is found.Mais the key for the change of the surface fluctuation pattern. The number of the hydrothermal wave is changed byQ, and the generation of hydrothermal wave which is near the crystal is determined byQ. With the enhancement of theQandMa, the flow pattern enventually changes into curved spoke patterns.

thermocapillary convection; Cz configuation; bidirectional temperature gradients; numerical simulation; surface tension

1001-9731(2016)12-12090-07

國家自然科學基金資助項目(51276203)

2016-01-04

2016-03-15 通訊作者：彭 嵐，E-mail：penglan@cqu.edu.cn

魏程星云 (1991-)，男，武漢人，在讀碩士，師承彭嵐教授，從事工程熱物理、熱毛細對流相關研究。

TK124

A

10.3969/j.issn.1001-9731.2016.12.014