CZ法硅單晶生長電阻網(wǎng)絡建模及MPC仿真分析

弋英民， 張　潼

(西安理工大學 自動化與信息工程學院, 陜西 西安 710048)

?

CZ法硅單晶生長電阻網(wǎng)絡建模及MPC仿真分析

弋英民， 張潼

(西安理工大學 自動化與信息工程學院, 陜西 西安 710048)

CZ法硅單晶生長系統(tǒng)具有強非線性及大滯后的特點。本文就該系統(tǒng)建立了傳熱傳質的低階模型，計算了該模型在初始條件下的穩(wěn)態(tài)解，并進行線性化處理及分析，根據(jù)系統(tǒng)的可控性及穩(wěn)定性確定控制策略，并設計控制器。結果表明：在低階模型中，采用Gevelber提出的簡化輻射角系數(shù)能較好的反映工業(yè)中晶體生長過程特性；系統(tǒng)線性化后的性能分析表明，僅將加熱器功率作為輸入時，系統(tǒng)是不可控不穩(wěn)定的，故采用模型預測控制算法作為控制策略，設計了無干擾無延遲的無約束模型預測控制器和無干擾有延遲的無約束模型預測控制器，仿真結果驗證了所設計的方法的有效性。

單晶硅; CZ法; 電阻網(wǎng)絡類比; 低階模型; 模型預測控制

單晶硅半導體硅材料是電子信息技術產業(yè)以及新能源產業(yè)的基礎材料，也是半導體器件和集成電路的重要原材料，因此在增大硅單晶直徑的同時，提高直拉單晶硅的質量仍是研究的主要內容[1]。而CZ法是工業(yè)化生產單晶硅的最主要方法，其生長過程中既有物質的傳輸，也有動量和熱量的傳輸[2]。坩堝內熔體的流動、傳熱及傳質現(xiàn)象直接影響單晶硅的質量。國內外主要從理論分析、實驗測定和數(shù)值模擬[3,4]三個方面對晶體生長過程進行研究[5]，Derby、Brown[6-7]、Hurle[8]等在晶體生長理論建模方面做出了大量貢獻；Brice[9]、Jones[10]、Seidl[11]等通過實驗研究分析了晶體生長過程中坩堝旋轉和雜質氧等因素對硅單晶質量的影響。但是實驗成本高、耗時長等原因制約了實驗測定方面的進一步研究。

本文對CZ法生長單晶硅過程中的等徑階段采用電路類比法建立基于傳熱傳質現(xiàn)象的低階模型，并以此作為系統(tǒng)分析和控制算法研究的基礎。

1　模　型

本文在熱量、質量傳遞基本定律的基礎上[12-15]，建立了具有兩個輸入量(Vp和Pin)、八個狀態(tài)量(Th、Tc、Tm、Ri、Hi、Re、Hm和φ)和兩個輸出量(Tm和Ri)的低階模型，根據(jù)加熱器、坩堝、熔體的輻射傳熱和傳導傳熱現(xiàn)象，采用電阻網(wǎng)絡類比的方法建立單晶硅生長系統(tǒng)的溫度狀態(tài)方程；主要在熔體-晶體分界面處，根據(jù)分界面的幾何形變建立幾何狀態(tài)方程。系統(tǒng)低階模型是一系列具有如下形式的非線性方程：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

式中，Vp是提拉速率，Pin是加熱器輸入功率，Th是加熱器溫度，Tc是坩堝溫度，Tm是熔體溫度，Ri是分界面半徑，Hi是分解面高度，Re是有效半徑，Hm是熔體高度，φ是接觸角。

2　建模的方法步驟

2.1建立晶體生長的溫度狀態(tài)方程

單晶爐內的熱量傳遞方向及單晶爐部分尺寸如圖1所示。溫度狀態(tài)方程涉及CZ法晶體生長系統(tǒng)中的大部分組分，即加熱器、坩堝和熔體的溫度。該狀態(tài)方程以輻射傳熱和傳導傳熱引起的能量平衡公式(9)為基礎。

Ein-Eout=Estored

(9)

式中，Ein是系統(tǒng)接收的能量，Eout是系統(tǒng)耗散的能量，Estored是系統(tǒng)存儲的能量。根據(jù)式(9)的系統(tǒng)能量平衡方程可以得到系統(tǒng)的溫度狀態(tài)方程。

對于加熱器：

(10)

對于坩堝：

(11)

對于熔體：

(12)

