直拉硅單晶生長過程中工藝參數對相變界面形態的影響?

張妮 劉丁 馮雪亮

1)(西安理工大學,晶體生長設備及系統集成國家地方聯合工程研究中心,西安 710048)2)(陜西省復雜系統控制與智能信息處理重點實驗室,西安 710048)(2018年2月7日收到;2018年7月28日收到修改稿)

為改善晶體相變界面形態,提高晶體品質,提出了一種融合浸入邊界法(immersed boundary method,IBM)和格子Boltzmann法(lattice Boltzmann method,LBM)的二維軸對稱浸入邊界熱格子Boltzmann模型來研究直拉法硅單晶生長中的相變問題.將相變界面視為浸沒邊界,用拉格朗日節點顯式追蹤相變界面;用LBM求解熔體中的流場和溫度分布;用有限差分法求解晶體中的溫度分布.實現了基于IB-LBM的動邊界晶體生長過程研究.得到了不同晶體生長工藝參數作用下的相變界面,并用相變界面位置偏差絕對值的均值和偏差的標準差來衡量界面的平坦度,得到平坦相變界面對應工藝參數的調整方法.研究表明,相變過程與晶體提拉速度、晶體旋轉參數和坩堝旋轉參數的相互作用有關,合理地配置晶體旋轉參數和坩堝旋轉參數的比值,能夠得到平坦的相變界面.

1 引 言

直拉法硅單晶生長中相變界面的形狀會直接影響晶體中的位錯密度大小以及剖面上電阻率的均勻性,因而一直受到學術界和產業界的關注[1].依據能量守恒原理可知相變界面形態由其上下兩側的溫度梯度決定,同時相變界面的動態變化也會影響熔體中的溫度場和速度場.相變與傳熱和流動之間的雙向耦合作用及其所具有的不確定性使晶體生長過程中相變問題的求解變得非常復雜和困難.因此,通過建立硅單晶生長模型,尋求復雜相變問題的求解方法,對研究硅單晶生長過程中不同工藝參數作用下的相變界面具有實際意義.

國內外學者對晶體生長過程中的相變、流動和傳熱問題進行了大量的研究,并取得了相應的研究成果[2?6].在流動與傳熱方面,格子Boltzmann法以其物理意義清晰、程序簡單易行等優點成為一種新興的研究方法并被廣泛使用.與傳統的流體動力學方法相比,格子Boltzmann法由于不需要離散Navier-Stokes(NS)方程從而避免了對流項離散帶來的數值不穩定問題[7].因此,針對流動與傳熱問題,本文采用格子Boltzmann方法求解NS方程.相變研究的核心主要在于相變界面的追蹤.在此方面,文獻[8]采用自適應網格法研究了自然對流作用下的相變過程.通過不斷建立貼體網格顯式追蹤相變界面,清晰地描述了相變的物理過程,但迭代過程時間成本高;文獻[9,10]用相場法分別研究了合金的凝固和鎵的融化過程.通過引入相場變量來區分液相和固相,即隱式追蹤相變界面.該方法用歐拉網格描述模擬區域,不需要建立貼體網格,從而節省了由于相變界面運動引起的網格重新劃分所需的時間,大大提高了計算效率.然而,文獻所描述的相變界面為糊狀區域,只有當網格劃分極為精細時,才能得到準確的相變界面;文獻[11,12]用熱焓法分別研究了強迫對流作用下水滴的凝固過程和自然對流作用下固體的融化過程.通過引入新變量——焓,來計算網格點的液相分數以判斷網格點的相態,不需要建立貼體網格,具有計算簡單且不需要顯式追蹤相變界面的優點,因此在相變研究領域得到廣泛應用.但也存在和相場法同樣的問題,即相變界面為糊狀區域.尤其當固體區域在外力作用下運動時,相變界面的糊狀區域可能會擴大,難以準確地獲得相變界面[13].文獻[14]用浸入邊界法求解流固耦合問題,該方法采用了兩套網格:歐拉網格和拉格朗日網格,歐拉網格只需要劃分一次,用歐拉網格描述模擬區域,用拉格朗日節點描述和顯式追蹤固-液界面.文獻[13]用浸入邊界法模擬了自然對流作用下移動固體的融化過程,結果表明,在解決具有移動邊界的相變問題時,該方法可以準確地跟蹤相變界面,不存在相場法和熱焓法中的糊狀邊界問題.在計算精度方面,文獻[13]的結果顯示浸入邊界法的計算精度高于熱焓法.在計算速度方面,文獻[8]的結果顯示,用自適應網格法模擬固-液相變過程時,網格生成所耗費的時間占總計算時間的20%—30%.浸入邊界法由于不需要建立貼體網格從而節省了大量的計算時間.因此,浸入邊界法以網格生成簡單、計算速度快且精度高的特點,在固-液相變領域展現出一定的優越性.目前,尚未見文獻報道如何使用該方法分析研究硅單晶生長這種不僅具有自然對流而且具有強迫對流的相變問題.本文采用浸入邊界法結合格子Boltzmann方法研究了晶體生長中的固-液相變和流動傳熱過程,并分析了工藝參數與相變界面形態之間的關系.

