CN1249272C - 单晶与熔液固液界面形状和单晶点缺陷分布的模拟方法 - Google Patents

Abstract

将单晶提拉装置中的热区模型化成网格结构；将各元件的物理性质数值输入到计算机中；根据加热器的发热量和各元件的热发射率获得各元件的表面温度分布；根据各元件的表面温度分布和热导率获得内部温度分布，从而获得考虑了对流的熔液内部温度分布；根据包括单晶三相点的等温线获得单晶与熔液之间的固液界面的形状；重复第三到第五步骤直到三相点成为单晶的熔点。单晶与熔液之间的固液界面形状的计算结果与实际测量结果彼此吻合得非常好。

Description

单晶与熔液固液界面形状 和单晶点缺陷分布的模拟方法

本发明涉及到对硅之类的单晶与熔液之间固液界面形状以及单晶中点缺陷的分布的计算机模拟方法，所述单晶用切克劳斯基(以下称为CZ)方法拉制。

作为这样一种模拟，如图7所示，已知一种常规方法，它借助于在用CZ方法拉制硅单晶4时，使用总体热传导模拟器，根据提拉装置1中的热区结构和硅单晶4的提拉速度而控制硅熔液2的导热性，来估计硅熔液2的内部温度分布，并利用计算机，从这一内部温度分布获得硅单晶4与硅熔液2之间的固液界面的形状。

还已知有另一种常规方法，用来从所述硅熔液2的内部温度分布获得硅单晶4的各个网格的坐标和温度，然后根据硅单晶4中的填隙硅原子和原子空位的扩散系数和边界条件来求解扩散方程，从而用计算机获得所述填隙硅原子和空位的密度分布。

在这些方法中，热区中的各个元件被网状分割并模型化成网格结构。确切地说，硅熔液2被分割成大约10mm的比较粗略的网格，以便缩短计算时间。

然而，在上述的常规方法中，由于没有考虑实际提拉装置中产生的熔融硅的对流，且熔融硅的网格比较粗略，故存在着固液界面形状的模拟结果与实际测量结果彼此大不相同，且填隙硅原子和空位的密度分布的模拟结果(图6(b))与实际测量结果(图6(e))也彼此差别很大的问题。

本发明的目的是提供一种模拟单晶与熔液之间固液界面形状的方法，其中单晶与熔液之间的固液界面形状的计算结果与实际测量结果彼此吻合得非常好。

本发明的另一个目的是提供一种模拟单晶中的点缺陷的分布的方法，其中单晶中的点缺陷分布的计算结果与实际测量结果彼此吻合得非常好。

如图1和2所示，本发明的第一情况的特征是一种利用计算机模拟单晶与熔液之间固液界面形状的方法，它包含：

第一步骤，将待要计算的单晶14的提拉装置11中的热区模型化成网格结构，

第二步骤，组合热区中各个元件的网格，并将对应于被组合的网格的各个元件的物理性质数值输入到计算机中，

第三步骤，根据加热器的发热量和各个元件的热发射率，获得各个元件的表面温度分布，

第四步骤，借助于根据各个元件的表面温度分布和热导率求解热传导方程，获得各个元件的内部温度分布，然后借助于同时求解在假设熔液12是湍流时得到的湍流模型方程以及Navier-Stokes方程，进一步获得考虑了对流的熔液12的内部温度分布，

第五步骤，根据包括单晶14的三相点S的等温线，获得单晶14与熔液12之间的固液界面的形状，以及

第六步骤，重复第三到第五步骤，直到三相点S成为单晶14的熔点，其中：

沿单晶14的径向紧邻单晶14下方的熔液12的某些或全部网格被设定为0.01-5.00mm，且

沿单晶14的纵向的熔液12的某些或全部网格被设定为0.01-5.00mm。

由于根据本发明第一情况的用来模拟单晶与熔液之间固液界面的方法考虑了熔液12的对流并设定了熔液12的比较精细的网格，故计算得到的单晶14与熔液12之间的固液界面形状，与实际测量结果吻合得非常好。

待要在第二步骤中赋予各个元件的物理性质数值，最好是各个元件的热导率、热发射率、粘滞率、热膨胀系数、密度和比热。

而且，湍流模型方程最好是方程(1)表示的k1模型方程，并使用0.4-0.6范围内的可选数值作为此模型方程的湍流参数C：

${\mathrm{\κ}}_t=\frac{c}{{\mathrm{Pr}}_t}\×\mathrm{\ρ}\×C\×d\sqrt{k}---\left(1\right)$

