使用反应性分子动力学方法模拟催化复合反应时人为提高气体压力的局限性

《International Journal of Heat and Mass Transfer》：Limitation of artificially elevating gaseous pressure in simulating catalytic recombination using reactive molecular dynamics method

　　

在高速飞行器的气 - 面相互作用分子动力学模拟中，通常采用人工提高气体压力的方法来提高计算效率，这可能导致相互作用机制或复合系数与飞行或地面试验条件下的情况产生偏差。为了揭示气体压力对表面催化复合的影响模式并寻找模拟期间提高气体压力的物理依据，本研究分析了α-石英上解离氧原子的非均相催化反应。统计结果表明，人工提高的气体压力改变了主导基元步骤，从吸附分子的热脱附转变为 Eley-Rideal 复合。同时，在极高压力条件下，由于气相原子碰撞的发生，速率系数随气体压力的变化规律显著不同。根据气体压力对速率系数的影响，定义了基于气体区域高度的微观克努森数，以阐明在气 - 面相互作用分子动力学模拟中提高气体压力的局限性。本研究加深了对高焓解离原子催化复合机制的理解，并为微观动力学演化结果应用于宏观热量和质量传输提供了质量保证。

1. 引言
解离原子的催化复合显著影响高速飞行器表面的热流 [1], [2], [3]，其中热防护材料作为催化剂加速了复合过程。寻找基元反应对于构建催化复合模型至关重要 [4], [5], [6]，例如原子吸附或脱附、气相原子与吸附原子之间的 Eley-Rideal (ER) 复合以及吸附原子之间的 Langmuir-Hinshelwood (LH) 复合。基于特定力场的分子动力学 (MD) 模拟有助于揭示反应网络，因为它具有原子级分辨率 [7], [8], [9]。量化热防护材料实际服役条件下的反应速率是模型构建的另一个关键方面 [10]。因此，在使用指定的统计方法从 MD 模拟结果中获得速率结果时，需要足够长的演化周期。

为了缩短气 - 面相互作用系统的演化时间，气体压力通常被人为提高到远高于边界层中的值。例如，Cozmuta [11] 在 MD 计算中使用了 10.0–100.0 atm 的压力以确保发生足够多的反应；Norman 团队 [12], [13] 在 100 atm 的压力下进行建模；Yang [14] 在其研究中采用了 10.0 atm 的压力以缩短计算时间，这高于实际压力；Cui 等人 [15] 通过确保不同模拟案例中具有相同的入射频率来保持恒定的压力。以往关于催化复合速率的数值研究很少考虑压力对反应的影响。热防护材料表面氧原子复合系数的地面试验通常发生在低于 1.0 atm（1.0 atm 为标准大气压）的条件下 [16]。在实际高速飞行期间，表面压力通常在 0.01 到 1.0 atm 之间 [11]。

然而，研究表明气体压力可能会影响原子的复合系数或基元反应的速率系数 [17]。Pidan 等人 [18] 测量了不同驻点压力下的氧原子复合系数，表明随着压力增加，复合系数下降。Cheung 等人 [19] 实验研究了氮原子复合系数随温度和压力的变化，尽管测试的压力范围相对较窄（10? Pa 量级）。在他们对表面上解离原子催化复合实验研究的综述中，Massuti-Ballester 等人 [20], [21] 得出结论，氮原子和氧原子的复合系数均表现出随压力增加而下降的趋势。然而，地面试验设施可达到的压力范围通常有限。Yang 等人 [22] 也通过实验方法获得了不同压力下的氧原子复合系数，并观察到系数随压力升高而下降。这些发现共同表明，解离原子的复合系数会受到气体压力条件的影响，但由于实验条件的限制，仅通过实验方法很难定量评估人为提高的气体压力对界面催化反应的影响。

此外，速率系数的压力依赖性可能会在不同的压力范围内发生变化。考虑到 MD 模拟涉及的尺寸极小，原子模型中的气体属于自由分子流 [23], [24]，这是利用气体动理学计算气 - 面碰撞频率的基础。在过高压力下，由于显著的气相碰撞开始发生，速率系数的推导变得无效。正如 Panerai 和 Chazot [17] 所指出的，微观模拟，特别是在低压条件下 [11] 的模拟，可能有助于揭示热防护材料与化学非平衡气体之间的相互作用机制和规律。

在本研究中，通过反应分子动力学 (RMD) 模拟和指定的统计方法，分析了气体压力对解离氧原子表面催化复合的影响，特别是对基元反应速率系数的影响。第 2 节介绍了 RMD 模拟设置和用于统计反应速率的统计方法。第 3 节分析了不同气体压力下的 RMD 模拟结果。相应地，第 4 节提出了在气 - 面相互作用 RMD 模拟中人为提高气体压力的局限性。最后，第 5 节给出了本研究的结论。