針對二維軸對稱硅單晶生長過程,本文提出了一種融合浸入邊界法和格子Boltzmann法的二維軸對稱浸入邊界熱格子Boltzmann模型,結合有限差分法,對伴有相變問題的硅單晶生長過程進行了研究.將相變界面視為浸入邊界,通過計算反饋力的形式修正流體邊界的速度與溫度.同時,根據界面能量守恒方程,建立相變界面運動模型,用拉格朗日節點顯式追蹤相變界面的位置;采用D2Q9模型構建熔體的密度演化方程、旋轉速度演化方程和溫度演化方程,用有限差分法求解晶體熱傳導方程,最后對晶體生長中的相變、流動和傳熱進行了深入的分析,得到直拉法硅單晶生長中工藝參數的選擇和調整方法.

2 控制方程

三維硅單晶生長簡化模型如圖1(a)所示,在穩態及準穩態情況下,可將三維晶體生長模型簡化為二維軸對稱晶體生長模型如圖1(b)所示.固-液相變界面處溫度恒等于相變溫度Tm,晶體半徑、坩堝半徑和坩堝高度之間滿足Rx=0.5Rc=0.5H的幾何尺寸關系.

圖1 直拉法硅單晶生長模型(a)及邊界條件(b)Fig.1.Model(a)and boundary conditions(b)of silicon single crystal in Czochralski method.

2.1 熔體對流與傳熱控制方程

為了用格子Boltzmann方法求解熔體中的流速分布和溫度分布,將圓柱坐標系下的晶體生長控制方程轉換到偽笛卡爾坐標系下,如下式所示:

其中(u,v)代表熔體流速在x和y方向的分量;w代表熔體繞y軸旋轉的速度;υ,ρ,α,g和β分別表示硅熔體的運動黏度、密度、熱擴散系數、重力加速度和熱膨脹系數.(2)和(3)式中fb=(fbx,fby)是相變界面產生的體積力項,具體形式見第3.2節;fg=gβ(T?T0)是由溫差引起的熱浮力.(5)式中gb是相變界面產生的能量力項,具體形式見第3.2節.熔體控制方程中fb和gb的引入,體現了相變界面對熔體流動和傳熱的作用.

2.2 相變控制方程

晶體生長中的相變界面即熔體與晶體的交界面,如(6)式所示的能量守恒方程,在單位時間內相變界面傳遞至晶體域的熱量等于通過熔體傳遞至相變界面的熱量和相變過程中釋放的潛熱之和:

λl(?T)l+ρsL(Us?Wpull)·n=λs(?T)s,(6)其中λs,λl分別為硅晶體熱傳導系數和硅熔體熱傳導系數;L為相變潛熱;Us為相變界面的運動速度;Wpull為晶體的提拉速度;(?T)s和(?T)l分別為晶體側和熔體側的溫度梯度.可以看出,相變界面的運動速度由晶體側的溫度梯度和熔體側的溫度梯度共同決定.

