张妮 刘丁 冯雪亮

1)(西安理工大学,晶体生长设备及系统集成国家地方联合工程研究中心,西安 710048)2)(陕西省复杂系统控制与智能信息处理重点实验室,西安 710048)(2018年2月7日收到;2018年7月28日收到修改稿)

为改善晶体相变界面形态,提高晶体品质,提出了一种融合浸入边界法(immersed boundary method,IBM)和格子Boltzmann法(lattice Boltzmann method,LBM)的二维轴对称浸入边界热格子Boltzmann模型来研究直拉法硅单晶生长中的相变问题.将相变界面视为浸没边界,用拉格朗日节点显式追踪相变界面;用LBM求解熔体中的流场和温度分布;用有限差分法求解晶体中的温度分布.实现了基于IB-LBM的动边界晶体生长过程研究.得到了不同晶体生长工艺参数作用下的相变界面,并用相变界面位置偏差绝对值的均值和偏差的标准差来衡量界面的平坦度,得到平坦相变界面对应工艺参数的调整方法.研究表明,相变过程与晶体提拉速度、晶体旋转参数和坩埚旋转参数的相互作用有关,合理地配置晶体旋转参数和坩埚旋转参数的比值,能够得到平坦的相变界面.

1 引 言

直拉法硅单晶生长中相变界面的形状会直接影响晶体中的位错密度大小以及剖面上电阻率的均匀性,因而一直受到学术界和产业界的关注[1].依据能量守恒原理可知相变界面形态由其上下两侧的温度梯度决定,同时相变界面的动态变化也会影响熔体中的温度场和速度场.相变与传热和流动之间的双向耦合作用及其所具有的不确定性使晶体生长过程中相变问题的求解变得非常复杂和困难.因此,通过建立硅单晶生长模型,寻求复杂相变问题的求解方法,对研究硅单晶生长过程中不同工艺参数作用下的相变界面具有实际意义.

国内外学者对晶体生长过程中的相变、流动和传热问题进行了大量的研究,并取得了相应的研究成果[2−6].在流动与传热方面,格子Boltzmann法以其物理意义清晰、程序简单易行等优点成为一种新兴的研究方法并被广泛使用.与传统的流体动力学方法相比,格子Boltzmann法由于不需要离散Navier-Stokes(NS)方程从而避免了对流项离散带来的数值不稳定问题[7].因此,针对流动与传热问题,本文采用格子Boltzmann方法求解NS方程.相变研究的核心主要在于相变界面的追踪.在此方面,文献[8]采用自适应网格法研究了自然对流作用下的相变过程.通过不断建立贴体网格显式追踪相变界面,清晰地描述了相变的物理过程,但迭代过程时间成本高;文献[9,10]用相场法分别研究了合金的凝固和镓的融化过程.通过引入相场变量来区分液相和固相,即隐式追踪相变界面.该方法用欧拉网格描述模拟区域,不需要建立贴体网格,从而节省了由于相变界面运动引起的网格重新划分所需的时间,大大提高了计算效率.然而,文献所描述的相变界面为糊状区域,只有当网格划分极为精细时,才能得到准确的相变界面;文献[11,12]用热焓法分别研究了强迫对流作用下水滴的凝固过程和自然对流作用下固体的融化过程.通过引入新变量——焓,来计算网格点的液相分数以判断网格点的相态,不需要建立贴体网格,具有计算简单且不需要显式追踪相变界面的优点,因此在相变研究领域得到广泛应用.但也存在和相场法同样的问题,即相变界面为糊状区域.尤其当固体区域在外力作用下运动时,相变界面的糊状区域可能会扩大,难以准确地获得相变界面[13].文献[14]用浸入边界法求解流固耦合问题,该方法采用了两套网格:欧拉网格和拉格朗日网格,欧拉网格只需要划分一次,用欧拉网格描述模拟区域,用拉格朗日节点描述和显式追踪固-液界面.文献[13]用浸入边界法模拟了自然对流作用下移动固体的融化过程,结果表明,在解决具有移动边界的相变问题时,该方法可以准确地跟踪相变界面,不存在相场法和热焓法中的糊状边界问题.在计算精度方面,文献[13]的结果显示浸入边界法的计算精度高于热焓法.在计算速度方面,文献[8]的结果显示,用自适应网格法模拟固-液相变过程时,网格生成所耗费的时间占总计算时间的20%—30%.浸入边界法由于不需要建立贴体网格从而节省了大量的计算时间.因此,浸入边界法以网格生成简单、计算速度快且精度高的特点,在固-液相变领域展现出一定的优越性.目前,尚未见文献报道如何使用该方法分析研究硅单晶生长这种不仅具有自然对流而且具有强迫对流的相变问题.本文采用浸入边界法结合格子Boltzmann方法研究了晶体生长中的固-液相变和流动传热过程,并分析了工艺参数与相变界面形态之间的关系.

