N-S 方程01 N-S01 宏观模型是把气体当作连续介质，关于它的描述是以流动参数，如速度、密度、压力和温度在空间和时间的变化

這是 http://202.112.143.25:8001/xwlw/document?RecordNo=195&ColumnName=FREETEXT_URL&MultiNo=0&issource=yes&type=bin&channelid=65195 的 HTML 檔。
G o o g l e 在網路漫遊時會自動將檔案轉換成 HTML 網頁。

 

Page 1

中图分类号：V211.3

论文编号：10006SY0505218

硕 士 学 位 论 文

高超声速飞行器稀薄气体效应

的数值模拟

作者姓名

王勃

学科专业

流体力学

指导教师

吴颂平　教授

培养院系

航空科学与工程学院

---

Page 2

Numerical Simulation of Rarefied Gas Effects

on Hypersonic Vehicle

A Dissertation Submitted for the Degree of Master

Candidate：Wang Bo

Supervisor：Prof. Wu Songping

School of Aeronautic Science and Engineering

Beihang University, Beijing, China

---

Page 3

中图分类号：V211.3

论文编号：10006SY0505218

硕 士 学 位 论 文

高超声速飞行器稀薄气体效应的

数值模拟

作者姓名

王 勃

申请学位级别 工学硕士

指导教师姓名 吴颂平

职 称

教授

学科专业

流体力学

研究方向

计算流体力学

学习时间自 2005 年 9 月 15 日 起至 2007 年 12 月 14 日止

论文提交日期 2007 年 12 月 7 日 论文答辩日期 2007 年 12 月 14 日

学位授予单位

学位授予日期 2008 年 月 日

---

Page 4

关于学位论文的独创性声明

本人郑重声明：所呈交的论文是本人在指导教师指导下独立进行研究工作所取得的

成果，论文中有关资料和数据是实事求是的。尽我所知，除文中已经加以标注和致谢外，

本论文不包含其他人已经发表或撰写的研究成果，也不包含本人或他人为获得北京航空

航天大学或其它教育机构的学位或学历证书而使用过的材料。与我一同工作的同志对研

究所做的任何贡献均已在论文中作出了明确的说明。

若有不实之处，本人愿意承担相关法律责任。

学位论文作者签名：

日期： 年 月 日

学位论文使用授权书

本人完全同意北京航空航天大学有权使用本学位论文（包括但不限于其印刷版和电

子版），使用方式包括但不限于：保留学位论文，按规定向国家有关部门（机构）送交学

位论文，以学术交流为目的赠送和交换学位论文，允许学位论文被查阅、借阅和复印，

将学位论文的全部或部分内容编入有关数据库进行检索，采用影印、缩印或其他复制手

段保存学位论文。

保密学位论文在解密后的使用授权同上。

学位论文作者签名：

日期： 年 月 日

指导教师签名：

日期： 年 月 日

---

Page 5

- I -

摘 要

飞行器在再入过程中要经历 70～130Km 的过程，此阶段的特点是高速、高温和低

密度，所以有两个方面的效应是必须要考虑的，即高温效应和稀薄效应。在这个阶段，

基于连续介质的控制方程近似成立，但物面处的边界条件必须进行修正。本论文的目的

在于研究稀薄效应和高温效应对流场以及其它流场参数的影响。

本文首先假设气体为完全气体，在引入速度滑移和温度跳跃边界条件的基础上，研

究了不同网格尺度、不同边界条件（滑移、无滑移）、不同 Kn 数和不同壁面协调系数对

流场结构，壁面处温度、速度以及热流等参数的影响，得到了一些有意义的结论。

接着本文考虑了气体的高温效应，建立了 5 组分 17 化学反应的化学非平衡模型，

针对球头和双椭球外形进行了数值模拟，研究了不同 Ma 数下流场的结构及流场参数的

变化，并讨论了气体模型对气动力、气动热的影响。

关键词：高超声速，稀薄，速度滑移，温度跳跃，Kn 数，化学非平衡，热流

---

Page 6

- II -

Abstract

In the reentry process, the hypersonic vehicle would go through a flight region about 70

to 130 kilometers high in altitude, which is characterized by its high velocity, high

temperature and low density. Two effects, the high temperature effect and the rarefied effect,

should be taken into account. In the flight region illustrated above, the governing equation

based on the Continuous Medium Hypothesis would be proper approximately with a

correction of the boundary condition on the wall. The objective of this paper is to investigate

the impact of the rarefied effect and the high temperature effect, mainly on flow field

characteristics and flow field parameters.

Firstly, based on the perfect gas hypothesis, coupled with a velocity slip and temperature

jump as the wall boundary condition, the effects that different grid scale, different boundary

condition (slip or nonslip boundary), different Knudsen Number and different wall harmony

coefficient made on the flow structure are investigated in this paper, as well as the impact that

wall temperature, velocity and heat transfer made on flow field. Some useful conclusions

were made.

Secondly, the high temperature effect of gas is considered and a five component

chemical nonequilibrium model is adopted to take into account the chemical reaction. By the

numerical tests of a semi-sphere and a double-ellipsoid, the flow structure and flow field

parameters change are investigated at different Mach Number , the aerodynamic and the

aerothermdynamic are discussed meanwhile.

Keywords:

Hypersonic, Rarefied, Velocity slip, Temperature jump, Knudsen number,

Chemical nonequilibrium, Heat transfer

---

Page 7

- III -

目 录

第一章 绪论····························································································································1

1.1 研究背景及意义..................................................................................................1

1.2 真实气体效应的特点和难点..............................................................................2

1.3 高超音速流动数值模拟研究..............................................................................5

1.3.1 化学非平衡流的数值模拟 ························································································5

1.3.2 稀薄气体效应的数值模拟 ························································································7

1.4 论文构成..............................................................................................................7

第二章 滑移边界条件概述·····································································································9

2.1 流动特性的分区..................................................................................................9

2.2 边界条件的推导................................................................................................ 10

2.2.1 速度滑移边界条件·································································································· 11

2.2.2 温度跳跃边界条件···································································································15

2.3 关于壁面协调系数............................................................................................ 18

2.4 向多组分非平衡流情况的推广........................................................................ 19

2.5 本章小结............................................................................................................ 20

