用于工艺优化的W–C–B涂层的热壁化学气相沉积（CVD）多物理场建模 安德烈·V·波利根科（Andrey V. Poligenko）、 叶夫根尼·A·鲁班（Evgeny A. Ruban）、 基里尔·M·奥西波夫（Kirill M. Osipov）、 安德烈·A·沙波伦科夫（Andrey A. Shaporenkov） 以及弗拉基米尔·V·杜希克（Vladimir V. Dushik）

《Ceramics》：Multiphysics Modeling of Hot-Wall CVD Deposition of W–C–B Coatings for Process Optimization Andrey V. Poligenko, Evgeny A. Ruban, Kirill M. Osipov, Andrey A. Shaporenkov and Vladimir V. Dushik

　　摘要
本研究开发了一个多物理场有限元模型，用于模拟在热壁管式化学气相沉积（CVD）反应器中，通过六氟化钨（WF6）、氢气（H2）和三甲胺硼烷（(CH3)3N:BH3）气体混合物，在550°C和5托压力下沉积W-C-B涂层的过程。本研究旨在深入理解反应物传输机制，并优化工艺参数，以获得目标钨碳化物或硼化物相。模拟是在COMSOL Multiphysics（版本6.1）中进行的，使用二维轴对称公式耦合了层流、热传递和多组分扩散，并考虑了反应器壁上的异相化学反应。模拟得到的反应物浓度空间分布表明，沿着反应器长度前驱体逐渐耗尽。将13种操作条件下的计算得到的B/W和C/W化学计量比与实验数据进行了比较，证实了在三甲胺硼烷流量较低时形成W和W-B相，在流量较高时则形成基于钨碳化物的涂层。此外，还通过参数扫描确定了合成钨硼化物的最佳参数范围。

1. 引言
现代数字建模方法可以大幅减少实验工作中寻找最佳参数的实证性和迭代性工作。有限元方法（FEM）就是其中之一。FEM将研究系统离散为在节点点相连的简单几何元素，然后使用插值函数表示每个元素内的场变量，以求解控制全局方程。FEM能够高效模拟非线性系统，包括考虑由温度或压力引起的材料属性变化，以及规定不同相之间物质的相互作用[1,2,3]。对于化学气相沉积（CVD）过程，FEM可以在单一耦合模型中计算热和质量传递、流体动力学和化学反应动力学[4,5,6]。这使得可以预测基底表面和反应器体积内的温度分布、反应物和产物的浓度分布以及局部膜生长速率——所有这些对于确保涂层均匀性都至关重要。CVD反应器或特定沉积过程的计算模型经常用于优化任务，其中数字孪生模型使用实验运行获得的实际参数进行校准。由此产生的数字副本可以系统地探索参数空间，快速识别满足目标合成条件的参数[7,8,9]。
CVD过程能够通过气相和异相边界中的化学转化合成各种材料，如粉末、涂层和单个颗粒[2,10,11]。CVD的一个关键要求是前驱体在操作温度和压力下必须在气相中保持稳定，并提供足够的蒸气压来驱动异相反应。作者之前的实验工作展示了使用六氟化钨（WF6）和三甲胺硼烷（(CH3)3N:BH3）进行三元W-C-B涂层的CVD合成[12]。发现沉积膜的相组成对前驱体进料速率非常敏感。这种依赖性表明反应器内硼和碳源（三甲胺硼烷—TMAB）的空间分布不均匀，这限制了异相反应动力学并直接决定了所得钨硼化物和碳化物的化学计量比。尽管有这些实证观察结果，但沉积过程的物理化学机制仍然知之甚少。到目前为止，还没有一个统一的多物理场模型通过耦合热和质量传递、流体动力学以及异相反应动力学来定量描述这一系统。
本研究的主要目的是提供对三元W-C-B涂层化学气相沉积（CVD）过程中反应物传输机制和反应器内发生的物理化学过程的更深入理解。此外，这项工作还研究了前驱体流速的优化，以细化选择性合成钨硼化物或钨碳化物相所需的操作条件。此外，CVD建模中的一个重要挑战是经常缺乏对物理界面和边界条件选择的明确解释，这严重限制了已发布模型的可重复性和适应性。因此，本文的第二个目标是提供一个全面的方法论框架。通过详细说明在软件实现之前计算无量纲相似性准则（如雷诺数、克努森数、佩克莱特数和达姆科勒数），本文旨在展示一种物理上严格的方法来简化模型和选择计算模块。这种方法可以很容易地被从事其他复杂CVD研究的研究人员所采用。

