直拉法硅晶體生長中單晶爐坩堝內(nèi)熔體的數(shù)值模擬

弋英民，張潼，劉丹

(西安理工大學(xué) 自動(dòng)化與信息工程學(xué)院，陜西 西安 710048)

大尺寸、高品質(zhì)硅單晶是制造微細(xì)納米級(jí)集成電路芯片的關(guān)鍵材料，目前我國全部依賴進(jìn)口。硅單晶爐作為生長電子級(jí)硅單晶材料的主要設(shè)備，其工藝復(fù)雜，技術(shù)集成度高，實(shí)現(xiàn)困難，在集成電路生產(chǎn)流程中處于首要地位。因此，針對(duì)制備單晶硅的各種基礎(chǔ)研究和仿真模擬成為了學(xué)者們研究的熱點(diǎn)。單晶硅按生長方法分有直拉法和區(qū)熔法等[1]，直拉法也稱CZ法，目前廣泛應(yīng)用于晶體制造領(lǐng)域。在CZ法制備單晶硅的過程中，影響晶體質(zhì)量的因素有坩堝內(nèi)的硅熔體流動(dòng)[2]、熔體內(nèi)雜質(zhì)分布[3]以及晶體生長面彎月面的形狀[4]等。使用計(jì)算機(jī)模擬晶體的生長環(huán)境在晶體制備的仿真中起到了重要的作用，為了衡量仿真算法的正確性，文獻(xiàn)[5]提出了Wheeler標(biāo)準(zhǔn)問題，Wheeler標(biāo)準(zhǔn)問題是驗(yàn)證仿真算法正確與否的檢驗(yàn)標(biāo)準(zhǔn)。

為了能夠模擬晶體的生長情況，一般采用專用軟件或者編程實(shí)現(xiàn)。文獻(xiàn)[6]和[7]采用FEMAG-CZ對(duì)單晶爐內(nèi)的熱場(chǎng)和氧濃度分布進(jìn)行了模擬；文獻(xiàn)[8]采用CGSIM對(duì)晶體中的微缺陷進(jìn)行了模擬研究；文獻(xiàn)[9]使用Fortran編程，利用Petrov-Galerkin有限元法實(shí)現(xiàn)了不同物性參數(shù)對(duì)晶體生長的影響的模擬；文獻(xiàn)[10]也使用通用編程軟件，模擬了晶體的生長過程，并通過了Wheeler標(biāo)準(zhǔn)的驗(yàn)證。然而，專用軟件需要大量的授權(quán)費(fèi)用，使用通用軟件編程需要較高的編程功底和理論水平，因此本文研究使用通用的有限元軟件來實(shí)現(xiàn)晶體生長相關(guān)的數(shù)值模擬。

本文采用Fluent 6.3，利用有限容積法結(jié)合三階迎風(fēng)QUICK離散格式模擬了單晶爐中坩堝內(nèi)熔體的流場(chǎng)和熱場(chǎng)，并用Wheeler標(biāo)準(zhǔn)檢驗(yàn)了算法的有效性。對(duì)坩堝直徑為600 mm的熔體進(jìn)行了自然對(duì)流下的流場(chǎng)、熱場(chǎng)以及雜質(zhì)場(chǎng)的數(shù)值模擬。

1 數(shù)值模型和物性參數(shù)

直拉法晶體生長的簡化結(jié)構(gòu)如圖1所示。

圖1 晶體生長結(jié)構(gòu)簡化圖

