φ300 mm直拉硅单晶生长过程中的变晶现象及其影响因素

张 晶,刘 丁

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

0 引 言

作为一种可持续发展的绿色行业，光伏产业正受到越来越多的关注。近年来，随着电池生产技术的提升，光伏产业不断朝着降低成本、减少硅材料的损耗和提高生产效率的方向发展。直拉(Czochralski, Cz)硅单晶拥有生产成本低，晶体质量和性能好的优点，目前已广泛应用于光伏产业领域[1-3]。在直拉法制备太阳能级硅单晶中，为满足行业需求，晶体生长热系统的设计朝着提升晶体生长速度和降低加热器功率的方向发展。其中，提高晶体生长速度主要是通过提高硅单晶的提拉速度来实现的。据统计，国际上生长200 mm硅单晶的提拉速度已超过1.8 mm/min，直径为300 mm硅单晶最高提拉速度已经达到1.5 mm/min[4]，硅单晶的生产效率有了明显提高。然而，过高的提拉速度会给晶体生长过程带来问题，主要体现在晶体等径生长阶段[5]。

等径生长阶段硅单晶以均匀直径生长，其生长时间占晶体生长全周期的80%以上，该阶段所生长的晶体直径近似相等、成分均匀、结构完整、缺陷较少，是晶体生长的主要阶段。等径精度越高，晶体的质量越好。然而，在拉晶实验中观察到，当提拉速度较快，等径生长的晶体没有按照圆柱形轴对称的晶格结构生长，出现扭曲生长或者螺旋生长的现象[6]，这种现象称为晶体的“扭晶”。扭晶的出现使晶体的部分晶格在特定方向产生塑性变形，并且变形区原子与未变形区原子在交界处仍维持紧密接触，即生长的晶体由单晶体变为多晶体，发生变晶[6]，极大地降低了晶体品质。由晶体生长实验获得扭晶后的硅单晶棒如图1所示，可以看到晶体直径呈现剧烈变化，光滑的晶体侧面呈现波浪状，晶体生长直径不均匀。在合理范围内调控提拉速度是解决扭晶现象的有效方法，然而降低提拉速度不仅会降低生产效率，增加生产成本，也使晶体内部缺陷分布发生变化，不利于晶体的生长。因此，深入研究扭晶现象的成因并避免这一现象发生，对直拉硅单晶工艺中寻求最大提拉速度并降低成本具有一定指导意义。

与上述扭晶现象类似的情况还发生于化合物晶体生长研究中。文献[7]对不同加热条件下的BGO晶体形状进行研究，结果表明降低熔体自由表面的热辐射有利于改善晶体形状。Okan等[8]解释了DyAlO3晶体在生长过程中比HoAlO3更容易扭晶。文献[9]表明在直拉法生长钙镁锆掺杂钆镓石榴石(SGGG)晶体过程中，提升坩埚位置有助于提高固液界面温度梯度。Kamada等[10]探讨了由于熔体对流不稳定所引起的Cz法生长Pr∶LuAG晶体时非对称温度分布现象。熔体对流的不稳定是影响晶体扭晶的一个原因，而建立熔体稳定性与晶体直径的变化关系十分困难。主要原因是熔体对流的不稳定直接影响熔体内温度分布，进而影响固液界面处的温度梯度，从而影响晶体生长直径的变化。

为解决此问题，本文采用有限元法[11]对硅单晶生长过程进行数值模拟，通过数值模拟和理论相结合的方法，分析提拉速度过大时晶体生长出现扭晶的成因[12]，建立晶体直径与熔体温度分布的关系。通过对影响因素的分析，提出改进工艺参数措施。通过仿真和实验验证该方法的有效性，为在晶体生长热场中设计最大提拉速度以及寻求最佳工艺条件提供参考。

1 数学模型

采用有限元数值模拟的方法研究硅熔体的温度分布和对流情况，是目前分析硅单晶典型生长过程的有效方法[13-16]。在仿真过程中建立单晶炉热场准稳态数学模型，为了方便数值计算，对其进行轴对称简化。对整个单晶炉进行二维(2D)建模，包括熔体、晶体、热屏、气氛、石英坩埚、石墨坩埚、保温套、加热器、炉体等，如图2所示。

(1)

