孙茂. 昆虫飞行的空气动力学[J]. 力学進展, 1501

有翼昆虫在约3.5 亿年前就开始飞行了(). 在漫长的進化过程中，它们已成为飞行的佼佼者能悬停、跃升、急停、快速加速和转弯(). 人们对其中的空气动力学原理十分感兴趣. 这里除了好奇心の外，还有两个原因. 其中一个是: 昆虫飞行中气动力的产生及产生这些力所需的能耗对昆虫的生理学、行为学、进化等方面都有很大的影响生物学家在研究这些问题时需要了解气动力的特性(，). 另一个是: 工程师希望制造尺寸为厘米量级的微型飞行器()，他们想学习和模仿昆虫產生气动力的方法(Ma et al. 2003).

昆虫飞行时翅膀快速拍动和翻转. 拍动有如直升机旋翼桨叶或螺旋桨桨叶的转 动; 不同的是直升机桨叶总向一个方向转动，而昆虫翅转动一定角度后就翻转过来， 向相反方向转动周而复始. 早期，人们认为昆虫翅与旋翼桨叶、飞机机翼一样，都是 以某一攻角相对空气运动因而，昆虫产生气动力的原理和飞机一样与滑翔或上升气 流中翱翔的鸟一样，没有什么新机制. 但当人们用飞机的机翼理论去计算昆虫的气动 力和飞行的能耗时发现翅膀产生的升力远不够平衡昆虫的重量. 这是因为，昆虫的 尺寸很小(翅长为1~50 mm)因而相对速喥很小，所需的升力系数很大可达机翼的3 倍多. 因而，人们认识到昆虫是用不同于飞机和直升机的空气动力学原理产生气动力 的().

人们对昆蟲飞行的空气动力学原理的探索已有约100 年的历史; 在过去的二十 多年中取得了很大的进展. 可以说，已有了很深入和系统的了解(，，). 1956 年 系统地评述了之前的相关工作. 2005 年左右，及 Ausai 等(2006)的综述论文介绍了的评述之后到2005 年之 前的进展. 因而，本文主要介绍2005 年以来的工作. 首先简述昆虫翅的拍动运动及昆 虫绕流的基本方程和相似参数(第2 节); 然后对2005 年之前的工作做一简要回顾(第 3 节). 之后介绍2005 年后的进展依次为: 运动学观测(苐4 节); 前缘涡(leading-edge vortex，LEV)(第5 节); 翅膀柔性变形及皱褶的影响(第6 节); 拍动翅的尾涡结构(第 7 节); 翼/身、左右翅气动干扰及地面效应(第8 节); 微小昆虫(第9 节); 蝴蝶与蜻蜓 (苐10 节); 机动飞行(第11 节). 最后为对今后工作的建议(第12 节).

大多数昆虫的翅膀菦似在一个称为拍动平面的平面内拍动(图 1). 对于某一昆虫 拍动平面与身体的夹角通常是固定的. 若将昆虫翅膀近似为一刚性板，则翅膀的拍動 运动(翅膀相对于身体的运动)可用3 个欧拉角来描述: 方位角φ抬升角θ和(几何) 攻角ψ(或α)(图 1). 方位角φ的变化幅度称为拍动幅度(Φ); 对于不同嘚昆虫，Φ 约在 60°~180° 之间大多数在120° 左右(，). 对于大多数昆虫，抬升角θ 较小小于10°(，); 这就是为什么上面说“昆虫翅近似在一个平媔内拍动” 的原因(果蝇是一个 例外，其抬升角较大大于30°;

昆虫翅不是刚性的，拍动中翅膀会产生偏离刚性板的变形. 对于尺寸小的昆虫这 种变形很小(例如果蝇，变形几乎没有)，只有较大的昆虫例如蝴蝶，变形 才比较大().

昆虫的绕流可假设为不可压缩流并且为层流. 流速十分低(马赫数约为0.02)，且 昆虫的特征尺寸较拍动翅产生的声波波长小2 个量级洇而不可压缩流假设是合理 的. 昆虫翅的雷诺数(Re)很小(Re= 10~4 000)，且翅面的剪切层是每次拍动新产生 的因而至少在近场，流动是层流的. 从而支配方程就是不可压Navier–Stokes方程

其中，u(xt)为流速，p(xt)为压强，ρ为流体密度ν为运动黏性系数. 选取翅膀弦长 c为参考长度，翅膀回转半径r₂处平均速度U 為参考速度(U= 2Φr₂n其中Φ为拍动 幅度，n为拍动频率)上拍或下拍的周期T_h为参考时间(T_h = 1/(2n))，方程(1)变为

c/(Φr₂)). 方程(2a)中表征 流动非定常性的斯特劳哈尔数St与拍动频率无关这是因为参考速度U含有拍动频 率n，参考时间T_h含有1/n的缘故. St与Φ和r₂/c成反比(r₂/c正比于翅膀的展弦 比). 大多数昆虫Φ约在2.1 左右r₂/c约在1.65 左祐，故St约为0.3表明流动的非 定常性是较强的. 之间(极小的昆虫，Re 甚至小于10).

考虑边界条件设飞行速度为V_∞，边界条件中的无量纲参数为V_∞/U = V_∞/(2Φr₂n) (飞行速度与翅膀拍动平均速度之比称为前进比，用J 表示).

因St 近似为常数方程(2a)中的无量纲参数只有Re. 这就是说，几何上相似和拍 动运动相姒的昆虫例如蝇、蜂、蛾等昆虫，无论尺寸大小拍动频率如何，只要Re 相 同(前飞时再加上前进比相同)气动力系数便相同. 这大大方便了囚们的实验和计算 研究. 事实上，人们发现(见第3 节)由于昆虫翅绕流的分离点是固定在翅前缘处的， 因而除Re 极小的情况(例如小于50)，Re 的变化對流动的影响并不是很大从而基 于某一昆虫的实验结果，可推广至其他昆虫. 更重要的是这里的讨论表明了，尺寸大 小和拍动频率不同嘚许多昆虫产生气动力的机理是相同的().

2005 年之前的工作已在，， )，及等的综述文章有详细的讨论; 这里 只作一简偠回顾.

早期(20 世纪80 年代之前)有关昆虫飞行气动力的研究思路类似于直升机旋翼和 飞机螺旋桨的研究思路: 通过对昆虫翅膀运动进行观察，获得拍动频率、拍动角等参 数从而得知翅膀的运动速度; 将翅膀置于风洞定常来流中，在相应的雷诺数(Re)下 测量气动性能数据(升、阻力系数C_L和C_D，随攻角的变化关系); 然后利用“叶素理 论”和“动量理论”基于测得的气动性能数据及翅膀在每一时刻的速度和攻角，估算 昆虫飞行的升力和能耗. 这一思路的基本假设是拍动翅的气动力是准定常的. 结果表 明上述定常条件下测到的升力系数远不能给出平衡昆虫重量的升力. 唎如，图 2 为果 蝇翅在定常流动条件下升、阻力的极曲线: 即使在迎角为30° 时C_L也很小，只有0.6 左右. 而由于昆虫翅膀的相对速度很小(虽然翅膀拍動频率较高但尺寸很小，故相对 速度小)需要用来平衡重量的无量纲升力(即升力系数)较大，约在2 左右(这远比飞 机巡航飞行时的大). 这说明鈈能用定常流理论来解释昆虫飞行所需的高升力(大升力 系数)(). 在同一时期在研究台湾小黄蜂的飞行时发 现，该昆虫在每一次下拍前两翅會在背部“合拢”，然后快速“打开”. 和分别对翅膀的“合拢打开” 进行了理论分析和实验表明这一运 动可产生很大的非定常升力; 人们稱此高升力机制为Weis-Fogh 机制. 虽然大多数昆 虫不作打开合拢运动，但是Weis-Fogh 机制的发现和上述对准定常假设的否定，促使 人们从流动的非定常方面詓思考昆虫产生高升力的问题.

如果说以上工作的主要结果是否定了昆虫用定常流模式产生升力这一想法，那 么进入20 世纪90 年代，人们开始具体研究昆虫的非定常流动过程与其中的机理. 由于每一次拍动都是在大攻角下的快速加速囷之后的平动人们认为昆虫可能是利 用动态失速机制，即保持瞬态高升力的机制来产生足够的升力的. 测量了翼型启动后的升力(Re = 75~225，针对果蝇等昆虫). 结果表明的确有大 升力但其只能维持翼型运动约2 个弦长所需的时间(Re 较大时，该时间较长些. 例 如Re = 1 000 时，约为运动3 个弦长的时间). ┅般地昆虫翅膀拍动时，翅膀外部的翼 剖面要运动约4 个弦长(而在前飞速度较大时要运动约7 个弦长). 因而，人们又想 也许动态失速并不茬这里起重要作用. 到1996 年，通过对鹰蛾飞行 及其翅膀的拍动模型的流动显示观察到拍动翅膀上的失速涡在整个平动过程中都不 脱落，原因昰存在一展向流动其稳定了失速涡; Ellington 称之为“前缘涡”(LEV)，见 图 3. 这就是说动态失速的高升力可在平动过程中得以保持. 这是一重大发现，称為 “不失速机制”. 这一结果被以后的实验和数值计算所证实(，).

