大尺寸碳化硅晶体生长热-质输运过程建模及数值仿真

卢嘉铮,张 辉,郑丽丽,马 远

(1.清华大学航天航空学院,北京 100084;2.清华大学工程物理系,北京 100084;3.中电化合物半导体有限公司,宁波 315336)

0 引 言

碳化硅(SiC)单晶衬底是制造新一代功率器件、射频器件的重要基础材料,但当前SiC衬底成本偏高,扩大SiC晶体直径并提高晶体品质可增加晶圆利用率,达到降低器件成本目的。目前国内业界主要采用中频感应加热的物理气相传输(physical vapor transport, PVT)法生产100 mm直径SiC单晶[1],150 mm单晶具备批量供应能力[2-3],200 mm衬底已研发成功[3]。国外公司已大量供应200 mm晶片[4-5]。另一方面,由于电阻加热方式可有效控制坩埚内局部温场,电阻加热式PVT系统逐渐崭露头角[6-7]。

PVT法的关键技术之一是设计合理的热-质输运环境,为生长低缺陷晶体提供稳定的热场和均匀的流场。实验表明:晶体生长界面中心低温,晶体表面微凸,晶体质量较高[8];预烧原料可增大粉料颗粒度[9]、减小原料在长晶过程中的形状变化,从而稳定坩埚内温场[10];原料掺铈(CeO2)能有效控制多型生长[11]。PVT坩埚是工作在2 000 ℃高温的封闭结构,实验手段难以监测坩埚内部热-质输运过程。计算机仿真是开展研究的重要办法,在模拟原料多孔结构变化[12-13]、热场设计[14]、预测成核点[15]、计算晶体热应力和位错密度[16-17],以及合成掺矾原料[18]等方面提供了理论指导。但针对电阻炉生长大直径晶体的全过程模拟研究较少,原料消耗、晶体形貌变化等过程的数学模型不清楚,晶体生长过程中特征量之间的相互作用不够清晰。

本团队前期研究了PVT法SiC单晶生长热场设计[14],为本研究提供了有力支撑。本研究针对电阻加热式PVT炉生长150 mm的SiC单晶开展热-质输运过程建模和晶体生长过程的数值模拟研究。首先建立碳化硅原料分解及其多孔结构演变、气相组分输运、能量输运、晶体生长的物理和数学模型,再研究不同原料温度分布(侧面高温、底部高温)对晶体生长形貌变化的影响规律,最后厘清晶体生长形貌变化与原料消耗、温场变化等过程之间的关系。

1 热-质输运数学模型

PVT工艺生长SiC晶体的主要物理过程是:中频感应线圈或电阻加热器产生热源,在坩埚内建立轴向和径向温度梯度,装于坩埚下部的原料被加热后分解,热解产生的气相组分在温度梯度的驱动下被输运至坩埚顶部的低温籽晶面发生再结晶生成单晶体。如图1所示,气相组分输运先后经过原料堆、坩埚上部的空腔区、籽晶面,每个区域内物理过程不同,因此将坩埚内部划分为相应三个区域,并依次建立数学模型。电阻加热式PVT工艺的详细物理过程和热力学过程参见前期工作[14]。

1.1 原料区控制方程

1.1.1 原料热解、碳化和再结晶

原料区是坩埚内装载SiC原料粉末的区域(见图1),一般处于2 500～2 600 K,其最高温度点位于侧面或底部,最低温度点一般在顶部。

高温区域的SiC粉末受热分解生成Si、SiC2、Si2C、SiC等气相组分和疏松石墨,气相组分在温度梯度的驱动下输运至低温的SiC粉末表面和籽晶面再结晶,即:高温的原料处于欠饱和热力学状态,被消耗;低温的原料为过饱和状态,SiC质量增加。

假设:

1)SiC原料为直径0.24 mm的球形颗粒,SiC颗粒石墨化后直径不变[19];

2)原料区为多孔介质区域;

