張 晶,劉 丁,杜燕軍

(西安理工大學晶體生長設備及系統(tǒng)集成國家地方聯(lián)合工程中心,陜西西安 710048)

1 引言

直拉法生長直徑均勻的硅單晶是目前半導體芯片材料制造的關鍵技術[1-2],也是衡量晶體品質的重要指標之一.在硅單晶生長過程中,通過提拉運動在相變界面(即熔體-晶體交界面)處形成一定的溫度差,依據晶體生長原理[3],相變界面溫度場的均勻性決定了該溫度差的穩(wěn)定性,是硅單晶均勻生長的根本保證.然而,相變溫度場受提拉速度、加熱器溫度、晶體旋轉和坩堝旋轉等工藝參量的影響且隨著熱能、動能和質量輸運的狀態(tài)而變化,呈現非均勻性特征.

當前,針對相變溫度場非均勻性問題,在晶體生長溫度建模的研究中,在表征固液界面徑向溫度梯度和軸向溫度梯度共同作用規(guī)律時,往往忽略了提拉速度動力學因素[4-6].然而,由拉速引起的結晶潛熱所引發(fā)的熱補充對固液界面軸向溫度梯度有著不可忽視的影響.Brown[7]和Lan[8]分別給出了詳細的晶體生長過程建模和動力學分析方法,合理的調整晶轉和堝轉工藝參數,能夠提高固液界面溫度的均勻性.然而,隨著硅單晶生長直徑增大,熱場尺寸增加,液面位置變化等因素仍對相變溫度場的均勻性產生影響,從而引起晶體直徑波動.與此同時,熔體內部熱傳輸機理復雜,晶體生長直徑控制系統(tǒng)所具有的強非線性、大時滯、不確定性特征,給控制器設計帶來困難[9-10].因此,從晶體內部熱傳輸和溫度分布入手研究其直徑控制問題是一有效途徑.Derby等人[11]通過研究質量和能量平衡關系建立了晶體的熱傳輸模型,該模型為時變空間域上的偏微分方程.Wang等人[12-13]研究了基于穩(wěn)定溫度場的最優(yōu)控制器,提出了時變空和間域的分布式系統(tǒng)的穩(wěn)定與控制模型.然而,該模型描述的晶體溫度隨時空變換且具有無窮個自由度,本質上是無限維的,難以獲得解析解.同時,受實際生產條件限制,生長溫度難以實時測量,大大增加了直徑控制難度.

針對這一問題,對無限維分布參數系統(tǒng)進行降維是分析、優(yōu)化和控制生長直徑的有效方法[14].傳統(tǒng)的降維方法如差分法和有限元法只能得到階數非常高的近似模型,并不適用于分布參數系統(tǒng)的快速仿真和控制器設計[15-17].因此,研究合理的降維方法來獲得較為準確的溫度場模型是十分重要的.針對上述問題,本文建立一種基于偏微分方程(partial differential equation,PDE)模型的溫度場優(yōu)化模型,該模型采用二階常微分方程(ordinary differential equation,ODE)描述了機械拉動臂引起的邊界演化,考慮生長速率波動的時變空間域上的拋物型PDE來描述晶體溫度的變化.基于PDE無限維空間上晶體熱傳輸的時變模型,采用譜方法和Galerkin方法獲得有限維近似模型,以此來簡化溫度調節(jié)和域運動控制問題的數值實現.采用線性二次型對改進后的溫度模型進行控制,對比分析了相變界面溫度分布和晶體生長直徑,從而驗證方法的有效性.

2 非均勻相變溫度場問題描述

2.1 問題分析

依據點缺陷在固液界面處的相變溫度場形成理論,相變溫度場的某時空臨界尺度由比值決定,V為晶體的提拉速度,G為固液界面靠近晶體側的軸向溫度梯度.當比值維持在0.134(1±10%)mm2/?C·min范圍時,晶體可在穩(wěn)定的相變溫度場中生長,形成無缺陷的完美硅單晶[18].當相變溫度場不穩(wěn)定波動超過理想范圍時,會引起晶體生長速率波動并反映在晶體等徑生長時直徑的變化上.實際中往往通過提拉速度來調控固液界面相變溫度場的均衡性,這是由于直徑能夠快速響應拉速調節(jié),而通過加熱器溫度調節(jié)晶體直徑的過程中,由加熱器傳至固液界面的熱傳輸通道存在較大滯后,因而調節(jié)響應慢.因此,由生長速率波動引起的提拉速度的波動是導致晶體直徑波動的根本原因.