式中，Ch是加熱器熱容，qhc是加熱器傳遞至坩堝的傳熱速率；Cc是坩堝熱容，qco是坩堝傳遞到周圍環(huán)境的傳熱速率，qcs是坩堝傳遞到熔體表面的傳熱速率，qcm是坩堝傳遞到熔體內部的傳熱速率；Cm是熔體變化的熱容，qso是熔體表面到周圍環(huán)境的傳熱速率，qm是熔體到固液分界面區(qū)域底部的傳熱速率，Rcru是坩堝的半徑；Hcru為坩堝的高度。

圖1　傳熱率及系統(tǒng)尺寸Fig.1　Heat transfer rates and dimensions

2.1.1輻射角系數(shù)計算

單晶爐內主要涉及的輻射角系數(shù)有：坩堝到熔體表面的輻射角系數(shù)、坩堝到周圍環(huán)境的輻射角系數(shù)、熔體表面到周圍環(huán)境的輻射角系數(shù)。

1) 基本模型：圖2為輻射角系數(shù)的角度近似。根據(jù)圖2所示的坩堝內角度近似關系，可以得出輻射角系數(shù)的近似結果。圖中，γ、ω分別為坩堝表面、熔體表面之間的夾角，H為暴露的坩堝壁高度。

圖2　輻射角系數(shù)的角度近似Fig.2　View factor angle approximation

2) 輻射增強模型：對兩個表面A1和A2之間輻射角系數(shù)的定義式進行精確求解：

(13)

式中，R是兩個表面之間的距離；θ1和θ2分別表示兩個表面上的微小面積dA1、dA2間各自的天頂角。

2.1.2傳熱速率計算

將坩堝內熔體以上部分看作一個圓柱形外殼。下表面為熔體表面，即表面1；上表面為關于坩堝內部環(huán)境的假象表面，即表面2；側表面為暴露坩堝壁的表面區(qū)域，即表面3。

1) 基本模型：認為在換熱過程中一次僅有兩個表面相互接觸[8]，并進行熱量傳遞。

采用電阻類比法對接觸表面進行分析，即表面的傳熱速率qj、表面間的傳熱速率qjk類比為電流，輻射能量Ebj和能量Jj類似電壓，輻射角系數(shù)和對應面積的相關項記為表面輻射熱阻Rj和空間輻射熱阻Rjk(j=1、2、3，k=1、2、3)。通過電路求解得出相關傳熱速率。

2) 輻射增強模型：認為換熱過程中三個表面之間的傳熱同時進行，同樣采用電阻類比法對接觸表面進行分析，如圖3所示。在節(jié)點1和節(jié)點3處根據(jù)基爾霍夫定律列出節(jié)點方程，計算求解相關傳熱速率。

圖3　圓柱體外殼電路類比Fig.3　Circuit analog of cylindrical enclosure

2.2建立固液分界面的幾何狀態(tài)方程

幾何狀態(tài)方程涉及CZ法晶體生長系統(tǒng)的幾何形變，尤其是熔體-晶體分界面處的形變，即分界面半徑、分界面高度、接觸角、有效半徑以及熔體高度。該狀態(tài)方程以熔體-晶體分界面處的傳熱傳質現(xiàn)象為基礎。

熔體高度：

(14)

分界面高度：

(15)

分界面半徑：

(16)

有效半徑：

(17)

τ=fmRi/Vp

(18)

分界面高度和接觸角：

(19)

(20)

(21)

Dφ=βRi(8Ri2cosφ+βcosφsinφ-sinφS1)

(22)

式中，ρs和ρl分別是硅固體和硅液體的密度，qx是晶體的傳熱速率，qi是從彎月板區(qū)域到固液分界面的傳熱速率，Hf是材料的融合比熱，φ0是非零接觸角，β是材料的拉普拉斯常數(shù)，fm是介于0.25～0.50間的比例參數(shù)。

3　模型驗證及仿真實驗

仿真中的狀態(tài)參數(shù)如表1、表2所示[9]，在Matlab中對基本模型和輻射增強模型中的輻射角系數(shù)進行計算仿真，結果如圖4所示。

表1　仿真中的常量值

表2　硅晶體生長參數(shù)

圖4　基本模型和輻射增強模型中輻射角系數(shù)比較Fig.4　Radiation angle coefficients in basic model and radiation enhancement model