3)有限体积法的网格尺寸大于颗粒尺寸,同一网格体内的颗粒处于同样的热力学状态;

4)球形颗粒呈紧密交错排列;

5)气相组分中Si、C摩尔比为1∶1,即总的SiC气相组分摩尔流量JA=2JSiC2。

在颗粒表面,根据上述物理过程,可用Hertz-Knudsen公式描述原料粉末表面产生的热解气相流量:

(1)

(2)

(3)

式中:ΔA为控制体(网格)内所有SiC颗粒的表面积;ΔV为控制体(网格)体积。ΔA由下式计算:

ΔA=αNpAp

(4)

式中:Np为控制体内颗粒个数;Ap为单个颗粒的表面积;α为修正系数。根据假设4),控制体内颗粒个数为:

(5)

式中:ε0为原料的初始孔隙率;r为颗粒原始半径。由SiC质量守恒与C质量守恒可得:

(6)

与

(7)

ε=1-γSiC-γC

(8)

定义原料消耗程度dg:

(9)

PSi∶PSi2C∶PSiC2≈4∶1∶2

(10)

根据上述组分所涉及反应的化学计量比和假设5),可近似认为:

(11)

(12)

通过式(1)～(7)可知,当颗粒表面SiC气相组分的分压P大于此处的平衡压力P*时,颗粒处于过饱和状态,气相组分在颗粒表面再结晶,SiC体积分数增加;反之SiC颗粒被消耗,SiC体积分数降低,C体积分数增加。由于假设5),异相反应的质量不守恒,这是合理的。原因是:虽然SiC晶体生长一般处于富硅气相,但晶体生长所需Si与C原子数是1∶1,在宏观传质层面,可忽略多余气相Si对晶体生长的影响。

1.1.2 气相组分输运

根据假设1)和2),气相组分在原料中的输运可视作气体在球形颗粒组成的多孔介质中流动和扩散。可用修正后的连续性方程、动量方程和组分方程描述:

(13)

(14)

(15)

式中:ρg为气相平均密度;μ为黏度系数;K为Darcy定律渗透系数(permeability);Deff为有效扩散系数。分别由下式计算:

(16)

(17)

式中:τ为多孔介质的迂曲率[21],其大小与多孔介质基体形状、孔隙率有关;Dij为气体的二元扩散系数。

1.1.3 能量守恒

SiC原料粉末中存在能量输运,颗粒间通过热辐射、热对流和固相导热进行热量传递(见图1),原料区能量守恒方程如下:

(18)

式中:(ρcp)eff是等效热容,为SiC原料、石墨和气相组分热容的体积平均值;keff是有效热导系数,为SiC颗粒导热、疏松石墨导热、气体导热和孔隙间辐射传热的体积平均值。分别由下式计算:

(ρcp)eff=γSiC(ρcp)SiC+γC(ρcp)C+ε(ρcp)g

(19)

(20)

式中:εp为SiC颗粒表面的辐射发射系数;(ρcp)SiC、(ρcp)C和(ρcp)g分别为SiC原料、石墨和气相组分的热容;kSiC、kC和kg分别为SiC原料、石墨和气相组分的热导率。

1.2 生长室中的热-质输运

气相组分从原料表面逸出进入坩埚上部的生长室空腔,在此区域内,气相组分自由流动扩散。忽略气相辐射,坩埚壁、晶体外表面和原料区上表面之间存在辐射传热。坩埚外的气相区域也存在相似物理过程,如电阻加热器与坩埚外壁间、坩埚与保温棉间的辐射传热,坩埚外氩气的流动扩散等。首先,坩埚内除原料区,对气相有质量守恒、动量守恒和组分方程:

(21)

(22)

(23)

对包含保温棉、坩埚和加热器的整个系统(除原料区)有能量守恒:

(24)

(25)

式中:δij是Kronecker delta;Fji为视角系数,表示从j面发射的能量被i面拦截的份额,用下式计算:

(26)