2.3 晶體熱傳導控制方程

晶體域只存在由溫度梯度產生的能量變化,滿足熱傳導方程,如下式所示:

3 基于浸入邊界法的軸對稱不可壓熱格子Boltzmann模型

3.1 軸對稱不可壓熱格子模型

格子Boltzmann方法通過求解介觀上描述熔體特征的演化方程來替代求解宏觀上復雜的NS方程.利用離散粒子的碰撞交換粒子節點之間的能量信息,利用離散粒子的遷移過程完成粒子節點在時間上的演化.本文采用數值穩定性和計算精度較高的D2Q9(2代表二維,9代表9個離散速度)模型[15]構建三個用來描述熔體的密度、旋轉速度和溫度的演化方程.

D2Q9模型中設置網格步長δx=1,時間步長δt=1,因此格子速度c=δx/δt≡1,9個方向的權值分別為:ω0=4/9,ω1,2,3,4=1/9,ω5,6,7,8=1/36.離散速度定義如下:

描述熔體密度分布的演化方程為

描述熔體旋轉速度的演化方程為

描述熔體溫度分布的演化方程為

這里τf,τh,τg代表無量綱松弛時間,滿足τf=3υ+0.5,τh=3υ+0.5,τg=3α+0.5;fi(x,t),hi(x,t)和gi(x,t)分別為t時刻格子點x=(x,y)處沿i方向的密度分布函數、旋轉速度分布函數和溫度分布函數;ψ(T,u,i)分別代表平衡態密度分布函數、平衡態旋轉速度分布函數和平衡態溫度分布函數,在D2Q9模型中,平衡態分布函數ψ(λ,u,i)=ωiλ[1+3ei·u/c2+4.5(ei·u)2/c4?1.5u2/c2].

流體的宏觀參數由下式確定:

3.2 相變界面作用力

在硅單晶生長過程中,相變界面不斷地吸收或釋放熱量,會對其周圍流體的速度和溫度分布產生影響.動量方程式(2)和(3)中fb項體現了相變界面對周圍流體產生的力作用,能量方程式(5)中gb項體現了相變界面對周圍流體產生的熱作用.體動量力fb和體能量力gb可根據沿浸沒邊界上的表面動量力Fs(sk,t)和表面能量力Gs(sk,t)求解,具體形式如下:

式中sk表示描述浸入邊界的第k個拉格朗日節點;X(sk,t+δt)表示第k個拉格朗日節點在t+δt時刻的歐拉坐標;是與狄拉克函數相關的平滑函數.狄拉克函數為

其中浸沒邊界的臨時速度Us和臨時溫度Ts可通過邊界周圍流體節點的速度以及溫度信息求解得到,如下式所示:

將(18)和(19)式代入(17)式中,便可得到相變界面對其周圍流體節點產生的體積力項fb和能量力項gb.最后依據(21)式對流體節點的信息進行修正:

其中u?(x,t+δt)和T?(x,t+δt)是未考慮相變作用時,由格子Boltzmann演化方程計算得到,如(16)式所示.可以看出用二維軸對稱浸入邊界熱格子Boltzmann模型求解伴有相變的流動與傳熱問題時,流體粒子的信息由標準格子Boltzmann演化方程和相變界面的作用力共同決定.

3.3 邊界處理

直拉法硅單晶生長存在四種邊界,如圖1(b)所示.本文采用非平衡態外推格式[16]處理坩堝壁,用對稱格式處理y軸和熔體自由表面.硅單晶生長界面隨時間推移不斷發生變化,屬于動態曲邊界問題,且邊界位置決定了晶體和熔體的幾何區域,其處理過程非常復雜.

從相變界面控制方程式(6)可以看出,相變界面的運動由其上下兩側的溫度梯度決定,而相變界面處溫度梯度的求解相當麻煩,但每一時間步長格子點釋放的熱量與界面能量力的積分相等[13],即