第三章 N-S 方程计算方法··································································································· 21

3.1 笛卡尔坐标系下的 N-S 方程............................................................................ 21

3.2 一般坐标系下的 N-S 方程................................................................................ 22

3.3 N-S 方程的无量纲化.......................................................................................... 24

3.4 时间离散............................................................................................................ 25

3.4.1 时间推进···················································································································25

3.4.2 时间步长的确定·······································································································27

3.5 空间离散............................................................................................................. 27

3.5.1 空间格式···················································································································28

3.5.2 熵修正函数···············································································································28

3.5.3 限制器函数···············································································································29

3.6 粘性项的离散..................................................................................................... 29

3.7 定解条件............................................................................................................ 30

---

Page 8

- IV -

3.7.1 初始条件··················································································································30

3.7.2 边界条件··················································································································30

3.8 本章小结............................................................................................................ 34

第四章 完全气体下计算结果与分析··················································································· 35

4.1 程序有效性验证................................................................................................ 35

4.1.1 二维圆柱有效性验证 ······························································································35

4.1.2 三维球头有效性验证 ······························································································36

4.1.3 大 Kn 数的验证·······································································································37

4.2 不同边界条件对流场的影响............................................................................. 39

4.3 Kn 数对流场的影响 ........................................................................................... 41

4.3.1 对圆柱绕流的模拟··································································································42

4.3.2 对球头绕流的模拟··································································································44

4.4 壁面协调系数的影响........................................................................................ 45

4.4 本章小结............................................................................................................ 47

第五章 化学非平衡流的控制方程····················································································· 48

5.1 控制方程............................................................................................................ 48

5.1.1 直角坐标系下的控制方程 ·······················································································48

5.1.2 热力学性质和状态方程 ··························································································49

5.1.3 化学反应源项··········································································································50

5.1.4 输运系数··················································································································50

5.2 无量纲化............................................................................................................ 51

5.3 坐标变换............................................................................................................ 52

5.4 数值方法............................................................................................................ 52

5.5 边界条件............................................................................................................ 53

5.6 本章小结............................................................................................................ 53

第六章 数值计算结果及讨论······························································································· 54

6.1 绕球头的化学非平衡流动数值模拟................................................................ 54

6.1.1 计算条件···················································································································54

6.1.2 计算结果及分析·······································································································54

6.2 绕双椭球外形的化学非平衡流动数值模拟.................................................... 57

---

Page 9

- V -

6.2.1 计算条件···················································································································57

6.2.2 计算结果分析··········································································································58

6.3 本章小结............................................................................................................ 61

结 论···································································································································· 62

参 考 文 献 ···················································································································· 63

附录 A 各组元的热力学特性数值······················································································· 65

附录 B 牛顿迭代公式求温度······························································································· 68

附录 C 大气组成··················································································································· 70

攻读硕士期间发表的论文····································································································· 71

致 谢·································································································································· 72

---

Page 10

---

Page 11

北京航空航天大学硕士论文

- 1 -

第一章 绪论

1.1 研究背景及意义

从载人航天飞行器成功以来至今己经有 40 年的历史，人类征服太空已经一步一步

地走向现实。40 年来计算机及计算技术突飞猛进，人们在设计新的飞行器型号再也不需

完全依靠地面实验和飞行实验，而是首先求助于数值模拟。尤其是现在尚没有能完全模

拟这种飞行条件的地面设备可提供地面实验，而飞行实验则耗资巨大而不能随意进行。

在发达国家，利用数值模拟进行航天飞行器的设计在 20 世纪 80 年代就已经很普遍了。

再入飞行器在再入过程中将经历 130～70km 高空区域，由于低密度粘性效应，对返

回舱的气动特性将发生相当大的影响，例如，稀薄气体的阻力系数为连续流区的 1-3 倍。

飞船返回舱在高空受到的阻力，尽管其量值不大，但由于飞行时间较长，阻力将改变其

飞行轨迹，影响飞行器落点精度。另外，伴随着高温，波后气体将发生离解甚至电离，

这对飞行器的气动力、气动热也有显著的影响，因此高温效应也是不得不考虑的。

随着飞行高度和速度的提高，空气动力的数值模拟也从宏观走向微观。宏观模型是

把气体当作连续介质，关于它的描述是以流动参数，如速度、密度、压力和温度在空间

和时间的变化为依据。N-S 方程就是连续介质的数学模型。微观模型或分子模型是把气

体当作无数单独的分子，提出其每一时刻的位置、速度和每一分子的状态，在这一层次

上的数学模型为 Boltzmann 方程。

理论上来说，所有气体流动问题，都可以通过微观气体模型，即求解 Boltzman 方

程得以解决，其中也包括连续介质问题，但分子模型对计算机的资源要求是很高的，因

此能够避免使用就应当尽量避免。经过理论计算和实际设计经验，大约以 70km 为界，

在这个高度以下的问题可以使用连续介质模型。

在连续介质模型中，当气流速度提高到高超声速时就不能将气体当作理想气体看

待，而必须考虑真实气体效应的影响，因此在用宏观模型的情况下又多出一个层次，这

就是带真实气体效应的连续介质模型。这一模型既把气体当作宏观模型又考虑由于流速

提高到高超声速时温度升高引起的微观分子结构上的变化。也就是说，在高超声速阶段

必须考虑真实气体的影响，即高温情况下分子的离解和电离。

在再入航天器的高超声速流动中，气体的稀薄效应也是不能忽略的。此时，传统的

物体表面上边界条件取无滑移条件将失效，而必须考虑更加真实的速度滑移和温度跳跃

边界条件。

---

Page 12

第一章 绪论

- 2 -

高超声速飞行器的气动、热力环境的研究采用理论分析、地面设备模拟、飞行实验

和数值计算作为主要手段来获取数据。但是，地面模拟和飞行实验将耗费大量的人力、

物力和财力。随着计算机技术和数值方法的迅速发展，为计算流体力学的发展创造了十

分有利的条件，因此，计算流体力学在高超声速飞行器设计中就占据了更加重要的地位。

1.2 真实气体效应的特点和难点

所谓的高超音速流动就是当飞行器的飞行马赫数增加到一定程度以后，某些物理现