2. 计算细节和模型描述
模拟中使用的沉积条件基于实验结果，总结在表1中。表1. 涂层沉积的工艺参数[12]。CVD反应器模型是使用COMSOL Multiphysics软件包开发的。图1展示了反应区的三维模型，主要反应器参数总结在表2中。图1. 反应器的三维模型：外部视图（a），样品在反应区内的位置（b）。表2. 模型中使用的主要CVD反应器参数。由于实际反应器由圆柱形管组成，因此使用了二维轴对称公式来开发数学模型，以减少所需的计算资源。
为了确定COMSOL Multiphysics中适当的物理模块集合，首先对气体流动特性和反应动力学进行了初步评估。这一评估是通过计算一系列无量纲数（相似性准则）来进行的。这些准则用于评估流动状态（雷诺数Re）、气体稀疏度（克努森数Kn）、对流传质与扩散传质的比率（佩克莱特数Pe）、自然对流的影响（格拉肖夫数Gr）以及发生反应的速率限制步骤（达姆科勒数Da）（表3）。表3. 调查过程（表1中的操作条件1-13）的计算无量纲相似性准则的值范围。以操作条件13（表1）为例，逐步展示了相似性准则的计算过程。根据定义，雷诺数（Re）表征系统中惯性力与粘性力的比率，公式如下：
(1)
其中ρ是流体密度，单位kg/m3；V是流速，单位m/s；D是系统的特征长度，单位m（反应器直径）；μ是混合物的动态粘度，单位Pa·s。
气体密度使用理想气体状态方程（2）计算：
(2)
其中P是系统压力，单位Pa；Mmix是混合物的摩尔质量平均值（根据组分的摩尔分数计算），单位kg/mol；R是通用气体常数，单位J/(mol·K)；T是过程温度，单位K。
流速V（m/s）是通过将标准条件下的气体混合物体积流量（转换为操作条件）除以反应器的横截面积来计算的，公式如下：
(3)
其中Pstd是标准条件下的压力，单位Pa；Tstd是标准条件下的温度，单位K；P是操作压力，单位Pa；T是操作温度，单位K；Qstd是标准条件下的气体混合物体积流量，单位m3/s；Q是操作条件下的气体混合物体积流量，单位m3/s；D是特征长度（反应器直径），单位m；r2是反应器半径的平方。
由于氢在所有操作条件下都是主要成分，因此混合物的动态粘度μ（4）被近似为1.84 × 10?? Pa·s。这个值代表550°C下纯氢的粘度，使用Sutherland常数（S）对于273 K下的氢进行外推[13,14]。由于较重成分（TMAB, WF6）的贡献只会增加而不是减少μ，因此基于混合物中最不粘稠和最轻的成分（氢）来估计粘度可以被视为一种保守的方法。
(4)
其中μ0是273 K下氢的参考粘度，单位Pa·s；S是氢的Sutherland常数（96.67 K，T0 = 273 K）；T0是对应于选定常数S的参考温度（273 K）；T是操作温度（对应于550°C），单位K。
获得的Re = 9.9的值明显低于湍流的临界阈值（Re > 2300），表明气体流动状态是层流的。
克努森数Kn（5）表征了我们系统中的气体稀疏程度。从实际角度来看，这是评估将介质视为连续体的有效性所必需的，在这种情况下，经典的流体力学原理和纳维-斯托克斯方程适用。
(5)
其中λ是分子的平均自由路径，单位m；D是特征长度（反应器直径），单位m。
分子的平均自由路径λ是使用公式（6）计算的：
(6)
其中KB是玻尔兹曼常数，单位J/K；T是沉积过程温度，单位K；是混合物（WF6/TMAB/H2）的有效平均分子直径，单位m；P是压力，单位Pa。
鉴于Kn = 0.00024，介质应被视为连续体，因为分子的平均自由路径比特征长度（反应器直径D = 0.1 m）小几个数量级。
佩克莱特数Pe通过描述体积气体流动（对流）与分子热运动（扩散）的比率来表征质量传输状态。由于研究系统由三种主要成分组成——WF6、TMAB和H2——因此使用公式（7）-（9）计算了每个成分的佩克莱特数：
(7)
(8)
(9)
其中V是气体流速，单位m/s；D是特征长度（反应器直径），单位m；
//是相互扩散系数，单位m2/s。
由于在所有操作条件下氢的体积分数超过60%，因此计算了H2-H2（自扩散）、WF6-H2和TMAB-H2对之间的相互扩散系数DAB。由于没有WF6和TMAB在氢中的参考扩散系数，因此使用273 K和1 atm下类似系统的数据进行了估算（表4）。例如，WF6-H2的扩散系数是通过基于质量减少的比例缩放SF6-H2的数据得出的[15]。TMAB-H2对的扩散系数是通过类比氢中的正丁烷来估算的，考虑到它们相似的分子直径（约5-6 ?）[15]。
表4. 估算的扩散系数[15]。
为了准确考虑研究过程的温度和压力，对表4中选定的扩散系数应用了缩放方法。缩放因子是使用公式（10）计算的：
(10)
其中T是操作过程温度，单位K；T0是已知扩散系数对应的温度，单位K；P0是已知扩散系数对应的压力，单位Pa；P是操作过程压力，单位Pa。
将指数1.75应用于操作温度与参考扩散温度的比率（ΔT），以便更准确地缩放扩散系数，因为扩散系数表现出非线性的温度依赖性。根据Fuller等人的经典工作[16]以及更近期的研究[17,18]，幂律指数1.75能够充分描述大多数实际气体的扩散系数温度依赖性。
H2-H2、WF6-H2和TMAB-H2对的扩散系数的缩放分别使用公式（11）-（13）进行：
(11)
(12)
(13)
使用公式（7）-（9）为操作条件13计算的佩克莱特数（Pe）表明，在模拟中必须考虑扩散和对流质量传输。
格拉肖夫数（Gr）表征了自然对流对质量传输的影响。平均而言，在操作条件1-13（表1）下，进口气体混合物的温度为60°C，这比反应器温度（550°C）低一个数量级。乍看之下，进气流与稳态反应器环境之间的这种显著温差应该会导致明显的自然对流以及气体入口附近可能形成的冷区。然而，需要注意的是，大约0.05米的气体输送管线也位于热区内。在相对较低的流速（V = 0.07–0.3 m/s）和较低的压力条件下，进气流有足够的时间加热到操作温度，进入反应区时已经预热，不会引起径向温度梯度。为了证实这一点，计算了条件13（表1）参数下气体管线“热”部分的传热单元数（NTU），其中管线长度Lline = 0.05 m，直径Dline = 0.02 m。计算是针对氢（H2）进行的，因为它是混合物中含量最多（在条件13下体积占比超过60%）且热容量最高的成分[19]。较重的成分（WF6、TMAB）具有明显较低的热容量，并且升温更快[19,20]。NTU代表换热器容量与气体流动热容量率的比值，使用公式（14）进行计算：
(14)
其中：h 是传热系数，单位为 W/(m2·K)；A 是气体管道的内表面积，单位为 m2；m 是气体的质量流量，单位为 kg/s；Cp 是气体的比热容，单位为 J/(kg·K)。
传热系数 h 是使用公式（15）计算的。H2 的热导率（θ = 0.253 W/(m·K)）是在 578 K 的温度下测得的，该温度介于气体入口温度（333 K）和气体管道壁温度（823 K）之间。努塞尔特数（Nu）用于表征从壁面传递热量的效率，在具有恒定壁温的圆管中，其值为 3.66[21]。
(15)
计算得到的 NTU（公式（14）的物理意义在于，理论上，通过气体管道“高温”部分的气体混合物有足够的时间与管壁温度达到近24次平衡。因此，在整个实验过程中，格拉肖夫数（Gr）将接近于零。这意味着自然对流的影响可以忽略不计，在建模过程中可以安全地不予考虑。
为了使用达姆科勒数（Da）来确定限速步骤，首先必须分析驱动过程的化学反应。在模拟 W–C–B 化合物沉积过程中，一个主要挑战是准确描述加热基底和反应器壁附近边界层内的异质反应。在我们之前的工作中[12]，我们提出了一种机制，该机制使用全局反应方程描述了碳化钨和硼化物的气相形成过程。根据这一反应方案，WF6 被氢（H2）还原生成固态钨，并产生氢氟酸（HF）气体作为副产品，该气体通过气流不断从反应区排出（反应 16）。同时，TMAB 及其热分解产物作为碳和硼的来源，用于合成碳化钨和硼化物（反应式 17（示意图，非化学计量方程），18，19）。
(16)
(17)
(18)
(19)
反应（16）已经得到了广泛研究，是成熟的基于氟化物的钨化学气相沉积（CVD）过程的基本组成部分[22,23,24,25]。通常，WF6 的还原过程是在动力学控制下进行的，这确保了基底上层的均匀厚度和优异的涂层一致性。然而，所建模过程的主要焦点在于碳化钨和硼化物的形成机制。已知 TMAB 在低至 200 °C 的温度下就开始积极分解为三甲胺和硼烷，到 370 °C 时分解基本完成[26]。因此，之前的 NTU 计算（公式（14）证实 TMAB 在进入反应区之前，在气流中已分解为硼烷（BH3）和三甲胺（(CH3)3N）。随后，三甲胺及其热分解产物作为碳源参与碳化钨相的形成，而硼烷则提供硼用于在基底表面沉积硼化钨。鉴于这两种物质在所研究过程的操作温度下会在基底上迅速分解[26,27]，可以合理假设硼和碳的沉积速率超过了反应物向基底输送的速率。因此，与从 WF6 沉积纯钨不同，碳化钨和硼化物的形成过程很可能受到扩散限制。对于我们的模型来说，这意味着 W–C–B 化合物的沉积可能处于动力学-扩散控制机制下。因此，与基于简化三组分系统（WF6、H2 和 TMAB）的佩克莱数（Pe）计算不同，准确确定限速步骤需要计算系统中所有主要异质反应的达姆科勒数（Da）：WF6 的还原以及钨硼化物（WxBy）和碳化物（WxCy）的形成。
达姆科勒数（Da）定义了表面化学反应速率与通过扩散传递反应物质量的速率之比。WF6 被氢还原的达姆科勒数是使用公式（20）计算的。为了估算反应速率（Rreac），将实验观察到的平均层生长速率（Vg = 1 μm/h）除以相应组分的摩尔体积（Vm）。为了评估扩散速率（Rdiff），使用了已知温度范围内或相关系统的尺度化方法来确定扩散系数（DWF6）——类似于之前描述的程序（10）——并根据公式（12）进行定义。
对于 WF6 还原反应（公式 16），在表 1 中指定的流量条件下（Regime 13），达姆科勒数（Da）由公式（20）给出：
(20)
其中：
Rreac 是表面化学反应速率，单位为 mol/(m2·s)；
Rdiff 是扩散速率（扩散通量），单位为 mol/(m2·s)；
Vg 是层生长速率，单位为 m/s；
Vm 是组分的摩尔体积，单位为 m3/mol；
DWF6 是 WF6 的扩散系数，单位为 m2/s；
C0 是体积浓度，单位为 mol/m3；
D/2 是特征尺寸（对应于从流动中心到壁面的特征扩散长度，即反应器半径），单位为 m。
达姆科勒数（Da）小于 1 表明反应（16）受到动力学限制。
接下来我们计算钨硼化物形成反应的达姆科勒数（Da）。前面提到的 BH3 高反应性源于其缺电子的性质（作为强路易斯酸）及其获取电子密度的强烈倾向。这可以通过二聚化或在活性金属表面（在 CVD 中为高温基底和反应区加热的壁面）上的快速化学吸附来实现。由于 BH3 在高温下几乎无法二聚化[27]，富电子的钨表面则成为理想的路易斯碱。路易斯酸-碱相互作用几乎是无障碍的过程，其中金属表面的电子迅速填充硼的空轨道，形成强的 W–B 化学吸附键。释放出的氢随后参与 WF6 的还原，促进氟从反应区的移除，生成氟化氢（HF）。对于所考虑的过程，这意味着硼从 BH3 吸附到钨表面的活化障_DIRS 接近于零。因此，根据阿伦尼乌斯方程，反应速率仅受分子碰撞频率的限制。因此，可以使用碰撞理论来估算这一速率。BH3 的表面动力学速率常数 k 使用公式（21）计算：
(21)
其中：
γ 是碰撞时的反应粘附概率（保守估计为 0.1 [25]）；
v 是 BH3 分子的平均热速度，单位为 m/s。
这些物质的质量传递率使用质量传递系数（公式 22）估算：
(22)
其中：
DWF6 是 BH3 的扩散系数，根据分子量（估算的分子量比例）进行缩放，单位为 m2/s；
D/2 是反应器的特征尺寸（半径），单位为 m。
达姆科勒数（Da）大于 1 表明硼的沉积过程受到扩散限制。因此，在模型中假设 BH3 的表面浓度为零。关于扩散系数的简化需要特别说明。诚然，将钨六氟化物和硼烷这样不同的分子之间的扩散系数进行缩放可能会引入系统误差，可能低估硼的真实扩散系数。尽管如此，公式（23）表明 BH3 在钨表面的化学吸附速率非常快，即使质量传递系数增加数倍，也不会改变关于限速步骤的定性结论。因此，在缺乏实验参考数据的情况下，根据格雷厄姆定律缩放扩散系数也能为该特定过程的工程建模提供足够可靠的结论。
与 BH3 不同，三甲胺（TMA，(CH3)3N）是一种更稳定的、价态饱和的分子。然而，在 550 °C 时，其稳定性显著降低。在较高温度下，甲基内的 C–H 键会减弱。此外，过渡金属是脱氢反应的有效催化剂[28,29,30]。因此，一旦钨表面催化脱除第一个氢原子，TMA 分子会转变为高度活性的自由基，并在基底上迅速分解为氢原子、氮原子和碳原子[30]。根据钨的温度程序反应光谱（TPRS）数据[30]，烷基胺的分解伴随着在 325–550 K 范围内的 H2 释放（在 350 K 和 450 K 处有明显的峰值）。对于 (CH3)3N，在低表面暴露量下，观察到完全且不可逆的分解，生成 H2 和 N2。因此，在 550 °C（823 K）时，TMA 的异质分解步骤相对于其向反应区的质量传递速率几乎是瞬时的。因此，使用碰撞理论估算 TMA 的达姆科勒数（Da）是合理的。TMA 的表面动力学速率常数 Ks(TMA) 使用公式（24）计算，类似于公式（21）：
(24)
其中：
γ 是碰撞时的反应粘附概率（保守估计为 0.1 [25]）；
v 是 TMA（三甲胺）分子的平均热速度，单位为 m/s。
TMA 的质量传递系数 βTMA 使用公式（25）给出：
(25)
其中：
DMA 的扩散系数基于分子量（估算的分子量比例）进行缩放，单位为 m2/s；
D/2 是反应器的特征尺寸（半径），单位为 m。
TMA 的达姆科勒数（Da）使用公式（26）计算：
(26)
达姆科勒数（Da）大于 1 表明硼的沉积过程受到扩散限制。因此，在模型中假设 BH3 的表面浓度为零。需要注意的是关于扩散系数的简化。虽然将钨六氟化物和硼烷这样不同的分子之间的扩散系数进行缩放可能会引入系统误差，但公式（23）表明 BH3 在钨表面的化学吸附速率非常快，即使质量传递系数增加数倍，也不会改变关于限速步骤的定性结论。因此，在没有实验参考数据的情况下，根据格雷厄姆定律缩放扩散系数可以为该特定过程的工程建模提供足够的可靠性。
与 BH3 不同，三甲胺（TMA，(CH3)3N）是一种更稳定的、价态饱和的分子。然而，在 550 °C 时，其稳定性相比标准条件大幅降低。在较高温度下，甲基内的 C–H 键会减弱。此外，过渡金属是脱氢反应的有效催化剂[28,29,30]。因此，一旦钨表面催化脱除了第一个氢原子，TMA 分子会转变为高度活性的自由基，并在基底上迅速分解为氢原子、氮原子和碳原子[30]。根据钨的温度程序反应光谱（TPRS）数据[30]，烷基胺的分解伴随着在 325–550 K 范围内的 H2 释放（在 350 K 和 450 K 处有明显的峰值）。对于 (CH3)3N，在低表面暴露量下，观察到完全且不可逆的分解，生成 H2 和 N2。因此，在 550 °C（823 K）时，TMA 的异质分解步骤相对于其向反应区的质量传递速率几乎是瞬时的。因此，使用碰撞理论估算 TMA 的达姆科勒数（Da）是合理的。TMA 的表面动力学速率常数 Ks(TMA) 使用公式（24）计算，类似于公式（21）：
(24)
其中：
γ 是碰撞时的反应粘附概率（保守估计为 0.1 [25]；
v 是 TMA（三甲胺）分子的平均热速度，单位为 m/s。
TMA 的质量传递系数 βTMA 使用公式（25）给出：
(25)
其中：
DMA 的扩散系数基于分子量（估算的分子量比例）进行缩放，单位为 m2/s；
D/2 是反应器的特征尺寸（半径），单位为 m。
TMA 的达姆科勒数（Da）使用公式（26）计算：
(26)
达姆科勒数（Da）大于 1 表明碳的沉积和随后的钨碳化物形成同样受到质量传递限制的支配——具体来说，是前体分子从气体主体到沉积表面的扩散速率。计算出的达姆科勒值（Da 和 ）之间的接近性表明硼和碳向表面的传递是同时进行的、按化学计量进行的。这一比例本质上是由前体 TMAB 分子的组成（C:B = 3:1）决定的。因此，在高 TMAB 流量下，这可能导致表面碳的过饱和，从而使钨碳化物抑制其他相的形成。这一效应在实验中得到了验证，操作条件 11–13 下沉积层的相组成发生了变化[12]。
值得注意的是，实验中未检测到硼或钨氮化物，这也反映在计算沉积模型中。如前所述，当 TMA 被吸附到涂有钨的基底上时，钨表面会主动脱氢 TMA 分子，将其转化为表面结合的、高度不稳定的自由基，这些自由基迅速分解为碳原子、氮原子和氢原子[30]。这些中间体的例子包括二甲氨基甲基（CH2N(CH3)2）和二甲氨基（(CH3)2N）自由基[28,30]。类似的现象也观察到其他胺（如甲胺、乙胺和三乙胺）[30,31]在加热的过渡金属表面上发生。在 TMA 分解过程中，氮很可能不会以自由原子形式进入气相，而是以化学吸附物种的形式留在金属表面上（N(ads)）。鉴于钨碳化物和硼化物的形成在热力学上更受青睐（图 2），这些表面结合的氮原子通过朗缪尔-欣谢尔伍德机制重新组合，随后以 N2 分子的形式释放。值得注意的是，沉积样品的XRD图谱显示了β-W相的存在[12]。这一相的特点是具有较宽的均匀性范围，并且已知由小半径的间隙原子（如氮原子）稳定[35]。氮原子的这种间隙掺入可能是导致未观察到离散氮化物相的额外原因。图2显示了碳化钨、氮化物和硼化物以及氮化硼的标准吉布斯自由能形成焓变（?G）的Ellingham图[19,36]。该图是使用Ellingham近似方法构建的，假设?H和?S是温度无关的常数。基于之前为建模过程所有区域计算的相似性标准（表3），在COMSOL Multiphysics环境中做出了关于所需计算模块集的合理决策。鉴于气体流的层流特性，因此采用了层流模块来描述其在反应器内的行为。计算出的努森数（Kn）表明处于连续流状态；因此，选择了浓缩物种传输界面来模拟质量传输机制。尽管初步估计表明在指定流速下出现显著温度梯度的可能性较低，但仍在多物理模型中加入了流体传热界面。这一 inclusion 允许准确应用热交换边界条件，从而防止了反应器边缘出现非物理的过热现象。由于TMAB在进入反应器前会分解为硼烷（BH3）和TMA，因此化学模块专门考虑了控制钨、硼和碳沉积到基底上的异相反应。