結合(6)式,可推導出相變界面的運動速度,

這種處理方式避免了相變界面處溫度梯度的求解.這里Cps是常壓下的晶體比熱;n=(nx,ny)是相變界面的運動法向量,如圖2所示,

因此,相變界面可依據上述理論進行動態演化.在演化過程中,一部分節點從流體狀態轉變為固體狀態,一部分節點由固體狀態轉變為流體狀態(如圖3中A點).對于前一種情況,只需將新產生的固體節點的分布函數置零;而對于新產生的流體節點,就必須進行分布函數的重構,本文采用鄰近插值法進行重構[17].另外,相變界面邊界處的流體粒子如圖3中B點,存在未知分布函數,本文采用Guo等[18]提出的非平衡態外推法重構曲邊界上流體粒子的未知分布函數.

圖2 相變界面法向量Fig.2.Normal vector of the interface.

圖3 曲邊界條件Fig.3.Curve boundary sketch.

3.4 無量綱參數

針對直拉法硅單晶生長過程,為了計算方便,本文將所有的有量綱參數進行無量綱化處理.涉及的無量綱參數有傅里葉數Fo,格拉斯霍夫數Gr,晶轉雷諾數Rex,堝轉雷諾數Rec,斯蒂芬數Ste和普朗特數Pr.它們的定義如下:

其中Fo數是用來描述非穩態熱傳導即分子擴散的無量綱數,可視作無量綱的時間;Gr數等于作用在流體上的浮力與黏性力之比,反映自然對流程度;Re數是判別黏性流體流動狀態的無量綱數;Ste數為相變貯能材料固相顯熱與相變潛熱之比;Pr數是由流體物性參數組成的無量綱數,反映流體物理性質對流動傳熱過程的影響;t為非穩態導熱過程所經歷的時間;?x,?c分別代表晶體旋轉角速度和坩堝旋轉角速度.在后面的計算中,用(Th?Tb),Rc和υ/Rc分別對溫度、長度和流體速度進行無量綱化處理;用速度的均方根誤差作為程序收斂條件,

4 算例驗證

為了驗證本文提出的融合了浸入邊界法和格子Boltzmann法的模型在伴有流動的相變問題中的正確性,本文以固-液相變基準測試問題——方腔熔化進行驗證.方腔內部初始化為固體狀態,設置左壁面恒為高溫Th,右壁面恒為低溫Tb,上下壁面均為絕熱狀態,無量綱參數Pr=0.02,Ra=2.5×104(Ra=Gr·Pr),Ste=0.01.模擬得到不同Fo數下方腔內溫度及流線分布,結果如圖4所示.

結果顯示,隨著Fo數的增大,流體區域不斷變大,熔化程度增強;同時,由于熱浮力的作用,方腔上部的熔化速度大于底部的熔化速度.將不同Fo數下相變界面提取出來,并與熱焓法[19]及自適應網格法[8]得到的結果進行對比,獲得了良好的一致性,如圖5所示,表明將本文模型應用在伴有流動的相變過程中是可行的.

圖4 不同Fo數下二維方腔熔化形態(溫度分布和流線分布) (a)Fo=4.0;(b)Fo=10.0;(c)Fo=20.0;(d)Fo=30.0Fig.4.Configuration of two-dimensional melting in square cavity(temperature distribution and streamlines):(a)Fo=4.0;(b)Fo=10.0;(c)Fo=20.0;(d)Fo=30.0.

圖5 不同相變處理方法下相變界面的比較Fig.5.Comparisons of the location of solid-liquid interface among different methods.

5 仿真研究與分析

首先,以二維軸對稱晶體生長模型為研究對象,驗證網格獨立性,確定最優網格數.模型邊界條件如圖1(b)所示,無量綱參數Pr=0.013,Gr=2.5×106,Ste=0.01.表1為穩態情況下不同網格數對應的熔體最大、最小流函數值,可以看出,當網格大小為90×90時,流函數最大值幾乎不會因為網格大小而發生變化,因此本文選用最優網格數90×90進行模擬.

