數值模擬晶棒成長之熱流分析

國

立

交

通

大

學

機械工程學系

碩

士

論

文

數值模擬晶棒成長之熱流分析

Numerical Simulations of Crystal Growth

研 究 生：陳柏瑋

數值模擬晶棒成長之熱流分析

Numerical Simulations of Crystal Growth

研 究 生：陳 柏 瑋 Student：Po-Wei Chen 指導教授：陳 慶 耀 Advisor：Ching-Yao Chen 國 立 交 通 大 學 機 械 工 程 學 系 碩 士 論 文 A Thesis

Submitted to Department of Mechanical Engineering College of Engineering

National Chiao Tung University in Partial Fulfillment of The Requirements

for the Degree of Master of Science

in

Mechanical Engineering July 2010

數值模擬晶棒成長之熱流分析

學生:陳柏瑋 指導教授:陳慶耀 國立交通大學機械工程學系碩士班 摘要 本研究模擬在 cusp 磁場下，變化磁場強度及變化電磁線圈的高度位 置對 Cz 長晶法之矽熔湯流場和溫度場的影響，同時也探討固液界面形狀 和矽晶棒拉速及其溫度梯度的比值(V/G 值) 變化。統御方程式包含動量 方程式、熱傳方程式以及電流流動方程式。紊流效應和因氣體流動而在氣 液界面上所產生的剪切應力效應均考慮於本研究中。本文採用商用軟體 CGsim (版本 9.02) 進行計算。從模擬的結果可以推測出磁線圈往上移動 較多時會導致矽熔湯流速減慢、坩堝中心處的溫度梯度上升、坩堝底部的 溫度升高以及矽熔湯流場的不穩定。坩鍋內低溫區的範圍會朝向固液界面 的方向縮小。此外，界面的變化差異和 V/G 值均會隨調整磁線圈的高度而 降低。 關鍵字：電腦模擬、熱傳、紊流對流、Cz 長晶法、矽晶棒

誌 謝

時光匆匆、稍縱即逝，兩年充實的碩士生涯即將劃上句點，又將是人 生下一個階段的開始，回首剛進研究所時的期待，以及現在收成的喜悅， 這一切都要感謝許多人對我的提攜與幫助。 這篇論文能夠順利完成，首先要感謝我的指導老師－陳慶耀教授。感 謝老師在求學期間悉心的指導，給予學生明確的方向，帶領我渡過研究中 關卡與挑戰，僅此表示萬分的誠摯與感謝。再來要感謝口試委員葉俊良教 授與許立傑教授在百忙之中撥空前來給予論文及研究上許多的建議與指 教，使得本論文得以更趨於臻善，在此致上誠摯的謝意。 此外，還要感謝奇軒、詒翔、宏元等諸位學長在研究上給予指導及提 供相當寶貴的意見，以及裕盛、佳芬、啟鳴，很高興大家有緣，兩年同在 一個實驗室，因為你們充實了我生活中的點點滴滴，還有謝謝學弟文良、 世聰、國安、杰民幫忙分擔實驗室的業務。有你們的陪伴，讓我的生涯更 加的歡樂、多采多姿。 最後，我要感謝永遠支持我的家人，爸爸、媽媽謝謝你們養育之恩與 栽培，仟玓、廷嘉、紓百謝謝你們的照顧與關心，好友 M4D1 的品維、宏 偉、垠元、建安謝謝你們的情義相挺，讓我能全力以付的完成學業，再次 獻上最真誠謝意，謝謝大家!

目 錄

摘要 ... I 誌 謝 ... II 目 錄 ... III 圖目錄 ... V 表目錄 ... VII 符號說明 ... VIII 第一章 緒 論 ... 1 1.1 前言 ... 1 1.2 柴氏單晶生長法簡介 ... 2 1.3 文獻回顧 ... 3 1.4 研究目的與動機 ... 5 第二章 物理模型與系統描述 ... 10 2.1 分析假設 ... 10 2.2 矽熔湯流動現象 ... 10 2.3 長晶過程熱傳遞行為 ... 13 2.4 影響長晶熱流場的因素 ... 13

第三章 數值方法 ... 19 3.1 數值方法 ... 19 3.2 固液界面計算 ... 21 3.3 收斂條件 ... 22 第四章 結果與討論 ... 27 4.1 模式驗證 ... 27 4.2 網格測試 ... 28 4.3 改變拉升速度之模擬 ... 28 4.3.1 拉升速度對矽熔湯流場和溫度場影響... 29 4.3.2 拉升速度對固液界面與 V/G 參數之影響... 30 4.4. 改變 Cusp 磁場強度之模擬 ... 31 4.4.1 Cusp 磁場強度對矽熔湯流場和溫度場影響 ... 32 4.4.2 Cusp 磁場強度對矽熔湯固液界面與 V/G 參數之影響 .... 32 4.5 改變 Cusp 磁線圈位置之模擬 ... 33 4.5.1 Cusp 磁線圈位置對矽熔湯流場和溫度場影響 ... 34 4.5.2 Cusp 磁線圈位置對固液界面變化與 V/G 參數 ... 34 第五章 結論 ... 56 參考文獻 ... 57

圖目錄

圖 1-1 典型 CZ 矽晶爐之示意圖 ... 7 圖 1-2 傳統柴氏矽晶爐 ... 8 圖 1-3 不同 Gr number 與固液界面關係圖 ... 8 圖 1-4 Purge tube 影響整體流場圖 ... 9 圖 1-5 400mm 晶棒固液界面形狀、V/G 值 ... 9 圖 2-1 柴氏長晶法熱傳示意圖 ... 17 圖 2-2 熔湯流動型態(a)垂直自然對流(b)水準自然對流(c)晶棒拉速影響 的強制對流(d)晶棒旋轉影響的強制對流(e)坩堝旋轉影響的強制對流.. 18 圖 2-3 水準磁場影響之熔湯流動 ... 18 圖 3-1 網格中心點及離散面相關位置示意圖 ... 25 圖 3-2 研究流程圖 ... 26 圖 4-1 模式驗證比較圖 ... 36 圖 4-2 CGSim 模組網格測試 ... 36 圖 4-3 Flow 模組網格測試 ... 37 圖 4-4 二維軸對稱穩態全域熱傳計算之網格與溫度分佈圖 ... 37 圖 4-5 二維穩態軸對稱流動模組網格與溫度分佈圖 ... 38

圖 4-8 拉升速度 0.45mm/min 之模擬，(a)溫度場 (b)流場 ... 41 圖 4-9 拉升速度 0.50mm/min 之模擬，(a)溫度場 (b)流場 ... 42 圖 4-10 不同拉升速度下之固液界面 ... 43 圖 4-11 不同拉升速度下之 V/G 值 ... 43 圖 4-12 cusp 磁場中心示意圖 ... 44 圖 4-13 cusp 磁場位置圖 ... 44 圖 4-14 B=0.00T 之模擬，(a)溫度場 (b)流場 ... 45 圖 4-15 B=0.05T 之模擬，(a)溫度場 (b)流場 ... 46 圖 4-16 B=0.10T 之模擬，(a)溫度場 (b)流場 ... 47 圖 4-17 不同磁場強度下之固液界面 ... 48 圖 4-18 不同磁場強度下之 V/G 值 ... 48 圖 14-19 cusp 磁線圈位置移動圖 ... 49 圖 4-20 B=0.05T dy=0cm 之模擬，(a)溫度場 (b)流場... 50 圖 4-21 B=0.05T dy=5cm 之模擬，(a)溫度場 (b)流場... 51 圖 4-22 B=0.05T dy=10cm 之模擬，(a)溫度場 (b)流場... 52 圖 4-23 B=0.05T dy=-5cm 之模擬，(a)溫度場 (b)流場... 53 圖 4-24 B=0.05T dy=-10cm 之模擬，(a)溫度場 (b)流場... 54 圖 4-25 不同高度的磁線圈位置之固液界面 ... 55 圖 4-26 不同高度的磁線圈位置之 V/G 值 ... 55