如前所述，该模型使用了二维轴对称几何形状进行开发（图3）。沉积过程在热壁上模拟，这代表了整个反应器的纵向截面。这种几何形状非常适合评估物种传输的性质，并捕捉气体混合物向反应器出口移动时扩散限制前驱体通量的可能耗尽情况。计算网格（图3）使用“Extra fine”质量设置生成，共包含1100个元素。为了准确解析反应器壁附近边界层内的传输现象（该层在模型中作为沉积表面），应用了边界层网格细化（图3b）。网格独立性研究显示，将元素数量翻倍对计算出的沉积速率的影响小于1%。这证实了解的网格独立性，并证明了所选元素数量的充分性。图3显示了模型的计算网格（a）以及反应器壁附近的边界层网格细化（b）。由于气体介质是浓缩的（努森数，Kn < 0.001），且多个反应受到扩散限制，因此在层流模块中应用了无滑移边界条件。这一条件假设固体反应器壁处的气体速度为零。对于低速粘性流动，这种假设在物理上是合理的，并确保了边界层内速度剖面和前驱体浓度梯度的准确计算。在反应器入口处，应用了入口速度边界条件来定义进入的气体流速。该参数根据表1提供的数据为每个操作区域（1-13）进行了指定。对于出口边界条件，即离开反应器的气体，操作过程压力被设定为667 Pa（5 Torr）。此外，还启用了“Suppress backflow”功能；该功能修正了出口处的压力场，以防止非物理性的涡流和气体重新循环进入反应器，从而提高了数值求解器的稳定性。参考压力被设定为0 Pa，并且物理流动模型被定义为可压缩流动，这对于在过程的低操作压力下准确考虑气体混合物密度变化至关重要。

