CN103668436A - 一种熔体热毛细对流过程模拟预测系统及方法 - Google Patents

Abstract

一种熔体热毛细对流过程模拟预测系统及方法，属于半导体材料技术领域，该系统包括熔体数据采集模块、热毛细对流模拟预测模块、模拟预测结果存储模块和结果图像显示模块；热毛细对流模拟预测模块包括网格划分子模块、热毛细对流计算子模块和界面识别子模块；该方法包括采集熔体液桥区域的几何参数、介质物性参数、熔体所在环境侧流体的物性参数和熔体热毛细对流过程初始参数；进行热毛细对流模拟预测；将速度场、温度场、流场和压力场和两相界面位置分类且按时间顺序存储并以曲线图、矢量图或云图形式分类显示。本发明可用于预测在微重力环境及磁场控制下晶体生长过程中内部热毛细对流发展，温度场和速度场的分布，为制备高质量单晶提供预测手段。

Description

一种熔体热毛细对流过程模拟预测系统及方法

技术领域

本发明属于半导体材料技术领域，具体涉及一种熔体热毛细对流过程模拟预测系统及方法。

背景技术

半导体晶体是半导体工业的主要基础材料，但是，天然的晶体无论在品种、数量和质量方面都远远达不到现代科学技术发展的需要，因此促进了人工晶体的快速发展。晶体材料性能的提高，晶体新品种的开发，在很大程度上依赖于晶体生长法的改进、发展与完善。

熔体生长法是目前工业生产半导体单晶最常用和最重要的生长晶体方法，在熔体生长晶体过程中，先是将固体加热熔化，然后在受控条件下，通过降温使熔体逐渐凝固成晶体。浮区晶体生长方法是熔体中生长晶体的最常用的方法之一。

浮区法晶体生长技术，其生长示意图如图1所示。在浮区法晶体生长中，利用辐射加热、射频加热或电子束加热等方法作为加热源，首先在籽晶与多晶棒之间形成熔区，然后自下而上移动加热源或者自上而下移动熔体，完成晶体生长。处于多晶柱和单晶柱之间一定长度的熔体完全依靠自身表面张力来维持熔区不破裂(如图2所示)，浮区法晶体生长过程中的熔区主要受两种力：表面张力和重力。溶区靠表面张力来维持不破裂，另一方面由于切向表面张力梯度的驱动会在熔体中形成表面张力对流，或者称为热毛细对流。另外，在重力环境下，温度梯度引起熔体密度的不均匀性会在熔体中形成热浮力对流。

熔区中的表面张力流和热浮力流对浮区法生长晶体的质量产生重要作用。它们会影响晶体固化界面附近的温度梯度和杂质分布，在晶体中形成宏观分凝现象。而且随着对流强度的增大，熔体对流将转变为依赖时间的振荡对流，从而导致晶体中化学组分的变化产生微观分凝现象。因此，在浮区法晶体生长过程中，研究熔体中表面张力对流和热浮力流成为提高晶体质量的主要途径。

重力环境下的浮区法晶体生长中，由于温度梯度而引起的热浮力流是形成宏观和微观不均匀性(生长条纹)的主要原因之一。微重力环境下，随着重力逐渐趋于零时，此时由重力引起的热浮力对流就几乎消失，从而就可以避免不稳定浮力流对晶体生长质量的影响。另一方面，重力环境下用浮区法进行晶体的生长需要足够大的表面张力来维持熔区的稳定，从而极大地限制了生长晶体的尺寸。但在空间微重力条件下，熔区中的流体静压力几乎完全消失，从而采用浮区法生长较大尺寸的晶体成为可能。

目前，人们对不同微重力水平下的半导体晶体生长进行的实验研究主要利用弹道火箭飞行、落塔及轨道空间站等手段。相关的实验结果表明，在微重力环境下生长的晶体，其完整性、杂质分布的宏观均匀性和微观均匀性都得到明显改善。但是，现有的微重力条件下的晶体生长有局限性，首先，受到目前航空航天发展水平的制约，空间实验的机会少，并且实验成本高。其次是在微重力环境下的浮区法晶体生长过程中，虽然热浮力对流已经极度地衰减，但是熔体自由表面上由不平衡的表面张力所驱动的表面张力流仍然存在，并在熔体流动中占主导地位。当熔体中的温度差大于一定数值后，熔体表面张力流将会产生失稳而引起晶体生长条纹的产生，不利于高质量晶体的形成，甚至可能导致晶体生长的失败。再次是空间残余重力以及“g-颤动”现象，这些也会给晶体生长带来一些消极影响。因此，人们采用各种手段(外加磁场、表面敷层、气流剪切)来抑制热毛细对流的产生，特别是使用外加磁场来提高硅单晶的质量。