本質上,相變溫度場不穩(wěn)定波動一方面來自提拉速度對溫度分布影響,即晶體不斷生長使得晶體邊界發(fā)生變化,影響固液界面處的溫度分布.另一方面熔體液位不斷下降導致熱輻射增強,影響了固液界面處溫度分布的均勻性.因此,建立晶體生長過程提拉動力學模型與熱傳輸模型時,需要進一步考慮晶體生長和熔體液位下降對相變溫度場的影響,從而改善相變溫度場的非均勻性.

2.2 提拉動力學模型

根據晶體生長原理,晶體生長過程滿足質量守恒定律和牛頓第二定律,如圖1所示.

圖1 直拉法晶體生長過程示意圖Fig.1 Schematic diagram of the Czochralski crystal growth progress

文獻[19]中假設固液界面平直,提出基于晶體生長提拉過程的動力學模型,該模型能夠很好的描述提拉作用對晶體直徑和晶體生長速率的影響.但是,在晶體生長過程中,固液界面處熱通量并不是一成不變的,會隨著溫度場的變化而發(fā)生改變,從而影響晶體生長直徑.由第2.1節(jié)分析可知,液面位置下降所引起的相變溫度場非均勻性是改變固液界面形狀的根本原因,如圖1所示,固液界面形狀凸向晶體.針對這一問題,通過建立考慮坩堝上升速率的提拉動力學模型,以保證液面位置的穩(wěn)定來改善相變溫度場的非均勻性,從而避免生長速率的波動.

在晶體生長過程中,根據質量守恒定律可得Mt=Ml+Mm+Ms,Mt為總投料量,Ml為熔體質量,Mm為彎月面的質量,Ms為晶體質量.假設在晶體生長過程中Mm恒定，則熔體減少的質量等于晶體增加的質量,即

其中k為堝跟比.

直拉法硅單晶生長過程中,相變界面處能量的傳遞和變化均直接影響界面溫度,從而影響生長界面溫度差.圖1給出了晶體生長過程中固液界面的熱量傳遞過程.根據固液界面的能量守恒方程,可得晶體生長速度與熱量的關系為

其中:Φs表示穿過固液界面進入晶體的熱量,Φl表示從熔體進入固液界面的熱量,Φh表示晶體結晶時釋放的潛熱,?H為晶體結晶時釋放的潛熱.上式表明,在單位時間內相變界面?zhèn)鬟f到晶體中的熱量等于通過熔體傳遞至相變界面的熱量和相變過程中釋放的潛熱之和.因此,希望相變界面處有穩(wěn)定熱量傳輸就意味著Φs?Φl在晶體生長過程中保持穩(wěn)定.

根據牛頓第二定律可得

聯(lián)立式(1)-(5),可得提拉動力學模型

2.3 溫度場優(yōu)化模型

根據提拉動力學模型得到了晶體邊界演化和晶體的生長速率,通過對晶體生長動力學特性與熱傳輸方式的分析,式(8)描述了在時變空間域上對流擴散過程的熱傳輸過程[16]:

其中:x(r,z,t)表示晶體的溫度,Pe為Peclet常數,kr表示熱傳導率,V表示晶體生長速度.忽略徑向生長速率,并通過適當的坐標變換,在柱坐標系下式(8)可以表示為

其中:k0=,Vz(t)表示晶體沿軸向的生長速度,A(r,z,t)為對x(r,z,t)的操作的算子.

如圖1所示,加熱器u(t)作用在晶體的邊界上,并且在晶體生長過程中,加熱器位置保持不變.晶體在固液界面處結晶,固液界面處的溫度恒為晶體的熔點xf,且在r=0處無熱量交換,并假設在其余邊界處熱通量為零,邊界條件可表示為