表1 不同網格數下熔體流函數的最大值和最小值Table 1.Maximum and minimum values of flow function of melt under different grids.

5.1 提拉速度對相變界面的影響

在Gr=2.5×106且不考慮晶體和坩堝旋轉的情況下,模擬不同晶體提拉速度Wpull作用下晶體生長中的相變問題,為了清晰地對比溫度分布與流動結構,對所得結果關于y軸做鏡像顯示,結果如圖6所示.

由圖6可知,提拉速度對熔體流動結構和溫度分布并未產生明顯的作用,但相變界面發生了較大的變化,且由圖7可明顯地看到,當Wpull增大至0.0003時,相變界面凸向熔體的情況得到了很大程度的改善.

圖6 不同Wpull作用下晶體和熔體的溫度分布(左側)和流線分布(右側) (a)Wpull=0;(b)Wpull=0.0001;(c)Wpull=0.0002;(d)Wpull=0.0003Fig.6.Temperature distribution(left)and streamlines(right)of crystal and melt under different Wpullvalues:(a)Wpull=0;(b)Wpull=0.0001;(c)Wpull=0.0002;(d)Wpull=0.0003.

圖7 不同Wpull作用下的相變界面形態Fig.7.Phase transition interface with different Wpullvalues.

5.2 晶體旋轉對相變界面的影響

晶體旋轉是實際晶體生產中為實現熱場均勻性而必不可少的工藝手段之一. 本節以Wpull=0.0001,Gr=2.5×106為例,研究不同晶轉雷諾數下晶體生長中的相變問題,結果如圖8和圖9所示.

由圖8可知,在不考慮晶體旋轉作用,即Rex=0時,熔體內部由熱浮力驅動的逆時針渦流占據,熱熔體在熱浮力的作用下沿坩堝側壁向上運動,然后沿自由表面向晶體側運動,帶動冷熔體向坩堝底部運動形成回流.當晶體旋轉后,在熱浮力和晶體旋轉的共同作用下,晶體下方產生了兩個順時針的小渦流,且隨著晶體旋轉雷諾數的增大,小渦流的強度增強.當Rex=4000時,晶體下方的熔體幾乎全部由順時針的渦流控制,該渦流帶動坩堝底部的熔體向上運動,促進熱量向晶體側的傳遞,使得相變界面下方的溫度梯度增大,溫度分布變得更加均勻.

圖9顯示,隨著晶體旋轉參數的增大,相變界面形態從最初的凸向熔體變為部分凸向晶體,表明晶體旋轉對改善固-液界面附近熱場的均勻性起主要作用,這與文獻[20]的結果一致,達到了調節晶體轉速來改善相變界面形態的目的.當Rex=4000時,相變界面位置與自由表面(y=1)之間的偏差變小,相變界面形態呈M形狀.

圖8 不同Rex數下晶體和熔體的溫度分布(左側)和流線分布(右側) (a)Rex=0;(b)Rex=2000;(c)Rex=3000;(d)Rex=4000Fig.8.Temperature distribution(left)and streamlines(right)of crystal and melt in different Rexnumbers:(a)Rex=0;(b)Rex=2000;(c)Rex=3000;(d)Rex=4000.

圖9 不同Rex數下相變界面Fig.9.Phase transition interface in different Rex numbers.

5.3 晶體坩堝反向旋轉對相變界面的影響

本節以Wpull=0.0001,Gr=2.5×106,Rex=2000為例,研究坩堝與晶體反向旋轉時晶體生長中的相變問題.

