CN108754599A - 一种基于有限元数值模拟的硅单晶生长温度控制方法 - Google Patents

Abstract

本发明公开了一种基于有限元数值模拟的硅单晶生长温度控制方法，首先进行参考晶体形状下的有限元建模，设置整个过程中晶体的参考形状和边界条件；以晶体冷却速度最大为目标，以加热器功率的可调范围和晶体内温度梯度的最大值为约束条件，利用优化方法离线优化计算出满足约束的晶体温度最优轨迹和加热器功率最优轨迹，并将其作为控制系统的参考输入；然后对晶体半径进行控制；最后实现对晶体内部温度轨迹的跟踪控制，本发明解决了现有技术中存在的硅单晶生长控制方法无法满足晶体温度控制、造成晶体位错缺陷的问题。

Description

一种基于有限元数值模拟的硅单晶生长温度控制方法

技术领域

本发明属于直拉硅单晶生长过程温度自动控制技术领域，具体涉及一种基于有限元数值模拟的硅单晶生长温度控制方法。

背景技术

硅单晶是最重要的半导体材料。90％以上的集成电路器件制作在硅单晶上。近年来极大规模集成电路工艺快速发展，特征线宽不断缩小，对晶体微观品质的要求也越来越高。位错是衡量晶体品质的重要指标之一，而温度梯度和热应力是形成晶体位错缺陷的主要原因。热应力产生的根源在于，晶体内部由于热胀冷缩特性而存在的温度梯度。温度梯度越大，热应力越大，当热应力超过临界值时，晶体出现位错缺陷，甚至导致晶体断裂，从而导致晶体特性不符合要求。因此，对生长中的晶体温度梯度进行及时有效地控制可从根源上避免晶体位错缺陷的产生，这也是该领域的技术难点和研究热点。目前的单晶生长控制系统大多仍停留在对晶体直径的控制，控制方法多为基于测量数据的比例-积分-微分(Proportion Integration Differentiation,简称PID)控制方法和基于数据驱动辨识系统模型的模型预测控制(Model Predictive Control，简称MPC)。PID控制属于无模型控制，只需要根据检测值和参考值的偏差即可对被控量进行控制。基于数据驱动的模型预测控制，是根据测量到的数据建立一个近似模型从而对被控量进行控制。这两种控制方法的前提都是被控变量是可检测变量。由于晶体内部的温度梯度是不可测变量，因此，现有的控制方法无法对晶体位错进行有效控制。针对晶体位错的控制，目前主要通过实验来观察工艺参数与晶体位错之间的关系，不仅实验成本高、周期长，工艺参数的调节也具有一定程度的盲目性。因此，提出一种可同时控制晶体直径和晶体内部温度、温度梯度的控制方法以满足对晶体尺寸和晶体品质的高需求是目前亟待解决的重要问题。

发明内容

本发明的目的是提供一种基于有限元数值模拟的硅单晶生长温度控制方法，解决了现有技术中存在的硅单晶生长控制方法无法满足晶体温度控制、造成晶体位错缺陷的问题。

本发明所采用的技术方案是，基于有限元数值模拟的硅单晶生长温度控制方法，具体按照以下步骤实施：

步骤1、参考晶体形状下的有限元建模，设置晶体生长过程中晶体的参考形状和边界条件；

步骤2、以晶体冷却速度最大为目标，以加热器功率的可调范围和晶体内温度梯度的最大值为约束条件，利用优化方法计算出满足约束的晶体温度最优轨迹和加热器功率最优轨迹，并将其作为控制系统的参考输入；

步骤3、对晶体半径进行控制：首先建立晶体生长过程中提拉动力学模型并进行求解，得到晶体的实时半径和晶体高度，然后对比实际半径与晶体参考半径之间的偏差，利用PID控制算法对晶体实际输出半径与高度进行控制；

步骤4、实现对晶体内部温度变化轨迹的跟踪控制：首先在步骤1建立的有限元模型中，指定监测点P的位置，然后以步骤2得到的加热器功率轨迹作为参考加热器功率，以步骤2得到的最优温度轨迹作为P点的参考温度，以步骤3得到的晶体半径和高度作为动态有限元模型的实时半径和实时高度，然后利用MPC控制算法对晶体温度进行实时控制；