同时，在理论研究方面，虽然有很多关于数值计算和理论分析的报道，但很难发现准确的稳定性分析和数值模拟。一个重要的原因是绝大多数研究采用了简化模型，不考虑熔区液桥自由面变形，将熔区液桥自由面假定为不可变形的柱面或者是由流体体积和Young-Laplace方程确定的不可变形曲面，影响了研究结果的准确性。根本原因是采用数值模拟方法来俘获液桥自由面微小迁移时，现有的各种运动界面跟踪法都不能满足微重力水平下的要求。

发明内容

针对现有技术存在的问题，本发明提供一种熔体热毛细对流过程模拟预测系统及方法。

本发明的技术方案是：

一种熔体热毛细对流过程模拟预测系统，包括：

熔体数据采集模块：用于采集熔体液桥区域的几何参数、介质的物性参数、熔体所在环境的物性参数和熔体热毛细对流过程初始状态参数，并传输至热毛细对流模拟预测模块；

所述熔体液桥区域的几何参数包括：液桥直径、液桥高度、环境测宽度、环境测高度；

所述熔体液桥区域的介质物性参数包括：熔区介质的密度、熔区介质动力粘度、熔区介质表面张力随温度变化系数、熔区介质比热容、熔区介质导热系数和熔区介质热膨胀系数；

所述熔体所在环境的物性参数包括：空气的密度、空气动力粘度、空气的比热容、空气导热系数和环境温度；

所述熔体热毛细对流过程初始状态参数包括：上盘温度、下盘温度和重力加速度；

热毛细对流模拟预测模块：包括网格划分子模块、热毛细对流计算子模块和界面识别子模块；

其中，网格划分子模块用于根据熔体数据采集模块采集的熔体液桥区域的几何参数，确定熔体液桥区域并划分成网格，形成计算区域；

热毛细对流计算子模块：用于对计算区域的速度场、温度场、压力场和流场进行模拟预测，并将结果分别传输至界面识别子模块和模拟预测结果存储模块；

界面识别子模块：用于根据计算区域的速度场、温度场、流场和压力场对两相界面位置进行计算，并将计算结果传输至模拟预测结果存储模块；

模拟预测结果存储模块：用于将计算区域的速度场、温度场、流场和压力场的计算结果和两相界面位置计算结果分类且按时间顺序存储，并传输至结果图像显示模块；

结果图像显示模块：用于将模拟预测结果存储模块传来的数据以曲线图、矢量图或云图的形式分类显示。

所述热毛细对流模拟预测模块还包括匀强磁场施加子模块，用于对网格划分子模块划分网格后形成的计算区域施加匀强磁场，使计算区域的速度场、温度场、流场和压力场发生变化，进而使熔区热毛细对流发生改变。

采用所述的熔体热毛细对流过程模拟预测系统进行熔体热毛细对流过程模拟预测的方法，包括以下步骤：

步骤1：采集熔体液桥区域的几何参数、介质物性参数、熔体所在环境侧流体的物性参数和熔体热毛细对流过程初始参数；

所述熔体液桥区域的几何参数包括：液桥直径、液桥高；

所述熔体液桥区域的介质物性参数包括：熔区介质的密度、熔区介质动力粘度、熔区介质表面张力随温度变化系数、熔区介质比热容、熔区介质导热系数和熔区介质热膨胀系数；

所述熔体所在环境的物性参数包括：空气的密度、空气动力粘度、空气的比热容、空气导热系数和环境温度；

所述熔体热毛细对流过程初始状态参数包括：上盘温度、下盘温度和重力加速度；

步骤2：根据熔体液桥区域的几何参数进行热毛细对流模拟预测；

步骤2.1：根据熔体数据采集模块采集的熔体液桥区域的几何参数，确定熔体液桥区域并划分成网格，形成计算区域；

步骤2.2：对计算区域的速度场、温度场、流场和压力场进行模拟预测；

步骤2.3：采用Level Set Method方法对两相界面位置进行计算；

步骤2.3.1：将当前时刻熔体的运动界面构造成Level Set函数的零等值面，对计算区域的任一点到零等值面的距离进行计算，若距离小于零，则该点位于熔体内部，若距离大于零，则该点位于熔体外部，若距离等于零，则该点为两相界面上的点；

步骤2.3.2：对两相界面位置进行修正：根据下一时刻的熔体面积与起始时刻的熔体面积之差及界面的曲率关系确定修正比例，对两相界面位置按修正比例进行修正；