其中，κt是熔液的湍流热导率，c是熔液的比热，Pr_t是普朗特数，ρ是熔液的密度，d是到盛放熔液的坩埚的内壁的距离，而k是熔液平均流速的可变分量的平方和。

如图2-4所示，本发明的第二情况是一种用计算机来模拟单晶点缺陷的分布的方法，它包含：

第一步骤，在单晶14已经被提拉装置11拉制到特定长度的情况下，将单晶14的提拉装置11中的热区模型化为网格结构，

第二步骤，组合热区中各个元件的各个网格，并将对应于被组合的网格的各个元件的物理性质数值、单晶14的被拉制的长度、以及对应于被拉制的长度的单晶14的提拉速度输入到计算机中，

第三步骤，根据加热器的发热量和各个元件的热发射率，获得各个元件的表面温度分布，

第四步骤，借助于根据各个元件的表面温度分布和热导率求解热传导方程，获得各个元件的内部温度分布，然后借助于同时求解在假设熔液12是湍流时得到的湍流模型方程以及Navier-Stokes方程，进一步获得考虑了对流的熔液12的内部温度分布，

第五步骤，根据包括单晶14的三相点S的等温线，获得单晶14与熔液12之间的固液界面的形状，

第六步骤，重复第三到第五步骤，直到三相点S成为单晶14的熔点，计算提拉装置内部的温度分布，获得单晶14的各个网格的坐标和温度，并将各个数据输入到计算机中，

第七步骤，随着阶段性改变单晶14的拉制长度，重复第一到第六步骤，计算提拉装置11中的温度分布，获得单晶14的各个网格的坐标和温度，并将各个数据输入到计算机中，

第八步骤，将单晶14的各个网格的坐标和温度数据、以及单晶14中的空位和填隙原子的扩散系数和边界条件，输入到计算机中，以及

第九步骤，根据单晶14的各个网格的坐标和温度以及空位和填隙原子的扩散系数和边界条件，求解扩散方程，从而获得单晶14被冷却之后的空位和填隙原子的密度分布。

在根据本发明第二情况的模拟单晶点缺陷的分布的方法中，由于单晶14的内部温度分布是在考虑到熔液12的对流的情况下获得的，而且单晶14中的点缺陷的分布是根据单晶14的这一内部温度分布并在考虑到单晶14中的点缺陷的扩散的情况下获得的，故单晶14中的点缺陷分布的计算数值与实际测量数值彼此吻合得非常好。

图1流程图示出了根据本发明第一实施例的用来模拟硅单晶与硅熔液之间的固液界面形状的方法。

图2是根据本发明第一和第二实施例的具有模型化为网格结构的硅熔液的硅单晶的提拉装置主要部分的剖面图。

图3流程图示出了根据本发明第二实施例的用来模拟硅单晶与硅熔液之间的固液界面形状的方法的第一部分。

图4流程图示出了根据本发明第二实施例的用来模拟硅单晶与硅熔液之间的固液界面形状的方法的第二部分。

图5是主要部分的正视图，示出了各根据实施例1、比较例1、和实际测量结果的硅单晶与硅熔液之间的固液界面的形状。

图6是纵向剖面图，示出了根据实施例2、比较例2(常规例子)、比较例3、比较例4、和实际测量结果的各个硅单晶中的填隙硅原子和空位的分布。

图7是根据常规例子以及比较例1和2的具有模型化为网格结构的硅熔液的硅单晶的提拉装置的主要部分的剖面图。

接着，参照附图来描述本发明的第一实施例。

如图2所示，硅单晶提拉装置11的工作室配备有石英坩埚13以盛放硅熔液12。此坩埚13通过石墨托和支持转轴(未示出)被连接到坩埚驱动装置，而坩埚驱动装置被组成来旋转和升降石英坩埚13。石英坩埚13的外表面被加热器(未示出)围绕，加热器与石英坩埚13保持固定的间隔，而此加热器被隔热筒(未示出)围绕。加热器对置于石英坩埚13中的高纯度多晶硅进行加热并使之熔化，以形成硅熔液12。未示出的圆筒状外壳被连接到工作室的顶部，且此外壳配备有提拉装置。提拉装置被组成来在旋转的情况下提拉硅单晶。

参照图1和2来描述这样组成的硅单晶提拉装置11中的硅单晶14与硅熔液12之间的固液界面形状的模拟方法。

