張妮 劉丁 馮雪亮

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

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

1 引 言

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

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

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

2 控制方程

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

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

2.1 熔體對(duì)流與傳熱控制方程

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

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

2.2 相變控制方程

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

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

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

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

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

3.1 軸對(duì)稱不可壓熱格子模型

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

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

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

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

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

這里τf,τh,τg代表無(wú)量綱松弛時(shí)間,滿足τf=3υ+0.5,τh=3υ+0.5,τg=3α+0.5;fi(x,t),hi(x,t)和gi(x,t)分別為t時(shí)刻格子點(diǎn)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 相變界面作用力

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

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

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

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

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

3.3 邊界處理

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

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

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

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

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

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

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

3.4 無(wú)量綱參數(shù)

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

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

4 算例驗(yàn)證

為了驗(yàn)證本文提出的融合了浸入邊界法和格子Boltzmann法的模型在伴有流動(dòng)的相變問(wèn)題中的正確性,本文以固-液相變基準(zhǔn)測(cè)試問(wèn)題——方腔熔化進(jìn)行驗(yàn)證.方腔內(nèi)部初始化為固體狀態(tài),設(shè)置左壁面恒為高溫Th,右壁面恒為低溫Tb,上下壁面均為絕熱狀態(tài),無(wú)量綱參數(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ū)域不斷變大,熔化程度增強(qiáng);同時(shí),由于熱浮力的作用,方腔上部的熔化速度大于底部的熔化速度.將不同F(xiàn)o數(shù)下相變界面提取出來(lái),并與熱焓法[19]及自適應(yīng)網(wǎng)格法[8]得到的結(jié)果進(jìn)行對(duì)比,獲得了良好的一致性,如圖5所示,表明將本文模型應(yīng)用在伴有流動(dòng)的相變過(guò)程中是可行的.

圖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 仿真研究與分析

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

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

5.1 提拉速度對(duì)相變界面的影響

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

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

圖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)對(duì)相變界面的影響

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

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

圖9顯示,隨著晶體旋轉(zhuǎn)參數(shù)的增大,相變界面形態(tài)從最初的凸向熔體變?yōu)椴糠滞瓜蚓w,表明晶體旋轉(zhuǎn)對(duì)改善固-液界面附近熱場(chǎng)的均勻性起主要作用,這與文獻(xiàn)[20]的結(jié)果一致,達(dá)到了調(diào)節(jié)晶體轉(zhuǎn)速來(lái)改善相變界面形態(tài)的目的.當(dāng)Rex=4000時(shí),相變界面位置與自由表面(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)對(duì)相變界面的影響

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

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

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

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

實(shí)際晶體生產(chǎn)中,為了減少位錯(cuò)的產(chǎn)生,相變界面應(yīng)該越平坦越好,而不是凸向熔體或凸向晶體.因?yàn)槠教沟南嘧兘缑姹硎揪w中徑向溫度梯度較小,從而保證晶體內(nèi)部熱應(yīng)力較小,避免由于熱應(yīng)力過(guò)大而產(chǎn)生的位錯(cuò)缺陷.以上的分析結(jié)果表明,晶體旋轉(zhuǎn)對(duì)改善相變界面附近熱場(chǎng)的均勻性起主要作用,坩堝旋轉(zhuǎn)作為一種輔助調(diào)節(jié)手段,可通過(guò)合理地配置晶體旋轉(zhuǎn)參數(shù)和坩堝旋轉(zhuǎn)參數(shù)實(shí)現(xiàn)對(duì)相變界面形態(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ù)對(duì)相變界面形態(tài)的影響時(shí),發(fā)現(xiàn)當(dāng)坩堝與晶體反向旋轉(zhuǎn)時(shí),在一組晶體旋轉(zhuǎn)參數(shù)下,總能找到一組相應(yīng)的坩堝旋轉(zhuǎn)參數(shù)使得相變界面形態(tài)趨于平坦;當(dāng)坩堝與晶體同向旋轉(zhuǎn)時(shí),相變界面形態(tài)隨坩堝旋轉(zhuǎn)參數(shù)變化的規(guī)律性不明確,同一組晶體旋轉(zhuǎn)參數(shù),可能存在兩組坩堝旋轉(zhuǎn)參數(shù)使得相變界面變得平坦.因此,考慮到工程應(yīng)用中也大多采用反向旋轉(zhuǎn)技術(shù),本節(jié)僅針對(duì)反向旋轉(zhuǎn)的情況討論.為了找到相變界面平坦時(shí)對(duì)應(yīng)的晶體和坩堝旋轉(zhuǎn)參數(shù),根據(jù)實(shí)際生產(chǎn)中工藝參數(shù)的調(diào)節(jié)范圍,本文對(duì)不同晶體旋轉(zhuǎn)和坩堝旋轉(zhuǎn)參數(shù)下的相變界面形態(tài)進(jìn)行了大量的仿真,并進(jìn)行了定量分析.

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

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

表2為反向旋轉(zhuǎn)時(shí)不同Rex數(shù)和不同Rec數(shù)配置下相變界面偏差絕對(duì)值的均值和偏差的標(biāo)準(zhǔn)差.將每行結(jié)果中均值和標(biāo)準(zhǔn)差最小時(shí)對(duì)應(yīng)的Rex數(shù)和Rec數(shù)提取出來(lái),發(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可分析出:當(dāng)Rex<2000時(shí),隨著Rex數(shù)的增大,要保持相變界面平坦,Rec數(shù)也應(yīng)該增大;當(dāng)Rex>2000時(shí),隨著Rex數(shù)的增大,要保持相變界面平坦,Rec數(shù)應(yīng)該減小.這是因?yàn)?當(dāng)Rex數(shù)較小時(shí),晶體旋轉(zhuǎn)主導(dǎo)的順時(shí)針渦流較小,尚無(wú)法傳遞足夠的熱量到相變界面以實(shí)現(xiàn)相變界面向晶體側(cè)移動(dòng),相變界面完全凸向熔體.因此,需要較強(qiáng)的坩堝旋轉(zhuǎn)作用,將晶體旋轉(zhuǎn)主導(dǎo)的順時(shí)針渦流向晶體一側(cè)擠壓,改變晶體下方的溫度梯度,從而改變相變界面形態(tài).當(dāng)Rex數(shù)較大時(shí),晶體旋轉(zhuǎn)主導(dǎo)的順時(shí)針渦流本身就較強(qiáng),甚至?xí)耆紦?jù)晶體下方,直接影響相變界面形態(tài),此時(shí),坩堝旋轉(zhuǎn)作為一種輔助的調(diào)節(jié)手段,僅需要較小的坩堝旋轉(zhuǎn)參數(shù),對(duì)晶體下方順時(shí)針渦流的位置進(jìn)行微調(diào),就可以很大程度上調(diào)節(jié)相變界面形態(tài).

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

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

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

6 結(jié) 論

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

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

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

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