象变得越来越重要的流动。这些现象包括高马赫数产生的高度非线性的流体动力学特性

和由于流动能量很大而引起的高温物理化学特性，如很薄的粘性激波层，很大的熵梯度，

粘性/无粘相互干扰和真实气体效应，其中真实气体效应对高超音速流动的影响最为显

著。

真实气体效应对高超音速飞行器的气动力有重要的影响，是高超音速流动的重要特

征之一，对于复杂外形的再入飞行器，影响更大。美国登月飞船 Apollo

2 号、4 号、6 号的飞行实验结果表明，在高马赫数范围其指挥舱的配平攻角比风洞试验

的预测值要大 2°～4°[1]。由于风洞试验无法完全模拟高温气体流动，所以对造成这种

差别的原因是否来源于真实气体效应，长期以来一直存在不同的看法[2]。八十年代以来，

计算流体力学的发展为澄清这个问题创造了条件。根据文献[3]的报导，美国用无化学反

应气体和有化学反应气体两种模型计算了 Apollo 飞船指挥舱的配平攻角，其差别和风洞

预测结果与飞行试验结果的差别一致，从而说明了研究真实气体效应，特别是研究有化

学反应的气体流动的重要性。在美国航天飞机 STS-1 的飞行中，为了保持其配平攻角为

40°，机身襟翼偏转了 16°，而风洞预测只需偏转 7°[4]。事后的分析表明，这一现象

就是由于真实气体效应引起迎风面和背风面压力分布变化而产生抬头力矩的结果[5] [6]

[7] 。文献[8]在风洞试验数据上加入了马赫数效应、真实气体效应和粘性效应的修正后，

就与飞行试验结果一致。这三种因素中，影响最大的是真实气体效应。文献[9]考虑了七

种组元求解薄层纳维－斯托克斯方程，其计算结果已比较接近飞行试验结果。

随着我国 863 计划航天技术领域预研工作的发展，近年来国内也越来越重视对真实

气体效应的研究，并已从早期的方法研究逐渐转向开展对各种外形航天飞行器含化学反

应绕流的数值模拟研究。对于高超声速钝体绕流，不仅物面附近的边界层是高温流动区

域，而且在整个激波层流动区域内气体温度都是很高的。当温度升到大约 800K，空气

分子的振动能被激发，空气的比热比不再是常数 1.4；当温度达到 2500K 时，空气中的

---

Page 13

北京航空航天大学硕士论文

- 3 -

氧分子开始分解；到 4000K 时，氧分子基本上完全离解，这时氮分子也开始分解；当温

度达到 9000K 时，氮分子也基本离解完毕，氧原子和氮原子开始电离，形成部分电离的

等离子体。同时，在 4000K 到 9000K 之间，由于离解气体的复合反应，可形成少量的

一氧化氮（NO）分子，其中一部分经电离形成 NO 离子和自由电子 e。概括来说，高温

真实气体效应就是指空气在高温下产生的振动能激发、分子离解、原子复合、组元之间

的化学反应，以及电离等现象对流动特性的影响。

图 1.1 空气分子振动激发、离解和电离的温度区分[10]

一般的高超声速飞行器均具有复杂的钝头外形结构。在高超声速飞行下，钝头前将

产生非常强的弓形激波，绕飞行器的空气受强激波以及飞行器摩擦等因素的加热而发生

离解和电离，此时将产生真实气体效应。真实气体效应包括三种流动类型:

（1）冻结流。其流动特征时间小于化学反应时间和振动松弛时间，即化学反应速

率为零、振动松弛时间趋于无穷的情况；

（2）平衡流。其流动特征时间远大于反应特征时间和振动弛豫时间。这表明流体

微元在流场中滞留的时间足够长，微元的热力学特性和化学特性有充分的时间自身调

节，达到局部热力学平衡和化学平衡状态；

（3）非平衡流。其流动特征时间和化学反应时间与振动松弛时间可比，需要考虑

有限的化学反应速率和振动松弛时间。

高温真实气体效应对高超声速飞行器产生重要的影响。真实气体效应特别是表而催

化效应对高超声速飞行器的气动加热将改变空气的热力学性质，不仅影响飞行器的气动

力系数和物而气动加热，而且有可能导致飞行器飞行中的“通讯中断”。

非平衡效应是真实气体效应中的主要问题。由于非平衡流动非常复杂，高温所引发

的分子输送、化学反应、振动激发、辐射等多种非平衡过程的物理化学机理尚未完全掌

---

Page 14

第一章 绪论

- 4 -

握；高超声速高温气体动力学问题，要求求解一组复杂的控制方程组。对于化学非平衡

问题，控制方程为带化学反应源项的 N-S 方程，总的连续方程由组元连续方程组代替；

对于热化学非平衡问题还包括振动能量方程，计算量相当巨大。

另外，高超声速再入飞行器再入过程中，从高空到低空要经历自由分子区、稀薄过

渡区、滑移区，然后经历连续区。判断不同流态的相似参数通常采用努森数（Knuslen），

其定义为Kn

L

l

=

，这里的l 是分子的平均自由程长度，与大气层的高度有关；L 为物

体的特征长度。不同区域对应着不同得处理方法，如下图。

图 1.2 以当地 Kn 数确定使用的模型[11]

由努森数的定义可以看出，绕飞行器表面周围的流动状态不仅与飞行高度有关，而

且还与飞行器的特征尺寸有关。在有些情况下，连续流的假设在较高的高度上仍然可以

成立，但由于稀薄效应影响造成的边界层背离经典薄边界层概念的现象，可能在较低的

高度就会变得很明显，稀薄气体效应对气动热交换率有较大的影响。研究结果表明[12]，

稀薄气体效应对气动热的影响趋势与稀薄度有关，随着稀薄度增大，稀薄效应增加，根

据边界层理论计算的气动加热率低于实际试验结果，而随着稀薄度的进一步增加，稀薄

气体效应的进一步增大，边界层理论计算值又出现大于试验结果的现象、即随稀薄度的

增大，根据边界层理论计算的高超声速飞行器表面的气动加热率，对实际情况先是估计

不足，继而是估计过量。

关于再入飞行器稀薄过渡区气动热问题，尽管从二十世纪六十年代起，就是许多研