假設(shè)坩堝和晶體為同軸圓柱，坩堝半徑為Rc，坩堝角速度為Ωc，坩堝內(nèi)熔體高度為H，晶體半徑為Rx(Rx

模擬中假定：①晶體等徑生長；②熔體流動(dòng)和換熱是軸對(duì)稱的；③熔體與晶體界面的自由表面為平面；④忽略熔體自由表面張力和旋轉(zhuǎn)時(shí)界面的變化；⑤熔體為軸對(duì)稱的不可壓縮的牛頓流體；⑥不考慮磁場(chǎng)；⑦熔體內(nèi)的流動(dòng)為層流；⑧坩堝和生長界面處的流動(dòng)滿足無滑移條件，浮升力滿足Bousinesq假設(shè)。除密度外，其它物性參數(shù)均為常數(shù)，熔體的流動(dòng)和傳熱滿足N-S方程和能量方程。

該模型引入無量綱變量：

(1)

式中帶*標(biāo)的變量為有量綱變量，(r,z)為質(zhì)點(diǎn)坐標(biāo)；u、v、w分別為質(zhì)點(diǎn)的軸向速度、徑向速度和周向速度；p為壓力，T為溫度；μ為熔體的粘性系數(shù)；ρ為熔體密度；g為重力加速度。

基于之前的假設(shè)，無量綱的控制方程為：

(2)

(3)

(4)

(5)

(6)

本文所采用熔體模型的邊界條件為：

(7)

相應(yīng)的，模型中的無量綱參數(shù)定義如下：

(8)

式中，α為熔體高度與半徑之比；β為晶體與熔體半徑之比；Pr為普朗特?cái)?shù)；Gr為格拉斯霍夫數(shù)；Rex和Rec分別是因晶體和坩堝旋轉(zhuǎn)產(chǎn)生的雷諾數(shù)。

仿真中采用的硅熔體的物性參數(shù)如表1所示[10]。

表1 物性參數(shù)表

2 Wheeler標(biāo)準(zhǔn)化問題

在1989年召開的北約晶體生長計(jì)算模擬大會(huì)上，大會(huì)主席Wheeler提出了四個(gè)測(cè)試問題，即Wheeler標(biāo)準(zhǔn)問題。這些標(biāo)準(zhǔn)問題包括四組情況：A組只考慮晶體旋轉(zhuǎn)；B組考慮晶體和坩堝共同旋轉(zhuǎn)；C組只考慮自然對(duì)流，晶體坩堝不旋轉(zhuǎn)；D組考慮自然對(duì)流和晶體旋轉(zhuǎn)，坩堝不旋轉(zhuǎn)。每組又分三個(gè)不同的算例，Wheeler標(biāo)準(zhǔn)問題的參數(shù)如表2所示。本文計(jì)算了表2中的12組算例。

表2 Wheeler標(biāo)準(zhǔn)問題參數(shù)

3 數(shù)值模擬及分析

本文使用Fluent 6.3[11]進(jìn)行數(shù)值模擬，采用層流模型，壓力和速度的解耦合用Fluent 6.3中的SIMPLE算法，擴(kuò)散項(xiàng)差分采用二階截差中心差分格式，對(duì)流項(xiàng)差分采用三階截差QUICK格式，壓力插值采用PRESTO！格式。對(duì)Wheeler標(biāo)準(zhǔn)問題的12個(gè)算例進(jìn)行了數(shù)值模擬。

算例A結(jié)果的流函數(shù)如圖2所示。在這組算例中，熔體內(nèi)部只有一個(gè)由晶體旋轉(zhuǎn)產(chǎn)生的一個(gè)順時(shí)針的渦，熔體沿軸線向上運(yùn)動(dòng)，直至晶體邊界，然后流向坩堝側(cè)壁，晶體附近的流體由于粘性的作用和晶體一起旋轉(zhuǎn)，產(chǎn)生從軸心指向坩堝壁的離心力，在離心力的作用下，熔體流向坩堝壁，遇到坩堝壁阻擋，被迫下流。轉(zhuǎn)速將渦流的中心移向坩堝側(cè)壁，偏向右上角，并發(fā)生形變。同時(shí)，速度最快的區(qū)域從左上角移至右下角，速度最小的區(qū)域在左下角擴(kuò)大。

圖2 算例A的流函數(shù)圖