首先，第一步骤将硅单晶提拉装置11的热区中的各个元件分割成网格，并将各个元件模型化为网格结构，所述各个元件是工作室、石英坩埚13、硅熔液12、硅单晶14、石墨托、隔热筒等。具体地说，所述热区中的各个元件的各个网格点的坐标数据被输入到计算机中。此时，沿硅单晶14的径向紧邻硅单晶14下方的硅熔液12的某些或全部网格(以下称为径向网格)，被设定为0.01-5.00mm，最好是0.25-1.00mm。而沿硅单晶14的纵向紧邻硅单晶14下方的硅熔液12的某些或全部网格(以下称为纵向网格)，被设定为0.01-5.00mm，最好是0.1-0.5mm。

将径向网格限制在0.01-5.00mm范围内的理由是，小于0.01mm的范围需要很长的计算时间，而超过5.00mm的范围使得计算不稳定，且即使重复计算也无法确切地确定固液界面的形状。而将纵向网格限制在0.01-5.00mm范围内的理由是，小于0.01mm的范围需要很长的计算时间，而超过5.00mm的范围使得固液界面形状的计算数值与实际测量数值不吻合。在限制某些径向网格在0.01-5.00mm范围内的情况下，靠近硅单晶14外边沿的硅单晶14紧邻下方的硅熔液12，最好被限制在所述范围内，而在限制某些纵向网格在0.01-5.00mm范围内的情况下，靠近硅熔液12表面和底部的硅熔液12，最好被限制在所述范围内。

接着，第二步骤对所述热区中的各个元件的各个网格进行组合，并将被组合的网格的各个元件的物理性质数值输入到计算机中。例如，当工作室由不锈钢制成时，不锈钢的热导率、热发射率、粘滞率、热膨胀系数、密度、和比热被输入到计算机中。并将稍后描述的湍流模型方程(1)的湍流参数C输入到计算机中。

第三步骤根据加热器的发热量和各个元件的热发射率，用计算机获得热区中的各个元件的表面温度分布。亦即，加热器任意设定的发热量被输入到计算机中，并利用计算机从各个元件的热发射率得到各个元件的表面温度分布。接着，第四步骤借助于根据热区中各个元件的表面温度分布和热导率求解热传导方程(2)，获得各个元件的内部温度分布。虽然此处为了简化描述而采用了xyz直角坐标系统，但在实际计算中却采用了柱面坐标系统。

$\mathrm{\ρc}\frac{\∂T}{\∂t}=\frac{\∂}{\∂t}\left({\mathrm{\λ}}_x\frac{\∂T}{\∂x}\right)+\frac{\∂}{\∂y}\left({\mathrm{\λ}}_y\frac{\∂T}{\∂y}\right)+\frac{\∂}{\∂z}\left({\mathrm{\λ}}_z\frac{\∂T}{\∂z}\right)+q\·\·\·\·\·\·\left(2\right)$

其中，p是各个元件的密度，c是各个元件的比热，T是各个元件的各个网格点的绝对温度，t是时间，λ_x、λ_y和λ_z是各个元件的热导率的x、y和z分量，而q是加热器的发热量。

另一方面，至于硅熔液12，用所述热传导方程(2)得到熔液12的内部温度分布，然后根据熔液12的这一内部温度分布，借助于用计算机同时求解在假设硅熔液12是湍流时得到的湍流模型方程(1)以及Navier-Stokes方程(3)-(5)，获得硅熔液12的内部流速分布。

${\mathrm{\κ}}_t=\frac{c}{{\mathrm{Pr}}_t}\×\mathrm{\ρ}\×C\×d\sqrt{k}\·\·\·\·\·\·\left(1\right)$

其中，κ_t是硅熔液12的湍流热导率，c是硅熔液12的比热，Pr_t是普朗特数，ρ是硅熔液12的密度，C是湍流参数，d是到从盛放硅熔液12的坩埚内壁算起的距离，而k是硅熔液12平均流速的可变分量的平方和。

$\frac{\∂u}{\∂t}+u\frac{\∂u}{\∂x}+v\frac{\∂u}{\∂y}+w\frac{\∂u}{\∂z}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂x}+({\mathrm{\υ}}_1+{\mathrm{\υ}}_t)(\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{\∂y^2}+\frac{{\∂}^{2}u}{\∂z^2})+\frac{F_x}{\mathrm{\ρ}}\·\·\·\·\·\·\left(3\right)$

