镱金属真空蒸发装置坩埚及束形结构如图1所示，由于镱金属本身蒸发温度较低，加热方式采用欧姆电阻加热。为控制蒸发束流发散角，坩埚出口处安装有蜂窝板、多层狭缝板等束流结构。

在该装置中，由坩埚内部物料蒸发产生的原子束流在蜂窝孔道结构中多次碰撞后束流形状得到初步约束，束流密度下降，之后发散角较大的束流沉积在狭缝板表面，仅少部分发散角较小的束流能够通过多层狭缝板结构，参与后续光电离等步骤。

本文所用的DSMC方法是稀薄气体动力学研究中常用的数值模拟方法之一，该方法用少量虚拟粒子代替真实粒子，实时追踪每个虚拟粒子信息，引入碰撞频率的概念将粒子间碰撞过程转变为随机过程，实现了粒子运动与碰撞过程的解耦，极大减小了计算压力。该方法广泛应用于稀薄气体动力学研究，可靠性已被多位学者证明^[14-17]。一般情况下，DSMC方法的计算流程一般可以分为以下八个运算步骤：

（1）给定温度、气体密度、粒子与物面作用条件等物理参量和虚拟粒子数量、网格尺寸、时间步长等模拟参量；

（2）给定系统初始状态下虚拟粒子信息，确定粒子运动的边界条件；

（3）计算单位时间步长内虚拟粒子的运动信息；

（4）判断虚拟粒子是否与边界发生作用，如果粒子与边界发生作用，优先计算作用后粒子的状态信息；

（5）统计仍留存在模拟区域内的粒子数量以及各网格内的粒子数量并进行排序；

（6）计算各网格内的碰撞过程，更新粒子碰撞后的速度信息；

（7）统计网格内虚拟粒子参数并计算所需宏观量（密度、温度等），判断其是否达到稳态，如未达到，重复步骤3-7；

（8）统计所需参量，判断参量在空间或时间上的数值波动是否达到要求，如没有，重复步骤3-8。

上述运算步骤的重点在于代表粒子与物面作用的步骤4和代表粒子间碰撞过程的步骤6。下面将对步骤4和步骤6中涉及的理论模型做简单介绍。

粒子与物面作用形式包括吸附和反射两种。吸附过程的处理相对简单，粒子沉积在表面后便不再参与后续束流传输过程。反射过程相对复杂，目前普遍应用的碰撞模型有完全漫反射模型、Maxwell类型边界条件、Cercignani、Lampis和Lord共同完善的CLL模型等。本文粒子与物面作用模型选取Maxwell类型边界条件，即认为来流的$\sigma (0 < \sigma < 1)$部分发生漫反射，而$(1 - \sigma )$部分发生镜面发射，相较于完全漫反射模型更符合情况，相较于CLL模型更容易实现。一般情况下认为$\sigma $的取值与来流速度、与法向的夹角、物面温度等参数无关^[14]。对于镱等低温下饱和蒸气压较高的金属，接近坩埚出口的壁面与坩埚温度差距较小，表面吸附的金属可能发生二次蒸发。考虑到蒸发后粒子速度分布与粒子发生漫反射时速度分布接近，可以将二次蒸发的粒子等效为发生漫反射的粒子，简化分析过程。

粒子间碰撞模型种类较多，有碰撞截面固定的硬球模型、碰撞截面与相对速度成比例关系的变径硬球（variable hard sphere, VHS）模型^[18]、适用于存在不同种分子间碰撞的变径软球（variable soft sphere, VSS）模型碰撞^[19]等。本文涉及的蒸发实验在高真空情况下进行，没有考虑背景气体的影响，模拟程序中碰撞截面$ \sigma\mathrm{_T} $的计算采用VHS模型，公式如下：

式中，${\sigma _{{\text{r,ref}}}} = {\text{π}} d_{{\mathrm{ref}}}^2$为参考温度${T_{{\text{ref}}}}$下的参考碰撞截面，${c_{{\text{r,ref}}}} = 4\sqrt {k{T_{{\text{ref}}}}/{\text{π}} m} $，为参考相对速度，$ c\mathrm{_r} $为粒子间相对速度。模拟程序中${T_{{\text{ref}}}} = 2000{\text{ K}}$，${d_{{\text{ref}}}} = 0.590{\text{ nm}}$，$\eta = 6.669$。

为确保本文所用程序的可靠性，选取金属钆二维平面蒸发过程进行模拟验证。如图2所示，对于坩埚长度L远大于蒸发液面宽度D的长槽型坩埚蒸发过程，可以将其视为无限长窄平面蒸发。由于该蒸发过程左右对称，模拟时只需要模拟1/2蒸发区域。