针对二维轴对称硅单晶生长过程,本文提出了一种融合浸入边界法和格子Boltzmann法的二维轴对称浸入边界热格子Boltzmann模型,结合有限差分法,对伴有相变问题的硅单晶生长过程进行了研究.将相变界面视为浸入边界,通过计算反馈力的形式修正流体边界的速度与温度.同时,根据界面能量守恒方程,建立相变界面运动模型,用拉格朗日节点显式追踪相变界面的位置;采用D2Q9模型构建熔体的密度演化方程、旋转速度演化方程和温度演化方程,用有限差分法求解晶体热传导方程,最后对晶体生长中的相变、流动和传热进行了深入的分析,得到直拉法硅单晶生长中工艺参数的选择和调整方法.

2 控制方程

三维硅单晶生长简化模型如图1(a)所示,在稳态及准稳态情况下,可将三维晶体生长模型简化为二维轴对称晶体生长模型如图1(b)所示.固-液相变界面处温度恒等于相变温度Tm,晶体半径、坩埚半径和坩埚高度之间满足Rx=0.5Rc=0.5H的几何尺寸关系.

图1 直拉法硅单晶生长模型(a)及边界条件(b)Fig.1.Model(a)and boundary conditions(b)of silicon single crystal in Czochralski method.

2.1 熔体对流与传热控制方程

为了用格子Boltzmann方法求解熔体中的流速分布和温度分布,将圆柱坐标系下的晶体生长控制方程转换到伪笛卡尔坐标系下,如下式所示:

其中(u,v)代表熔体流速在x和y方向的分量;w代表熔体绕y轴旋转的速度;υ,ρ,α,g和β分别表示硅熔体的运动黏度、密度、热扩散系数、重力加速度和热膨胀系数.(2)和(3)式中fb=(fbx,fby)是相变界面产生的体积力项,具体形式见第3.2节;fg=gβ(T−T0)是由温差引起的热浮力.(5)式中gb是相变界面产生的能量力项,具体形式见第3.2节.熔体控制方程中fb和gb的引入,体现了相变界面对熔体流动和传热的作用.

2.2 相变控制方程

晶体生长中的相变界面即熔体与晶体的交界面,如(6)式所示的能量守恒方程,在单位时间内相变界面传递至晶体域的热量等于通过熔体传递至相变界面的热量和相变过程中释放的潜热之和:

λl(∇T)l+ρsL(Us−Wpull)·n=λs(∇T)s,(6)其中λs,λl分别为硅晶体热传导系数和硅熔体热传导系数;L为相变潜热;Us为相变界面的运动速度;Wpull为晶体的提拉速度;(∇T)s和(∇T)l分别为晶体侧和熔体侧的温度梯度.可以看出,相变界面的运动速度由晶体侧的温度梯度和熔体侧的温度梯度共同决定.

2.3 晶体热传导控制方程

晶体域只存在由温度梯度产生的能量变化,满足热传导方程,如下式所示:

3 基于浸入边界法的轴对称不可压热格子Boltzmann模型

3.1 轴对称不可压热格子模型

格子Boltzmann方法通过求解介观上描述熔体特征的演化方程来替代求解宏观上复杂的NS方程.利用离散粒子的碰撞交换粒子节点之间的能量信息,利用离散粒子的迁移过程完成粒子节点在时间上的演化.本文采用数值稳定性和计算精度较高的D2Q9(2代表二维,9代表9个离散速度)模型[15]构建三个用来描述熔体的密度、旋转速度和温度的演化方程.

D2Q9模型中设置网格步长δx=1,时间步长δt=1,因此格子速度c=δx/δt≡1,9个方向的权值分别为:ω0=4/9,ω1,2,3,4=1/9,ω5,6,7,8=1/36.离散速度定义如下:

描述熔体密度分布的演化方程为

描述熔体旋转速度的演化方程为

描述熔体温度分布的演化方程为