式中:A为微元面i的面积;R为i面与j面心连线长度;θi为R与i面外法向量夹角。参与辐射的面包括晶体表面、生长室坩埚内壁、原料区顶部边界在内的所有气-固边界。

1.3 晶体生长界面

在晶体生长界面上存在Stefan流,组分边界条件采用Hertz-Knudsen模型,对SiC气相组分和氩气分别有:

(27)

JAr=0

(28)

式中:J为摩尔流;ci为组分i的摩尔浓度;vn为法向速度。忽略晶体的径向生长,晶体在法向的生长速率G为:

(29)

式中:Mcrystal为SiC晶体的摩尔质量;ρcrystal为SiC晶体密度。采用动网格技术,在每个时间步前根据上一步计算得到的G值调整晶体生长界面的节点位置,模拟晶体生长界面形状变化和晶体增厚,并重构长晶界面两侧的体网格。

2 数值模拟研究

本节展开对电阻加热式150 mm晶体生长过程的数值模拟研究,探究不同原料温度分布条件对晶体生长的影响规律,梳理晶体生长与原料和热场变化之间的耦合关系。

2.1 算例说明

2.1.1 几何模型与计算说明

计算所用几何模型是对前期工作[14]的改进,如图2所示,系统结构呈圆柱形,半径约570 mm,高约1 600 mm。布置侧面和底部共2个电阻式加热器。

图2 计算域几何结构示意图Fig.2 Schematic diagram of computational domain geometry

为简化计算,采用轴对称计算域,对称轴边界条件为温度0梯度、压力0梯度、轴向速度0梯度和径向速度为0,其余外边界为320 K定温、速度为0、压力0梯度。坩埚内氩气压力300 Pa。采用商业软件ICEM对计算域划分非结构网格,网格尺度为0.5 mm。采用商业软件Fluent对方程进行离散和求解,使用Dynamic Mesh对晶体生长界面网格及其附近体网格进行重构。物性参数和几何尺寸参见表1。

表1 物性和参数[22-24]Table 1 Material properties and parameters for computation[22-24]

2.1.2 算例设置

150 mm晶体生长炉尺寸较大,原料内温度、气流在径向上的不均匀分布影响晶体生长形貌。因此设置2个算例,模拟在典型原料温度分布下的晶体生长过程,探究原料温度对晶体生长形貌变化的影响规律,分析坩埚内关键参量间的耦合关系。算例1仅打开侧面加热器,算例2同时打开侧面和底部加热器,算例模拟时长90 h,算例设计说明如表2所示。

表2 算例设计说明Table 2 Case description

2.2 原料温度分布对晶体形貌影响

2.2.1 侧面高温

算例1模拟的工况是只用侧面加热器对坩埚进行加热,初始时刻坩埚内温场如图3(a)所示,原料最高温度约2 556 K,位于侧面中部。原料顶部和底部是低温区,温度2 475～2 485 K。晶体生长界面中心点温度约2 410 K,边缘约2 440 K。图3(b)展示初始时刻坩埚内的流场,高温区原料热解生成的气相产物在温度梯度的驱动下被输运至低温的原料底部、原料顶部和晶体生长界面。对比流场和多孔结构变化,与文献[19,25]符合度较高。

图3 算例1热/流场结果图。初始时刻坩埚内温度云图(单位:K)(a)和流线与速度矢量图(b);30 h(c)、60 h(d)和90 h(e)的晶体形貌、流场、原料含量和孔隙率,左半图展示流线与SiC原料体积变化,右半图呈现速度矢量与孔隙率(初始值0.65)Fig.3 Heat/flow field results of case 1. At the initial moment, the temperature distribution (unit: K) (a), and streamline and velocity inside the crucible (b); the crystal morphology, flow field, feedstock content and porosity after 30 h (c), 60 h (d) and 90 h (e), the left half of the graph shows the streamlines and the volume change of the SiC source material, and the right half of the graph presents the velocity vector and porosity (initial value 0.65)