2. 原子模拟与统计方法
2.1. 氧 - 二氧化硅界面反应系统的 RMD 模拟
基于 ReaxFF 力场的 RMD 方法能够模拟化学键的形成和断裂，并已被广泛用于研究化学反应流与热防护材料之间的气 - 面相互作用效应。本工作采用 ReaxFF_GSiSiO 力场 [25] 模拟α-石英上解离氧原子的催化复合。所有 RMD 模拟均使用 LAMMPS 软件包 [26] 进行。催化复合反应系统的原子模型如图 1 所示。模拟盒尺寸为 73.74 ? (x) × 72.37 ? (y) × 26.0 ? (z)，其中 x 和 y 维度对应于α-石英晶格的周期性。下半部分 (0.0 ? ≤ z ≤ 13.0 ?) 填充了α-石英晶格。顶部 1.0 ? (25.0 ? ≤ z ≤ 26.0 ?) 作为原子注入区。在 0.0 ? ≤ z ≤ 2.0 ? 区域内的四层原子被固定以保持体晶格结构，而在 2.0 ? ≤ z ≤ 5.0 ? 区域内的五层原子使用 Langevin 恒温器进行热化。z = 5.0 ? 以上的原子可以自由反应。使用 Verlet 算法进行时间积分，时间步长为 0.5 fs。

图 1. RMD 模拟的原子模型。

总模拟持续时间随压力变化：在 1000 K 时，模拟时间从 80 atm 时的约 10 ns 到 0.625 atm 时的 100 ns 以上。在所有情况下，系统演化直至表面吸附物种群达到稳态，确保有足够的反应统计量进行可靠的速率拟合。在气 - 面相互作用之前，系统在 NVT 动力学下松弛了 40.0 ps，然后在气 - 面相互作用模拟期间采用 NVE 动力学。

入射气体原子在与表面恒温器相同的温度下生成；未对气 - 面界面建模温度梯度，这与气体 - 气体碰撞可忽略不计的自由分子流区域一致。入射原子氧在注入区生成，其速度遵循气体温度 ??下的麦克斯韦 - 玻尔兹曼分布。三个笛卡尔速度分量 (????, ????, ????) 按照 Garcia 和 Wagner [27] 的方法生成，其中 ???? 和 ???? 取自高斯正态分布，而 ???? 取自通量加权分布：
$$ \begin{cases} v_x = \sqrt{\frac{k_B T}{m}} R_G \\ v_y = \sqrt{\frac{k_B T}{m}} R_G \\ v_z = -\sqrt{\frac{2 k_B T}{m}} \sqrt{-\ln R_U} \end{cases} \tag{1} $$
其中 ??B 是玻尔兹曼常数，?? 是氧的原子质量，??G 是高斯正态随机数，??U 是 (0, 1] 中的均匀随机数。

注入是一个确定性的恒定速率过程。单位时间内要插入的原子数由入射频率决定：
$$ r_{in} = \frac{n_{AO} \bar{V} A}{4} \tag{2} $$
根据气体动理学使用原子氧数密度 ??AO、平均热速度 ?? 和表面积 ??计算得出。原子在注入区内的随机位置插入，其速度矢量指向表面。对于本研究中考虑的所有温度 (800–1400 K) 和压力 (0.625–80.0 atm)，均应用相同的注入方案，仅 ?? 和 ??AO 的值随预定的热力学条件变化。

基于该系统的模拟发现界面有 10 个基元反应步骤，包括气相氧原子 O(g) 在 Si 位点 (*) 上吸附形成 O(ad) 的 Ad1、氧原子在 O(ad) 上吸附形成 O2(ad) 的 Ad2、O(g) 和 O(ad) 之间的 ER1 复合、O(g) 和 O(bridge) 之间的 ER2 复合、O(g) 和 O2(ad) 之间的 ER3 复合、O(ad) 和 O(ad) 之间的 LH1 复合、O(ad) 和 O2(ad) 之间的 LH2 复合、O2(ad) 和 O2(ad) 之间的 LH3 复合、整个 O2(ad) 分子脱附的 De1，以及从 O2(ad) 脱附一个氧原子的 De3。这些基元反应步骤的表达式列于表 1 中。

表 1. 基元反应步骤的反应路径。

| 编号 | 基元反应步骤 | 反应路径 |
| :--- | :--- | :--- |
| 1 | Ad1 | ??(??????)+(*)→??(????) |
| 2 | Ad2 | ??(??????)+??(????)→??2(????) |
| 3 | ER1 | ??(??????)+??(????)→??2(??????)+(*) |
| 4 | ER2 | ??(??????)+??(????????????)→??2(??????)+(*) |
| 5 | ER3 | ??(??????)+??2(????)→??2(??????)+??(????) |
| 6 | LH1 | ??(????)+??(????)→??2(??????) |
| 7 | LH2 | ??(????)+??2(????)→??2(??????)+??(????) |
| 8 | LH3 | ??2(????)+??2(????)→??2(??????)+2??(????) |
| 9 | De1 | ??2(????)→??2(??????)+(*) |
| 10 | De3 | ??2(????)→??(??????)+??(????) |