步骤2.3.3：返回步骤2.3.1，对下一时刻两相界面位置进行模拟计算；

步骤3：将计算区域的速度场、温度场、流场和压力场的计算结果和两相界面位置计算结果分类且按时间顺序存储；

步骤4：将计算区域的速度场、温度场、流场、压力场数据以及两相界面位置，以曲线图、矢量图或云图的形式分类显示。

所述步骤2.1中划分网格后形成的计算区域施加匀强磁场，使计算区域的速度场、温度场、流场和压力场发生变化，进而使热毛细对流发生改变。

有益效果：

本发明提供了熔体热毛细对流过程模拟预测系统及方法，该系统是基于Level Set法的模拟重力和微重力环境下非等温熔区液桥热毛细对流数值模拟系统平台，可用于预测在微重力环境及磁场控制下，晶体生长过程中内部热毛细对流的发展，温度场和速度场的分布。利用本系统进行大量的数值实验，可以从中摸索规律特性，减少空间实验的盲目性。并且，可以建立真实空间浮区法晶体生长环境下的晶体生长模型，为制备高质量单晶提供必要的预测手段。同时也为研究载人空间站环境下非等温液桥热毛细对流振荡特性实验建立前期理论依据，本系统能够满足高分辨率以及俘获任意微小尺度自由面迁移的需要，对于研究大尺度熔区液桥热角区温度、速度和自由面振荡特性都具有建设性的意义。

附图说明

图1是熔体生长晶体的浮区法示意图；

图2是热毛细对流示意图；

图3是本发明具体实施方式的熔体热毛细对流过程模拟预测系统结构示意图；

图4是本发明具体实施方式的熔体热毛细对流过程模拟预测方法流程图；

图5是本发明具体实施方式的交错网格(a)主控制容积(b)u控制容积(c)v控制容积；

图6是本发明具体实施方式的熔区内部流函数等值线随时间的变化示意图；

图7是Level Set Method方法在界面识别计算中介质区域及其边界示意图；

图8是本发明具体实施方式的均匀磁场磁力线示意图，(a)轴向均匀磁场，(b)横向均匀磁场；

图9是本发明具体实施方式的流函数等值线随轴向磁场强度的变化(g_v＝0.6g₀、ΔT＝1355℃)；

图10是本发明具体实施方式的不同轴向磁场强度下液桥表面变形(g_v＝0.6g₀、ΔT＝1355℃)；

图11是本发明具体实施方式的流函数等值线随横向磁场强度的变化(g_v＝0.6g₀、ΔT＝1355℃)；

图12是本发明具体实施方式的不同横向磁场强度下液桥表面变形(g_v＝0.6g₀、ΔT＝1355℃)。

具体实施方式

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

如图3所示，本实施方式的熔体热毛细对流过程模拟预测系统，包括：

熔体数据采集模块：用于采集熔体液桥区域的几何参数、介质的物性参数、熔体所在环境的物性参数和熔体热毛细对流过程初始状态参数，并传输至热毛细对流模拟预测模块；

所述熔体液桥区域的几何参数包括：液桥直径、液桥高度、环境测宽度、环境测高度；

所述熔体液桥区域的介质物性参数包括：熔区介质的密度、熔区介质动力粘度、熔区介质表面张力随温度变化系数、熔区介质比热容、熔区介质导热系数和熔区介质热膨胀系数；

所述熔体所在环境的物性参数包括：空气的密度、空气动力粘度、空气的比热容、空气导热系数和环境温度；

所述熔体热毛细对流过程初始状态参数包括：上盘温度、下盘温度和重力加速度；

热毛细对流模拟预测模块：包括网格划分子模块、热毛细对流计算子模块和界面识别子模块；

其中，网格划分子模块用于根据熔体数据采集模块采集的熔体液桥区域的几何参数，确定熔体液桥区域并划分成网格，形成计算区域；

热毛细对流计算子模块：用于对计算区域的速度场、温度场、压力场和流场进行模拟预测，并将结果分别传输至界面识别子模块和模拟预测结果存储模块；

界面识别子模块：用于根据计算区域的速度场、温度场、流场和压力场对两相界面位置进行计算，并将计算结果传输至模拟预测结果存储模块；

模拟预测结果存储模块：用于将计算区域的速度场、温度场、流场和压力场的计算结果和两相界面位置计算结果分类且按时间顺序存储，并传输至结果图像显示模块；

结果图像显示模块：用于将模拟预测结果存储模块传来的数据以曲线图、矢量图或云图的形式分类显示。