在本組算例中，隨著晶體轉(zhuǎn)動(dòng)速度的增加，流函數(shù)變大，這說明渦流強(qiáng)度增加。算例A1、A2、A3的流函數(shù)最大值分別為8.79E-5、1.18E-4、1.47E-4。

算例B結(jié)果的流函數(shù)如圖3所示。在算例B中，熔體內(nèi)部出現(xiàn)兩個(gè)方向相反的渦流，左上角的順時(shí)針渦流是由晶體旋轉(zhuǎn)形成的，右下角的逆時(shí)針渦流是由坩堝旋轉(zhuǎn)造成的；隨著晶體旋轉(zhuǎn)速度和坩堝旋轉(zhuǎn)速度的同時(shí)增大，雖然晶體轉(zhuǎn)速遠(yuǎn)大于坩堝轉(zhuǎn)速，但是右下角的渦流仍主導(dǎo)著整個(gè)流場(chǎng)；晶體生長的理想條件是熔體區(qū)域只有一個(gè)渦流，顯然，這樣的情況是不利于晶體生長的。

圖3 算例B的流函數(shù)圖

在本組算例中，隨著晶體轉(zhuǎn)動(dòng)速度和坩堝轉(zhuǎn)動(dòng)速度的增加，流函數(shù)變大，說明渦流強(qiáng)度顯著增加。算例B1、B2、B3的流函數(shù)最大值分別為8.25E-7、1.86E-6、3.92E-5。加入晶體和坩堝旋轉(zhuǎn)是為了平衡自然對(duì)流產(chǎn)生的渦流，但晶體和坩堝的轉(zhuǎn)速過大，產(chǎn)生的渦流強(qiáng)度過大，也對(duì)晶體等徑生長帶來擾動(dòng)，因此晶體和坩堝的轉(zhuǎn)速應(yīng)控制在一定范圍之內(nèi)。

算例C中只有浮升力引起的自然對(duì)流，熔體內(nèi)部由一個(gè)逆時(shí)針的渦占據(jù)，隨著Gr數(shù)的增大，渦流也在變大，最大速度移向坩堝側(cè)壁，渦流的中心幾乎沒有移動(dòng)位置。算例C結(jié)果的流函數(shù)圖如圖4所示。Gr數(shù)在108時(shí)，層流模型無法得到收斂的解，證明層流模型不適用，Gr數(shù)在108以上的數(shù)量級(jí)時(shí)，應(yīng)該使用湍流模型。

圖4 算例C的流函數(shù)圖

算例D結(jié)果的流函數(shù)圖如圖5所示。算例D相比算例C加上了晶體旋轉(zhuǎn)，結(jié)果加上晶體旋轉(zhuǎn)的自然對(duì)流和只有自然對(duì)流的情況，流函數(shù)不管是在形狀上，還是在數(shù)值上都沒有什么變化，和C組的幾乎一樣。

圖5 算例D的流函數(shù)圖

通過Wheeler標(biāo)準(zhǔn)問題檢驗(yàn)后，證明此方法可用，于是，本文將此方法應(yīng)用于大直徑硅晶體，對(duì)坩堝半徑為300 mm內(nèi)的硅熔體，進(jìn)行自然對(duì)流下的流場(chǎng)、熱場(chǎng)和氧濃度的數(shù)值模擬。熔體高度與半徑之比α=0.6；晶體與熔體半徑之比β=0.4。

圖6為坩堝半徑為300 mm時(shí)，硅熔體在自然對(duì)流下的流函數(shù)圖?？梢钥闯觯骱瘮?shù)的基本形狀和層流模型的基本相同，是充滿整個(gè)熔體區(qū)域的一個(gè)逆時(shí)針大渦，但是渦心更偏向坩堝側(cè)壁，此渦是由于坩堝底部和側(cè)壁的溫度比熔體高，坩堝壁附近區(qū)域的熔體受熱后，密度變小，在浮升力的作用下，沿著側(cè)壁上升，然后在自由表面和生長界面處受冷，密度變大，下沉至坩堝底部，從而形成完整的渦流。