式中：ρl为硅熔体的密度；v表示速度矢量；μ为熔体的粘滞系数；β为温度引起的体膨胀系数；T为温度；下标s表示固体；g为重力加速度；P为压力；式(2)为连续性方程:

(2)

熔体的热传输方程如式(3)所示：

(3)

式中：c为热容量；kl为硅熔体的热传导系数。晶体与熔体固液界面处的热平衡满足能量守恒方程，如公式(4)，表示流出固液界面的热量Q3等于流入界面的热量Q2与界面处相变所释放的潜热Q1之和：

(4)

式中：ks为硅单晶热传导系数；ρs为硅单晶的密度；A为硅单晶横截面积；n为固液界面的法线向量；(∂T/∂n)l和(∂T/∂n)s分别表示相变界面附近熔体的温度梯度和相变界面附近晶体的温度梯度；ΔH为晶体潜热。

所建立的局部三维(3D)模型包括晶体、熔体、石英坩埚和石墨坩埚，对于三维局部模型来说，边界条件的设定影响仿真的精度，其中需要设定的边界条件主要包括：熔体自由表面、晶体外表面、石英坩埚内表面、石英坩埚顶表面、石墨坩埚顶表面。由于晶体的外表面与熔体自由表面温度分布受诸多因素影响，边界条件不确定。原因是晶体外表面和熔体自由表面的温度分布一方面受热辐射影响，另一方面受氩气流动影响[17]。在热系统热屏的作用下，热辐射及氩气流向改变，因此在晶体外表面和熔体自由表面处采用第一类边界条件进行温度约束是不恰当的。为了解决此问题，建立热流密度方程来描述晶体外表面和熔体自由表面由于受到氩气吹拂和辐射产生的热损耗。

对晶体外表面处建立的热流密度边界条件如式(5)所示:

(5)

式中q′s为：

q′s=qout,k-qin,k=σεT4-εqin,k,qin,k=sumj=1～N(Fkjqout,j)

(6)

熔体自由表面处的热损耗如式(7)所示:

(7)

式中q′l为：

q′l=qout,k-qin,k=σεT4-εqin,k
qin,k=sumj=1～N(Fkjqout,j)

(8)

式(5)和式(7)左侧表示热流密度边界条件，δ[T(x)-T0(x)]1.25表示因气体对流而产生的热损耗，∂T/∂r，∂T/∂z为径向温度梯度和轴向温度梯度，q′s表示晶体表面由辐射而产生的热损耗，q′l表示熔体液面由辐射而产生的热损耗，ra和l分别为晶体直径和熔体深度，qout,k是第k个表面上由于热辐射向外散发的热量，qin,k是第k个表面上由于外部热辐射而吸收的热量。式(6)和式(8)与辐射模型相关，在本文中考虑晶体表面与熔体自由表面的自身辐射热流、反射热流和入射辐射热流。δ为对流的热损耗系数，σ为Stefan-Boltzmann常数，ε为辐射系数，Fkj为k、j两个表面之间的角系数[18]。

石英坩埚内表面、石英坩埚顶表面和石墨坩埚顶表面的边界条件同晶体，均采用热流密度边界条件。此外，在仿真过程中，密度采用Boussinesq模型近似，该简化下，控制方程中的密度项只在重力项中考虑密度变化，而其他项中的密度取常数，即平均值ρ0。可近似表示为如公式(9)所示：

Δρ=ρ0(1-βΔT)

(9)

式中：Δρ为密度；β为温度引起的体膨胀系数；ΔT为温差。

在熔体自由表面采用自由滑移壁面，其他壁面均采用非滑移。

方程组(1)～(4)采用有限元方法离散，ELEMENTS数量为8 237，NODES数量为17 993。热场的数值计算以及对固液界面耦合边界处理已发表在文献[19-21]中，压力-速度修正采用SIMPLEC算法。对于几何形状规则的区域采用结构化网格进行离散，对于氩气流动域这样不规则的复杂区域采用三角形非结构化网格离散。determinant方法是有限元分析中评价网格质量常用的方法之一，对其网格划分的雅可比值进行比较，该比值基于雅可比矩阵的行列式。本文采用determinant方法检测网格质量，在分析过程中，采用雅可比矩阵将单元矩阵从理论形状转换为实际形状。理想各单元的雅可比值为1，该值离1越远表示单元网格质量越差。对上述离散的网格进行评价，所得雅可比值最小处为0.885，最大处为0.907，表明网格划分效果较好。