ER 复合、LH 复合和 De1 脱附直接形成气相分子 O2(gas)，因此它们对总复合系数有贡献。基元路径 ?? 的复合系数 ??O,?? 定义为通过路径 ?? 复合成分子氧的入射原子氧的比例：
$$ \gamma_{O,i} = \frac{2 r_{O2,i}}{r_{in}} \tag{3} $$
其中 ??O2,?? 是通过路径 ?? 形成 O2 分子的速率（通过对累计反应事件与时间进行线性拟合获得），??in 是原子氧的入射频率 (式 (2))。系数 2 考虑到每个复合事件消耗两个氧原子这一事实。因此，?? 是一个介于 0 和 1 之间的无量纲效率，代表入射氧原子通过特定路径发生复合的概率。

此外，O(bridge) 稳定地结合在 SiO2 表面的两个 Si 原子上，难以与撞击的气相氧原子复合 [12], [13], [28], [29]，导致在氧原子催化复合的一些建模研究中忽略了 ER2 复合和 O(bridge) [30], [31], [32]。当前研究列出了 ER2 复合的贡献，但不会考虑 O(bridge) 的形成和消耗。

2.2. 基元反应速率的统计方法
根据上述基元步骤的特征，开发了一种自动跟踪基元反应的算法 [33]。该算法包括三个主要过程。首先，识别表面吸附物并记录最近 ??t 个时间步的信息。然后，捕获当前时间步吸附物或复合分子（即反应事件）的形成。最后，追溯过去 ??t 个时间步中反应物的状态，以匹配数据库中列出的基元步骤特征。通过这种方式，可以找到每个反应事件的路径。RMD 模拟提供了大量的位置、速度和键合数据。统计算法有助于从这些结果中提取两类重要信息：任意时刻的吸附物数密度和累计反应事件数。吸附物数密度稳定表明反应系统动态稳定，而累计反应事件数随时间的斜率意味着反应速率。

表 1 中提到的所有基元反应步骤的速率系数都可以从气相原子的入射频率、吸附物数密度和反应速率推导出来。区分复合系数 ?? 和基元步骤的速率系数 ???? 很重要。速率系数 ????（单位：s?1）表征基元反应的固有动力学——每单位时间每反应位点的反应事件数（对于诸如脱附等单分子步骤）或每反应对的反应事件数（对于诸如 LH 复合等双分子步骤）。虽然 ?? 既取决于固有反应动力学又取决于气动环境（压力、入射通量），但 ???? 仅反映反应物理学，并且在克努森流条件下，与入射通量无关。这种区分至关重要：?? 是从所有基元步骤的耦合动力学中涌现出的系统级可观测量，而 ???? 是单个步骤的基本参数。

速率系数与这些物理量之间的具体关系因化学反应表达式而异。

对于气-吸附体反应，即吸附和ER复合，反应步骤??的速率系数是入射频率??in与平均反应概率???????的乘积（4）：{??Ad,??=??in·???????Ad,??，??ER,??=??in·???????ER,??}，其中??in根据气体动理论，由平均热速度???????、氧原子数密度??????和表面积??计算得出（5）：??in=??????·???????·??/4，而???????可通过气-表面相互作用体系的RMD模拟统计获得，如式（6）所示：?????????=????/(??in·??????,??)。在式(6)中，????为基元反应步骤??的反应速率，??????,??表示作为步骤??反应物的表面吸附体数量，从每个时间步的RMD轨迹中统计得出。具体而言：对于ER1复合反应[O(g)+O(ad)→O?(g)+(*)]，??????,??为顶吸附氧原子数??O(ad)；对于ER3复合反应[O(g)+O?(ad)→O?(g)+O(ad)]，??????,??为顶吸附氧分子数??O?(ad)；对于Ad1吸附反应[O(g)+(*)→O(ad)]，??????,??为空吸附位点数??(*)；对于Ad2吸附反应[O(g)+O(ad)→O?(ad)]，??????,??为O(ad)吸附体数。该数量通过以下描述的基础上键阶跟踪算法在每个输出时间步识别表面吸附体而统计获得。

表面吸附体利用ReaxFF反应力场提供的键阶信息进行区分。在每个输出时间步，通过ReaxFF键阶函数从原子间距离评估所有原子对之间的键阶。分类标准如下：

•O(ad)——顶吸附氧原子：一个单独的氧原子，与恰好一个表面Si原子形成Si-O键阶超过键合阈值（键阶>0.3），且与其他任何氧原子之间无O-O键。

•O?(ad)——顶吸附氧分子：一对通过O-O键连接的氧原子（键阶>0.3），且其中至少一个原子还与表面Si原子成键。

•O(bridge)——桥吸附氧原子：一个单独的氧原子，与两个不同的表面Si原子形成Si-O键。

•(*)——空吸附位点：配位不足的表面Si原子。

