MOLECULAR DYNAMICS STUDY OF GAS-SURFACE INTERACTIONS IN UPPER ATMOSPHERE
-
摘要: 海拔100 ~ 300 km的上层大气层内, 主要成分之一的原子氧化学活性高, 容易在超低轨飞行器表面材料上发生吸附和侵蚀现象, 显著改变壁面的物理化学状态和气固相互作用特征, 进而影响飞行器的气动性能. 由于时空尺度差异和成本限制, 目前的地面试验手段和模拟方法很难复现上层大气层真实空间环境的气固相互作用过程. 针对这一难题, 文章设计了一种分子动力学等效模拟方案, 让气体分子在短时间内高频入射材料表面, 来反映真实环境中长时间低频率的气固相互作用过程, 从而研究上层大气层不同海拔高度的空间环境中, 原子氧在典型航天材料壁面的吸附/侵蚀状态演化过程以及动量适应系数的变化规律. 研究表明, 在满足一定等效原则时, 这种模拟方案能够合理地反映相关的真实物理图像, 展示了洁净壁面从少量吸附到大量吸附原子氧, 再到原子氧侵蚀壁面表层结构形成松散氧化物的3阶段演化过程, 同时也给出了高速入射的气体分子在相应壁面状态下的动量适应特征. 另外, 本文也研究了上层大气层另一主要组分N2的存在对原子氧吸附、侵蚀和适应过程的影响. 这一研究可加深人们对上层大气层气固相互作用过程的认识, 并为新一代飞行器的气动设计提供参考.Abstract: In the upper atmosphere 100 ~ 300 km above sea level, the atomic oxygen, one of the main constituents, is chemically active and can be easily adsorbed and even corrode surfaces of vehicles flying in very low earth orbits. This phenomenon will significantly change the physical and chemical properties of surfaces as well as the gas-surface interactions characteristics, and then affect the aerodynamic performance of vehicles. Due to the cost limitations and huge distinction in space-time scales, it is difficult for current ground tests and simulation methods to duplicate the gas-surface interaction processes in the real space environment of the upper atmosphere. To solve this hard problem, an equivalent method of molecular dynamics simulation is designed in this study, which allows gas molecules to impinge on the surface at a high frequency in a short time to capture the main features of long-time but low-frequency gas-surface interaction processes in the real environment. Thus, the evolution of oxygen atoms' adsorption/corrosion state on the surface of a typical aero-space material is explored, and the momentum accommodation processes at different altitudes are analyzed. The results indicate that when certain equivalence principles are satisfied, this simulation method can reasonably reveal the real physical pictures, showing three evolution stages of the surface state, i.e., (1) adsorbing a few oxygen atoms, (2) adsorbing a large number of oxygen atoms, and (3) corroded by oxygen atoms and covered by a loose oxide. At the same time, the momentum accommodation characteristics of high-speed incident gas molecules on the corresponding surface states are also investigated. Moreover, the existence of another main component N2 in the upper atmosphere is also considered and its effects on the adsorption, corrosion and accommodation of oxygen atoms are evaluated. This study can deepen the understanding of gas-surface interactions in the upper atmosphere and provide instruction for the aerodynamic design of the upper atmosphere vehicles.
-
引 言
上层大气层(海拔100 ~ 300 km[1])飞行器的设计是中国科协遴选的2023年度十大前沿科学问题之一, 也是世界范围内航空航天领域研究的一个新热点. 上层大气层飞行器主要飞行在超低地球轨道(very low earth orbit, VLEO)上, 在军事和民用领域都有明显的应用优势, 目前美国、日本及欧洲的相关组织都在探索这一空域[2-4], 但由于上层大气层中仍然存在非常稀薄的大气, 对飞行器造成的气动阻力影响显著, 现有的飞行器, 例如超低轨卫星, 无法维持长时间在轨飞行. 在这种近自由分子流动中, 气体分子间碰撞频率极低, 高速气体分子与飞行器表面的相互作用对飞行器气动力/热的影响凸显[5]. 靳旭红等[6]在最近的综述文章中对上层大气层飞行器在气动设计问题方面面临的几项技术挑战进行了详细地介绍, 其中的气-固相互作用(gas-surface interaction, GSI)建模是影响阻力系数的关键因素. 前人的不少研究[1,7], 已经关注了GSI模型对上层大气层飞行器的气动减阻设计和机动性能的影响. 另外, 不少学者也在研究适用于此类新型飞行器的吸气式电推进系统[8-11], 其进气道设计同样依赖于对气固相互作用现象的深入认识. 这些研究大都采用了相对简单的GSI模型和假设的模型参数, 但由于GSI过程对许多影响因素(包括气体和壁面的物理化学性质、热力学状态以及气体在壁面吸附和壁面粗糙程度随时间的演化过程等)非常敏感, GSI模型参数的准确获取仍然是本领域的一大难题. 为了更精确地设计新型飞行器的气动外型和推进系统, 亟需深入理解高速气体分子, 特别是原子氧(atomic oxygen, AO), 与航天材料的相互作用机理, 确定合适的GSI模型及参数范围.
工程上常采用直接模拟蒙特卡罗(direct simulation Monte Carlo, DSMC)方法结合散射核模型来研究稀薄大气环境中飞行器的气动力/热问题. 被广泛使用的散射核模型包括Maxwell模型[12]和Cercignani-Lampis-Lord (CLL)模型[13-15]等, 能够对GSI过程中复杂的物理化学过程进行唯象简化描述, 并通过适应系数(accommodation coefficient, AC)和反应系数等参数, 综合反映气体分子与固体分子之间以及壁面势场区域内气体分子之间的相互作用, 甚至还包括了壁面微观结构的影响. 其中适应系数
$$ {{AC}} = \frac{{\overline {{\varPhi _{\text{i}}}} - \overline {{\varPhi _{\text{r}}}} }}{{\overline {{\varPhi _{\text{i}}}} - \overline {{\varPhi _{\text{s}}}} }} $$ (1) 由入射与反射的物理量通量定义, 其中$ \varPhi $表示某个物理量的通量, $ \text { i } $, $ \rm{r} $和$ \rm{s} $分别表示入射、反射以及与壁面达到热平衡时的状态. 要想应用散射核模型, 其中的适应系数等参数需要由实验或者分子动力学(molecular dynamics, MD)模拟等手段获得.
已有不少研究者[16-19]对GSI过程中适应系数的获取方式进行了讨论或综述, 本文作者[19]在前期也讨论了MD模拟方案对适应系数获取的影响. 对于自由分子流状态的来流, 通常认为气体分子间不存在相互作用, 只需多次重复模拟单个分子入射到壁面的过程并进行统计学分析, 就可以计算出适应系数的值[18]. 但是, 这种简化的模拟方案并不适用于上层大气环境. 上层大气受太阳辐射的影响, 大量氧气分子离解为AO[20]. AO的化学性质活泼, 很容易在飞行器表面吸附, 甚至侵蚀表面材料. 因而, 在此研究中不能再简单地考虑单个分子与理想光洁壁面的碰撞, 还需要考虑吸附的气体分子以及侵蚀造成的微观粗糙度对后续气体分子散射过程的影响. 另外, 文献[21-25]表明, AO的吸附/侵蚀现象会受到轨道高度与太阳活动的周期性影响, 进而导致上层大气层气固相互作用现象更加复杂化. 上述过程的时间尺度大, 远超MD模拟方法的模拟能力. 并且, 这种极端环境在地面实验中也无法完全复现, 而飞行试验也会受到成本和测试设备等多方面限制[26-29].
目前对上层大气区域GSI过程的研究, 更多是直接给定一系列不同的适应系数, 研究飞行器升阻力或吸气推进系统性能随之变化的规律. 例如Mehta等[30]研究了不同几何外形的飞行器在不同适应系数条件下的阻力系数, 发现预测结果与所采用的GSI模型以及适应系数的取值密切相关; Romano等[11]对比发现镜面反射($ {AC} = 0 $)与漫反射($ {AC} = 1 $)条件下吸气式电推进系统吸气效率存在显著差异; Rapisarda等[8]研究了适应系数对吸气式电推进系统不同气体组分的吸气效率的影响. 最近几年, 国内在该方面也开展了不少相关研究. 沈清等[1]的研究发现, 平板外形的飞行器要想获得大于1的升阻比, 适应系数需低于$ 0.2 $; 张俊等[7]提出了一种基于表面适应系数的卫星减阻设计新方法, 其计算结果表明, 卫星迎风面适应系数较大而侧面适应系数较小时, 能够实现有效减阻; 靳旭红等[31]在适应系数典型取值下对吸气式电推进系统的可行性进行了系统分析; 周靖云等[32]的研究表明, 气-固相互作用会显著影响吸气式电推进系统的进气道压缩和收集性能. 值得注意的是, 上述研究结果要进一步应用于实际飞行器设计之中, 还必须通过实验或者MD模拟获得实际问题中材料表面的适应系数. 过去研究中, 这些数据往往是从超低轨卫星飞行试验阻力数据反推出来的[23,33-34], 但这种近似估算的做法会受到飞行器结构、大气环境等因素影响, 且忽略了复杂物理化学过程的时变效应, 其应用范围值得进一步检验.
为了深入认识上层大气层GSI现象及其对气动力精确预测的影响, 本文采用MD方法研究了AO在典型航天材料TiO2上的吸附和侵蚀过程. 在前期研究[19]基础上, 本文基于AO入射数目一致和壁面充分弛豫两个原则, 提出了一种等效的MD模拟方案. 该方案通过模拟短时间内气体分子高频率入射壁面, 来反映真实空间环境中长时间低频率的分子入射情况. 研究结果展示了不同条件下材料壁面受AO吸附和侵蚀的3个特征状态, 并给出了各状态下气体分子动量适应规律的分析.
值得说明的是, MD模拟结果的可信度会受到所选取的分子作用势函数的影响, 但对于涉及复杂物理化学现象的分子作用势函数构建, 目前并没有成熟数据库可用. 特别是第一宇宙速度下GSI现象的MD模拟, 国内外相关研究尚处于起步阶段. 在此现实背景之下, 本文选取了该领域相对可靠的分子作用势模型, 并与有限的文献数据进行了比较验证, 所得结果可以认为是对上层大气层复杂GSI现象的初步探索.
1. 分子动力学模型
1.1 物理模型
本文采用开源软件LAMMPS[35]和OVITO[36]研究AO在金红石型TiO2(110) 表面的散射过程. 金红石相TiO2结构稳定且耐高温, 是航天器上常用的保护性涂料之一[29,37-38], 其 (110) 切面是表面能最低、结构最稳定的切面[39-40]. 为了确保模拟方案能够全面反映GSI中涉及的物理化学过程, 这里选择了Van Duin等[41-42]开发的反应力场(reactive force field, ReaxFF)来模拟原子之间的作用力. 反应力场曾被用于研究AO与聚酰亚胺(Kapton)、无定形SiO2等材料之间的相互作用, 并与飞行实验结果相近[43-44]. 以下所有模拟中的反应力场参数来自Kim等[45], 可以描述Ti, O和N等元素组成的材料和相关化学反应, 常被用于金红石型TiO2(110)表面的模拟研究[46-47].
为了验证反应力场数据的可靠性, 本文复现了Sokolović等[48]关于O2/O与TiO2相互作用的一个典型算例. 文献采用密度泛函理论(density functional theory, DFT)计算了尖端为O原子的探针靠近吸附了O原子的TiO2表面时, 二者之间的受力曲线, 并与原子力显微镜(atomic force microscopy, AFM)实验数据进行了比较, 如图1(a)所示, 其中${F_{{\text{ts}}}}$为探针尖端受到的作用力, ${Z_{{\text{ts}}}}$表示尖端O原子垂直壁面方向的坐标. 如图1(b)所示, 本文基于反应力场的MD模拟结果与文献中试验数据符合更好, 可用于本研究中对AO与TiO2表面之间的动量传递的模拟. 另外, Tian等[49] 采用DFT计算N2在TiO2表面的吸附能为0.15 eV, 而Lee等[50]的DFT计算结果则为0.04 eV, 其中的差异和DFT计算中的设置不同有关. 本文计算了如图1(c)所示N2在TiO2表面的吸附能为0.21 eV, 与Tian等[49]的结果较为接近. 总的来看, 关于O和N在TiO2表面吸附的研究很少见, 且不同文献结果之间存在一定散度. 作为初步探索, 本文采用的MD势函数模型足以获得AO和N2两种气体与TiO2之间相互作用的一些规律. 当然, 未来要想得到更精确的结果, 则需要针对AO, N2与TiO2的各种相互作用情形进行精确的量子化学计算并对反应力场的参数进行有针对性的修正.
![]()
图 1 (a)本文MD模拟原子力探针受力结果与文献[48]中AFM测量和DFT结果的对比; (b)本文MD模拟原子力探针算例截图; (c)本文MD模拟N2在TiO2表面吸附的截图Figure 1. (a) Comparison of present MD simulation of AFM tip force with measurement data and DFT result in Ref. [48]; (b) Screenshot of present MD simulation of AFM tip; (c) Screenshot of present MD simulation of N2 adsorption on TiO2 surface本文作者最近在一个关于GSI过程的研究中[19], 详细对比讨论了单个入射(single scattering, SS)与连续入射(continual scattering, CS)的MD模拟方案. SS方案重复单个气体分子入射到壁面的过程并进行统计学分析, 这种方案不存在气体分子之间的相互作用及宏观流动; CS方案通过控制入射频率模拟气体稀薄程度的变化, 让气体分子连续入射到壁面, 这种方案包含了壁面吸附气体分子以及势场范围内气体分子之间相互作用的影响. 研究发现, 即使在高稀薄环境中, 壁面吸附的气体分子也会对后续入射分子GSI过程产生影响.
根据常用的大气模型NRLMSIS 2.0[51], 如图2所示, 上层大气层气体分子数密度为${10^{13}}\sim {10^{18}}{{\text{ m}}^{ - 3}}$, 其中AO数密度${n_{{\text{AO}}}}$ = ${10^{13}}\sim {10^{17}}{{\text{ m}}^{ - 3}}$, 对于图3中所示面积${A_{{\text{sur}}}} = 1495.53$ Å2的壁面, 连续两个AO入射到壁面的时间间隔为$ 10^{4} \sim 10^{8} \mathrm{ ~ ns} $, 远长于一般分子动力学模拟的规模$(\sim 100{\text{ ns}})$[52]. 因此, 想要通过CS方案直接模拟上层大气环境是不现实的, 需要结合物理过程的特征分析提出一种等效模拟的新方案. 考虑到连续两次入射的间隔时间尺度极大, 可以认为每个时刻的壁面状态只与此前已经入射和吸附的AO数目有关, 因此, 可以在模拟中适当提高AO入射的频率. 但是, AO连续入射的频率也不能过高, 否则壁面分子状态没有足够时间弛豫, 模拟结果将明显偏离真实情况, 这一点在1.2小节模拟设置中有详细说明. 基于以上讨论, 本文确立了等效模拟方案的两条原则如下.
![]()
![]()
图 3 AO与金红石型TiO2(110)表面相互作用的模拟过程示意图Figure 3. Schematic diagram illustrating the simulation process of AO interacting with the rutile TiO2(110) surface(1) AO入射数目一致原则. 根据AO入射数目一致, 壁面在AO来流中的暴露时长$t$与AO数密度$n$之间需满足$nt = {\mathrm{const}}.$ 由此可以建立高频AO入射模拟中的物理时间与真实飞行时间之间的对应关系, 进而得到随时间演化的模拟壁面状态与暴露在不同AO密度环境中的真实壁面状态之间的映射关系.
(2)壁面充分弛豫原则. 即模拟中的连续两次AO入射的时间间隔$\Delta t$原则上应大于壁面弛豫特征时间$\tau $. 这就要求AO入射频率不能过高, 以使壁面充分弛豫, 达到较为稳定的状态.
基于以上两条原则, 本文把前文所述的CS和SS模拟方案进行结合, 首先通过CS方案快速获得一系列不同吸附和侵蚀程度的壁面状态, 再选择某密度和某时刻对应的壁面状态, 进行SS方案的AO入射模拟, 整个流程如图3所示, 其中来流入射的AO被显示为深红色, 以区别于壁面材料中原有的O原子.
另外, 从图2可以注意到, 海拔$200{\text{ km}}$以下, N2才是大气的主要成分. 即使N2不如AO活性高, 不易吸附和反应, 但N2的存在是否会对AO在壁面的吸附和侵蚀产生影响也是值得探究的问题. 在本文第3部分, 也对海拔$200{\text{ km}}$以下N2/AO混合来流条件的GSI问题开展了讨论.
1.2 模拟设置
上层大气飞行器的飞行速度一般在第一宇宙速度附近, 因此, 所有模拟中的来流气体分子速度被设为$ {v_{\text{g}}} = 7800{\text{ m/s}} $. 在低轨卫星等航天器设计中, 其内部温度一般维持在$ 300{\text{ K}} $左右[21], 以保证精密元器件的正常运行. 为了尽可能反映实际情况, 所有模拟中的壁面进行了如图4所示的分层处理. 固定层是为了确保在整个模拟中壁面不会由于受到气体分子的碰撞而发生位移. 温控层通过Langevin热浴将温度稳定在$ {T_{\text{s}}} = 300{\text{ K}} $. 过渡层通过虚拟弹簧力连接固定层与温控层, 这是文献中常见的处理方法[53]. 直接作用层(active layer)则采用微正则系综(NVE), 只受到原子之间势函数作用, 能够更真实地反映壁面原子与来流气体分子之间的动力学过程.
图4中的蓝色虚线以下代表壁面附近的势场作用范围. 对于每个入射的气体分子, 记录其进入与离开势场范围时的速度$ {{\boldsymbol{v}}_{\text{i}}} = \left( {v_{{x}}^{\text{i}},v_{{y}}^{\text{i}},v_{{z}}^{\text{i}}} \right) $与$ {{\boldsymbol{v}}_{\text{r}}} = \left( {v_{{x}}^{\text{r}},v_{{y}}^{\text{r}},v_{{z}}^{\text{r}}} \right) $. 若有气体分子在势场范围内滞留50 ps以上, 认为其被壁面捕获并达到完全适应, 对应适应系数为1.0. 对大量相同模拟条件下的分子散射算例进行统计学处理, 可以获得这一模拟条件下的适应系数. 在入射平面内计算入射与反射分子的切向速度分量
$$ v_{\text{t}}^{\text{i}} = \sqrt {v{{_{{x}}^{\text{i}}}^2} + v{{_{{y}}^{\text{i}}}^2}} $$ (2) $$ v_{\text{t}}^{\text{r}} = \frac{{v_{{x}}^{\text{i}}v_{{x}}^{\text{r}} + v_{{y}}^{\text{i}}v_{{y}}^{\text{r}}}}{{\sqrt {v{{_{{x}}^{\text{i}}}^2} + v{{_{{y}}^{\text{i}}}^2}} }} $$ (3) 并假设有$ \xi $比例的来流分子被壁面捕获, 就可以计算切向动量适应系数(tangential momentum accommodation coefficient, TMAC)[19,54]
$$ {{TMAC}} = \frac{{\overline {mv_{{t}}^{\text{i}}} - \overline {mv_{{t}}^{\text{r}}} }}{{\overline {mv_{{t}}^{\text{i}}} }} \cdot (1 - \xi ) + \xi $$ (4) 法向动量适应系数(normal momentum accommodation coefficient, NMAC)在文献中有多种计算方式[16], 考虑到工程中对阻力系数的估计通常选择适应系数在$ 0\sim 1 $的范围, 一般采用下式计算
$$ {{NMAC}} = \frac{{\overline {\left| {mv_{{z}}^{\text{i}}} \right|} - \overline {\left| {mv_{{z}}^{\text{r}}} \right|} }}{{\overline {\left| {mv_{{z}}^{\text{i}}} \right|} - \sqrt {\text{π} m{k_{\text{B}}}{T_{\text{s}}}/2} }} \cdot (1 - \xi ) + \xi $$ (5) 在高频连续入射的等效方案中, 入射频率是非常重要的参数. 为了与上层大气环境的气体分子数密度进行直观的对比, 下面将入射频率换算成等效的来流数密度进行衡量. 理想情况下, 单位面积与单位时间内气体分子在壁面的碰撞数为$ n{v_{\text{g}}} $, 其中$ n $代表数密度. CS模拟方案中, 如果面积$ {A_{{\text{sur}}}} $壁面上每隔$ \Delta t $时间受到一次分子入射, 则其等效的来流数密度为
$$ {n_{{\text{CS}}}} = \frac{1}{{{v_{\text{g}}} \cdot {A_{{\text{sur}}}} \cdot \Delta t}} $$ (6) 连续的高能AO入射会对壁面局部造成剧烈的加热效应, 而温度急剧上升势必会导致壁面更容易被侵蚀和破坏. 这一现象在不少文献[28,55-57]中都有指出, 这些研究中的等效数密度在$ {10^{24}}\sim 10^{27}{{\text{ m}}^{ - 3}} $范围内, 远比上层大气环境稠密. 但由于上层大气非常稀薄, 每个气体分子入射后, 壁面都有足够长的弛豫时间, 实际飞行器表面温度不会显著升高. 实际环境中的壁面状态弛豫时间与许多因素有关, 在MD模拟中, 还需要考虑模拟温控设置的影响, 很难给出一个明确的弛豫时间结果. 因此, 本文直接监控图4所示的直接作用层的温度变化, 来检验模拟中时间步长和气体分子入射频率设置的合理性. 图5显示了在不同模拟时间步长及等效气体分子数密度下, 壁面温度随着入射分子数增加而升高的程度. 对于任一时间步长, 等效来流数密度越低, 壁面温度都越接近$ {T_{\text{s}}} = 300{\text{ K}} $. 当CS方案的等效来流数密度在$ {n_{{\text{CS}}}} = 9 .0\times {10^{23}}{{\text{ m}}^{ - 3}} $以下时, 壁面温度在AO大量入射之后仅有稍许升高.
![]()
图 5 不同时间步长以及等效分子数密度下, CS方案对壁面的加热效应随入射分子数目的变化Figure 5. The heating effects on the surface in CS approach vary with the incident AO number under different time step settings and equivalent number densities在不涉及化学反应的高速气体GSI研究中[19], 壁面温度的稍许变化对适应过程的影响不大. 当考虑化学反应时, 温度升高会促进化学反应的发生, AO越容易被壁面捕获. 为了精准捕捉化学键的断裂与生成, 通常将MD模拟时间步长设为$ {\text{0}}{\text{.25 fs}} $. 然而上层大气环境中, 壁面侵蚀状态显著变化的时间尺度却很长, 即使图5中的等效分子数密度都远高于上层大气的真实密度, MD模拟的时间成本依然很高. 事实上, 采用相同的时间步长, 并模拟相同的分子入射数目, $ {n_{{\text{CS}}}} $ = $ 9.0 \times {10^{22}}{{\text{ m}}^{ - 3}} $模拟时间比$ 9.0 \times {10^{23}}{{\text{ m}}^{ - 3}} $情况高了两个数量级. 在等效分子数密度和相同分子入射数目条件下, 时间步长$ {\text{1}}{\text{.0 fs}} $的模拟成本是$ {\text{0}}{\text{.25 fs}} $情况的4倍. 严格来说, 通过某时刻AO所具有的平动能是否大于其所处位置的势能可以判断此时的AO是否处于被壁面吸附的状态[58-59], 进而可以统计在壁面势场范围内吸附的AO数目并计算其面密度${\sigma _{{\text{AO}}}}$. 不过在实际操作中发现, 直接统计滞留在壁面势场范围内的AO数目基本可以代表吸附AO的数目, 这也是一些文献中计算来流分子吸附时长或判断吸附状态的简化方法[60-61]. 图6显示了CS方案设置不同时间步长以及等效分子数密度条件下, 随着入射分子数目增加, 壁面上吸附分子的面密度${\sigma _{{\text{AO}}}}$的变化趋势. 对比图6中相同$ {n_{{\text{CS}}}} $条件下${\sigma _{{\text{AO}}}}$的数据, 并结合图5中各$ {n_{{\text{CS}}}} $对应的温度差异, 可以发现时间步长对AO吸附程度的影响远远小于壁温变化的影响. 因此, 综合考虑计算精度和效率之后, 本文选择了$ {n_{{\text{CS}}}} = 9.0 \times {10^{23}}{{\text{ m}}^{ - 3}} $(每两次分子入射的时间间隔为$ \Delta t = 9545{\text{ fs}} $)与时间步长$ {\text{1}}{\text{.0 fs}} $进行初步研究.
![]()
图 6 不同时间步长以及等效分子数密度下, 吸附AO的面密度${\sigma _{{\text{AO}}}}$随入射分子数目的变化Figure 6. The surface densities ${\sigma _{{\text{AO}}}}$ vary with the incident AO number under different time step settings and equivalent number densities2. 结果与讨论
2.1 高速AO在壁面的吸附/侵蚀状态
金红石型的TiO2中每个Ti原子周围有6个配位O原子, 每个O原子周围有3个配位Ti原子. 理想的TiO2(110)表面O原子的配位环境不变, 而Ti原子则一半是五配位的Ti(5), 一半是四配位的Ti(4)[40]. 在极为特殊的条件下, 外界的O原子可以诱使TiO2(110)表面重构成更稳定的表面TiO2(110) (1 × 1)[62]. 而在本文的模拟中, 来流的AO也会按照一定的排布吸附在壁面Ti原子配位不饱和的地方. 图3中所示的壁面最表层有156个Ti原子和156个O原子. 若所有配位缺失位点均吸附一个AO, 则至多可能吸附234个AO. 不过实际上, AO倾向于按照近似于TiO2(110) (1 × 1)的结构吸附, 如图7所示, Ti(5)会再吸附一个AO以形成TiO2晶体中的八面体结构, Ti(4)则会两两吸附一个AO, 总共吸附156个AO. 据此可以估计, 吸附的AO完整铺满壁面一层时的面密度$ {\sigma _{\text{0}}} $约$ 1.04\times {10^{15}} \sim 1.56 \times {10^{15}}{\text{ c}}{{\text{m}}^{ - 2}} $. 考虑到实际吸附状态没有那么完美, 而且除了填补配位缺失, AO还可能与壁面原有的O原子发生置换, 因此上面预估的AO吸附量仅作参考.
随着模拟进程发展, AO连续不断入射到壁面, 从模拟结果可发现AO吸附状态经历了3个阶段: 从少量吸附到大量吸附, 最后发展为松散的氧化层. 整个演化过程如图8从右向左所示. 同时, 根据AO入射数目一致的原则, 图8中图像从左至右也等效于相同壁面在不同海拔AO环境中暴露$ {\text{15 s}} $后的吸附/侵蚀情况. 从洁净壁面到阶段Ⅰ, 入射AO很容易和壁面原子直接接触并被吸附. 当壁面大半区域被AO覆盖, 壁面达到阶段Ⅱ, 入射AO与壁面原子的直接接触受到吸附AO的阻碍, 反而主要与吸附的AO发生碰撞, 甚至发生复合反应形成O2. 当壁面被过量AO覆盖后, 壁面状态转入阶段Ⅲ, 逐渐增多的AO导致表层Ti原子发生位移, 一同混合形成松散的氧化层, 致使晶体结构受到破坏. 若采用前文估计的参考面密度$ {\sigma _{\text{0}}} $作进一步定量分析, 阶段Ⅰ转化为阶段Ⅱ时对应的临界面密度近似为$ {\sigma _{\text{0}}}/2 $, 阶段Ⅱ转化为阶段Ⅲ时的临界面密度则近似为$ 2{\sigma _{\text{0}}} $. 图6显示, 这样通过物理分析得到的壁面状态演化过程与实际模拟出的${\sigma _{{\text{AO}}}}$曲线变化情况相符.
![]()
图 8 壁面吸附AO的状态演化过程以及与海拔的对应关系Figure 8. The evolution process of AO adsorption on surface and its corresponding relationship with altitude将3个阶段的吸附速率分别提取出来, 可以近似预测给定数密度AO来流下的壁面在某时刻的吸附速率. 由于阶段Ⅰ表面活性位点密集排布, 可以近似认为入射AO都受到化学吸附. 类似现象在He等[63]关于SiO2材料催化特性的MD模拟也有显示, 在模拟的初始阶段, 晶体表面密集的活性位点将导致高能AO几乎全部被壁面吸附. 假定入射AO全部被壁面捕获, 理论预测出的吸附面密度曲线如图9中红色虚线所示, 虚线斜率即为吸附速率, 由红色文字标注. 理论预测与MD模拟初始阶段的结果(数据点)完全吻合, 表明了上述物理分析的合理性. 对于阶段Ⅱ与阶段Ⅲ, 可根据模拟结果拟合出等效数密度$ {n_{{\text{CS}}}} $条件下的吸附速率, 在图9中分别由蓝色与绿色虚线和文字标出. 用$ \eta _{{\text{AO}}}^{{\text{ori}}} $表示模拟出的吸附速率, 若实际飞行中来流AO数密度为$ {n_{{\text{real}}}} $, 则某时刻的吸附速率$ \eta _{{\text{AO}}}^{{\text{real}}} $为
![]()
图 9 吸附过程3个阶段中的AO面密度及吸附速率Figure 9. The surface densities of AO and the corresponding adsorption rates in three stages$$ \eta _{{\text{AO}}}^{{\text{real}}} = \eta _{{\text{AO}}}^{{\text{ori}}} \cdot \frac{{{n_{{\text{real}}}}}}{{{n_{{\text{CS}}}}}} $$ (7) 其中, $\eta _{{\text{AO}}}^{{\text{ori}}}$需要根据1.1小节提出的AO入射数目一致原则来确定. 据此, 一次CS方案的模拟可以换算出一系列不同海拔不同时刻的壁面状态, 大大降低计算成本. 真实飞行条件的AO数密度越低, 则能够等效的真实飞行时间越长.
表1显示了壁面暴露在不同海拔的AO数密度中一些典型时刻的吸附速率, 其值基本在$ {10^{13}}\sim {10^{14}}{\text{ c}}{{\text{m}}^{ - 2}} \cdot {\text{mi}}{{\text{n}}^{ - 1}} $范围内, 可以与Wood[64]和Moe等[65]分别从地面与飞行实验得到的数据进行比较. Wood等[64]分别在金属Au和Ag表面做地面实验, 金属Ag的AO吸附速率保持为$ 3.3 \times {10^{14}}{\text{ c}}{{\text{m}}^{ - 2}} \cdot {\text{mi}}{{\text{n}}^{ - 1}} $, 而金属Au表面的AO吸附速率分为两个阶段, 在入射$ {\text{5 min}} $之内为$ 1.7 \times {10^{14}}{\text{ c}}{{\text{m}}^{ - 2}} \cdot {\text{mi}}{{\text{n}}^{ - 1}} $, 在$ {\text{5 min}} $之后壁面吸附状态趋于稳定. 同样是在金属Au表面, Moe等[65]的飞行实验受到太阳活跃周期的影响, 吸附速率也随之变化, 在近地点(海拔$ {\text{420 km}} $)的吸附速率为$ 7.7 \times {10^{11}}{\text{ c}}{{\text{m}}^{ - 2}} \cdot {\text{mi}}{{\text{n}}^{ - 1}} $, 更趋近于Wood地面实验$ {\text{5 min}} $后的结果. 这是因为他们飞行器的在轨运行时间长, 金属Au表面已被AO大面积覆盖, 进一步吸附AO的能力剧烈下降. 这一现象对应于阶段Ⅲ的壁面状态, 不同的是金属Au很难被氧化, 因此吸附速率会变得极小, 而金属Ag的吸附速率与TiO2更为相近, 二者均容易被氧化, 并且被AO覆盖后会继续形成氧化物. 3种材料吸附能力的差异, 定性反映了其抗氧化性能高低, 与文献[29]中根据体积损失率推算出的结果相符.
表 1 不同海拔H不同时刻AO在壁面的吸附速率$ {\boldsymbol{\eta}} _{{\bf{AO}}}^{{\bf{real}}} $ (单位: ${{\boldsymbol{10}}^{{\boldsymbol{14}}}}{\bf{\;c}}{{\bf{m}}^{ {\boldsymbol{- 2}}}} {\boldsymbol{\cdot}} {{\bf{min}} ^{ {\boldsymbol{- 1}}}}) $Table 1. AO adsorption rates $\eta _{{\text{AO}}}^{{\text{real}}}$ on the surface at different altitudes H and times $({\text{unit: }}{10^{14}}{\text{ c}}{{\text{m}}^{ - 2}} \cdot {\min ^{ - 1}}) $H/km ${n_{{\text{AO}}}}$/m−3 15 s 1 min 2 min 3 min 4 min 200 34.9 33.3 — — — — 250 12.9 46.1 — — — — 300 5.2 18.6 5.0 5.0 5.0 — 350 2.2 7.8 7.8 2.1 2.1 2.1 400 0.93 3.3 3.3 3.3 3.3 0.89 本文MD模拟中的阶段Ⅰ时间极短, 它对应于理想光洁表面的初始吸附阶段, 在实验环境中很难复现. 考虑到MD模拟与实验在时空尺度上的巨大差异, 以及真实壁面材料和来流环境的复杂性, 可以认为本文研究已经捕捉到了上层大气层AO吸附/侵蚀材料的主要物理过程, 并推算出了定量合理的吸附速率. 并且, 本文在模拟中也观察到极少Ti原子以TiOx化合物的形式脱离壁面. 因而本文提出的方法, 原则上可进一步获得材料体积损失率[29]等重要数据, 但这需要在未来开展更长时间和更大规模的模拟.
2.2 不同吸附/侵蚀状态壁面上AO的动量适应
如图10所示, 观察AO在壁面附近的运动轨迹, 它除了被壁面捕获之外, 还可能从壁面直接解吸离开, 或者复合反应形成O2离开. 若气体分子被壁面捕获, 其动量会与壁面达到热平衡, 可以认为动量完全适应, 这与式(4)和式(5)所示适应系数常规的计算方式相同; 若来流AO在壁面形成了O2才离开, 则需要计算整个O2向壁面传递的动量, 以准确计算飞行器表面所受的气动力. 根据适应系数的定义式(1), 如果$\xi_{1}$和$ \xi_{2} $分别表示来流AO被壁面捕获以及形成O2的占比, 则需要对适应系数的计算进行如下修正
![]()
图 10 AO入射壁面可能的运动轨迹Figure 10. The possible trajectories of AO interacting with the surface$$ \begin{split} & {{TMA}}{{{C}}_{{\text{modified}}}} = \frac{{\overline {mv_{{t}}^{\text{i}}} - \overline {mv_{{t}}^{\text{r}}} }}{{\overline {mv_{{t}}^{\text{i}}} }} \cdot \left( {1 - {\xi _1} - {\xi _2}} \right) +\\ &\qquad \frac{{\overline {mv_{{t}}^{\text{i}}} - \overline {2mv_{{t}}^{\text{r}}} }}{{\overline {mv_{{t}}^{\text{i}}} }} \cdot {\xi _2} + {\xi _1} \end{split}$$ (8) $$ {{NMA}}{{{C}}_{{\text{modified}}}} = \frac{{\overline {\left| {mv_{{z}}^{\text{i}}} \right|} - \varPhi _{{n}}^{\text{r}}}}{{\overline {\left| {mv_{{z}}^{\text{i}}} \right|} - \sqrt {\text{π} m{k_{\text{B}}}{T_{\text{s}}}/2} }} $$ (9) 其中
$$ \varPhi _{{n}}^{\text{r}} = \overline {\left| {mv_{{z}}^{\text{r}}} \right|} \cdot \left( {1 - {\xi _1} - {\xi _2}} \right) + \overline {\left| {2mv_{{z}}^{\text{r}}} \right|} \cdot {\xi _2} + {\xi _1} \cdot \sqrt {\text{π} m{k_{\text{B}}}{T_{\text{s}}}/2} $$ 图8显示了从CS方案获得的一系列暴露在不同海拔AO环境$ 15\;{\mathrm{s}} $后的壁面吸附/侵蚀状态. 海拔较高(来流AO数密度较低)时, 少量AO吸附在壁面, 壁面状态处于2.1小节中所述的阶段Ⅰ, 入射AO容易被壁面捕获或与壁面吸附的AO形成O2. 大部分吸附AO沿着壁面晶格规律排布, 形成一道道“棱”, 使表面变得不平整. 在海拔$ 300 \sim 350 \mathrm{ ~ km} $, 将出现两种现象: 一方面, 壁面上大部分区域被AO覆盖, 减少了壁面可以吸附AO的活性位点, 增大了入射AO与吸附AO形成O2的可能性; 另一方面, 大量的吸附填补了表面的“棱”间凹陷, 使得壁面的“粗糙度”下降, 并阻碍了入射AO与壁面的直接作用. 在海拔$ { 200 \sim 250\;{\mathrm{km}} } $, 壁面上吸附的氧原子堆叠起来, 进一步的侵蚀使得表面的Ti原子显著位移, 甚至有极少量的TiOx化合物离开壁面. 堆叠的氧化层结构松散, 有明显的孔洞结构, 粗糙程度高, 会极大地促进入射AO被壁面捕获的可能.
为了佐证上述观点, 针对图8所示的6个壁面, 用SS方案模拟AO以0°及45°角度入射, 每组算例重复3000次以上并统计相应壁面条件下的TMAC, NMAC以及各类运动轨迹的比例, 结果见图11. 采用式(4)和式(5)常规方式计算的TMAC与NMAC排除了形成O2的算例, 而采用式(8)和式(9)计算的TMACmodified与NMACmodified则考虑O2对动量计算的贡献. 对比常规与修正的适应系数, 可以直观看出O2散射在适应过程中所起的作用. 由图12可见, 洁净壁面极易捕获AO, 在前文中已经说明这与壁面上密集排布的活性位点有关. 虽然现实中很难保持壁面完全洁净, 但这种理想状态可以作为一种极限参考. 如前文所述, 在海拔较高时, 大部分入射AO会被壁面捕获, 充分地向壁面传递动量; 在海拔较低时, 结构松散、凹凸不平的表面也容易吸附AO, 因而接近完全适应. 在海拔$ 300 \sim 350 \mathrm{ ~ km} $时, 前文提到的两种现象对适应过程分别起到促进和阻碍的作用, 二者相互抵消之后的综合表现为阻碍作用, 使TMACmodified与NMACmodified稍有降低, 但如果不考虑O2散射的影响, 则只有上述的阻碍作用存在, 适应系数TMAC与NMAC显著降低. 因而, 当有较大比例的入射AO形成O2时, 其对动量适应的影响不可忽略.
![]()
图 11 不同吸附/侵蚀状态下壁面的AO动量适应系数(AC)Figure 11. Momentum ACs of AO on surfaces at different adsorption/corrosion states![]()
图 12 不同吸附/侵蚀状态下壁面上AO各类运动轨迹占比Figure 12. The proportion of various trajectories of AO at different adsorption/corrosion states从图12也可以发现, 入射角为45°时, AO在洁净壁面更容易被吸附, 在阶段Ⅰ与阶段Ⅱ的壁面更容易形成O2. 这是因为法向动量分量较低时, 入射AO在壁面势场范围运动时间更长, 因而更容易被吸附, 而切向运动则促使入射AO与壁面上更多的吸附AO相互作用, 二者均对适应过程有促进作用. 但入射AO受到壁面吸附时, 在水平方向上依然有一定动量, 需要一段时间达到完全适应, 而处于阶段Ⅱ的壁面上吸附的AO可能阻碍入射AO与壁面直接接触, 其扰动效应也可能导致后续入射AO更容易解吸附. 因而, 切向运动的存在同时也对处于阶段Ⅱ的壁面上的适应过程表现出一定的阻碍作用.
若壁面吸附情况处于阶段Ⅲ, 则粗糙度成为影响适应的主要因素. 阶段Ⅰ和阶段Ⅱ表面的“粗糙度”由吸附的AO形成, 基本不高于一层原子, 但阶段Ⅲ的粗糙度凹陷与空洞结构尺度较大, 使气体分子与壁面相互作用的有效面积增加. 并且, 粗糙表面由Ti和O共同构成, 相较于AO, 较重的Ti原子更难离开壁面. 0°入射的AO有足够的能量打入松散的氧化层并被困于孔洞中, 但45°入射的AO法向动量较低, 分子更倾向于吸附在表面. 上文所述切向运动对形成O2的促进效果被表面粗糙结构弱化了, 但阻碍作用却依然存在. 也就是说, 尽管粗糙壁面本身为入射AO提供了更多相互作用面积而促进适应过程, 但凹凸不平的壁面对具有切向运动的吸附AO也有扰动, 使AO在完全适应之前就已经解吸附. 因此, 对比0°和45°的NMACmodified可以发现, 入射角45°时, AO在阶段一壁面上被吸附和形成O2的比例高, 因而动量适应程度更高; 处于阶段Ⅲ的壁面会促使切向运动的AO解吸附, 综合表现为阻碍作用, 其动量适应程度低于0°入射角的情况; 处于阶段Ⅱ的壁面, 入射角变化对适应过程的促进和阻碍作用相当, 综合影响不显著.
2.3 N2/AO混合气体在壁面吸附和侵蚀的特征
海拔低于200 km时, 大气占比最高的成分为N2, 高速的N2入射到壁面上也可能发生吸附和侵蚀现象, 甚至一部分N2会发生离解. 不过总的来说, N2的化学性质没有AO活泼, 因此在壁面附近吸附和反应的N2占比很小. 图13显示了不同海拔不同混合比例(以粒子数之比${N_{{{\text{N}}_{\text{2}}}}}:{N_{{\text{AO}}}}$表征)的N2/AO混合物连续入射壁面的吸附和侵蚀情况, 其中每个阶段不同海拔处的情况, 保持入射AO数目相同, N2是在此基础上附加的, 以对比显示N2组分存在与否的影响. 图14则显示了前文所述AO吸附3个阶段中, 混合比例对壁面吸附情况的影响. 少量的N2离解为N原子, 并占据表面的活性位点, 使得AO吸附的面密度${\sigma _{{\text{AO}}}}$相应地下降. 当N2和AO混合比接近10时, N原子与AO在壁面吸附的面密度达到一致. 阶段Ⅰ与阶段Ⅱ的壁面吸附状态与吸附位点相关, 而离解的N原子比AO的配位能力强, 因此壁面总的吸附粒子面密度${\sigma _{{\text{total}}}}$会随着入射N2数目增多而先降低后增加. 阶段Ⅲ的壁面粗糙程度高, N原子的配位能力对AO吸附干扰小, N原子与AO均容易被壁面吸附, 因此${\sigma _{{\text{total}}}}$随着入射N2数目增多而显著增加. 总的来说, 海拔200 km以下的N2/AO混合入射对${\sigma _{{\text{total}}}}$的影响达不到量级变化的程度. 在海拔200 km及以上, N2组分比例很低, 可以合理推断, 极少量的N2掺混对壁面总的吸附状态影响不大. 但是, 由于N原子一旦吸附, 就比AO更难脱离壁面, 其对后续入射AO的动量适应过程可能会产生一定影响, 详见下文讨论.
![]()
图 13 不同海拔和不同N2/AO混合比气体在壁面的吸附/侵蚀Figure 13. Adsorption/corrosion features of surfaces exposed in different N2/AO gas mixtures at different altitudes![]()
图 14 不同海拔和不同N2/AO混合比气体环境中的壁面吸附面密度Figure 14. Surface densities of adsorbed atoms in different N2/AO gas mixtures at different altitudes由于高速的N2在壁面上发生吸附和侵蚀的概率极小, 很难充分传递其携带的动量, 因此N2/AO混合气体会比纯AO在壁面的动量适应程度低. 以图13(a) 处于阶段一的4个典型壁面为例, 可以定量研究N2/AO混合气体的适应规律. N2与AO均以$ {v_{\text{g}}} = 7800{\text{ m/s}} $的速度垂直入射壁面, 壁面温度仍然设为$ {T_{\text{s}}} = 300{\text{ K}} $, N2的转动角速度根据海拔$ {\text{100}} $, $ {\text{150}} $和$ {\text{200 km}} $处的温度($ {\text{187}}{\text{.7}} $, $ {\text{716}}{\text{.0}} $和$ {\text{920}}{\text{.8 K}} $)进行设定. 除了前文讨论过的纯AO入射壁面上的3种轨迹类型, 此时AO还可能与壁面吸附的N原子形成氮氧化合物(最主要是NO). 另外, 除了在壁面发生吸附-解吸或被壁面捕获, N2还可能在与壁面相互作用过程中离解, 而N原子在壁面上可能的轨迹包括散射离开壁面、被壁面捕获、与壁面的O原子或者吸附AO形成氮氧化合物、与吸附N原子形成N2等. 综上, N2/AO混合气体入射的气固相互作用过程非常复杂, 因此这里对其中占比极小(见图15)的散射轨迹进行了简化处理. 如果N原子与其他原子结合形成了任何一种分子并离开壁面, 则认为N原子以所形成的分子的速度散射, 忽略与其结合的其他原子的贡献. 因此, 法向动量适应系数被修正为
![]()
图 15 N2/AO混合气体入射和纯AO入射条件下粒子各种轨迹的占比Figure 15. Proportions of various particle trajectories on surfaces exposed in N2/AO gas mixtures and pure AO gas$$ {{NMA}}{{{C}}_{{{\text{N}}_{\text{2}}}{\text{+AO}}}} = \frac{{\displaystyle\sum\limits_{{j} = 1}^{{N_{{\text{total}}}}} {{{\left| {mv_{{z}}^{\text{i}}} \right|}_{j}}} - \varPhi _{{n}}^{\text{r}}}}{{\displaystyle\sum\limits_{{j} = 1}^{{N_{{\text{total}}}}} {{{\left| {mv_{{z}}^{\text{i}}} \right|}_{j}}} - \displaystyle\sum\limits_{{j} = 1}^{{N_{{\text{total}}}}} {\sqrt {\text{π} {m_{j}}{k_{\text{B}}}{T_{\text{s}}}/2} } }} $$ (10) 其中
$$ \begin{split} & \varPhi _{{n}}^{\text{r}} = \sum\limits_{{j} = 1}^{{N_{{\text{desorbed}}}}} {{{\left| {mv_{{z}}^{\text{r}}} \right|}_{{j}}}} + \sum\limits_{{j} = 1}^{{N_{{\text{NO}}}}} {{{\left| {\left( {{m_{\text{N}}} + {m_{\text{O}}}} \right)v_{{z}}^{\text{r}}} \right|}_{{j}}}} + \sum\limits_{{j} = 1}^{{N_{{{\text{O}}_{\text{2}}}}}} {{{\left| {2{m_{\text{O}}}v_{{z}}^{\text{r}}} \right|}_{{j}}}} + \\ &\qquad \sum\limits_{{j} = 1}^{{N_{{\text{2N}}}}} {{{\left| {{m_{\text{N}}}\left( {v_{{z}}^{{\text{r,1}}} + v_{{z}}^{{\text{r,2}}}} \right)} \right|}_{{j}}}} + \sum\limits_{{j} = 1}^{{N_{{\text{trapped}}}}} {\sqrt {\text{π} {m_{{j}}}{k_{\text{B}}}{T_{\text{s}}}/2} } \end{split}$$ 式中, $ m $, ${m_{\text{O}}}$和${m_{\text{N}}}$分别表示具体某种原子或分子、O原子和N原子的质量. Ntotal, Ndesorbed和Ntrapped分别表示总算例数、入射粒子被壁面吸附后又解吸和被壁面捕获的算例数; $ {N_{{\text{NO}}}} $和${N_{{{\text{O}}_2}}}$分别表示AO在壁面附近形成NO和O2的算例数; $ {N_{2{\text{N}}}} $表示N2离解且两个N原子(由上标1, 2区分)轨迹不同的算例数. 由于忽略了与N原子形成分子的其他粒子所携带的动量, 上式计算的适应系数偏大. N原子能够带走的其他粒子质量${m_{{\text{err}}}}$在$0\sim 2{m_{\text{O}}}$之间, 散射速度为$ v_{{z}}^{\text{r}} $, 由此可以估计出式(10)的计算方法造成的误差, 也显示在图16中.
![]()
图 16 不同海拔的N2/AO混合气体在吸附阶段一壁面上的NMAC结果Figure 16. The NMAC results of N2/AO gas mixtures impinging on surfaces at the first adsorption stage at different altitudes图16显示了不同海拔高度, 图13(a)中处于阶段一的4个壁面上的NMAC. 显然, 由于高速N2难以与壁面完全适应, 所以海拔越低, N2占比越高, 混合气体等效的${{NMA}}{{{C}}_{{{\text{N}}_{\text{2}}}{\text{+AO}}}}$越小. 这一结果表明, 对气动力进行工程计算时, 应该针对来流中N2的占比选择合适的等效${{NMA}}{{{C}}_{{{\text{N}}_{\text{2}}}{\text{+AO}}}}$值, 或者将不同组分分开考虑.
若仅统计其中AO入射对应的NMACmodified, 可以发现NMACmodified受N2占比高低的影响并不大, 但与纯AO环境计算所得的NMACmodified仍有稍许差异. 结合图15(b)中AO的吸附-解吸轨迹占比变化可知, 如果壁面上有N原子占据部分位点, 将导致AO被壁面捕获或者形成分子的过程受到一定程度阻碍, 因而更多的AO解吸并离开壁面, 使得法向动量适应程度偏低. 因此, 即使壁面吸附粒子总量处于同一水平, 但极少量的N原子吸附在壁面, 也可能会对后续入射AO的动量适应产生影响. 需要注意的是, 由于计算规模限制, 图16显示的是尚处于气体吸附初始阶段的壁面上的动量适应特征, 能够反映在大气环境中暴露很短时间的情况. 但若要反映飞行器长时间运行后达到吸附/侵蚀饱和的壁面特性, 未来还需要进一步研究.
3. 结 论
本文采用物理建模和分子动力学方法, 对上层大气层气固相互作用进行了初步分析和模拟, 得到以下结论.
(1)鉴于上层大气极度稀薄, 并且其主要成分之一AO化学性质活泼, 本文基于AO入射数目一致原则与壁面充分弛豫原则提出了一种等效的MD模拟方案, 通过短时间内让气体分子高频率入射到壁面, 来等效反映实际情况中长时间低频率的气固相互作用. 与文献中飞行和地面实验数据的对比表明, 这种等效方案模拟出的壁面受AO吸附/侵蚀的主要特征是合理的.
(2)对海拔$ {\text{200 km}} $以上区域, 本文分析了壁面的吸附/侵蚀情况随AO入射数目增加的演化过程. 洁净的壁面会经过3个阶段的演化, 从少量吸附AO到大量吸附AO, 最后吸附的AO会破坏壁面表层结构并形成松散氧化物. 在不同阶段, 影响动量适应过程的主要因素也不同, 适应的规律会随之变化. 其中, 一部分入射AO形成O2的现象对动量适应的促进作用不可忽略.
(3)对海拔$ {{100\sim 200\;{\mathrm{km}}}} $区域, 本文模拟了不同比例的N2/AO混合气体在壁面吸附和侵蚀的情况, 发现少量N2在壁面上离解并占据活性位点, 对后续入射AO在壁面的吸附和适应有一定影响, 在工程计算中需加以考虑.
上层大气层气-固相互作用现象非常复杂, 开展实验和计算研究难度很大, 本文的探索性研究结果对上层大气层飞行器气动设计和吸气式电推进系统设计都有一定指导意义, 目前尚未见到类似研究. 由于MD模拟结果精度依赖于势函数的选择, 本文某些定量结论还需要结合未来更精确计算或实验结果进行讨论. 另外, 高速AO对壁面的加热需要足够长的弛豫时间来冷却, 目前这种等效模拟方案的时间成本仍然较高. 对更长时间尺度气-固相互作用现象的研究, 例如材料质量损失等, 未来仍有很长的路要走.
-
![]()
图 1 (a)本文MD模拟原子力探针受力结果与文献[48]中AFM测量和DFT结果的对比; (b)本文MD模拟原子力探针算例截图; (c)本文MD模拟N2在TiO2表面吸附的截图
Figure 1. (a) Comparison of present MD simulation of AFM tip force with measurement data and DFT result in Ref. [48]; (b) Screenshot of present MD simulation of AFM tip; (c) Screenshot of present MD simulation of N2 adsorption on TiO2 surface
![]()
![]()
图 3 AO与金红石型TiO2(110)表面相互作用的模拟过程示意图
Figure 3. Schematic diagram illustrating the simulation process of AO interacting with the rutile TiO2(110) surface
![]()
图 5 不同时间步长以及等效分子数密度下, CS方案对壁面的加热效应随入射分子数目的变化
Figure 5. The heating effects on the surface in CS approach vary with the incident AO number under different time step settings and equivalent number densities
![]()
图 6 不同时间步长以及等效分子数密度下, 吸附AO的面密度${\sigma _{{\text{AO}}}}$随入射分子数目的变化
Figure 6. The surface densities ${\sigma _{{\text{AO}}}}$ vary with the incident AO number under different time step settings and equivalent number densities
![]()
图 8 壁面吸附AO的状态演化过程以及与海拔的对应关系
Figure 8. The evolution process of AO adsorption on surface and its corresponding relationship with altitude
![]()
图 9 吸附过程3个阶段中的AO面密度及吸附速率
Figure 9. The surface densities of AO and the corresponding adsorption rates in three stages
![]()
图 10 AO入射壁面可能的运动轨迹
Figure 10. The possible trajectories of AO interacting with the surface
![]()
图 11 不同吸附/侵蚀状态下壁面的AO动量适应系数(AC)
Figure 11. Momentum ACs of AO on surfaces at different adsorption/corrosion states
![]()
图 12 不同吸附/侵蚀状态下壁面上AO各类运动轨迹占比
Figure 12. The proportion of various trajectories of AO at different adsorption/corrosion states
![]()
图 13 不同海拔和不同N2/AO混合比气体在壁面的吸附/侵蚀
Figure 13. Adsorption/corrosion features of surfaces exposed in different N2/AO gas mixtures at different altitudes
![]()
图 14 不同海拔和不同N2/AO混合比气体环境中的壁面吸附面密度
Figure 14. Surface densities of adsorbed atoms in different N2/AO gas mixtures at different altitudes
![]()
图 15 N2/AO混合气体入射和纯AO入射条件下粒子各种轨迹的占比
Figure 15. Proportions of various particle trajectories on surfaces exposed in N2/AO gas mixtures and pure AO gas
![]()
图 16 不同海拔的N2/AO混合气体在吸附阶段一壁面上的NMAC结果
Figure 16. The NMAC results of N2/AO gas mixtures impinging on surfaces at the first adsorption stage at different altitudes
表 1 不同海拔H不同时刻AO在壁面的吸附速率$ {\boldsymbol{\eta}} _{{\bf{AO}}}^{{\bf{real}}} $ (单位: ${{\boldsymbol{10}}^{{\boldsymbol{14}}}}{\bf{\;c}}{{\bf{m}}^{ {\boldsymbol{- 2}}}} {\boldsymbol{\cdot}} {{\bf{min}} ^{ {\boldsymbol{- 1}}}}) $
Table 1 AO adsorption rates $\eta _{{\text{AO}}}^{{\text{real}}}$ on the surface at different altitudes H and times $({\text{unit: }}{10^{14}}{\text{ c}}{{\text{m}}^{ - 2}} \cdot {\min ^{ - 1}}) $
H/km ${n_{{\text{AO}}}}$/m−3 15 s 1 min 2 min 3 min 4 min 200 34.9 33.3 — — — — 250 12.9 46.1 — — — — 300 5.2 18.6 5.0 5.0 5.0 — 350 2.2 7.8 7.8 2.1 2.1 2.1 400 0.93 3.3 3.3 3.3 3.3 0.89 
下载: 导出CSV
-
[1] 沈清, 黄飞, 程晓丽等. 飞行器上层大气层空气动力特性探讨. 气体物理, 2021, 6(1): 1-9 (Shen Qing, Huang Fei, Cheng Xiaoli, et al. On characteristics of upper atmosphere aerodynamics of flying vehicles. Physics of Gases, 2021, 6(1): 1-9 (in Chinese) Shen Qing, Huang Fei, Cheng Xiaoli, et al. On characteristics of upper atmosphere aerodynamics of flying vehicles. Physics of Gases, 2021, 6(1): 1-9 (in Chinese)
[2] 唐绍锋, 张静. 世界主要空天飞行器研制情况及未来发展趋势. 国际太空, 2017, 466: 30-37 (Tang Shaofeng, Zhang Jing. Development status and trend of the world’s major aerospace vehicles. Space International, 2017, 466: 30-37 (in Chinese) doi: 10.3969/j.issn.1009-2366.2017.10.007 Tang Shaofeng, Zhang Jing. Development status and trend of the world’s major aerospace vehicles. Space International, 2017, 466: 30-37 (in Chinese) doi: 10.3969/j.issn.1009-2366.2017.10.007
[3] Sasaki M. Mission for super-low earth orbit, SLATS. Space Research Today, 2017, 198: 10-18
[4] Crisp NH, Roberts PCE, Livadiotti S, et al. In-orbit aerodynamic coefficient measurements using SOAR (satellite for orbital aerodynamics research). Acta Astronautica, 2021, 180: 85-99
[5] Traub C, Herdrich GH, Fasoulas S. Influence of energy accommodation on a robust spacecraft rendezvous maneuver using differential aerodynamic forces. CEAS Space Journal, 2020, 12(1): 43-63 doi: 10.1007/s12567-019-00258-8
[6] 靳旭红, 黄飞, 张俊等. 上层大气层飞行器研究进展及气动技术挑战. 航空学报, 2024, 46(1): 130254 (Jin Xuhong, Huang Fei, Zhang Jun, et al. Spacecraft in the upper atmosphere: Research development and aerodynamic challenges. Acta Aeronautica et Astronautica Sinica, 2024, 46(1): 130254 (in Chinese) doi: 10.7527/S1000-6893.2024.30254 Jin Xuhong, Huang Fei, Zhang Jun, et al. Spacecraft in the upper atmosphere: Research development and aerodynamic challenges. Acta Aeronautica et Astronautica Sinica, 2024, 46(1): 130254 (in Chinese) doi: 10.7527/S1000-6893.2024.30254
[7] 张俊, 蒋亦凡, 陈松等. 超低轨卫星气动阻力计算与减阻设计研究综述. 航空学报, 2024, 45(20): 029796 (Zhang Jun, Jiang Yifan, Chen Song, et al. Research progress in aerodynamic drag calculation and reduction design for very low earth orbit satellites. Acta Aeronautica et Astronautica Sinica, 2024, 45(20): 029796 (in Chinese) Zhang Jun, Jiang Yifan, Chen Song, et al. Research progress in aerodynamic drag calculation and reduction design for very low earth orbit satellites. Acta Aeronautica et Astronautica Sinica, 2024, 45(20): 029796 (in Chinese)
[8] Rapisarda C, Roberts PCE, Smith KL. Design and optimisation of a passive atmosphere-breathing electric propulsion (ABEP) intake. Acta Astronautica, 2023, 202(12): 77-93
[9] Zheng P, Wu J, Zhang Y, et al. Design and optimization of vacuum intake for atmosphere-breathing electric propulsion (ABEP) system. Vacuum, 2022, 195: 110652 doi: 10.1016/j.vacuum.2021.110652
[10] 岳宇贤, 耿金越, 冯冠华等. 基于吸气式电推进的超低轨飞行系统分析. 力学学报, 2023, 55(12): 2835-2843 (Yue Yuxian, Geng Jinyue, Feng Guanhua, et al. Analysis of very-low earth orbit flight system based on air-breathing electric propulsion technology. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(12): 2835-2843 (in Chinese) Yue Yuxian, Geng Jinyue, Feng Guanhua, et al. Analysis of very-low earth orbit flight system based on air-breathing electric propulsion technology. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(12): 2835-2843 (in Chinese)
[11] Romano F, Espinosa-Orozco J, Pfeiffer M, et al. Intake design for an atmosphere-breathing electric propulsion system (ABEP). Acta Astronautica, 2021, 187: 225-235
[12] Maxwell JC VII. On stresses in rarified gases arising from inequalities of temperature. Philosophical Transactions of the Royal Society of London, 1879, 170: 231-256
[13] Cercignani C, Lampis M. Kinetic models for gas-surface interactions. Transport Theory and Statistical Physics, 1971, 1(2): 101-114 doi: 10.1080/00411457108231440
[14] Lord RG. Some extensions to the Cercignani-Lampis gas-surface scattering kernel. Physics of Fluids A: Fluid Dynamics, 1991, 3(4): 706-710 doi: 10.1063/1.858076
[15] Lord RG. Some further extensions of the Cercignani-Lampis gas-surface interaction model. Physics of Fluids, 1995, 7(5): 1159-1161 doi: 10.1063/1.868557
[16] Cao BY, Sun J, Chen M, et al. Molecular momentum transport at fluid-solid interfaces in MEMS/NEMS: A review. International Journal of Molecular Sciences, 2009, 10(11): 4638-4706 doi: 10.3390/ijms10114638
[17] Agrawal A, Prabhu SV. Survey on measurement of tangential momentum accommodation coefficient. Journal of Vacuum Science & Technology A : Vacuum, Surfaces, and Films, 2008, 26(4): 634-645
[18] Nejad SM, Nedea S, Frijns A, et al. The influence of gas-wall and gas-gas interactions on the accommodation coefficients for rarefied gases: a molecular dynamics study. Micromachines, 2020, 11(3): 319 doi: 10.3390/mi11030319
[19] Tao R, Wang Z. Gas-surface interaction features under effects of gas-gas molecules interaction in high-speed flows. Chinese Journal of Aeronautics, 2024, 37(5): 228-242 doi: 10.1016/j.cja.2024.01.038
[20] Mahmoud WM, Elfiky D, Robaa SM, et al. Atomic oxygen in low earth orbits, a retrospective review study. Kosmìčna Nauka ì Tehnologìâ, 2023, 29(2): 32-44
[21] Moe K, Moe MM. Gas-surface interactions in low-earth orbit. AIP Conference Proceedings, 2011, 1333(1): 1313-1318
[22] Moe K, Moe MM. Gas-surface interactions and satellite drag coefficients. Planetary and Space Science, 2005, 53(8): 793-801 doi: 10.1016/j.pss.2005.03.005
[23] Pardini C, Anselmo L, Moe K, et al. Drag and energy accommodation coefficients during sunspot maximum. Advances in Space Research, 2010, 45(5): 638-650 doi: 10.1016/j.asr.2009.08.034
[24] Pilinski MD. Dynamic gas-surface interaction modeling for satellite aerodynamic computations. [PhD Thesis]. University of Colorado at Boulder, 2011
[25] 翟睿琼, 姜海富, 田东波等. 空间站原子氧环境仿真研究. 装备环境工程, 2014, 11(3): 35-39 (Zhai Ruiqiong, Jiang Haifu, Tian Dongbo, et al. Simulation of atomic oxygen interaction with space station. Equipment Environmental Engineering, 2014, 11(3): 35-39 (in Chinese) Zhai Ruiqiong, Jiang Haifu, Tian Dongbo, et al. Simulation of atomic oxygen interaction with space station. Equipment Environmental Engineering, 2014, 11(3): 35-39 (in Chinese)
[26] Norman P, Schwartzentruber TE, Leverentz H, et al. The structure of silica surfaces exposed to atomic oxygen. Journal of Physical Chemistry C, 2013, 117(18): 9311-9321 doi: 10.1021/jp4019525
[27] Miller SKR, Banks BA, Waters DL. Investigation into the differences between atomic oxygen erosion yields of materials in ground-based facilities and LEO. High Performance Polymers, 2008, 20(4-5): 523-534 doi: 10.1177/0954008308089711
[28] Qiao S, Jiang L, Jiang H, et al. Reactive molecular dynamics simulation on the disintegration of Kapton and Upilex-S during atomic oxygen impact. Frontiers in Materials, 2023, 10: 1-13
[29] Banks BA, Rutledge SK, Paulsen PE, et al. Simulation of the low earth orbital atomic oxygen interaction with materials by means of an oxygen ion beam. NASA Technical Memorandum, 1989, 10: 197-198
[30] Mehta PM, Walker A, McLaughlin CA, et al. Comparing physical drag coefficients computed using different gas-surface interaction models. Journal of Spacecraft and Rockets, 2014, 51(3): 873-883 doi: 10.2514/1.A32566
[31] 靳旭红, 程晓丽, 沈清等. 吸气式电推进系统进气道气体流动数值分析. 中国科学: 物理学 力学 天文学, 2024, 54(3): 234712 (Jin Xuhong, Cheng Xiaoli, Shen Qing, et al. Numerical analysis of inlet flows in an atmosphere-breathing electric propulsion system. Scientia Sinica: Physica, Mechanica et Astronomica, 2024, 54(3): 234712 (in Chinese) doi: 10.1360/SSPMA-2023-0212 Jin Xuhong, Cheng Xiaoli, Shen Qing, et al. Numerical analysis of inlet flows in an atmosphere-breathing electric propulsion system. Scientia Sinica: Physica, Mechanica et Astronomica, 2024, 54(3): 234712 (in Chinese) doi: 10.1360/SSPMA-2023-0212
[32] 周靖云, 靳旭红, 程晓丽等. 气固相互作用对吸气式电推进系统进气道性能的影响. 清华大学学报 (自然科学版), 2024, 64(9): 1536-1546 (Zhou Jingyun, Jin Xuhong, Cheng Xiaoli, et al. Effects of gas-surface interaction on the inlet performance of an atmosphere-breathing electric propulsion system. Journal of Tsinghua University (Science and Technology), 2024, 64(9): 1536-1546 (in Chinese) Zhou Jingyun, Jin Xuhong, Cheng Xiaoli, et al. Effects of gas-surface interaction on the inlet performance of an atmosphere-breathing electric propulsion system. Journal of Tsinghua University (Science and Technology), 2024, 64(9): 1536-1546 (in Chinese)
[33] Karr GR. A study of effects of the gas-surface interaction on spinning convex bodies with application to satellite experiments. [PhD Thesis]. University of Illinois, 1969
[34] Moe K. Absolute atmospheric densities determined from the spin and orbital decays of explorer VI. Planetary and Space Science, 1966, 14(11): 1065-1075 doi: 10.1016/0032-0633(66)90022-5
[35] Thompson AP, Aktulga HM, Berger R, et al. LAMMPS—a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Computer Physics Communications, 2022, 271: 108171 doi: 10.1016/j.cpc.2021.108171
[36] Stukowski A. Visualization and analysis of atomistic simulation data with OVITO—the open visualization tool. Modelling and Simulation in Materials Science and Engineering, 2010, 18(1): 015012 doi: 10.1088/0965-0393/18/1/015012
[37] Reddy MR. Effect of low earth orbit atomic oxygen on spacecraft materials. Journal of Materials Science, 1995, 30(2): 281-307 doi: 10.1007/BF00354389
[38] Smith CTG, Delkowki M, Anguita JV, et al. Complete atomic oxygen and UV protection for polymer and composite materials in a low earth orbit. ACS Applied Materials and Interfaces, 2021, 13(5): 6670-6677 doi: 10.1021/acsami.0c21552
[39] Perron H, Domain C, Roques J, et al. Optimisation of accurate rutile TiO2 (110), (100), (101) and (001) surface models from periodic DFT calculations. Theoretical Chemistry Accounts, 2007, 117(4): 565-574 doi: 10.1007/s00214-006-0189-y
[40] 郭玉宝, 杨儒. 金红石型TiO2(110)表面原子结构和电子结构的理论研究. 北京化工大学学报, 2004, 31(5): 64-68 (Guo Yubao, Yang Ru. Theoretical study of atomic structure and electronic structure of nanosized rutile TiO2(110) surface. Journal of Beijing University of Chemical Technology, 2004, 31(5): 64-68 (in Chinese) Guo Yubao, Yang Ru. Theoretical study of atomic structure and electronic structure of nanosized rutile TiO2(110) surface. Journal of Beijing University of Chemical Technology, 2004, 31(5): 64-68 (in Chinese)
[41] Van Duin ACT, Dasgupta S, Lorant F, et al. ReaxFF: A reactive force field for hydrocarbons. Journal of Physical Chemistry A, 2001, 105(41): 9396-9409 doi: 10.1021/jp004368u
[42] Aktulga HM, Fogarty JC, Pandit SA, et al. Parallel reactive molecular dynamics: Numerical methods and algorithmic techniques. Parallel Computing, 2012, 38(4-5): 245-259 doi: 10.1016/j.parco.2011.08.005
[43] Rahnamoun A, Van Duin ACT. Reactive molecular dynamics simulation on the disintegration of Kapton, POSS polyimide, amorphous silica, and teflon during atomic oxygen impact using the reaxff reactive force-field method. Journal of Physical Chemistry A, 2014, 118(15): 2780-2787 doi: 10.1021/jp4121029
[44] Groh KK De, Banks BA, Mccarthy CE, et al. MISSE Peace polymers atomic oxygen erosion results//2006 MISSE Post-Retrieval Conference, 2006. NASA/TM-2006-214482
[45] Kim SYY, Van Duin ACTT. Simulation of titanium metal/titanium dioxide etching with chlorine and hydrogen chloride gases using the ReaxFF reactive force field. Journal of Physical Chemistry A, 2013, 117(27): 5655-5663 doi: 10.1021/jp4031943
[46] Monti S, Li C, Ågren H, et al. Dropping a droplet of cysteine molecules on a rutile (110) interface: Reactive versus nonreactive classical molecular dynamics simulations. Journal of Physical Chemistry C, 2015, 119(12): 6703-6712 doi: 10.1021/acs.jpcc.5b00932
[47] Barcaro G, Sementa L, Carravetta V, et al. Experimental and theoretical elucidation of catalytic pathways in TiO2-initiated prebiotic polymerization. Physical Chemistry Chemical Physics, 2019, 21(10): 5435-5447 doi: 10.1039/C9CP00167K
[48] Sokolović I, Reticcioli M, Čalkovský M, et al. Resolving the adsorption of molecular O2 on the rutile TiO2(110) surface by noncontact atomic force microscopy. Proceedings of the National Academy of Sciences of the United States of America, 2020, 117(26): 14827-14837
[49] Tian N, Comer BM, Medford AJ. Screening and discovery of metal compound active sites for strong and selective adsorption of N2 in air. Chem Sus Chem, 2023, 16(22): e202300948 doi: 10.1002/cssc.202300948
[50] Lee J, Liu X, Kumar A, et al. Phase-selective active sites on ordered/disordered titanium dioxide enable exceptional photocatalytic ammonia synthesis. Chemical Science, 2021, 12(28): 9619-9629 doi: 10.1039/D1SC03223B
[51] Emmert JT, Drob DP, Picone JM, et al. NRLMSIS 2.0: a whole-atmosphere empirical model of temperature and neutral species densities. Earth and Space Science, 2021, 8(3): e2020EA001321 doi: 10.1029/2020EA001321
[52] Norman P, Schwartzentruber T, Cozmuta I. A computational chemistry methodology for developing an oxygen-silica finite rate catalytic model for hypersonic flows//42nd AIAA Thermophysics Conference, Reston, Virigina: American Institute of Aeronautics and Astronautics, 2011: 1-17
[53] Minkowycz WJ, Sparrow EM. Advances in Numerical Heat Transfer, Vol. 2. CRC Press, 2000
[54] Reinhold J, Veltzke T, Wells B, et al. Molecular dynamics simulations on scattering of single Ar, N2, and CO2 molecules on realistic surfaces. Computers and Fluids, 2014, 97: 31-39 doi: 10.1016/j.compfluid.2014.03.024
[55] Fu Y, Peng F, Zhang C, et al. Multi-faceted simulation of atomic oxygen erosion of deorbit sail for cleaning space debris. Acta Astronautica, 2021, 187(3888): 61-69
[56] Cui Z, Zhao J, Yao G, et al. Molecular insight of the interface evolution of silicon carbide under hyperthermal atomic oxygen impact. Physics of Fluids, 2022, 34(5): 052101 doi: 10.1063/5.0092005
[57] Cui Z, Zhao J, He L, et al. A reactive molecular dynamics study of hyperthermal atomic oxygen erosion mechanisms for graphene sheets. Physics of Fluids, 2020, 32(11): 112110 doi: 10.1063/5.0030749
[58] Filinov A, Bonitz M, Loffhagen D. Microscopic modeling of gas-surface scattering I. A combined molecular dynamics-rate equation approach. Plasma Sources Science and Technology, 2018, 27(6): 064003
[59] Filinov A, Bonitz M, Loffhagen D. Microscopic modeling of gas-surface scattering: II. Application to argon atom adsorption on a platinum (111) surface. Plasma Sources Science and Technology, 2018, 27(6): 064002
[60] Murray VJ, Zhou L, Xu C, et al. Scattering dynamics of glycine, H2O, and CO2 on highly oriented pyrolytic graphite. The Journal of Physical Chemistry C, 2019, 123(6): 3605-3621 doi: 10.1021/acs.jpcc.8b11293
[61] Rutigliano M, Pirani F. On the influence of rotational motion of oxygen molecules on the scattering from graphite surfaces. The Journal of Physical Chemistry C, 2019, 123(18): 11752-11762 doi: 10.1021/acs.jpcc.9b01966
[62] Li M, Hebenstreit W, Diebold U, et al. Oxygen-induced restructuring of rutile TiO2(110): Formation mechanism, atomic models, and influence on surface chemistry. Faraday Discussions, 1999, 114(110): 245-258
[63] He L, Cui Z, Sun X, et al. Sensitivity analysis of the catalysis recombination mechanism on nanoscale silica surfaces. Nanomaterials, 2022, 12(14): 1-21
[64] Wood BJ. The rate and mechanism of interaction of oxygen atoms and hydrogen atoms with silver and gold. Journal of Physical Chemistry, 1971, 75(14): 2186-21895 doi: 10.1021/j100683a021
[65] Moe K, Moe MM, Wallace SD. Improved satellite drag coefficient calculations from orbital measurements of energy accommodation. Journal of Spacecraft and Rockets, 1998, 35(3): 266-272 doi: 10.2514/2.3350
计量
- 文章访问数: 231
- HTML全文浏览量: 97
- PDF下载量: 63