2 结果与讨论

在实际拉晶过程中，提拉速度越大越有利于加快生产进程、提高生产效率。然而，提拉速度受到晶体生长环境的种种约束，由能量守恒定理可知[22]，晶体生长速度与固液界面温度梯度和结晶潜热有关，由式(4)可知晶体生长速度为：

(10)

式中：Gl和Gs分别表示相变界面附近熔体的温度梯度和相变界面附近晶体的温度梯度；由式(10)可知，当Gl=0时，晶体的生长速度可以达到极限值[22]：

(11)

在受控的直拉晶体生长中，提拉速度决定着晶体的轴向生长速度，即提拉速度的快慢决定了晶体生长的快慢[23]。但是由于晶体生长环境热场复杂，熔体处的温度梯度不能满足Gl=0的条件，晶体生长速度不可能取得极限值。并且过大的提拉速度会导致晶体扭晶，由式(10)可知，晶体能够获得的最大生长速度取决于晶体和熔体中的温度梯度的大小，所以对不同提拉速度下熔体内温度分布进行分析，以寻找扭晶出现的原因和最大提拉速度。为了分析这一问题，仿真过程以西安理工大学某型号单晶炉24英寸(1英寸=2.54 cm)热场为模型。在等径阶段，测量得到晶体发生扭晶时的工艺参数l=180 mm，φ=300 mm，其中，l为晶体长度，φ为晶体直径。针对不同提拉速度下的晶体生长过程中的熔体温度、流形分布等参数进行分析，模拟过程所采用的晶体生长过程控制参数为：晶体直径300 mm，晶体长度180 mm，提拉速度0.3～1.2 mm/min，晶体转速4 r/min，坩埚转速8 r/min，氩气流量0.001 m3/s，炉内压力约1 500 Pa。

由图3可知，当提拉速度为0.45 mm/min时，熔体自由表面温度均高于固液界面温度(1 685 K)。随着提拉速度逐渐增大，熔体自由表面的平均温度逐渐降低。当提拉速度为0.5 mm/min时，熔体自由表面靠近三相点处的温度达到晶体的熔点；当提拉速度为1.2 mm/min时，熔体自由表面靠近三相点处温度为1 680 K，说明在熔体自由表面存在比硅熔点(1 685 K)低的温度区域。由式(10)可知，当提拉速度加快时，为了维持晶体生长的相对稳定性，结晶速率也相对提高，导致晶体固液界面附近轴向温度梯度增大。而热系统无法提供足够的热量，于是在熔体的自由表面三相点附近的区域出现了过冷区。出现过冷区的区域开始结晶，影响晶体生长直径变化，而持续的直径变化就会使晶体产生扭晶现象。观察图3固液界面形状，提拉速度由0.45 mm/min升高至1.2 mm/mm的过程中，固液界面形状逐渐凸向晶体。这是由于熔体向晶棒侧传热效率增强，晶体侧壁生长速度快于中心位置，生长界面从凸向熔体变得凸向晶体。

提拉速度分别为0.45 mm/min和1.2 mm/min时，二维仿真和三维仿真下熔体自由表面的径向温度梯度如图4所示。观察图4可得，在同一仿真方法下，不同提拉速度得到的熔体自由表面径向温度梯度分布趋势相似。但是对比同一提拉速度下不同仿真方法得到的结果，当提拉速度为0.45 mm/min时，不论是二维仿真还是三维仿真，熔体自由表面处的温度都大于晶体熔点，表示没有过冷区的产生。但是当提拉速度为1.2 mm/min时，三维仿真结果出现了低于晶体熔点的过冷区，即提拉速度过大时熔体自由表面会产生过冷区。而二维仿真的结果中熔体自由表面没有出现低于晶体熔点的温度，即没有过冷区的产生，与图3中提拉速度为1.2 mm/min时产生过冷区的现象不符。该结果表明了三维仿真的必要性，通过三维仿真，在提拉速度过大时产生的过冷区才可以在熔体自由表面的径向温度梯度分布中被明显地观察到。