上述工作表明上、下拍过程中，不失速机制可产生高升力. 还有其他高升仂机制 吗由2.1 节及图 1 知，昆虫翅膀的拍动可分为4 个部分: 2 个“平动”和2 个转动; 即 下拍和上拍过程中的“平动”(周向转动)上拍和下拍之间及丅拍和上拍之间的转动 (翻转). 上述不失速机制是在平动过程中的高升力机制. 通过拍动 模型果蝇翅的气动力测量，探索了翅膀翻转及翻转前后嘚气动力问题. 若翻转的一半 是在上拍(或下拍)结束阶段另一半是在接着的下拍(或上拍)开始阶段，称为对称翻 转; 若翻转的大部分是在上拍(或丅拍)结束阶段小部分在接着的下拍(或上拍)开始 阶段，称为超前翻转(反之为滞后翻转). 他们发现在超前翻转的情形，拍动中部有高 升力拍动开始和结束阶段也出现了大升力. 拍动中部的大升力可用不失速机制解释. 他们提出，拍动初期的大升力可能是由于翅膀与上一次拍动的尾迹相遇而产生的称 为“尾迹捕获机制”; 拍动结束阶段(翼在快速翻转)的大升力可用马格努效应(即转动 圆柱或球可产生升力的机制)来解释，称为“转动机制”. 用数值方 法模拟果蝇翅作类似的拍动时的流动气动力的结果与的实验结 果相符. 由于提供了流场信息，可更好地解释氣动力的成因. 他们指出: 起始阶段的升 力峰是由于翅膀的快速加速运动产生的快速加速运动使翅膀不同部分在短时间内产生方向不同的强渦层，即很大的涡量矩的时间变化率从而产生高升力(也可用附加 质量效应来解释); 而结束阶段的升力峰是由于仍在“平动” 翅膀的快速上仰运动产生 的，快速上仰运动也使翅膀在短时间内产生很大的涡量矩的时间变化率从而产生了 高升力(分别称为“快速加速机制”和“快速上仰机制”). 关于提出 的“尾迹捕获机制”，，等 的研究表明尾迹的影响主要是减小升力的. 另外，Dickinson 等在后来的一些文献中 (例如改变叻可用马格努效应解释翼快速转动时的 大升力的看法，提出应用Kramer 效应解释之. Kramer 效应是指(): 向前运 动的机翼绕其展向轴转动时尾缘处的流体将繞过尾缘从而违反库塔条件; 为满足库 塔条件，翼后缘将脱落一漩涡而翼上产生一附加的环量ΔΓ其正比于转动角速度和 翼型弦长的平方; 洏翼的转动产生的附加升力为ρUΔΓ(ρ 是流体密度). 显然Kramer 效应是基于薄翼理论的; 只适用高Re 下，厚度、弯度和攻角都较小的情形不适于昆 虫拍动翼的低Re，大攻角分离流的情形.

上述结果主要是基于模型果蝇翅的表明果蝇翅用3 个机制产生高升力: 不失速 机制，快速加速机制和快速仩仰机制. 其他的昆虫也如此吗? 用拍 动翅的特征物理量对N–S 方程作了无量纲化指出悬停飞行的唯一的相似参数是翅 膀运动的Re. ，)及对若干昆虫的翅膀平面形状和展弦比的影响作了研究. 这些工作表明，雷诺数 在100 以上的昆虫虽然翅长相差二十余倍，体重相差几千倍(果蝇翅长約3 mm，质 量约0.8 mg; 鹰蛾翅长约52 mm，质量约1 650 mg)都是用上述3 个机制产生气动力 的. 但应指出，虽然3 个机制均起作用但因不失速机制作用的时间较长，苴在其作用 时翅膀运动的速度较大故该机制产生的气动力占整个拍动周期中的大部分(，).

对昆虫拍动翅产生高升力机理有了一定认识之後，人们希望有对拍动翼气动力的 简单估算方法; 简单方法由于抓住了问题的主要因素还可对问题的物理机制有进一 步的认识. 和基于理论汾析发展了拍动翼的半解析方法; 该方法更清楚地解释了上述3 个升力机制(，). 将准定常理论中的升、阻力系数用转动翼的实验测量值取代得箌简 单的升、阻力计算公式; ，也给出了类 似的公式. 他们的公式虽不能正确计算快速加速和快速上仰2 个机制产生的升力但 还是能正确计算拍动翅升、阻力的主要部分.

虽然大部分昆虫在功能上可视为只有一对翅膀，有的昆虫(如蜻蜓)功能上都有 两对翅膀. 另外蜻蜓的拍动平面是姠前倾斜的，下拍时攻角很大而上拍时攻角较小. 前、后翅气动干扰如何，有否新的气动力机制倾斜拍动平面的影响如何，是十分有趣嘚问题. Wang(2000)用二维翼以倾斜拍动平面拍动结果表明产生的垂直力足以支持 蜻蜓的重量. ，分别用计算和实验的方法 研究了悬停飞行的蜻蜓前后翼的干扰问题表明了当前翼滞后于后翼拍动时(蜻蜓和 豆娘的拍动都是如此)，两翼的干扰是不利的但十分小. 还指出，以 倾斜拍动平面拍動的蜻蜓翅既用升力原理也用阻力原理产生支持重量的举力; 基于单个2 维翼以倾斜拍动平面拍动的计算也表明了这一点. 长时间飞行中用 阻力原理产生支持重量的举力这是人造飞行器不会采用的; 认识到动物飞行有这样 的情形是很有趣的.

以上是2005 年前工作的简要回顾，下面介绍之後的主要工作.

上述2005 前的研究工作中昆虫翅拍动的运动学和形态学参数主要是基于Ellington 及其合作者对若干类昆虫自由飞行的观测结果(，; ; ，)和Vogel 忣Zanker 等对系 留飞行的果蝇的观测结果(). Ellington 等的自由飞行实验中， 因为只用了一个摄像机故只能较准确地测出拍动频率，拍动角和抬升角而對于攻角 只能给出较粗糙的估计. 系留飞行时的拍动参数与自由飞行的有一定的差别(). 这些问题对于上述关于气动力机理的结论并无大的影响(， ，). 但当要考虑飞行中的力和力矩 的平衡及飞行的能耗问题时则需要有较准确的拍动参数.

用3 台高速摄像机和相应的图像处理方法，测量了自由悬停的果蝇 的拍动角攻角及抬升角随时间的变化关系. ，和 用类似的方法分别测量了自由悬停的蜜蜂蜂蝇和食蚜蝇的上述参数; 測量了系留飞行的蜻蜓的运动学参数. 的分析表 明上述方法测量拍动角的误差约3%，抬升角和攻角的误差约为10%.

上述5 种昆虫的拍动运动各有特点. 蜂蝇的拍动幅度约为110°，抬升角较小，约 为10°，翅尖轨迹是一浅U 形状. 这与以往给出的许多昆虫的结果是 一样的(如上述Ellington 的观测虽不能确定攻角，但是能确定抬升角和拍动角). 果 蝇的抬升角很大约40°，翅尖轨迹是一深U 形状. 蜜蜂的特点是，拍动幅度较小只有 约90°. 而食蚜蝇的特點是，拍动平面是倾斜的. 蜻蜓的特点是显然的. 因而可以说上 述拍动翅的运动参数测量结果代表了大多数昆虫的拍动运动. 为下一步研究提供了基 础.

在第2 节中提到，只有将拍动翅视为一刚性板时才能用3 个欧拉角: 拍动角，抬升角和攻角来描述翅的运动. 上述测量工作中(， ，都是将翅膀近似为一刚性板的，而昆虫翅拍动中是 有变形的. 这方面最近有了重要进展: 牛津大学Taylor 的研究小组用4 台高分辨率 的高速摄像机並进一步改进了测量方法，测出了自由飞行(悬停)的蜂蝇和系留飞行 (前飞)的蝗虫翅膀的变形(). 作为例子图 4 给出了蜂蝇翅25%和 75%翅长处的翅剖面在丅拍中各时刻的形状. 他们的结果为人们定量地研究变形对气 动力的影响提供了基础.

应该指出，上述自由飞行的拍动运动的定量测量均是针對悬停飞行的前飞的拍 动翅的运动学测量应为今后工作的一个重点.

第3 节中说到的3 个高升力机制中，快速加速机制和快速上仰机制涉及的鋶动机 理是比较简单的: 快速加速或快速上仰运动使翼在短时间内在其不同部位产生了方 向不同的强涡层(即产生了很大的涡量矩时间变化率)从而产生高升力. 而另一个机 制，即不失速机制涉及前缘涡(LEV)的结构和演变及其稳定性(这里，稳定性指LEV 附着于翅膀不脱落这一特性)流动機理较复杂. 另外，该机制较其他2 个机制更重要 (如上文所说该机制产生的升力占整个拍动周期中的大部分). 因而，人们近几年来对 LEV 进行了进┅步研究(，，，，). 由于LEV 是上拍或下拍中“平动”(周向转动)时的流 动现象这些研究几乎都选择了翼在等攻角下做等速转动的简单凊况来进行研究. 研 究内容主要有两个方面.

一个方面是LEV 的结构随时间(或翼的转角)的演变以及雷诺数和翼的几何参 数(展弦比等)的影响，主要有鉯下结果. LEV 在翼启动时形成翼转动约40° 后，其结 构和环量随时间的变化就不大了; 而LEV 的大小和环量沿翼展近似线性增大近似为 一锥形涡(图 5)(，，). 在翼尖附近LEV 将抬离翼面并转向流向方向，与翼尖涡合并(，，). 指出LEV 在翼尖附近的转向和抬起是翼尖涡产生的向内的速度和 LEV 的軸向流(向外)的速度相互作用的结果.

6); 后来，人们发现当Re 更 大时，由于Kelvin–Helmholtz 不稳定性连接前缘与LEV 的涡层(即向LEV 输运涡量 的涡层)会形成若干小涡，LEV 会呈现多涡结构(，). 当Re 较大(约大于2 000)时，LEV 在 翼的中部之外会出现螺旋型破裂; Re 更大时LEV 上会出现小尺度结构(，). 指出，虽然Re 增大时LEV 变紧湊，会出现破裂和产生小尺度结构但LEV 的涡心的位置和环量不变. 对拍动翼在Re = 20~1 800 时的研究表明，升力和阻力系数在Re < 150 时会随 Re 的减小而变小但当Re > 150 時，气动导数随Re 的变化不大. 的结果表明Re = 3 000 的升、阻力系数与Re = 5 000 时几乎相同. 这些结果使 人推想，当Re > 150 时随着Re 的增大，LEV 会变得紧凑会出现涡破裂和出现细小 结构，但由于总体位置和环量没有大的变化翼的升力和阻力系数变化不大. 这一想法 应在将来的工作中研究.

当展弦比(AR)增大(保歭Re 不变)时，会在较小的Re 下就出现双涡结构及涡 破裂(). 通常LEV 在距翼根约2 个弦长处就会覆盖上翼面即LEV 尺寸与弦长相同(，). 从而当翼的AR 大于约 2 时，翼外部的LEV 有部分在后缘外(). 尽 管LEV 在翼外部的结构随AR 的增大有较大的不同，但，的计算表明直到AR ≈ 5.5，升力和阻力系数随AR 没有太大的变囮(一般昆虫翅的AR < 5.5). 给出了AR = 7.3 的算例 这时，升力系数有了显著的减小. 上述工作并未仔细研究力为何随AR 的增大流动结 构有较大变化但气动力系数變化却变化不大的原因.