表目錄

符號說明

文字符號 B _{磁通量密度} p c _比熱 eff D _{有效擴散係數} d _坩堝直徑 g _重力常數 h _熔湯深度 j _感應電流 m _分子量 0 P _參考壓力 g R 理想氣體常數 T _溫度 希臘符號

β

_{熱膨脹係數} 0

ρ

參考質量密度

ρ

_質量密度

τ

應力張量 eff

μ

有效黏滯係數

ν

_{動力黏滯係數}

λ

熱擴散係數 eff

λ

有效熱擴散係數 s

ω

晶棒轉速 c

ω

坩堝轉速

γ

_表面張力 Φ _磁勢能 σ _電導率 無因次參數 Gr 葛拉秀夫數( 3_/ 2 TgL

β

Δ

ν

) Pr 普朗特數(

ν α

/ ) Re 雷諾數( 2_/ sr k

ω

υ

) Ta 泰勒數₍₂ 2_{/ )}2 ch k

ω

υ

Ma 馬若哥尼數(Δ

_ρυ

Td (_∂∂

γ

))

第一章 緒 論

1.1 前言

現今環保意識的抬頭與能源枯竭造成能源危機，基於環保與節能需 求之發展趨勢，取之不盡，用之不竭的太陽能源又開始受到大家的重 視。在太陽能轉換為光能方面，大部分利用太陽能板把光能轉會為電 能，矽系的太陽能電池材料，主要可以分為單晶矽、多晶矽、非晶矽3 大類，單晶矽的太陽能電池因轉換效率較多晶矽、非晶矽太能能電池 高，以目前條件來估算，增加1%的轉換效率，就可以降低 5~7%的生產 成本，因此開發高效率的單晶矽太陽能電池為結晶矽太陽能電池的發展 方向。 由於單晶矽晶圓的需求量越來越大，因此需要大尺寸的晶棒生產來 滿足需求及降低生產成本，而矽單晶必須由特定之晶體生長法來生長出 單晶晶棒，再將晶棒經由切割研磨等製程，加工成矽晶圓作為太陽能基 板。這些長晶法中以柴氏單晶生長法(Czochralski single crystal growth method)是目前單晶生長法中較快速，經濟且大尺寸之生長方法。大尺寸 矽單晶不只是依靠傳統的CZ法製程，更需要有精密的控制設備與良好的 長晶爐熱場設計，才能達到生長出高品質的矽單晶圓。 本研究使用Basic CGSim模擬軟體來進行柴氏法生長單晶矽之氣 場、熱場及流場三場耦合的模擬研究與分析。也將探討長晶過程中磁場 及熔湯熱流場分佈情況及各項參數的影響，作為將來開發及生長高品質

1.2 柴氏單晶生長法簡介

在1950 年，貝爾實驗室(Bell Lab.)的 Teal 及 Little [1]兩人將 Czochralski [2]方法應用在生長鍺及矽單晶上，雖然這種拉單晶的方 法，1951 年 Buckley 發表”Crystal Growth”一書之後，有關提拉方式的 長晶法被稱為CZ(Czochralski)法 [3]，但嚴格講起來，這方法應該被 稱為Teal-Little 法，因為當初 Czochralski 只是在研究固液界面的結晶 速度，而Teal 及 Little 則是在研究拉單晶的方法。在電子元件材料超 過98%都是使用單晶矽。在工業界上，有兩種方法可用來生長單晶 矽。其中CZ 法佔了約 85%，其他部分則是由浮融法(Float Zone)生長 之。由CZ 法生長出的矽晶，主要用來生產低功率的積體電路元件， 例如:DRAM、SDRAM、ASIC 等。而 FZ 法生長出的矽晶，則主要用 在高功率的電子元件，例如:電晶體、閘流晶體管等。CZ 法所以比 FZ 法更普遍被半導體工業採用，主要在於它的高含氧量提供了強化晶片 及gettering 兩種優點;另外一個原因是 CZ 法比 FZ 法更易生產出大尺 寸的單晶棒。 在一般業界使用柴氏矽晶爐體內的主要元件如圖 1-1 所示分為矽 晶爐體、加熱構件、坩堝與上升旋轉機構及控制系統等。柴氏法是屬 於熔融液生長法(Melting Growth)的一種，熔融液生長過程牽涉固體與 液體的相變化過程，熔融液在固定的條件下沿晶種方向凝固，通過固 液界面的移動逐漸形成規則性的原子堆積，而得到高品質的單晶結 構。藉由晶種或自發成核的繼續生長，固液界面需處於過冷的狀態，

1.3 文獻回顧

單晶矽可藉由CZ 法獲得，但在 1950 年當時，生長過程中差排 (Dislocation)卻往往會出現在晶體中。在 1959 年，Dash [5]發明瞭一 種可以完全消除差排的方法(Dash Technique)，這種零差排 (Dislocation-free)的方法，才使得生長大尺寸晶棒成為可能。目前晶棒 的尺寸已由1950 年代的 1 吋增加到現在的 12 吋，預計在未來的 10 年內會開始生產18 吋的晶棒。1960 年以後，以柴氏法生長矽單晶已 經陸陸續續的建立穩定的機制。柴氏法在熔湯傳導機制及晶體的組成 之探討、固液界面形狀對流場的影響、Hot zone 的設計改善等，都顯 著的改善晶體在品質上的要求，在生長控制上也朝向全自動化生長。 隨著晶棒的尺寸增大，熔湯內流場相較越來越複雜，使得控制參數更 為困難，所以在近二十年來紛紛有許多學者投入研究柴氏法長晶系統 的熱流數值模擬，模擬熱場的傳遞模式及坩堝內矽熔湯流動行為，有 助於利用數值模擬接近真實的晶體生長熱場分佈情形。 早在 1986 年時，Derby 與 Brown [6，7]利用有限元素法模擬柴氏 法生長矽單晶時的爐體內的熱傳行為之外，主要著重在固液界面形狀 與熔湯內的流動彼此間的關聯。同一年，Dupret 等人 [8]針對柴氏矽 晶爐體作全域熱傳遞行為的分析(Global Heat Transfer)，使用有限元素 法分析柴氏矽爐體的整體熱場分佈。

而在 1997 年，Dupret 等人 [9，10]利用暫態的方式來模擬柴氏 晶體生長系統，在此文獻中提高了計算時的速度，同時也更詳細考慮 系統內的熱傳遞現象，包括有熱傳導及熱輻射等，在這些暫態模擬的 文獻之中皆提出一個重要的觀念，就是晶體生長的系統是一個暫態的