执行步骤1和步骤2得到整个晶体生长过程中晶体监测点的最优温度变化轨迹和最优加热器功率轨迹，将其作为真实晶体生长中的设定轨迹；重复步骤3得到晶体生长过程中任意时刻晶体的真实半径，在PID控制器的作用下，晶体半径与设定半径趋于一致；重复步骤4得到晶体生长过程中任意时刻晶体内部的温度分布和最大温度梯度，在MPC控制器的作用下，模拟得到的监测点P处的温度将与设定轨迹趋于一致，同时加热器功率也得到确定，至此，基于有限元数值模拟的晶体温度控制过程结束。

本发明的特点还在于，

步骤1晶体生长的有限元建模，具体按照以下步骤实施：

步骤1.1、晶体热传导方程如下式所示：

式中▽表示圆柱坐标系下的空间梯度，r表示晶体半径，z表示晶体晶体高度，T(r,z,t)表示t时刻晶体中(r,z)处的温度，Pe＝ρC_pv₀R_cruc./k_s为Peclet常数，ρ、v₀、R_cruc、k_s、k_r、C_p分别为密度、提拉速度、坩埚半径、热传导系数、热传导率和比热容，V(r,z,t)为晶体区域的速度，T(r,z,t)为晶体区域的温度；

步骤1.2、用变量r和z分别描述晶体的半径和高度，用表示晶体边界上指向晶体内部的法向量，设置晶体的边界条件如下：

当z＝0时，T|_z＝0＝T_f＝1411.85℃；

当0＜z＜z₁时，

当z₁≤z≤z₂时，

当z₂＜z＜z_top时，

当z＝z_top时，

步骤1.3、根据式(1)用有限元法求解晶体中的温度分布：首先将晶体区域划分成Np个网格单元，然后对公式(1)进行展开，结合设置的对应边界条件，将热传导模型转换为有限元模型，如式(2)所示

式中A_k、B_k是有限元模型中的时变矩阵，是与晶体表面第一类边界条件相关的矩阵，X_k代表晶体内部k时刻所有网格单元的温度矩阵，Y_k代表k时刻所有网格单元的径向温度梯度矩阵，u_k代表k时刻的加热器功率值，C_k是将温度转化为温度梯度的输出矩阵。

步骤2计算最优的温度变化轨迹和加热器功率轨迹，具体如下：

目标函数及约束条件如下：

subject to:

X_k+1＝A_kX_k+B_ku_k+B_Ck (3)

u_min≤u_k≤u_max

|Y_k|≤G_max

式中u_min、u_max分别为加热器功率可调范围的上下限，G_max为晶体内部的最大温度梯度约束，Q和R均为正定矩阵，利用优化算法对上述优化问题进行求解，得到整个晶体生长过程中加热器的最优功率轨迹和晶体的最优温度轨迹，优化算法的最大迭代次数为Np。

步骤3对晶体实施半径控制，具体如下：

晶体提拉动力学模型如下式所示：

式中，r(t)、z(t)分别表示t时刻晶体的实时半径和实时高度，V_c(t)、ρ分别表示晶体的体积增长速率、坩埚半径、晶体体积和晶体密度，F_ext(t)表示晶体的牵引力，也是晶体直径控制的控制量，根据式(4)得出晶体的实时半径r(t)和实时高度z(t)，将实时半径与步骤1中设置的晶体参考半径r_d(t)比较，根据两者之间的偏差e(t)＝r(t)-r_d(t)，采用PID算法对控制量F_ext(t)进行调整，调整规则如下：

其中，V_c0为晶体的初始体积，K_P和K_I分别为控制器的比例调节系数和积分调节系数，将式(5)代入式(4)计算出PID控制作用下晶体的实时半径和实时高度，为方便计算，假设晶体体积增长速率C_growth为常量，将晶体体积增长速率的变化作为一种扰动d施加在C_growth上。

步骤4对晶体内部温度实施轨迹跟踪控制，具体如下：

步骤4.1、模型预测：建立有限元模型如式(6)所示，根据该模型预测出晶体内部的温度及温度梯度，包括P点的温度：