圖6 硅熔體流函數(shù)圖

圖7為坩堝半徑為300 mm時(shí)，硅熔體在自然對(duì)流下的等溫圖。低缺陷晶體生長的溫度條件是在整個(gè)晶體生長過程中固液交界面保持平坦，這要求等溫線近似垂直于生長方向，即等溫線比較平坦。圖7中，靠近晶體生長面的等溫線比較平坦，隨著遠(yuǎn)離晶體生長界面，等溫線變得更為陡峭。自然對(duì)流的湍流會(huì)加劇溫度的波動(dòng)，而溫度的波動(dòng)是產(chǎn)生生長條紋的重要原因，因此要獲得高品質(zhì)單晶硅，控制其溫度波動(dòng)是關(guān)鍵。

圖7 硅熔體等溫圖

圖8顯示的是在自然對(duì)流作用下坩堝半徑為300 mm的硅熔體中的氧原子濃度分布圖。可以看出，在熔體內(nèi)部的大部分區(qū)域，氧原子濃度都很小。并且，由于渦流在自然對(duì)流的作用下，沿坩堝壁經(jīng)由自由表面再流向?qū)ΨQ軸中心，因此一部分氧雜質(zhì)也會(huì)從自由表面揮發(fā)出去，在圖中就可以看出大部分的氧雜質(zhì)聚集在坩堝底部和側(cè)壁。如果熔體內(nèi)的流動(dòng)加強(qiáng)，那么氧原子雜質(zhì)有可能會(huì)來不及揮發(fā)而被帶入熔體，從而形成熔體內(nèi)的雜質(zhì)，因此抑制自然對(duì)流的流動(dòng)強(qiáng)度對(duì)降低熔體內(nèi)的氧雜質(zhì)至關(guān)重要。

圖8 硅熔體氧濃度圖

根據(jù)本研究結(jié)論，通過在直徑300 mm硅單晶爐修改拉晶的SOP表有關(guān)參數(shù)，進(jìn)行不同SOP溫度參數(shù)下的拉晶實(shí)驗(yàn)如圖9所示。

圖9 在不同拉晶SOP參數(shù)下的晶棒

由圖9可以看出，在晶體生長界面處溫度波動(dòng)大，會(huì)導(dǎo)致生長處的晶體無法成晶，甚至斷裂。而當(dāng)控制晶體生長界面(彎月面)溫度為硅熔體成晶溫度1 420℃±0.1℃時(shí)，可以長出完美的硅單晶。

4 結(jié) 論

本文通過Fluent6.3對(duì)Wheeler標(biāo)準(zhǔn)問題中的12個(gè)算例進(jìn)行了數(shù)值計(jì)算，結(jié)果表明，使用Fluent計(jì)算的結(jié)果與文獻(xiàn)[10]中的基本一致。本文的仿真實(shí)驗(yàn)結(jié)果表明，對(duì)于直徑600 mm坩堝必須控制合適的晶轉(zhuǎn)和堝轉(zhuǎn)范圍，以保證晶體生長，控制溫度的波動(dòng)可以抑制生長條紋，抑制自然對(duì)流的可以降低熔體內(nèi)的氧雜質(zhì)含量。本研究下一步將在自然對(duì)流下增加晶體和坩堝旋轉(zhuǎn)，并增加間壁對(duì)氧濃度的抑制方面的研究。

參考文獻(xiàn)：

[1] 蔣榮華,肖順珍. 半導(dǎo)體硅材料的進(jìn)展與發(fā)展趨勢(shì)[J].四川有色金屬, 2000,3:1-7.

Jiang Ronghua,Xiao Shunzhen. Progress and development trend of semiconductor Si materials[J]. Sichuan Nonferrous Metals,2000,3:1-7.

[2] 宇慧平,隋允康,張峰翊,等. 紊流模型模擬分析旋轉(zhuǎn)對(duì)提拉大直徑單晶硅的影響[J]. 人工晶體學(xué)報(bào), 2004,33(5): 835-840.