$\frac{\∂v}{\∂t}+u\frac{\∂v}{\∂x}+v\frac{\∂v}{\∂y}+w\frac{\∂v}{\∂z}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂y}+({\mathrm{\υ}}_1+{\mathrm{\υ}}_t)(\frac{{\∂}^{2}v}{{\∂x}^2}+\frac{{\∂}^{2}v}{\∂y^2}+\frac{{\∂}^{2}v}{\∂z^2})+\frac{F_y}{\mathrm{\ρ}}\·\·\·\·\·\·\left(4\right)$

$\frac{\∂w}{\∂t}+u\frac{\∂w}{\∂x}+v\frac{\∂w}{\∂y}+w\frac{\∂w}{\∂z}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂z}+({\mathrm{\υ}}_1+{\mathrm{\υ}}_t)(\frac{{\∂}^{2}w}{{\∂x}^2}+\frac{{\∂}^{2}w}{\∂y^2}+\frac{{\∂}^{2}w}{\∂z^2})+\frac{F_z}{\mathrm{\ρ}}\·\·\·\·\·\·\left(5\right)$

其中u、v和w分别是硅熔液12的各个网格点处流速的x、y和z分量，υ₁是硅熔液12的分子运动粘滞率(物理性质数值)，υ_t是硅熔液12的湍流效应提供的运动粘滞率系数，而F_x、F_y和F_z分别是作用在硅熔液12上的体积力的x、y和z分量。

所述湍流模型方程(1)被称为k1模型方程，且最好使用0.4-0.6范围内的可选数值作为此模型方程的湍流参数C。将湍流参数C限制在0.4-0.6范围内的理由是，小于0.4或大于0.6的范围有以下缺点，会使计算得到的界面形状与实际测量的结果不吻合。所述Navier-Stokes方程(3)-(5)是假设硅熔液12是不可压缩的粘滞率恒定的流体时的运动方程。

借助于用计算机根据上述得到的硅熔液12的内部流速求解热能方程(6)，还获得了考虑到对流的硅熔液12的内部温度分布。

$u\frac{\∂T}{\∂x}+v\frac{\∂T}{\∂y}+w\frac{\∂T}{\∂z}=-\frac{1}{\mathrm{\ρc}}({\mathrm{\κ}}_1+{\mathrm{\κ}}_t)(\frac{{\∂}^{2}T}{{\∂x}^2}+\frac{{\∂}^{2}T}{\∂y^2}+\frac{{\∂}^{2}T}{\∂z^2})\·\·\·\·\·\·\left(6\right)$

其中u、v和w分别是硅熔液12的各个网格点处的流速的x、y和z分量，T是硅熔液12的各个网格点处的绝对温度，ρ是硅熔液12的密度，c是硅熔液12的比热，κ₁是分子热导率(物理性质数值)，而κ_t是用方程(1)计算得到的湍流热导率。

接着，第五步骤根据包括图2的点S所示的硅的三相点S(固液气的三接点)的等温线，用计算机得到硅单晶14与硅熔液12之间的固液界面的形状。而且，改变(逐渐提高)待要输入到计算机中的加热器的发热量，此方法重复所述第三到第五步骤，直至三相点成为单晶14的熔点。

以这种方法得到的硅单晶14与硅熔液12之间的固液界面的形状，与实际测量的结果几乎重合。结果，本发明得到的固液界面的形状能够被用作计算的基础来估计考虑到拉制硅单晶14时点缺陷的扩散的情况下晶体内部的点缺陷的分布。

图2-4示出了本发明的第二实施例。

在此实施例中，以相同于所述第一实施例的方式使用图2所示的硅单晶提拉装置11。

参照图2-4描述了模拟此硅单晶提拉装置11中的硅单晶14中的点缺陷的分布的方法。

首先，在硅单晶14已经被拉制到特定长度L₁(例如100mm)的情况下，第一步骤将硅单晶提拉装置11的热区中的各个元件分割成网格，并将各个元件模型化为网格结构，所述各个元件是工作室、石英坩埚13、硅熔液12、硅单晶14、石墨托、隔热筒等。具体地说，以相同于第一实施例的方式，所述热区中的各个元件的各个网格点的坐标数据被输入到计算机中。此时，以相同于第一实施例的方式，分别设定硅熔液12的径向网格和纵向网格。

第二步骤对所述热区中的各个元件的各个网格进行组合，并将被组合的各个网格的各个元件的物理性质数值输入到计算机中。例如，当工作室由不锈钢制成时，不锈钢的热导率、热发射率、粘滞率、热膨胀系数、密度、和比热被输入到计算机中。并将硅单晶14的已经拉制的长度、对应于这一拉制长度的硅单晶14的提拉速度、以及第一实施例所述的湍流模型方程(1)的湍流参数C，输入到计算机中。