步骤4.2、滚动优化：使将晶体实际温度X_k与参考温度轨迹之间的误差平方和以及实际加热器功率u_k与参考加热器功率轨迹之间的误差平方求和最小作为优化目标，以监测点温度梯度不超过临界值G_max为约束，同时对加热器功率加以约束，目标函数和约束条件设置如式(7)所示，用优化算法对该优化问题进行求解，优化得到一组最佳的控制量u，即加热器功率：

subject to:

y_k＝C'_kX_k

u_min≤u_k≤u_max

|Y_k|≤G_max

式中，A'_k和B'_k分别表示根据步骤3计算出的晶体实时半径和实时高度推导出的有限元时变矩阵，Nc是优化算法的最大迭代次数；

步骤4.3、反馈校正：将每次滚动优化得到的第一个控制量施加在被控对象上，得到被控对象监测点P处的实时温度，将监测点P处的温度X^J反馈到MPC控制器，与MPC的预测值X^Y进行对比，根据两者的偏差对MPC的预测值进行修正，如式(8)所示，其中ε为校正系数：

步骤4.4、再重新进行步骤4.1～步骤4.3，预测模型在每一次迭代过程中都将被更新以更加准确的预测晶体内的温度和温度梯度分布。

本发明的有益效果是，基于有限元数值模拟的硅单晶生长温度控制方法，根据晶体热传导方程，建立参考晶体形状下的晶体热传递过程的动态有限元数值模型，以晶体冷却速度最大为目标建立目标函数，根据晶体温度梯度要求和加热器功率的范围确定约束条件，利用优化方法离线优化计算出满足约束的晶体温度最优轨迹和加热器功率最优轨迹，并将其作为控制系统的参考输入；建立晶体生长提拉动力学模型，采用PID控制方法实时控制晶体半径和晶体高度的变化，使得实际晶体形状跟随参考晶体形状变化；建立考虑了热动力学模型影响的动态有限元模型，并对晶体形状加入扰动，然后将其作为被控对象，采用MPC控制方法对被控对象进行参考温度轨迹跟踪控制，使得被控对象的温度跟踪参考输入，同时保证被控对象的温度梯度满足约束条件，避免了晶体位错缺陷的产生。以上三个步骤可以实现对晶体半径和晶体温度的同步控制，不仅满足了晶体外形尺寸的要求，还保证了晶体内部的温度梯度符合生长工艺要求，避免了晶体位错缺陷的产生，突破了传统晶体生长控制方法无法满足晶体品质需求的局限。这种把数值模拟和控制方法结合起来的策略，不仅提高了计算效率，而且，工艺操作人员可不用通过实验逐次试参数的方法就能确定在不产生晶体位错缺陷的前提下加热器的最优功率轨迹和对应晶体内部最优温度轨迹。

附图说明

图1是本发明基于有限元数值模拟的硅单晶生长温度控制方法的控制框图；

图2是本发明基于有限元数值模拟的硅单晶生长温度控制方法中参考晶体几何模型、边界条件设置及温度监测点设置示意图；

图3是本发明基于有限元数值模拟的硅单晶生长温度控制方法中在获得晶体温度参考轨迹之后对晶体温度和直径同时控制的流程图；

图4是本发明基于有限元数值模拟的硅单晶生长温度控制方法中PID控制作用下的晶体的半径曲线图；

图5是本发明基于有限元数值模拟的硅单晶生长温度控制方法中PID控制作用下的晶体的高度曲线图；

图6(a)是本发明基于有限元数值模拟的硅单晶生长温度控制方法中计算得到的t＝13.6min放肩阶段晶体的温度分布结果图；

图6(b)是本发明基于有限元数值模拟的硅单晶生长温度控制方法中计算得到的t＝40min放肩阶段晶体的温度分布结果图；

图6(c)是本发明基于有限元数值模拟的硅单晶生长温度控制方法中计算得到的t＝192min等径阶段晶体的温度分布结果图；

图6(d)是本发明基于有限元数值模拟的硅单晶生长温度控制方法中计算得到的t＝256min等径阶段晶体的温度分布结果图；

图6(e)是本发明基于有限元数值模拟的硅单晶生长温度控制方法中计算得到的t＝352min等径阶段晶体的温度分布结果图；

图6(f)是本发明基于有限元数值模拟的硅单晶生长温度控制方法中计算得到的t＝400min收尾阶段晶体的温度分布结果图；