该分类在每个输出时间步应用，生成吸附体种群的时间序列，用于稳态检测、反应速率提取以及反应概率和速率系数的计算。

LH复合已被证明发生在最近邻吸附体之间[32]，[33]，[34]，因此需要根据RMD模拟假设LH复合可以发生时的最大距离??LH。考虑到计算表面上最近邻吸附体的平均数量，LH复合的速率系数可推断为（7）：

{??LH1 = 2·??·??LH1/(??O(ad)2·π·??LH12)
??LH2 = ??·??LH2/(??O(ad)·??O?(ad)·π·??LH22)
??LH3 = 2·??·??LH3/(??O?(ad)2·π·??LH32)}

热脱附反应的速率系数也由反应速率及其反应物的数密度计算得出（8）：{??De1 = ??De1/??O?(ad)；??De3 = ??De3/??O?(ad)}

有效截断半径??LH从RMD模拟中观察到的LH复合事件中估计。对于每个此类事件，记录复合反应前两个反应吸附体之间的距离；在所有条件下，这些距离均落在最近邻吸附体间距范围内。采用代表性值??LH = 2.5 ?。由于??LH作为乘数常数进入LH速率系数（??LH,?? ∝ 1/??LH,??2），因此??LH,??的压力依赖趋势——即核心关注量——对??LH的具体选择不敏感。

除有效LH半径??LH,??外，式(7)和(8)中的速率系数均通过自动反应跟踪算法[33]对RMD轨迹后处理获得。该算法通过ReaxFF键阶识别表面吸附体，检测反应事件（成键或脱附），并追溯参与原子在反应前的状态，将每个事件分配给相应的基元路径。反应速率通过对稳态采样窗口内累计事件计数进行线性回归获得，反应物数密度在同一窗口内取平均值。所有速率系数的单位均为s?1。

式(4)、(5)、(6)、(7)、(8)表明，所有基元反应步骤的速率系数均可从RMD模拟的统计结果中提取。值得一提的是，只有在气体流动处于Knudsen区域（即自由分子流）时，气体-表面碰撞频率才可通过气体动理论按式(5)计算，此时气体原子仅与固体表面相互作用而不与其他气体原子相互作用。在RMD模拟中，反应体系的特征长度为纳米量级，因此只要气体不过于稠密，即处于自由分子流区域。

3. 气相压力对催化复合的影响

3.1. 氧原子入射频率

为研究氧原子压力??O对不同温度下非均相催化复合的影响，分别在800 K、1000 K、1200 K和1400 K下进行压力效应分析。对于每个温度条件，基于ReaxFF力场在8个??O值下进行了RMD计算：0.625 atm、1.25 atm、2.5 atm、5.0 atm、10.0 atm、20.0 atm、40.0 atm和80.0 atm（其中1.0 atm = 101,325.0 Pa）。在每个温度和压力条件下，经过足够的气-固原子碰撞后——以表面吸附体数量达到动态稳态且基元反应计数与时间呈强线性关系为标志——对基元反应速率进行拟合。随后，计算表面吸附体覆盖率和基元步骤的速率系数，以实现对??O对反应机制影响的定量分析。

式(5)表明入射频率由???????和??????决定，它们与??和??O相关，如式(9)所示。在特定温度下，??in与??O成正比。使用当前方法在不同温度和压力条件下计算的??in值如图2所示。

(9) ??????? = √(8????/(π??))
(10) ??AO = ??A·??O/(??·??)

下载：下载高分辨率图片 (70KB) 下载：下载全尺寸图片
图2. 不同温度和压力条件下的入射频率。

3.2. 气相压力对复合系数的影响

图3(a)～(d)分别展示了四种温度条件下氧原子复合系数??O随压力??O的变化。不同基元反应引起的复合系数以不同颜色显示。图3表明，各温度条件下的总复合系数??O,total随??O变化相对较小。当??O增加到一定程度后，随着压力进一步升高，??O,total略有下降。此外，表面温度越高，高压条件下??O,total的下降越显著。尽管??O,total变化相对较小，但各基元路径的复合系数受气相压力显著影响。特别是起主导作用的ER3复合和De1分解，分别随??O增加呈递增和递减趋势。

下载：下载高分辨率图片 (318KB) 下载：下载全尺寸图片
图3. 气相压力对基元复合系数的影响。

两种关键中间产物的吸附氧原子覆盖率??O(ad)和氧分子覆盖率??O?(ad)随??O的变化如图4所示。在四种温度条件下，??O(ad)和??O?(ad)随??O的变化模式高度一致。当气相压力从0.625 atm增加到10.0 atm时，??O(ad)和??O?(ad)均显著增加；当??O进一步增加时，??O(ad)和??O?(ad)的增长趋势减缓。在1400 K下，当??O增加到80.0 atm时，??O(ad)开始呈下降趋势。这是因为随着??O增加，气-固原子之间的碰撞和吸附反应加速，导致表面吸附体生成速率更高。然而，ER复合、LH复合和热脱附反应对吸附体的消耗速率不能同比增加，使得高??O条件下表面吸附体的数密度得以累积。