根据上述仿真结果，可以根据熔体自由表面的温度分布预测出晶体生长热系统最大提拉速度，将熔体自由表面的温度与提拉速度之间的关系作为判断热系统最大稳定提拉速度的依据。通过该方法可以判断晶体生长到各个阶段所能获得的最大提拉速度，从而避免扭晶现象的发生，设计出更为优良的热系统。

对不同提拉速度下熔体自由表面过冷区和固液界面形状进行仿真，结果如表1所示。由表1可知，随着晶体提拉速度增加，熔体自由表面附近的过冷区逐渐扩大，固液界面形状从凸向熔体逐渐变为凹向熔体。当V>0.45 mm/min时过冷区出现，并随着提拉速度的增加而逐渐扩大，固液界面形状越凹向熔体；当V≤0.45 mm/min时，熔体自由表面不存在过冷区，随着提拉速度的减小固液界面形状越凸向熔体。于是，可以判定在晶体生长等径长度为180 mm时可以获得的最大提拉速度为0.45 mm/min。当提拉速度大于0.45 mm/min时，熔体自由表面产生过冷区，过冷区的持续存在使晶体直径猛然增大，晶体内部温度梯度分布剧烈变化，产生较大的热应力。当热应力超过硅的临界应力后引起位错，导致晶体变晶断苞产生扭晶，等径生长失败。

表1 不同提拉速度下的熔体自由表面过冷区和固液界面形状仿真结果Table 1 Simulation results of free melt surface supercooling zone and the solid-liquid interface shape at different values of pulling speed

根据上述分析，过大的提拉速度使得晶体产生扭晶现象，导致晶体变晶，拉晶失败。为了得到结构完整的合格晶体，提拉速度不能超过各个生长阶段的临界值，该临界值取决于热系统中材料的性质和生长工艺参数。由表1可知，增大晶体提拉速度，固液界面的形状逐渐凸向晶体，即凹向熔体。其原因是根据固液界面能量守恒定理，流入固液界面的热量等于流出固液界面的热量与晶体结晶时释放的潜热之和，如式(4)所示。当提拉速度增大，为了维持晶体的稳定生长，固液界面结晶速率也相应增大，导致单位时间内由于相变产生的结晶潜热增多，进而影响固液界面形状。并且晶体生长过程中固液界面的形状与熔体热传输均匀性密切相关，过大的提拉速度导致熔体温度分布不均匀，也影响固液界面形状，降低生长质量。因此，通过判定固液界面形状是获得晶体生长最大提拉速度的又一判据。合理优化热系统的温度梯度并且获得晶体生长过程中的最大生长速度及其临界值的判断条件，为设计高效热场提供依据，对实际拉晶实验有重大指导意义。

3 改进措施

由上述分析可知，过大的提拉速度产生扭晶，是因为提拉速度过大时熔体自由表面产生了过冷区。当提拉速度较高时，固液界面前端的结晶速率也相应提高，由固液界面热量平衡关系可得，晶体固液界面附近轴向温度梯度增大，而热系统无法提供足够的热量，于是在熔体自由表面三相点附近的区域出现了过冷区，晶体生长固液界面形状由凸向熔体变得凹向熔体，使得晶体扭晶。所以扭晶现象与熔体自由表面三相点附近过冷区的出现和固液界面形状有关，通过改善三相点附近温度分布和固液界面形状，可以避免扭晶的产生。在晶体生长工艺参数中，晶体旋转会产生强迫对流，影响熔体自由表面温度分布和固液界面的温度分布，适当的提高晶体旋转的速度，可以改善三相点附近温度分布以及固液界面形状[24]。

为了验证不同晶体旋转速度对熔体自由表面温度分布和固液界面形状的影响，通过数值模拟得到了不同转速下晶体和坩埚内的温度分布曲线。图5为坩埚旋转速度为0，晶体旋转速度(ωs)分别为0 r/min、5 r/min和10 r/min时晶体和坩埚内的温度分布。观察仿真结果可得，晶体旋转使熔体温度分布逐渐均匀。随着晶体旋转速度的增加，熔体内部等温线向晶体下方移动，熔体内部强迫对流增强，使得晶体下方高温熔体运动增强，熔体内部温度分布更加均匀。并且随着晶体旋转速度的增加，晶体固液界面的形状得到改善，固液界面变得更加平坦，更有利于晶体生长。