究和设计者非常感兴趣的课题，且已进行了大量的理论分析和实验研究。但迄今为止，

对稀薄过渡流高超声速气动加热的计算，还仅限于驻点区和沿飞行器迎风中心线的最大

热流分布。对有攻角三维流动情况下沿整个再入飞行器表面(沿纵向和周向)的热流分布

还没有合适的计算方法。新一代高超声速再入弹头的发展趋势是小型化，因而对气动热

---

Page 15

北京航空航天大学硕士论文

- 5 -

的预测提出了新的问题，一是对热防护和气动热的计算提出了更高的要求，要求气动热

的预测要更为准确，二是由于小型化弹头的特征尺寸比较小，绕其周围的流动到较低的

高度(约 40 公里)仍然为稀薄过渡流，可能从 40 公里到 80 公里都属于稀薄过渡流，也就

是说，小型化弹头的再入过程中，在较大的高度跨度范围内都是处于稀薄过渡流区，因

而，准确地计算稀薄过渡流的稀薄气动热，对小型化弹头的防热设计是非常重要的。对

航天飞机、载人飞船以及新一代跨大气层飞行器来说，其峰值加热区都是在 60 至 80 公

里的高空，都是处于稀薄过渡流区，因而，对稀薄过渡区的气动热同样需要进行认真地

试验和理论研究。虽然，由于电子计算机和数值模拟能力的极大提高，可以通过结合滑

移边界条件直接数值模拟求解 N-S 方程、直接模拟蒙特卡罗方法(DSMC)或者结合考虑

激波和壁面滑移、温度跳跃等稀薄效应求解粘性激波层(VSL)方程等方法来计算稀薄过

渡流三维情况下的气动热问题，但所需要的计算机计算时间和费用，对飞行器设计的实

际工程需要来说是不适宜的。

1.3 高超音速流动数值模拟研究

高超音速气流穿过激波后，高温所引发的分子输运、化学反应及辐射传热等多种非

平衡过程的物理化学机理还没有完全研究清楚，这是目前数值模拟的困难之一。由于新

一代先进天地往返运输系统的外形都非常复杂，其周围流场中存在强激波、激波波系干

扰、激波/粘性边界层干扰等复杂现象，这些都要求数值计算格式具有高分辨率和高精度。

而对于化学非平衡流动，除去以上的困难之外，组元连续方程中的化学反应源项进一步

增加了计算的难度。当组元的化学反应时间尺度和流体流动的时间尺度相差几个数量级

时，就会导致刚性问题。这类问题会使得数值模拟计算的效率下降，甚至导致某些算法

无法收敛到定常解。另外，高超音速飞行器表面热流的准确计算也是亟待解决的问题之

一。

1.3.1 化学非平衡流的数值模拟

为满足先进天地往返运输系统发展的需求，人们逐步从模拟高超音速完全气体发展

为模拟高超音速化学反应流动等非平衡流动。计算机技术的迅速发展则为此提供了物质

基础。对于高超音速化学反应流动的模拟也分为平衡流和非平衡流两类：平衡流的计算

采用自由能极小化的方法，在不涉及具体反应模型的情况下，对达到平衡态的流动进行

模拟。非平衡流的计算则采用有限速率化学反应模型，直接模拟化学反应的非平衡过程。

与高超音速完全气体流动的数值模拟相比，除求解大型非线性方程组所遇到的困难

---

Page 16

第一章 绪论

- 6 -

外,化学非平衡反应流动组元连续方程中的化学反应源项进一步增加了计算的难度，如何

解决方程中出现的众多不同量级的时间尺度等问题，提高计算效率和精度，成为数值模

拟中的主要困难。

因此，高超音速化学非平衡流数值模拟多采用稳定性更好的隐式格式，以放宽对时

间步长的限制。Bussing 和 Murman 提出了点隐式格式[13]，将化学反应源项做了隐式处

理。该格式的实质是通过预调节(Pre-conditioning)矩阵将时间步长统一到伪时间步长上，

使每一个反应均在各自特征时间上进行；Yee 和 Shinn 将隐式 TVD 格式运用于化学非平

衡流[14]。当采用隐式格式方法求解时，随着化学组元个数的增加，也增加了矩阵的阶数

和计算的复杂性。对于三维问题的隐式方法，由于需要做大量的分块矩阵求逆，计算量

更加庞大，因此 LU 算法[15]不需要分块矩阵求逆的特点受到了人们的重视，也被用于化

学非平衡流的数值模拟。但是由于化学反应源项的存在，要想完全避免分块矩阵求逆，

还需要采用对角化处理，对角化处理同时还是降低刚性的有效方法。对于按坐标方向的

近似因子分解，Pulliam 和 Chaussee 首先提出了对角化形式[16]，Eberhardt 和 Imlay 给出

了化学反应源项 Jacobi 矩阵的对角化形式[17]。由于刚性问题随着马赫数的增加而变得更

加严重，不同格式的性能在不同马赫数下会有不同的表现。国内对于高超音速化学非平

衡流场 N-S 方程数值模拟计算开始于上世纪 90 年代初期，中国空气动力发展研究中心、

中国科学院、北京空气动力研究所、北京航空航天大学等单位都相继开展了卓有成效的

工作。黎作武，张涵信利用杂交通量分裂法建立了隐式 NND 格式，并求解了完全气体

和化学非平衡气体 N-S 方程[18]。瞿章华等人开展了化学非平衡流隐式向量算法研究[19]。

董维中在他的博士论文[20]中，对轴对称 N-S 方程和 11 组元的离解电离空气绕流高超声

速钝锥体流场进行了单温模型，双温模型和三温模型的数值计算，分析了温度模型和化

学模型对流场特性的影响。传统 Harten-Yee 型的 TVD 格式用于计算钝头体高超音速绕

流时，由于所采用的熵修正函数无法提供足够的数值粘性捕捉驻点前的激波，导致前驻

点流场中出现压力集中的问题。李椿萱，褚以华[21]提出了一类 non-MUSCL 的 Harten-Yee

型的 TVD 格式，解决了这一问题。该文献采用此格式计算了绕圆柱的二维粘性高超音

速化学非平衡流动。李椿萱，周华[22]以 Harten-Yee 型 TVD 格式为基础，构造了计算化