这里τf,τh,τg代表无量纲松弛时间,满足τf=3υ+0.5,τh=3υ+0.5,τg=3α+0.5;fi(x,t),hi(x,t)和gi(x,t)分别为t时刻格子点x=(x,y)处沿i方向的密度分布函数、旋转速度分布函数和温度分布函数;ψ(T,u,i)分别代表平衡态密度分布函数、平衡态旋转速度分布函数和平衡态温度分布函数,在D2Q9模型中,平衡态分布函数ψ(λ,u,i)=ωiλ[1+3ei·u/c2+4.5(ei·u)2/c4−1.5u2/c2].

流体的宏观参数由下式确定:

3.2 相变界面作用力

在硅单晶生长过程中,相变界面不断地吸收或释放热量,会对其周围流体的速度和温度分布产生影响.动量方程式(2)和(3)中fb项体现了相变界面对周围流体产生的力作用,能量方程式(5)中gb项体现了相变界面对周围流体产生的热作用.体动量力fb和体能量力gb可根据沿浸没边界上的表面动量力Fs(sk,t)和表面能量力Gs(sk,t)求解,具体形式如下:

式中sk表示描述浸入边界的第k个拉格朗日节点;X(sk,t+δt)表示第k个拉格朗日节点在t+δt时刻的欧拉坐标;是与狄拉克函数相关的平滑函数.狄拉克函数为

其中浸没边界的临时速度Us和临时温度Ts可通过边界周围流体节点的速度以及温度信息求解得到,如下式所示:

将(18)和(19)式代入(17)式中,便可得到相变界面对其周围流体节点产生的体积力项fb和能量力项gb.最后依据(21)式对流体节点的信息进行修正:

其中u∗(x,t+δt)和T∗(x,t+δt)是未考虑相变作用时,由格子Boltzmann演化方程计算得到,如(16)式所示.可以看出用二维轴对称浸入边界热格子Boltzmann模型求解伴有相变的流动与传热问题时,流体粒子的信息由标准格子Boltzmann演化方程和相变界面的作用力共同决定.

3.3 边界处理

直拉法硅单晶生长存在四种边界,如图1(b)所示.本文采用非平衡态外推格式[16]处理坩埚壁,用对称格式处理y轴和熔体自由表面.硅单晶生长界面随时间推移不断发生变化,属于动态曲边界问题,且边界位置决定了晶体和熔体的几何区域,其处理过程非常复杂.

从相变界面控制方程式(6)可以看出,相变界面的运动由其上下两侧的温度梯度决定,而相变界面处温度梯度的求解相当麻烦,但每一时间步长格子点释放的热量与界面能量力的积分相等[13],即

结合(6)式,可推导出相变界面的运动速度,

这种处理方式避免了相变界面处温度梯度的求解.这里Cps是常压下的晶体比热;n=(nx,ny)是相变界面的运动法向量,如图2所示,

因此,相变界面可依据上述理论进行动态演化.在演化过程中,一部分节点从流体状态转变为固体状态,一部分节点由固体状态转变为流体状态(如图3中A点).对于前一种情况,只需将新产生的固体节点的分布函数置零;而对于新产生的流体节点,就必须进行分布函数的重构,本文采用邻近插值法进行重构[17].另外,相变界面边界处的流体粒子如图3中B点,存在未知分布函数,本文采用Guo等[18]提出的非平衡态外推法重构曲边界上流体粒子的未知分布函数.

图2 相变界面法向量Fig.2.Normal vector of the interface.

图3 曲边界条件Fig.3.Curve boundary sketch.

3.4 无量纲参数

针对直拉法硅单晶生长过程,为了计算方便,本文将所有的有量纲参数进行无量纲化处理.涉及的无量纲参数有傅里叶数Fo,格拉斯霍夫数Gr,晶转雷诺数Rex,埚转雷诺数Rec,斯蒂芬数Ste和普朗特数Pr.它们的定义如下:

其中Fo数是用来描述非稳态热传导即分子扩散的无量纲数,可视作无量纲的时间;Gr数等于作用在流体上的浮力与黏性力之比,反映自然对流程度;Re数是判别黏性流体流动状态的无量纲数;Ste数为相变贮能材料固相显热与相变潜热之比;Pr数是由流体物性参数组成的无量纲数,反映流体物理性质对流动传热过程的影响;t为非稳态导热过程所经历的时间;Ωx,Ωc分别代表晶体旋转角速度和坩埚旋转角速度.在后面的计算中,用(Th−Tb),Rc和υ/Rc分别对温度、长度和流体速度进行无量纲化处理;用速度的均方根误差作为程序收敛条件,