在流体传热界面中，对反应器壁施加了823 K（550 °C）的恒定温度边界条件。由于氢载气的高热导率以及预热区的传热单元数（方程（14）的计算结果，进入的气体混合物几乎在进入瞬间就与壁达到了热平衡。由异相表面反应产生的热量被认为与对流热通量相比可以忽略不计。因此，传热模拟确认主要反应区在550 °C的近似等温条件下运行。模型的参考温度被设定为标准室温293 K（20 °C）。在控制质量传输的浓缩物种传输界面中，选择了菲克定律扩散模型。分配给反应物的扩散系数对应于之前为相似性标准计算的值。因此，该界面以质量分数而不是摩尔浓度来求解传输方程。由于所有质量分数之和严格限制为一，这种方法本质上确保了整个计算域内的总质量守恒，并在计算反应区内的气体混合物密度和物种浓度时最小化了数值误差。对于扩散限制的前驱体（TMA和BH3，其中Da > 1），在热反应器壁上应用了质量分数=0的边界条件。这反映了解离反应的高速率以及这些分子在异相相互作用下的完全分解。相反，对于动力学控制的WF6还原反应（Da < 1），使用了通用内向通量边界条件。在这种情况下，质量通量由一阶动力学速率方程控制。由于文献中没有关于所研究过程条件下还原反应动力学参数的直接测量数据，因此通过从实验观察到的Vg = 1 μm/h的涂层生长率反推确定了WF6还原的表面反应速率常数（Ks）。表5总结了为沉积物种定义的关键边界条件。表5列出了计算系统的所有组件及其根据物质状态分配的属性。表6列出了计算模型的各个组件。这个耦合的模块系统使用分离的求解器进行求解。求解过程分为三个顺序计算的步骤：流体流动、质量分数和热传递。研究类型被定义为稳态。