下载：下载高分辨率图片 (222KB) 下载：下载全尺寸图片
图4. 气相压力对表面覆盖率的影响。

在消耗表面吸附体的反应中，ER复合速率直接受氧原子入射和碰撞频率的影响，因为气相氧原子作为必要反应物参与反应。因此，ER复合引起的复合系数随??O增加而增大（如图3所示），这意味着ER复合对吸附体的消耗速率随??O增加而显著增大。然而，LH复合和热脱附不受氧原子入射频率的直接影响，尤其是占比较大比例的De1反应，随??O增加而减小。即，它们对吸附体的消耗速率随??O增加而减小。因此，在热脱附主导的低??O区域，吸附体覆盖率随??O增加而急剧上升，而在ER复合逐渐主导的高??O条件下，??O(ad)和??O?(ad)的增加则趋于平缓。

在上述四种温度条件下，ER复合、LH复合和热脱附引起的??O占总??O,total的比例??随??O的变化如图5所示。??ER随??O增加而逐渐增大，主要由于氧原子入射频率的增加。??LH随??O增加先升后降，但整体变化相对较小。热脱附的贡献随??O增加而显著减小。

下载：下载高分辨率图片 (362KB) 下载：下载全尺寸图片
图5. 气相压力对三条反应路径占比的影响。

需要注意的是，LH路径的复合系数定义为??LH = 2??LH/??in。随着气相压力增加，LH反应速率??LH（由更高表面覆盖率驱动）和入射频率??in（由更高气体密度驱动）同时增加。因此??LH随压力近似不变反映了两个量的比例缩放关系，而非LH动力学对表面覆盖率的无敏感性：LH复合的绝对速率确实随压力增加，与吸附体种群数的预期依赖关系一致。

图5还表明，在800 K温度下，在计算??O范围内，ER复合始终占主导地位。在1000 K温度下，当??O降至0.625 atm时，De1脱附效应与ER复合效应相当。而在1200 K和1400 K温度下，当??O降至某一值时，De1脱附超过ER复合并成为主导反应路径。因此，??O对主导反应路径的影响受计算系统温度条件的影响。当温度条件足以触发的脱附反应时，??O降低将导致主导反应路径从ER复合转移至De1脱附。换言之，人为提高??O以改善一定数量总反应事件的计算效率，将导致主导反应路径出现偏差。

3.3. 气相压力对速率系数的影响

3.3.1. 气相与吸附相之间的反应

基元复合系数表征原子复合形成分子的概率，由反应本身的速率系数和气动力环境决定，是反应体系综合演化的结果。速率系数表征反应的本质特征，是建立界面反应动力学模型的关键参数。除温度和表面材料外，研究表明速率系数还受气相压力影响[17]，[35]。因此，研究气相压力对速率系数的影响规律具有重要意义。基于第2.2节所述的反应速率系数获取方法，采用RMD演化结果的稳定阶段进行基元反应速率系数的计算。

在四种温度条件下，Ad1反应的平均吸附概率???????Ad1随??O的变化如图6(a)所示。随着??O增加，???????Ad1减小。图6的横坐标为??O^(1/2)。这是因为在上述所有计算条件的RMD计算结果中，???????Ad1与??O^(1/2)呈良好的线性关系，且对于每个温度条件，???????Ad1与??O^(1/2)之间的线性关系可分为斜率不同的两个区域。为客观确定拐点，测试了四个候选压力——2.5、5.0、10.0和20.0 atm，每个候选压力两侧至少保留三个数据点。在每个候选压力??O,brk处，对反应概率??关于??O^(1/2)进行约束双段线性拟合，要求两段在??O,brk处相交。

对于每次拟合，拟合质量通过实测概率值与拟合值之间的差平方和进行评估，即对所有八个数据点求和∑(?)2。选取使该和最小的候选压力作为最佳断点。

Ad1的最优断点结果分别为5.0 atm（800 K）、2.5 atm（1000 K）、10.0 atm（1200 K）和5.0 atm（1400 K）。对于Ad1，两段拟合始终比单一直线更好地描述数据，确认了两段行为是数据的真实特征。在小于O,brk的范围内，Ad1随1/2O减小的趋势更为明显。当O大于O,brk时，下降趋势减缓。图6(b)展示了1000 K温度下，高低O范围内Ad1与1/2O的线性拟合结果，两者的斜率差异显著。

值得注意的是，在低压区域（O≤10.0 atm），1400 K时的Ad1吸附概率超过1200 K时的吸附概率（图6(a)）。在低压下，1400 K时空吸附位点数(*)低于1200 K，这导致尽管绝对吸附效率Ad1随温度升高而增大，表面单位位点的概率仍然被放大。在高压下（O≥20.0 atm），1400 K时的(*)超过1200 K时，常规趋势得以恢复。因此，这种反转是归一化方式依赖温度所致，而非物理吸附倾向的真实逆转。