在坩埚转速为0，晶体旋转速度分别为0 r/min、5 r/min和10 r/min时固液界面附近硅熔体径向温度分布曲线如图6所示。由图6可以观察到，不论晶体旋转速度多大，固液界面附近熔体温度分布呈现两侧温度高，中心轴处温度低的现象，符合晶体生长工艺。在晶体旋转速度ωs=0 r/min时，相变界面附近熔体温度低于ωs=5 r/min和ωs=10 r/min。在晶体旋转速度ωs=5 r/min时，强迫对流使得晶体下方熔体运动，从而实现热量传递，使得相变界面附近熔体温度升高。在晶体旋转速度ωs=10 r/min时，熔体温度升高更为明显，固液界面相变温度场温度分布均匀性得到改善。

基于上述仿真结果，增加晶体旋转速度可以增强晶体下方熔体流动，改善熔体内温度分布和固液界面形状，使熔体自由表面温度上升，改善扭晶产生的过冷区现象。在提拉速度为0.8 mm/min时，分析不同晶体旋转速度下熔体中的对流形态及温度分布，如图7所示。由图3分析，当提拉速度为1.2 mm/min时，熔体自由表面出现了低于晶体熔点1 685 K的过冷区，该过冷区的温度为1 680 K，与晶体熔点温差为5 K，所以图7显示了等温线间隔为5 K左右的温度分布情况。观察可知，随着晶体旋转速度的增加，固液界面形状更加凸向晶体，熔体自由表面的过冷区有所改善。

由图7可知，晶体旋转能够控制固液界面的形状，当晶体旋转速度为8 r/min时，如图7(a)所示，与图3(c)中晶体旋转速度为4 r/min进行对比，固液界面形状平坦，是因为晶体旋转速度较低时，加热器从坩埚外部加热，固液界面上的温度梯度呈现出中间低，四周高的温度分布。当晶体旋转速度增加，晶体下方传热加快，使得固液界面下方的温度梯度增加，弥补了中间区域较低的温度梯度。当晶体旋转速度为12 r/min时，固液界面形状过于凸向晶体，容易产生较大的热应力，不利于晶体生长，如图7(b)所示。因此，在提拉速度过大时，适当增加晶体旋转速度可以使得固液界面稳定生长，避免熔体自由表面产生过冷区和固液界面形状发生变化导致晶体扭晶，有利于硅单晶生长。

按照上述理论分析和有限元数值模拟结果，在晶体生长国家及地方联合工程中心进行了12英寸硅单晶生长实验，图8(a)所示为12英寸直拉硅单晶生长炉实验平台。为了验证提拉速度、晶体旋转速度对直拉硅单晶生长的影响，设置其余工艺参数相同，在不同的晶体长度、提拉速度和晶体旋转速度下，比较了6组晶体生长。其中，L是晶体长度，V是提拉速度，ωs是晶体旋转速度。实验中设置的工艺参数如表2所示。

表2 拉晶实验参数Table 2 Experimental parameters of pulling crystal

通过上述方法，获得拉晶实验中直径均匀的晶体，如图8(b)所示。结果表明，较大的提拉速度下适当的提高晶体旋转速度可维持晶体稳定生长，获得直径均匀的硅单晶。

4 结 论

本文针对直拉硅单晶实验中提拉速度过大时观察到的扭晶现象，通过有限元数值模拟得到熔体的温度分布和热分布，并结合晶体生长理论分析了扭晶的成因，提出一种最大提拉速度估计的方法，并且通过实验验证了该方法的准确性。过高的提拉速度会影响熔体内温度分布的均匀性，主要表现为熔体自由表面出现过冷区。随着提拉速度的增加，熔体自由表面过冷区面积逐渐扩大，导致晶体生长直径不均匀产生扭晶现象，极大地降低了晶体生长品质。基于理论分析以及仿真和多次拉晶实验验证，提出了提高晶体旋转速度的工艺措施，在维持较高的提拉速度的同时能够有效避免扭晶现象，从而获得直径均匀的硅单晶，提高晶体生长效率。这一方法的提出对实际晶体生长工艺参数设定具有一定的指导作用。