二氧化硅表面催化复合机理的原子尺度揭示

2023-10-24贺立超

科技资讯 2023年19期

贺立超

(北京航空航天大学北京 100191)

高超声速飞行器的发展，对于未来国家安全、政治和经济利益有着巨大的战略意义和潜在价值[1]。在高速飞行情况下，激波层和边界层内的气体温度急剧升高，离解后的O、N原子在表面热防护材料的催化效应作用下发生复合反应并释放热量，因而不同的壁面催化特性会极大地影响飞行器表面热流[2]。

采用数值计算方法预测高超声速飞行器表面热环境时，由于材料表面的催化过程较复杂，多数情况下假设材料表面处催化边界条件为完全催化条件或者完全非催化条件[3]。数值计算表明：在飞行器再入过程中，完全催化作用的材料表面热流约为完全非催化表面热流的3 倍。准确预测催化反应速率，可以为高超声速飞行器表面热防护材料的设计提供依据。随着高性能计算系统的飞速发展，反应分子动力学模拟方法已经广泛应用于不同材料的物理性质的研究中，成为了研究表面催化反应过程的新方式[4-5]。

SiO2等硅基复合材料具有强度大、高温稳定性以及化学稳定性好的特点，被广泛应用在高超声速飞行器的热防护中[6-7]。虽然在先前的研究中，已经有研究者通过RMD研究提高了对SiO2表面与O原子碰撞、吸附、催化等过程的认识，但是由于表面附近的环境条件相互耦合，O 原子与SiO2表面的反应过程的影响因素仍需要进一步详细讨论。因而，本文对SiO2表面温度、入射O原子动能这两种单一条件对α-SiO2表面与O原子催化重组反应机制的影响进行研究，以此揭示单一因素对壁面的影响。

1 气固界面反应模型构建

本模拟使用了ReaxFFSiO反应力场，该力场计算结果与DFT 吻合得很好，能较好地再现多种SiO2表面结构，可以应用于大规模Si-O界面的分子动力学模拟。

为了研究O 原子与SiO2表面反应过程，首先建立了一个45 Å×45 Å×240 Å的模拟框，以距离模拟框x、y 方向的底面30 Å 高度的平面为中心，建立了一块沿着（001）表面切割的α-SiO2平板，裸露在平板表面的原子为Si原子。α-SiO2沿z轴方向的尺寸为12 Å，沿x轴和y轴方向的尺寸为45 Å×45 Å，共使用了五层SiO2。为了避免模拟过程中SiO2平板受到O原子的冲击向下移动或变形，距离平板底面0 Å～1 Å 范围内的原子被固定住，距离平板底面1 Å～10 Å 范围内的原子使用Langevin恒温器控制温度，Langevin恒温器的阻尼常数为10 fs，距离平板底面10 Å～12 Å 范围内的原子自由移动。

为了保证O 原子持续与SiO2表面碰撞，使用了通量边界条件。为避免在初始阶段SiO2表面与O原子相互作用，在SiO2表面上方15 Å的平面上，每隔1 ps随机位置产生一个O原子。设立距离表面15 Å以内的区域为反应区域，高于表面15 Å 的区域设立为统计区域，如图1所示。O原子z轴方向速度指向SiO2表面，以保证所有O 原子与表面相撞。每隔1 ps，对位于统计区域的原子统计种类和数目后删除。一方面，可以避免与SiO2表面碰撞后的原子在模型内大量堆积，在长时间的模拟中形成高密度气态区域，改变表面附近环境条件；另一方面，统计区域内原子密度的疏解也可以减少与SiO2表面反应后的产物在统计区域再次反应，进一步精确种类统计结果。

图1 通量边界模型的结构

使用两种环境条件来分析不同条件对O 原子与SiO2表面反应过程：SiO2表面温度的影响和入射O原子初始动能的影响。模拟过程对系统整体使用NVE 系综，在模拟开始前，对SiO2平板进行能量最小化，并使用Nose/Hoover 恒温器对SiO2平板进行50 ps 的充分弛豫，仿真模拟时长均为1 ns，时间步长为0.25 fs。在对SiO2表面温度影响的研究中，O 原子入射方向垂直于SiO2表面，携带动能为0.05 eV。通过改变Langevin 恒温器的温度来获取不同的表面温度，温度改变范围500～2 100 K。在对O原子入射时携带的动能影响的研究中，O 原子入射方向垂直于SiO2表面，SiO2表面通过Langevin 恒温器稳定为1 000 K。O原子入射时的动能变化范围从0.01～4 eV。

2 结果与分析

2.1 表面温度的影响

通过对SiO2表面进行温度控制，笔者对SiO2表面温度的影响进行了研究。如图2所示，经过50 ps的弛豫后，系统温度随着模拟的运行保持稳定；而系统势能在平板弛豫过程中，不同温度下的势能均得到了释放。随着系统模拟的进一步运行，不同温度下的势能都有所升高，但是整体升高幅度较小。这是由于原子和表面碰撞后，碰撞、吸附和复合产生的热能经过热传导，通过Langevin 恒温器重新进行了调整，而离开表面的碰撞后的产物，进入统计区域后被删除，因而本模型很好地维持了系统温度和势能的稳定平衡。