与Ad1反应类似，Ad2反应的平均吸附概率Ad2也随O的增加而减小，如图7(a)所示。采用相同的断点选择程序，Ad2的最优断点在800 K、1000 K和1200 K时为5.0 atm，在1400 K时为10.0 atm。Ad2在低气压范围内的变化幅度较大，而在高O范围内的变化幅度较小。与Ad1吸附不同，Ad2在四个温度条件下的值相近，表明Ad2主要受压力影响，而非温度。

图8(a)显示，在不同气压条件下，ER1复合反应的平均反应概率ER1随O的增加而逐渐增大。采用相同的断点选择程序（在候选压力2.5、5.0、10.0和20.0 atm上最小化∑(?)2），ER1的最优断点在800 K和1000 K时为2.5 atm，在1200 K和1400 K时为10.0 atm。在所有四条气体-吸附物反应路径中，两段拟合对ER1的提升幅度最大。在O小于O,brk的范围内，ER1与1/2O呈近似线性增长。当O继续增加时，ER1以较小的斜率线性增长。例如，在1200 K温度下，以10.0 atm为边界点，ER1与1/2O的分段线性拟合结果如图8(b)所示。

图9(a)展示了ER3复合反应的平均反应概率ER3随气压的变化规律。四个温度条件下的ER3值均随O的增加而减小，且四个温度条件下的ER3-1/2O关系相似，表明ER3主要受O影响，而温度影响相对较弱。采用相同的断点选择程序，ER3在四个温度下的最优断点均为5.0 atm，如图9(b)所示。以5.0 atm为边界点，ER3与1/2O呈现两段具有不同斜率的线性关系，且仍保持ER3在低压条件下随O变化更为剧烈的特征。以1000 K的计算结果为例，图9(b)展示了ER3与1/2O的分段线性拟合结果，说明了线性拟合结果斜率之间的差异。

Ad1吸附、Ad2吸附、ER1复合和ER3复合均涉及气相原子与吸附相原子或位点之间的相互作用。这些步骤的平均相互作用概率与气压呈现显著关系，即i;与1/2O呈线性关系，且在低压和高压区域内的斜率差异显著。其原因是，在RMD模拟中，反应概率等其他结果基于局部克努森流的假设推导得出，这意味着系统中仅有气固碰撞。

当气压较低时，气相入射原子能够充分与吸附位点或吸附粒子碰撞和相互作用，这属于克努森流的范畴。随着气压升高，气相原子的入射频率增加，原子模型气相区域内的原子数量增多，导致气相原子间的相互作用逐渐凸显，并偏离克努森流。这种偏离体现在i随O的变化逐渐偏离低压条件下的趋势。因此，出现了图6至图9所示的分段线性函数规律。

在地面高焓试验中，很难创造高压条件，O通常低于1.0×105 Pa量级。在真实飞行状态中，激波解离后的原子分压通常也小于1.0 atm。因此，试验和飞行条件下热防护材料与来流原子流的相互作用均在克努森流有效的压力范围内，可以通过扩展低压范围内的微观模拟结果来进行表征和建模。

3.3.2 吸附相之间的反应

除了气相与吸附相之间的四种基本反应路径外，LH1、LH2和LH3复合反应属于吸附相原子之间的相互作用，而De1和De3是吸附物热运动导致的脱附效应。它们的反应速率不直接受入射原子频率的影响。然而，吸附相原子的热运动受入射原子及其运动的影响，因此有必要分析O对其反应速率系数的影响。

对3.3.1节描述的相同断点选择程序应用于本节的所有五个基本步骤（在候选压力2.5、5.0、10.0和20.0 atm上最小化∑(???)2）。对于LH1和De3，两段拟合比单一直线有明显的提升，最优断点为10.0–20.0 atm。对于LH2、LH3和De1，两段拟合在最优断点两侧的斜率几乎相同，表明这些速率系数在整个压力范围内与1/2O基本呈单一线性关系。因此，采用O,brk=10.0 atm来展示LH2、LH3和De1的两段拟合，而LH1和De3则采用最优断点。

RMD计算结果表明，在四个温度条件下，LH1复合反应的速率系数LH1随O的增加而增大，如图10(a)所示。在O小于O,brk的范围内，LH1与1/2O具有良好的线性关系；当O大于O,brk时，LH1仍与1/2O呈线性相关，但斜率增大。以1000 K的结果为例，两个压力区域的拟合结果如图10(b)所示。因此，当气态偏离克努森流区域，即接近连续流区域时，气压对LH1复合反应速率系数的影响模式发生改变。