图7是本发明基于有限元数值模拟的硅单晶生长温度控制方法中从引晶到收尾整个过程中P点的参考温度轨迹和MPC控制作用下的温度轨迹曲线图；

图8是本发明基于有限元数值模拟的硅单晶生长温度控制方法中从引晶到收尾整个过程监测点温度梯度的变化过程图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明一种基于有限元数值模拟的硅单晶生长温度控制方法，晶体生长主要包括放肩、转肩、等径和收尾四个阶段，控制框图如图1所示，具体按照以下步骤实施：

步骤1、参考晶体形状下的有限元建模，设置整个过程中晶体的参考形状和边界条件，具体按照以下步骤实施：

步骤1.1、晶体热传导方程如下式所示：

式中表示圆柱坐标系下的空间梯度，r表示晶体半径，z表示晶体晶体高度，T(r,z,t)表示t时刻晶体中(r,z)处的温度，Pe＝ρC_pv₀R_cruc./k_s为Peclet常数，ρ、v₀、R_cruc、k_s、k_r、C_p分别为密度、提拉速度、坩埚半径、热传导系数、热传导率和比热容，V(r,z,t)为晶体区域的速度，T(r,z,t)为晶体区域的温度；

步骤1.2、用变量r和z分别描述晶体的半径和高度，用表示晶体边界上指向晶体内部的法向量，设置晶体的边界条件如下：

当z＝0时，T|_z＝0＝T_f＝1411.85℃；

当0＜z＜z₁时，

当z₁≤z≤z₂时，

当z₂＜z＜z_top时，

当z＝z_top时，

步骤1.3、根据式(1)用有限元法求解晶体中的温度分布：首先将晶体区域划分成Np个网格单元，然后对公式(1)进行展开，结合设置的对应边界条件，将热传导模型转换为有限元模型，如式(2)所示

式中A_k、B_k是有限元模型中的时变矩阵，是与晶体表面第一类边界条件相关的矩阵，X_k代表晶体内部k时刻所有网格单元的温度矩阵，Y_k代表k时刻所有网格单元的径向温度梯度矩阵，u_k代表k时刻的加热器功率值，C_k是将温度转化为温度梯度的输出矩阵；

步骤2、以晶体冷却速度最大为目标，以加热器功率的可调范围和晶体内温度梯度的最大值为约束条件，利用优化方法离线优化计算出满足约束的晶体温度最优轨迹和加热器功率最优轨迹，并将其作为控制系统的参考输入，具体如下：

目标函数及约束条件如下：

subject to:

X_k+1＝A_kX_k+B_ku_k+B_Ck (3)

u_min≤u_k≤u_max

|Y_k|≤G_max

式中u_min、u_max分别为加热器功率可调范围的上下限，G_max为晶体内部的最大温度梯度约束，Q和R均为正定矩阵，利用优化算法对上述优化问题进行求解，得到整个晶体生长过程中加热器的最优功率轨迹和晶体的最优温度轨迹，优化算法的最大迭代次数为Np；

步骤3、对晶体半径进行控制：首先建立晶体生长过程中的提拉动力学模型并进行求解，得到晶体的实时半径和晶体高度的变化，然后对比实际半径与晶体参考半径之间的偏差，利用PID控制算法对晶体实际输出半径与高度进行控制，具体如下：

晶体提拉动力学模型如下式所示：

式中，r(t)、z(t)分别表示t时刻晶体的实时半径和实时高度，V_c(t)、ρ分别表示晶体的体积增长速率、坩埚半径、晶体体积和晶体密度，F_ext(t)表示晶体的牵引力，也是晶体直径控制的控制量，根据式(4)得出晶体的实时半径r(t)和实时高度z(t)，将实时半径与步骤1中设置的晶体参考半径r_d(t)比较，根据两者之间的偏差e(t)＝r(t)-r_d(t)，采用PID算法对控制量F_ext(t)进行调整，调整规则如下：

其中，V_c0为晶体的初始体积，K_P和K_I分别为控制器的比例调节系数和积分调节系数，将式(5)代入式(4)计算出PID控制作用下晶体的实时半径和实时高度，为方便计算，假设晶体体积增长速率C_growth为常量，将晶体体积增长速率的变化作为一种扰动d施加在C_growth上；