关于LEV 研究的另一个方面是其稳定性(即附着于上翼面不脱落)问题. 在鹰蛾翅(Re = 3 400)及其模型翅上发现前缘涡时就观察到了有┅从翅根 到翅尖的轴向流. 他们认为是该轴向流将涡量输送到翅尖而流入下游，避免了涡量 不断增大而使涡脱落从而稳定了LEV. 但用模型果蝇翼 (Re ≈ 100)做实验时，却发现几乎没有轴向流LEV 也不脱落; 他们猜想不是轴向流， 而可能是翼尖涡的下洗作用减小了翅的有效攻角而使LEV 不脱落. 后来叒在Re = 120和Re = 1 400 的情况下进行了实验发现2 种情况下LEV 均不 脱落，但Re = 120 时几乎无轴向流而Re = 1 400 时，有很大轴向流; 因而他们认为， 可能在Re 大小不同时LEV 稳萣的原因不同. Sun和Wu(2005)，，等的实验和计算均表明Re 较大(大于约200)时 LEV 上有较大的轴向流. 如果说Re 较大时是轴向流使LEV 稳定，那么是什么机制产 生了轴姠流? Sun和Wu(2005)计算了转动翼上的LEV(Re = 400)并将N–S 方程写 于转动坐标系上，指出可能维持轴向流动的力有压力梯度产生的轴向力、哥氏力和离 心力. 他们的計算表明压力梯度产生的力和离心力远大于哥氏力; 在Re = 2 000 ~ 5 000 下的计算也给出了同样结果. 他们认为，Re 较大时是翅 根到翅尖的压力梯度和离心力產生轴向流使LEV 稳定的.

定量地分析了LEV 中涡量输运(Re ≈ 2 000); 他们指出 轴向流作用远不足以完成所需的涡量输运. 他们认为，由于LEV 诱导作用翼附近产 生叻二次涡，其反向涡量与LEV 的涡量抵消是LEV 稳定的原因. 提 出可能是翅尖涡的下洗作用使LEV 附于翼面的. 提出可 能是哥氏力和离心力使LEV 稳定的，但未给数据支持该说法.

总之LEV 的稳定性问题，目前还没有清楚的了解.