以相同于第一实施例的第三步骤的方式，第三步骤根据加热器的发热量和各个元件的热发射率，用计算机获得热区中的各个元件的表面温度分布，并以相同于第一实施例的第四步骤的方式，第四步骤借助于利用计算机，根据热区中各个元件的表面温度分布和热导率，求解热传导方程(2)，而获得各个元件的内部温度分布。

$\mathrm{\ρc}\frac{\∂T}{\∂t}=\frac{\∂}{\∂t}\left({\mathrm{\λ}}_x\frac{\∂T}{\∂x}\right)+\frac{\∂}{\∂y}\left({\mathrm{\λ}}_y\frac{\∂T}{\∂y}\right)+\frac{\∂}{\∂z}\left({\mathrm{\λ}}_z\frac{\∂T}{\∂z}\right)+q\·\·\·\·\·\·\left(2\right)$

另一方面，至于硅熔液12，以相同于第一实施例的方式，用所述热传导方程(2)得到熔液12的内部温度分布，然后根据熔液12的这一内部温度分布，借助于用计算机同时求解在假设硅熔液12是湍流时得到的湍流模型方程(1)以及Navier-Stokes方程(3)-(5)，获得硅熔液12的内部流速分布。

${\mathrm{\κ}}_t=\frac{c}{{\mathrm{Pr}}_t}\×\mathrm{\ρ}\×C\×d\sqrt{k}\·\·\·\·\·\·\left(1\right)$

$\frac{\∂u}{\∂t}+u\frac{\∂u}{\∂x}+v\frac{\∂u}{\∂y}+w\frac{\∂u}{\∂z}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂x}+({\mathrm{\υ}}_1+{\mathrm{\υ}}_t)(\frac{{\∂}^{2}u}{{\∂x}^2}+\frac{{\∂}^{2}u}{\∂y^2}+\frac{{\∂}^{2}u}{\∂z^2})+\frac{F_x}{\mathrm{\ρ}}\·\·\·\·\·\·\left(3\right)$

$\frac{\∂v}{\∂t}+u\frac{\∂v}{\∂x}+v\frac{\∂v}{\∂y}+w\frac{\∂v}{\∂z}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂y}+({\mathrm{\υ}}_1+{\mathrm{\υ}}_t)(\frac{{\∂}^{2}v}{{\∂x}^2}+\frac{{\∂}^{2}v}{\∂y^2}+\frac{{\∂}^{2}v}{\∂z^2})+\frac{F_y}{\mathrm{\ρ}}\·\·\·\·\·\·\left(4\right)$

$\frac{\∂w}{\∂t}+u\frac{\∂w}{\∂x}+v\frac{\∂w}{\∂y}+w\frac{\∂w}{\∂z}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂z}+({\mathrm{\υ}}_1+{\mathrm{\υ}}_t)(\frac{{\∂}^{2}w}{{\∂x}^2}+\frac{{\∂}^{2}w}{\∂y^2}+\frac{{\∂}^{2}w}{\∂z^2})+\frac{F_z}{\mathrm{\ρ}}\·\·\·\·\·\·\left(5\right)$

借助于用计算机根据上述得到的硅熔液12的内部流速分布求解热能方程(6)，还获得了考虑到硅熔液12的对流的硅熔液12的内部温度分布。

$u\frac{\∂T}{\∂x}+v\frac{\∂T}{\∂y}+w\frac{\∂T}{\∂z}=-\frac{1}{\mathrm{\ρc}}({\mathrm{\κ}}_1+{\mathrm{\κ}}_t)(\frac{{\∂}^{2}T}{{\∂x}^2}+\frac{{\∂}^{2}T}{\∂y^2}+\frac{{\∂}^{2}T}{\∂z^2})\·\·\·\·\·\·\left(6\right)$

接着，以相同于第一实施例的第五步骤的方式，第五步骤根据包括图2的点S所示的硅的三相点S(固液气的三接点)的等温线，用计算机得到单晶14与硅熔液12之间的固液界面的形状。第六步骤改变待要输入到计算机中的加热器的发热量，亦即逐渐提高加热器的发热量；并重复所述第三到第五步骤，直至三相点成为硅单晶14的熔点，然后计算提拉装置11内部的温度分布，获得硅单晶14的各个网格的坐标和温度，并使计算机存储这些数据。