步骤4、实现对晶体内部温度轨迹的跟踪控制：首先在步骤1建立的有限元模型中，指定监测点P的位置，然后以步骤2得到的加热器功率轨迹作为参考加热器功率，以步骤2得到的最优温度轨迹作为P点的参考温度，以步骤3得到的晶体半径和高度作为动态有限元模型的实时半径，然后利用MPC控制算法对晶体温度进行实时控制；

执行步骤1和步骤2得到整个晶体生长过程中晶体监测点的最优温度轨迹和最优加热器功率轨迹，将其作为真实晶体生长中的设定轨迹；重复步骤3得到晶体生长过程中任意时刻晶体的真实半径，在PID控制器的作用下，晶体半径与设定半径趋于一致；重复步骤4得到晶体生长过程中任意时刻晶体内部的温度分布和最大温度梯度，在MPC控制器的作用下，模拟得到的监测点P处的温度将与设定轨迹趋于一致，同时加热器功率也得到确定，至此，基于有限元数值模拟的晶体温度控制过程结束，具体如下：

步骤4.1、模型预测：建立有限元模型如式(6)所示，根据该模型预测出晶体内部的温度及温度梯度，包括P点的温度：

步骤4.2、滚动优化：使将晶体实际温度X_k与参考温度轨迹之间的误差平方和以及实际加热器功率u_k与参考加热器功率轨迹之间的误差平方求和最小作为优化目标，以监测点温度梯度不超过临界值G_max为约束，同时对加热器功率加以约束，目标函数和约束条件设置如式(7)所示，用优化算法对该优化问题进行求解，优化得到一组最佳的控制量u，即加热器功率：

subjectto:

X_k+1＝A'_kX_k+B'_ku_k+B_ck (7)

y_k＝C'_kX_k

u_min≤u_k≤u_max

|Y_k|≤G_max

式中，A'_k和B'_k分别表示根据步骤3计算出的晶体实时半径和实时高度推导出的有限元时变矩阵，Nc是优化算法的最大迭代次数；

步骤4.3、反馈校正：将每次滚动优化得到的第一个控制量施加在被控对象上，得到被控对象监测点P处的实时温度，将监测点P处的温度X^J反馈到MPC控制器，与MPC的预测值X^Y进行对比，根据两者的偏差对MPC的预测值进行修正，如式(8)所示，其中ε为校正系数：

步骤4.4、再重新进行步骤4.1～步骤4.3，预测模型在每一次迭代过程中都将被更新以更加准确的预测晶体内的温度和温度梯度分布。

根据步骤1，结合晶体生长过程，设定晶体参考形状如图2所示，晶体总高度z_top为20cm，等径半径r_x为5cm，晶体的生长速度为5cm/h，建立如式(2)所示有限元模型，考虑到数值模拟的实时性，本发明建立的数值模型暂且只考虑了晶体部分。加热器位置处于z₁和z₂之间，z₁＝1cm，z₂＝6cm。考虑加热器作用下，设置晶体侧壁z₁～z₂之间为第二类边界条件，设置晶体顶部和底部温度为第一类边界条件，晶体顶部温度T_top设置为676.85℃,底部温度T_f为晶体熔点温度1411.85℃。在提拉动力学模型中，坩埚半径设为12.5cm，体积增长速率C_growth＝0.75cm³/s。一般而言，靠近晶体生长界面处的晶体表面的热应力最大，因此，在选择晶体温度控制监测点时，本发明选择靠近晶体生长界面且位于晶体侧壁处的P点(本发明中选择距离生长界面0.25处的晶体表面的点)为监测点，如图2所示，并将P(r,z)点设置为温度控制系统的监测点。

借助matlab平台对上述有限元模型进行求解,晶体内部温度动态变化的迭代时间间隔为0.8s。

如图3所示，根据步骤2，设置目标函数和相应的约束条件，利用整数规划求解该优化问题,得到加热器功率的最优轨迹u_r和P点温度的最优轨迹T_r。在本发明中加热器功率的调节范围为0～3W/m²，温度梯度在放肩和收尾阶段最大值不得超过140，在等径阶段最大值不得超过50。

根据步骤3，依据提拉动力学模型求解出晶体的实时半径r(t)，根据其与参考半径r_d(t)之间的偏差，利用PID控制器对半径进行实时控制，控制参数设置为K_P＝1.1,K_I＝0.6。得到PID控制器作用下的晶体的半径R和对应的晶体高度分别如图4和图5所示。