圖4(a)說明了基本模型中坩堝到熔體表面的輻射角系數(shù)Fcs1和輻射增強模型中坩堝到熔體表面的輻射角系數(shù)Fcs2隨熔體高度Hm的變化趨勢；圖4(b)則說明了基本模型中熔體表面到周圍環(huán)境的輻射角系數(shù)Fso1和輻射增強模型中熔體表面到周圍環(huán)境輻射角系數(shù)Fso2隨熔體高度Hm的變化趨勢。從圖中可以看出，這兩種類型的輻射角系數(shù)在數(shù)值上并不相同，但在趨勢上是相同的，都是隨著熔體高度的增加而遞增。其中，Gevelber等人已驗證了基本模型在CZ法晶體生長系統(tǒng)中應用的正確性[9]，因此這兩種模型的輻射角系數(shù)的計算方法都能夠用于CZ法晶體生長系統(tǒng)的研究。

k1Tc4+k2Tc+k3=0

(23)

從而得到Tc的解軌跡表達式，k1、k2、k3為該方程的系數(shù)。同理，依次推導Th和Pin的解軌跡表達式。利用maple軟件進行解軌跡的計算，結果如圖5所示。

圖5　初始條件下系統(tǒng)開環(huán)解軌跡Fig.5　Trajectory of open loop solution under initial conditions

在表1和表2所示的初始條件下(即：Ri=2cm，Tm=1 740K，Vp=2cm/h)，坩堝溫度、加熱器溫度、輸入功率隨熔體高度Hm的變化軌跡如圖5(a)、(b)、(c)所示，在圖中分別用Tc0、Th0、Pin0代表初始條件下的變化軌跡。Hm隨t的變化軌跡如圖5(d)所示。由于在實際工業(yè)生產中，坩堝內熔體的高度是不斷下降的，因此觀察圖5所示的溫度變化軌跡時，應從坐標軸的橫軸由右向左看，從圖中可以看出，單晶爐內的主要組分的溫度在晶體生長過程中的變化都是隨著熔體高度的下降而逐漸升高的。此外，Hm隨時間的增加而線性遞減，其變化非常小，一小時約下降0.5 mm，因此可以認為在一定的時間范圍內，熔體高度是不變的，這也就意味著，在一定程度上單晶爐內的溫度是恒定的。

圖7為坩堝溫度隨熔體高度的變化關系。從圖中可以看出其趨勢是相同的，僅在數(shù)值上有所不同，其中輻射增強模型的溫度值更高一些。這一現(xiàn)象完全符合實際物理特性，由于輻射增強模型每次輻射換熱時同時考慮到三個表面之間的相互作用，這與基本模型相比，捕獲了更多的輻射換熱損失。

圖6　熔體溫度Tm和qi/qx的關系Fig.6　The relation between Tm and qi/qx

圖7　坩堝溫度隨熔體高度的變化關系Fig.7　The relation between crucible temperature and melting height

4　系統(tǒng)分析

本研究所建立的系統(tǒng)方程是高度耦合的、非線性的以及時變的，這意味著對系統(tǒng)的分析會是在數(shù)學上的加強分析，而且控制器的設計會被復雜性所限制。因此，控制系統(tǒng)開發(fā)的第一步是利用線性化理論對該系統(tǒng)進行線性化。

4.1系統(tǒng)線性化

根據(jù)所建立狀態(tài)方程的復雜性以及易于實現(xiàn)的程度，本研究采用泰勒級數(shù)逼近法對狀態(tài)方程線性化。系統(tǒng)模型是具有如下形式的一組方程式：

(24)

y(t)=g(x(t),u(t))

(25)

式中，x(t)是一個8×1的狀態(tài)向量(Th、Tc、Tm、Ri、Hi、Re、Hm和φ)，u(t)是一個2×1的輸入向量(Vp和Pin)，y(t)是一個關于衡量參數(shù)(Ri和Tm)的2×1輸出向量。

根據(jù)泰勒級數(shù)展開式可以將系統(tǒng)線性化為標準狀態(tài)空間模型，系數(shù)矩陣A(t)和B(t)分別是f對x(t)和u(t)求偏導的函數(shù)行列式，即：

A(t)=

CZ法晶體生長系統(tǒng)模型的輸出方程為：

(26)

4.2系統(tǒng)特征值

系統(tǒng)的特征值定義為特征方程∏(λ)的根。表3給出了在不同的熔體高度Hm下兩種模型的線性化系統(tǒng)的特征值。由于Hm每下降1cm對應大約20小時的晶體生長，因此，特征值隨時間的增加幾乎無變化。