圖10和圖11為坩堝與晶體反向旋轉作用下的結果.隨著坩堝旋轉參數的增大,雖然整個熔體中的溫度分布基本沒變,但相變界面兩側的局部溫度分布發生了非常大的變化,隨著坩堝旋轉參數的增大,晶體下方的溫度梯度增大.從圖11可看到,當Rec=?400時,相變界面形變量最小,相變界面最平坦;當Rec=?500時,相變界面呈W形狀,這是由于坩堝旋轉增大了逆時針渦流的強度,將晶體下方的渦流向晶體側擠壓,因此更多的熱量在該渦流的作用下被傳遞到相變界面處,相變界面向晶體側移動,最終該渦流處的相變界面凸向晶體.坩堝與晶體的反向旋轉使得原本凸向熔體的相變界面逐漸凸向晶體,與文獻[4]的結果一致.

圖10 不同Rec數下晶體和熔體的溫度分布(左側)和流線分布(右側) (a)Rec=?200;(b)Rec=?300;(c)Rec=?400;(d)Rec=?500Fig.10.Temperature distribution(left)and streamlines(right)of crystal and melt in different Recnumbers:(a)Rec=?200;(b)Rec=?300;(c)Rec=?400;(d)Rec=?500.

圖11 不同Rec數下的相變界面Fig.11.Phase transition interface in different Recnumbers.

5.4 晶體坩堝同向旋轉對相變界面的影響

本節以Wpull=0.0001,Gr=2.5×106,Rex=3000為例,研究坩堝與晶體同向旋轉時晶體生長中的相變問題.

圖12和圖13為坩堝與晶體同向旋轉作用下的結果.與反向旋轉作用時的結果相比(圖10),除了熱浮力和旋轉作用驅動的逆時針渦流外,在坩堝底部中心區域產生了一個順時針渦流,隨著坩堝旋轉參數的增大,坩堝底部中心區域的順時針渦流逐漸變大,并控制了對稱軸周圍的熔體流動.這是因為坩堝和晶體的同向旋轉,帶動了大部分熔體繞y軸做類似于剛體的旋轉運動,隨著坩堝旋轉作用的增強,熔體流動速度增大,而坩堝底部中心區域的熔體無法達到坩堝的轉速,從而在該區域形成順時針的渦胞.另外,從圖13可見,在Rec數增大的過程中,相變界面形態從最初的平坦變為凸向熔體,最后又出現凸向晶體的趨勢,其形態的演化與Rec數的變化不具有明確的規律性.

圖12 不同Rec數下晶體和熔體的溫度分布(左側)和流線分布(右側) (a)Rec=100;(b)Rec=200;(c)Rec=300;(d)Rec=400Fig.12.Temperature distribution(left)and streamlines(right)of crystal and melt in different Recnumbers:(a)Rec=100;(b)Rec=200;(c)Rec=300;(d)Rec=400.

實際晶體生產中,為了減少位錯的產生,相變界面應該越平坦越好,而不是凸向熔體或凸向晶體.因為平坦的相變界面表示晶體中徑向溫度梯度較小,從而保證晶體內部熱應力較小,避免由于熱應力過大而產生的位錯缺陷.以上的分析結果表明,晶體旋轉對改善相變界面附近熱場的均勻性起主要作用,坩堝旋轉作為一種輔助調節手段,可通過合理地配置晶體旋轉參數和坩堝旋轉參數實現對相變界面形態的改善.

圖13 不同Rec數下的相變界面Fig.13.Phase transition interface in different Recnumbers.

5.5 平坦相變界面下晶堝轉參數之間的函數關系

在模擬晶體旋轉和坩堝旋轉參數對相變界面形態的影響時,發現當坩堝與晶體反向旋轉時,在一組晶體旋轉參數下,總能找到一組相應的坩堝旋轉參數使得相變界面形態趨于平坦;當坩堝與晶體同向旋轉時,相變界面形態隨坩堝旋轉參數變化的規律性不明確,同一組晶體旋轉參數,可能存在兩組坩堝旋轉參數使得相變界面變得平坦.因此,考慮到工程應用中也大多采用反向旋轉技術,本節僅針對反向旋轉的情況討論.為了找到相變界面平坦時對應的晶體和坩堝旋轉參數,根據實際生產中工藝參數的調節范圍,本文對不同晶體旋轉和坩堝旋轉參數下的相變界面形態進行了大量的仿真,并進行了定量分析.