如图3所示，模拟区域尺寸为$ 30\ \mathrm{cm}\times40\text{ cm} $，上边界和右边界取自由边界条件，左边界取对称边界条件，下壁面取完全漫反射条件，壁面温度${T_{\text{w}}} = 3414{\text{ K}}$，蒸发源半宽度$R = 5{\text{ cm}}$。蒸发源表面温度$ T\mathrm{_{\mathrm{s}}}=3414\text{ K} $，Knudsen数$Kn = 0.04$，粒子碰撞模型选取VHS模型。

网格划分方面，根据前人模拟经验，网格尺度应取网格内粒子平均自由程的1/3左右，时间步长应保证单位时间步内粒子运动距离小于网格尺度的同时远小于粒子平均碰撞时间，网格内粒子数保持在30左右^[20]。由于金属蒸发过程粒子数密度变化较为剧烈，本文所用模拟程序中网格划分以蒸发源处碰撞最剧烈的区域为基础，保证接近蒸发源处网格尺度在该处粒子平均自由程的1/3~1/2，远离蒸发源处网格尺度适当增大；时间步长以确保蒸发源处粒子单个时间步运动距离小于网格尺度为准；单个网格粒子数除蒸发源表面和距蒸发源过远的区域外基本保持在20~80水平，后续模拟过程网格划分结果类似。

模拟结果和用DS2 V模拟程序获得的结果^[21]对比如图4~图7所示，所有参数均已归一化处理。其中，原子数密度的归一化因子是蒸发液面表面的原子数密度，速度参量的归一化因子是蒸发温度下的最可几速度$ v_{\mathrm{m}}=\sqrt{2kT_{\mathrm{s}}/mR} $，温度的归一化因子是蒸发液面的温度。

可以看到，除蒸发源出口处部分网格点由于网格尺寸结果有较大差异外，整体上本文所用程序的模拟结果与DS2 V程序模拟结果比较接近，远离蒸发源处原子速度分布的相对误差不超过13%，温度分布的相对误差不超过5%，且模拟结果比较符合二维平面蒸发过程的理论描述，在一定程度上验证了本文所用程序的可靠性。

对于图1所用蒸发装置，针对束形结构尺寸参数、入口束流参数、壁面作用参数等影响最终出射束流速度和密度分布的主要参数，分别开展束流在蜂窝孔道和狭缝板间传输过程的数值模拟，结果如下。

蜂窝板由若干尺寸相同的孔道组成，孔径略小于实际蒸发温度（700-820 K）下镱原子的平均自由程，且各孔道入口条件近似相同，因此束流在蜂窝板中的传输过程可以近似为在单个孔道中的传输过程，采用一维DSMC方法进行数值模拟。

如前文所言，影响出口束流特性的主要参数包括束形结构尺寸参数、入口束流参数、壁面作用参数。对于束流在蜂窝板中的传输过程，束形结构尺寸参数为孔道半径与长度，入口束流速度分布认为服从蒸发温度下的Maxwell分布，角分布服从余弦分布，束流密度分布近似为均匀分布，可变参量为密度大小。考虑到实际应用中蜂窝板温度与坩埚温度接近，可认为蜂窝孔道表面不发生金属沉积，壁面条件选取Maxwell边界条件，可变参量为漫反射比率σ（孔道长度较短，可近似认为孔壁物性条件均匀，σ取值一致）。

模拟区域分别选取$ \mathit{\Phi}2\ \mathrm{mm}\times\mathit{\Phi}10\text{ mm} $、$ \mathit{\Phi}2\ \mathrm{mm} \times \mathit{\Phi}20\text{ mm} $、$ \mathit{\Phi}2\ \mathrm{mm}\times\mathit{\Phi}30\text{ mm} $、$ \mathit{\Phi}2\ \mathrm{mm}\times\mathit{\Phi}40\text{ mm} $的圆柱形孔道；蒸发源侧粒子均匀分布，出口处取自由边界条件，从两侧离开模拟区域的粒子不再返回模拟区域；模拟过程不考虑粒子间相互作用，壁面条件取$\sigma = 1$的情况。模拟结果如图8所示。

图8（a）中，数据点为DSMC方法计算值，实线为Clausing等针对无粒子间碰撞、壁面条件为漫反射条件下束流在长圆柱形孔道中运动过程得到的解析解^[22]，形式如下：

式中，$\beta = 4R/3L$，R和L分别为孔道半径和长度。不同长径比条件下数值解与解析解间相关系数分别为0.988（$L/R = 10$）、0.986（$L/R = 20$）、0.996（$L/R = 40$）。