表3　系統(tǒng)開環(huán)特征值

4.3系統(tǒng)穩(wěn)定性

系統(tǒng)的特征值結果為該系統(tǒng)是不穩(wěn)定系統(tǒng)提供了數(shù)據(jù)支持。從表3中可以看出，線性化的模型具有右半平面的特征值，這表明系統(tǒng)是不穩(wěn)定的。

4.4系統(tǒng)可控性

為了確定僅在加熱器功率作用下的系統(tǒng)是否可控，通過計算可控矩陣的秩進行判斷。

在本文所建立的系統(tǒng)中，系統(tǒng)階數(shù)是8，采用符號計算求解矩陣A的秩，結果為6，可以得出：僅在加熱器功率作為輸入的情況下系統(tǒng)是不可控的。由于本文所建立的線性化模型是時變的、高度耦合的復雜模型，因此不能采用常規(guī)PID算法設計控制器，相比之下采用易調整、能夠系統(tǒng)解決約束和時延等實際問題的模型預測控制算法作為控制策略。

5　CZ法晶體生長系統(tǒng)的模型預測控制

模型預測控制(MPC)是一種基于模型的控制算法，包括經(jīng)典模型預測控制和綜合模型預測控制，利用內部模型的狀態(tài)或輸出預測，同時應用了有限時域的滾動計算思想和反饋及預測校正，最后采用了對某個系統(tǒng)性能指標的最優(yōu)化計算以確定在一個控制時域內的最優(yōu)控制序列。由于其魯棒性強、抗干擾性強，而且能夠在優(yōu)化控制理論的框架內很好地處理系統(tǒng)的控制約束等特點，非常適用于實際工業(yè)過程的控制，而CZ法晶體生長過程本身便是一個典型的工業(yè)過程。

CZ法晶體生長過程的閉環(huán)控制系統(tǒng)初步設計如圖8所示，其中r是參考軌跡，u是輸入，y是輸出，w、z是干擾。對于直拉法晶體生長系統(tǒng)，u是系統(tǒng)輸入，即：加熱器功率Pin和提拉速率Vp；y是系統(tǒng)輸出，即熔體溫度Tm和分界面半徑Ri；w是加熱器功率Pin上的干擾；z是分界面半徑Ri上的干擾。

圖8　閉環(huán)控制系統(tǒng)Fig.8　Closed loop control system

根據(jù)標稱無延遲線性系統(tǒng)(即無干擾且不存在內部模型和被控系統(tǒng)的不匹配)設計無約束MPC控制器來研究系統(tǒng)性能。系統(tǒng)的激勵給定為階躍信號。

通過對預測時域、控制時域選擇不同的取值，進行仿真實驗，確定合適的參數(shù)。

當設定熔體溫度Tm為5K的階躍激勵、晶體半徑Ri設定為常數(shù)0時，仿真結果如圖10所示，這一結果也在可接受范圍內。

圖9　Ri以0.5mm的步長變化時的系統(tǒng)響應Fig.9　Response to a 0.5mmstep change in Ri

圖10　Tm以5K的步長變化時的系統(tǒng)響應Fig.10　Response to a 5Kstep change in Tm

設計有延遲系統(tǒng)的控制器需要考慮在標稱系統(tǒng)中加入時延，即考慮輸入功率Pin到系統(tǒng)輸出的純時延，初步設定為60s。將預測時域P和控制時域M默認值下的無延遲結果與有延遲結果以及增加P和M取值的有延遲結果進行比較。當P和M取默認值且有延遲時，閉環(huán)系統(tǒng)無法執(zhí)行并且是不穩(wěn)定的；當P和M取值增加時，系統(tǒng)性能得以恢復。最后得出的仿真結果如圖11所示，在可接受范圍內，認為控制器性能良好。

圖11　Ri以0.5mm的步長變化時的系統(tǒng)響應，時延60 sFig.11　Response to a 0.5mmstep change in Ri, delay=60 s

6　結　語

本文基于傳熱傳質基本原理建立了CZ法晶體生長等徑階段的低階模型，并通過驗證開環(huán)特征根的方式證明了模型的有效性，采用模型預測控制算法針對系統(tǒng)有無延遲的情況分別設計了相對應的無約束無干擾控制器，在可接受范圍內，認為控制器的控制效果良好。本研究下一步將在控制器的設計中考慮約束和干擾的問題，進行更結合實際工業(yè)問題的研究。