為了準確地衡量相變界面的平坦度,首先定義相變界面位置與熔體自由表面位置(y=1)之間的偏差為ek=X(sk,1)?1.用偏差絕對值的均值和偏差的標準差

來描述相變界面位置的偏差大小和波動程度,其中N是拉格朗日節點的個數,X(sk,1)表示拉格朗日節點sk的縱坐標.均值越小表明相變界面越靠近y=1,標準差越小表明相變界面的波動越小.因此,在實際晶體生長中,要求均值和標準差均越小越好.

表2為反向旋轉時不同Rex數和不同Rec數配置下相變界面偏差絕對值的均值和偏差的標準差.將每行結果中均值和標準差最小時對應的Rex數和Rec數提取出來,發現它們之間存在一定的關系,如圖14所示.為了定量地衡量Rex數和Rec數之間的關系,定義Re′x=Rex/1000,利用曲線擬合找到圖14中Rex數和Rec數之間的函數關系分別如下:

從表2和圖14可分析出:當Rex<2000時,隨著Rex數的增大,要保持相變界面平坦,Rec數也應該增大;當Rex>2000時,隨著Rex數的增大,要保持相變界面平坦,Rec數應該減小.這是因為,當Rex數較小時,晶體旋轉主導的順時針渦流較小,尚無法傳遞足夠的熱量到相變界面以實現相變界面向晶體側移動,相變界面完全凸向熔體.因此,需要較強的坩堝旋轉作用,將晶體旋轉主導的順時針渦流向晶體一側擠壓,改變晶體下方的溫度梯度,從而改變相變界面形態.當Rex數較大時,晶體旋轉主導的順時針渦流本身就較強,甚至會完全占據晶體下方,直接影響相變界面形態,此時,坩堝旋轉作為一種輔助的調節手段,僅需要較小的坩堝旋轉參數,對晶體下方順時針渦流的位置進行微調,就可以很大程度上調節相變界面形態.

另外,圖14顯示當晶轉雷諾數小于2000時,最小均值和最小標準差對應的晶轉與晶堝轉之比的關系基本相同,因此根據圖14中的曲線調節堝轉雷諾數就能保證相變界面偏差絕對值的均值最小且偏差的標準差最小.當晶轉雷諾數大于2000時,最小均值和最小偏差對應的晶堝轉參數之間的關系表現出明顯的差異.此時,需要結合對相變界面偏差以及波動程度的不同要求來選擇合適的晶堝轉參數.圖14可作為實際晶體生產中晶體旋轉參數和坩堝旋轉參數調節的參考.

圖14 最佳相變界面對應的晶堝轉關系Fig.14.Relationship of Rexand Recin optimal interface.

表2 反向旋轉時不同晶-堝轉作用下相變界面位置的標準差和相對誤差Table 2.Standard deviation and relative error of interface in different Rexnumbers and Recnumbers.

6 結 論

針對直拉法硅單晶生長中的相變、流動和傳熱問題,本文提出了一種融合浸入邊界法和格子Boltzmann法的二維軸對稱浸入邊界熱格子Boltzmann模型.用固-液相變基準測試問題驗證了所提模型的正確性.利用該模型分析了晶體生長中不同工藝參數對相變界面形態的影響,結果表明:

1)適當地提高晶體提拉速度能有效改善相變界面凸向熔體的問題;

2)在只有晶體旋轉作用時,雖然相變界面偏差較小但界面波動情況未得到改善;

3)在晶體-坩堝協同旋轉作用下,無論是反向旋轉還是同向旋轉,相變界面的偏差和波動均能得到有效調節,且在保持相變界面平坦的情況下,反向旋轉時晶體旋轉參數與晶堝旋轉參數比之間滿足一定的函數關系,對調整和優化工藝參數具有指導意義.