从图3可知,在原料表面附近,速度大小沿径向(A至B)逐渐变大,这是由于原料侧面(B至D)温度梯度大,气流驱动力大,A至C点温度梯度小,驱动力小,此现象在后续过程中进一步加剧。SiC气相组分在A点附近再结晶,该区域SiC颗粒体积在90 h增幅约30%,孔隙率从0.65降至0.55,气流通道缩小。而D至B点区域原料不断消耗,孔隙率增加,气流通道逐渐打开,导致A点附近气流速度不断减小,B点附近气流速度不断增大,气流不均匀性凸显。原料底部再结晶程度比顶部更严重,90 h后该区域SiC固相体积增加了70%,造成气源浪费,原料使用率低。

径向分布不均匀的气流量导致晶体中心区域的生长速率比5/8半径处慢,由图3(d)可见生长界面呈“W”形。根据图4所示,平均晶体生长速率先增加后降低,拐点位于40 h附近;A点的生长速率从0.22 mm/h不断降低至0.15 mm/h;前期B点生长速率高于A点,先从0.24 mm/h增至0.25 mm/h,在40 h后迅速降低,最终减至0。因此,中前期B点晶体厚度比A点大,即界面呈“W”形状,后期A点晶体厚度增加较快,最终A/B点晶体厚度相近,生长界面形状平整。

图4 算例1晶体生长速率曲线。(a)晶体生长界面上的A和B点(A点位于界面中心,B点在处);(b)A/B点的晶体生长速率(GA、GB)、晶体厚度(TA、TB)和长晶界面平均生长速率随时间变化图Fig.4 Crystal growth rate curve of case 1. The crystal growth interface taking point A at the center of the interface and point B at 5/8 R; (b) plots of the crystal growth rate (GA, GB), crystal thickness (TA, TB) at point A/B and average growth rate as a function of time

综上,对于原料内最高温度点位于侧面的情况,原料侧面被消耗使气流通道打开,同时原料顶部再结晶使气流通道缩小,造成贴壁附近气流量大、中心区域气流量小,导致中前期籽晶面中心的生长速率低于边缘的生长速率。后期晶体中心生长速率大于晶体边缘生长速率,中心的凹陷被补齐,此现象是PVT法生长晶体的共性问题,在后文2.3小节中会进一步讨论。侧面气流通道的扩张为晶体生长提供较稳定的气源,长晶界面平均生长速率变化较平缓。

2.2.2 底部高温

算例2用侧面加热器和底部加热器同时加热坩埚。通过图5可知,初始时刻,原料底部是高温,顶部是低温,且原料区径向温度梯度小。在此原料温度分布情况下,底部的SiC原料首先热解,气相组分沿负温度梯度方向输运至籽晶面,气流速度在径向上的变化小,气流均匀。值得注意,算例2与算例1的总加热功率相等,但算例2中坩埚内温度更高,因为本研究使用的几何模型中,底部加热器离坩埚更近,且系统底部的保温棉更厚,使坩埚底部受热更充分。

由图5可知,前期晶体从中心区域开始生长,长晶界面中心凸起,边缘平缓。中后期长晶界面形状几乎无变化,各处生长速率相近,晶体表面近似做“平移”运动。该现象与算例1后期的晶体中心区域补长类似。

原料底部的SiC原料持续消耗,无再结晶现象,气相组分从原料区下方流向上方的低温区,并在原料顶部低温区域再结晶。气相组分均匀流出原料表面,气流速度在径向上变化小。但随时间增加,未热解原料减少,且原料顶部气流通道在90 h缩小60%以上,导致气流速度逐渐减小,气流量降低。理论上算例2的原料使用率更高,晶体更厚,晶体形貌更符合需求。