学非平衡流动的有限体积法，该算法为了解决刚性问题，对化学反应源项采用点隐式格

式，并对源项的 Jacobi 矩阵进行缩并，提高了计算效率。谢锦睿 [23]通过对 TVD 格式中

几类常用的熵修正函数的研究，提出了具有更高稳定性和适用于热流计算的改进型熵修

正函数，并研究了化学非平衡流下法向网格的划分原则。欧阳水吾，谢中强[24]详细介绍

---

Page 17

北京航空航天大学硕士论文

- 7 -

了模拟高超音速化学非平衡流动的数值方法。

1.3.2 稀薄气体效应的数值模拟

对于稀薄气体动力学的研究是从 19 世纪 Maxwell 等人的工作开始的，但当时的研

究仅限于速度很低的情况。第二次世界大战后与空间探索引起的高空高速飞行的蓬勃发

展相适应，稀薄气体动力学走上了前台，取得了显著的进展。钱学森提出了稀薄气体动

力学流动领域的划分[25]，为现代稀薄气体动力学的研究与发展作出了开创性的工作。

迄今为止，对于流动的数值模拟主要是通过求解基于连续性假设的 NS 方程和基于

微观分子模型的直接模拟蒙特卡罗(DSMC)方法得到的。这两种方法在研究高超声速稀

薄气体流动中互为补充。

DSMC 已经成为研究稀薄效应的最有效的手段之一。理论上来说，所有流动都可以

由 DSMC 方法求解，但是该方法受网格划分、时间步长的限制，需要耗费大量的计算

机资源，特别是当流动接近连续区时，DSMC 需要耗费更大的计算机资源，因此，从本

质上来说，DSMC 方法在大 Kn 数下才是可行的。

另一方面，理论和实践都已经证明，在连续区求解 NS 方程是非常有效的。但是，

随着 Kn 数的增大，N-S 方程就不再准确了。Maxwell[26]证实了在分子平均自由程与流动

的特征长度相当的滑流区无滑移边界条件将不再成立，而应当对边界条件作一定的修正

[27] [28]。

对于滑移区的流动，很多人作了研究。X.Zhong[29] [30] ，Keon-Young Yun [31]等人通

过求解 Burnett 方程研究了滑移区的流动，F.W.Vogenitz[32]等人运用 DSMC 方法模拟了

滑移区的流场，A.C.Jain，H.G.Woods[33]等人通过求解粘性激波层方程得到了类似的结果。

1.4 论文构成

本文的主要目的在于考察稀薄气体效应及高温效应对流场和流动参数的影响。采用

有限差分方法，空间离散采用 Harten-Yee 的二阶迎风 TVD 格式，时间迭代采用隐式

LUSGS 方法，数值求解 N-S 方程。本为分别采用完全气体模型和化学非平衡模型，结

合滑移边界条件，分别对多个外形进行了数值模拟。论文所完成的工作包括：

1．查阅并推导了滑移、跳跃边界条件，并将其转化为本文所需要的形式。

2．研究了实验室原有的完全气体及化学非平衡流计算程序，并针对本论文研究内

容对其进行了补充和完善。

3．由代数方法生成网格，采用改进的程序，分别对圆柱、圆球、双椭球进行了一

---

Page 18

第一章 绪论

- 8 -

系列的数值模拟，着重研究了不同网格尺度、不同边界条件（滑移、无滑移）、不同 Kn

数、不同壁面协调系数和不同 Ma 数对流场结构，壁面处温度、速度以及热流等参数的

影响，得到了一些有意义的结论。

4．综合考虑稀薄气体效应和高温效应，并对球头和双椭球进行了初步的数值模拟。

论文的第一章是对本文的研究背景和意义、真实气体效应的特点和难点以及真实气

体的研究现状的叙述。第二章主要对滑移边界条件进行了简要的推导和总结，并归纳出

本文需要的边界条件。第三章对完全气体的控制方程进行了转换，并对空间离散和时间

推进方法进行了介绍，最后给出了所需的边界条件，包括物面边界、远场边界、奇性轴

条件和对称边界。第四章给出了完全气体下的计算结果，并对结果进行了分析。第五章

给出化学非平衡流场的数学模型，包括流场控制方程以及相关的热力学性质、输运系数。

第六章给出了化学非平衡下的一些计算结果，并对结果进行了分析。

---

Page 19

北京航空航天大学硕士论文

- 9 -

第二章 滑移边界条件概述

2.1 流动特性的分区

在气体绕流问题中，假定流动问题的特征长度为 L，分子运动的特征长度取为热运

动的平均自由程λ 。现假定

L

λ << ，于是当一个气体分子 A 入射到物体表面后，从平均

意义上说反射分子将在λ 距离内与另一个分子 B 相碰撞。因为

L

λ << ，且已假定碰撞分

子的散射按立体角均匀分布，因此碰撞分子各有一半的机会再次入射到物体表面上。假

定

L

λ << 意味着在物体表面附近存在大数目的气体分子，入射分子与反射分子之间碰撞

数目更是巨大。气体不断的把入射分子带到壁面附近，气体分子之间的碰撞导致气体分

子在物体表面附近的聚集。显然，在物体表面附近聚集层的厚度l 与分子的热运动平均

自由程λ 具有同一量级。通常将这一层称为 Knudsen(努森)层。

当λ 相对 L 趋近零时，努森层的厚度l 也趋于零，意味着入射分子被物体表面完全

吸附。于是气体分子微观特征量的集体表现，与物体表面分子微观特征量的集体表现一

致，也就是在物体表面上边界条件是无滑移边界条件和等温壁条件。当λ 相对与 L 增大，

形成努森层，若仍保持

L

λ << ，那么在宏观意义上考察这一流动时，可以不必顾及努森

层的存在，但是物体表面上的附着壁和等温壁边界条件已不再成立，必须考虑入射到物

体表面上的气体分子与壁面分子之间动量和能量的交换，也就是速度滑移和温度跃跳效

应。以上两种情形，由于

L

λ << ，处于物体边界以外的区域中气体密度是相对较大的，

经典气体动力学的连续介质模型假设仍然成立。因此控制气体流动的方程是经典的流动