根据步骤4，以加热器功率u^r和P点的温度轨迹T^r作为MPC控制器的参考轨迹；根据图4得到的晶体实时半径r建立新的有限元模型，并对晶体半径加入扰动，(扰动大小为0.01～0.09之间的随机数)；利用MPC控制器实现对新的有限元模型中P点温度的轨迹跟踪控制，其中预测步长为10，控制步长为5，式(6)中所示的校正系数为0.01。图6(a)～图6(f)分别展示了晶体从引晶到收尾过程中t＝13.6min，t＝40min，t＝192min，t＝256min，t＝352min和t＝400min时刻晶体内部的温度分布结果。图7显示了从引晶到收尾整个过程中P点的参考温度轨迹和控制作用下的温度轨迹，可以看出，MPC控制器作用下P点的温度较好的跟踪了参考温度轨迹。图8为从引晶到收尾整个过程中监测点温度梯度的变化过程，可以看出，在等径阶段，P点的径向温度梯度几乎为0，在放肩和收尾阶段，温度梯度相对较大，但也满足了本文设定的最大温度梯度约束。

本发明使工艺操作人员不仅可以有效控制晶体尺寸，还可以对晶体内部温度进行有效控制，从而保证晶体品质。在仿真阶段就可以确定晶体生长从放肩到收尾过程中加热器功率参数的轨迹。并且，随着计算机运算速度的提高，可以将有限元模型从目前只考虑了晶体部分扩展到考虑晶体、熔体甚至弯月面部分，使建立的有限元数值模型更加完整，更加贴合真实的生长过程，使优化计算得到的加热器最优功率轨迹能更准确的调整晶体内部的温度分布。所以，本发明克服了传统控制方法无法对晶体内部温度进行控制的不足，弥补了需通过实验调节参数这种实验法来调节晶体品质的高成本和长周期的缺点，提供了一个可以同时控制晶体直径与晶体品质的方法，具有高度的实际应用价值。

Claims (5)

1.一种基于有限元数值模拟的硅单晶生长温度控制方法，其特征在于，具体按照以下步骤实施：

步骤1、参考晶体形状下的有限元建模，设置晶体生长过程中晶体的参考形状和边界条件；

步骤2、以晶体冷却速度最大为目标，以加热器功率的可调范围和晶体内温度梯度的最大值为约束条件，利用优化方法计算出满足约束的晶体温度最优轨迹和加热器功率最优轨迹，并将其作为控制系统的参考输入；

步骤3、对晶体半径进行控制：首先建立晶体生长过程中提拉动力学模型并进行求解，得到晶体的实时半径和晶体高度，然后对比实际半径与晶体参考半径之间的偏差，利用PID控制算法对晶体实际输出半径与高度进行控制；

步骤4、实现对晶体内部温度变化轨迹的跟踪控制：首先在步骤1建立的有限元模型中，指定监测点P的位置，然后以步骤2得到的加热器功率轨迹作为参考加热器功率，以步骤2得到的最优温度轨迹作为P点的参考温度，以步骤3得到的晶体半径和高度作为动态有限元模型的实时半径和实时高度，然后利用MPC控制算法对晶体温度进行实时控制；

执行步骤1和步骤2得到整个晶体生长过程中晶体监测点的最优温度变化轨迹和最优加热器功率轨迹，将其作为真实晶体生长中的设定轨迹；重复步骤3得到晶体生长过程中任意时刻晶体的真实半径，在PID控制器的作用下，晶体半径与设定半径趋于一致；重复步骤4得到晶体生长过程中任意时刻晶体内部的温度分布和最大温度梯度，在MPC控制器的作用下，模拟得到的监测点P处的温度将与设定轨迹趋于一致，同时加热器功率也得到确定，至此，基于有限元数值模拟的晶体温度控制过程结束。

2.根据权利要求1所述的一种基于有限元数值模拟的硅单晶生长温度控制方法，其特征在于，所述步骤1具体按照以下步骤实施：

步骤1.1、晶体热传导方程如下式所示：