图5 算例2热/流场结果图。初始时刻的坩埚内温场(a)和流场(b);30 h(c)、60 h(d)和90 h(e)的形貌,左半图展示坩埚内SiC原料体积变化量、气流的流线,右半图呈现孔隙率变化量、速度矢量Fig.5 Heat/flow field results of case 2. The temperature field (a) and the flow field (b) in the crucible at the initial moment of case 2. Crystal morphology after 30 h (c), 60 h (d) and 90 h (e), the left half of the figures show the volume change of the SiC source material in the crucible and the streamline of the gas flow, and the right half of the figures show the porosity change and the velocity vector

2.3 晶体生长、原料演变与热场、流场变化等过程间的耦合关系

晶体生长界面的温度分布是影响界面形状的重要因素,由2.2小节已知晶体形貌受原料温度分布、气流均匀性影响,长晶界面形状动态变化是多因素综合作用的结果。本节分别对原料和晶体等特征区域温度、形状结构变化带来的影响展开讨论。

2.3.1 原料区温度和孔隙结构变化的影响

原料区温度是决定晶体生长界面温度的主要因素之一,也是晶体生长的气源。电阻加热器通过热辐射加热坩埚外壁,再由坩埚导热将热量传递给SiC原料,因此原料内最高温度点总是贴近坩埚壁。高温区的原料热解后不断生成疏松石墨,形成较大热阻,以算例2为例,原料区平均温度从初始时刻的2 551 K降低至90 h时的2 522 K,如图6(d)所示。原料温度下降导致晶体生长界面温度降低,加剧了原料低温区的再结晶速度,使原料顶部孔隙率降低、气流通道缩小(见图6(e)),抵达晶体生长界面的气流速度减小(见图5)。

原料顶部再结晶使SiC固相体积分数从0.35增至0.57,增幅60%。顶部区域导热性能增强,A点附近温度梯度减小。原料底部a点附近由于SiC消耗生成疏松石墨,导热性能减弱,该处温度梯度变大。但由于加热器功率、坩埚外环境无变化,a点处温度稳定在2 598 K左右。若用类似算例2温场进行实际长晶试验,上述现象可能造成测量误差,即坩埚底部的温度测量值变化较小,但原料内实际温度可能已低于设计值。

2.3.2 晶体生长界面温度与形貌动态变化

长晶界面与原料表面存在较强辐射传热,面对面辐射传热强度主要受各表面温度和面与面距离影响。晶体增厚,晶体表面与原料表面间距缩短,晶体生长界面温度与原料表面温度、晶体厚度成正相关。以算例2结果为例(见图6(d)),在0～20 h,原料温度略降(2 551～2 548 K),长晶界面均温升高(2 438～2 445 K),即该时段内,晶体增厚对长晶界面温度影响更大。20～40 h,原料温度持续降低,二者对晶体表面温度的影响相互抵消,长晶界面均温保持在2 443 K左右。40 h后,原料温度降低导致长晶界面均温从40 h时的2 443 K降至90 h时的2 424 K。

图6 算例2典型时刻温场和原料体积分数变化。30 h(a)、60 h(b)和90 h(c)坩埚内温度分布;(d)原料区平均温度(Tc)、晶体生长界面平均温度(Ts)和原料底部a点(见图(a))温度(Tc_max)随时间变化情况;(e)原料顶部A点(见图(a))SiC原料的体积分数、孔隙率随时间变化图Fig.6 Typical temperature field and source material volume fraction changes in case 2. The temperature distributions in the crucible at 30 h (a), 60 h (b) and 90 h (c); (d) average temperature in the feed zone (Tc), average temperature of the growth interface (Ts) and the temperature (Tc_max) at point ‘a’ (see Fig(a)) changes with time; (e) volume fraction and porosity of the SiC source material at point A (see Fig(a)) change with time

长晶界面形状主要由界面温度控制[8]。如图7所示,晶体从中心开始生长,中心生长速率大,边缘生长速率小,长晶界面向等温线形状趋近。长晶界面径向温度梯度不断减小,90 h晶体表面形状与2 426 K等温线近乎平行。晶体附近等温线形状是径向温度梯度的体现,主要由系统顶部的散热与加热设计决定。