LH2复合反应的速率系数LH2在当前O范围内随1/2O呈良好的线性增长关系，且斜率没有显著变化，如图11(a)所示。以O=10 atm分别对LH2和1/2O进行线性拟合，两者的斜率差异较小。图11(b)展示了800 K条件下计算结果的分段线性拟合结果，其斜率变化相对显著，表明气体碰撞对LH2复合反应的影响相对较小。

图12(a)展示了气压对LH3复合反应速率系数LH3的影响。LH3也随O的增加而增大。两段拟合在最优断点两侧也产生了几乎相同的斜率，表明速率系数在整个压力范围内与1/2O基本呈单一线性关系。上述分析表明，高O条件导致LH3~1/2O的变化规律相对于克努森流区域略有偏差。

图13(a)表明，De1脱附反应的速率系数De1随O的增加呈稳定上升趋势。在计算的O范围内，四个温度条件下的De1均与1/2O保持良好的线性关系。以1000 K的结果为例，分别在O≤10.0 atm和O>10.0 atm的气压范围内对De1与1/2O进行线性拟合，拟合曲线如图13(b)所示。两条拟合曲线的斜率极为接近，表明De1脱附（即吸附的O2分子整体脱附）的速率系数受气体碰撞效应的影响较小。

De3脱附反应的速率系数De3也随O的增加而增大，且De3随O增加呈现两个显著不同的增长趋势，如图14(a)所示。在O≤10.0 atm的范围内，De3随1/2O的线性增长较为平缓，而在O>10.0 atm的范围内，增长趋势变得更加陡峭。图14(b)中1000 K的拟合结果充分展示了高低O范围内线性拟合结果的显著差异。图14表明，De3脱附（即吸附的O2分子中一个原子的脱附）的速率系数受气体碰撞效应的影响很大。

O对非均相催化反应速率的影响模式表明，O显著影响基本反应路径对总复合系数的贡献。此外，O对基本反应速率系数的影响揭示了人为提高气压的局限性。

4. 人为提高气压的局限性

4.1. 微观克努森数的定义

3.3.1节气相与吸附相之间的反应和3.3.2节吸附相之间的反应的分析表明，在特定温度条件下，当O足够低时，吸附反应和ER复合反应的平均反应概率均与1/2O近似呈线性关系。在此压力范围内，LH复合和热脱附的速率系数也与1/2O呈现稳定的线性关系。可以看出，随着O的降低，基本反应步骤的速率系数随压力变化的趋势趋于稳定，表明反应系统中气固原子碰撞逐渐接近克努森流区域。

在克努森流条件下，可以定量表征基本步骤的速率系数与O之间的关系。

因此，满足克努森流区条件的原子压力下的计算结果可推广至较低压力条件，因为在那些条件下 RMD 的计算成本过高而不可接受。也就是说，在人为提高的压力条件（但不违反克努森流区）下的速率系数和复合系数的计算结果，可以外推至实验条件或实际飞行条件。为了表征人为提高原子压力的限制，即在克努森流区范围内处理气 - 面相互作用，定义了无量纲数 $Kn_h$ 如下：
(11) $Kn_h = \lambda/h$
其中 $\lambda$ 是气体原子的平均自由程，$h$ 是原子模型中气相区域的高度。$Kn_h$ 可视为克努森数的微观形式。$\lambda$ 决定了气体原子之间的平均距离，而 $h$ 代表了气体原子与表面原子之间的最短距离。当前的 RMD 模拟中采用 $h=13$ ?，该值大于 ReaxFF-GSISiO 力场的截断距离。较高的 $Kn_h$ 值表示气体更为稀薄，这使得气体原子与表面原子之间在进行任何气相碰撞之前能够进行更广泛的相互作用。为确保碰撞仅发生在解离气体与固体表面之间，$Kn_h$ 应足够大。也就是说，当使用 RMD 方法定量表征反应速率并建立界面处的气 - 固相互作用模型时，气相压力（或原子入射频率）的选择应满足 $Kn_h \ge Kn_{ch}$ 的约束条件，从而能够通过外推获得较低压力条件下的速率系数。$Kn_{ch}$ 是纯气 - 面碰撞条件被破坏时的临界克努森数。平均自由程 $\lambda$ 与气体温度 $T$ 和压力 $p_O$ 相关
(12) $\lambda = k_B T / (\sqrt{2} \pi d^2 p_O)$
其中 $k_B$ 是玻尔兹曼常数，$d$ 是有效原子直径。因此，针对克努森流区的限制条件 $Kn_h > Kn_{ch}$ 可表示为
(13) $p_O \le k_B T / (Kn_{ch} \sqrt{2} \pi d^2 h)$
这表明气体压力可被人为提高的程度受到温度、模型尺寸和临界克努森数的限制。