热毛细对流模拟预测模块还包括匀强磁场施加子模块，用于对网格划分子模块划分网格后形成的计算区域施加匀强磁场，使计算区域的速度场、温度场、流场和压力场发生变化，进而使熔区热毛细对流发生改变。

采用上述的熔体热毛细对流过程模拟预测系统进行熔体热毛细对流过程模拟预测的方法，流程如图4所示，包括以下步骤：

步骤1：采集熔体液桥区域的几何参数、介质物性参数、熔体所在环境侧流体的物性参数和熔体热毛细对流过程初始参数；

熔体液桥区域的几何参数包括：液桥直径D＝5mm、液桥高H＝2.5mm，纵横比为A＝1.0(H/R)；

熔体液桥区域的介质物性参数包括：熔区介质的密度ρ_l＝2500kg/m³、熔区介质动力粘度μ_l＝7.8×10^-5kg/(m·s)、熔区介质表面张力随温度变化系数σ_T＝1.87×10^-2、熔区介质比热容C_p＝0.703KJ/(kg·℃)、熔区介质导热系数K＝64W/m·K和熔区介质热膨胀系数1；

介质(熔体)所在环境的物性参数包括：空气的密度ρ_g＝1.205kg/m³、空气动力粘度μ_g＝1.80×10^-5kg/(m·s)、空气的比热容C_p＝1.005KJ/(kg·℃)、空气导热系数K＝0.02624W/m·K和环境温度T_s＝298K；

熔体热毛细对流过程初始状态参数包括：上盘温度T_b＝1685K、下盘温度T_t＝330K和重力加速度g＝9.8m/s²，无量纲时间间隔量t＝0.0005；

步骤2：根据熔体液桥区域的几何参数进行热毛细对流模拟预测；

步骤2.1：根据熔体数据采集模块采集的熔体液桥区域的几何参数，确定熔体液桥区域并划分成网格，形成计算区域，计算区域横向节点数81和轴向节点数41；

步骤2.2：对计算区域的速度场、温度场、流场和压力场进行模拟预测：用SIMPLE算法求解N-S方程，得到计算区域的速度场、流场和压力场；

步骤2-2-1：假设一个压力场，记为p^*；

步骤2-2-2：利用压力场p^*求解动量离散方程，得出该压力场对应的横向速度u^*、纵向速度v^*；

步骤2-2-3：该利用质量守恒方程来改进压力场p^*，使改进后的压力场相对应的速度场能满足连续性方程，用p′，u′和v′分别表示压力场的修正值、横向速度的修正值及纵向速度的修正值；

步骤2-2-4：以p^*+p′、u^*+u′、v^*+v′作为本次迭代的解并据此开始下一次的迭代计算；

以直角坐标系中二维稳态流动为例，控制方程为：

$\frac{\∂\left(\mathrm{\ρu}\right)}{\∂x}+\frac{\∂\left(\mathrm{\ρv}\right)}{\∂y}=0---\left(1\right)$

$\frac{\∂\left(\mathrm{\ρuu}\right)}{\∂x}+\frac{\∂\left(\mathrm{\ρuv}\right)}{\∂y}=\frac{\∂}{\∂x}\left(\mathrm{\μ}\frac{\∂u}{\∂x}\right)+\frac{\∂}{\∂y}\left(\mathrm{\μ}\frac{\∂u}{\∂y}\right)-\frac{\∂p}{\∂x}+S_x---\left(2\right)$

$\frac{\∂\left(\mathrm{\ρvu}\right)}{\∂x}+\frac{\∂\left(\mathrm{\ρvv}\right)}{\∂y}=\frac{\∂}{\∂x}\left(\mathrm{\μ}\frac{\∂v}{\∂x}\right)+\frac{\∂}{\∂y}\left(\mathrm{\μ}\frac{\∂v}{\∂y}\right)-\frac{\∂p}{\∂y}+S_y---\left(3\right)$

其中：ρ密度项，μ粘度项，S_x和S_y分别为源项。

对(1)、(2)、(3)式进行离散化，同时采用交错网格，如图5所示，离散后的连续性方程如下：

(ρuΔy)_e-(ρuΔy)_w+(ρuΔx)_n-(ρuΔx)_s＝0          (4)

离散后的动量方程分别为：

$a_e{u}_e={\mathrm{\Σa}}_{\mathrm{nb}}^u{u}_{\mathrm{nb}}+b^u+(p_P-p_E)A_e---\left(5\right)$

$a_n{u}_n={\mathrm{\Σa}}_{\mathrm{nb}}^v{u}_{\mathrm{nb}}+b^v+(p_P-p_N)A_n---\left(6\right)$

其中，

分别是与u_e、u_n有关的速度修正项；b^u、b^v为不包含压力在内的源项中的常数部分；A_e，A_n对应网格的截面积。