式中▽表示圆柱坐标系下的空间梯度，r表示晶体半径，z表示晶体晶体高度，T(r,z,t)表示t时刻晶体中(r,z)处的温度，Pe＝ρC_pv₀R_cruc./k_s为Peclet常数，ρ、v₀、R_cruc、k_s、k_r、C_p分别为密度、提拉速度、坩埚半径、热传导系数、热传导率和比热容，V(r,z,t)为晶体区域的速度，T(r,z,t)为晶体区域的温度；

步骤1.2、用变量r和z分别描述晶体的半径和高度，用表示晶体边界上指向晶体内部的法向量，设置晶体的边界条件如下：

当z＝0时，T|_z＝0＝T_f＝1411.85℃；

当0＜z＜z₁时，

当z₁≤z≤z₂时，

当z₂＜z＜z_top时，

当z＝z_top时，

步骤1.3、根据式(1)用有限元法求解晶体中的温度分布：首先将晶体区域划分成Np个网格单元，然后对公式(1)进行展开，结合设置的对应边界条件，将热传导模型转换为有限元模型，如式(2)所示

式中A_k、B_k是有限元模型中的时变矩阵，是与晶体表面第一类边界条件相关的矩阵，X_k代表晶体内部k时刻所有网格单元的温度矩阵，Y_k代表k时刻所有网格单元的径向温度梯度矩阵，u_k代表k时刻的加热器功率值，C_k是将温度转化为温度梯度的输出矩阵。

3.根据权利要求2所述的一种基于有限元数值模拟的硅单晶生长温度控制方法，其特征在于，所述步骤2具体如下：

目标函数及约束条件如下：

式中u_min、u_max分别为加热器功率可调范围的上下限，G_max为晶体内部的最大温度梯度约束，Q和R均为正定矩阵，利用优化算法对上述优化问题进行求解，得到整个晶体生长过程中加热器的最优功率轨迹和晶体的最优温度轨迹，优化算法的最大迭代次数为Np。

4.根据权利要求3所述的一种基于有限元数值模拟的硅单晶生长温度控制方法，其特征在于，所述步骤3具体如下：

晶体提拉动力学模型如下式所示：

式中，r(t)、z(t)分别表示t时刻晶体的实时半径和实时高度，V_c(t)、ρ分别表示晶体的体积增长速率、坩埚半径、晶体体积和晶体密度，F_ext(t)表示晶体的牵引力，也是晶体直径控制的控制量，根据式(4)得出晶体的实时半径r(t)和实时高度z(t)，将实时半径与步骤1中设置的晶体参考半径r_d(t)比较，根据两者之间的偏差e(t)＝r(t)-r_d(t)，采用PI算法对控制量F_ext(t)进行调整，调整规则如下：

其中，V_c0为晶体的初始体积，K_P和K_I分别为控制器的比例调节系数和积分调节系数，将式(5)代入式(4)计算出PID控制作用下晶体的实时半径和实时高度，为方便计算，假设晶体体积增长速率C_growth为常量，将晶体体积增长速率的变化作为一种扰动d施加在C_growth上。

5.根据权利要求4所述的一种基于有限元数值模拟的硅单晶生长温度控制方法，其特征在于，所述步骤4具体如下：

步骤4.1、模型预测：建立有限元模型如式(6)所示，根据该模型预测出晶体内部的温度及温度梯度，包括P点的温度：

步骤4.2、滚动优化：使将晶体实际温度X_k与参考温度轨迹之间的误差平方和以及实际加热器功率u_k与参考加热器功率轨迹之间的误差平方求和最小作为优化目标，以监测点温度梯度不超过临界值G_max为约束，同时对加热器功率加以约束，目标函数和约束条件设置如式(7)所示，用优化算法对该优化问题进行求解，优化得到一组最佳的控制量u，即加热器功率：

式中，A'_k和B'_k分别表示根据步骤3计算出的晶体实时半径和实时高度推导出的有限元时变矩阵；

步骤4.3、反馈校正：将每次滚动优化得到的第一个控制量施加在被控对象上，得到被控对象监测点P处的实时温度，将监测点P处的温度X^J反馈到MPC控制器，与MPC的预测值X^Y进行对比，根据两者的偏差对MPC的预测值进行修正，如式(8)所示，其中ε为校正系数：

步骤4.4、再重新进行步骤4.1～步骤4.3，预测模型在每一次迭代过程中都将被更新以更加准确的预测晶体内的温度和温度梯度分布。