图8（a）表明，孔道长径比越大，出口束流发散角越小。进一步地，利用该解析解考虑孔道对束流约束作用的极限情况，即当孔道长径比无穷大时，式（2）中$ p = \dfrac{L}{{2R}}\tan \theta < 1 $的情况占比接近0，此时出口束流通量角分布$I(\theta ) = 2\pi \sin \theta \cdot T(\theta )$（未进行归一化处理）可近似表示为：

式（3）即为长圆柱形孔道束流发散角最佳分布形式，此时出口束流平均发散角25°，最可几发散角25.23°，$\theta < {34^ \circ }$的范围内集中了约50%的流量，$\theta < {56^ \circ }$的范围内集中了约80%的流量，相较于入口束流速度角分布$I(\theta ) = \sin 2\theta $的情况（平均发散角、最可几发散角均为45°，$\theta < {45^ \circ }$的范围内集中了约50%的流量，$\theta < {63^ \circ }$的范围内集中了约80%的流量）有一定提升。

图8（b）中，完全漫反射条件下，孔道透过率与孔道长径比近似服从反比关系。对此作如下解释：

如图9所示，假设长度为${L_1}$的孔道透过率为α，长度为${L_2}$的孔道透过率为β。对于完全漫反射壁面条件，粒子碰撞后速度方向随机。考虑孔道足够长的情况，粒子必须与孔壁发生多次碰撞才能离开孔道，则对于长度为${L_1} + {L_2}$的孔道，单个粒子从孔道入口向右运动第一次到达$ {L}_{1} $位置处的概率为α，从左侧离开孔道的概率为$1 - \alpha $，粒子继续运动到达${L_1} + {L_2}$位置离开孔道的概率近似为$\alpha \beta $，向右运动重新回到$ {L}_{1} $位置处的概率近似为$(1 - \beta )\alpha $。如此往复，可得长度为${L_1} + {L_2}$的孔道透过率为$\alpha \beta /(\alpha + \beta - \alpha \beta )$。

假设长度为L的孔道透过率为α，长度为$kL$的孔道透过率为${a_k}$，则长度为$(k + 1)L$的孔道透过率表达式如下：

对等式左右取倒数求解可得：

当孔道长径比较大时，孔道透过率α与长径比近似为反比关系。

模拟区域选取$ \mathit{\Phi}2\ \mathrm{mm}\times\mathit{\Phi}30\text{ mm} $的圆柱形孔道，改变壁面作用参数σ取值，其余条件与2.1.1小节相同，计算结果如图10所示。

模拟结果表明，Maxwell边界条件下，漫反射占比越高，孔道透过率越低，出口束流发散角越小，当漫反射占比达到1时，二者达到极限，出口束流速度角分布近似服从式（1）。

对该结果做如下解释：漫反射与镜面反射两种反射方式中粒子速度方向变化如图11所示。粒子发生镜面反射垂直于孔道方向速度反向，平行于孔道方向的速度不变，几乎不影响束流发散角与孔道透过率；粒子发生漫反射时离开壁面的粒子速度方向满足余弦分布，大小满足壁面温度下的Maxwell速率分布，在不考虑粒子间碰撞的情况下，原本不能直接离开孔道的粒子发生漫反射后有一定几率以更小的角度直接离开孔道，因此长孔道中漫反射对束形有一定的约束作用，但由于漫反射有概率使粒子沿孔道方向速度反向，因此漫反射占比过高孔道透过率会有所下降。

在金属蒸发过程模拟中，初始束流密度一般通过金属蒸发饱和蒸气压经验公式确定，仅与蒸发温度相关。程序中所用镱饱和蒸气压的经验公式如下所示：

该式由623 K~931 K温度范围内的157个蒸发数据拟合得到的，式中压强单位为mmHg，温度单位为K^[23]。

这里为避免引入温度变量，直接给定液面处的Knudsen数，利用Knudsen数定义（原子平均自由程λ与流动区域特征长度L之比）确定λ，再由λ的表达式（7）求解初始束流密度。

蒸发源表面Kn取0.05、1.00、5.00、25.00，分别对应$\lambda \ll R$、$\lambda \approx R$、$R < \lambda < L$和$\lambda > L$四种经典情况，模拟区域选取$ \mathit{\Phi}2\ \mathrm{mm}\times\mathit{\Phi}10\text{ mm} $的圆柱形孔道，考虑粒子间相互作用，其余条件与2.1.1小节相同，计算结果如图12所示。