为了验证模型，将计算出的前驱体到达基底表面的摩尔比率与在相同沉积条件下实验获得的涂层相组成进行了比较[12]。表7展示了硼和碳与钨的模拟化学计量比与实验结果的比较分析。在系列2（区域5-10）中，实验获得了最平衡的相组成（两相W + WB混合物），模型预测的B/W比率在0.2-1.0范围内。平均而言，这低于WB的1:1化学计量比。这种明显差异可以归因于表面反应的动力学差异。BH3在钨表面的分解是一个几乎无障碍的过程[30]，硼原子有很大的机会掺入钨晶格。相比之下，WF6的还原通过一系列中间阶段进行（WFx），并受到动力学限制[37,38]。因此，钨在实际涂层中的沉积效率可能低于模型预测的值，使得实验观察到的固相化学计量比向硼富集方向偏移。这解释了为什么即使模型预测的B/W通量比小于1，实验中仍然观察到化学计量的WB的形成。在系列1（区域1-4）中，计算出的B/W比率（0.8-3.6）有利于形成较高的硼化物；实验中确实检测到了W2B5、β-W和纯W相。值得注意的是，模型预测的C/W比率显著高于系列2。这与实验中出现的β-W相直接相关。如前所述，亚稳态的β-W相具有A15型晶体结构，由小间隙杂质原子稳定[35]。我们的模型显示，相对于钨的通量，碳和氮（来自TMAB分解）有过量。这种杂质过量可能促进了β-W相的稳定，从而阻碍了进一步的硼化物形成。在系列3（区域11-13）中，实验中的相组成向碳化钨（WC, W2C）发生了显著转变。模型通过为这些区域计算的极高C/W比率（9.3-21.8）解释了这一转变。尽管模型预测的B/W比率（3.1-7.0）理论上足以生成较高的硼化物，如W2B5甚至WB4，但碳的过量促进了热力学稳定的碳化物的竞争性形成，有效地阻止了硼化物的生成。表面碳物种的高浓度很可能屏蔽了硼吸附所需的活性位点，从而将主导反应路径从硼化转变为碳化。