其中:r(t)表示晶體的半徑,l(t)表示晶體的高度,u(t)表示加熱器輸入.

由式(8)可知,晶體的溫度是相對時間和空間位置的函數,具有無窮維自由度,且往往不能通過有限個狀態(tài)進行求解.因此,需要對該模型進行降維.

3 晶體生長溫度場降維

針對無限維分布參數系統(tǒng)的求解問題,采用有限維的常微分方程描述的系統(tǒng)來近似無限維分布參數系統(tǒng).首先需要確定熱傳輸作用的有效邊界,即晶體直徑和長度的變化.在晶體生長過程中,總是期望晶體直徑恒定,因此,控制晶體直徑尤其重要.本文在提拉動力學模型和直徑控制的基礎上確立熱傳輸邊界,進一步對熱傳輸模型進行降維.

3.1 熱傳輸邊界的確立

根據式(6)-(7)所示的晶體提拉動力學模型,對輸入輸出進行線性化,可定義Fext(t)為

3.2 PDE模型降維方法

基于上述方法確立的邊界條件對PDE模型進行降維.譜方法在整個空間域上選擇全局和正交的空間基函數,相比較于有限單元法的空間離散能夠得到維數更低的近似模型.譜方法適用于一類能夠進行快慢變量分離的系統(tǒng),快慢變量是指空間基函數在空間頻率域中的排列順序往往是由慢到快,由于快變量對系統(tǒng)的貢獻較小,因此,選擇慢變空間基函數即可滿足降維的要求.快慢變量的區(qū)分通過Galerkin截斷準則來判定.

譜方法對模型降維的基本步驟為

1) 選取空間基函數,將系統(tǒng)變量在空間基函數上展開.

空間微分算子的特征函數常用作空間基函數[14],對于算子A,其與特征值和特征函數之間關系為

式中:λ為特征值,Φ為特征函數,可通過式(14)對算子A分離變量進行求解[19].

將系統(tǒng)變量x(r,z,t)在空間基函數上展開,即

對于式(9)-(10)所描述的系統(tǒng),輸入控制u(t)加在邊界上,并假設輸入控制作用函數為

該函數代表輸入的作用范圍,其中ε1>0,通過采用Drica delta函數將邊界控制系統(tǒng)轉化為與其等價的分布控制系統(tǒng),則輸入算子B(t)可以定義為

結合式(9)和式(15)可以將熱傳輸模型寫為

2) 采用Galerkin方法得到描述展開系數的無限維常微分方程組.

Galerkin方法中選取與空間基函數相同的權重函數使系統(tǒng)殘差最小,獲得展開系數的無限維常微分方程組,記為

其中:mn=1,2,···,∞,λmn表示算子A的第mn個特征值.

3) 通過有限維截斷與空間基函數綜合得到系統(tǒng)的低維近似模型.

為了對無窮維常微分方程降維,需根據Galerkin截斷理論確定慢變空間基函數的個數.Galerkin截斷準則[20]如下.

對于分布參數系統(tǒng),Reλj表示的是系統(tǒng)空間算子特征值λj的實部,如果0 ≥Reλ1≥Reλ2≥···≥Reλj≥···,且存在k使之滿足

3.3 溫度控制

對于系統(tǒng)的最優(yōu)控制,線性二次型性能指標能夠更好的兼顧系統(tǒng)的狀態(tài)以及輸入能量消耗.因此,將線性二次型最優(yōu)控制應用于拋物PDE系統(tǒng)低維模型,結合晶體域運動來簡化晶體溫度調節(jié),得到系統(tǒng)的最優(yōu)輸入與溫度分布.設計一個最優(yōu)輸入u(t),使如下的性能指標函數達到最小值:

其中:Q和R為非負定對稱和正定對稱的系數矩陣,a為系統(tǒng)的狀態(tài).為使性能指標J最小,最優(yōu)輸入取為

其中P(t)是Riccati方程的解.

4 結果與分析