佈，還是要考慮成暫態的系統較為合理。 在晶棒與矽熔湯交接處的介面，形狀的凹形(concave)或凸形 (convex)會造成缺陷、差排、氣泡產生，所以固液界面現象在長晶過 程中視為很重要的影響，在1996 年，Wu [11]等人分別提出自然對流、 晶棒旋轉、坩堝旋轉造成固液界面形狀的變化，也解釋當自然對流增 加時，固液界面會產生凸型，如圖 1-3，而增加晶棒轉速會產生凹型。 在2000 年時，Brown 學者等人 [12]模擬柴氏生長矽單晶時使用 有限元素法(FEM)計算熱傳模式和有限容積法(FVM)計算矽熔湯內紊 流模式，使用此種方式分割整體熱場與熔湯流場，使得較複雜的熔湯 流場得到更高的精確度，也使計算結果同時考慮熔湯內的紊流行為， 去分析紊流對於熔湯的熱場分佈所造成的影響，更進一步考慮紊流下 熔湯內的氧析出物的分佈是否有所改變。 在 2001 年，Tankano 學者等人 [13]利用 FEMAG 軟體，模擬拉升 不同晶棒直徑與熱應力的關係，並加以解釋當晶棒直徑越大時，熱應 力也會跟著增加，而影響到固液界面的形狀，造成性質的不穩定。此 外，模擬在氬氣流道上加裝一Purge tube 可以減少渦流的的產生，如 圖1-4。 在矽熔湯中不穩定的熱對流容易導致固液界面的不穩定性及雜 質分析的微觀偏析現象(growth striation)。由於大部分的金屬及半導體

利用柴氏長晶法(Czochralski method)生產晶棒的過程中，當晶棒 直徑提升至300公厘或更大尺寸450公厘時，坩堝中矽熔湯自然對流效 應的影響是工程師所必需注意的地方 [15]。若使用不同類型的穩態 直流電 (direct current, DC) 磁場 (如垂直式、水平式或者cusp式) 於 長晶製程中，可以改變矽熔湯的流動與對流熱質傳現象的特性。例 如，當磁場強度B大於0.014T時，矽熔湯的流動會形成層流 [16]；當 晶棒長度超過300公厘時，若不加磁場效應則矽熔湯流動會變成紊流 [17]。因此，使用磁場來控制矽熔湯內的自然對流效應再進而探討最 佳的固液界面形狀以及V/G值是ㄧ種有效率的方法(V為晶棒拉速，G 為 晶 棒 溫 度 梯 度)。在2007年，Kalaev [23]學者模擬大尺寸晶棒 400mm，施加橫向磁場和cusp磁場來有效抑制熔湯內的自然對流，探 討有無開起磁場之流場、溫度場、固液界面形狀及V/G的好壞，如圖 1-5。 然而，長晶設備及其實驗和磁場裝置的開發成本相當昂貴，而且 在實際製程中固液界面也不易觀測。在過去數十年間，數值模擬已變 成一種用來降低實驗成本以及提供相關物理現象與定性分析的利 器。因此，許多關於使用磁場影響矽熔湯流動與溫度分佈進而預測晶 棒品質的穩態與非穩態研究已經發表出來 [18-32]。

1.4 研究目的與動機

柴氏長晶法用在生長均勻大直徑矽單晶的研究已經有許多學者 累積了相當多的經驗，在晶棒尺寸與品質上也相當顯著的進步。由 於半導體產業蓬勃發展，單晶矽晶圓的需求量也日益增加，所以生

擬方法，幫助我們釐清完整的生長機制。 根據摩爾定律的推算，下一代的矽晶爐應由原先晶棒直徑 300mm 到直徑450mm，而目前文獻中尚未有針對柴氏長晶法 450mm 晶棒生 長過程的模擬研究，所以這方面有相當大的空間做深入探討與研究。 此外，我們選用有限元素法為基礎的Basic CGSim 軟體，來進行柴氏 長晶生長矽單晶之氣場、熱場、流場，三場耦合的模擬與研究分析， 來接近實際物理模式之模擬。長晶過程繁雜且耗費成本，加上複雜的 爐體外型對於熱流場分佈會造成相當大的影響，所以事前的模擬與研 究分析有助於未來矽單晶生長製程上有更好的建議。

圖 1-2 傳統柴氏矽晶爐

圖 1-4 Purge tube 影響整體流場圖

第二章 物理模型與系統描述

單晶生長過程中，熱流場是一個相當繁雜的系統，由長晶系統整 體分析而言，由加熱器以熱輻射的方式加熱坩堝，坩堝再經由熱傳導 將熱傳至矽熔湯中，矽熔湯因溫度差使矽熔湯的密度產生變化，因而 進行自然對流的流動行為。因為矽晶爐內大部分時間均處於高溫狀況 下，所以爐體內坩堝、支撐元件及絕熱層的材質及熱傳遞行為的種種 原因均須考慮。材料性質上的研究需由實際的長晶過程而得，費時且 影響晶體生長的因素還要做進一步的探討。所以站在實務的觀點，我 們研究軟體與數值方法的限制下，必須考量到模擬的簡化，使其在複 雜的系統中得以在適當的假設下加以簡化而不失真，如此才能使模擬 在既經濟又有效率的情況下，得到晶棒生長過程的相關分析。

2.1 分析假設

(1)熱傳過程左右對稱，設為二維軸對稱系統 (2)工作流體設定為不可壓縮流體 (3)長晶過程為穩態分析

2.2 矽熔湯流動現象

柴氏長晶系統中的流動現象如下，如圖2-2、圖 2-3: (1)自然對流(Natural Convection) 一般柴氏系統裡，熱源是由坩堝外側的加熱器所提供，造成熔湯

Gr= 3_/ 2 Tgd

β

Δ

ν

(2-1) 其中:

β

:熱膨脹係數

Δ

T

:溫度差 g:重力常數 d:坩堝直徑 根據公式所示，Gr 值越大，自然對流的程度越明顯，Gr 值越小， 矽熔湯流場流動現象趨於緩和。本研究的重力常數=9.81 _/ 2 m s ，熱膨 脹係數=

_{1.43 10 (1/ )}

4

K

−

×

，坩堝直徑=0.962m，矽熔湯內溫度差約 15K~30K 之間，故 Gr 值約

_{1.83 10 ~ 3.66 10}

×

9

×

9 ，相較於Gr 的臨界 值

10

5高出許多，故矽熔湯流動方式很容易出現不穩定性。 Pr=

ν α

/

(2-2) 其中:

ν

:運動黏滯係數

α

:熱擴散係數 根據熱傳理論，Pr 值越大，熔湯內的熱對流(heat convection)對溫 度場之影響越顯著，其矽熔湯內的熱傳遞現象主要由熱對流(heat convection)所主導；反之若 Pr 值越小則熔湯內熱傳遞現象由熱傳導 (heat conduction)主導。本研究所採用矽熔湯的運動黏滯係數為 6 2

3.2 10

×

−

m

/

s

，熱擴散係數為

8.338 10

×

−4，所計算出的 Pr 值為

(2)強迫對流-晶軸旋轉影響 晶軸旋轉會使緊臨固液界面下的熔湯往上流動，再藉著離心力往 外側流動而造成強制對流區。由晶軸旋轉引起的熔湯擾動程度，可由 無因次的 Reynolds 常數表示 2

Re

=

ω

s

r

/

υ

(2-3) (3)強迫對流-坩堝旋轉的影響 坩堝的旋轉將使得坩堝外側的熔湯往中心流動。由坩堝旋轉引起 的對流程度可由無因次的Taylor 常數表示 2 2 (2 _c / ) Ta=

ω

h

υ

(2-4) 其中:

ω

_c

=

坩堝轉速 h=熔湯深度 (4)表面張力對流 由液面的溫度梯度，所造成表面張力的差異，而引起的對流型 態，稱為表面張力對流。其對流程度大小可由無因次的Marangoni 常 數表示 ( ) Td Ma T

γ

ρυ

Δ ∂ = ∂ (2-5) 其中:

γ

∂

由上式可知， d

Bo

∝

d g

2 ，也就是說在地表上較大的長晶系統主 要受自然對流主導。而表面張力對流則在低重力狀態及小長晶系統， 才會突顯其重要性。

2.3 長晶過程熱傳遞行為