改进后的压力场与速度场也满足本次迭代的动量离散方程，即线性化的动量方程，将u^**＝u^*+u′，v^**＝v^*+v′，p^**＝p^*+p′代入动量方程，可得：

$a_e(u_e^{*}+u_e^{\′})={\mathrm{\Σa}}_{\mathrm{nb}}^u(u_{\mathrm{nb}}^{*}+u_{\mathrm{nb}}^{\′})+b+[(p_P^{*}+p_P^{\′})-(p_E^{*}+p_E^{\′})]A_e---\left(7\right)$

注意到u^*，v^*是根据p^*之值从这一离散方程解出的，因而它们满足：

$a_e{u}_e^{*}={\mathrm{\Σa}}_{\mathrm{nb}}^u{u}_{\mathrm{nb}}^{*}+b+(p_P^{*}-p_E^{*})A_e---\left(8\right)$

假定由源项构成的b的值保持不变，于是将式(6)、式(7)相减得：

$a_e{u}_e^{\′}={\mathrm{\Σa}}_{\mathrm{nb}}^u{u}_{\mathrm{nb}}^{\′}+(p_P^{\′}-p_E^{\′})A_e---\left(9\right)$

如果直接按(9)式来确定速度修正值，将导致十分复杂的计算，因此，在上述两个影响因素中压力修正的直接影响是主要的，四周邻点速度修正值的影响可近似地不予考虑，假设在

中系数

于是得到速度修正方程：

a_eu′_e＝(p′_P-p′_E)A_e              (10)

或

$u_e^{\′}=\frac{A_e}{a_e}(p_P^{\′}-p_E^{\′})=d_e(p_P^{\′}-p_E^{\′})---\left(11\right)$

类似可得

$v_n^{\′}=\frac{A_n}{a_n}(p_P^{\′}-p_N^{\′})=d_n(p_P^{\′}-p_N^{\′})---\left(12\right)$

于是，改进后的速度为：

u＝u^*+d_e(p′_P-p′_E)           (13)

v＝v^*+d_n(p′_P-p′_N)           (14)

而改进后的压力为：

p＝p^*+p′              (15)

压力修正值p′应当满足的条件是：根据p′而改进的速度场能满足连续性方程。把(11)式代人连续性方程的离散形式，即可获得能满足上述条件的p′的代数方程，可得：

a_pp_P＝a_Ep_E+a_Wp_W+a_Np_N+a_Sp_S+b             (16)

其中：a_E＝ρ_ed_eΔy、a_W＝ρ_wd_wΔy、a_N＝ρ_nd_nΔy、a_S＝ρ_sd_sΔy、a_P＝a_E+a_W+a_N+a_S、b＝(ρu^*Δy)_w-(ρu^*Δy)_e+(ρu^*Δy)_s-(ρu^*Δy)_n。

求解能量方程，对计算区域的温度场进行模拟预测。其中，Ma为马格尼数，Re为雷诺数，

Pr为普朗特数，

Θ表示无量纲过余温度，

熔区内部流函数等值线随时间的变化如图6所示，从图中可以看出，起始时刻涡心位于两相界面附近，从t＝7.0到t＝10.0时刻，流函数等值线逐渐向自由面靠近，这也反映了此时间段内涡心的运动趋势。在表面张力梯度的作用下，熔区液桥的表面发生变形，熔区内部产生回流，形成一对胞元。同时，周围空气中由于界面上粘性力的作用也产生一对胞元流。在微重力环境中，由于浮力对流几乎消失，浮区法晶体生长中的主要对流成为由不平衡表面张力引起的表面张力流。表面张力对流与晶体生长的质量有着密切的关系。振荡表面张力流或者紊流是微重力环境下浮区法晶体生长过程中条纹的来源之一。

步骤2.3：采用Level Set Method方法对两相界面位置进行计算；

分析依赖于时间、位置、界面几何性质以及外部物理特性的速度场作用下的界面并发运动；

步骤2.3.1：将当前时刻熔体的运动界面构造成Level Set函数的零等值面，对计算区域的任一点到零等值面的距离进行计算，若距离小于零，则该点位于熔体内部，若距离大于零，则该点位于熔体外部，若距离等于零，则该点为两相界面上的点；

步骤2.3.1.1：将熔体的运动界面看作函数

的零等值面：在每一个时刻t，构造函数

使得在任意时刻，运动界面Γ(t)恰是

的零等值面，即：

$\mathrm{\Γ}\left(t\right)=\{\stackrel{\‾}{x}\∈\mathrm{\Ω}:\mathrm{\φ}(\stackrel{\‾}{x},t)=0\}---\left(17\right)$