Yu Huiping,Sui Yunkang,Zhang Fengyi, et al. Numerical simulation on the effect of rotation in a czochralski silicon crystal growth with a turbulence model[J]. Journal of Synthetic Crystal, 2004,33(5): 835-840.

[3] Ren Bingyan,Zhao Long,Zhao Xiuling, et al. Effects of argon gas flow rate and guide shell on oxygen concentration in czochralski silicon growth [J]. Rare Metals, 2006, 25(1): 7-10.

[4] 滕冉,戴小林,肖清華,等. 大直徑硅單晶生長過程中固/液界面形狀及熔體流動(dòng)的數(shù)值分析[J]. 人工晶體學(xué)報(bào),2013,42(4): 611-615.

Teng Ran,Dai Xiaolin,Xiao Qinghua, et al. Numerical simulation of melt/crystal interface and melt flow during large diameter single silicon crystal growth[J]. Journal of Synthetic Crystal, 2013,42(4): 611-615.

[5] Wheeler A A. Four test problems for the numerical simulation of flow in czochraski crystal growth [J]. Journal of Crystal Growth,1990,102(4):691-695.

[6] 王慶,張婷曼. 基于FEMAG-CZ熱場(chǎng)仿真軟件的單晶爐熱場(chǎng)分析[J]. 計(jì)算機(jī)光盤軟件與應(yīng)用, 2011,17: 65-66.

Wang Qing,Zhang Tingman. Single crystal furnace thermal field analysis on FEMAG-CZ thermal field simulation software[J]. Computer CD Software and Applications, 2011,17: 65-66.

[7] 常麟,周旗鋼,戴小林,等. CUSP磁場(chǎng)對(duì)直拉硅單晶氧濃度分布影響的數(shù)值模擬[J].稀有金屬,2011,35(6): 909-915.

Chang Lin,Zhou Qigang,Dai Xiaolin, et al. Numerical simulation of CUSP magnetic field on oxygen concentration distribution in CZ-Si crystal growth[J]. Chinese Journal of Rare Metals, 2011,35(6): 909-915.

[8] 曾慶凱,關(guān)小軍,潘忠奔,等. Ф400mm直拉硅單晶生長過程中氧濃度對(duì)微缺陷影響的數(shù)值模擬[J]. 人工晶體學(xué)報(bào),2011,40(5):1150-1156.

Zeng Qingkai,Guan Xiaojun,Pan Zhongben, et al. Simulation of microdefect with different oxygen concentration in Ф400 mm CZ silicon crystal growth[J]. Journal of Synthetic Crystal,2011,40(5): 1150-1156.

[9] 李友榮,魏東海,余長軍,等. 物性參數(shù)對(duì)硅單晶Czochralski生長過程的影響[J]. 熱科學(xué)與技術(shù),2006,5(4): 351-355.

Li Yourong,Wei Donghai,Yu Changjun, et al. Effects of thermophysical properties on silicon single crystal czochraski growth processes[J]. Journal of Thermal Science and Technology,2006,5(4): 351-355.

[10] 宇慧平. MCZ大直徑單晶體生長的數(shù)值模擬及控制參數(shù)優(yōu)化[D]. 北京:北京工業(yè)大學(xué),2005.

Yu Huiping. The numerical simulation of the large diameter MCZ crystal growth and the optimization of the operational parameters[D]. Beijing:Beijing University of Technology,2005.

[11] 李國棟,哈岸英,鐘小彥,等. 基于FLUENT 的滲流場(chǎng)數(shù)值模擬分析[J]. 西安理工大學(xué)學(xué)報(bào),2011,27(3): 317-320.

Li Gudong,Ha Anying,Zhong Xiaoyan, et al. Numerical simulation of seepage field based on FLUENT[J]. Journal of Xi’an University of Technology, 2011,27(3): 317-320.