柴式法生長矽單晶的系統內包含著各項元件間複雜的熱傳遞行 為，由圖 2-1 中為爐體內有熱傳導、熱對流、熱輻射的熱傳遞行為 示意圖，在長晶過程中包括所有的熱傳遞現象，有固體內部的熱傳 導，如坩堝與支撐座內部的熱傳遞行為。而熱對流行為會發生在矽熔 湯中，矽熔湯內部會因壁面溫度比中心高，底端溫度比液面高，使溫 度梯度造成密度變化而產生熱對流。而矽晶爐大部份熱源由加熱器經 由熱輻射到石墨坩堝壁面藉此加熱矽熔湯，因長晶過程爐體內為高溫 狀態所以熱輻射現像是影響整個熱場分佈很重要的一環，其中包括晶 柱表面、熔湯自由液面及固體元件間所產生的熱輻射效應，皆是影響 爐體內熱場分佈的重要因素。

2.4 影響長晶熱流場的因素

晶棒與矽熔湯之間固液界面的平緩程度，影響著晶棒品質的好 壞，因此固液界面除了與晶棒的拉升速度有關，並且與軸向熱通量相 關，所以靠近固液界面的矽熔湯軸向溫度梯度，也間接反映熔體流動 的穩定性，同時也影響了晶體生長的品質，因此調整加熱功率及控制 適合晶體生長的熱場環境是晶體生長的重要參數。

2.5 統御方程式

( u) 0 t ρ _ρ→ ∂ _{+ ∇ ⋅ =} ∂ (2-7) 0 g p m R T ρ= _(2-8) 其中:

u