力学方程，即纳维－斯托克斯方程。随着λ 相对 L 继续增大，努森层的厚度与 L 可相比

拟的时候，努森层与整个流动区域已浑然一体，不可能再单独考察努森层的流动状态。

这已流动状态称为过渡流动。当λ 相对 L 趋于无穷大时，气体分子之间不存在相互碰撞，

碰撞仅仅出现在入射的气体分子与物体表面之间，这一流动状态通常称为自由分子流。

Knudsen（努森）数表达式为

/

Kn

L

λ

=

。其中λ 是气体粒子的平均自由程长度，而

L是问题中物体的特征尺寸，如飞行体的直径或长度，边界层厚度，激波厚度等，按照

Knudsen 数的数值范围，气体的传热与流动问题可以分类如下[34]：

1.

3

10

Kn

−

<

称为连续介质（continuum）流区。在这一区域中，气体分子之间的碰撞

频率远大于气体分子与物体之间的碰撞频率。Navier-Stokes 方程、Fourier 热传导定律质

适用，物面处的无速度滑移和无温度跳跃假设成立。这是传统气体动力学的研究领域。

2． 3

1

10

10

Kn

−

−

<

<

，为速度滑移和温度跳跃（velocity slip and temperature jump）区。

---

Page 20

第二章 滑移边界条件概述

- 10 -

此时气体粒子之间的碰撞频率仍然比气体粒子与物体表面之间的碰撞频率高很多，但稀

薄气体效应（rarefaction effect）已经不能忽略。此区域中，在流场内 Navier-Stokes 方程、

Fourier 热传导关系式仍然成立，由于气体稀薄程度增加，边界条件需考虑气体与固体交

界处的速度滑移和温度跳跃。

3.

1

10

10

Kn

− <

< 为过渡（transition）区。在这一区域，分子的平均自由程只与物体

的特征长度 L 具有同一量级。气体分子之间的碰撞和气体分子与物体表面的碰撞对流动

具有同等重要的影响。连续介质假设不再成立，这一区域内的传热和流动问题的严格处

理是非常困难的，往往要从包括碰撞积分项 Boltzmann（玻尔兹曼）方程出发，采用数

值算法（如 Monte Carlo 算法）计算。工程计算中则常用一些半经验性的关系式。

4.

10

Kn >

为自由分子(free-molecule)区。这一区域内，分子平均自由程远远大于物

体的特征长度值，可以忽略气体分子之间的碰撞而仅仅考虑和气体分子与物体表面的相

互作用。这一流动特征导致问题处理的极大简化。该区域内传热与流动问题的求解需依

靠动力论方法。

2.2 边界条件的推导

根据Schaaf和Chambre[35]以及Beskok和Karniadakis[36]等人的分类，当 3

1

10

10

Kn

−

−

<

<

时，气体流动与换热将进入速度滑移和温度跳跃区，该区的特点是:气体分子的平均自由

程与流道特征尺寸相比仍然不大，但己经不能忽略不计，因而在远离边界的中心区域，

动量和能量传输过程中气体分子之间的碰撞仍起支配作用；但是此时的 Knudsen 数又不

是小得可以忽略不计，因而气体的稀薄效应应加以考虑，因此在紧靠壁面的、厚度约为

气体分子平均自由程长度的气体层（Knudsen 层）内的动量和能量传输过程必须用分子

动力论来处理，从而导出边界上将出现的速度滑移和温度跳跃条件。对滑移流区的流动

与传热问题，采用的处理方法是：在主流中心区内仍采用 Navier-Stokes 方程和 Fourier

定律，但边界条件则采用动力论方法导出的壁面处的滑移速度和跳跃温度修正，这样的

求解结果可以满足工程计算的精度要求。鉴于本论文主要采用滑移模型研究流动与换热

问题，因此首先就一阶和高阶滑移边界条件及温度跳跃边界条件对前人的工作做简要介

绍。更详细的分析参见文献[34] [37]。

---

Page 21

北京航空航天大学硕士论文

- 11 -

2.2.1 速度滑移边界条件

2.2.1.1 速度滑移边界条件的导出

图 2.1 壁面附近气体滑移速度分析简图

滑移区壁面处的滑移速度表达式是首先由Maxwell通过对Knudsen层内动量交换的

分析导出的。设气体仅在 x 方向上具有宏观速度( x

U ) ,且壁而附近既存在速度梯度又存

在温度梯度。在此种情况下，气体分子的速度分布函数不再是 Maxwell 分布，气体也不

再是平衡气体。但当气体只是稍稍偏离平衡气体时，气体分子的速度分布函数可以在

Maxwell 分布函数的基础上加上小的摄动项予以修正。因此此时在 Knudsen 层外缘处向

壁面运动的气体分子的速度分布函数可表示为：

(

)

2 2 2

exp

0

3/2

2

/

2

/

v v v

n

x y z

s

f

T ms

T ms

κ

πκ

⎡

⎤

+ +

⎢

⎥

−=

−

⎢

⎥

⎢

⎥

⎣

⎦

(

)

(

)

( )

x

2

y

z

x

s

y

s

T

T

T

u

v u

v

v

v u v

x sx

y

z

y

T /m

m

r k

轾

ïï

抖

抖

ï

犏

ý

?

+

+

-

-

犏

ï

抖

抖

犏

ïï

臌

（2.1）

因此，从 Knudsen 层外缘入射到物体表面的分子通量可计算如下：

0

s

y

x

y

z

s

T

v dv dv dv

n

2 m

k

y

p

¥

¥

- - ?

= -

=

ò 蝌

（2.2）

而从 Knudsen 层外缘向壁面输运的切向（x 向）动量可由下式计算

( )

0

y

x

x

y

z

v mv f dv dv dv

¥

¥

-

-

---

Page 22

第二章 滑移边界条件概述

- 12 -

( )

s

x

s

s

s

s

s

1

1

T

T

U

n

mu

k

2

5 2 T /m

2 m

x

y

k

m

pk

p

骣

骣

÷

ç

抖 ÷

ç

÷

ç

÷

÷

ç

=

-

+

ç

÷

÷

ç

÷

ç

ç

÷

桫

ç

抖

桫

（2.3）

图 2.2 镜反射、漫反射示意图