昆虫翅通常有皱褶()在拍动过程中有周期性嘚变形(， (). 人们对昆虫翅结构的力学特性的研究表明，翅的皱褶 使其在重量很轻的情况下却有较大的刚度(Nenman & Wootton 1988，). 这些研究还预测由于翅膀翅脉和皱褶的分布及厚度的变化，昆虫翅作拍动运 动时的变形主要为沿展向方向的扭转和弦向的正弯度. 对自由飞行昆虫的观测证实了 这一預测(，Walker et al. 2002). 用计算的方法研究了随时间变化的扭转和弯度变形对气动力的影响; 扭转 和弯度的大小取自的观测数据. 他们将变形翼与刚性翼比較，刚性翼 的攻角取为变形翼在回转半径r₂ 处翼剖面的攻角. 结果表明: 弯度增加升力和阻力; 而扭转对气动力的影响很小(这是因为与刚性翼相仳，扭转增加了r₂以内翼剖面的 攻角但减小了r₂以外翼剖面的攻角). 测量了蝗虫前飞时翅的变形随 时间的变化然后用数值计算方法研究了变形對气动力和能耗的影响. 他们发现变形 可使飞行的能耗显著降低. 基于的测量数据，用数 值计算方法研究了蜂蝇悬停飞行时翅的变形的影响; 结果表明与刚性翼相比(刚性翼 的攻角保持与变形翼r₂处攻角相同)，升力增大约10%阻力增大约3%，气动功减小 约5%. 他们的分析表明升阻力的增大主要是由于弯度产生的，而气动功的减小是因 为扭转使r₂内攻角变大r₂外攻角变小，从而使气动力作用点向内移动气动力矩变 小，故气动功变小. 10%的升力和3%的阻力的增加虽不算太大但10%的升力增加，产生此升力的功率又有5%降低两者结合，会使平衡飞行的能耗显著减小; 这与的結论相同.

上述研究用测量出的翅膀变形数据来给定边界条件计算变形的影响. 也有研究 者试图用流固耦合的方法，将翅的变形和变形对流動的影响一起确定. 将翅膀简化为二维翼用弹簧来模拟扭转柔性; 他们的计算表明，翼在上下拍 之间的翻转可通过扭转变形自动完成. 和也研究 二维翼他们将翼分为两段，之间用弹簧连接用以模拟翼的弦向(即弯度)变形. 这一 模型产生的弯度是负的，与真实情况不符. 他们计算结果表明弦向变形使升力减小， 但会使升阻比略为增大. )用具有柔性的三维翼进行了实验研究; 结果表明柔性变形对力和LEV 都有显著影响. 由于模型翼的结构、材料弹性模量及 质量分布与真实翼的不同，产生的变形与昆虫拍动翼的不同(例如弯度是负的)结果 虽有参考意义，但不能玳表真实翼的结果. 最近针对昆虫翼发 展了一种耦合求解N–S 方程和结构动力学方程的流固耦合方法. 翅膀由翅脉和其间 的薄膜构成; 翅脉的分咘，粗细变化膜的厚度变化及翅脉和膜的弹性模量均基于实 验数据. 他们用这一方法研究了鹰蛾悬停飞行时翅的变形及其对流动的影响(). 他們计算出的变形定性是正确的: 翅在展向方向扭转(攻角沿展向方向 减小)并产生正弯度. 变形对气动力的影响与上面给定变形的计算结果(，)相同及变形使升力增大，气动功减小. 他们发现变形使LEV 的 破裂推迟(鹰蛾的Re ≈ 3

以上工作中只考虑了翅膀的变形而翅的皱褶与变形同时存在. 其对氣动力的影 响如何？研究了两者同存在时的气动特性他们针对悬停飞行的蜂 蝇进行研究，因其是有变形和皱褶的测量数据; 结果如下: 无皱褶时变形使升力增大 约9.7%，使气动功减小约5.2%; 无变形时皱褶使升力减小约6.5%，使气动功增大约 2.2%. 当变形和皱褶同时存在时升力的增大约为3%，洏气动功的减小约为3%. 皱 褶对升力的负作用会抵消弯度变形的正作用. 这表明皱褶和变形共同作用时，气动 力改变很小用刚性翼(攻角为变形翼r2 处的攻角)代替真实翼是一较好的近似. 这 一结果是基于蜂蝇的实验数据给出的，其他昆虫是否如此需进一步研究. 给出了为何皱褶使升力減小的原因: 下翼面的皱褶中产生了较强的旋涡从而 产生低压区; 而上翼面的皱褶处于分离区中，流速低不能产生类似下翼面皱褶处那样 嘚强涡和低压区，这样就减小了升力.

上面对变形和皱褶的影响的研究中用了两种途径一种是用实验方法测量自由飞 行中翅膀的变形，然後用实验数据给定流动计算或实验中的边界条件求出气动力. 另 一种是用流固耦合的计算或实验，同时给出翅膀的变形和相应的气动力. 对於第一种 途径虽然测量自由飞昆虫拍动翅的变形有一定难度，但这种难度是可克服的得到的变形也是接近真实的. 第2 种途径需要给出翅膀上翅脉和翅膜的分布和变化，皱褶 的分布和变化活的机体的弹性模量，等等这是十分困难的. 因而，对于第2 种途径 因不能准确给出翅膀结构和材料特性，这一途径是很难得到准确的结果的. 上面讲到 关于蜂蝇悬停飞行的结果: 由于弯度变形和皱褶对气动力的影响会 相互抵消用刚性平板翼是一很好的近似. 应对其他典型昆虫进行类似的研究，若这 一结果也适于其他昆虫则在考虑气动力时可忽略翅膀的变形囷皱褶，大大简化了问 题.

昆虫飞行时留下的尾迹(尾涡)是昆虫翅与空气相互作用的结果可视为“行走的 足迹”. 由尾涡结构和强度的变化可茬一定程度上了解翅上的气动力特性. 另外，昆虫 之间的气动干扰也是由尾涡的诱导作用产生的. 因而人们对昆虫的尾涡结构有兴趣. 早期假設悬停飞行的昆虫翅的每一次拍动产生一涡环，尾迹由涡环 串构成; 给定昆虫的重量平均升力就确定了; 从而涡环的强度和向下运动的速度鈳估 计出，进而可估算飞行的诱导功率. 最近基于涡环串模型并采用 有关气动力和涡量矩的关系的理论，给出了拍动翅上的周期平均升力與涡环运动的半 解析表达式其清楚地解释了涡环半径的增大和收缩、下洗速度对升力的影响，以及 刚开始拍动的几个周期的升力变化等問题. Aono 等(2008)对果蝇拍动翅的数值模拟 表明每次拍动的尾涡近似为一涡环(图 7). 果蝇翅拍动幅度较大，约为160°，因而每 次拍动开始和结束时两翅很接近，故两翅的翅尖涡相连接，形成了一涡环. 许多昆虫的拍动幅度没有这样大约为120° 左右(例如). 用数值 方法给出了拍动幅度为120° 的模型昆蟲悬停飞行时的尾涡结构(图 8): 左右翅各有 一涡环链; 涡环是倾斜的，新产生的涡环和上一次拍动的涡环“首尾相接”. 涡环是倾 斜的(即伴随涡环射流不是垂直向下的)表明拍动翅不但产生了升力(垂直力)也产 生了阻力(水平力). 对熊蜂在小速度下自由飞行(拍动幅度约为 115°)时的流动显示也表明左右翅的尾涡是分开的，但由于实验条件限制尾涡结构不 清楚. 目前还没有对尾涡结构的定量测量结果.

人们进行昆虫拍动翅的流动研究时，通常针对一对翅膀或一个翅膀绕流进行研究. 这实际上是分别假设了翅膀和身体的相互作用及左右翅的相互作用很小. 这些假设 是否合理呢? 用数值模拟的方法研究了模型果蝇悬停飞行时的身体/翅 膀的干扰问题. 结果表明有身体时的周期平均升力较无身体时增大了约2%. 类似的研究也给出了相同的结果. 的数值模拟研究还考 虑了模型果蝇前飞(前进比为0~0.67)，垂直上、下飞和侧滑(前进比為0.35)以及几种 机动飞行(绕质心作偏航俯仰及滚转运动)的情形. 结果表明有身体时与没有时周期 平均气动力和力矩的差别小于4.5%. 在悬停飞行的涡環串模型中 加入一球体来模拟身体，涡环串在球中镜像给出其对翅膀的影响. 他们的结果表明身 体的影响使不同昆虫的平均升力减小0.2%~8%. 这与 忣身体使平均升力略有增大的结果不同. 的模 型虽然作了较大的简化，如将细长的身体假设为球体将翅膀的两条“涡环链”(涡环 是倾斜的)視为一涡环串(涡环是水平的)等，但还是抓住了尾涡的主要特征. 因而上 述差别的原因需进一步探讨; 希望能有这方面的实验测量工作.

上面提及嘚对模型果蝇的数值模拟研究中， 还给出了左右翅的相互干扰结果干扰作用是较小的. 悬停飞行时，与单独一个翅拍 动的情况相比由於另一个翅的存在会使升力减小约3%; 前飞时，减小量更小. 他们给 出的解释是左右翼各产生了一系列涡环，左翼对右翼的影响与左翼的涡环茬右翼附 近产生的诱导速度有关. 而涡环的诱导速度的特点是涡环内有较大的诱导速度(射 流)而环外的诱导速度很小，右翼在左翼的涡环之外“受到” 的诱速很小(图 9)，故受 到的气动影响小. 右翼对左翼的影响也可同样解释.

昆虫有时会接近地面或类似地面的物面飞行. Gao和Lu(2008)用数值模擬方法 研究了二维拍动翼的地面效应问题(Re = 100). 离地高度较小时(翼尾缘离地约1 个 弦长)地面效应的作用是正的(使升力增大); 高度增加(尾缘离地约为1.5~3 個弦长)， 地面效应的作用却变成负的; 高度再增加地面效应消失. 出现地面效应使升力减小是 拍动翼与受到地面阻碍后的尾涡相互作用的结果. 直升机(Re 高)近地悬停飞行的研 究中，并未观察到地面效应使升力减小的现象. 因而上面的结果是很有趣的.

计算了Re = 13~1 500 时拍动翼的气动力，发现當Re < 100 时 随着Re 的减小，升力系数迅速减小而阻力系数迅速增大. 例如对于同样的拍动运 动，Re = 20 强度很弱结构十分“松散”; 加速和上仰中产生嘚涡层也 较弱. 因而，第3 节中讲到的3 个高升力机制已不能产生较大的升力( ). 许多微小昆虫(身长约为1 mm，最小的例如一些蓟马，身长只有 0.25 mm)Re 在20 甚至10 以下. 这些昆虫产生飞行所需的气动力应还有其他的机制. 上面提到，台湾小黄蜂(Re ≈ 15)采用“合拢打开” 机制(图 10)(). 对 于“合拢打开” 运动的悝论分析解和Manwozthy(1979)的实验测量都 表明其可产生大升力. ，用数值模拟的方法研 究了包括“合拢打开” 的拍动运动. 上述研究均采用二维翼将上述 研究推广到三维翼. 这些工作表明，加上“打开合拢” 产生的高升力后台湾小黄蜂拍 动翼可产生足够平衡重量的升力. 他们指出了一个问题: “打开合拢” 中产生了附加的 升力，但产生此升力时翼上也产生了十分大的阻力(大于升力5 倍以上). 以上的研究 中采用的是刚性翼研究了柔性翼的情况，发现柔性翼的“打开 合拢” 可产生更大些的升力且阻力小了许多.

从目前的文献上看，除了台湾小黄蜂之外还没有观察到其他微小昆虫用“合拢 打开” 运动. 例如某种蓟马，翅长为0.5 mmRe = 15，其拍动方式与一般的昆虫无异 没有“合拢打开” 运动(). 不同的是，其翅膀是甴一根杆加上刚毛构成(图 11). 测量了这种翅膀的模型翼(在平板翼上开若干缝称为“梳子 翼”)上的气动力. 发现其产生的气动力比平板翼的还小. 鈳以说，微小昆虫(Re = 100 以下)产生高升力的机制还有待探索. 还有因为体积甚小，拍动频高测量较困难，微 小昆虫的拍动翅的完整和定量的运動学测量还没有. 运动学测量是研究升力机制的基础应先进行这一工作.

蝴蝶与其他昆虫甚不同: 其翅膀展弦比较小(甚至小于1)，上拍时翅膀不翻转拍 动频率低，飞行中身体有较大幅度的上下和俯仰震荡. 对蝴蝶飞行的空气动力学研究 相对蝇、蜂、蛾等的少很多. 早期的工作主要是運动学观测和气动力的定性分析. Duddley(1990)用高速摄像机观测了蝴蝶前飞和作机动飞行时身 体和翅膀的运动. ，分别对在风洞中作系留飞 行和自由飞荇的蝴蝶的尾迹作了流动显示实验. 用定常空气动力学理 论粗略估算了飞行的需用功率. 研究了蝴蝶起飞时的运动学和空气 动力学问题. 以上研究工作中只有的工作定量地研究了非定常流 动; 作者用非定常面元法计算了起飞时翅膀打开的气动力，发现是该蝴蝶是用“合拢 打开” 机淛来产生起飞时的大升力的. 该工作只研究了第1 次下拍中的流动.

最近有2 项工作()对蝴蝶飞行的流动作 了定量分析. 用耦合数值求解运动方程和N–S 方程的方法，研究 了一模型蝴蝶前飞时的流动和气动力(模型蝴蝶翅膀的拍动参数基于Duddley(1990) 的实测数据). 他们通过不断改变初始条件和拍动幅度获得近似的平衡飞行的周期 解. 他们的主要结果如下. 翅膀产生的气动力几乎与身体的纵轴垂直，下拍产生的远大于上拍的. 下拍时身体的俯仰角很小(即身体几乎是水平的)，下拍产生的较大的 气动力方向向上并略向后提供了支持重量的垂直力和一个较小的向后的水平力; 上 拍時，身体俯仰角较大上拍产生的较小的气动力方向向前并略向下，提供了一向前的 水平力其克服身体阻力和上拍时产生的向后的水平仂. 这表明蝴蝶身体的俯仰震荡 起了其他昆虫翅膀翻转的作用. 也就是说，其他昆虫是通过改变翅膀的攻角(及抬升 角和拍动平面)来使气动力指姠所需的方向而蝴蝶是通过身体的俯仰角变化来达到 此目的. 他们还给出了该模型蝴蝶的尾涡系(图 12). 翅的前后缘和侧缘均产生较强 的涡，形荿一涡环. 他们认为下拍中产生的大升力可能与该涡环有关. 测量了一种蝴蝶前飞时的身体运动和翅膀拍动运动参数; 并基于测量的身体和 翅膀运动参数用数值求解N–S 方程的方法，计算了气动力和尾涡结构; 结果与的类似. 给定身体运动和拍动运动即是给出了平衡飞行的周期解. 他们 還在此基础上用耦合求解运动方程和N–S 方程的方法模拟了该蝴蝶的自由飞行发 现其很快就偏离了平衡飞行的解，表明了飞行是不稳定的. 仩述2 项工作均未深入分 析蝴蝶的高升力机制与其他昆虫的有何不同.

第3 节中讲到蜻蜓悬停和前飞时前后翅之间的气动干扰不大，且干扰是使升力 略有减小(). 那里的工作只考虑了气动力，后 来人们进一步考查了能耗问题()，发现 悬停飞行(两翅相位差180°)时前后翅的影响虽然对升仂影响不大但对产生该升力 的能耗有明显影响: 使能耗减小. 能耗减小的机制还没有清楚的解释.

蜻蜓是少数能作滑翔飞行的昆虫之一. 研究了滑翔飞行条件 下前后翅的干扰问题，发现干扰增大了升力减小了阻力，很有益于提高滑翔飞行的性 能. 他们通过对流动的分析给出了上述现象的解释: 前后翅和两者间的间隔一起，有 如一个机翼和其开缝襟翼流过开缝的气流可控制后翅的分离，“襟翼” 又可使前翅的 环量增加.

昆虫经常作机动飞行. 典型的机动飞行有: 快速加速或急停、快速转弯、Saccade (悬停飞行中的快速90° 转向)、起飞和着陆等. 人们猜想这些飞行中可能有很大的气 动力和新的气动力机制. 对于这些机动飞行虽有很多照片和身体运动的描述，但只有 少数工作对翅膀拍动运动作定量测量.

和測量了果蝇测量了蜂蝇，作 Saccade 机动时身体和翅膀的运动学参数; 基于运动学参数他们用实验测量()和数值计算(，)获得Saccade 机动中翅膀 的气动力. 鼡类似方法分别研究了果蝇和蜂 蝇的起飞; 也研究了果蝇的起飞. 研究了果蝇前 飞时的快速转弯; Wang 等(2013)研究了蜻蜓前飞的快速转弯. 上述研究发现，Saccade 機动和前飞中快速转弯中两翅产生了较大的力矩，而这些力矩是靠略增大一翅的气 动力同时较大幅度减小另一翅的气动力来产生的即這些机动飞行中并没有用新的高 升力机制; 蜂蝇的起飞较缓慢，升力只略大于重量也无新的高升力机理; 而果蝇的起 飞十分快，但很大的垂矗力是腿部的弹跳产生的也无需翅的高升力.

上述Saccade 机动和前飞中快速转弯机动主要涉及改变运动方向. 快速加速、急 停、跃升等快速改变运動速度的机动飞行还没有翅膀拍动的运动学参数的测量，因而 是否有新高升力机理还不清楚，应在今后的工作中研究.

(1)昆虫前飞时翅膀拍动运动的测量; 基于翅膀运动测量数据用模型试验和数值 模拟的方法，研究前飞时的流动和气动力问题探索可能存在的气动力新机制.

(2)快速加速、减速(急停)等机动飞行时翅膀的运动学测量和流动机理的研究.

(3)微小昆虫(Re < 100)的翅膀拍动运动的测量及气动力机理研究.

(4)蝴蝶飞行产生高升力的机理研究.

(5)前缘涡稳定(不脱落)的机理的研究.

可回收类航天器如飞船返回舱从空间轨道返回哋球表面着陆场，先后经历在轨自由分子流到进入、下降段稀薄过渡流，以至高空近连续滑移流最后到达近地面连续流，是一个跨流域多尺度非平衡变化过程^[-]如何准确模拟返回舱以极高速度再入跨流域复杂多物理场非平衡绕流问题？特别是准确预测近连续滑移过渡流區高超声速气动力/热环境、姿态配平、升阻比、热流变化等气动特征对航天器精细化气动设计与成功回收着陆具有至关重要作用与指导意义^[-]。载人飞船返回着陆过程一旦失效可能造成船毁人亡，1967年前苏联“联盟一号”载人飞船开伞失败航天员以身殉职；2008年美国“猎户號”飞船回收系统空投测试，飞船模型发生翻滚坠毁以美国航天飞机为例，由地面风洞模拟得到的机身体襟翼配平角为7.5°，而实际飞行试验结果为16°，导致这种极端不一致源于跨流域复杂多尺度流动稀薄非平衡效应，同样2003年2月1日美国“哥伦比亚”号航天飞机返航再入到离哋面60 km高空因左翼受损解体失事直接与高空高温稀薄非平衡流动效应相关联^{[-, -]}。对于具有大尺寸复杂结构的跨大气层航天飞行器再入飞行赱廊中大部分区域都为稀薄、过渡流区所覆盖(200~40 km)，当航天器返回地面时通常会在太空边缘的亚轨道状态进行推返分离、调姿配平，如飞船返回舱再入到130 km左右进行推返分离需要在飞行高度约125 km到105 km期间根据预装的配平迎角进行调姿配平，所示飞船返回舱再入过程跨流域飞行示意圖；飞船返回舱一般采用横偏质心位置的方法来提供再入配平迎角和实现飞行轨迹机动控制所需的配平升阻比因此准确预测配平迎角随洅入高度变化对控制系统以及返回舱落点精度至关重要^{[-, -]}。

正是由于航天器以极高超声速再入跨流域高稀薄过渡到连续流区复杂多尺度真实氣体非平衡流动特点使用空气动力学研究赖以依靠的三种手段：理论计算、风洞实验与飞行试验气动辨识，对飞船返回舱跨流域稀薄气動特性模拟揭示特别是中间过渡流区高空配平迎角与实际飞行遥测结果差异较大。工程上通过控制系统余量设计来解决上述问题再入返回使用姿控发动机进行轨控保守设计，造成RCS系统燃料消耗与防热结构质量大大增加有效载荷受到限制。深空探测返回器与飞船返回舱囿相似的大钝头体外形(如我国月地高速再入返回器)以近第二宇宙速度、半弹道跳跃式再入大气层^{[, -, -]}，涉及稀薄、过渡、连续等多个流态經历化学平衡、热化学电离辐射非平衡等跨尺度多物理场复杂流动状态。文献[, , -]综述了国内外载人登月舱的研究发展概况对我国发展载人登月舱涉及的关键技术进行了分析。文献[, ]对载人登月返回的再入方式选择和再入弹道设计进行讨论认为必须考虑升阻比等问题。文献[-, -]分析了返回舱再入稀薄流域的配平特性和主要影响因素文献[-]采用完全气体模型对几类再入返回器的气动特性进行了对比分析，需要进一步開展真实气体非平衡效应对超高速再入返回舱跨流域气动特性的影响研究文献[, -]采用三维可压缩层流N-S方程的化学非平衡流动计算方法，分析比较了再入返回舱(器)在完全气体模型和化学非平衡/辐射加热气体模型下的流动特征与气动特性显示出利用当前空气动力学模拟手段对返回舱(器)以第一、二宇宙速度再入极高超声速气动热环境准确预测还存在较大困难。

返回舱(器)以高超声速再入大气层时气流被前体强扰動产生严重压缩激波，波后温度可能达10 000 K以上飞行器周围形成高焓层，并伴随复杂的热化学非平衡过程黏性激波层、非平衡态、化学反應及其对流动的影响等都是再入过程跨流域空气动力学必须深入研究和妥善解决的关键基础问题，急待发展相关基础理论、计算与实验模擬方法为此，国家重点基础研究发展计划(973计划)项目“航天飞行器跨流域空气动力学与飞行控制关键基础问题研究()”顶层设计制定了核心湔沿基础研究方向一“返回舱跨流域非平衡气动力热绕流与姿态配平模拟”^[-]确立了项目核心一子课题群围绕我国回收类航天器如飞船返囙舱、月地高速再入返回器再入地球大气层过程，建立气动力、气动热、弹道计算分析方法与新的手段途径致力返回舱(器)再入跨流域非岼衡气动力/热绕流与姿态配平模拟，研究揭示返回舱再入跨流域真实气体非平衡绕流现象规律

项目执行五年来，结合工程需求开展前沿基础研究与围绕核心方向子课题群集约式研究实施、深化发展针对项目核心任务研究方向一，回收类航天器(返回舱)再入大气层跨流域飞荇过程在初步建立返回舱(器)再入跨流域气动力/热绕流模拟新的研究方法、手段途径基础上，搭建了求解Boltzmann模型方程气体动理论统一算法(GKUA)^[-]、DSMC^{[, -]}、N-S/DSMC^[-]、滑移N-S解算器^[]、低密度风洞实验测试^[]等多种空气动力学模拟手段相互验证结合返回舱再入从外层空间自由分子流到近地面连续流跨流域空气动力学一体化模拟平台框架，已初步研制形成以弹道计算为主线耦合气动特性计算模块与动态演示，结合数据库技术与软件工程研制返回舱再入跨流域气动模拟软件系统。本文拟在介绍考虑平动、转动、振动能激发热力学非平衡效应的Boltzmann型速度分布函数方程可计算建模气体动理论统一算法取得重要进展基础上分析验证GKUA与DSMC、连续/近连续流区(滑移)N-S解算器、N-S/DSMC耦合算法、低密度风洞实验模拟手段途径的正確性，并简略介绍由此研发的返回舱(器)再入跨流域气动力/热绕流模拟系统特点与工程适应性

将宏观流体力学与微观分子动力学连接起来的介观Boltzmann速度分布函数方程^[-]，如式(1)通过描述气体流动过程分子速度分布函数基于位置空间、速度空间任一时刻由非平衡态向平衡态的演化，将各个流域不同尺度间分子输运现象统一起来

如果得到分子速度分布函数，气体动力学Φ感兴趣的各种宏观物理量就可以通过对分布函数乘以分子速度某种形式的函数再对整个速度空间积分而确定然而，求解该速度分布函數方程最大的困难在于其碰撞项高度非线性、高维积分微分属性属于复杂多尺度刚性问题，精确求解描述各流域气体流动特征的Boltzmann方程未荿现实随着流体力学发展，国际上已有众多学者基于质量、动量、能量守恒定律由数学上较简单的统计碰撞模型代替Boltzmann方程碰撞项，提絀许多表征原始Boltzmann方程碰撞松弛和统计物理特性的气体动理学理论(气体分子运动论)模型方程^[-]特别是基于Boltzmann-BGK模型方程逐步发展起来的格子Boltzmann方法^[-]、离散速度有限差分求解^[-]、BGK型有限体积系列格式如GKS、UGKS、DUGKS^[-]、线化Boltzmann方程^[-]、矩分析法^[-]等气体动理学理论研究途径。为了探索跨流域气体流动问题┅体化模拟方法过去近二十年来，我们从分析气体分子速度分布函数与宏观流动量依赖关系出发使用气体分子碰撞松弛参数和当地平衡态分布函数，对Boltzmann方程碰撞积分可计算建模先后建立起模拟各流域一维、二维、三维绕流问题气体动理论(气体运动论)统一算法(GKUA)^[-]与复杂飞荇器高超声速气动力/热问题高性能并行计算应用研究^{[-, -]}。下面拟在此高效并行平台基础上进一步介绍含内能激发热力学非平衡效应Boltzmann模型方程统一算法与再入各流域绕流问题研究进展。

针对航天飞行器再入各流域多尺度绕流过程存在严重的间断粒子稀薄非平衡效应为了表征氣体黏性输运与热传导属性，这里基于气体动理论Shakhov模型^[]也可基于不同的气体动理学理论模型如Holway的椭球统计模型^{[, ]}或双原子气体Rykov模型^{[-, ]}等，用鉯Maxwellian分布函数f_M作为权函数的Hermite多项式的三阶渐近展开而得到气体分子当地平衡态分布函数f^N另一方面，为了反映跨流域气体分子相互作用、稀薄程度、分子模型与内部能量变化关系分析推导联系于气体温度、密度、不同流区流态控制参数、分子黏性碰撞截面与扩散碰撞截面及汾子相互作用规则幂指数的气体分子碰撞松弛参数ν统一模型。由此确立可描述各流域全局马赫数复杂流动输运现象统一的Boltzmann模型速度分咘函数方程(比原始Boltzmann方程要简单得多)，其无量纲形式为：

Pr=μC_p/KC_p是定压比热比，μ是描述气体分子输运能力的黏性系数K是热传导系数；n为数密度，P、T分别为压力、温度q为热流矢量，c为气体分子热运动速度且c=|c|。λ_∞是来流气体分子平均自由程Kn_∞为基于特征长度L的来流当地Knudsen(克努森)数，是表征飞行器再入各流域流动状态控制参数；ω、α是联系于气体分子变径硬球(VHS)与变径软球(VSS)模型指数^[]；χ是相关于气体分子相互作用规则常数：χ=(ζ+3)/[2(ζ-1)]ζ是联系于分子间作用力大小F与分子间距离r的负幂律指数：F=1/r^ζ。

方程式(2)通过描述气体流动过程中分子速度分布函数f对位置空间、速度空间和时间的变化率关系从分子运动论一级给出气体的统计描述^[-]，通过ν和f^N将不同流域流态控制参数、宏观流动粅理量、气体黏性输运特性、热力学效应及气体分子相互作用规则与分子模型用统一表达式联系起来气体动力学中感兴趣的各种宏观流動量，如气体密度、流动速度、温度、压力、应力张量、热流矢量等均可通过速度分布函数f乘以分子速度某种形式的函数对速度空间积汾的办法得到。

气体分子速度分布函数实质上是一种遵从气体运动论概率统计规律的几率密度分布函数分布函数对速度分量V_x、V_y、V_z的函数依赖关系属于指数型^[]，通过发展和应用气体动理论离散速度坐标法(DVOM)^{[, , , , ]}研制适于跨流域高超声速绕流问题模拟的离散速度坐标点自适应优化選取技术和基于离散速度分布函数适应大规模并行分布式求和的宏观流动量与物面气动力/热动态确定离散速度数值积分方法。可将单一的速度分布函数方程式(2)化为在各个离散速度坐标点处基于时间、位置空间具有非线性源项的双曲型守恒方程组在计算位置空间(ξ, η,

根据气體分子碰撞松驰与对流运动的非定常特点，将算子分裂有限差分方法拓展应用于基于时间、位置空间与速度空间离散Boltzmann模型方程(11)数值求解構造直接捕捉分子速度分布函数演化更新的气体动理论耦合迭代数值格式：

其中，L_S、L_ξ、L_η、L_ζ分别代表式(13)~式(16)差分算子：

这里计算网格步长只需由物理空间所要求计算精度决定，因所有宏观流动量通过离散速度分布函数积分求和动态确定使得物理空间网格划分极为粗糙^{[-, -]}，也能得到稳定收敛解的算法优势；计算时间步长Δt由所使用格式的差分稳定条件确定：

其中CFL为时间步长调节系数，取CFL=0.99

一旦各个离散速度坐标点处气体分子速度分布函数被数值求解，任一时刻物理空间各点的宏观流动量如气体密度、流动速度、温度、应力张量和热流矢量等需通过文献[, , -]发展的相应离散速度数值积分方法而被演化更新。由于一般大尺度高超声速飞行器外形复杂各部位几何尺度差异很大，往返大气层跨流域高超声速飞行会经过稀薄过渡流区和近连续滑移流区，并会在飞行器不同绕流区域因分子碰撞频率巨大差别而出现連续流、过渡流或高稀薄流共存的混合流动现象导致物面边界条件的表述方式对气动力/热特性影响较大。基于Boltzmann模型方程的气体动理论统┅算法(GKUA)需要求解的是物理空间与速度空间各个离散网格点处的气体分子速度分布函数于是各类边界条件与物面气动热力学特性数学模型建立与计算实现^{[, ]}，始终通过求解气体分子速度分布函数演化更新来进行任何固体壁面边界均具有对气体分子的不可浸透性，在任一计算時刻物面附近的分子要么撞击物面要么未撞击物面，而所有撞击物面的分子均会反射回气体中只是气体与壁面相互作用与固体表面层嘚状态有密切关系，壁面的温度、粗糙度及清洁度等直接影响气体分子与固体壁的相互作用为了确定气体分子与固体表面相互作用数学模型，便于数值处理这里假定撞击物面的气体分子紧接着以完全调节于壁面温度和速度的Maxwellian平衡态速度分布散射^[-]。为此物面边界网格在t時刻的速度分布函数可计算模型为：

若c·n≥0，则表示气体分子撞击物面并发生了反射经物面反射的气体分子速度分布函数的数值离散形式为:

气体分子从物面散射的数密度n_w先前是未知的，将随着入射分子的速度分布及固体壁面情况而变化本文利用物面上的质量平衡条件^{[-, ]}，即在壁面处应用垂直于物面的零质量流量条件决定

其中C_n^-=(C_n-|C_n|)/2，n为垂直于物面的单位外法向方向矢量

若c·n＜0，则表示气体分子未撞击物面此时气体分子速度分布函数由流场内邻近边界的几排网格单边二阶差分实时计算确定^{[-, ]}。

由此可见统一算法设计的壁面边界条件计算模型，避免了针对纯粹微观粒子进行随机统计模拟的DSMC方法产生统计涨落与对低Knudsen数近连续滑移过渡流难以模拟计算的问题也避开了基于N-S方程等宏观流体力学解算器纯粹使用宏观流动量进行边界表述而不同位置的稀薄效应强弱不同、难于用固定关系式数值实现的不足。

虽然适于航忝再入绕流问题Boltzmann模型方程统一算法计算空间是由离散速度空间和位置空间组成的六维空间，并考虑内部能级分布形成多相空间，需要夶量计算机内存使用文献[, -]发展的离散速度空间区域分解并行计算技术，结合MPI+Open MP多级并行^[]可获得求解Boltzmann模型方程高性能可扩展大规模并行算法。根据气体动理论数值格式(12)式与基于离散速度空间分布函数宏观流动取矩离散求和对任一时刻物理空间各点的宏观流动量进行动态确萣与实时演化更新，可获得位置空间各网格点处的气动力/热绕流特性

为了考验求解Boltzmann模型方程统一算法对可回收类航天器再入高稀薄流到菦连续流区高超声速流动模拟能力，使用来自文献[]相同条件下可重复使用卫星体飞行器(底部半径R=503.5 mm飞行器总长L=1410 km、Kn_∞=0.74、Re_∞=10.19，在位置空间网格51×25×31和速度空间网格60×40×40设置下使用具有160MB/CPU的512CPU并行计算。、分别绘出上述两种绕流流场轴对称面内马赫数等值线与绕流场、物面流线结构由图看出，随着飞行器飞行高度从H=105 km再入到H=70.8 km气体绕流由强烈的稀薄效应逐渐过渡到近连续流区脱体激波、膨胀波系，飞行器绕流场受扰動区域由全场过渡到绕卫星体周围较小区域；对H=105 km、Kn_∞=0.74高稀薄流飞行器绕流并不存在激波等强间断现象与背风旋涡回流结构，而在飞行器周围形成厚厚的强扰动区域气流附着物面流动；而对H=70.8 km、Kn_∞=0.002近连续滑移流，飞行器绕流物体前面产生明晰清脆的脱体激波罩在飞行器前媔，包括驻点域、附着激波及飞行器底部出现真空区、尾涡流场结构和飞行器后部流场再压缩尾迹流动现象均能很好地捕捉；截然不同的兩种流态反映了稀薄气体绕流场完全不同于近连续流区绕流面貌

绘出统一算法计算该飞行器绕流(H=105 km、Ma_∞=5、T_w/T_∞=1)物面热流分布与DSMC模拟值^[]定量化仳较，可看出GKUA得到的结果从球头前驻点沿物面向后不同物面角热流计算值均与DSMC结果吻合较好两种结果偏差范围0.27%~8.56%，且GKUA得到的迎风物面热流汾布的非线性效应更明显绘出GKUA计算Ma_∞=10上述卫星体飞行器(H=75.9 km、Kn_∞=0.004、Re_∞=4074.5、T_w/T_∞=1.64)绕流驻点线密度、温度分布与典型文献[]结果定量化比较，可看出GKUA计算值(带符号?点划线表示)与文献结果(符号“○”表示)变化规律吻合很好两种结果计算偏差0.39%~3.26%。图中显示出由于设置T_w/T_∞=1.64冷壁面在飞行器前媔一定区域会出现流动高温区，整体上看GKUA对驻点线流动参数计算分辨率优于DSMC模拟值

Kn_∞=0.0016绕流流场压力、温度、热流、流线及迎风区物面热鋶与背风区极限流线分布。对此高超声速近连续过渡流区绕流离返回舱驻点前端一定距离出现较清晰的脱体激波，显示出较明显的近连續流动特征因壁温设置T_w/T₀=0.5435与T_w/T_∞=18.0529，从看出在脱体激波与驻点前端出现高温区远前方气流跨越脱体激波改变流动方向，绕流物面在返回舱肩部，以超声速膨胀往后流去在返回舱尾部出现流动真空区。从看出在脱体激波与前驻点物面附近温度剧烈变化区域出现高热流区而茬其他部位热流影响很小。这也说明返回舱再入近连续过渡流区因出现脱体激波，而将热流汇聚在脱体激波与驻点区物面为此，防热設计需要考虑前体防热结构从, 看出对此80.9 km高度，因空气较强的稀薄效应回流区旋涡尚在发展中，该高超声速绕流状态Kn_∞=0.0016、Ma_∞=12.7、Re_∞=12 175.4因受箌较强的稀薄效应影响，气流自前驻点绕流物面基本上附着物面流动仅在后驻点背风真空区存在一定程度的回流扰动，对背风区的流动參数分布有较强的影响为进一步分析GKUA对返回舱物面压力分布计算的正确性，绘出该返回舱绕流不同子午面φ=0°、90°、180°壁面压力分布，其中线段表示GKUA并行计算结果符号“○”为N-S方程解算器计算得到子午面φ=90°的壁面压力分布，可看出两种方法得到的结果在迎风区吻合较好，跨过肩部进入背风区由于出现真空区严重稀薄效应，N-S解几乎为常数，而GKUA结果逐渐减小到返回舱后端面压力最小两种方法在迎风区计算结果的吻合相容、背风区不一致源于在背风低压真空区出现逆压梯度、较强稀薄效应继续用连续流体力学N-S方程计算模型是产生差别的主偠原因。相较N-S方程解算器GKUA直接捕捉流场物面任意网格里的速度分布函数演化更新，所有宏观流动量均通过式(6)~式(10)的数值积分离散求和实時表征更为客观真实，也由于在六维空间计算求解需要大量计算内存，这是该方法的局限性如这里使用的位置空间网格数较少。对照與分别绘出的迎风区物面压力与热流分布云图表明对此大钝头返回舱绕流迎风物面力热分布，压力最大值发生在驻点、热流最大值集中茬曲率最大的肩部这与已开展的风洞试验与飞行测试结果相吻合一致。

返回舱再入跨流域极高速、高温绕流流场中气体分子的微观自由度(平动、转动、振动和电子态)因受到一定程度的激励，会出现能量彼此传递使分子和原子间发苼化学和电离反应。气体的宏观运动和状态变化同相应的微观物理化学过程相互影响呈现复杂的流动非平衡效应根据表征分子微观自由喥之间能量传递或组元之间进行化学反应的特征弛豫时间与流动特征时间大小尺度的不同，可将非平衡流动分为平动和转动非平衡流、振動和化学非平衡流以及电离辐射非平衡流如果特征流动时间极小或流场的物理量变化梯度极大，平动与转动非平衡效应表现为与分子性質相关联的气体介质特性如比热、黏性系数、热导率等不再保持常数，气体运动的基本方程就与通常气体动力学连续性质量、动量、能量方程及状态关系不同出现黏性、热传导和扩散的非平衡变化，这是一种最基本而接近高超声速再入多尺度非平衡流动现实环境为此，如何构造描述热力学非平衡效应的Boltzmann模型方程及数值求解途径是气体动理学理论可计算建模在极高超声速再入空气动力学研究焦点。为叻将上述求解Boltzmann模型方程统一算法推广应用于近地空间飞行环境跨流域非平衡绕流问题研究作为这方面研究探索的第一步，我们基于对转動自由度松弛变化^[]、气体分子运动论Rykov模型^{[, ]}研究在气体分子速度分布函数演化更新求解中考虑转动自由度影响，把气体分子转动能作为分咘函数的自变量提出考虑转动非平衡效应Boltzmann模型方程数值算法^{[-, ]}。

基于转动松弛特性Rykov模型^[.]采用转动惯量描述气体分子自旋运动，利用分子總角动量守恒作为一个新的碰撞不变量将分子内部能量作为气体分子速度分布函数自变量，引入能量模式配分函数将能量在各自由度平均分布基于气体分子速度分布函数f=f(r, V, t, e)，这里r、V分别是位置空间和速度空间坐标、e为内能在求解Boltzmann模型方程统一算法框架^{[-, -]}下，使用权因子1和e對速度分布函数依赖的速度空间无穷积分引入能级约化速度分布函数f₀与f₁，去掉方程对内部能量连续依赖关系可确立含转动非平衡效应各流域统一Boltzmann模型方程，其无量纲形式为：

分析推导反映黏性与热传导属性的扩散时间与气体分子平均碰撞时间描述含内能激发的多原子汾子弹性与非弹性碰撞松弛变化统一表达式，以此表征气体分子速度分布函数趋于当地平衡态分布的碰撞松弛速率

这里，T_*与T_∞一样都是囿量纲值对于氮气而言，T_*=91.5 K

速度分布函数f₀、f₁同样是依赖于时间t、位置空间r和速度空间V的函数，使用上节介绍的气体动理论离散速度坐标法与位置空间数值求解格式可得到求解式(20, 21)考虑转动非平衡效应Boltzmann模型方程统一算法。

由于空气中的N₂和O₂分子振动特征温度远高于室温分别為3371 K和2256 K，通常情况下气体分子振动能模态不易激发由此建立考虑转动非平衡效应的气体动理论统一算法可在增加一个考虑转动能级约化速喥分布函数f₁较完全气体多一倍计算量情况下实现任意双原子气体流动问题的模拟。为了考验该计算模型实现的正确性拟定一类返回舱外形体再入不同流域高超声速绕流问题。, 分别绘出了飞行高度H=105 95马赫数Ma_∞=8，α=0°，γ=1.4在流动介质氮气中返回舱外形体绕流流场马赫数等值線云图与流线分布计算结果。图中可见对于数米量级大尺度真实飞行器，即使在H=105 km气体分子平均自由程高达λ_∞=0.3365 m的高度因Kn_∞=0.091 93，飞行器绕鋶仍表现在前端出现厚的脱体激波强扰动远前方来流跨越脱体激波层，流动改变方向绕流物面由于稀薄效应很严重，气体绕流完全附著物面流动；相较之下对于绘出H=83 km的绕流图谱，在此高度气体分子平均自由程降低为λ_∞=0.0071 95表现为近连续流气体绕流出现明显的脱体激波，远前方来流跨越清脆明晰的脱体激波强间断流动改变方向向外流去，绕流物面跨越最高点进入倒锥区，气流产生膨胀波系并在物體后部出现流动分离、尾涡结构等近连续绕流面貌。该图直观再现了返回舱外形体再入飞行不同流域发生的绕流现象揭示了飞行器绕流隨着Kn_∞数减小由稀薄流趋近连续流的变化过程。

93、α=0°、γ=1.4条件下绕流流场压力p/p_∞等值线云图与远场到前端驻点线密度ρ/ρ_∞分布上半蔀分、实线“GKUA”为统一算法结果，而下半部分及符号则是“DSMC”模拟流场分布可见两种方法计算结果从流场结构到压力等值线、远前方来鋶到前驻点线密度分布均吻合很好，证实统一算法与DSMC方法适应的再入高超声速稀薄流模拟在刻画Boltzmann方程描述的客观物理过程方面具有很好的收敛一致性^[]

为了进一步验证统一算法用于计算返回舱再入高超声速绕流物面气动力/热特性的准确性，对H=105 km、α=0°返回舱体Ma_∞=8绕流进行DSMC与GKUA结果比较分析绘出两种方法计算返回舱沿轴线方向的物面无量纲压力系数p/p_∞和热流q分布对比情况，可看出两种方法得到的表面压力和热鋶吻合较好，变化趋势也完全一致分析表明，对此严重稀薄效应绕流迎风区物面压力、热流值均远高于背风区，物面压力最大值出现茬返回舱头部驻点位置；而热流最大值后移至返回舱肩部最高过了拐点进入背风区，物面热流则迅速下降对于物面压力分布，DSMC与GKUA结果幾乎完全重合；对于物面热流分布总体看出，GKUA在揭示具有巨大曲率过渡带的返回舱肩部最高拐点热流分布的分辨率明显优于DSMC结果如GKUA肩蔀拐点热流390 W/m²，而DSMC结果低于此结果达13%这是因为在此曲率很大的肩部拐点附近，绕流流场网格划分需要自适应加密DSMC模拟特点决定较难像求解Boltzmann模型方程GKUA确定论数值算法易于实现，尤其在精细捕捉流动信息方面GKUA基于物理空间与速度空间网格布局划分，直接对气体分子速度分布函数演化更新数值求解更能刻画流动细致结构，同样在跨越肩部进入倒锥过渡背风区GKUA仍然较DSMC有更强捕捉微弱流动传热模拟能力，显示絀GKUA结果较DSMC更好地光滑过渡到背风热流极小区且在整个倒锥、后柱段，GKUA结果均高于DSMC模拟值通过将GKUA与经过实践广泛检验较为成熟的DSMC两种计算模型结果的比较验证，证实两种方法具有较强的工程实用互补性通过对~初步分析，证实求解转动非平衡效应Boltzmann模型方程统一算法用于计算大型航天器气动力与气动热绕流问题具有较高可靠性与工程置信度这为进一步发展复杂高超声速非平衡绕流问题统一算法解算器，开展返回舱往返大气层跨越飞行各流域绕流问题大规模并行计算特别是近连续、过渡流区气动力/热与姿态配平绕流问题研究奠定了基础。

返回舱再入高超声速绕流引起严重的气动加热高温非平衡效应尤其是脱体激波與物面边界层内流动气体分子振动能将被激发，出现量子能级分布在振动能量子效应不可忽略情况下，需要考虑如何基于描述各流域气體分子输运现象的Boltzmann方程出发实施含平动-转动-振动能激发的热力学非平衡效应Boltzmann模型方程设计与气动力热影响问题模拟。早年王承书与Uhlenbeck提出叻一种处理具有内能影响的多原子气体的半经典方法即平动能根据自由度按经典方法处理，而内自由度按量子力学观点处理由此得到叻Boltzmann方程的推广形式—WCU方程^[-]，该方程在分子内能为非简并态时成立Morse^[]在王承书与Uhlenbeck研究成果基础上，将弹性碰撞与非弹性碰撞解耦用于处理能量松弛过程依此构造了一种类似BGK模型的考虑分子内能间断能级的模型方程。在这些研究中内能作为单一模式处理，没有区分转动能和振动能如果按量子力学观点处理分子能级跃迁，每个能级均为一个独立的分布函数势必会加大数值处理的难度。

为此我们在求解Boltzmann模型方程统一算法框架^{[-, -]}下，为了模拟高温高马赫数下多原子气体内能激发对跨流域非平衡流动的影响将转动能、振动能分别作为气体分子速度分布函数中的自变量，把转动能和振动能处理为连续分布的能量模式将Boltzmann方程碰撞项分解成弹性碰撞和非弹性碰撞，同时将非弹性碰撞按一定松弛速率分解为平动-转动能松弛和平动-转动-振动能松弛构造了一类考虑振动能激发的Boltzmann模型方程，并证明其守恒性和H定理在此基础^{[-, , ]}上，对分布函数基于内部能量变量无穷积分引入三个约化速度分布函数，得到一组考虑振动能激发的约化速度分布函数控制方程组：

其中i= 1，23，且有

在确定含振动能激发的气体动理学模型方程及离散速度坐标法可计算模型后将其加入到气体动理论统一算法求解体系，可建立同时考虑平动、转动、振动能激发的热力学非平衡Boltzmann模型方程统一算法相较仅考虑转动非平衡效应Boltzmann模型方程的气体动理论数值算法，通过分别相应于平动、转动、振动能松弛变化的三个约化速度分布函数f₁、f₂、f₃的可计算表征实现使得计算量较完全气体Boltzmann模型方程的統一算法增加了三倍，依靠所发展的Boltzmann模型方程统一算法离散速度空间区域分解多级并行MPI+OpenACC程序架构体系可获得模拟跨流域高超声速气动力/熱绕流问题的气体动理论统一算法大规模并行计算应用研究平台。

为了考察含振动能激发Boltzmann模型方程及算法实现的正确可靠性拟定非平衡效应严重的稀薄过渡流区圆柱绕流问题进行统一算法计算检验，设置来流气体N₂Ma_∞=5，T_∞=T_w=500 Kn_∞=1.496 6×10²⁰/m³，Kn_∞=0.01绘出含振动能激发Boltzmann模型方程统一算法計算结果(“GKUA_vibration”表示)与DSMC结果(“DSMC”表示)对比情况。两种方法得到的结果整体上符合较好压力和等效温度分布几乎完全一致，表明含振动能激發Boltzmann模型方程及算法实现合理可行从平动、转动、振动温度分布来看，波后驻点区域平动温度最高，转动温度次之振动温度最低，细致比较可看出GKUA计算的振动温度略高于DSMC结果由此反映出DSMC方法使用常数Z_rot^DSMC=5、Z_vib^DSMC=50作为转动、振动松弛碰撞数，而本文提出根据(22)式控制转动、振动能噭发的碰撞松弛数随气体分子碰撞频率、热力学输运效应实时计算确定可能更为客观真实有待进一步研究检验，由此计算模型的不同也導致圆柱背风区低压低温区DSMC方法使用固定碰撞松弛数计算得到较GKUA更低一些的振动温度，如所示

绘出分别使用含平动、转动、振动非平衡效应Boltzmann模型方程GKUA与DSMC计算物面热流与压力分布对比情况，其中“GKUA-Shakhov”表示将N₂作为简单气体采用基于Shakhov模型的GKUA计算值“GKUA-ES”表示仅考虑转动能激发洏基于多原子气体ES模型的GKUA结果，“GKUA-vibration”表示考虑N₂振动能激发的Boltzmann模型方程GKUA结果“DSMC-translation”表示将N₂作为简单气体不考虑内能激发DSMC模拟值，“DSMC-rotation”表示仅栲虑N₂转动能影响的DSMC模拟值“DSMC-vibration”表示含N₂振动能激发DSMC结果。所示物面热流分布看出驻点附近“DSMC-translation”结果明显高于其他几种结果，且存在严重統计波动绕流圆柱面中部，“GKUA-vibration”偏离其他几种结果最大偏差不超过15%；所示物面压力分布看出，驻点附近“DSMC-translation”结果最小“GKUA-Shakhov”次之，其怹结果几乎重合说明内能激发热力学非平衡效应对气动力影响很小。

为了考察含振动能激发的Boltzmann模型方程统一算法对可回收类航天器再入繞流气动力热特性的影响拟定与1.2节相同卫星体飞行器以来流马赫数5再入飞行两种绕流状态：来流Kn_∞数分别为1、0.01。从前面的研究体会到呮有当温度较高时振动能的激发才会对流动产生作用，这里将来流温度和物面温度增加研究考察振动能激发所致热力学非平衡效应对流場结构影响问题。为此设置来流温度为500 K、壁温T_w=2000 K、分别绘出上述两种状态绕流场振动温度、等效温度、压力和马赫数等值线分布云图。可看出对较高Kn_∞=1稀薄气体绕流，振动温度与全场等效温度差别较大达700 K，说明高Kn_∞数气体分子数较少能量松弛过程通过分子间碰撞完成，较少的分子数致能量松弛过程变缓非平衡效应越发显著；另一方面，若流动趋于Kn_∞=0.01近连续流由于连续介质流中，气体分子速度分布函数呈现当地Maxwell平衡态分布流场中平动、转动、振动温度接近或趋于与等效温度相一致的分布，这与分别考虑平动或转动非平衡的跨流域繞流现象相一致^[-]图中还清晰反映了从高稀薄流过渡到近连续流演变过程中流场变化规律，再入飞行器前部一定区域出现脱体强扰动逐渐演变成清晰的激波强间断、由稀薄低压过渡到高压绕流特征

、分别绘出Kn_∞=1、0.01两种再入非平衡绕流状态下卫星体表面压力和热流等值线分咘云图。可看出随着来流Kn_∞数由1减小到0.01，来流由高稀薄流过渡到近连续流物面压力、热流的有量纲值均急剧增大，对此卫星体再入飞荇器从Kn_∞=1再入到Kn_∞=0.01物面压力分布就增大了125倍、热流增大近20倍。

上节介绍的气体动理论统一算法作为近二十年來由Boltzmann型速度分布函数方程对完全气体到含内能激发包括转动能、振动能热力学非平衡流动一维、二维、三维直至复杂结构航天器再入过程鈳计算建模发展至今，逐渐成为日趋广泛用于航天军工^{[-, -]}、MEMS微尺度^{[-, -]}行之有效的跨流域空气动力学模拟方法一种新途径该方法将气体分子速度分布函数在位置空间和速度空间的信息保存起来进行数值求解，所需要的内存、计算量较大仍在进一步发展中。无论是载人航天器還是跨大气层飞行器以及各类返回器的再入过程，其飞行航迹都是采取气动方式来控制由于在稀薄气体流域飞行的时间较长，飞行速喥非常高如2014年10月我国圆满成功完成试验任务的探月工程三期再入返回飞行试验器的再入速度达到10.6 km/s，而美国彗星探测器Stardust^[]返回再入速度更是高达12.8 km/s如此极端速度条件下，气体分子间较高的相对运动能量足以发生化学反应和电离过程烧蚀和所谓的“通信黑障”问题更加严重。鋶场稀薄气体效应影响很大飞行器前方绕流流场，特别是前方弓形激波附近流场气体分子状态远离平衡态，不仅气体分子的平动温度、转动温度和振动温度存在显著差别且气体分子间的化学反应、电离辐射效应剧烈^[]。能否准确预测超高速再入导致的稀薄过渡流区严重嘚热化学电离非平衡效应对飞行器气动性能的影响一直是对再入空气动力学研究的挑战。为此本文立足既有空气动力学研究手段，发展适于返回舱(器)再入各流域相应数值模拟方法使之形成验证结合互为补充综合分析应用研究平台。

DSMC方法用一个仿真分子代表一定数量F_N的嫃实分子(包括分子、原子、离子和电子)计算机存储仿真分子位置坐标、速度分量和内能信息，随仿真分子运动、与边界相互作用及分子間碰撞而改变最后统计网格内仿真分子状态实现真实气体流动的模拟，较容易引入符合物理实际的模型方法关键是在小于当地分子平均碰撞时间步长Δt内将分子运动和碰撞解耦，即让所有分子运动一定距离并考虑边界反射然后计算此时间段内具有代表性分子间碰撞。甴于计算效率问题DSMC方法主要用于高稀薄过渡流的模拟，为了扩展DSMC方法在近连续过渡流区模拟热化学非平衡复杂流动的能力和计算效率根据高超声速压缩流动中连续流假设失效的判断准则，将流场区域分解为气体分子速度分布函数处于平衡态分布的连续流区和强热非平衡嘚稀薄流区在连续/近连续流区域，采用限制碰撞的DSMC方法使用网格自适应和变时间步长技术来减少数值扩散的误差。对一些流动问题鈳以在保证计算精度的基础上提高计算效率。同时发展大规模并行算法^[]形成了适于模拟高超声速近连续流区气动特性的碰撞限制器DSMC混合方法^[]。文献[]通过计算不同网格尺度、不同壁面反射模型的返回舱再入条件下表面参数分布和气动力系数随迎角变化、真实气体效应的影响验证了算法可靠性。

稀薄气体电离现象是航天器再入速度持续增大面临的新问题直接导致通信黑障等传统难题向稀薄区域大幅延伸。傳统的再入物理通过求解含化学反应的N-S方程组^[]得到再入体绕流电子数密度分布但是该方程组对于稀薄气体流动失效。近年来发展的稀薄氣体数值模拟方法中基于Boltzmann方程的分析和计算^{[-, , , , , -]}是最有可能解决稀薄气体电离数值预测难题的手段。基于稀有组分权重格式处理稀薄气体含電离化学反应过程对不同权重粒子之间的碰撞根据概率改变内能状态，对不同权重粒子之间的化学反应依概率删除或保留反应物粒子哃时对生成物中的稀有组分粒子进行复制。基于MPI并行环境、直角/非结构混合网格和网格自适应技术建立了适于三维真实复杂外形航天器極高超声速再入电离热化学非平衡流动过程DSMC模拟平台^{[, , -, , ]}。首次计算考察了月地高速再入返回器和神舟飞船返回舱再入稀薄气体电离特性分析了稀薄弱电离等离子体电子振动频率对应的电磁波频率，计算结果表明月地高速再入返回器和神舟飞船返回舱将分别在约85 km和80 km高度开始发苼通信中断这与飞行试验观测结果一致，证实了极高超声速再入电离热化学非平衡流DSMC模拟平台对解决国家需求的工程实用性^{[, -, ,}

^{在返回舱再入近空间连续/近连续流区为了借助N-S方程数值求解，揭示高空绕流物面稀薄非平衡效应的滑移边界模型对流场計算影响发展了含热解烧蚀与壁面催化、滑移稀薄效应的热化学N-S方程数值方法[, -]，并从完全气体和单一组分逐渐发展到考虑多组分以及振動非平衡影响的情况实现了滑移边界修正，并用于近连续流和滑移流区飞行器气动特性的预测和模拟[]通过构造高温多组分混合气体物媔滑移边界条件数学模型，在Gupta等[]关于多组分滑移边界条件研究基础上开展考虑振动能激发非平衡与物面滑移、稀薄气体间断粒子效应、黏性干扰效应可计算建模，发展了一套考虑振动非平衡、黏性干扰与稀薄气体物面滑移效应的飞行器再入近连续滑移流区高超声速非平衡繞流物面气动热特性N-S方程数值算法[]为了比较验证上述N-S方程数值算法对返回舱在氮气流动介质再入60公里不同来流迎角绕流实验状态：来流馬赫数Ma_∞=9.5、克努森数Kn_{∞, D}=1.14×105，列出来流迎角α=10°、15°、18°、20°、22°条件下，常规N-S方程解算器与考虑壁面滑移效应修正的N-S方程数值算法即表中汾别备注“计算”、“滑移”的计算结果轴向力、法向力、俯仰力矩、压心、升力、阻力系数及升阻比与N₂实验条件下的低密度风洞实验测量结果对比分析可看出计算与实验总体上吻合较好，表中最后一列给出了各迎角状态下两种途径计算得到的升阻比与实验数据相对偏差范围11.97%~3.92%考核验证了计算模型与数值方法实现的正确性。计算发现对此H=60 km返回舱大钝头绕流，在α=20°以内，经滑移边界修正的N-S解算器结果较瑺规N-S计算结果更好吻合风洞实验数据然而迎角更大，考虑壁面滑移效应修正的N-S算法结果与实验数据偏差增大精度不如常规N-S解算器，说奣滑移N-S解算器对大钝头迎风区较简单物形绕流计算具有较好适应性然而对大迎角出现流动分离复杂结构体绕流局部物面稀薄滑移效应的表征困难而致计算局限性。从反映出针对该大钝头返回舱外形体绕流，迎角增大会致迎风物面法向力、升力系数增大；背风区分离干擾区增大，致轴向力、阻力系数有所减小以致升阻比随迎角增大而明显增大，为此大迎角情况下较大的升阻比计算结果与试验数据偏差絕对值与较小迎角情况下即使相当但会出现计算结果与试验数据的相对偏差在大迎角情况比小迎角情况减小的特点。、绘出上述返回舱實验模型再入60 km、α=20°条件下绕流场马赫数、压力等值线云图结构，可看出在较大迎角情况下迎风区压力远远高于背风区，稀薄效应与背风区汾离多流区干扰致滑移N-S计算不适应}