步骤2.3.1.2：值满足在Γ(0)附近为法向单调，在Γ(t)上为零，一般取

为

点到界面Γ(0)的符号距离，即

$\mathrm{\φ}(\stackrel{\‾}{x},0)=\left\{\begin{array}{cc}d(\stackrel{\‾}{x},\mathrm{\Γ}\left(0\right))& \stackrel{\‾}{x}\∈{\mathrm{\Ω}}_1\\ 0& \stackrel{\‾}{x}\∈\mathrm{\Γ}\left(0\right)\\ -d(\stackrel{\‾}{x},\mathrm{\Γ}\left(0\right))& \stackrel{\‾}{x}\∈{\mathrm{\Ω}}_2\end{array}\right.---\left(18\right)$

Ω₁、Γ(0)、Ω₂分别表示两种介质区域及其边界(如图7所示)，

表示点

点到界面Γ(0)的距离。

步骤2.3.1.3：在任意时刻，函数φ的零等值面就是运动界面，从而有：

$\frac{\mathrm{d\φ}}{\mathrm{dt}}=\frac{\∂\mathrm{\φ}}{\∂t}+\stackrel{\‾}{V}\·\▿\mathrm{\φ}=0---\left(19\right)$

对于具体问题方程(19)有具体的形式，在自由追踪或二相流问题中，主场物理量控制方程一般是某种具体形式的N-S方程，此时方程(19)就是：

φ_t+uφ_x+vφ_y＝0                (20)

其中，

是流体速度。

步骤2.3.2：对两相界面位置进行修正：根据下一时刻的熔体面积与起始时刻的熔体面积之差及界面的曲率关系确定修正比例，对两相界面位置按修正比例进行修正；

采用重新初始化手段改造

使其重新成为

到界面Γ(t)的符号距离。

步骤2.3.2.1：设在时刻t求得Level Set函数值φ₀，重新构造函数

使之满足下面两个条件：

(a)、

满足(17)，即是符号距离函数；

(b)、

与φ₀有相同的零等值面，即满足式(18)。

步骤2.3.2.2：在满足以上两个条件基础上，通过求解初值问题的稳态解来实现

为满足(19)的符号距离函数；

初值问题如下：

$\mathrm{sig}{n}_r\left({\mathrm{\φ}}_0\right)=\frac{{\mathrm{\φ}}_0}{\sqrt{{{\mathrm{\φ}}_0}^2+{\mathrm{\ϵ}}^2}}---\left(21\right)$

步骤2.3.2.3：跟踪面积损失以重新修正Level Set函数值，采用一个变化的比例来与其损失的面积进行适度的比例修正，具体修正方法是求解下式至稳定状态：

$\frac{\∂\mathrm{\φ}}{\∂\mathrm{\τ}}+(A_0-A\left(t\right))f\left(\mathrm{\κ}\left(\mathrm{\φ}\right)\right)|\▿\mathrm{\φ}|=0---\left(22\right)$

其中，A₀为初始时刻(t)＝0的总质量，A(t)为时刻t对应于

的总质量，f＝(κ(φ))为关于曲率的函数，可以有不同的表达形式；

步骤2.3.3：返回步骤2.3.1，对下一时刻两相界面位置进行模拟计算；

步骤3：将计算区域的速度场、温度场、流场和压力场的计算结果和两相界面位置计算结果分类且按时间顺序存储；

步骤4：将计算区域的速度场、温度场、流场、压力场数据以及两相界面位置，以曲线图、矢量图或云图的形式分类显示。

步骤2.1中划分网格后形成的计算区域施加匀强磁场，使计算区域的速度场、温度场、流场和压力场发生变化，进而使热毛细对流发生改变。

采用匀强磁场控制熔区内高Pr数流体热毛细对流。均匀磁场包括轴向均匀磁场和横向均匀磁场，其示意图如图8所示，轴向均匀磁场方向如图8(a)所示，磁场方向为Z轴正方向；横向磁场磁力线与生长轴垂直，横向均匀磁场方向如图8(b)所示，磁场方向为X轴正方向。重复步骤1的操作，但在步骤1-5中操作人员需要加入额外的轴向磁场强度B分别为0T、0.1T、0.2T、0.3T、0.4T、0.5T，对应的Ha数分别为0、16、32、48、64、80(Ha为哈德曼数，是液态金属磁流体流动中反应磁场大小影响的磁流体力学相似准则数。反映导电流体在磁场中运动时所受的洛伦兹力和粘性力的相对大小)。

重复步骤2.2～步骤2.3，操作人员根据加入匀强磁场后的模拟预测结果，绘制出流函数等值线随轴向磁场强度的变化图以及不同轴向磁场强度下液桥表面变形图。流函数等值线随轴向磁场强度的变化如图9所示，不同轴向磁场强度下液桥表面变形图如图10所示。

图9图中所取的时刻为无量纲时刻t＝500，这个时刻液桥中的流动基本已经达到平衡状态，液桥中流动机构及温度的变化随着时间的增加变化极小。随着轴向磁场强度的增加流函数等值线图不断向自由表面聚集，同时自由表面的涡越来越扁，这是由于轴向磁场的洛伦兹力的方向与径向平行，磁场对于熔区内的热毛细对流有抑制作用，能够抑制流体的径向流动，而且随着磁场强度的不断增加其抑制作用不断增强。而对轴向流动的影响较小。由图8可以发现随着磁场强度的增加在液桥内部又出现了一个涡，这是因为轴向均匀磁场的方向平行于液桥的轴向，对与磁场方向垂直的径向的流体流动具有明显的抑制作用，而对轴向的流体流动的抑制作用不明显，从而自由表面附近的流体的径向流动并不能有效渗透到熔区中心轴附近区域就在熔区近自由表面附近区域回流，从而形成了自由表面附近的涡流。

图10说明在同一时刻随着轴向磁场强度的增加液桥的表面变形幅度也在增加，在靠近热端附近区域内液桥的径向变形随着磁场强度的增加而逐渐向外即为正值，在靠近冷端附近区域内液桥的径向变形向内即为负值且随着磁场强度的增加变形幅度也不断的增大，图中的变形形状与实验相符。从中可以看到熔区液桥表面在微重力环境中因其热毛细对流的产生确实存在着变形。

操作人员根据需要，绘制出计算区域内流场随磁场强度的变化图，如图11所示。不同磁场强度下液桥表面变形图，如图12所示。

图11是在t＝300时刻，不同横向磁场强度下液桥及周围空气的流场速度矢量图。从图中可以发现同一时刻随着磁场强度的增强对于液桥内部的影响不像轴向磁场那样会在熔区内部形成一个近似静止的区域。在液桥整个内部仍然存在着对流，只是对流的速度随着磁场强度的增强而逐渐减小。熔区中热毛细对流的中心也会随着磁场强度的增强而逐渐向液桥表面靠近。

图12显示了在t＝450时刻熔区界面变形随横向磁场强度的变化，可以看到在同一时刻随着磁场强度的增加液桥的表面变形幅度逐渐减小，这是因为在浮区法的晶体生长中，横向静态磁场垂直于晶体生长轴向，可以较好地抑制熔体表面张力对流。随着横向磁场强度的增强晶体生长界面形状趋于平坦，结果与实验结果吻合较好。

综上所述，通过平台的预测分析，可以得出以下结论：随着时间的发展，熔区液桥中的热毛细对流会逐渐进入稳定状态，热毛细对流的中心位置从液桥的热端下移到液桥的中部，并保持稳定，液桥中的热毛细对流主要集中在液桥表面附近，在液桥内部的流动基本是静止的。均匀轴向磁场对液桥内热毛细对流具有抑制作用，这种抑制作用随着磁场强度的增强不断增强，但轴向均匀磁场对自由表面附近的强烈对流抑制作用不明显，形成了自由表面局部对流强烈的区域和液桥内部对流较弱的近似静止区域。在磁场强度0.5T下熔区热端自由面的变形向外，且随着时间的演变其变形程度稍微增大，到稳定状态后维持不变。在液桥冷端自由面的变形向内，随着时间的发展其变形程度逐渐减小，并且变形最大值所处的位置逐渐向冷端移动。总体上，热端自由表面变形幅度小于冷端自由表面变形幅度。

加入横向磁场后，会对熔区液桥内部的热毛细对流产生一定的抑制作用，使整个熔区液桥内部热毛细对流强度减弱；由于横向磁场对热毛细对流中的轴向速度有很好的抑制作用，液桥内部的温度分布因此就会更均匀，呈现热传递型温度分布，受热毛细对流的影响比较小。同时，随着磁场强度的增强液桥自由表面变形幅度减小，在磁场强度为0.5T时液桥自由面变形幅度约为无磁场时的2倍，表明随着横向磁场强度的增大，实际晶体生长界面形状将趋于平坦。

本实施方式的系统及方法克服了以往绝大多数预测手段没有同时考虑周围气相及熔区自由面变形，而是将熔区自由面作为预测域的外边界或者以Young-Laplace方程定义固定表面曲率等造成的影响。

Claims (4)

1.一种熔体热毛细对流过程模拟预测系统，其特征在于，包括：

熔体数据采集模块：用于采集熔体液桥区域的几何参数、介质的物性参数、熔体所在环境的物性参数和熔体热毛细对流过程初始状态参数，并传输至热毛细对流模拟预测模块；

所述熔体液桥区域的几何参数包括：液桥直径、液桥高度、环境测宽度、环境测高度；

所述熔体液桥区域的介质物性参数包括：熔区介质的密度、熔区介质动力粘度、熔区介质表面张力随温度变化系数、熔区介质比热容、熔区介质导热系数和熔区介质热膨胀系数；

所述熔体所在环境的物性参数包括：空气的密度、空气动力粘度、空气的比热容、空气导热系数和环境温度；

所述熔体热毛细对流过程初始状态参数包括：上盘温度、下盘温度和重力加速度；

热毛细对流模拟预测模块：包括网格划分子模块、热毛细对流计算子模块和界面识别子模块；

其中，网格划分子模块用于根据熔体数据采集模块采集的熔体液桥区域的几何参数，确定熔体液桥区域并划分成网格，形成计算区域；

热毛细对流计算子模块：用于对计算区域的速度场、温度场、压力场和流场进行模拟预测，并将结果分别传输至界面识别子模块和模拟预测结果存储模块；

界面识别子模块：用于根据计算区域的速度场、温度场、流场和压力场对两相界面位置进行计算，并将计算结果传输至模拟预测结果存储模块；

模拟预测结果存储模块：用于将计算区域的速度场、温度场、流场和压力场的计算结果和两相界面位置计算结果分类且按时间顺序存储，并传输至结果图像显示模块；

结果图像显示模块：用于将模拟预测结果存储模块传来的数据以曲线图、矢量图或云图的形式分类显示。

2.根据权利要求1所述的熔体热毛细对流过程模拟预测系统，其特征在于：所述热毛细对流模拟预测模块还包括匀强磁场施加子模块，用于对网格划分子模块划分网格后形成的计算区域施加匀强磁场，使计算区域的速度场、温度场、流场和压力场发生变化，进而使熔区热毛细对流发生改变。

3.采用权利要求1所述的熔体热毛细对流过程模拟预测系统进行熔体热毛细对流过程模拟预测的方法，其特征在于：包括以下步骤：

步骤1：采集熔体液桥区域的几何参数、介质物性参数、熔体所在环境侧流体的物性参数和熔体热毛细对流过程初始参数；

所述熔体液桥区域的几何参数包括：液桥直径、液桥高；

所述熔体液桥区域的介质物性参数包括：熔区介质的密度、熔区介质动力粘度、熔区介质表面张力随温度变化系数、熔区介质比热容、熔区介质导热系数和熔区介质热膨胀系数；

所述熔体所在环境的物性参数包括：空气的密度、空气动力粘度、空气的比热容、空气导热系数和环境温度；

所述熔体热毛细对流过程初始状态参数包括：上盘温度、下盘温度和重力加速度；

步骤2：根据熔体液桥区域的几何参数进行热毛细对流模拟预测；

步骤2.1：根据熔体数据采集模块采集的熔体液桥区域的几何参数，确定熔体液桥区域并划分成网格，形成计算区域；

步骤2.2：对计算区域的速度场、温度场、流场和压力场进行模拟预测；

步骤2.3：采用Level Set Method方法对两相界面位置进行计算；

步骤2.3.1：将当前时刻熔体的运动界面构造成Level Set函数的零等值面，对计算区域的任一点到零等值面的距离进行计算，若距离小于零，则该点位于熔体内部，若距离大于零，则该点位于熔体外部，若距离等于零，则该点为两相界面上的点；

步骤2.3.2：对两相界面位置进行修正：根据下一时刻的熔体面积与起始时刻的熔体面积之差及界面的曲率关系确定修正比例，对两相界面位置按修正比例进行修正；

步骤2.3.3：返回步骤2.3.1，对下一时刻两相界面位置进行模拟计算；

步骤3：将计算区域的速度场、温度场、流场和压力场的计算结果和两相界面位置计算结果分类且按时间顺序存储；

步骤4：将计算区域的速度场、温度场、流场、压力场数据以及两相界面位置，以曲线图、矢量图或云图的形式分类显示。

4.根据权利要求3所述的熔体热毛细对流过程模拟预测的方法，其特征在于：所述步骤2.1中划分网格后形成的计算区域施加匀强磁场，使计算区域的速度场、温度场、流场和压力场发生变化，进而使热毛细对流发生改变。

Images