若从 Knudsen 层外缘入射到表面的分子中，漫反射分子份额为 v

s ，其余部分为镜反

射。由于镜反射分子切向动量与入射分子相同，因此这部分分子对 Knudsen 层中分子的

动量交换没有贡献，因此入射分子输运的切向动量为

( )

s

x

s

s

v

s

s

s

T

1

T

1

U

n

mu

k

2 m

x

2

y

5 2 T /m

k

s

m

p

pk

轾

骣

骣

犏

抖

÷

ç

÷

ç

÷

犏

-

+÷ ç

ç

÷

÷

ç

÷

ç

÷

ç

犏

桫抖

桫

臌

（2.4）

另一方面，在表面漫反射的分子对动量交换没有贡献，即由壁面向外漫反射出去的

分子在切向方向的动量为零，因此式（2.4）即为 Knudsen 层中气体分子交换的净切向

动量。又因为 Knudsen 层内净切向动量交换应与表面的粘性应力相等，故

( )

s

x

x

s

s

v

s

s

s

s

T

1

T

1

U

U

n

mu

k

2 m

x

2

y

y

5 2 T /m

k

s

m

m

p

pk

轾

骣

骣

÷

骣

犏

抖

÷

÷

ç

ç

÷

ç

÷

÷

犏

-

+

=

÷

ç

ç

ç

÷

÷

÷

ç

ç

÷

ç

÷

÷

珑

犏

桫抖

桫

桫

臌

（2.5）

其中：

s

s

s

1

2 T

v

mn

2

m

k

m

r l

l

p

=

=

，

s

s

v

5

3 j

2 T

k

vc

4

4

m

m

k

k

r

l

rl

p

骣

֍֍

=

֍ ֍

桫

＋

＝

则从（2.5）式可最终导出滑移速度的表达式为：

v

s

w

v

s

s

s

2

3 j

u

T

u

u

4

n

T

s

m

l

s

r

t

骣

骣

÷

-

+ ç

抖 ÷

÷

ç

ç

÷

÷

ç

ç

-

=

+÷

÷

ç 麋

ç

÷

桫

÷

ç

抖

桫

u

s

s

s

3 j

u

T

4

n

T

m

x

r

t

骣

骣

÷

+ ç

抖 ÷

÷

ç

ç

÷

÷

ç

ç

=

+÷

÷

ç 麋

ç

÷

桫

÷

ç

抖

桫

（2.6）

式中 s

u 为气体在近壁面处的滑移速度， w

u 表示壁面切向运动速度， U/ y

抖

是切向

速度

x

U 在法向的梯度， T/ x

抖

是温度在切向的梯度，n 是壁面的法线方向， u

x 为速度

滑移系数(coefficient of slip)，具有长度量纲，表示为

v

v

2

u

s

l

s

x

-

=

（2.7）

---

Page 23

北京航空航天大学硕士论文

- 13 -

式中 l 为气体分子的平均自由程，由下式计算

16μ

=

5ρ 2πRT

l

（2.8）

v

s

称为切向动量协调系数（tangential momentum accommodation coefficient），反映

的是气体分子在壁面处的漫反射分数，事实上根据其定义：

i

r

v

i

t t

s

t

-

=

（2.9）

也反映了入射气体分子与壁面之间动量交换的完善程度，式中 i

t

和 r

t

分别表示入

射分子与反射分子所具有的平均动量。在工程应用中经常经验性的认为分子在固体表面

的切向动量协调系数（TMAC）为 1，但是大多数学者测量得到的切向动量协调系数是

小于 1 的，其对固体表面的环境条件比较敏感，如气体的种类、表面材料的种类、表面

的粗糙度、表面的氧化和吸附物等有关，在工程应用中常将动量协调系数取在 0.75～0.85

之间，如图 2.3。

很明显，当

v

1

s = 即气体分子在表面完全漫反射时，滑移系数 u

x l

= ；而当 v

0

s = 时，

即气体分子在表面完全镜反射时，式（2.6）无意义，此时由于气体分子与表面之间没有

净切向动量交换，表面无切应力，因而表面附近气体无速度梯度，从而表面处的气体滑

移速度等于气体中的宏观速度的平行与表面方向的分量。

图 2.3 切向动量协调系数的实验结果

---

Page 24

第二章 滑移边界条件概述

- 14 -

2.2.1.2 一阶滑移边界条件

当忽略壁面附近气体切向温度梯度时，式（2.6）简化为

v

v

s

2

U

U

u us w

u

n

n s

s

l

x

s

骣

骣

-

抖　鼢

珑

鼢

珑

-

=

鼢

珑　鼢

珑桫

桫

抖

=

（2.10）

这就是所谓的一阶滑移边界条件。

2.2.1.3 耦合滑移边界条件

当气体在壁面附近存在温度梯度时，将引起气体在近壁面处产生由低温向高温方向

的流动，称之为热蠕流(thermal creep)，它会影响壁面处的滑移速度。气体因温度梯度存

在而产生附加的宏观速度为

s

s

3 j

T

u

4

T x

m

r

骣

+

÷

ç ¶÷

ç

÷

=

ç

÷

ç

÷

ç ¶

桫

（2.11）

热蠕流将影响到壁面处的滑移速度，因此当壁面处既存在速度梯度又存在温度梯度

时，应综合考虑两者的作用，则此时滑移边界条件变为：

s

w

u

s

s

s

3 j

u

T

u

u

4

n

T x

m

x

r

骣

骣

÷

+ ç

抖 ÷

÷

ç

ç

÷

÷

ç

ç

-

=

+÷

÷

ç 麋

ç

÷

桫

÷

ç

抖

桫

（2.12）

其中 j 气体分子内部自由度数目，对于常温下的单原子气体分子，j＝0，对于内部

仅有转动自由度的双原子气体分子，j＝2。

2.2.1.4 高阶滑移边界条件

图 2.4 高阶滑移边界条件的分析模型

Beskok 等人[38]在研究滑移流动问题时提出了高阶滑移边界条件的概念，如图所示。

---

Page 25

北京航空航天大学硕士论文

- 15 -

假设在距离壁面一个气体平均自由程高度的地方存在一个假想的层面，壁面附近的气体

分子中一半来自该假想层面，且具有平均的切向速度 u