该结果表明气体稀薄程度提高，出口束流发散角减小，孔道透过率下降。这是由于气体稀薄程度提高导致粒子间速度交换频率下降，与壁面作用后速度反向的粒子通过粒子间碰撞调整速度方向的概率下降，速度分布逐渐偏离Maxwell分布，同时孔道通过率下降。但当气体稀薄程度足够高时，孔道透过率又略有提升，这一点在其他类似的模拟工作中也有所体现^[24]。

束流在狭缝板间传输过程采用二维DSMC模拟方法，模拟区域如图13所示，为实际蒸发区域的一半。与蜂窝孔道通过漫反射约束来流不同，狭缝板主要通过吸附发散角较大的粒子约束来流，狭缝板表面吸附率为1.0时，出口束流发散角最小，密度最低。但如前文所言，镱金属蒸发温度较低，狭缝板与坩埚温差较小，狭缝板、收集板表面并不能认为是全吸收边界条件。此外，实验过程中为避免狭缝被沉积物堵塞，通常采用多层狭缝板结构，喷口支板、各层狭缝板表面温度不同，壁面作用参数也会有所不同。为简化处理，壁面作用形式考虑反射与吸附两种情况，反射模型选取Maxwell反射条件，根据类似研究经验，将σ取值固定为0.85 ^[25]，仅考虑吸附率的变化。其余边界中侧壁与蒸发源距离较远，温度相对较低，取全吸收边界条件，y=0处取对称边界条件。

双层狭缝板原始结构参数设定如下：模拟区域整体尺寸为$ 200\ \mathrm{mm}\times110\text{ mm} $，蒸发源尺寸$2{L_0}$为4 mm，位于喷口支板下方10 mm处，蒸发温度为810 K，蒸发速率设定为50.0 g/h，源表面金属原子均匀分布，速度大小服从Maxwell分布，方向服从余弦分布；喷口支板处狭缝宽度$2L$为20 mm，表面温度640 K，吸附率0.5，反射中漫反射占比0.85；狭缝板1尺寸$2{L_1}$、狭缝板2尺寸$2{L_2}$均为8 mm，表面温度分别为610 K、580 K，吸附率均取0.50，反射中漫反射占比0.85；收集板表面温度为550 K，吸附率0.50，反射中漫反射占比0.85；喷口支板与狭缝板1间距${d_1}$、狭缝板1与狭缝板2间距${d_2}$分别为20 mm、30 mm，狭缝板2与收集板间距${d_3}$分别为20 mm、30 mm、60 mm；侧壁面采用全吸收边界条件。

改变各级狭缝板尺寸参数，模拟结果如表1、表2所示。

上表中，${\alpha _0} \sim {\alpha _3}$分别代表喷口支板、狭缝板1、狭缝板2、收集板表面沉积物占比，${R_0}$为收集板上的镀斑尺寸，用以表征狭缝板2处出口束流的发散程度，$(0,{R_0})$范围内包含了约80%的沉积量。

模拟结果表明，提高单级狭缝板尺寸，该级狭缝板和上一级表面沉积量下降，下一级狭缝板和收集板表面沉积量上升，镀斑尺寸增大，表征出口束流发散角增大，该级狭缝板尺寸越小，该现象越显著，缩小狭缝板间距效果相同。此外，改变狭缝板2相关参数相较于改变狭缝板1对出口束流影响更显著，在蒸发装置可接受调整范围内改变狭缝尺寸相较于改变板间距对出口束流影响更显著。可以推断，对于多层狭缝板结构，出口束流受最后一级狭缝的约束最为显著。狭缝调节应优先考虑最后一级狭缝的大小和位置，前级狭缝板应当通过合理设计分担多余的金属沉积，避免实际蒸发过程中最后一级狭缝发生堵塞。

受表面温度分布、金属沉积量分布等多种因素的影响，实际情况下狭缝板表面吸附率的取值相当复杂，这里仅考虑单级狭缝板表面吸附率取值相同且不随时间变化的简单情况，对吸附率与出口束流特性间的关系做定性讨论，具体模拟条件除各级吸附率有所改变外与2.2.1小节模拟条件相同，模拟结果如表3所示。

表中${\gamma _0}$、${\gamma _1}$、${\gamma _2}$、${\gamma _3}$口支板、狭缝板1、狭缝板2、收集板表面吸附率。表3中，提高单级狭缝板吸附率，该级狭缝板表面金属沉积量明显上升，相邻狭缝板和收集板表面金属沉积量下降，同时出口束流发散角减小，且各级吸附率对出口束流发散角的影响从大到小排序依次为狭缝板2、狭缝板1、喷口支板；统一提高所有狭缝板吸附率，狭缝板1和收集板表面金属沉积量明显上升，喷口支板表面金属沉积量下降，出口束流发散角减小。该结论表明，对于壁面吸附现象明显的真空蒸发装置，出口束流中有相当一部分是经过多次与装置壁面的反射作用后离开束流装置的，且各级狭缝板吸附率越低，经过多次反射离开束流装置的比例越高，这部分束流的存在对出口束流密度、发散角都有较大影响。此外，注意到壁面吸附率中对收集板表面镀斑尺寸影响最大的是收集板吸附率，后续镱金属蒸发装置设计过程中应尽量避免收集板表面粒子反射对出口束流特性的干扰。