接着，第七步骤使硅单晶14的被拉制长度L₁增加δ(例如50mm)，重复所述第一到第六步骤，然后计算提拉装置11内部的温度分布，得到硅单晶14的各个网格的坐标和温度，并使计算机存储这些数据。此第七步骤一直被执行到硅单晶14的被拉制的长度L₁达到长度L₂(L₂＞L₁)。当硅单晶14的被拉制的长度L₁达到长度L₂时，此方法进行到第八步骤，将单晶14的各个网格的坐标和温度数据、以及硅单晶14中的空位和填隙硅原子的扩散系数和边界条件，输入到计算机中。而且，此方法根据空位和填隙硅原子的这些扩散系数和边界条件，求解扩散方程，从而获得单晶14已经被冷却之后的空位和填隙硅原子的密度分布。

具体地说，填隙硅原子的密度C_i的计算表达式由下列方程(7)表示，而空位的密度C_v的计算表达式由下列方程(8)表示。在方程(7)和(8)中，为了计算密度C_i和C_v随时间变化的状态，假设在晶体的侧面和顶面中以及在固液界面中保持填隙硅原子与空位之间的热平衡。

$C_i^e=K_1\exp (-\frac{E_i}{\mathrm{kT}})\·\·\·\·\·\·\left(7\right)$

$C_v^e=K_2\exp (-\frac{E_v}{\mathrm{kT}})\·\·\·\·\·\·\left(8\right)$

其中K₁和K₂是常数，E_i和E_v分别是填隙硅原子与空位的形成能量，上标e是平衡量，k是波尔兹曼常数，而T是绝对温度。

所述平衡量方程被时间微分，分别成为关于填隙硅原子和空位的下列方程(9)和(10)。

$\frac{{\mathrm{dC}}_i}{\mathrm{dt}}=-\▿(D_i\▿C_i)+\▿(\frac{E_i^t}{{\mathrm{kT}}^2}{D}_i{C}_i\▿T)+k_{\mathrm{iv}}\left(T\right)[C_i{C}_v-C_e^i\left(T\right)C_e^v\left(T\right)]\·\·\·\·\·\·\left(9\right)$

$\frac{{\mathrm{dC}}_v}{\mathrm{dt}}=-\▿(D_v\▿C_v)+\▿(\frac{E_v^t}{{\mathrm{kT}}^2}{D}_v{C}_v\▿T)+k_{\mathrm{iv}}\left(T\right)[C_i{C}_v-C_e^i\left(T\right)C_e^v\left(T\right)]\·\·\·\·\·\·\left(10\right)$

方程(9)和(10)的各个右侧的第一项D_i和D_v分别表示具有方程(11)和(12)所示扩散系数的Fickian扩散。

$D_i=d_i\exp (-\frac{{\mathrm{\ΔE}}_i}{\mathrm{kT}})\·\·\·\·\·\·\left(11\right)$

$D_v=d_v\exp (-\frac{{\mathrm{\ΔE}}_v}{\mathrm{kT}})\·\·\·\·\·\·\left(12\right)$

其中ΔE_i和ΔE_v分别是填隙硅原子和空位的激活能，而d_i和d_v分别是常数。方程(9)和(10)的各个右侧的第二项的E_i ^t和E_v ^t是热扩散提供的填隙硅原子和空位的激活能，而方程(9)和(10)的各个右侧的第三项的各个k_iv分别是成对的填隙硅原子与空位的复合常数。

以这种方式得到的硅单晶14的点缺陷的分布，与实际测量结果几乎重合。结果，有可能在设计提拉装置11阶段，估计提拉装置11拉制的硅单晶14中的点缺陷的分布，且反之，有可能在其设计阶段检验提拉装置11的结构，以便获得被拉制的硅单晶14中的点缺陷的所希望的分布。

在所述第一和第二实施例中，已经使用了硅单晶，但也可以使用GaAs单晶、InP单晶、ZnS单晶或ZnSe单晶。

[实施例]

下面详细地描述本发明的实施例和比较例。

<实施例1>

如图2所示，利用基于图1流程图的模拟方法，获得了从盛放在石英坩埚13中的硅熔液12拉制直径为6英寸的硅单晶14的情况下的硅单晶14与硅熔液12之间的固液界面的形状。亦即，硅单晶提拉装置11中的热区被模型化为网格结构。此处，沿硅单晶14的径向紧邻硅单晶14下方的硅熔液12的网格，被设定为0.75mm，而沿硅单晶14的径向不紧邻硅单晶14下方的硅熔液12的网格，被设定为1-5mm。且沿硅单晶14的纵向的硅熔液12的网格，被设定为0.25-5mm。而且，0.45被用作湍流模型方程的湍流参数C。

