直拉硅单晶非均匀相变温度场最优控制

张 晶,刘 丁,杜燕军

(西安理工大学晶体生长设备及系统集成国家地方联合工程中心,陕西西安 710048)

1 引言

直拉法生长直径均匀的硅单晶是目前半导体芯片材料制造的关键技术[1-2],也是衡量晶体品质的重要指标之一.在硅单晶生长过程中,通过提拉运动在相变界面(即熔体-晶体交界面)处形成一定的温度差,依据晶体生长原理[3],相变界面温度场的均匀性决定了该温度差的稳定性,是硅单晶均匀生长的根本保证.然而,相变温度场受提拉速度、加热器温度、晶体旋转和坩埚旋转等工艺参量的影响且随着热能、动能和质量输运的状态而变化,呈现非均匀性特征.

当前,针对相变温度场非均匀性问题,在晶体生长温度建模的研究中,在表征固液界面径向温度梯度和轴向温度梯度共同作用规律时,往往忽略了提拉速度动力学因素[4-6].然而,由拉速引起的结晶潜热所引发的热补充对固液界面轴向温度梯度有着不可忽视的影响.Brown[7]和Lan[8]分别给出了详细的晶体生长过程建模和动力学分析方法,合理的调整晶转和埚转工艺参数,能够提高固液界面温度的均匀性.然而,随着硅单晶生长直径增大,热场尺寸增加,液面位置变化等因素仍对相变温度场的均匀性产生影响,从而引起晶体直径波动.与此同时,熔体内部热传输机理复杂,晶体生长直径控制系统所具有的强非线性、大时滞、不确定性特征,给控制器设计带来困难[9-10].因此,从晶体内部热传输和温度分布入手研究其直径控制问题是一有效途径.Derby等人[11]通过研究质量和能量平衡关系建立了晶体的热传输模型,该模型为时变空间域上的偏微分方程.Wang等人[12-13]研究了基于稳定温度场的最优控制器,提出了时变空和间域的分布式系统的稳定与控制模型.然而,该模型描述的晶体温度随时空变换且具有无穷个自由度,本质上是无限维的,难以获得解析解.同时,受实际生产条件限制,生长温度难以实时测量,大大增加了直径控制难度.

针对这一问题,对无限维分布参数系统进行降维是分析、优化和控制生长直径的有效方法[14].传统的降维方法如差分法和有限元法只能得到阶数非常高的近似模型,并不适用于分布参数系统的快速仿真和控制器设计[15-17].因此,研究合理的降维方法来获得较为准确的温度场模型是十分重要的.针对上述问题,本文建立一种基于偏微分方程(partial differential equation,PDE)模型的温度场优化模型,该模型采用二阶常微分方程(ordinary differential equation,ODE)描述了机械拉动臂引起的边界演化,考虑生长速率波动的时变空间域上的抛物型PDE来描述晶体温度的变化.基于PDE无限维空间上晶体热传输的时变模型,采用谱方法和Galerkin方法获得有限维近似模型,以此来简化温度调节和域运动控制问题的数值实现.采用线性二次型对改进后的温度模型进行控制,对比分析了相变界面温度分布和晶体生长直径,从而验证方法的有效性.

2 非均匀相变温度场问题描述

2.1 问题分析

依据点缺陷在固液界面处的相变温度场形成理论,相变温度场的某时空临界尺度由比值决定,V为晶体的提拉速度,G为固液界面靠近晶体侧的轴向温度梯度.当比值维持在0.134(1±10%)mm2/◦C·min范围时,晶体可在稳定的相变温度场中生长,形成无缺陷的完美硅单晶[18].当相变温度场不稳定波动超过理想范围时,会引起晶体生长速率波动并反映在晶体等径生长时直径的变化上.实际中往往通过提拉速度来调控固液界面相变温度场的均衡性,这是由于直径能够快速响应拉速调节,而通过加热器温度调节晶体直径的过程中,由加热器传至固液界面的热传输通道存在较大滞后,因而调节响应慢.因此,由生长速率波动引起的提拉速度的波动是导致晶体直径波动的根本原因.

本质上,相变温度场不稳定波动一方面来自提拉速度对温度分布影响,即晶体不断生长使得晶体边界发生变化,影响固液界面处的温度分布.另一方面熔体液位不断下降导致热辐射增强,影响了固液界面处温度分布的均匀性.因此,建立晶体生长过程提拉动力学模型与热传输模型时,需要进一步考虑晶体生长和熔体液位下降对相变温度场的影响,从而改善相变温度场的非均匀性.

2.2 提拉动力学模型

根据晶体生长原理,晶体生长过程满足质量守恒定律和牛顿第二定律,如图1所示.

图1 直拉法晶体生长过程示意图Fig.1 Schematic diagram of the Czochralski crystal growth progress