l ，而另一半来自壁面的反射。

根据 Schaaf 和 Chambre[35]的有关结论可知气体在近壁面处的速度为：

(

)v

v w

s

1

u

u

1

u

u

2 l

l

s

s

轾

=

+ -

+

犏

臌

（2.13）

将式(2.4)中的 u

l 按内法线方向 n 在壁面处作 Taylor 展开，整理得

2 2

3 3

v

2

3

s

w

v

2

u

u

u

u

u

n 2

6

n

n

s

l

l

l

s

轾

-

抖

犏

-

=

+

+

+ L

犏¶

抖

臌

（2.14）

上式就是滑移边界条件得高阶表达式。高阶滑移边界条件得引入，降低了由于 Kn

数增加而引起得滑移速度增加的幅度。

2.2.2 温度跳跃边界条件

2.2.2.1 温度跳跃边界条件的导出

图 2.5 温度跳跃区壁面附近的温度分布

采用与滑移速度导出类似的方法，可以导出壁面处得温度跳跃。假设气体无宏观流

动且表面附近仅有法线方向的温度梯度，在此种情况下，Knudsen 层外缘处向物体表面

运动的气体分子的速度分布可表示为：

(

)

2 2 2

1

1

0

2

5

/

/

v v v

k

T

x y z

f

f

vy

T m

y

s

T ms

κ

ρ κ

⎫

⎧

⎡

⎤

⎪

+ +

−

−

⎪

⎪

⎛

⎞

∂

⎢

⎥

⎬

=

+

× −

⎨

⎜

⎟

⎢

⎥

∂ ⎪

⎝

⎠

⎪

⎢

⎥

⎣

⎦

⎪

⎩

⎭

（2.15）

其中， Ts 为 Knudsen 层外缘处气体的温度， 0f

−

是 Maxwell 速度分布函数：

(

)

2

2

2

exp

0

3/2

2

/

2

/

v

v

v

n

x

y

z

s

f

T m

s

T m

s

κ

πκ

⎡

⎤

+

+

⎢

⎥

−

=

−

⎢

⎥

⎢

⎥

⎣

⎦

（2.16）

---

Page 26

第二章 滑移边界条件概述

- 16 -

因此，从 Knudsen 层外缘入射到表面的分子通量可由下式计算：

0

T

s

v f dv dv dv

n

y

x y z

s 2 m

k

y

p

ゥ

-

= -

=

蝌

- ? ?

（2.17）

而入射到表面的气体分子所携带的能量则由下式计算：

(

)

0

2

2

2

y

x

y

z

s

x

y

z

1

j

E

v m v v v

T f dv dv dv

2

2

k

¥

¥

-

-

- - ?

轾

犏

= -

+ +

+

犏

臌

ò 蝌

(

)

0

2

2

2

y

x

y

z

x

y

z

s

j

1

T

v

m v v v f dv dv dv

2

2

y k

¥

¥

-

- - ?

骣

轾

÷

ç

÷

ç

犏

=

-

+ +

÷

ç

÷

ç桫

犏

臌

ò 蝌

s

s

s

1

j

T

2 T

T

k

2

2

y

y k

k

骣

骣

÷

ç ¶

÷

ç

÷

ç

÷

÷

ç

=

+

+

ç

÷

÷

ç

÷

ç

ç

÷

桫

ç ¶

桫

（2.18）

此时设从表面漫反射出去的分子的温度为

w

T ，则漫反射分子所携带的能量为：

(

)

(

)

2

2

2

x

y

z

2

2

2

y

x

y

z

x

w

2

y

z

w

w

0

v v v

1

j

E

v

m v v v

exp

dv dv dv

T

2

2 T /m

2

2

T /m

y

y k

k

p k

¥

ゥ

+

-

-

轾

骣

+ +

轾

÷

ç

犏

÷

ç

犏

=

+ +

-

+

÷

犏

ç

÷

ç桫

犏

犏

臌

臌

ò 蝌

w

w

j

2 T

T

2

y k

k

骣

÷

ç=

+

÷

ç

÷

ç桫

（2.19）

入射分子所携带的能量与反射分子所携带的能量之差即为 Knudsen 层中的净能量

交换，而净能量交换则等于以热传导方式传向壁面方向的热流。

为了反应入射的气体分子与壁面之间能量交换的完善程度，定义热协调系数

T

s (thermal accommodation coefficient):

i

r

T

i

w

E E

E E

s

-

=

-

（2.20）

式中

i

E 和 r

E 分别表示入射于物面的分子与从物面反射回气体流场的分子所具有的

平均能量，

w

E 表示当反射分子的温度等于实际壁温 w

T 时气体分子所具有的平均能量。

显然，当入射分子与壁面间的能量交换完全充分，即入射分子所携带的能量全部交给壁

面时，

r

w

E E= ，此时， T

1

s = ；而当入射分子与壁面之间完全不发生能量交换时， r

i

E E= ，

因此

T

0

s = 。由此可见，热协调系数 T

s 的实际数值总是处于 0 与 1 之间，具体数值取决

与气体的种类、温度、压力、壁面的材料、光洁度及清洁度等因素。对于常用的工程材

料表面，通常取

T

0.8

s =

。

由于热协调系数

T

s 的引入，从而可得到能量交换的等式为：

把气体当作连续介质，关于它的描述是以流动参数，如速度、密度、压力和温度在空间

 

和时间的变化为依据。N-S 方程就是连续介质的数学模型。微观模型或分子模型是把气

 

体当作无数单独的分子，提出其每一时刻的位置、速度和每一分子的状态，在这一层次

 

上的数学模型为 Boltzmann 方程

和时间的变化为依据。N-S 方程就是连续介质的数学模型。微观模型或分子模型是把气

 

体当作无数单独的分子，提出其每一时刻的位置、速度和每一分子的状态，在这一层次

 

上的数学模型为 Boltzmann 方程

把气体当作连续介质，关于它的描述是以流动参数，如速度、密度、压力和温度在空间

和时间的变化为依据。N-S 方程就是连续介质的数学模型。微观模型或分子模型是把气

体当作无数单独的分子，提出其每一时刻的位置、速度和每一分子的状态，在这一层次

上的数学模型为 Boltzmann 方程