4 算例验证

为了验证本文提出的融合了浸入边界法和格子Boltzmann法的模型在伴有流动的相变问题中的正确性,本文以固-液相变基准测试问题——方腔熔化进行验证.方腔内部初始化为固体状态,设置左壁面恒为高温Th,右壁面恒为低温Tb,上下壁面均为绝热状态,无量纲参数Pr=0.02,Ra=2.5×104(Ra=Gr·Pr),Ste=0.01.模拟得到不同Fo数下方腔内温度及流线分布,结果如图4所示.

结果显示,随着Fo数的增大,流体区域不断变大,熔化程度增强;同时,由于热浮力的作用,方腔上部的熔化速度大于底部的熔化速度.将不同Fo数下相变界面提取出来,并与热焓法[19]及自适应网格法[8]得到的结果进行对比,获得了良好的一致性,如图5所示,表明将本文模型应用在伴有流动的相变过程中是可行的.

图4 不同Fo数下二维方腔熔化形态(温度分布和流线分布) (a)Fo=4.0;(b)Fo=10.0;(c)Fo=20.0;(d)Fo=30.0Fig.4.Configuration of two-dimensional melting in square cavity(temperature distribution and streamlines):(a)Fo=4.0;(b)Fo=10.0;(c)Fo=20.0;(d)Fo=30.0.

图5 不同相变处理方法下相变界面的比较Fig.5.Comparisons of the location of solid-liquid interface among different methods.

5 仿真研究与分析

首先,以二维轴对称晶体生长模型为研究对象,验证网格独立性,确定最优网格数.模型边界条件如图1(b)所示,无量纲参数Pr=0.013,Gr=2.5×106,Ste=0.01.表1为稳态情况下不同网格数对应的熔体最大、最小流函数值,可以看出,当网格大小为90×90时,流函数最大值几乎不会因为网格大小而发生变化,因此本文选用最优网格数90×90进行模拟.

表1 不同网格数下熔体流函数的最大值和最小值Table 1.Maximum and minimum values of flow function of melt under different grids.

5.1 提拉速度对相变界面的影响

在Gr=2.5×106且不考虑晶体和坩埚旋转的情况下,模拟不同晶体提拉速度Wpull作用下晶体生长中的相变问题,为了清晰地对比温度分布与流动结构,对所得结果关于y轴做镜像显示,结果如图6所示.

由图6可知,提拉速度对熔体流动结构和温度分布并未产生明显的作用,但相变界面发生了较大的变化,且由图7可明显地看到,当Wpull增大至0.0003时,相变界面凸向熔体的情况得到了很大程度的改善.

图6 不同Wpull作用下晶体和熔体的温度分布(左侧)和流线分布(右侧) (a)Wpull=0;(b)Wpull=0.0001;(c)Wpull=0.0002;(d)Wpull=0.0003Fig.6.Temperature distribution(left)and streamlines(right)of crystal and melt under different Wpullvalues:(a)Wpull=0;(b)Wpull=0.0001;(c)Wpull=0.0002;(d)Wpull=0.0003.

图7 不同Wpull作用下的相变界面形态Fig.7.Phase transition interface with different Wpullvalues.

5.2 晶体旋转对相变界面的影响

晶体旋转是实际晶体生产中为实现热场均匀性而必不可少的工艺手段之一. 本节以Wpull=0.0001,Gr=2.5×106为例,研究不同晶转雷诺数下晶体生长中的相变问题,结果如图8和图9所示.

由图8可知,在不考虑晶体旋转作用,即Rex=0时,熔体内部由热浮力驱动的逆时针涡流占据,热熔体在热浮力的作用下沿坩埚侧壁向上运动,然后沿自由表面向晶体侧运动,带动冷熔体向坩埚底部运动形成回流.当晶体旋转后,在热浮力和晶体旋转的共同作用下,晶体下方产生了两个顺时针的小涡流,且随着晶体旋转雷诺数的增大,小涡流的强度增强.当Rex=4000时,晶体下方的熔体几乎全部由顺时针的涡流控制,该涡流带动坩埚底部的熔体向上运动,促进热量向晶体侧的传递,使得相变界面下方的温度梯度增大,温度分布变得更加均匀.