→ :坩堝內矽熔湯的流速 0 P :參考壓力 g R =8314 k kmol J ⋅ 動量方程式: 0 ( ) ( ) ( ) u u u p g j B t ρ→ → _ρ→ _{τ ρ ρ} → → → ∂ _{+ ⋅∇ = −∇ + ∇ ⋅ + − + ×} ∂ (2-9) 2 ( ) 3 j i ij eff eff ij j i u u u x x τ =μ ∂ +∂ − μ δ ∇ ⋅→ ∂ ∂ (2-10) ( ) j σ u B → → → = −∇Φ + × (2-11) t molecular eff μ μ μ = + _(2-12) ) ( → → × ∇ ⋅ = ΔΦ B u (2-13) ij

δ =1，for i=j ;δij=0，for i≠j (2-14)

其中:

B:磁通量密度

∫

∫

Ω Ω Ω = ∞ → π λ λ 4 0 ( , ) ) (r I r d d q_rad (2-16) t t eff Pr μ λ λ = + _(2-17) 擴散方程式: i i eff i i i S D u t ρ ϕ ϕ ϕ ϕ ρϕ + ∇ ⋅ ∇ = ⋅ ∇ + ∂ ∂ → ) ( ) ( ) ( , (2-18) P i i u i i S S S_ϕ = _φ +φ _φ (2-19) 其中: eff i D_ϕ, :有效擴散係數 紊流動能方程式(one-equation): __ t k j k j dk k P dt x x

ν

ν

ε

σ

⎡_{⎛ ⎞} ⎤ ∂ ∂ = _⎢⎜ + _{⎟ ⎥}+ − ∂ _{⎣⎝ ⎠}∂ _⎦ (2-20) 其中: ₂ 2 k t P =

ν

S :來源項 (2-21) 3/ 2 __ _k l F_ε

ε

= ⋅ :耗散項 (2-22) F_ε = f(Re )_{y (2-23)} Re_y k y

ν

= (2-24)

2.6 邊界條件

在系統內的邊界設定如下: 1.流場邊界條件: (1) 矽晶爐上方入口端為氣體進口(inlet)

(3) 爐體固體表面皆為無滑移邊界(no-slip condition) 2.熱場邊界條件: 中心軸設定為軸對稱；其矽晶爐內的矽熔湯及晶柱表面、石英坩堝、 石墨坩堝、固體元件表面等區塊的熱傳遞方式如下所示: (1)爐體表面皆設為絕熱(insulation)狀態 (2)固液界面處定溫 1685K，爐體外部定溫 300K (3)熱傳經固體元件表面邊界條件為 1 2 ( T) ( T) n n

λ

∂ =

λ

∂ ∂ ∂ (2-25) 1 2 ( )T =( )T (2-26) (4)熱傳經石墨坩堝與加熱器的邊界條件為 4 ( ) ( ) in

opaque transparent b rad w rad rad

T T T Q n n

λ

∂ =

λ

∂ +

σ ε

−

ε

∂ ∂ (2-27) 8 2 4 5.67 10 /( ) b W m K

σ

_{= ⋅} − _{:波茲曼常數} : rad

ε

放射率 : w T 壁面溫度 in rad Q :進輻射熱通量 (5)熱傳經矽熔湯或晶棒與空氣接觸表面和邊界條件為 1 2 1 2 ( ) ( ) ( in out) ( in out)

rad rad rad rad

T T Q Q Q Q n n

λ

∂ =

λ

∂ + − − − ∂ ∂ (2-28) 1,2 1,2 ( in ) in ( , ) rad k Q q r z d λ λ

λ

=

_∫

；( out)_1,2 out_1,2( , ) rad k Q q r z d λ λ

λ

≠∞ =

_∫

_(2-29)

圖 2-2 熔湯流動型態(a)垂直自然對流(b)水準自然對流(c)晶棒拉速影 響的強制對流(d)晶棒旋轉影響的強制對流(e)坩堝旋轉影響的強制對

第三章 數值方法

3.1 數值方法

本研究利用數值方法，求解矽晶爐內熱流變化，此數值方法使用 CGSim 長晶軟體為計算平台，計算紊流模式及不可壓縮流模式，解 二維軸對稱穩態Navier-Stokes equations，並使用 CGSim 內 Basic CGSim 模組利用有限元素法(finite-element)計算矽晶爐內全域熱場， 再由Flow Module 模組帶入 Basic CGSim 的熱場邊界條件利用有限體 積法(finite-volume)計算坩堝內矽熔湯區域的熱場與流場，物理變數定 義於每個控制容積中心做積分而求得，如圖 3-1，研究流程圖如圖 3-2。 將方程式(2-7)、(2-9)、(2-15)、(2-18)積分後得: 0 V S dV d S t

u

ρ

_ρ

→ → ∂ _{+ ⋅ =} ∂

∫

∫

(3-1) ( _i) i i ij j V S V S u dV u d S p I dV I d S t

u

ρ

_ρ

→ →

_τ

→ ∂ _{+ ⋅ = − ∇ ⋅ +} ∂

∫

∫

∫

∫

0 ( ) i i u S V V g dV j B dV S dV

ρ ρ

→ → +

_∫

− + ×

_∫

+

_∫

(3-2) ( _p ) p eff T V S S V c T dV C T d S T dS S dV t

u

ρ

ρ

→ →

λ

∂ + ⋅ = − ∇ ⋅ + ∂

∫

∫

∫

∫

(3-3) ( ) 0 S V d S B u dV → → → ∇Φ ⋅ + ⋅ ∇ × =

∫

∫

(3-4) ( ) i i i i eff i V S S V dV d S D d S S dV t

u

ϕ ϕ

ρϕ

_ρϕ

→ →

_ρ

_ϕ

→ ∂ _{+ ⋅ = − ∇ ⋅ +} ∂

∫

∫

∫

∫

(3-5) 有限體積法的守恆方程式（conservation equation）可表式如下：

(

ρφ

)

div

(

ρ φ

v

_φ

grad

( ))

φ

s

_φ

∂

₊

JJK

_{− Γ}

₌

(3-6)

數， 和 分別為擴散係數（diffusion）和來源係數（source）。 將(3-6)式積分後可得：

(

( ))

P j r j V S V

d

dV

v

grad

d S

s dV

dt

∫

ρφ

+

∑

∫

ρ φ

− Γ

φ

φ

⋅

=

∫

φ

JJG

JG

(3-7) T1 T2 T3 根據(3-7)式，對各項作離散： 1. 擴散項(diffusion terms) 2 ( r )j ( ( ) )j j j j j j j T ≈

∑

ρ φ

vJJG JG⋅S −

∑

Γ_φgrad

φ

⋅SJG ≡

∑

C −

∑

D (3-8) 擴散項用中點近似法可下列來表示： [ ( ) { ( ) ( ) } ] , l l D f grad S f grad _dPN j≈Γ

φ

j j N

φ

−

φ

P +

φ

⋅ − j

φ

⋅ j JG JG (3-9) 其中 為幾何因數， 為P 到 N 之向量， 為面擴散係數。 2. 來源項(source terms) 3 1 2 P

T

≈

s

−

s

φ

_(3-10) 3. 對流項(convective terms)

對流項以 Upwind Differencing Scheme 處理，根據(3-3)式的推 導，對流項可由下列表示：

,

F

0

P

j

插得到。總和上述，(3-7)式可表示為： 0 ( ) ( ) 0 n j V V F t ρ ρ δ − _{+ =}

∑

(3-13) 或 1 n o P P m m P P m A

φ

=

∑

A

φ

φ + +s B

φ

(3-14) 和 2 n P m m P m A ≡

∑

A

φ

+s +B (3-15)

3.2 固液界面計算

在矽熔湯凝固成矽晶棒之生長速度如下: 1 ( ) ( )

crys melt in out in out

crys crys melt rad rad crys rad rad melt x crys T _T V Q Q Q Q n

ρ

H

λ

n

λ

n ∂ ⎛ ∂ ⎞ = _⎜ − + − − − _⎟ Δ _⎝ ∂ ∂ _⎠ (3-16)

( )

_/ / , in inc

rad _{crys melt} crys melt

Q q d λ λ κ

λ

≠∞ =

_∫

(3-17)

( )

_/ / , out out

rad _{crys melt} crys melt

Q q d λ λ κ

λ

≠∞ =

_∫

(3-18) crys λ 、λmelt為矽晶棒與矽熔湯熱傳導係數，n 為沿固液界面的法向量。 相對生長速度: * relativ

crys crys crys

V =V −V (3-19)

*

crys

V

為固液界面上的平均生長速度

3.3 收斂條件

本研究計算溫場及流場判斷收斂條件是以下列兩項: 1. 每個時間步階(time step)的疊代次數 2. 每次疊代後所有控制容積內的殘值改變率小於給定值，即: k ₍ n n 0 0₎ P P P P C_φ =

∑

B

φ

− B

φ

<給定值

表(3-1) 材料參數 材料 狀態 性質 參數 單位 Conductivity 0.01+2.5×10−5T W/m．K Heat capacity 521 J/kg．K Molar 40 g/mol Dynamic viscosity T 8 6 _{5.365 10} 10 46 . 8 × − + × − 2 12 10 682 . 8 × − T − Pa．S Emissivity 0.8 氬 氣體 Crucible pressure 4000 Pa Conductivity Conductivity 0.874, T≦300K 0.2025, 300K≦T≦1000K 0.3690, 1000K≦T≦1500K 0.9638, 1500K≦T≦2300K W/m．K W/m．K 鐵 固體 Emissivity 0.8

Conductivity Graphite function W/m．K

石墨 固體 Emissivity 0.8 Conductivity 10 W/m．K 鉬 固體 Emissivity 0.18 Conductivity 4 W/m．K Emissivity 0.85 Heat capacity 1232 J/kg．K 石英 固體 Density(RHO) 2650 Kg/m3 Conductivity 15 W/m．K 鋼 固體 Emissivity 0.45 Conductivity 2 5 10 89 . 2 0943 . 0 89 . 98 − T + × − T W/m．K Emissivity 0.9016-0.00026208T Heat capacity 1000 J/kg．K 矽單晶 晶體 Density(RHO) 2330 Kg/m3 Conductivity 132 W/m．K Emissivity 0.3 Heat capacity 860 J/kg．K 矽 液體

Latent heat 1800000 J/kg

Surface tension 0.7835 N/m

Wetting angle 11 deg

Dynamic viscosity 4 10 225 . 8 _× − Pa．S

由Baic CGSim建立幾 何尺寸 定義邊界條件 設定網格 資料收集與處 理 設定演算法、控制參數 與相關設定 是否有錯誤或 計算時產生錯 誤 儲存計算 結果 後處理分 析計算結 果 結束模擬 計算 是 否 擷取熱場尺寸，由 Flow Module建立幾 何尺寸 設定網格 Basic GSim熱場 邊界條件 設定演算法、控制參數 與相關設定 由Baic CGSim建立幾 何尺寸 定義邊界條件 設定網格 資料收集與處 理 設定演算法、控制參數 與相關設定 是否有錯誤或 計算時產生錯 誤 儲存計算 結果 後處理分 析計算結 果 結束模擬 計算 是 否 擷取熱場尺寸，由 Flow Module建立幾 何尺寸 設定網格 Basic GSim熱場 邊界條件 設定演算法、控制參數 與相關設定 由Baic CGSim建立幾 何尺寸 定義邊界條件 設定網格 資料收集與處 理 設定演算法、控制參數 與相關設定 由Baic CGSim建立幾 何尺寸 定義邊界條件 定義邊界條件 設定網格 設定網格 資料收集與處 理 資料收集與處 理 設定演算法、控制參數 與相關設定 設定演算法、控制參數 與相關設定 是否有錯誤或 計算時產生錯 誤 儲存計算 結果 後處理分 析計算結 果 結束模擬 計算 是 否 是否有錯誤或 計算時產生錯 誤 是否有錯誤或 計算時產生錯 誤 儲存計算 結果 後處理分 析計算結 果 結束模擬 計算 儲存計算 結果 儲存計算 結果 後處理分 析計算結 果 後處理分 析計算結 果 結束模擬 計算 結束模擬 計算 是 否 擷取熱場尺寸，由 Flow Module建立幾 何尺寸 設定網格 Basic GSim熱場 邊界條件 設定演算法、控制參數 與相關設定 擷取熱場尺寸，由 Flow Module建立幾 何尺寸 擷取熱場尺寸，由 Flow Module建立幾 何尺寸 設定網格 設定網格 Basic GSim熱場 邊界條件 Basic GSim熱場 邊界條件 設定演算法、控制參數 與相關設定 設定演算法、控制參數 與相關設定 圖3-2 研究流程圖

第四章 結果與討論

本研究利用商用軟體 CGSim (版本 9.02)進行計算，其中 Basic CGSim 模組使用有限元素法(finite element method)針對整個矽晶爐 2 維軸對稱穩態全域(Global domain)進行計算，而 Flow 模組使用有限 體積法(finite volume method)對矽晶爐內矽晶棒、矽熔湯、石英坩堝、

石墨坩堝、輻射隔絕器等五個部份做為2 維軸對稱穩態區域(Local

domain)計算。方程式使用 2 階差分(second order approximations)針對 2 維軸對稱穩態全域進行離散。使用 3 階差分和 2 階差分針對統御方 程式(governing equations)中的對流項與傳導項進行離散。2 維穩態全 域熱傳(Basic CGSim 模組)的計算結果，可當成 2 維穩態軸對稱流動 模組中(Flow 模組)的熱邊界條件。收斂條件的相對殘值為 10-4。 本模擬之矽晶爐是使用 38 英吋坩堝，晶棒直徑為 450 公厘在晶 身長度為 1.34 公尺之尺寸進行穩態軸對稱模擬；長晶拉速固定為每 分鐘 0.5 公厘；坩堝和晶棒轉速分別為 3 rpm (順時鐘) 和 5 rpm (逆 時鐘)。爐體壓力為 0.024 bar，其氬氣流量為 355 slpm。 本研究針對大尺寸晶棒(直徑 450 公厘)之矽晶爐溫度分佈、流場 固液界面形狀、V/G 值等製程參數做深入研究並比較分析模擬的結 果，提供未來在生長大尺寸晶棒之生長製程參數上建議。若有需要做 爐體熱流場的改善或加熱效率最佳化之研究，也可做為實際最佳化設 計之參考依據。

4.1 模式驗證

米，及坩堝轉速5rpm (逆時針)、晶棒轉速 3rpm (順時針)、拉升速度 固定每分鐘2 厘米下進行模擬，網格數為 8685，利用穩態軸對稱 Flow 模組搭配五種不同的紊流模式，分別為Chien k-ε模式，標準 k-ε模 式，Spalart - Almaras 模式，one equation 模式(Flow 模組)，one equation 模式(Basic CGSim 模組)，所計算出的固液界面高度與文獻[33]的實驗 數據和模擬結果進行比較，如圖4-1 所示。經由觀察後可以發現，由 Flow 模組之 Chien k-ε 模式下計算得到的趨勢與實驗數據相當吻合， 也再次應證與文獻 [33]學者所使用的 Chien k-ε 紊流模式有較佳的固 液界面。因此，本研究選擇此紊流模式進行計算，4.2 節(含)之後的 結果均以此模式為基礎。

4.2 網格測試

在全域計算之 Basic CGSim 模組使用 17823、35708、51195 三種 網格數進行測試，由圖 4-2 所示，網格數為 17823 之固液界面模擬結 果偏離了35708 與 51195 之固液界面曲線，且網格數為 35708 與 51195 之固液界面模擬結果相當接近，為了有效率的計算及合理誤差範圍內 模擬結果為目的，故 Basic CGSim 模組採用網格數為 35708。在區域 計算之Flow 模組使用 7355、12235、35740 三種網格數進行測試，由 圖4-3 所示，12235 與 35740 網格之固液界面較為完全重合，故 Flow 模組採用較少網格數為12235。圖 4-4 為 2 維軸對稱穩態全域熱傳計 算之網格示意圖與溫度分佈圖，網格數為35708;圖 4-5 為 2 維穩態軸 對稱Flow 模組網格示意圖與溫度分佈圖，網格數為 12235。

增加時，晶棒容易形成空隙型的點缺陷(vacancy)，反之則易形成插入 型點缺陷(interstitial)。這些點缺陷在晶棒的冷卻過程中會聚集形成所 謂的D-defect 及插排環(dislocation loop)，而影響晶片在 IC 製程上的 應用良率。因此本研究採用0.25、0.35、0.45 及 0.50 mm/min 之拉升 速度探討其結果之影響。

4.3.1 拉升速度對矽熔湯流場和溫度場影響

熱的傳輸靠著三種主要模式，亦即輻射(radiation)，對流(convection) 及傳導(conduction)。在 CZ 法裡，熔湯是藉著由因石墨加熱器由熱輻 射的方式加熱坩堝，經由熱傳導加熱矽熔湯，而熔湯內部的熱傳輸則 是主要靠著對流，晶棒內部的熱傳輸是主要靠著傳導。由圖 4-6(a)至 圖4-9(a)溫度場所示，因石墨加熱器加熱矽熔湯，使得靠近石英坩堝 外側的矽熔湯溫度較熔湯中心溫度高，矽熔湯底部溫度較上方高，拉 升速度0.25、0.35、0.45 mm/min 之溫度梯度方向主要朝左上方(高溫 至低溫)，當拉升速度增加到 0.50mm/min，溫度梯度變化為向水平方 向，與拉升速度0.25、0.35、0.45 mm/min 比較，高溫區由矽熔湯右 下方上移至外側，矽熔湯下方的低溫區範圍增大。 由圖4-6(b) 至圖 4-9(b) 流場所示，拉升速度由 0.25 至 0.50 mm/min，矽熔湯內流場皆由三個渦旋(vortex)構成，可以發現矽熔湯 中心有一最大且較強烈的渦旋主導流場，在靠近坩堝壁附近則有較弱 的渦旋，而中間渦旋過渡在強烈渦漩與較弱渦旋之間。在坩堝中心 處，矽熔湯因為自然對流的關係而往下流動且其速度較快，相對來 說，坩堝壁面附近的矽熔湯流速較緩慢。當拉升速度由0.25mm/min、 0.35mm/min、0.45mm/min 增加至 0.5mm/min 時，因晶棒下方的低溫

4.3.2 拉升速度對固液界面與 V/G 參數之影響

在晶棒生長過程中固液界面之形狀極為重要，因晶體生長時的固 液界面形狀，和晶軸方向及溫度梯度有關。但在實際的長晶裡，由於 熱環境不斷在變化，很難去維持長時間的平面狀界面。在定性上(假 設晶棒的熱輻射、矽熔湯液面的熱輻射及晶棒拉升速度固定)，固液 界面形狀可由考慮熱平衡的關係式去了解，如方程式(3-16)中晶棒中 熱傳導通量大於矽熔湯之熱傳導通量時，固液界面形狀變成凸狀 (convex)，同理，凹形(concave)固液界面意味著過多的熱源進入固液 界面，而提升等溫線(isotherm)的位置。所以當固液界面形狀不是平 面狀(planer)時，在凝固過程中會使晶棒內產生熱應力。當熱應力小 於彈性應力時，這些熱應力在晶棒冷卻過程中消失掉。若固液界面形 狀太凸或太凹時，熱應力可能大於彈性應力，於是晶棒內便會產生差 排(dislocation)而失去結晶的完美性。 所以較平坦的固液界面被視為改善製程參數的首要目標，要產生 平面形的固液界面就必須由矽熔湯凝固成晶棒的徑向生長速度與晶 軸旋轉上升的拉升速度相近，故可藉由調整拉升速度來配合生長速度 之變化。由圖4-10 所示，拉升速度為 0.25mm/min、0.35mm/min 時， 因拉升速度小於生長速度，使其固液界面產生凸形，因凸形的固液界 面與水平面下方較高溫的矽熔湯容易產生反應使晶棒發生差排，凸形 與凹形固液界面相比下，凸形是製程中最不願意樂見的形狀，因此拉

晶棒拉速與晶棒軸向溫度梯度的比值為 V/G 值，此臨界值約 1.5 x 10-3 cm2 K-1 min-1左右;當 V/G 值高過臨界值，晶棒容易產生空隙型點 缺陷(low self-vacancies concentration)，當 V/G 值低於臨界值，晶棒容 易產生插入型點缺陷(self-interstitials concentration)，對於此兩種缺

陷，若矽晶圓中發生空位型點缺陷對於IC 元件的不良影響遠大於插

入型點缺陷，故V/G 值需大於且接近臨界值會有較好的矽晶棒品質。

由圖4-11 所示，拉升速度為 0.25 mm/min、0.35 mm/min、0.45 mm/min 在晶棒外側的V/G 值低於臨界值，而拉升速度為 0.50 mm/min，V/G 值均大於臨界值，故拉升速度為0.50 mm/min 時有較佳的晶棒品質。

4.4. 改變 Cusp 磁場強度之模擬

在矽晶生長中，熔湯裡的熱對流對晶棒的品質影響很大，特別是 一些不穩定的熱對流，容易導致固液界面的不穩定性及雜質分析的微 觀偏析現象(growth striation)。由於大部分的金屬及半導體熔液，都具 有高導磁性，所以在施以適當的磁場之下，可以抑制自然對流的程 度，避免紊流的發生。

MCZ 矽單晶生長法(Magnetic Field Applied Czochralski Method) 法，依據磁場方向來區分有橫型(horizontal)、縱型(vertical)及 cusp 型 三種，因橫型與縱型磁場使得晶棒雜質有分佈不均及石英坩堝溶解量 之優缺點，而使用Cusp 磁場如圖 4-12 Cusp 磁場示意圖所示，由上 下兩電磁線圈組成，兩磁線圈分別施予反向之磁場，使得在固液界面 處，因磁場縱向分量上下抵銷，使得晶棒是在磁場徑向分量為零的狀 況下生長，所以在雜質濃度在晶棒的徑向分佈均勻度得以維持，靠近 晶棒下方的熔湯是處於低強度的磁場中，所以該處熔湯仍能均勻攪

因此本研究使用Cusp 磁場來探討磁場強度與未開起磁場之溫度場與 流場及液固界面之好壞，圖4-13 為電磁線圈與與矽晶爐之相對位置 圖，其線圈位置參考文獻 [15]。

4.4.1 Cusp 磁場強度對矽熔湯流場和溫度場影響

本研究在石英坩堝與矽熔湯液面的接觸點設為滿足磁場強度之 參考點(r=48.5cm,dy=0cm)，分別以磁場強度 0.00T(未開磁場)、0.05T、 0.10T(開啟磁場)下，進行溫度場、流場與液固界面之整體探討。當考 慮Cusp 磁場效應時，溫度場如圖 4-14(a)至圖 4-16(a)所示，在未開磁

場時，矽熔湯中心有較大的低溫區，當磁場強度為0.05T 時，低溫區 漸漸的縮小，當磁場強度為0.10T 時，低溫區域縮至最小。而流場流 動形式如圖4-14(b)至 4-16(b)所示，在未開啟磁場時，矽熔湯流場由 三個較複雜的渦旋所構成流動，其中以晶棒正下方的較強烈渦流所主 導整個流場，當開啟磁場強度為0.05T、0.10T 時，使得流場與未開 啟磁場之相比下有明顯的差異，施予Cusp 磁場後在坩堝下方以垂直 磁場為主，抑制水平方向流動，故流動以上下為主。在熔湯液面以水 平磁場為主，抑制垂直方向流動，故流動以水平為主，所以矽熔湯流 場就變為一個較簡單的渦流為主的流動形式，這樣的流動形式可以減 少一些不穩定的熱流，在晶棒底部的熔湯往下流動的情形依舊可見， 但是其流動速度因為施加Cusp 磁場所產生的羅倫茲力 (Lorentz force) 抑制自然對流程度之關係而減緩，因開啟磁場也使得較快速的流動區 域由坩堝中心整個區域縮小至晶棒正下方區域。因此，較慢的流速降

在不同磁場強度下對固液界面的形狀如圖4-17 所示。由磁場強度 0.05T 至 0.10T，因為增加 cusp 磁場強度，減緩了矽熔湯內的自然對 流強度，因而縮小晶棒正下方的低溫區，減小了此區域軸向溫度梯 度，進而增加了生長速度，使開啟磁場後的固液界面較為開啟磁場之 固液界面為平坦，平坦的固液界面對於晶棒品質有較佳的穩定性。雖 然開啟磁場後的固液界面降低了界面的峰值，但卻增加了曲率半徑， 而曲率半徑對晶棒品質好壞目前沒有相關文獻或實驗數據做參考，關 於這部分需做更深入的探討。 開啟磁場與未開啟磁場之V/G 參數由圖 4-18 所示，未開啟磁場 之V/G 值在徑向上都有超過臨界值，而開啟磁場之 V/G 值在半徑大 於0.19 公尺卻小於臨界值，容易造成較不理想的插入型點缺陷，故 開啟磁場可效降低自然對流即產生較平坦的固液界面，但V/G 值卻 較未開磁場來得低。

4.5 改變 Cusp 磁線圈位置之模擬

許多關於使用磁場影響矽熔湯流動與溫度分佈進而預測晶棒品 質的穩態與非穩態研究相當完整，不過在Cusp 磁場中不同高度位置 之磁線圈(magnet coils)對矽熔湯流場與溫度場的影響尚未被進一步 研究且實際生長情況Cusp 磁場中心也未剛好在矽熔湯液面上，故本 研究藉由模擬不同磁線圈位置如圖4-19 所示，在磁場強度 B=0.05T 下，分別將線圈上移0、5 和 10 公分(dy=0cm、5cm 和 10cm)與下移 0、 5、10 公分(dy=0cm、-5cm 和-10cm)分別進行探討矽熔湯之流場、溫 度場、液固界面與V/G 值之好壞，提供較好的磁線圈位置供製程參 考使用。

4.5.1 Cusp 磁線圈位置對矽熔湯流場和溫度場影響

圖4-20 至圖 4-22 表示在 B = 0.05 T 下，磁線圈位置為上移 0、5 和10 公分 (dy = 0 cm、5 cm 和 10 cm) 時對熔湯流場與溫度場的影 響。可以觀察到，晶棒下方的熔湯流速隨著磁線圈高度增加而降低。 在dy = 0 和 5 公分時，有一較大渦旋存在於熔湯內，對於 dy = 10 公 分的情況下，熔湯內存在著一較小與一較大的渦旋。因此，由此結果 可以推測出當磁線圈往上移動越多，會導致熔湯流場的不穩定。對於 溫度場而言，由於羅倫茲力增加使得低溫區範圍朝向固液界面的方向 收縮，導致坩堝底部的溫度上升，同時也導致坩堝中心處的溫度梯度 增加。圖4-20、圖 4-23、圖 4-24 表示在 B = 0.05 T 下，磁線圈位置 為下移 0、5 和 10 公分 (dy = 0 cm、-5 cm 和-10 cm) 時對熔湯流場 與溫度場的影響。由圖可以得知，三種情況下均產生一較大的渦旋。 磁線圈往下移動越多，在坩堝中心處的速度也逐漸增加。對於溫度場 來說，三種情況下之低溫區的範圍是差不多的。由於向上之羅倫茲力 的影響，使得固液界面下附近的低溫區寬度變窄。此外，也可發現坩 堝壁面上的高溫區變薄。

4.5.2 Cusp 磁線圈位置對固液界面變化與 V/G 參數

在 B = 0.05 T 下，不同高度的磁線圈對固液界面的影響如圖 4-25 所示。從圖4-25 中可以發現，對所有情形而言，固液界面變化趨勢 均呈現往晶棒方向突起的現象。dy = 10 公分的情形使得固液界面的 變化最陡峭。dy = -10 公分的情形下可發現固液界面有鞍形 (saddle

溫度梯度的比值為V/G 值，此值約 1.5 x 10-3 cm2 K-1 min-1左右；通常 藉由此值可以判斷出矽晶棒是否有空隙型點缺陷 (low self-vacancies concentration) 的情形。從圖 4-26 可以發現到，dy = 0、5、-5 和-10 公分下的V/G 值很接近理想值，但是，dy = 0、-5 和-10 公分下的固 液界面差異的變化均大於dy = 5 公分下的結果。因此，依據本研究之 條件且經由比較固液界面差異變化以及V/G 值，本文建議將磁線圈 往上移動5 公分可得到較佳固液界面與 V/G 值。 此外，由圖 4-26 也可以再次得知，dy = 0 時的 V/G 值最大且移動磁線圈的高度可以降 低V/G 值。因此，固液界面的變化以及 V/G 均可以透過移動磁線圈 的高度而降低。

圖4-1 模式驗證比較圖

(a)

(b)

(a)

(b)

(a)

(b)

(a)

(b)

圖4-10 不同拉升速度下之固液界面

(a)

(b)

(a)

(b)

(a)

(b)

(a)

(b)

(a)

(b)

(a)

(b)

(a)

(b)

(a)

(b)

圖 4-25 不同高度的磁線圈位置之固液界面

第五章 結論

本研究利用CGsim 模擬軟體成功的模擬晶身長度和直徑分別為 500 公厘和 450 公厘的幾何尺寸，變化拉升速度、磁場強度及磁線圈 的位置來探討對38 英吋石英坩堝內之矽熔湯流場與溫度場的影響。 經由晶身長度為100 公厘之模擬與實驗數據比較後，發現 Chien k-ε 紊流模式較接近實驗數據，故本文採用 Chien k-ε 紊流模式進行計算。 首先進行模擬改變拉升速度對矽熔湯內流場與溫度場的影響，當 降低拉升速度有助於晶棒下方流速的減緩及縮小低溫區範圍，進而降 低了自然對流的強度，但較緩慢的拉升速度會使得固液界面產生凹面 (concave)，容易造成晶棒產生插排(dislocation)的發生，故調整較佳的 拉升速度影響著晶棒品質的好壞。 接著模擬考慮cusp 磁場效應時，增加磁場使得在流場方面晶棒底 部矽熔湯流速有明顯減緩，而整體流場由未開磁場之三個渦旋的複雜 流動形式轉變成單一渦旋主導整體流場的型態，而增加磁場強度也有 助於效降低液固界面高低差。 當磁線圈往上移動較多時會導致矽熔湯的流速減慢、坩堝中心處 的溫度梯度上升、坩堝底部的溫度升高以及矽熔湯流場的不穩定。低 溫區範圍朝向固液界面的方向收縮。此外，磁線圈往下移動對低溫區 的範圍影響不明顯。本文也針對固液界面的變化與V/G 值進行探討。 當考慮磁場效應且不移動磁線圈的情況下，固液界面差異變化以及

參考文獻

[1] G..K. Teal and J.B. Little, “Growth of Germanium Single Crystals Containing p-n Junctions,” Phys. Rev., 78, 647(1950).

[2] J.Z. Czochralski, Phys. Chem., 92, 219(1918).

[3] H.E. Buckley, Crystal Growth., John Wiley and Sons Inc., New York, (1951). [4] H. J. Scheel and T. Fukuda, “The Development of Crystal Growth Technology,” Crystal Growth Technology., 3-14(2003).

[5] W.C. Dash, J. of Applied Physics, 30, 4, 459(1959).

[6] J.J. Derby and R.A. Brown, “Thermal-capillary analysis of Czochralski and liquid encapsulated crystal growth I . Simulation,” Journal of Crystal Growth, 74,

605-624(1986).

[7] J.J. Derby and R.A. Brown, “Thermal-capillary analysis of Czochralski and liquid encapsulated crystal growth II. Processing strategies,” Journal of Crystal Growth, 75, 227-240(1986).

[8] F. Dupret, et al. “Numerical calculation of the global heat transfer in a Czochralski furnace,” Journal of Crystal Growth, 79, 84-91(1986).

[9] N. Van den Bogaert and F. Dupret, “Dynamic global simulation of the Czochralski process I. Principles of the method,” Journal of Crystal Growth, 171, 65-76(1997). [10] N. Van den Bogaert and F. Dupret, “Dynamic global simulation of the

Czochralski process II. Analysis of the growth of a germanium crystal,” Journal of Crystal Growth, 171, 77-93(1997).

[11] Xiao-Bo Wu and Xu Geng, ”Fundamental study of crystal/melt interface shape change in Czochralski crystal growth,” Journal of Crystal Growth, 169,

786-794(1996).

[12] A. Lipchin and R. A. Brown, “Hybrid finite-volume/finite-element simulation of heat transfer and melt turbulence in Czochralski crystal growth of silicon,” Journal of Crystal Growth, 216, 192-203(2000).

[13] K. Takano, et al., “Global simulation of the CZ silicon crystal growth up to 400mm in diameter,” Journal of Crystal Growth, 229, 26-30(2001).

[14] Y.R. Li, et al., “Global simulation of silicon Czochralski furnace in an axial magnetic field,” Journal of Heat and Mass Transfer, 46, 2887-2898(2003).

[15] Shiraishi, Y., et al., “Growth of silicon crystal with a diameter of 400 mm and weight of 400 kg,” Journal of Crystal Growth, 229, 17-21(2001).

[16] Gorbunov, et al., “Melt flow instability and vortex structures in Czochralski growth under steady magnetic fields,” Energy Conversion and Management, 43,

130-133(2000).

[18] Klein, et al., “Numerical simulations of the influence of a traveling magnetic field”, Magnetohydrodynamics 45, 557-567(2009).

[19] Krauze, A., et al., “Unsteady 3D LES modeling of turbulent melt flow with AC travelling EM fields”, Magnetohydrodynamics 45, 605-611(2009).

[20] Klein, et al., “Numerical simulation of Czochralski crystal growth under the influence of a traveling magnetic field generated by an internal heater-magnet module (HMM)”, Journal of Crystal Growth, 310, 1523-1532(2001).

[21] Rudolph, P., “Travelling magnetic fields applied to bulk crystal growth from the melt: The step from basic research to industrial scale”, Journal of Crystal Growth, 310, 1298-1306(2008).

[22] Lechner, C., et al., “Development of a software for the numerical simulation of VCz growth under the influence of a traveling magnetic field”, Journal of Crystal Growth, 303, 161-164(2008).

[23] Kalaev, V. V., “Combined effect of DC magnetic fields and free surface stresses on the melt flow and crystallization front formation during 400 mm diameter Si Cz crystal growth”, Journal of Crystal Growth, 303, 203-210(2007).

[24] Kishida, Y., et al., “Geostrophic turbulence in CZ silicon melt under CUSP magnetic field”, Journal of Crystal Growth, 273, 329-339(2005).

[25] Munakata, T., et al., “Effect of high frequency magnetic field on CZ silicon melt convection”, International Journal of Heat and Mass Transfer, 47, 4525-4533(2004). [26] Shiraishi, et al., “Prediction of solid–liquid interface shape during CZ Si crystal growth using experimental and global simulation”, Journal of Crystal Growth, 266, 28-33(2004).

[27] Krauze, A., et al, “Numerical 2D modelling of turbulent melt flow in CZ system with dynamic magnetic fields”, Journal of Crystal Growth, 266, 40-47(2004).

[28] Ivanov, N. G., et al, “Analysis of magnetic field effect on 3D melt flow in CZ Si growth”, Journal of Crystal Growth, 250, 183-188(2003).

[29] Vizman, D., et al., “3D numerical simulation and experimental investigations of melt flow in an Si Czochralski melt under the influence of a cusp-magnetic field”, Journal of Crystal Growth, 236, 545-550(2002).

AC and CUSP magnetic fields and its verification in a laboratory facility”, Journal of Crystal Growth 230, 81-91(2001).

[33] Lukanin, D. P., et al., “Advances in the simulation of heat transfer and prediction of the melt-crystal interface shape in silicon CZ growth”, Journal of Crystal Growth, 266, 20-27 (2004) .