張妮 劉丁 馮雪亮

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

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

1 引 言

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

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

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

2 控制方程

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

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

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

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

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

2.2 相變控制方程

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

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

2.3 晶體熱傳導(dǎo)控制方程

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

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

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

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

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

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

描述熔體旋轉(zhuǎn)速度的演化方程為

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

這里τ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方向的密度分布函數(shù)、旋轉(zhuǎn)速度分布函數(shù)和溫度分布函數(shù);ψ(T,u,i)分別代表平衡態(tài)密度分布函數(shù)、平衡態(tài)旋轉(zhuǎn)速度分布函數(shù)和平衡態(tài)溫度分布函數(shù),在D2Q9模型中,平衡態(tài)分布函數(shù)ψ(λ,u,i)=ωiλ[1+3ei·u/c2+4.5(ei·u)2/c4?1.5u2/c2].

流體的宏觀參數(shù)由下式確定:

3.2 相變界面作用力

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

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

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

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

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

3.3 邊界處理

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

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

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

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

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

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

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

3.4 無量綱參數(shù)

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

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

4 算例驗證

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

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

圖4 不同F(xiàn)o數(shù)下二維方腔熔化形態(tài)(溫度分布和流線分布) (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 仿真研究與分析

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

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

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

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

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

圖6 不同Wpull作用下晶體和熔體的溫度分布(左側(cè))和流線分布(右側(cè)) (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作用下的相變界面形態(tài)Fig.7.Phase transition interface with different Wpullvalues.

5.2 晶體旋轉(zhuǎn)對相變界面的影響

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

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

圖9顯示,隨著晶體旋轉(zhuǎn)參數(shù)的增大,相變界面形態(tài)從最初的凸向熔體變?yōu)椴糠滞瓜蚓w,表明晶體旋轉(zhuǎn)對改善固-液界面附近熱場的均勻性起主要作用,這與文獻[20]的結(jié)果一致,達到了調(diào)節(jié)晶體轉(zhuǎn)速來改善相變界面形態(tài)的目的.當Rex=4000時,相變界面位置與自由表面(y=1)之間的偏差變小,相變界面形態(tài)呈M形狀.

圖8 不同Rex數(shù)下晶體和熔體的溫度分布(左側(cè))和流線分布(右側(cè)) (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數(shù)下相變界面Fig.9.Phase transition interface in different Rex numbers.

5.3 晶體坩堝反向旋轉(zhuǎn)對相變界面的影響

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

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

圖10 不同Rec數(shù)下晶體和熔體的溫度分布(左側(cè))和流線分布(右側(cè)) (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數(shù)下的相變界面Fig.11.Phase transition interface in different Recnumbers.

5.4 晶體坩堝同向旋轉(zhuǎn)對相變界面的影響

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

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

圖12 不同Rec數(shù)下晶體和熔體的溫度分布(左側(cè))和流線分布(右側(cè)) (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.

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

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

5.5 平坦相變界面下晶堝轉(zhuǎn)參數(shù)之間的函數(shù)關(guān)系

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

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

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

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

從表2和圖14可分析出:當Rex<2000時,隨著Rex數(shù)的增大,要保持相變界面平坦,Rec數(shù)也應(yīng)該增大;當Rex>2000時,隨著Rex數(shù)的增大,要保持相變界面平坦,Rec數(shù)應(yīng)該減小.這是因為,當Rex數(shù)較小時,晶體旋轉(zhuǎn)主導(dǎo)的順時針渦流較小,尚無法傳遞足夠的熱量到相變界面以實現(xiàn)相變界面向晶體側(cè)移動,相變界面完全凸向熔體.因此,需要較強的坩堝旋轉(zhuǎn)作用,將晶體旋轉(zhuǎn)主導(dǎo)的順時針渦流向晶體一側(cè)擠壓,改變晶體下方的溫度梯度,從而改變相變界面形態(tài).當Rex數(shù)較大時,晶體旋轉(zhuǎn)主導(dǎo)的順時針渦流本身就較強,甚至?xí)耆紦?jù)晶體下方,直接影響相變界面形態(tài),此時,坩堝旋轉(zhuǎn)作為一種輔助的調(diào)節(jié)手段,僅需要較小的坩堝旋轉(zhuǎn)參數(shù),對晶體下方順時針渦流的位置進行微調(diào),就可以很大程度上調(diào)節(jié)相變界面形態(tài).

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

圖14 最佳相變界面對應(yīng)的晶堝轉(zhuǎn)關(guān)系Fig.14.Relationship of Rexand Recin optimal interface.

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

6 結(jié) 論

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

1)適當?shù)靥岣呔w提拉速度能有效改善相變界面凸向熔體的問題;

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

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