图9显示,随着晶体旋转参数的增大,相变界面形态从最初的凸向熔体变为部分凸向晶体,表明晶体旋转对改善固-液界面附近热场的均匀性起主要作用,这与文献[20]的结果一致,达到了调节晶体转速来改善相变界面形态的目的.当Rex=4000时,相变界面位置与自由表面(y=1)之间的偏差变小,相变界面形态呈M形状.

图8 不同Rex数下晶体和熔体的温度分布(左侧)和流线分布(右侧) (a)Rex=0;(b)Rex=2000;(c)Rex=3000;(d)Rex=4000Fig.8.Temperature distribution(left)and streamlines(right)of crystal and melt in different Rexnumbers:(a)Rex=0;(b)Rex=2000;(c)Rex=3000;(d)Rex=4000.

图9 不同Rex数下相变界面Fig.9.Phase transition interface in different Rex numbers.

5.3 晶体坩埚反向旋转对相变界面的影响

本节以Wpull=0.0001,Gr=2.5×106,Rex=2000为例,研究坩埚与晶体反向旋转时晶体生长中的相变问题.

图10和图11为坩埚与晶体反向旋转作用下的结果.随着坩埚旋转参数的增大,虽然整个熔体中的温度分布基本没变,但相变界面两侧的局部温度分布发生了非常大的变化,随着坩埚旋转参数的增大,晶体下方的温度梯度增大.从图11可看到,当Rec=−400时,相变界面形变量最小,相变界面最平坦;当Rec=−500时,相变界面呈W形状,这是由于坩埚旋转增大了逆时针涡流的强度,将晶体下方的涡流向晶体侧挤压,因此更多的热量在该涡流的作用下被传递到相变界面处,相变界面向晶体侧移动,最终该涡流处的相变界面凸向晶体.坩埚与晶体的反向旋转使得原本凸向熔体的相变界面逐渐凸向晶体,与文献[4]的结果一致.

图10 不同Rec数下晶体和熔体的温度分布(左侧)和流线分布(右侧) (a)Rec=−200;(b)Rec=−300;(c)Rec=−400;(d)Rec=−500Fig.10.Temperature distribution(left)and streamlines(right)of crystal and melt in different Recnumbers:(a)Rec=−200;(b)Rec=−300;(c)Rec=−400;(d)Rec=−500.

图11 不同Rec数下的相变界面Fig.11.Phase transition interface in different Recnumbers.

5.4 晶体坩埚同向旋转对相变界面的影响

本节以Wpull=0.0001,Gr=2.5×106,Rex=3000为例,研究坩埚与晶体同向旋转时晶体生长中的相变问题.

图12和图13为坩埚与晶体同向旋转作用下的结果.与反向旋转作用时的结果相比(图10),除了热浮力和旋转作用驱动的逆时针涡流外,在坩埚底部中心区域产生了一个顺时针涡流,随着坩埚旋转参数的增大,坩埚底部中心区域的顺时针涡流逐渐变大,并控制了对称轴周围的熔体流动.这是因为坩埚和晶体的同向旋转,带动了大部分熔体绕y轴做类似于刚体的旋转运动,随着坩埚旋转作用的增强,熔体流动速度增大,而坩埚底部中心区域的熔体无法达到坩埚的转速,从而在该区域形成顺时针的涡胞.另外,从图13可见,在Rec数增大的过程中,相变界面形态从最初的平坦变为凸向熔体,最后又出现凸向晶体的趋势,其形态的演化与Rec数的变化不具有明确的规律性.

图12 不同Rec数下晶体和熔体的温度分布(左侧)和流线分布(右侧) (a)Rec=100;(b)Rec=200;(c)Rec=300;(d)Rec=400Fig.12.Temperature distribution(left)and streamlines(right)of crystal and melt in different Recnumbers:(a)Rec=100;(b)Rec=200;(c)Rec=300;(d)Rec=400.

实际晶体生产中,为了减少位错的产生,相变界面应该越平坦越好,而不是凸向熔体或凸向晶体.因为平坦的相变界面表示晶体中径向温度梯度较小,从而保证晶体内部热应力较小,避免由于热应力过大而产生的位错缺陷.以上的分析结果表明,晶体旋转对改善相变界面附近热场的均匀性起主要作用,坩埚旋转作为一种辅助调节手段,可通过合理地配置晶体旋转参数和坩埚旋转参数实现对相变界面形态的改善.

图13 不同Rec数下的相变界面Fig.13.Phase transition interface in different Recnumbers.

5.5 平坦相变界面下晶埚转参数之间的函数关系