$Kn_h$ 对气相区域高度 $h = 13$ ? 的依赖性具有明确的物理基础。用于从 RMD 轨迹中提取速率系数的统计方法（见 2.2 节）假设，入射气体原子从注入区运动至表面期间，在中间空间内不发生气 - 气碰撞。只有在气体处于自由分子流区（其特征是平均自由程 $\lambda$ 远大于特征输运长度尺度 $h$）时，纯气 - 面碰撞的假设才有效。随着 $h$ 的增加，在给定气体密度下，入射原子在到达表面之前与其他气相原子发生碰撞的概率会增加。因此，纯气 - 面碰撞假设将在相应更低的压力下失效。这反映在 $Kn_h$ 准则中：较大的 $h$ 会在给定压力下降低 $Kn_h$ 值，从而缩小了自由分子流假设成立的压力范围。本研究采用的 $h = 13$ ? 值是物理上相关的长度尺度：它代表了气 - 面相互作用发生之前气相输运所跨越的距离。若采用较大的 $h$ 进行模拟，在相同的 $Kn_h$ 准则下，其压力上限界限将更为严格，但这不会改变 $Kn_h$ 主导自由分子流区与过渡流区之间转变这一基本结论。因此，$Kn_h$ 对 $h$ 的依赖性并非建模伪影，而是连接模拟几何与连续介质尺度克努森数的真实物理特征。

4.2 临界克努森数的获取
为确定 $Kn_{ch}$ 的值，考虑了 Ad2 吸附及其逆反应 De3 脱附，以及 ER3 复合反应。Ad2 吸附、De3 脱附和 ER3 复合的速率系数结果见图 15。图 15(a) 描绘了不同温度下平均反应概率 $\bar{P}_{Ad2}$ 随微观克努森数 $Kn_h$ 的变化情况。在较高的 $Kn_h$ 值下，$\bar{P}_{Ad2}$ 呈现平滑变化，而当进入较低的 $Kn_h$ 区域时，则转变为显著且急剧的变化。这种过渡行为发生在 $Kn_h$ 约为 100.0–200.0 的范围内，该范围已在图 15(a) 中予以高亮显示。
下载：下载高分辨率图像 (214KB) 下载：下载全尺寸图像
图 15. 关键基元步骤的速率系数随克努森数的变化。
图 15(b) 和图 15(c) 中观察到了相似的趋势，表明 $Kn_{ch}$ 的值位于 100.0 和 200.0 之间。以 $Kn_{ch}=100.0$ 为例，在 800、1000、1200 和 1400 K 的温度下，相应的气体压力分别为 8.6、10.7、12.9 和 15.0 atm。由 $Kn_{ch}$ 推导出的过渡压力值与 3.3 节中的分段点接近。气体压力的降低标志着流体状态从以密集原子碰撞为主过渡到以稀疏原子碰撞为特征的状态。当前研究结果揭示了临界微观克努森数 $Kn_{ch}$ 的存在。当 $Kn_h$ 超过该临界值时，基元反应步骤的速率系数与气体压力之间会出现稳定且可解析量化的关系。基于目前的计算结果，该临界 $Kn_{ch}$ 位于 100.0–200.0 范围内，近似值为 150.0。通过施加 $Kn_h \ge Kn_{ch}$ 的约束条件，人为提高气体压力可以在保持克努森流区的同时，提高基于 ReaxFF 的 MD 模拟的计算效率。

5. 结论
本研究采用 RMD 模拟和反应识别统计方法，在不同温度和压力条件下，探讨了气体压力对 α-石英表面上解离氧原子催化复合的影响。此外，拟合了基元反应速率系数对气体压力的依赖关系，从而提出了 RMD 模拟中人为提高气体压力的上限。主要结论如下：
气体压力显著影响基元反应路径对总复合系数的贡献，改变表面反应的演变过程及覆盖率，从而导致界面主导反应路径的变化。在高温下，随着气体压力的增加，主导路径从 ER3 复合过渡为 De1 脱附。
在较高气体压力下进行的分子模拟可能导致基元反应速率系数的压力依赖关系偏离实验或飞行条件下观察到的行为。过高的气体压力会加剧气相相互作用，从而导致无法维持纯气 - 面碰撞条件。因此，分子模拟中应避免采用该范围内的压力设置。
基于基元反应步骤速率系数随气体压力的变化规律，提出了气 - 固非均相催化模拟中人为提高压力的上限约束条件。该约束利用基于气相区域高度的微观克努森数，将计算压力限制在速率系数与压力关系符合真实低压条件下规律的范围之内。
该约束的建立明确界定了在气 - 面相互作用分子模拟中通过提高气体压力来提升计算效率的有效范围。通过在此有效范围内采用最优条件，该方法有助于微观计算数据的积累，并促进其在介观和宏观尺度上的应用。

CRediT 作者贡献声明
Wei Dong：监督、构思。Yanxia Du：监督、构思。Qiheng Chen：写作 – 审阅与编辑、数据整理。Xiaofeng Yang：写作 – 审阅与编辑、构思。Qin Li：写作 – 初稿、方法论、数据整理。
Qin Li|Xiaofeng Yang|Qiheng Chen|Yanxia Du|Wei Dong
中国空气动力学国家重点实验室