长晶界面形状在一定程度上受气流均匀性的影响。如图3、图4所示,算例1中气流不均匀,5/8半径处生长速率大于晶体中心,晶体生长界面先呈“W”形。后期温度变化调节生长速率,中心厚度被补上,90 h界面平坦,已消除“W”形。若继续模拟,晶体形貌将与算例2类似。

原料区为晶体生长提供气源。一是原料温度高则热解速率快,晶体生长速率大;二是剩余SiC原料越多,能提供的气源越多,晶体生长越快。由图8可知,晶体生长速率与原料区平均温度、剩余SiC原料质量接近成正比。因此在估算生长速率时,应考虑原料温度和剩余质量带来的影响。

图7 算例2长晶界面上的气相组分浓度与生长速率变化情况。(a)0、20、40、60和80 h晶体生长界面温度分布;(b) 20、40、60和80 h的晶体生长界面形状,与80/90 h时的2 426 K等温线形状;(c)晶体生长界面中心点A与边缘点B位置;(d)A/B点晶体生长速率(GA、GB)、平均生长速率温度(TA、TB)随时间变化折图线;(e)A/B点SiC2气相组分浓度(cA、cB)和平衡浓度随时间变化线图(平衡浓度由平衡分压换算获得,晶体生长速率与SiC2浓度和平衡浓度的差值Δc成正比)Fig.7 Variations of gas species concentration and growth rate at the growth interface in case 2. (a) Temperature distribution at the interface of crystal growth at 0, 20, 40, 60 and 80 h; (b) shape of the growth front at 20, 40, 60 and 80 h, and the shape of the isotherm of 2 426 K; (c) position of the center point A and the edge point B on the crystal growth interface; (d) growth rate (GA, GB) at point A/B, the average growth rate and the temperature(TA, TB) versus time; (e) species concentration (cA, cB) of SiC2 and equilibrium concentrations versus time at point A/B. Equilibrium concentration is obtained by conversion of equilibrium pressure, the crystal growth rate is proportional to the difference (Δc) between the SiC2 concentration and the equilibrium concentration

图8 原料温度及其剩余量与长晶速率关系。(a)算例1、2晶体平均生长速率(G1、G2)随原料平均温度的变化;(b)平均生长速率随剩余SiC原料质量的变化Fig.8 Graph of source temperature and its residual amount versus growth rate. Variation of the average crystal growth rate (G1, G2) with the average temperature of source materials (a) and quality of the remaining SiC source materials (b)

3 结 论

对150 mm晶体生长全过程的数值仿真结果表明,本研究提出的电阻加热PVT法SiC晶体生长热-质输运模型能较准确地模拟原料消耗、再结晶、晶体生长界面形貌和热场变化,流场和晶体生长的仿真趋势与文献吻合。若进一步利用实验数据对此模型进行修正,能广泛用于预测电阻加热式PVT炉中SiC晶体生长形貌、原料和热场变化情况。

仅当侧面加热器工作时,原料区径向温度梯度较大使气流不均匀,气相组分主要沿坩埚壁附近流出原料区,导致长晶界面四周的过饱和度比界面中心大,晶体表面生长成“W”形。但生长后期,晶体中心生长速率变大,凹陷被补上。当底部、侧面加热器同时开启时,能获得轴向温梯大、径向温梯小的原料区温度分布状态,有利于气相组分均匀流动,晶体生长呈表面微凸的理想形状。

晶体表面形状主要受表面的温度控制,晶体表面上各处温度变化引起气相组分平衡浓度变化,进而改变组分的过饱和度,影响生长速率。长晶界面的低温区域逐步升温,生长速率从快到慢,高温区域逐渐降温,生长速率从慢到快,晶体生长界面最终形成等温线形状,此时各处生长速率相近。工程实际中需估算晶体生长速率,晶体平均生长速率大小和原料区温度、剩余原料质量均成正相关。