文献[19]中假设固液界面平直,提出基于晶体生长提拉过程的动力学模型,该模型能够很好的描述提拉作用对晶体直径和晶体生长速率的影响.但是,在晶体生长过程中,固液界面处热通量并不是一成不变的,会随着温度场的变化而发生改变,从而影响晶体生长直径.由第2.1节分析可知,液面位置下降所引起的相变温度场非均匀性是改变固液界面形状的根本原因,如图1所示,固液界面形状凸向晶体.针对这一问题,通过建立考虑坩埚上升速率的提拉动力学模型,以保证液面位置的稳定来改善相变温度场的非均匀性,从而避免生长速率的波动.

在晶体生长过程中,根据质量守恒定律可得Mt=Ml+Mm+Ms,Mt为总投料量,Ml为熔体质量,Mm为弯月面的质量,Ms为晶体质量.假设在晶体生长过程中Mm恒定，则熔体减少的质量等于晶体增加的质量,即

其中k为埚跟比.

直拉法硅单晶生长过程中,相变界面处能量的传递和变化均直接影响界面温度,从而影响生长界面温度差.图1给出了晶体生长过程中固液界面的热量传递过程.根据固液界面的能量守恒方程,可得晶体生长速度与热量的关系为

其中:Φs表示穿过固液界面进入晶体的热量,Φl表示从熔体进入固液界面的热量,Φh表示晶体结晶时释放的潜热,∆H为晶体结晶时释放的潜热.上式表明,在单位时间内相变界面传递到晶体中的热量等于通过熔体传递至相变界面的热量和相变过程中释放的潜热之和.因此,希望相变界面处有稳定热量传输就意味着Φs−Φl在晶体生长过程中保持稳定.

根据牛顿第二定律可得

联立式(1)-(5),可得提拉动力学模型

2.3 温度场优化模型

根据提拉动力学模型得到了晶体边界演化和晶体的生长速率,通过对晶体生长动力学特性与热传输方式的分析,式(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)可知,晶体的温度是相对时间和空间位置的函数,具有无穷维自由度,且往往不能通过有限个状态进行求解.因此,需要对该模型进行降维.

3 晶体生长温度场降维

针对无限维分布参数系统的求解问题,采用有限维的常微分方程描述的系统来近似无限维分布参数系统.首先需要确定热传输作用的有效边界,即晶体直径和长度的变化.在晶体生长过程中,总是期望晶体直径恒定,因此,控制晶体直径尤其重要.本文在提拉动力学模型和直径控制的基础上确立热传输边界,进一步对热传输模型进行降维.

3.1 热传输边界的确立

根据式(6)-(7)所示的晶体提拉动力学模型,对输入输出进行线性化,可定义Fext(t)为

3.2 PDE模型降维方法

基于上述方法确立的边界条件对PDE模型进行降维.谱方法在整个空间域上选择全局和正交的空间基函数,相比较于有限单元法的空间离散能够得到维数更低的近似模型.谱方法适用于一类能够进行快慢变量分离的系统,快慢变量是指空间基函数在空间频率域中的排列顺序往往是由慢到快,由于快变量对系统的贡献较小,因此,选择慢变空间基函数即可满足降维的要求.快慢变量的区分通过Galerkin截断准则来判定.

谱方法对模型降维的基本步骤为

1) 选取空间基函数,将系统变量在空间基函数上展开.

空间微分算子的特征函数常用作空间基函数[14],对于算子A,其与特征值和特征函数之间关系为

式中:λ为特征值,Φ为特征函数,可通过式(14)对算子A分离变量进行求解[19].

将系统变量x(r,z,t)在空间基函数上展开,即

对于式(9)-(10)所描述的系统,输入控制u(t)加在边界上,并假设输入控制作用函数为

该函数代表输入的作用范围,其中ε1>0,通过采用Drica delta函数将边界控制系统转化为与其等价的分布控制系统,则输入算子B(t)可以定义为

结合式(9)和式(15)可以将热传输模型写为

2) 采用Galerkin方法得到描述展开系数的无限维常微分方程组.

Galerkin方法中选取与空间基函数相同的权重函数使系统残差最小,获得展开系数的无限维常微分方程组,记为

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

3) 通过有限维截断与空间基函数综合得到系统的低维近似模型.

为了对无穷维常微分方程降维,需根据Galerkin截断理论确定慢变空间基函数的个数.Galerkin截断准则[20]如下.

对于分布参数系统,Reλj表示的是系统空间算子特征值λj的实部,如果0 ≥Reλ1≥Reλ2≥···≥Reλj≥···,且存在k使之满足

3.3 温度控制

对于系统的最优控制,线性二次型性能指标能够更好的兼顾系统的状态以及输入能量消耗.因此,将线性二次型最优控制应用于抛物PDE系统低维模型,结合晶体域运动来简化晶体温度调节,得到系统的最优输入与温度分布.设计一个最优输入u(t),使如下的性能指标函数达到最小值:

其中:Q和R为非负定对称和正定对称的系数矩阵,a为系统的状态.为使性能指标J最小,最优输入取为

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

4 结果与分析