实际蒸发过程中坩埚蒸发速率的改变是通过提高加热温度、扩大喷口尺寸等方式实现的，这里为避免引入其他变量，直接改变了坩埚蒸发量，未对蒸发温度、喷口尺寸等参量进行调整，其余参数与2.2.1小节相同，模拟结果如表4所示。

表4中m喷口的处质量流量，代表不同的入口束流密度，喷口尺寸恒为$ 4\ \mathrm{mm}\times380\text{ mm} $，蒸发源温度恒为810 K。表4中，随入口束流密度的上升，喷口支板和狭缝板1沉积量下降，狭缝板2、收集板沉积量上升，出口束流发散角变大。该现象原因与2.1.3小节类似，粒子间碰撞破坏了束流运动的单向性。该结果也说明蒸发量大于25 g/h束流在狭缝板间传输过程模拟时粒子间的碰撞作用不可忽略。此外，注意到蒸发量小于100 g/h缝板表面沉积量变化并不显著，说明此时粒子间碰撞主要集中在狭缝中心区域。需要注意的是，随着入口束流密度的提高，出口束流发散角变大，同时狭缝板表面金属沉积速率提高，长时间蒸发使得狭缝发生堵塞的时间点提前，因此应当尽量避免通过提高入口束流密度的方式提高出口束流密度。

上述模拟结果表明应尽量避免通过提高入口束流密度的方式提高出口束流密度，因此，蒸发装置束形器设计时应考虑在出口束流发散角满足需求时束形器透过率保持在较高水平。考虑到蜂窝板结构对束流发散角的优化作用有限，并对束形器透过率有一定影响（长径比L/R=20的孔道出口束流平均发散角约37°，孔道透过率约11%），且蜂窝板内壁温度与蒸发温度相差较大时易发生堵塞。因此，束形器设计过程中仅保留坩埚出口处的蜂窝板。

多级狭缝板结构中，影响束流特性的主要因素包括狭缝板结构参数、壁面吸附率、入口束流密度等。其中狭缝板结构参数与入口束流密度属于可测量参量，而壁面吸附率取值由狭缝板温度分布、金属沉积速率等多种因素决定。壁面吸附率与上述因素间具体数值关系尚不确定，且实际蒸发过程中无法测量该参量，严重阻碍了束形器参数优化工作。

为解决该问题，狭缝板结构优化过程中，通过调整各级狭缝板壁面吸附率取值，确保结构、蒸发量、蒸发温度等参数比较接近的模拟结果与实验结果相对吻合，再以该壁面吸附率取值计算出口束流特性满足需求时狭缝板结构、蒸发量、蒸发温度参数的调整范围，为实际参数调整提供参考。当各参数变更较大时，需重复上述过程重新确定各级狭缝板壁面参数。

基于DSMC方法对蜂窝孔道和狭缝板两种常用的束形结构进行了初步分析讨论，主要结论如下：

（1）柱形蜂窝孔道长径比越大，粒子与壁面作用时漫反射占比越高，入口束流密度越低，出口束流发散角越小，孔道透过率越低，且孔道出口束流平均发散角最小为22.5°，但当入口束流密度足够低（$Kn > 10$）时，孔道透过率又略有提升；

（2）多级狭缝板结构中最后一级狭缝板对束流约束作用最为显著，减小狭缝尺寸，提高板间距，增大吸附率，减少蒸发量等操作均能使出口束流发散角和狭缝板整体透过率下降。

基于上述分析结果，汇总了部分束形器优化过程的部分经验方法，如多级狭缝板结构中除最后一级外其余狭缝板主要作用为分担多余的金属沉积；相较于改变狭缝板间距改变狭缝尺寸对束流特性影响更显著；应尽量避免通过提高入口束流密度的方式提高出口束流密度等。

最后，总结了束形器优化方法，即通过实验结果确定壁面吸附率，再由得到的壁面吸附率确定束形器参数调整范围。

后续工作中，将继续完善模拟程序，进一步优化壁面吸附率确定方法，确保模拟工作与实际蒸发过程的统一。