图2 系统温度和势能随时间的变化

如图3所示，在模拟初始阶段，O原子大量吸附在SiO2表面，随着模拟过程继续运行，较低温度下O原子在SiO2表面吸附数目逐渐饱和，而较高温度下O 原子仍在继续附着。一方面，由于较高的表面温度为吸附、复合等反应过程提供了足够的活化能；另一方面，则是随着温度升高，α-SiO2表面逐渐呈现无定形SiO2的形态，提供了更多的反应位点供O原子吸附。

图3 不同温度下的表面O原子吸附特性

在O 原子附着数目逐渐稳定后，更多与SiO2表面碰撞后的产物进入了统计区，进入统计区的O 原子数目初期呈现快速增长的趋势。由于进入统计区域的O原子统计后被删除，因而大幅度减少了气相内的反应，因而O 原子数目不会因为气相内的消耗而减少，O2分子增长速率也相对稳定，由于O3分子键长更长，键能较弱，因而比O2分子不稳定，生成数目较少。通量边界模型删除统计区域的原子，减少气相内反应产物的大量堆积，同样可以减少O3、O4以及其他不稳定的中间分子产物出现在统计区域内的可能性。由图4（a）可知，在相同时刻，随着SiO2表面温度的升高，产生的O2数目随之增长。在获取O2生成数目后，通过下式获取了O2在SiO2表面的催化重组系数。

图4 不同温度下催化重组反应特性

不同时刻温度下的催化重组系数如图4（b）所示，随着表面附着O原子数、O2生成数的稳定，催化重组系数随着模拟的进行逐步达到平衡。催化重组系数随温度升高而升高。

Eley Rideal（E-R）重组机制和Langmuir Hinshelwood（L-H）重组机制是O原子的重组过程占据主导地位的两种机制。如图5（a）所示，≡Si-O·位点和扩散而来的离解原子相撞，形成≡Si-O2，之后Si与O之间的化学键断裂，最终结合成O2分子，这是E-R重组过程。≡Si-O·位点得到了相邻的≡Si-O2的末端O原子，重新形成了一个≡Si-O2，之后Si与O之间的化学键断裂，并通过解吸附，O2分子离开表面，这是L-H 重组过程。通过自编程序的方式，通过对反应过程的轨迹文件进行后处理，对形成O2分子的两个O原子的运动过程进行轨迹追溯，并对不同反应机制下产生的分子分别进行统计，最终获取了不同温度下通过E-R 重机制或L-H 重组机制催化复合形成的O2的比例。由于受到O原子连续不断的轰击，E-R重组机制始终占据着较大比例，在500 K时，通过E-R重组机制重组形成的O2所占比例高达83%，随着温度升高，O原子可以得到更多的活化能，L-H重组机制占据的比例随之有所上升，通过E-R 重组机制重组形成的O2接近70%，具体详见图5（b）。

图5 不同温度下重组机制

2.2 O原子入射动能的影响

笔者通过改变O原子在入射时携带的动能（0.01～4 eV），以此来研究来流O 原子的动能对碰撞、催化重组过程的影响。随着O 原子入射动能的提高，可以为O 原子在表面位点吸附或重组提供的活化能随之增加，因而同一时刻下O 原子在SiO2表面吸附数目也随之迅速升高，具体详见图6。O 原子入射动能0.05 eV的表面附着数目远远小于入射动能1 eV的表面附着数目。当入射动能超过0.5 eV 后，活化能逐渐不再是限制表面附着数目的关键因素，SiO2表面可以提供的位点逐渐趋于饱和，吸附数目不再继续增长。

图6 不同入射动能下的表面O原子吸附特性

如图7（a）所示，在初始阶段，入射原子大量附着在SiO2表面，随着表面位点的减少，O原子吸附数目增长开始放缓。随着表面分布的O 原子越来越密集，相邻位点的O 原子会更容易复合，大量新入射的O 原子和附着在表面的O 原子碰撞的概率也会增加。如图7（b）所示，随着反应的进行，O2分子开始迅速生成。当入射动能较高时，有较多O3进入统计区域，这是由于入射原子携带的大量动能使相对不稳定的O3更容易脱离SiO2表面的吸附进入气相中。

图7 O原子入射动能1 eV时SiO2表面附着数和产物变化

由图8（a）所示，当入射动能小于1 eV 时，随着O原子入射动能增加，重组生成的O2也随之增加，这是由于随着入射动能增加，SiO2表面O 原子附着更加密集，可以用于O 原子表面催化重组的活化能也进一步增加，O原子与位点碰撞概率与催化重组概率均提升。由图8（b）所示，可见当入射O 原子动能小于1 eV 时，随着O 原子入射动能的增加，催化重组系数也随之升高。当入射O 原子动能大于1 eV 时，随着O 原子入射动能增加，O2生成数、催化重组系数都随之下降。这主要是由于较高的动能会使入射的O原子弹离表面或使吸附在SiO2表面的O 原子直接解吸附，而不是继续发生催化重组反应。因此，O原子与表面位点碰撞时，过高的动能可以降低E-R催化重组反应效率。