本次仿真實驗數據來自本中心12英寸硅單晶爐的實際拉晶試驗,記錄了整個晶體生長過程中直徑,生長速率和加熱器功率等數據,采樣時間為2 s.選取放肩和等徑階段10000個直徑數據作為晶體的參考直徑,其中前1500個數據為放肩階段,其余數據為等徑階段.仿真過程中各初始狀態(tài)取值參考實際晶體生長數據,如表1所示,且模型中堝跟比設置參考實際晶體生長中堝跟比.

表1 仿真過程參數選取Table 1 Parameter selection of simulation process

在晶體生長過程中,固液界面熱通量時刻變化.在不考慮坩堝上升時,C值的變化可以看作溫度場變化對固液界面熱通量的影響,如圖2中C1所示.在考慮坩堝上升時,由圖2中C2給出了在仿真階段C的變化曲線.

圖2 晶體生長過程固液界面C值變化Fig.2 Change of C at solid-liquid interface during crystal growth process

在晶體生長的不同階段,固液界面熱通量變化不同,反映在C值上.放肩階段是使晶體直徑平滑放大至預設直徑的過程.由圖2可知,在放肩前期,C值較高,隨著晶體進入等徑階段,C值基本保持恒定.原因是在放肩前期,生長界面距離晶棒頂部距離較小,大量的熱通過提拉桿流出,固液界面軸向溫度梯度較高,C值較大.隨著直徑不斷增大,在放肩階段液面變化微小可忽略不計,假設流入固液界面熱量不變,由式(4)可知,Φs逐漸減小,C值逐漸減小.在等徑階段,固液界面熱通量相對穩(wěn)定,晶體保持穩(wěn)定生長,符合晶體生長工藝要求,驗證了模型的正確性.由圖2可知,C2值波動較小,說明相變溫度場波動較小,在等徑階段維持液面位置不變有利于改善固液界面相變溫度場的均勻性.因此,在晶體生長過程中考慮坩堝上升速率能夠有效地減少固液界面下降引起的相變溫度場變化.

針對改進前和改進后兩種模型對晶體直徑控制進行仿真,仿真過程設置控制增益K=0.1,可得晶體生長速率變化曲線如圖3所示.

圖3 晶體生長速率Fig.3 Crystal growth rate

圖3給出了晶體在放肩和等徑階段生長速率的變化曲線.由C值分析可得,在放肩階段,提拉速度較高,生長速率近似于提拉速率,因此,在放肩階段晶體生長速率較高.由放肩至等徑的過程需要進行轉肩工藝,即需提高拉速抑制晶體直徑放大,如圖3所示,此處拉速明顯升高.將優(yōu)化后的生長速率進行對比,由圖3可知,改進前模型中C值波動較大,為了能夠獲得穩(wěn)定的直徑,須通過調節(jié)拉速來減小直徑波動,即導致晶體生長速率波動較大;而改進后模型中C值較為穩(wěn)定,能夠獲得穩(wěn)定的生長速率.

圖4給出了改進前和改進后模型直徑控制結果,并給出了與期望直徑之間的誤差.

由圖4可知,改進后模型能夠更好穩(wěn)定晶體的直徑.通過計算均方根誤差對結果進行精確評估,改進前模型求得直徑的均方根誤差為0.1636,改進后模型求得直徑的均方根誤差為0.0423,可以得出,考慮坩堝上升能夠減小相變溫度場變化,有利于實現晶體直徑控制.

圖4 晶體半徑變化及誤差對比Fig.4 Change of crystal radius and deviation comparison

為了進一步驗證液面位置對晶體生長過程的影響,在加熱溫度不變條件下對晶體直徑控制進行仿真,如圖5所示.在不考慮坩堝上升的情況下,由于液面位置下降,坩堝壁的裸露面積不斷增加,形成大量的熱輻射,使晶體生長熱量無法及時散出,即Φs減小,導致晶體直徑減小.實際中需降低拉速來進行補償,不但容易引起拉速的波動,同時降低了效率.由圖5可知,在考慮堝上升的情況下能夠獲得期望的晶體直徑.

圖5 兩種條件下晶體的生長半徑對比Fig.5 Comparison of the growth radius of crystals under two conditions