因此，开发的模型充分捕捉了W-C-B系统中相形成的竞争性。它可以用来预测相区域的边界，并确定实现目标沉积涂层相组成所需的最优前驱体流速。为了进一步支持所提出的处理窗口，计算出的化学计量比与已建立的平衡热力学数据进行了关联。虽然目前文献中缺乏低温区域（550 °C）的完整三元W-C-B相图，但二元子系统提供了一个可靠的比较基准。根据Duschanek和Rogl[39]以及Kurlov和Gusev[40]评估的经典W-B相图，计算出的局部原子比率与观察到的固体相吻合。在低TMAB流速下，计算出的B/W比率对应于纯W、W2B和WB平衡共存的热力学区域。相反，较高前驱体进料时C/W比的急剧增加将局部表面化学计量比推向了碳化钨的稳定域，从而解释了实验中观察到的硼化物抑制和向碳化物主导涂层的相转变。尽管低温CVD是一个非平衡动力学过程，但动力学建模的气相/表面化学计量比与热力学相图之间的这种相关性验证了所建立处理窗口的预测能力。图4使用区域13作为代表性示例，展示了反应体积内前驱体（WF6、BH3、TMA和H2）和反应产物（HF）的摩尔浓度分布（mol/m3）：W、B、C（a）和H2及HF（b）的主要来源。从浓度剖面可以看出，随着气体混合物沿反应器移动（从底部到顶部），反应物耗尽区逐渐扩大。这种耗尽在不同组分中的性质不同，因为前驱体的消耗受到它们之间相互作用以及与加热表面的相互作用的动力学和质量扩散的双重影响。虽然这种现象可能会影响涂层的局部化学计量比（图5），但实际实验中的样品位于反应器的中心区域（0.07至0.09 m），在那里扩散驱动的耗尽效应最小。图5显示了反应区内的化学计量比：B/W（a）、C/W（b）。为了量化这一观察结果，并对整个实验参数空间内的耦合传输模型进行独立验证，我们提取了在三种代表性边界条件下的计算沉积速率沿反应器壁的轴向分布：第1阶段（最小前体流量）、第5阶段（中等）和第13阶段（最大前体流量）。这些分布图显示在图6中。由于前体物质的消耗，沿着整个反应器长度，沉积通量逐渐减小，特别是在第1阶段。然而，模型预测在样品放置区域（0.07–0.09米范围内）内沉积速率保持稳定（如图6中高亮显示）。在这个区域内，所有评估阶段的沉积速率变化仅为3–6%。这种空间一致性与实验观察到的均匀涂层厚度非常吻合，为传输模块提供了可靠的独立验证。图6显示了沿反应器壁计算的轴向沉积速率分布。此外，图4中的氟化氢（HF）浓度分布也值得注意，因为它在接近反应器出口的边界层中积累。高浓度的氟化氢可能会对形成的涂层产生蚀刻作用，并根据勒夏特列原理通过吸附的氟原子毒化活性表面位点，从而使还原反应平衡向左移动，从而抑制钨的沉积[22]。然而，在当前模型框架内，假设沉积反应是不可逆的并且会完全进行，因此没有考虑这种效应。前体消耗率的差异导致B/W和C/W比值沿反应器长度呈下降趋势。为了在扩展的基底上获得成分均匀的涂层，我们必须建立能够最小化边界层厚度并平衡反应物分布的流动条件。