[1] 徐赟瑜. 直拉法單晶硅熔體內熱量動量及質量輸運的數(shù)值模擬[D]. 重慶:重慶大學,2009:5.

XU Yunyu. Numerical simulation of heat, momentum and mass transport in the melt of Czochralski silicon single crystal[D]. Chongqing：Chongqing University,2009:5.

[2] DERBY J J, BROWN R A. Thermal-capillary analysis of Czochralski and liquid encapsulated Czochralski crystal growth II. processing strategies[J]. Journal of Crystal Growth，1986，75(2):227-240.

[3]DERBY J J, BROWN R A. Thermal-capillary analysis of Czochralski and liquid encapsulated Czochralski crystal growth: I.Simulation[J]. Journal of Crystal Growth，1986，74(3):605-624.

[4]DONALD T J H. Crystal pulling from the melt[M]. New York：Springer Verlag，1993.

[5]BRICE J C , BRUTON T M, HILL O F,et al. The Czochralski growth of Bi12SiO20crystals[J]. Journal of Crystal Growth, 1974,24-25:429-431.

[6]JONES D W. An experimental model of the flow in Czochralski growth[J]. Journal of Crystal Growth,1983, 61(2): 235-244.

[7]SEIDL A, MARTEN R, MTILLER G. In situ investigation of oxygen distribution and transport in Czochralski silicon melts by electrochemical solid ionic sensors[J]. Journal of Crystal Growth. 1996,166(1-4):680-684.

[8]GEVELBER M, George S. Dynamics and control of the Czochralski process I: Modeling and dynamic characterization[J].Journal of Crystal Growth, 1987, 84(4): 647-668.

[9]GEVELBER M, George S. Dynamics and control of the Czochralski process II:Objectives and control structure design[J].Journal of Crystal Growth,1988, 91(1-2): 199-217.

[10] GEVELBER M, George S. Dynamics and control of the Czochralski process III:Interface dynamics and control requirements[J]. Journal of Crystal Growth,1994, 139(3-4): 271-285.

[11]GEVELBER M, George S. Dynamics and control of the Czochralski process IV:Control structure design for interface shape control and performance evaluation[J]. Journal of Crystal Growth, 1994, 139(3-4): 286-301.

[12]SCHEEL H J. Historical aspects of crystal growth technology[J]. Journal of Crystal Growth, 2000, 211(1):1-12.

[13]蔣榮華,肖順珍. 半導體硅材料的進展與發(fā)展趨勢[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.

[14]張尚中.提拉法晶體生長全局數(shù)值模擬[D]. 重慶:重慶大學資源與環(huán)境科學學院，2009:4.

ZHANG Shangzhong. Global simulation of Czochralski crystal growth[D]. Chongqing: College of Resource and Environmental Science of Chongqing University，2009:4.

[15]MARTINEZ D M. Modeling and control of the Czochralski crystal growth process[D]. Texas: Texas A&M University,2002.

(責任編輯周蓓)

Modeling based on resistive network and simulation of MPC method on Czochralski growth process

YI Yingmin， ZHANG Tong

(School of Automation and Information Engineering, Xi’an University of Technology, Xi’an 710048,China)

For strong nonlinear with large delay on CZ growth process of silicon, the model based on the heat and mass transfer phenomenon is developed. The steady state solutions in the initial conditions are calculated, linearization and analysis of system model are conducted, with control algorithm and designing determined and the controller based on the controllability and stability of the system designed. It is concluded that the simplifications in radiation view factors given by Gevelber can reflect the characteristics of crystal growth in industry. Moreover, the analysis of system performance indicates that the system is uncontrollable and unstable on the premise of heater power as the input, thus using the model predictive control algorithm as the control strategy. Then the unconstrained MPC is designed for the nominal, delay-free, linearized system and the unconstrained MPC for the nominal, delay, linearized system, with the result that verifies the validity of these methods.

single-crystal silicon; the Czochralski process; the resistive network analogy; the lower order model; model predictive control

1006-4710(2016)02-0142-07

10.19322/j.cnki.issn.1006-4710.2016.02.003

2015-05-21

國家重點基礎研究發(fā)展計劃973項目資助(2014CB360508)；國家自然科學基金資助項目(61533014);陜西省教育廳科研資助項目(16JS069)

弋英民,男,教授,博士,研究方向為晶體生長模擬仿真研究。E-mail: yiym@xaut.edu.cn

O782+.5

A