在模拟晶体旋转和坩埚旋转参数对相变界面形态的影响时,发现当坩埚与晶体反向旋转时,在一组晶体旋转参数下,总能找到一组相应的坩埚旋转参数使得相变界面形态趋于平坦;当坩埚与晶体同向旋转时,相变界面形态随坩埚旋转参数变化的规律性不明确,同一组晶体旋转参数,可能存在两组坩埚旋转参数使得相变界面变得平坦.因此,考虑到工程应用中也大多采用反向旋转技术,本节仅针对反向旋转的情况讨论.为了找到相变界面平坦时对应的晶体和坩埚旋转参数,根据实际生产中工艺参数的调节范围,本文对不同晶体旋转和坩埚旋转参数下的相变界面形态进行了大量的仿真,并进行了定量分析.

为了准确地衡量相变界面的平坦度,首先定义相变界面位置与熔体自由表面位置(y=1)之间的偏差为ek=X(sk,1)−1.用偏差绝对值的均值和偏差的标准差

来描述相变界面位置的偏差大小和波动程度,其中N是拉格朗日节点的个数,X(sk,1)表示拉格朗日节点sk的纵坐标.均值越小表明相变界面越靠近y=1,标准差越小表明相变界面的波动越小.因此,在实际晶体生长中,要求均值和标准差均越小越好.

表2为反向旋转时不同Rex数和不同Rec数配置下相变界面偏差绝对值的均值和偏差的标准差.将每行结果中均值和标准差最小时对应的Rex数和Rec数提取出来,发现它们之间存在一定的关系,如图14所示.为了定量地衡量Rex数和Rec数之间的关系,定义Re′x=Rex/1000,利用曲线拟合找到图14中Rex数和Rec数之间的函数关系分别如下:

从表2和图14可分析出:当Rex<2000时,随着Rex数的增大,要保持相变界面平坦,Rec数也应该增大;当Rex>2000时,随着Rex数的增大,要保持相变界面平坦,Rec数应该减小.这是因为,当Rex数较小时,晶体旋转主导的顺时针涡流较小,尚无法传递足够的热量到相变界面以实现相变界面向晶体侧移动,相变界面完全凸向熔体.因此,需要较强的坩埚旋转作用,将晶体旋转主导的顺时针涡流向晶体一侧挤压,改变晶体下方的温度梯度,从而改变相变界面形态.当Rex数较大时,晶体旋转主导的顺时针涡流本身就较强,甚至会完全占据晶体下方,直接影响相变界面形态,此时,坩埚旋转作为一种辅助的调节手段,仅需要较小的坩埚旋转参数,对晶体下方顺时针涡流的位置进行微调,就可以很大程度上调节相变界面形态.

另外,图14显示当晶转雷诺数小于2000时,最小均值和最小标准差对应的晶转与晶埚转之比的关系基本相同,因此根据图14中的曲线调节埚转雷诺数就能保证相变界面偏差绝对值的均值最小且偏差的标准差最小.当晶转雷诺数大于2000时,最小均值和最小偏差对应的晶埚转参数之间的关系表现出明显的差异.此时,需要结合对相变界面偏差以及波动程度的不同要求来选择合适的晶埚转参数.图14可作为实际晶体生产中晶体旋转参数和坩埚旋转参数调节的参考.

图14 最佳相变界面对应的晶埚转关系Fig.14.Relationship of Rexand Recin optimal interface.

表2 反向旋转时不同晶-埚转作用下相变界面位置的标准差和相对误差Table 2.Standard deviation and relative error of interface in different Rexnumbers and Recnumbers.

6 结 论

针对直拉法硅单晶生长中的相变、流动和传热问题,本文提出了一种融合浸入边界法和格子Boltzmann法的二维轴对称浸入边界热格子Boltzmann模型.用固-液相变基准测试问题验证了所提模型的正确性.利用该模型分析了晶体生长中不同工艺参数对相变界面形态的影响,结果表明:

1)适当地提高晶体提拉速度能有效改善相变界面凸向熔体的问题;

2)在只有晶体旋转作用时,虽然相变界面偏差较小但界面波动情况未得到改善;

3)在晶体-坩埚协同旋转作用下,无论是反向旋转还是同向旋转,相变界面的偏差和波动均能得到有效调节,且在保持相变界面平坦的情况下,反向旋转时晶体旋转参数与晶埚旋转参数比之间满足一定的函数关系,对调整和优化工艺参数具有指导意义.