实验数据和模型数据都表明，在W–C–B系统中可以区分出三种特征性的相形成阶段，这些阶段由基底表面的局部前体比率决定：
第1阶段：β-W相的形成。这一阶段发生在低前体流量下（Q(WF6) = 1 L/h且Q(TMAB) < 0.6 L/h）。在这些流量下，计算出的B/W比值相对较高（0.8–3.6），理论上应该促进硼化物的形成。然而，WF6的低绝对流量显著降低了钨的生长速率。这导致表面动力学发生改变：TMAB分解产生的含碳和氮物种被生长中的金属层动态捕获，稳定了亚稳态的β-W结构。因此，这一阶段受到动力学饥饿的影响，杂质驱动的稳定性超过了热力学相平衡。
第2阶段：硼化物的竞争性生长。这一阶段在前体钨流量增加（Q(WF6) = 3 L/h）和中等TMAB进料速率（Q(TMAB) = 0.14–0.68 L/h）时建立。在这些条件下，钨的通量足以抑制N和C杂质的影响，同时B/W比值保持在0.2–1.0范围内。这一阶段的关键驱动因素是硼吸附相对于碳的动力学优势。尽管C/W比值可以超过1（0.6–3.2），这通常表明有利于碳化钨的生长，但由于BH3物种的无障碍吸附，形成强W–B键在动力学上更为有利。因此，钨单硼化物（WB）的形成占主导地位，从而抑制了β-W和碳化钨的形成。
第3阶段：碳化物的主导。这一阶段在前体TMAB流量较高（Q(TMAB) > 1.7 L/h）时出现。在这个阶段，计算出的C/W比值达到临界值（>9），在表面产生了大量的碳。高浓度的碳可能会饱和活性表面位点，从而阻止含硼前体的吸附。从热力学角度来看，碳化钨（WC, W2C）是稳定相，因此在高碳活性条件下它们的形成成为主要反应途径。此外，显著的碳过剩与实验观察到的颗粒细化和涂层结构非晶化倾向相关。在研究碳化钨的过程中也报告过类似的现象[41,42,43]。