通過上述結果獲得了溫度模型邊界條件,將提拉動力學模型得到的晶體直徑與高度變化應用于偏微分方程算子A的特征值計算中,可得到特征值變化.表2為徑向特征函數對應的特征值,不隨時間變化.圖6為軸向特征函數對應的特征值,其隨時間變化,這里僅給出前5項特征值的變化.

表2 部分徑向特征值Table 2 Partial radial eigenvalues

圖6 部分軸向特征值Fig.6 Partial axial eigenvalues

根據Galerkin截斷準則可知,當算子A的特征值滿足截斷條件時,可以將分布參數系統(tǒng)進行快慢分離.一般,當=0.1,即λk十倍于λ1時,就可截斷前k階低維慢變量部分來近似該無限維系統(tǒng).由表2可得.徑向特征值滿足0 ≥Reλ1≥Reλ2≥···,

因此,截取到前4個特征值即可滿足降維要求.同理,根據圖6所示,軸向特征函數對應的特征值只需截取前4個即可滿足降維要求，即式(19)中M=4,N=4,因此,將分布參數系統(tǒng)模型降維到16維.

根據特征值的截斷數,在圖7-8分別給出了徑向和軸向前4個特征函數的變化曲線,其作為空間基函數來實現模型降維.從圖7-8可以看出,特征函數之間頻率變化明顯,具有良好的快慢分離特性,適于模型降維.

圖7 前4個徑向特征函數Fig.7 The first 4 radial eigenfunctions

圖8 前4個軸向特征函數Fig.8 The first 4 axial eigenfunctions

為了穩(wěn)定固液界面附近的溫度分布,對降維后的模型進行線性二次型最優(yōu)控制,圖9僅給出了前4維時間系數amn(t)的變化曲線:

最優(yōu)控制得到系統(tǒng)的輸入變化曲線如圖10所示.

圖11給出了降維后固液界面處晶體徑向溫度隨時間的變化.從中選取不同時刻的溫度分布加以分析,如圖12所示.圖12中初始時刻晶體徑向溫度存在波動,通過對系統(tǒng)進行最優(yōu)控制后,晶體固液界面附近溫度逐漸穩(wěn)定.

圖9 時間系數amn(t)的變化曲線Fig.9 Change curve of time coefficient amn(t)

圖10 輸入u(t)變化曲線Fig.10 Change curve of input u(t)

圖11中溫度差表示與凝固點溫度的差值,即可通過溫度差變化來表示溫度分布變化.從圖11中可以看出,固液界面處溫度分布隨時間增長而更加平穩(wěn).圖12中對模型改進前、后得到的溫度差分布比較可以得出,改進前與改進后模型中固液界面徑向溫度梯度的差值逐漸增大,是由于隨著晶體的不斷生長,液面位置不斷下降,使得相變界面的非均勻性增強,且隨著時間推移表現更為明顯.改進后模型能夠大幅度提高徑向溫度分布均勻性,得到的固液界面徑向溫度梯度較小,能夠獲得更加平坦的固液界面形狀.在線性二次型最優(yōu)控制作用下,改進前后模型的固液界面徑向溫度梯度在逐漸減小,最終獲得穩(wěn)定的相變界面溫度分布.在這一過程中,改進后模型溫度梯度低于改進前模型,進一步驗證了改進模型與降維方法的有效性.

圖11 固液界面處溫度分布Fig.11 Temperature distribution at solid-liquid interface

圖12 不同時刻溫度分布曲線Fig.12 Temperature distribution curve at different time

5 結論

直拉法晶體生長相變界面溫度分布非均勻性導致晶體生長直徑不均勻.針對這一問題本文提出了基于PDE模型的溫度場最優(yōu)控制策略.建立了考慮坩堝上升的改進提拉動力學模型,確定晶體域邊界的演化,從而使液位下降導致的相變溫度場非均勻性得到改善.采用譜方法選擇特征函數作為空間基函數,Galerkin將時間系數展開為無限維常微分方程組,并根據截斷準則得到低維近似模型,解決了PDE模型無法直接求解的問題.通過反饋線性化和線性二次型來求解域運動以及域運動單向耦合的溫度分布問題.仿真實驗結果表明,改進模型使得相變界面徑向溫度分布更加均勻,有利于直拉硅單晶穩(wěn)定生長.