本次仿真实验数据来自本中心12英寸硅单晶炉的实际拉晶试验,记录了整个晶体生长过程中直径,生长速率和加热器功率等数据,采样时间为2 s.选取放肩和等径阶段10000个直径数据作为晶体的参考直径,其中前1500个数据为放肩阶段,其余数据为等径阶段.仿真过程中各初始状态取值参考实际晶体生长数据,如表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值逐渐减小.在等径阶段,固液界面热通量相对稳定,晶体保持稳定生长,符合晶体生长工艺要求,验证了模型的正确性.由图2可知,C2值波动较小,说明相变温度场波动较小,在等径阶段维持液面位置不变有利于改善固液界面相变温度场的均匀性.因此,在晶体生长过程中考虑坩埚上升速率能够有效地减少固液界面下降引起的相变温度场变化.

针对改进前和改进后两种模型对晶体直径控制进行仿真,仿真过程设置控制增益K=0.1,可得晶体生长速率变化曲线如图3所示.

图3 晶体生长速率Fig.3 Crystal growth rate

图3给出了晶体在放肩和等径阶段生长速率的变化曲线.由C值分析可得,在放肩阶段,提拉速度较高,生长速率近似于提拉速率,因此,在放肩阶段晶体生长速率较高.由放肩至等径的过程需要进行转肩工艺,即需提高拉速抑制晶体直径放大,如图3所示,此处拉速明显升高.将优化后的生长速率进行对比,由图3可知,改进前模型中C值波动较大,为了能够获得稳定的直径,须通过调节拉速来减小直径波动,即导致晶体生长速率波动较大;而改进后模型中C值较为稳定,能够获得稳定的生长速率.

图4给出了改进前和改进后模型直径控制结果,并给出了与期望直径之间的误差.

由图4可知,改进后模型能够更好稳定晶体的直径.通过计算均方根误差对结果进行精确评估,改进前模型求得直径的均方根误差为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的特征值满足截断条件时,可以将分布参数系统进行快慢分离.一般,当=0.1,即λk十倍于λ1时,就可截断前k阶低维慢变量部分来近似该无限维系统.由表2可得.径向特征值满足0 ≥Reλ1≥Reλ2≥···,

因此,截取到前4个特征值即可满足降维要求.同理,根据图6所示,轴向特征函数对应的特征值只需截取前4个即可满足降维要求，即式(19)中M=4,N=4,因此,将分布参数系统模型降维到16维.

根据特征值的截断数,在图7-8分别给出了径向和轴向前4个特征函数的变化曲线,其作为空间基函数来实现模型降维.从图7-8可以看出,特征函数之间频率变化明显,具有良好的快慢分离特性,适于模型降维.

图7 前4个径向特征函数Fig.7 The first 4 radial eigenfunctions

图8 前4个轴向特征函数Fig.8 The first 4 axial eigenfunctions

为了稳定固液界面附近的温度分布,对降维后的模型进行线性二次型最优控制,图9仅给出了前4维时间系数amn(t)的变化曲线:

最优控制得到系统的输入变化曲线如图10所示.

图11给出了降维后固液界面处晶体径向温度随时间的变化.从中选取不同时刻的温度分布加以分析,如图12所示.图12中初始时刻晶体径向温度存在波动,通过对系统进行最优控制后,晶体固液界面附近温度逐渐稳定.

图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中可以看出,固液界面处温度分布随时间增长而更加平稳.图12中对模型改进前、后得到的温度差分布比较可以得出,改进前与改进后模型中固液界面径向温度梯度的差值逐渐增大,是由于随着晶体的不断生长,液面位置不断下降,使得相变界面的非均匀性增强,且随着时间推移表现更为明显.改进后模型能够大幅度提高径向温度分布均匀性,得到的固液界面径向温度梯度较小,能够获得更加平坦的固液界面形状.在线性二次型最优控制作用下,改进前后模型的固液界面径向温度梯度在逐渐减小,最终获得稳定的相变界面温度分布.在这一过程中,改进后模型温度梯度低于改进前模型,进一步验证了改进模型与降维方法的有效性.

图11 固液界面处温度分布Fig.11 Temperature distribution at solid-liquid interface

图12 不同时刻温度分布曲线Fig.12 Temperature distribution curve at different time

5 结论

直拉法晶体生长相变界面温度分布非均匀性导致晶体生长直径不均匀.针对这一问题本文提出了基于PDE模型的温度场最优控制策略.建立了考虑坩埚上升的改进提拉动力学模型,确定晶体域边界的演化,从而使液位下降导致的相变温度场非均匀性得到改善.采用谱方法选择特征函数作为空间基函数,Galerkin将时间系数展开为无限维常微分方程组,并根据截断准则得到低维近似模型,解决了PDE模型无法直接求解的问题.通过反馈线性化和线性二次型来求解域运动以及域运动单向耦合的温度分布问题.仿真实验结果表明,改进模型使得相变界面径向温度分布更加均匀,有利于直拉硅单晶稳定生长.