由于碳化钨的化学气相沉积（CVD）合成已经被广泛研究，并且通常使用更容易获得的碳源如甲烷（CH4）、丙烷（C3H8）或乙炔（C2H2）[44,45,46,47]进行，因此从WF6和TMAB合成碳化钨的实际应用有限。相反，能够在相对较低的温度（550 °C）下合成稳定的钨硼化物（WB）相具有极大的兴趣。为了准确界定加工窗口并确定从WF6–TMAB混合物中沉积钨硼化物的最佳条件，解决了一个多目标优化问题。入口前体流量被选为控制变量（方程（30）：
X = Q(WF6), Q(TMAB) （30）
其中：
X是过程控制变量的向量；
Q(WF6)是六氟化钨的流量，单位为L/h；
Q(TMAB)是TMAB的流量，单位为L/h。

优化标准是基于模型推导出的B/W和C/W比值（表8）建立的。表8显示了优化标准。优化问题被构建为最小化一个二次损失函数L(X)，该函数限制了涂层成分与目标WB相的化学计量比的偏差，同时考虑了两个主要条件：基底上的B/W和C/W比值。由于我们之前的研究表明，作为硼源的BH3和作为碳源的TMA在气相中都受到类似的扩散限制，这意味着基底上的B和C的比例将由气态前体（TMAB）的化学计量比决定，即B:C ≈ 1:3。因此，方程（31）的第一项对偏离WB相最佳化学计量比（1:1）的系统进行惩罚，而第二项对偏离化学计量碳比值的情况进行惩罚，后者本质上是前体性质的固有特征。因此，目标函数最小化了相对于目标化学计量平衡点的方差（31）。

优化过程受到钨前体流量下限约束，以避免形成亚稳态的β-W相（32），根据实验观察，在低钨沉积速率下，该相会发生杂质驱动的稳定。WF6流量的搜索空间下限设置为Qcrit = 2.0 L/h。这个值是根据实验数据得出的：当Q(WF6) ≤ 1.0 L/h时，在沉积的涂层中检测到了不希望出现的亚稳态β-W相，而当Q(WF6) ≥ 3.0 L/h时则没有观察到这种现象。选择Qcrit = 2.0 L/h作为这个不确定性区间的中点，从而为稳定生长的无β-W钨硼化物提供了安全边际。

优化问题使用COMSOL Multiphysics软件包中的参数扫描（Parametric Sweep）方法在以下参数范围内解决：Q(WF6) = 2–3 L/h和Q(TMAB) = 0.1–3.2 L/h。优化问题解决方案的图形解释以CVD相图的形式展示在图7中。图7显示了前体流量优化的CVD相图。对目标函数L(X)响应面的分析表明，全局最小值（最优值）在参数Q(WF6) = 2 L/h和Q(TMAB) = 0.48 L/h时达到。在这一点上，计算出的B/W比值约为1，C/W比值最小（约为3.1），这对应于最有利于WB相形成的条件。在Q(WF6) = 2–3.0 L/h和Q(TMAB) = 0.48–0.69 L/h的范围内观察到接近最优解的簇（局部最小值，用绿色标出）。因此，建议将这个特定范围作为钨硼化物（WB）沉积的加工窗口，因为它在保持目标相成分的同时，能够抵抗流量波动的影响。

**结论**
在本研究中，对使用WF6和三甲基胺硼烷（TMAB）前体的W–C–B涂层的CVD过程进行了全面的多物理场建模。开发的2D轴对称模型将耦合热传递和质量传递与表面动力学结合起来，确定了工艺参数与涂层相组成之间的关系。研究表明，前体向基底的传输处于混合扩散-动力学阶段。这可能导致反应气体混合物在反应器内的流径上反应物消耗不均匀，这一因素在放大工艺规模时必须加以考虑。将计算出的B/W和C/W摩尔比值与实验数据进行了比较，显示出很强的相关性。该模型准确再现了从亚稳态β-W相到钨硼化物，再到碳化钨的转变过程。

基于多目标优化问题的解决方案，构建了一个沉积阶段图，界定了钨单硼化物（WB）合成的加工窗口。结果表明，使用TMAB作为B和C的唯一前体会对合成涂层的相纯度造成限制，因为前体的化学计量比是固定的（C:B = 3:1）。在高TMAB流量下，碳的过量是抑制硼化物生长的主要因素。为了扩大硼化物的稳定性窗口并可靠地合成富硼相（W2B5, WB4），建议考虑替代的无碳硼源，如二硼烷（B2H6）和三氯化硼（BCl3）。