<比较例1>

如图7所示，利用常规模拟方法，获得了从盛放在石英坩埚3中的硅熔液2拉制直径为6英寸的硅单晶4的情况下的硅单晶4与硅熔液2之间的固液界面的形状。亦即，硅单晶提拉装置1中的热区被模型化为网格结构。沿硅单晶4的径向的硅熔液2的网格，被设定为10mm，且沿硅单晶4的纵向的硅熔液2的网格，被设定为10mm。没有考虑硅熔液2的对流(未使用湍流模型方程和Navier-Stokes方程的方程组)。除了上述条件之外，用相同于实施例1的方式，用计算机执行模拟。

<比较测试1和评估>

利用各个实施例1和比较例1的模拟方法，获得了硅单晶与硅熔液之间的固液界面的形状。结果示于图5。

如从图5显见，发现用比较例1的模拟方法得到的固液界面的形状(点划线所示)与实际测量结果(实线所示)相差很大，而用实施例1的模拟方法得到的固液界面的形状(虚线所示)几乎与实际测量结果重合。

<实施例2>

如图2所示，利用基于图3和4的流程图的模拟方法，获得了从盛放在石英坩埚13中的硅熔液12拉制直径为6英寸的硅单晶14的情况下的硅单晶14中的点缺陷的分布。

亦即，硅单晶提拉装置11中的热区被模型化为网格结构。沿硅单晶14的径向紧邻硅单晶14下方的硅熔液12的网格，被设定为0.75mm，而沿硅单晶14的径向不紧邻硅单晶14下方的硅熔液12的网格，被设定为1-5mm。且沿硅单晶14的纵向的硅熔液12的网格，被设定为0.25-5mm。而且，0.45被用作湍流模型方程的湍流参数C。在此条件下，获得了考虑到硅熔液12的对流的硅单晶14的内部温度分布，且根据硅单晶14的这一内部温度分布并考虑到硅单晶14中的点缺陷的扩散，获得了硅单晶14中的点缺陷的分布。

<比较例2>

如图7所示，利用常规模拟方法，获得了从盛放在石英坩埚3中的硅熔液2拉制直径为6英寸的硅单晶4的情况下的硅单晶4中的点缺陷的分布。亦即，硅单晶提拉装置1中的热区被模型化为网格结构。此处，沿硅单晶4的径向的硅熔液2的网格，被设定为10mm，且沿硅单晶4的纵向的硅熔液2的网格，被设定为10mm。而且没有考虑硅熔液2的对流(未使用湍流模型方程和Navier-Stokes方程的方程组)。除了上述条件之外，以相同于实施例1的方式，用计算机执行模拟。

<比较例3>

除了考虑到硅熔液的对流但未考虑硅单晶中的点缺陷的扩散之外，以相同于实施例1的方式，用计算机执行模拟。

<比较例4>

除了未考虑硅熔液的对流也未考虑硅单晶中的点缺陷的扩散之外，以相同于实施例1的方式，用计算机执行模拟。

<比较测试2和评估>

利用各个实施例2和比较例2-4的模拟方法，获得了硅单晶中的点缺陷的分布。模拟结果与图6(e)中的硅单晶中点缺陷分布的实际测量结果一起，被示于图6(a)-6(d)。

如从图6显见，发现用各个比较例1-3的模拟方法得到的硅单晶中的点缺陷的分布(图6(b)-6(d))与实际测量结果(图6(e))差别很大，而用实施例1的模拟方法得到的硅单晶中的点缺陷的分布(图6(a))几乎与实际测量结果重合。

如上所述，由于根据本发明的方法将模型化为网格结构的热区中的被组合的网格的各个元件的物理性质数值输入到计算机中，根据加热器的发热量和各个元件的热发射率得到各个元件的表面温度分布，根据各个元件的表面温度分布和热导率得到各个元件的内部温度分布，然后进一步得到考虑对流时的熔液的内部温度分布，根据包括单晶的三相点的等温线得到单晶与熔液之间的固液界面的形状，并重复上述各个步骤直至三相点成为单晶的熔点，并将熔液的网格限制在特定的范围内，故计算得到的单晶与熔液之间的固液界面的形状，与实际测量结果吻合得非常好。结果，用本发明的模拟方法得到的固液界面的形状能够被用作计算基础，以便估计考虑到其在拉制硅单晶时的扩散的硅单晶中的点缺陷的分布。

而且，根据获得考虑到熔液对流的单晶中的内部温度分布以及基于单晶中的这一内部温度分布并考虑到单晶中点缺陷的扩散而获得点缺陷的分布的方法，单晶中的点缺陷的分布的计算结果和实际测量结果彼此吻合得非常好。结果，有可能在设计提拉装置的阶段，对被提拉装置拉制的硅单晶中的点缺陷的分布进行估计，且反之，有可能在其设计阶段，检验提拉装置的结构，以便获得被拉制的硅单晶中的点缺陷的所希望的分布。

Claims (3)

1.一种利用计算机来模拟单晶与熔液之间固液界面形状的方法，它包含：

第一步骤，将待要计算的单晶的提拉装置中的热区模型化成网格结构，

第二步骤，组合所述热区中各个元件的网格，并将对应于所述被组合的网格的所述各个元件的由热导率、热发射率、粘滞率、热膨胀系数、密度和比热组成的物理性质数值输入到计算机中，

第三步骤，根据加热器的发热量和所述各个元件的热发射率，获得所述各个元件的表面温度分布，

第四步骤，借助于根据所述各个元件的表面温度分布和热导率，求解热传导方程，获得所述各个元件的内部温度分布，然后借助于同时求解在假设熔液是湍流时得到的湍流模型方程以及Navier-Stokes方程而进一步获得考虑了对流的熔液的内部温度分布，

第五步骤，根据包括所述单晶的三相点的等温线，获得所述单晶与所述熔液之间的固液界面的形状，以及

第六步骤，重复所述第三到第五步骤，直到所述三相点成为所述单晶的熔点，其中：

沿所述单晶的径向紧邻所述单晶下方的所述熔液的某些或全部网格被设定为0.01-5.00mm，且

沿所述单晶的纵向的所述熔液的某些或全部网格被设定为0.01-5.00mm。

2.根据权利要求1的方法，其中：

湍流模型方程是下列方程(1)表示的k1模型方程，并使用0.4-0.6范围内的可选数值作为此模型方程的湍流参数C：

${\mathrm{\&kappa;}}_t=\frac{c}{P{r}_t}\&times;\mathrm{\&rho;}\&times;C\&times;d\sqrt{k}\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\left(1\right)$

其中，K_t是熔液的湍流热导率，c是熔液的比热，Pr_t是普朗特数，ρ是熔液的密度，d是从盛放熔液的坩埚内壁算起的距离，而k是熔液平均流速的各个可变分量的平方和。

3.一种用计算机来模拟单晶中点缺陷分布的方法，它包含：

第一步骤，在单晶已经被提拉装置拉制到特定长度的情况下，将单晶提拉装置中的热区模型化为网格结构，

第二步骤，组合所述热区中各个元件的网格，并将对应于所述被组合的网格的所述各个元件的物理性质数值、所述单晶的被拉制的长度、以及对应于所述被拉制的长度的所述单晶的提拉速度输入到所述计算机中，

第三步骤，根据加热器的发热量和所述各个元件的热发射率，获得所述各个元件的表面温度分布，

第四步骤，借助于根据所述各个元件的表面温度分布和热导率求解热传导方程，获得所述各个元件的内部温度分布，然后借助于同时求解在假设熔液是湍流时得到的湍流模型方程以及Navier-Stokes方程，而进一步获得考虑了对流的熔液的内部温度分布，

第五步骤，根据包括所述单晶的三相点的等温线，获得所述单晶与所述熔液之间的固液界面的形状，

第六步骤，重复所述第三到第五步骤，直到所述三相点成为所述单晶的熔点，计算所述提拉装置内部的温度分布，获得所述单晶网格的坐标和温度，并将各个数据提供到所述计算机中，

第七步骤，随着阶段性改变所述单晶的被拉制长度，重复所述第一到第六步骤，计算所述提拉装置中的温度分布，获得所述单晶网格的坐标和温度，并将各个数据输入到所述计算机中，

第八步骤，将所述单晶网格的坐标和温度数据、以及所述单晶中的空位和填隙原子的扩散系数和边界条件，输入到所述计算机中，以及

第九步骤，根据所述单晶网格的坐标和温度以及所述空位和所述填隙原子的扩散系数和边界条件，求解扩散方程，从而获得所述单晶被冷却之后的所述空位和所述填隙原子的密度分布。

Images