分子气体稀薄效应的动理学建模
2022-05-10曾嘉楠
曾嘉楠，李 琪，吴 雷
(南方科技大学 力学与航空航天工程系，深圳 518055)
0 引 言
经典流体力学与空气动力学建立在连续介质假设之上，即宏观微元中需包含足够多的微观分子。与之相反，稀薄气体动力学描述的则是在此假设不成立时气体的流动与热力学等性质。克努森数Kn是稀薄气体研究中一个重要的无量纲参数，用来表征气体的稀薄程度，通常定义为气体分子平均自由程与流动特征长度的比值。因而，当分子平均自由程较大（即气体密度与压强较小）时，或在流动特征长度极小（即微纳米尺度下），气体流动均可呈现出显著的稀薄效应。
由此，不难发现，稀薄气体流动广泛存在于航空航天、工业技术和能源环境等重要的国防、民生领域，如返回舱再入大气层以及临近空间飞行器高空域飞行时的绕流[1-5]、微机电系统中的气体流动[6-9]以及页岩气开采中甲烷在复杂多孔介质中的流动[10]等。相关科学技术的发展对我国国防、国民经济和社会发展意义重大，如下所述：
1）临近空间飞行器既可以避免绝大部分的地面攻击，同时也能够有效实施对地攻击和对航天器的打击，因此临近空间是进行空中军事活动的理想区域，发展潜力极大，我国和不少其他国家已经将其纳入信息化武器装备体系[11]。在70 km以上的飞行空域，复杂构型的高超声速飞行器表面存在极为复杂的流动，因而需要建立新的流动模型和计算方法，来预测复杂的稀薄气体效应、高温真实气体效应以及解决激波/膨胀波与黏性作用的耦合干扰等问题。深入研究上述实际应用中的稀薄气体流动规律是优化飞行器设计的关键所在，具有重要意义。
2）机电系统微小集成化是提高能源利用效率的重要手段，因为微小化可大幅度减少传热传质阻力从而显著提高工质的相关输运能力。近年来，微流动系统在电子、化学、生物等领域已经取得了巨大的发展并得到了广泛的应用，但是微尺度下的空气动力学理论和模拟方法仍不完善，还存在许多亟待解决的问题[12]：如在电容式膜片压力计中，当被测气体的温度与压力计的温度不同时，稀薄气体自动从低温区向高温区蠕动，从而导致压力测量偏差，需要考虑稀薄气体效应进行修正[13]；工作在低压风洞的测量设备也需要考虑相应的稀薄气体效应修正以提高测量精度。
3）随着自然资源逐渐枯竭和环境破坏愈发严重，能源转型和能源利用效率提升成为各国科技发展的重要目标。相较于传统化石能源，以甲烷为主要成分的页岩气，因其具有弥补常规能源储量及产量、减少温室气体排放的巨大潜力，在世界能源供给中变得愈发重要。值得一提的是，我国页岩气储备位列全球第一，但复杂的地质条件要求对其输运机制有更清晰深入的理解，才能实现经济有效的开采；虽然页岩气埋藏在地表深处，气体压力大，但是页岩孔隙直径在纳米尺度，因而同样需要考虑稀薄气体效应[14]。
图1列举了典型稀薄气体流动中纳维-斯托克斯方程的数值预测结果与实验测量存在明显差距的例子。由此可以看出，基于连续性假设得到的流体力学方程不再适用于此类问题，其根本原因在于，在稀薄效应较强的区域内，由宏观量的一阶导数构成的线性本构关系不再准确。另一方面，直接模拟所有气体分子运动与相互作用的分子动力学方法，提供了一个更为准确地模拟稀薄气体流动的途径。然而，因其庞大的计算量，该方法仅能用于极小时空尺度下的气体输运问题。因此，对于更为广泛存在的介观尺度下稀薄气体流动，如何建立准确的动理论模型方程并有效求解，则是目前最为关切与亟待解决的两个问题。
图1 典型稀薄气体流动及纳维-斯托克斯方程的预测误差：第一列图显示，当马赫数大于2时，纳维-斯托克斯方程低估正激波的厚度[15]；第二列图显示，当声波振动频率接近气体分子碰撞频率时，纳维-斯托克斯方程低估声波接收器的声强[16-17]；第三列图显示，低压下纳维-斯托克斯方程可能远远低估微纳米通道的质量流量[18]Fig. 1 Typical rarefied gas flows and the incapability of Navier-Stokes equations. Figures in the first column show that, when the Mach number is larger than 2, the Navier-Stokes equations underpredict the thickness of normal shock wave[15]. Figures in the second column show that the Navier-Stokes equations underpredict the sound pressure at the receiver[16-17]. Figures in the last column show that the Navier-Stokes equations might underestimate the mass flow rate significantly[18]
（1）研究现状。
一般来说，传统的纳维-斯托克斯方程只能描述连续流区（Kn<0.001）的气体流动。随着克努森数的增加，气体流动将分别经历滑移流区（0 .001
作为最为重要的气体动理论方程，玻尔兹曼方程描述了单原子气体的速度分布函数在时间与空间中的演化规律。玻尔兹曼方程的求解，在航空航天、微机电系统和页岩气开发等领域中都发挥着重要作用。其中，基于统计随机方法的DSMC计算，是模拟稀薄气体流动的重要工具[24]。对于单原子气体的模拟，DSMC方法已在数学上被证明在粒子数趋于无穷的条件下，收敛于玻尔兹曼方程的解[25-26]。玻尔兹曼方程中气体分子的速度分布函数是定义在平均自由程尺度上，由于DSMC方法将气体分子自由运动与相互碰撞这两个过程进行了解耦，因此该方法的模拟中要求空间网格尺寸小于气体平均自由程，且时间步长小于平均碰撞时间。而这直接导致了DSMC方法对近连续流模拟的低效性。此外，由于流场的宏观量信息需要通过统计抽样得到，因此DSMC模拟中需进行大量系综平均去除噪声，特别地，当宏观物理量自身数值相对于其统计微观量的涨落较小时，例如低速流动或较小温差的系统，DSMC方法的时间成本将变得难以负担。
相比于DSMC等统计模拟方法，确定性数值方法在近几十年得到快速发展并展现出强大优势，我国学者在此方面贡献良多。但由于玻尔兹曼方程碰撞项的复杂性，除了少数方法可以求解玻尔兹曼方程外[27-31]，大部分确定性方法求解的是简化模型方程，如BGK模型[32]、ES-BGK模型[33-34]和Shakhov模型[35]等。中国空气动力研究与发展中心的李志辉和张涵信[36]最早应用气体动理论统一数值算法求解模型方程，并建立了三维复杂飞行器高马赫数绕流问题的统一算法大规模并行计算研究平台，已应用于返回舱回收和天宫一号离轨、再入、解体等航天工程[37-39]。香港科技大学徐昆首创的统一气体动理论格式（unified gas kinetic scheme，UGKS）耦合处理了模型方程的迁移项和碰撞项，相对于传统离散速度方法而言，在网格离散尺度远大于分子平均自由程时依然可以在连续流时恢复纳维-斯托克斯方程的解，从而在计算多尺度问题时可以采用更粗粒化的时空离散[40]。华中科技大学的郭照立、西北工业大学的钟诚文、新加坡国立大学舒昌等课题组在UGKS的基础上做了进一步发展[41-43]。最近，苏微和吴雷等提出的合成迭代算法（general synthetic iterative scheme，GSIS）耦合了宏观方程以实现全流场信息交换，具有良好的“渐近保持”和“快速收敛”性质，不但可以使用大空间网格和时间步长，而且可以在短短几十步迭代内找到稳态解，极大地提高了数值精度和效率[44-47]。更重要的是，该方法对碰撞项的具体形式没有要求，可以直接针对玻尔兹曼方程进行快速求解。
（2）存在的一些问题。
稀薄气体研究主要涉及动理学建模和多尺度计算两个方面，本文仅关注分子气体动理论模型中存在的一些问题。众所周知，目前大部分气体动理论方程是针对单原子气体提出的[48-50]，而实际气体大多是分子气体，即一个气体分子包含两个或两个以上的原子，例如占据地球大气绝大部分比例的N2和O2，占据火星大气绝大部分比例的CO2等。相对于仅具有平动自由度的单原子气体，分子气体中还存在转动与振动的内能形式，并且各类型自由度的能量可以相互转化。对于常见气体，其转动能在室温下已充分激发，而振动能也在上千度高温下逐渐激发。随着温度的继续升高，分子气体的离解与电离也将逐渐发生，与此同时，化学反应也通常存在于此。此时，高温真实气体效应显著。而单原子玻尔兹曼方程无法包含这些过程，因此需要更复杂的模型方程来给出相应的模拟。目前，大部分以此为目标的建模工作是针对连续流提出的，例如对高超声速飞行器周围气体的模拟。由于分子气体的热流由平动能量和内部能量的输运与相互转换共同作用导致，因此需要通过额外的能量方程来描述高温气体能量输运，并且保证能量交换速率和总热导率与实验数值一样，此类连续流下的模拟方法就可以较为准确地模拟高温真实气体效应，国内外学者在这方面已经做了大量工作方面已经开展了大量工作[51]。
然而，对于临近空间飞行器等类似问题而言，稀薄气体效应也非常显著[19]，因此有必要将稀薄气体效应和高温真实气体效应结合起来考虑。对此，王承书和乌伦贝克建立的WCU方程[52]给出了相对完备的描述，该方程以经典力学方式描述平动能量，而以量子力学方式描述内部（转动与振动）能量，并对每个能级下的速度分布函数建立方程。显然，该方法使得方程组自由度过于庞大，致使其求解变得相当困难从而难以应用于实际问题。因此建立相对简化的模型方程尤为必要。与此同时，由于各类自由度的特征碰撞时间不同，导致稀薄气体效应的程度也不同，因此动理论建模难度也大为增加。
对于分子气体的动理论建模， 如何在模型方程中准确反应气体的各个输运系数是首要问题。相对于单原子气体，由于更多内部自由度的存在，分子气体中出现了更多的输运过程及相应参数。其一为体积黏性，其表征了气体在膨胀压缩的做功过程中能量与内部自由度转换的阻力，这一阻力来源于有限时间的能量弛豫过程；其二为转动、振动热导率，即内部自由度同时参与了热流的传递，与热流的弛豫速率紧密相关。实际上，在分子气体动理论的建模中，正确的体积黏性是容易实现的，只需保证正确的内部能量（转动、振动等）和平动能量的弛豫速率即可。但是，实现正确的平动、转动和振动热导率却非常困难。本文以仅有转动自由度的分子气体为例说明问题所在。历史上，Eucken首先将气体总热导率 κ与剪切黏性 µ的关系表示为：
其中：m是 分子质量；kB是玻尔兹曼常数；dr是分子的转动自由度的数量，如N2在常温下具有两个转动自由度，则dr= 2； κtr和 κrot分别是气体的平动和转动热导率；ftr、frot和feu分别是平动、转动和总Eucken因子。对于单原子气体，ftr非常接近2.5，但对于分子气体，ftr和frot的数值很难确定。最初Eucken从简单分子动理论出发，认为：
后来，人们意识到气体内部热导率与分子扩散有关，将内部Eucken因子改写为：
其中 ρ为气体密度，D为气体自扩散系数，Sc为斯密特数。最后，Mason和Monchick根据王承书和乌伦贝克建立的WCU方程[52]，推导出Eucken因子的表达式[53]，但与实验数据比较仍然存在差距[54]。由此可知，目前并不存在简明准确的仅依赖于气体参数与其他输运系数的各个Eucken因子的表达关系，因此各热导率分量均需直接与实验数据进行匹配，才能保证相应输运系数的准确性。然而，虽然可以比较容易地测得总热导率，但从实验获取热导率各分量本身便是相当有难度的任务。
WCU方程的复杂性限制了其在实际问题中的应用，而DSMC方法则相对容易地从单原子气体模拟推广至分子气体的输运。即便如此，对于极为复杂的分子气体碰撞过程，DSMC方法依然采用了适当简化的模型来处理，从而提高计算效率。目前通用的分子碰撞模型最早由Borgnakke和Larsen提出[55]。此后，人们提出各种各样的改进，主要目的均在于正确实现分子各能级间的能量交换速率，从而获得正确的分子气体的体积黏性[56-58]。然而，当前的DSMC方法并没有包含可以单独调节热导率的相应机制与参数。与此同时，在实际中使用DSMC模拟分子气体流动时，该碰撞模型对应得到的平动和转动热导率是否正确尚未可知。因此，有必要对DSMC模拟分子气体的相应模型进行检查，这里的检查有双重意义：1）检查DSMC模拟得到的总热导率是否与实验测量数据一致，若不一致，则DSMC即便在近连续流下也会给出错误解。2）即使DSMC能给出正确的总热导率，也还需要检查ftr和frot的数值是否与实际气体一致。
在近连续流条件下，能量在平动自由度与内部自由度的分配趋于一致，因而平动温度与转动和振动温度的差别极小，此时保证总体热导率的正确性，即可得到相应正确的结果。因而，在稀薄效应较强的区域，保证各个热导率分量的正确性就显得尤为关键。例如，爱因斯坦等很早就发现，对于热蠕动等稀薄气体流动，温度梯度引起的气体流动只与ftr有关[59]。最近李琪和吴雷等发现，目前流行的、公认的DSMC并不一定能做到以上两点[54,60]，因此DSMC在分子气体流动中的模拟结果并不可靠。同样地，对于确定性方法，人们在原有单原子气体模型方程的基础上，提出了各种描述分子气体稀薄效应的动理论简化模型方程，如ES-BGK模型[33-34]、Rykov模型[61]等。与DSMC方法相似的是，这类模型对于分子气体的拓展同样专注于多种能量相互转换的速率（即气体体积黏性的正确性），而忽略了热流弛豫过程的速率（即热导率分量的正确性）。具体而言，除了Rykov模型可以独立调整平动和内部Eucken因子以对应真实参数外，其他模型只能调整总热导率，因此不适合用于研究热蠕动等稀薄气体现象。
（3）本文目的和结构。
分子气体稀薄效应的准确刻画是稀薄气体动理学和高温气体动力学研究的前沿和热点，在这个背景下，本文首先介绍分子气体动理论建模的思路和最新研究进展，并讨论现有的分子气体模型方程在稀薄气体流动中的适用性与准确性，为今后多组分分子气体和化学反应的非平衡流动的建模和模拟奠定坚实的理论基础。以下第1节介绍单原子气体玻尔兹曼方程的基本性质、输运系数的严格推导以及动理学建模的基本思路。特别地，我们指出，几乎所有的简化模型都可以从Gross-Jackson模型非线性化而来。第2节介绍王承书和乌伦贝克提出的WCU方程，推导分子气体的体积黏性和平动、转动热导率系数，并指出DSMC碰撞模型的缺陷。第4节介绍WCU方程的线性化Hanson-Morse模型以及主流的非线性模型，并在经典的稀薄气体流动中测试模型的精度。第4节和第5节分别讨论热弛豫速率的不确定性对稀薄气体流动的影响，以及从分子动力学模拟和实验中减小甚至消除热弛豫速率不确定性的方法。最后做总结和展望。
1 单原子气体玻尔兹曼方程及其简化模型
在气体动理论中，单原子气体在相空间中的概率密度由速度分布函数f(t,x,v)表 示，是时间t、空间坐标x=(x1,x2,x3) 和 分子速度v=(v1,v2,v3)的函数。定义f(t,x,v)dxdv是 体积为 dxdv的相空间上的分子数，则气体的数密度n、 宏观速度u、 温度T、压力张量pij和热流q可以分别通过对速度分布函数求矩得到。
其中，c=v−u是气体分子热运动速度，即气体分子速度与当地宏观速度的矢量差，而c是热运动速率。注意本文其他矢量也如此表示：黑体表示矢量，斜体表示矢量的模。定义应力偏量为：
其中，p=nkBT， δij为克罗内克函数。
对于稀疏气体（分子间距远远大于分子直径），在外部施加的加速度a=(a1,a2,a3)的作用下，速度分布函数的演化由玻尔兹曼方程描述：
方程（6）左边的三项分别表示速度分布函数在时间上的变化、在速度作用下在物理空间的变化以及在外力作用下在速度空间的变化，右边表示使得分布函数趋于平衡态的气体分子的碰撞过程。在玻尔兹曼方程中，两体碰撞的形式为：
其中，v和v∗分 别是碰撞前两个分子的速度，而v′和v′∗是碰撞后的速度。由于碰撞前后两分子的距离足够远，因此其相互作用可以忽略不计。根据动量和能量守恒定律，碰撞前后速度关系如下：碰撞示意图见图2。其中，碰撞前两分子的相对速度为vr=v−v∗， 碰撞后的相对速度为v′−v′∗。 Ω为定义在单位球空间的矢量，其与碰撞后的相对速度同方向。于是相对速度的偏转角 θ与碰撞前相对速度满足如下关系：
图2 (左)两体碰撞前后的速度分布：由于动量和能量守恒，碰撞前后的相对速度分布在球体上并且通过球心；(右上)中心力场作用下的经典两体碰撞示意图，其中b 为 瞄准距离，k为沿两分子间最短距离方向的单位矢量；(右下) 直径为σ的硬球分子的两体碰撞Fig. 2 (Left) Velocity redistribution after the binary collision.Due to the conservation of momentum and energy, the pre- and post-collision relative velocities fall in a sphere and pass through the sphere center. (Top right) Classical binary collision between molecules with central force, where b is the aiming distance, andk is the unit vector along the minimum distance between two colliding molecules. (Bottom right) Binary collision between HSmolecules of the diameter σ
最后，碰撞核B(θ,vr)是碰撞偏转角度和相对速度的函数，具体形式取决于分子间的作用力。
1.1 偏转角和微分散射截面
假设分子间通过中心力场作用，作用势 φ(r)已知，其中r为分子间距，则偏转角可以通过经典力学和量子力学两种方式求解。若气体温度不低（如氦气的温度高于100 K），两种方式得到的输运系数（如黏性和热导率）相同[62]。这里仅介绍经典力学的计算结果，如图2所示，偏转角可以表示为：
其中，W=b/r为瞄准距离b与 分子间距r的比值，积分上限W1对应于最短分子间距，即式（10）中括号内表达式等于零的方程的正根。
在气体动理论中，经常考虑如下形式的逆幂律分子作用势：
其中，K表征分子间相互作用的强度。从式（10）可知，偏转角只与变量s有 关，即 θ=Θ(s) ：
微分散射截面定义为：
而碰撞核为：
当式（11）中 η=5时，即为麦克斯韦分子，此时碰撞核与相对碰撞速度无关，记为对于硬球分子模型，可看作式（11）中 η=∞。从图2可以看出，偏转角可通过如下公式确定：b=σcos(θ/2) ，其中 σ为硬球直径。因此微分散射截面为 σD=σ2/4， 而碰撞核为B(θ,vr)=σ2vr/4。对于一般的气体，其行为介于麦克斯韦分子和硬球分子之间。
1.2 熵增原理
对于任意关于分子速度的函数Ψ (v)，对玻尔兹曼碰撞项（式（7））在速度空间积分，具有如下对称性：
其 中 Δ[Ψ]=Ψ(v∗)+Ψ(v)−Ψ(v′∗)−Ψ(v′)。 若 Ψ满 足则称其为碰撞不变量。根据质量、动量和能量守恒，可知Ψ =1、v、v2为碰撞不变量，而碰撞不变量的线性组合也是碰撞不变量。
定义H函数为：
H是与气体系统的熵相关的标量。在无外力的情况下，H的时间导数可写为：
式中，n为 系统表面微元 dS的外法线方向。对于孤立系统，上式右端第一项为零；利用方程（15），右端项第二项可写为：
因为碰撞核B非负，且对于任意的两个正整数a和b， 不等式 (a−b)(lna−lnb)≥0恒成立，所以有：
式（17）表明孤立系统的H函数不会减小，这就是著名的玻尔兹曼的熵增原理。H随时间单调增加，但存在有限的上界，上界对应式（17）中等号成立时的平衡态，即 lnf(v∗)+lnf(v)=lnf(v′∗)+lnf(v′)。 因此， lnf也是碰撞不变量，可以表示为五个基本碰撞不变量的线性组合： lnf=α1+α2·v+α3v2；给定系统的密度、速度和温度，参数α1、α3和α2可以唯一确定。此时，麦克斯韦平衡态速度分布函数为：
1.3 线性化玻尔兹曼方程
线性化玻尔兹曼方程在气体动理论中具有重要地位。第一，线性化玻尔兹曼碰撞项的本征值与本征函数，不仅在渐近展开推导纳维-斯托克斯方程的过程中至关重要，而且是发展简化动理学模型方程的源头理论。第二，在许多微机电系统中，气体的压力梯度与温度梯度非常小，因此使用线性化方程可以高效准确地模拟微流动。第三，虽然玻尔兹曼方程是单粒子速度分布函数确定性演化的平均方程，但是在某些问题中（例如：第5.2节中的瑞利-布里渊散射）可以通过线化玻尔兹曼方程研究粒子涨落带来的影响。
通常将速度分布函数在全局平衡态
下展开为：
其中扰动量 φ满足约束条件 |φ|≪1。不考虑外力项，只保留 φ的一次项，玻尔兹曼方程（6）可线性化为如下形式：
这里的速度和空间坐标分别通过最概然速率vm和特征尺寸 ℓ进行无量纲化，时间用 ℓ/vm归一化；无量纲分子速度为：
其中：
为气体分子的最概然速率，即与麦克斯韦速率分布的极大值对应的速率。
对于麦克斯韦分子，线性化玻尔兹曼碰撞项的本征值与本征函数为[52]：
式中，Pl(x) 是勒让德多项式，是索南多项式，为 ξ方向的球谐函数。
前三项本征值为 λ00=λ01=λ10=0，分别代表三大守恒律。另外两个非常重要的本征值为：
这两个本征值决定了剪切黏性系数和热导率的大小。若归一化的速度 ξ 在z轴 的投影为 ξz=ξcosθ（其中 θ为极距角），与上述本征值对应的本征函数为：
其与气体的分子数密度、速度、温度、压力和热流的扰动密切相关：
式（26）中，下标中的尖括号表示该张量是无迹张量，即ξ〈iξj〉=ξiξj−ξ2δij/3。为避免使用过多符号，在线性化问题中，本文使用与原物理量相同的符号表示归一化的扰动物理量，如这里的扰动数密度n应理解为偏离参考数密度n0的那部分，再除以参考数密度；扰动温度T为偏离参考温度T0的部分，再除以参考温度；在全局平衡态下，参考速度、应力偏量和热流全部为零，因此上式中它们分别以T0下 的最概然速率vm、n0kBT0和n0kBT0vm归一化。
1.4 输运系数和弛豫率：表象与本质
玻尔兹曼方程和纳维-斯托克斯方程分别在介观和宏观尺度上描述气体系统的演化，它们之间的关系一直是数学家和力学家研究的重要课题[63]。对玻尔兹曼方程求守恒量m、mc、mc2/2的矩，可以得到如下宏观方程组：
其中E=3kBT/2m+u2/2。但是该方程组并不封闭，因为压力张量和热流的表达式不能用密度、速度和温度表示。有三种主要方法尝试封闭宏观方程组，它们分别是希尔伯特展开法[64]、Chapman-Enskog展开法[21,65-66]和矩方法[22,67]。这里介绍最常用的第二种方法。
在Chapman和Enskog[65-66]方法中，首先将玻尔兹曼方程改写为：
然后将分布函数和玻尔兹曼碰撞项展开为 ε的幂级数：
值得注意的是，ε是一个小的形式参数，主要用于监测展开的阶数。展开完成后，取 ε=1。
同时，压力和热流也相应展开如下：
但是，碰撞不变量对应的宏观量CM={ρ,u,T}仅由分布函数的零阶展开决定，即：
而各个高阶展开对CM的贡献均为零，这称为兼容性条件。
将式（30）代入式（27），可知时间导数亦可表示为ε级数[67]，
且对于零阶近似，可得到：
从式（28）和式（29）易知Q(f(0))=0。因此，速度分布函数的零阶展开为：
而速度分布函数的一阶近似f(1)=Feqφ满足以下方程[67]：
其中，方程的推导过程中用到
以及式（33）。注意式（21）中的feq应改写为Feq。
积分方程（35）的解可分为齐次部分与非齐次部分。齐次部分满足J(φ)=0， 这说明 φ必然是碰撞不变量的线性组合，而兼容性条件要求所有线性组合的系数都为零。非齐次部分满足如下形式：
其中，A(ξ) 和B(ξ)的解满足：
一旦得到A(ξ)和B(ξ)，通过一阶展开就可以恢复牛顿黏性定理与傅里叶热传导定理：
且剪切黏性与热导率分别为：
一般将A(ξ) 和B(ξ)展开为索南多项式级数：
然后通过式（37）和索南多项式的正交性求解输运系数。对于麦克斯韦分子，取第一项展开可得到准确的剪 切 黏 性 和 热 导 率，即 取A(ξ)=a1(ξ2−5/2)和B(ξ)=b0；而对于其他相互作用势，仅考虑第一项展开得到的输运系数的相对误差也不超过2%。定义等效黏度截面为：
则剪切黏性为：
而热导率为：
从而单原子气体的普朗特数为：
对于逆幂律分子间相互作用势（式(11)），可以得出气体黏性和温度的关系：
其中 ω是黏性温度幂指数。对于麦克斯韦分子和硬球分子，ω分别取1和0.5；对于其他气体，ω一般介于0.5和1之间。
将式（38）代回式（27），即可获得纳维-斯托克斯方程。若继续计算Chapman-Enskog展开的高阶近似，可以得到Burnett、super-Burnett等宏观方程组[21]。但是可能由于应力和热流与宏观量CM展开的不匹配，导致高阶方程组的线性不稳定性。同时，由于边界条件数量随着方程组阶数升高而增加，但对附加的边界条件提法缺少充分的研究，因此高阶方程组很少被实际应用。另一方面，即使这些高阶方程组在某些情况下能得到精确解，解的精度也不一定随着展开阶数的增加而增加[17]。
值得一提的是，对于麦克斯韦分子，在空间均匀系统中，本征值 λ02和 λ11与应力偏量和热流的弛豫速率相关，而这些弛豫速率决定着剪切黏性系数和热传导的大小：
可以认为这些弛豫过程是本质，而黏性系数和热导率则为表象：在稀薄气体流动中，本构关系式（38）和等效输运系数会随着克努森数的改变而改变，但一旦分子作用势确定，弛豫系数就会确定不变。
利用弛豫时间和输运系数的关系，可以大大简化模型方程中输运系数的推导过程。即，如果应力偏量的弛豫时间为τ，则剪切黏性为pτ。如果热流的弛豫时间为 τ/A， 则普朗特数为A。
1.5 模型方程
玻尔兹曼方程碰撞项的复杂性给数值求解带来极大挑战[31]，因此，学者建议使用简化的碰撞模型来代替。在简化玻尔兹曼方程的碰撞项时，一般认为需要满足以下几点：
1）动理学模型必须保证质量、动量、能量守恒；
2）理想气体的局部平衡态速度分布函数为麦克斯韦分布；
3）从模型方程导出的黏性系数与热导率应与玻尔兹曼方程导出的结果保持一致；
4）满足熵增定理，当且仅当孤立系统处于平衡态时熵增为零，并达到其最大值。
其中，前两项是基本要求，第三项要求在连续流区域通过Chapman-Enskog展开恢复出纳维-斯托克斯方程。而第四项仅为从物理角度出发的考量，并不能保证模型方程的精度：如果一个动理学模型方程与玻尔兹曼方程完全一致，那么熵增速率也应该相同，然而这通常是不可能做到的。实际上，大多数时候，仅在线性化情况下满足熵增定理的Shakhov模型的精度优于完全满足熵增定理的ES-BGK模型。
由于玻尔兹曼碰撞项在形式上可记为Q=Q+−f/τ，因此，动理学模型方程的碰撞项一般采取如下形式：
其中：fr为参考速度分布函数；τ为特征碰撞弛豫时间，表征速度分布函数趋向参考分布函数的快慢，为简单计，通常假设与分子速度无关。
1.5.1 BGK模型
BGK模型是最著名的气体动理学简化模型，由Bhatnagar、Gross和Krook提出[32]，其参考速度分布函数是局部麦克斯韦分布：
易见此模型满足三大守恒律，且给出了在平衡态下的正确解f=Feq。在空间均匀系统中，由于 lnFeq是碰撞不变量的线性组合，有，所以可以证明BGK模型满足熵增原理：
空间均匀系统中应力偏量和热流的弛豫时间均为τ，因此根据式（46）随后段落的描述，可知剪切黏性为 µ=pτ，而普朗特数为1。然而，通过玻尔兹曼方程推导得到的普朗特数正确值为 2/3。这说明BGK模型不可能同时得到正确的剪切黏性和热导率。对于黏性主导的流动，一般选取 τ以恢复剪切黏性；对于热传导主导的流动，则以恢复热导率为目标来选取τ。
1.5.2 ES-BGK模型
Holway[33]提出了ES-BGK模型以修正BGK模型中错误的普朗特数，其中参考速度分布函数是在给定质量、动量、能量和应力张量条件下，通过最大化熵函数式（16）而来：
其中：
若b=0， 二阶张量 λij中仅有主对角元素非零，模型恢复为BGK模型。
容易证明，应力偏量的弛豫时间为 τ/(1−b)，而热流的弛豫时间为τ。因此，剪切黏性和普朗特数分别为：
为确保弛豫时间为正，则有b<1。 另外，有界的充要条件是b≥−1/2。因此，ES-BGK模型的普朗特数范围为 [2/3,+∞)。为了恢复单原子气体的普朗特数2/3，取
可以看出，当碰撞项为零时，ES-BGK方程给出f=，当且仅当 σij=0时，速度分布函数退化为麦克斯韦平衡态分布。而这并不违背动理学建模的第二个要求，因为当局部平衡态为f=时，应力偏量的演化满足如下形式 ∂σij/∂t=−pσij/µ，即应力偏量将随时间推进逐渐衰减至零。因此，ES-BGK模型中趋于平衡态的过程可以解释为：当速度分布函数朝着椭球分布（式（50））松弛，椭球分布也朝着麦克斯韦分布演化，最终实现f=Feq。
Andries等[34]证明了ES-BGK模型满足熵增定理，使之成为唯一满足熵增定理且能正确恢复各输运系数的动理学模型，因此，ES-BGK模型受到广泛关注。
1.5.3 Shakhov模型
不同于ES-BGK模型在参考速度分布函数中引入应力修正项，Shakhov模型通过在参考速度分布函数中引入热流修正项来恢复正确的输运系数：
此模型同样满足三大守恒律。从偏应力和热流的弛豫时间可以看出，该模型的黏性为µ =pτ，且能实现正确的普朗特数。此外，由于热流弛豫过程满足∂q/∂t=−(2p/3µ)q， 在速度分布函数向frS靠近的过程中，参考速度分布函数也将逐渐演化为麦克斯韦分布，从而使Shakhov模型在平衡态时满足f=Feq。
线性化的Shakhov模型满足熵增定理。而对于非线性流动，熵增定理未被证明但也无从证伪；另外，速度分布函数可能出现非物理的负值。尽管如此，Shakhov模型仍能在许多问题中给出准确结果，因此得到了广泛应用。
1.5.4 Gross-Jackson模型及其非线性化
针对麦克斯韦分子，根据线化玻尔兹曼碰撞项的本征值与本征函数，Gross和Jackson提出的以下形式的动理学模型去逼近线性化玻尔兹曼碰撞项J(φ)：
式中， λst是 一个任意负数， λnl是线性化玻尔兹曼项的本征值（式（24）），θ′为v和v′的夹角。该模型的物理意义在于保证模型方程和线性化玻尔兹曼方程各阶矩的弛豫率相同。
Gross和Jackson将L(h)的N阶近似表示为：
即保留了所有阶数小于或等于N的多项式，并且取λst=λ0N，因此该截断模型的最大弛豫率为 λ0N。
在线性化情况下，即把速度分布函数在全局平衡态下根据式（20）展开，BGK、ES-BGK和Shakhov模型有如下形式：
其中扰动宏观物理量的定义可参考式（26）。它们都可以看作是Gross-Jackson模型的特例。当N=2时，从式（56）可以得到=LBGK；令N=3，根据式（24）有：
从而式（56）给出：
虽然，这既不符合线性化ES-BGK模型，也不符合线性化Shakhov模型，但是可将其视为这两个模型的线性组合。另一方面，若给定约束 2n+l≤3，当λst=−2p/3µ时，式（55）恰好展开为线性化ES-BGK模型；而当 λst=−p/µ时，式（55）展开为线性化Shakhov模型。
因此，通过Gross-Jackson模型和线性化BGK、ES-BGK、Shakhov模型的关系，可以反推出各种形式的非线性动理学模型。例如，陈松泽等将ES-BGK模型与Shakhov模型线性组合为一个复合模型，通过改变两模型所占的比例和普朗特数，复合模型能够分别恢复到BGK、ES-BGK和Shakhov模型式[68]。类似地，吴雷在博士论文也提出一个带有自由参数b的复合模型[69]：
该模型实际上是在Gross-Jackson模型（55）中取2n+l≤N=3和 λst=1/τ，并通过将如下的线性化碰撞项变为非线性碰撞项得来：
需要注意到，含有应力偏量的项既可以被吸收到指数函数里使参考分布函数变成类似式（50）的形式，也可以放在外面（此时对b的 限制仅为b<1），两者对数值计算的结果影响不大[69]。通过将Gross-Jackson模型中的其他阶本征函数非线性化为相应的形式，可以构造更加准确的动理学模型。
2 分子气体的特性
玻尔兹曼方程是仅针对单原子气体提出的，为模拟具有两个或多个原子的分子结构的气体的动力学行为，则需要在气体动理论层面重新建模。对于分子气体，两体碰撞保持动量守恒，但是能量却可以在各形式自由度和能级之间进行交换。因此，相比于仅具有剪切黏性的单原子气体，分子气体还存在体积黏性。另外，除了由分子平动引起的平动热导率外，还有由分子转动和振动引起的内部热导率。复杂的物理过程和更多的输运参数都给分子气体动理论建模带来挑战。
分子的内能可以表示为转动能量、振动能量以及电离能之和，分子气体内部自由度导致了更加复杂的碰撞项。王承书和乌伦贝克从量子力学角度出发建立了WCU方程，定义fi(t,x,v) 为 第i个 能量为ei的能级的速度分布函数，并对每个能级的分布函数都写出一个玻尔兹曼方程，其中迁移项保持不变，而碰撞项Qi是玻尔兹曼碰撞项的唯象扩展，
碰撞后的分子速度为：
碰撞核和跃迁概率的具体表达式是非常复杂的。这里以N2为例，在只存在转动自由度的情况下，对WCU方程的计算复杂性进行说明。在Lennard-Jones作用势下，转动能可表达为：
而转动能级的简并度为：
对照分子动力学模拟数据[70-71]，将跃迁概率拟合为：
其中，P0为一常数，显然，如果总能级数为N，式（61）的计算量为单原子玻尔兹曼碰撞项（式（7））的N4倍，计算时间成本巨大。
根据碰撞过程中平动动能是否守恒，WCU方程中的碰撞项可以分为弹性碰撞与非弹性碰撞。若处在i能 级和j能级的两个分子发生碰撞，碰撞后仍然能分别处在i能 级和j能级，则为弹性碰撞，否则为非弹性碰撞：
对于弹性碰撞和非弹性碰撞，可以分别引入两个等效的弛豫时间[72]：
而将WCU碰撞项形式地改写为：
这两个等效的碰撞时间，表征着能量在各种自由度或能级间的弛豫速度，因而与下面将讨论的分子气体输运系数、体积黏性、平动/转动热导率等，有着直接的关联。
2.1 体积黏性
体积黏性表征了气体无剪切膨胀或压缩时所受的阻力，对于稀疏气体，这种阻力来自于气体分子平动动能与内能的交换过程：在膨胀或者压缩过程中，压力做功会立刻改变平动动能，但在一定的延迟后才会通过非弹性碰撞来改变内能。为方便起见，在下述推导中，假设只存在一个转动能级。
引入非弹性碰撞所对应的平动能与转动能交换的弛豫时间 τr， 在空间均匀系统中，转动温度Tr的弛豫过程由Jeans-Landau方程描述：
式中，Tt和Tr分别为平动和转动温度，而平衡态温度T可表达为：
忽略剪切黏性和热传导的影响，通过膨胀压缩过程中能量守恒可得到：
式中， D/Dt为物质导数，推导得到：
将Tt−Tr基 于小量 τr展 开为Tt−Tr=Ttτr+O(τ2r)，可知式（67）左侧的时间导数项较右侧为更高阶小量，于是得到：
这说明，内能弛豫的影响将压力从原始压力p0=nkBT改变为平动压力pt=nkBTt。结合式（65）可得：
式中速度散度前的系数即为体积黏性 µb。在此，可以引入无量纲转动碰撞数Z：
其中由分子平动引起的平均碰撞时间定义为τ=µ/n0kBTt。于是体积黏性与剪切黏性的比值有如下关系：
从式（70）和式（71）可知，体积黏性正比于能量交换时间。对于单原子气体，若不考虑电离，不存在内部能量的交换，因此体积黏性为零。而对于一些分子气体，例如CO2，体积黏性可以是剪切黏性的上千倍[73]。但是，需要强调的是，式（68）成立的条件是τr很小，即能量弛豫的时间远远小于系统的特征时间。如果系统特征频率 ϖ过高，则等效的体积黏性一般与系统频率有如下近似关系[74-77]：
也就是说，在 ϖτr→∞极限下，体积黏性趋向零而不是无穷大，这是因为能量交换时间过长以至于系统无法感知，此时可以认为内部自由度是冻结的。因此，与式（46）中提到的剪切黏性和热导率的弛豫速率是本质而剪切黏性和热导率是表象一样，这里Jeans-Landau方程里的能量交换也是本质，而体积黏性只是连续流下的一个表象。
2.2 热导率
分子气体的热流由两部分组成，一部分是由分子的平动引起的动能传递，一部分是由分子的扩散引起的内能的传递。由于气体分子迁移过程中存在动能和内能交换，这一非弹性碰撞改变了能量在不同自由度下的分布，从而间接改变了动能传递与内能随扩散传递的结果。因而，平动热流和转动热流的弛豫过程是相互耦合的。于是，在空间均匀系统中，类似于式（46），平动热流qt和 转动热流qr的弛豫过程一般可以写成如下形式：
而通过Chapman-Enskog展开得到的连续流下的平动热导率 κtr和 转动热导率 κrot可表示为：
将方程（74）代入方程（1），就可以通过热弛豫系数求解如下方程，得到对应的Eucken系数：
基于WCU方程（61），Manson和Manchick通过修正的Chapman-Enskog展开推导出四个无量纲热弛豫系数的具体形式：
其中D′是有效扩散系数，如果极性分子发生强共振碰撞，该系数可能与自扩散系数D有明显差异，例如HCl的有效扩散系数满足 ρD′/µ=0.89， 而 ρD/µ=1.33。根据式（75），忽略1 /Z的高阶项，可将平动Eucken因子和转动Eucken因子写为[53]：注意到，ftr随着转动碰撞数Z的减小而减小，即内能交换越频繁，平动热导率降低。这是因为在分子气体的碰撞过程中，部分平动动能转化为内能，当D′增加，Z减小时，能量转化效率增强。这在物理上是合理的。
Mason与Monchick的解析解依然存在的一个问题，即无法保证通过式（76）推导而来的无量纲热弛豫系数A与实际气体参数保持一致。其原因在于，在黏性和扩散系数确定的条件下，转动碰撞数Z作为公式中唯一用来调节Eucken因子的参数，只能使得总热导率可能达到与实验数据相近的数值，但无法分别独立调整平动与内能Eucken因子。与此同时，由于Z也是决定分子气体体积黏性的唯一调节参数，参考式（71），则总热导率与体积黏性的准确性通常无法同时被满足。
2.3 DSMC中的热弛豫速率
对于单原子气体，在粒子数和空间网格足够多的情况下，DSMC已被证明收敛于玻尔兹曼方程的解[25]。实际上，由于DSMC中两体碰撞后粒子速度的更新方式与玻尔兹曼方程一致，式（43）在可接受的误差范围内恒成立，即DSMC可确保平动Eucken因子ftr非 常接近 2.5。但对于分子气体，DSMC不一定能够准确恢复热导率。在DSMC的碰撞模型中，可变软球模型的提出是为了获得准确的剪切黏性与自扩散系数；唯象的Larsen-Borgnakke模型[55]可以通过恢复能量交换速率[56-58]从而准确获取体积黏性。然而迄今为止，尚无研究证明DSMC能够在分子气体模拟中正确恢复热导率。
我们通过提取DSMC中的热弛豫速率A来 分析其恢复热导率的能力。具体而言，在Larsen-Borgnakke模型中选取不同的转动碰撞数，对非极性气体N2与极性气体HCl进行研究。参照Bird书中的公式（4.50）和公式（5.67），DSMC中非弹性碰撞概率 Λ与本文中定义的转动碰撞数Z通过Jeans-Landau方程（64）联系在一起[24]：
式中，ω为黏性温度幂指数，见式（45）；α是可变软球模型中与斯密特数相关的参数：
使用开源软件SPARTA在空间均匀系统中提取DSMC方法中的热弛豫系数矩阵A。将一百万个DSMC模拟粒子均匀地分布在边长为10 nm的立方体内，所有的边界设置为周期性边界条件。气体数密度为n0=2.69×1025m−3， 温度为T0=300 K。如图3（a）所示，在模拟初始时刻，沿x轴正方向运动的粒子的速度分布函数为200 K温度下的麦克斯韦分布，沿x轴负方向运动的粒子速度分布函数为 400 K温度下的麦克斯韦分布。类似地，沿x轴正方向运动的粒子的转动能满足200 K温度下的麦克斯韦分布，而沿反方向运动的粒子转动能满足400 K温度下的麦克斯韦分布，见图3（b）。这样的设置是为了建立初始平动和转动热流分布。随后记录平动和转动热流随时间的演化，直至其在系统中按照熵增原理达到热平衡后消失。图3（c）展示的是通过100次系综平均得到的平滑计算结果。
图3 从DSMC模拟中提取热流的弛豫速率：(a)初始状态用于激发平动热流的速度分布函数，横坐标由最概然速率v m归一化； (b)初始状态用于转动热流的转动能分布函数，横坐标由k BT0归一化；(c, d)氮气平动热流和转动热流的弛豫过程，模拟参数为：式(78)中的非弹性碰撞概率Λ =0.25 ， 施密特数S c=1/1.33， 转动自由度d r=2 ， 温度黏性指数ω =0.74. 图中数据来源于文献[60]Fig. 3 Extraction of heat flux relaxation rates from DSMC simulations. The initial distribution of (a) molecular velocity and (b) rotational energy of nitrogen molecules in DSMC, where the abscissas are normalized by and k BT0, respectively. (c, d) The evolution of heat fluxes and their time derivatives in Nitrogen. Λ =0.25 in Eq. (78), the Schmidt number is S c=1/1.33, the rotational degree of freedom is d r=2 , and the viscosity index is ω =0.74. The figures are modified from figure 1 in Ref. [60]
图3（d）给出N2的热流变化率的演化过程，可以看出平动热流的时间导数随时间先增加后减小，这可以解释为平动热流与转动热流之间存在强耦合效应，即在平动热流衰减的同时得到了转动热流的补偿，且根据式（73）可知，Atr必定小于零。
图4显示通过最小二乘法求解线性回归问题（式73）得到的热流弛豫速率A。可以看出，在相同的转动碰撞数Z和斯密特数下，N2与HCl气体的热弛豫速率几乎相同。另外，当Z→∞ 时 ，有Att=2/3，Arr=Sc，Atr=Art=0，可以得到与式（3）一致的结果。图中虚线为Mason和Manchick的结果（式（76）），由于忽略了1 /Z的高阶项，其结果在1 /Z较大时与DSMC的结果出现明显偏差。
图4 根据理论公式(76)从DSMC模拟中提取的热弛豫速率. 施密特数为S c=1/1.33， 转动自由度为d r=2; 对N2和HCl，温度黏性指数ω 分别取0.74和1，图中数据来源于文献[60]Fig. 4 The extracted rates A in the relaxation of heat fluxes from the DSMC simulation, as per Eq. (76). The Schmidt number is Sc=1/1.33, and the rotational degree of freedom is d r=2; the viscosity index for Nitrogen and hydrogen chloride is ω = 0.74 and 1,respectively. The figures are modified from figure 1 in Ref. [60]
图5给出根据式（74）和式（75）计算得到的Eucken式因子，并与常温下的实验结果对比。对于N2，当施密特数取1 /1.33时 ，通过选择特定的碰撞数1 /Λ≈3.5，DSMC可以恢复实验测量得到的总热导率，但由于缺少相应实验数据，我们无法分辨其是否能正确恢复平动与转动热导率分量。若此温度下N2的体积黏性要求 1 /Λ≉3.5，则DSMC给出的总热导率必然存在误差。对于极性气体HCl，实验中常温下测得的总Eucken因子约为1 .7，但是从图5（a）中可见，如果取施密特数为1 /1.33，从DSMC中导出的总Eucken因子最小只能达到1 .8。因此，无论如何选择碰撞数，DSMC都不能恢复总热导率。
图5 室温下Eucken因子的解析解(77)与DSMC数值解的对比. DSMC中使用可变软球模型，自扩散系数 D取值为 D'，菱形、正方形、三角形分别表示平动、内能和总Eucken因子，实心点为解析解，空心点为DSMC结果，虚线表示300 K温度下实验测量的总Eucken因子，图中数据来源于文献[54]Fig. 5 The thermal conductivity obtained from the analytical solution (77) and the DSMC at room temperature, as a function of the inelastic collision probability Λ. The variable-soft-sphere model is used and the self-diffusion coefficient D in DSMC takes the value of D' . Filled diamonds, squares and triangles represent the translational, internal and total Eucken factors from Eq. (77),respectively, while the open symbols are the corresponding results from DSMC. Dashed lines show the total Eucken factor obtained from experiments at a temperature of 300 K.The figures are from Ref. [54]
从式（77）可得，当气体种类确定时，转动碰撞数决定了体积黏性，因此可以通过调整DSMC中的扩散系数以恢复总热导率，见图5（b）。然而，即便在体积黏性与总热导率均与实验值相同的情况下，方程（77）的近似解与使用Larsen-Borgnakke模型的DSMC方法得到的平动热导率ftr并不一致，且两者的准确性都尚不确定。对于HCl气体，极性分子的共振相互作用[51]导致 ρD′/µ从 1 .33降 低至 0.89，为了恢复实验测量的总热导率，从图5（b）可知，DSMC的非弹性碰撞概率约等于1，达到了平动能与内能交换的极限，因此并不是一个满足物理要求的值。若选择合理的非弹性碰撞概率，则将要求DSMC中的有效扩散系数远小于自扩散系数。这将带来另一个问题：对于多组分气体，扩散系数与热导率同等重要，而Larsen-Borgnakke模型难以同时正确恢复这两个输运系数，亟需改进。
3 分子气体的动理论模型及精度
WCU方程的每一个内能能级都对应一个速度分布函数，因而计算代价巨大。对于强激波问题，N2的能级数可达到70，给数值求解带来了极大的困难[71]。因此对WCU方程进行简化是非常必要的，而简化的基本原则与单原子气体相同，即需要满足必要的守恒律和平衡态分布，还需要尽量使各物理量的输运系数与WCU方程保持一致，如扩散系数、剪切黏性、体积黏性、热导率及其在各自由度下的分量等。更接近本质的说法是，各输运系数对应的物理量的弛豫速率与实际气体保持一致。学者们提出了许多描述分子气体稀薄效应的动理学模型方程，如Rykov模型和ES-BGK模型等。与在单原子气体中Gross-Jackson模型可以视为所有模型的源头类似，在分子气体中可以认为Hanson-Morse模型[78]是所有动理论模型的源头。
3.1 Hanson-Morse线性化模型
此模型在线性化情况下逼近WCU方程（61）的碰撞项。这种构造方法类似于单原子气体的Gross-Jackson模型[79]，即通过对玻尔兹曼碰撞项进行正交多项式展开，并保证各多项式对应的宏观量的弛豫速率与WCU方程一致而获得。对于分子气体，关于分子速度的多项式仍由麦克斯韦分子的本征函数（24）给出；关于内部能量的多项式的前两项取为1和ε=ei−e¯， 其中e¯为平均内部能量[52]。
分子气体的平衡态速度分布函数可表示为fi=oiFeq，其中
表征了内能为ei的分子比例。在线性化问题中，通常考虑宏观速度为零的全局平衡态，并将速度分布函数写为如下形式：
其中每个内部能级对应的扰动分布函数 φi满足|φi|≪1。注意到为与上文的单原子线性化模型保持一致，这里的速度和空间坐标分别通过最概然速率vm和 特征尺寸 ℓ 进行无量纲化，时间通过 ℓ/vm归一化，加速度通过v2m/ℓ无 量纲化，内能ei通 过kBT0无量纲化。归一化后的扰动物理量可以通过对扰动分布函数求矩得到：
若加速度为小量，WCU方程的线性化形式为：
其中Ji的形式与式（21）类似，不再赘述。根据Hanson-Morse模型[78]，该线性化碰撞项可近似为[78]：
Hanson-Morse模型假设各内部能级的弛豫时间相同，则与式（76）对应的弛豫率为：
在自发瑞利-布里渊散射中，如果模型中包含应力偏量 σαβ，则该模型称为Tenti-S7模型[80-81]，其缺点在于不能得到实验散射光谱。Tenti-S6模型通过移除应力偏量 σαβ，即，令模型中的 σαβ=0，获得与实验一致的散射光谱。实际上，如果不考虑内部能级，Tenti-S7和Tenti-S6模型分别退化为单原子ES-BGK模型和Shakhov模型，而在大多数情况下Shakhov模型的精度高于ES-BGK模型。
当在Hanson-Morse模型中移除应力偏量，为了复现应力偏量、能量弛豫、热流的弛豫过程，式（84）中的J030应 该改成 − 1/τ 。同时，将系数J011修改为以下形式[80,82]，从而保证正确恢复热导率。
3.2 Hanson-Morse模型的约化及非线性化
为了减小计算量，引入如下两个约化速度分布函数：
消去内能自由度，则Hanson-Morse模型可以简化为：
其中：
至此，注意到单原子气体中Gross-Jackson模型与BGK、ES-BGK、Shakhov等模型的关系，可以通过类似于式（60）中展示的非线性化程序构造出分子气体的动理学模型。即，首先写出扰动量g和r对应的速度分布函数，
并定义宏观量如下：
总热流定义为q=qt+qr，平动压力、转动压力和总压力分别定义为pt=nkTt、pr=nkTr和p=nkT。 注意此处G可以线性化为G=feq(1+g)， 而R应线性化为R=(drkBT/2)feq(1+g+r)。其次，把速度分布函数的演化过程形式地写为：
最后，根据式（88）和式（89）中的碰撞项，类比式（60），写出四个参考分布函数Gt、Gr、Rt和Rr的形式。
例如，考虑Hanson-Morse模型中不出现应力偏量。式（88）中的第一行和第二行即为线性化Shakhov模型的碰撞项，可以非线性化为：
结合非线性模型（式92）中的非弹性碰撞项(Gr−Gt)/Zτ，可以选取
以在线性化情况下退化为式（88）中的后两项，其中能量交换由Feq(T)和Feq(Tt)的差距而来，热流松弛由q′决定。
3.3 非线性模型
然而，历史并没有选择非线性化路线，而是由许多学者各自提出了非线性动理学模型[33-34,50,61,83-87]。下面先介绍主要的Rykov模型和ES-BGK模型，并揭示与Hanson-Morse模型的内在联系，最后介绍吴模型。
3.3.1 Rykov模型
原始的Rykov模型[61]仅适用于无振动模态的双原子气体，但是构造思想可以直接拓展至多原子分子气体建模[83]。在该模型中，四个参考速度分布函数定义为：
选取弹性碰撞弛豫时间 τ=µ/pt以及适当的转动碰撞数Z，可保证恢复正确的剪切黏性和体积黏性。四个热流弛豫系数为：
因此通过式（75）可以得到Eucken因子的表达式如下：
这说明，在Rykov模型中，可通过调整自由参数 ω0和ω1修改Eucken系数。
当Z→∞，该模型退化为单原子气体的Shakhov模型。如果将其线性化，可以得到类似于第3.2节中的约化Hanson-Morse模型。不同之处在于空间均匀系统中，Rykov模型中平动和转动热流的弛豫不会相互干扰，这是因为式（74）中Atr=Art=0。
3.3.2 ES-BGK模型
Holway通过最大熵原理提出了ES-BGK模型，该模型保留了标准BGK模型的数学简便性，并且能够恢复正确的普朗特数，同时满足熵增定理[33-34]。速度分布函数的演化方程为：
其中，弛豫温度Trel定义为转动温度与总温度的加权平均值：
当Z→∞，转动和平动能量交换消失，因此模型方程退化为单原子气体的ES-BGK方程。
取弹性碰撞弛豫时间 τ=µ/ptPr及适当的转动碰撞数Z，则分子气体的剪切和体积黏性自动满足要求。其中普朗特数可以通过调整b的取值得到：
为使参考分布函数有界且在物理上有意义，要求Z≥1和−1/2≤b≤1。 因此，普朗特数范围为2 /3≤Pr≤+∞。对于大部分双原子分子来说，普朗特数大约为 5 /7，因此可以保证要求。
与Rykov模型一样，ES-BGK模型的平动热流和转动热流的弛豫过程互不干涉；从其弛豫率可知，Eucken因子可以表达为：
注意到，对于多原子ES-BGK模型，平动Eucken因子与转动Eucken因子的比值始终为 5/3，并且Eucken因子完全由普朗特数决定。
3.3.3 吴模型
上述气体动理学模型的弛豫时间 τ都与分子速度无关，与玻尔兹曼方程的弛豫时间不相符。因此，在正激波问题中，Rykov模型计算得到的结果可能会出现波前温升的误差，且随着马赫数的增加，这一偏离显著增长。由此，为了提高模型方程准确性，吴雷等首先将Rykov模型中的弹性碰撞项 (Gt−G)/τ替换为单原子玻尔兹曼方程的碰撞项Q(G)，以使弹性碰撞弛豫时间是分子速度的函数[83]；然后，又考虑引入更为准确的热流弛豫过程[60]，从而在相同参数条件下，得到更为符合DSMC的模拟结果。
吴模型的演化方程为：
式中，R′t=(dr/2)kBTr(τQ+G)， 弛豫时间为 τ=µ/pt，对应的参考速度分布函数为：
为了恢复热流弛豫过程（式（73）），选取q′和q′′如下：
当Z→∞，吴模型完全退化为单原子玻尔兹曼方程。另一方面，如果将玻尔兹曼碰撞项换回Shakhov碰撞模型，在线性化的情况下，吴模型则完全与第3.2节中的约化Hanson-Morse模型一致。
3.4 模型方程精度比较
本节以正激波结构、库特流、热蠕动、二维热驱动问题为例，比较不同模型描述稀薄气体流动的精度。模拟中，分子气体为N2，转动自由度dr=2，模型中的相关参数选为Z=2.6671、ω=0.74，施密特数Sc=1/1.33。 吴模型中的热弛豫系数矩阵A来自第2.3节中的DSMC模拟，从而保证吴模型的Eucken因子为feu=1.993、ftr=2.365和frot=1.436，与DSMC保持一致；此外，进一步假设这些参数不随温度变化。Rykov模型可以通过式（97）调整 ω0和 ω1恢复正确的Eucken系数。从式（102）可知，ES-BGK模型只有一个自由度来恢复总Eucken系数，即ES-BGK模型没有分辨平动Eucken系数和转动Eucken系数的机制。幸运的是，对于N2，ES-BGK模型给出的平动和转动热导率与DSMC非常接近，见表1中关于Eucken参数的汇总。再次强调，ES-BGK和Rykov模型中平动和转动热流的弛豫过程是相互独立的。
表1 不同动理学模型和DSMC模拟中所使用的Eucken因子与热弛豫速率. N2的Eucken因子为feu=1.993Table 1 Eucken factors and heat flux relaxation rates used in different gas kinetic models and DSMC. The total Eucken factor of the Nitrogen at room temperature is feu=1.993
引入参考密度n0、 参考温度T0、最概然速率参考压力n0kBT0、 参考热流n0kBT0vm。宏观量通过参考值进行无量纲化，且无量纲克努森数定义为：
用于表征流动的稀薄程度，其中 ℓ为流动特征长度。
在以下各研究问题中，气体在固体表面的反射均用完全漫反射作为边界条件，即气体反射速度分布为满足固体表面宏观量的麦克斯韦分布。在数值解法上，使用离散速度法求解动理学模型，快速谱方法求解玻尔兹曼碰撞项[31]，开源软件SPARTA进行DSMC模拟。
3.4.1 正激波
图6给出马赫数为7的一维N2正激波的结构。空间坐标用激波上游的分子平均自由程对其进行无量纲化，密度与温度通过 (Q−Qu)/(Qd−Qu) 归 一化，其中Q=n、Tt、Tr，下标u和d分别代表上游和下游的物理量。所有动理学模型的密度分布均与DSMC计算结果大体上相符，只有ES-BGK模型给出的密度在波前稍低一些。Shakhov模型和ES-BGK模型的平动温度在波前存在“过早温升”现象，其中ES-BGK的平动与转动温度分布偏离最大，这是因为这两个模型方程的碰撞时间与分子速度无关。而吴模型和DSMC中的碰撞时间都是分子速度的单调递增函数，因而得到最为接近的结果。图6（c、d）给出了压力与热流分布的计算结果，可见吴模型与DSMC符合得更好，而其他模型均存在显著误差，尤其是热流分布的对比。
图6 不同动理学模型预测的马赫数为7的氮气正激波的密度、温度、应力偏量σ22和热流的对比Fig. 6 Comparisons in the density, temperature, stress, and heat flux between different kinetic models for a normal shock in Nitrogen gas with Mach number 7
3.4.2 库特流动
平板库特流动的上下两板温度相同，处在x2=L/2的 上板速度是vw=vm， 处在x2=−L/2的下板速度相反。图7给出各模型和DSMC的对比结果。当Kn=0.1时，所有动理学模型都给出了良好的密度分布。ES-BGK模型高估了平动温度和平动热流，却低估了转动热流。Rykov模型给出良好的温度和转动热流，却高估了平动热流。吴模型与DSMC的对比结果非常理想。当Kn=1时，吴模型同样能给出与DSMC较为一致的结果。Rykov模型方程的密度分布与DSMC的结果稍有差别，在两板中心区域略低。ES-BGK模型的转动温度也略低于其他模型。
图7 平板库特流动中的密度、温度和垂直于流动方向的热流分布. 第一行和第二行的克努森数分别为K n=0.1、1； 基于对称性，另一半区域− 0.5≤x1≤0没有显示Fig. 7 Profiles of the density, temperature, and heat flux (perpendicular to the plates) in the planar Couette flow of the Knudsen number Kn = 0.1 (first row) and 1 (second row). Due to symmetry, the other half region − 0.5≤x1≤0 is not shown
3.4.3 麦克斯韦妖驱动的微流动
考察处在两平行静止的无限长平板间的分子气体，每个气体分子受到与其速度相关的水平外力的作用：
在此例中，取 2a0L/v2m=0.0718， 以保证系统仅略微偏离平衡态。可以看出，当分子速率大于时，其加速度为正，反之则为负。因此该加速度可以看出是一种特殊形式的麦克斯韦妖，见图8。图8给出克努森数分别为0.1和1的模拟结果对比。吴模型与DSMC非常符合，速度剖面与热流分布的最大误差仅为 2%。值得注意的是，ES-BGK模型和Rykov模型的转动热流值远大于吴模型和DSMC的结果。其原因在于，这两个模型方程中的平动热流弛豫与转动热流弛豫相解耦，即Atr=Art=0，即便热导率相近，但热弛豫系数的差别仍可导致显著的宏观量的差异。这正说明了采用正确的热弛豫速率的重要性。
图8 麦克斯韦妖驱动下的稀薄气体热蠕动, 其中速度和热流均沿 x 1方向. 第二行和第三行图对应的克努森数分别为Kn = 0.1和1，基于对称性，0 .5≤x2≤1的区域没有显示Fig. 8 Profiles of the vertical velocity and heat flux in the thermal creep flow driven by the Maxwell demon, where the Knudsen number in the second and third rows is Kn = 0.1 and 1, respectively. Due to symmetry, the other half region 0 .5≤x2≤1 is not shown
3.4.4 封闭方腔中的热蠕动
最后，考察二维封闭方腔内的热蠕动问题。方腔宽度与高度比值为5，固定左板和右板的温度分别为200 K和400 K，而上下两板的温度在此间线性分布。由于此问题以x2=0.5为轴具有对称性，实际计算区域为方腔的下半部分。取参考温度为300 K，参考长度为方腔高度，模拟克努森数为Kn=0.6的气体流动。
图9中给出各个动理学模型方程与DSMC结果的比较。吴模型的热流云图与速度剖面均与DSMC结果相符，而ES-BGK模型热流与速度剖面的解与DSMC的解有非常大的偏离，并且ES-BGK模型计算得到的中心涡更偏向右侧的高温区，且结构范围更小。另外，Rykov模型在预测热流方面存在一定误差。
4 热弛豫速率不确定度量化
为了准确地预测稀薄分子气体的非平衡行为，气体动理学的建模过程应实现正确的剪切黏性、体积黏性、平动热导率、转动热导率。分子气体的热导率比剪切黏性和体积黏性更加复杂，实验通常只能确定总热导率，但是难以分辨其平动与转动分量。虽然在连续流极限下，气体动力学行为取决于气体黏性与总热导率，但是在稀薄流区，分子气体的平动与内能热导率可能对非平衡效应有各自不同的影响。
上文的算例表明，热弛豫速率在稀薄气体流动中也起着重要作用。尽管在DSMC和其他分子气体动理学模型中[33,72,88-89]，已有大量关于温度弛豫效应（与体积黏性相关）的研究工作，却从未有过关于热弛豫速率矩阵A对稀薄气体非平衡效应的影响的相关研究。通过实验可以直接测量总热导率，并且在某些情况下，还可以提取平动热导率[54,59,90-91]。然而，根据式（73），弛豫系数矩阵A中至少有两个元素不能确定，即不同的系数矩阵A可以对应一致的热导率及其分量。因此有必要量化在输运系数都一致的前提下，由系数A的变化所带来的宏观量的不确定性。
在DSMC中，因碰撞模型的原因，固定体积黏性后，热弛豫系数也随之确定。在吴模型中，可以通过独立修改Aij以量化热弛豫系数的影响。上文已验证吴模型的准确性，这里使用吴模型来量化这种不确定性。首先，在约束ftr和frot不变的前提下修改Aij，以研究热弛豫速率的影响。其次，在保持总热导率不变前提下，通过改变ftr和frot来量化热导率分量的影响。
4.1 正激波
给定ftr和frot，即体系中所有的热导率分量都被确定，若此时其余的输运系数也保持不变，则称此系统为宏观意义上的确定系统。然而，不同的Aij的选择将可能导致结果的不确定性。
考察一维正激波问题，Atr、Art的取值范围分别为[−5/(6Z),0]和 [−1/(3Z),0]， 而Att和Arr的值根据式（73）调整，以保持ftr和frot不变。在DSMC中非对角线元素 为Atr=−0.201,Art=−0.059，本 例 中，给 定Z=2.6671，Art和Atr的 最小值分别为− 0.3124 和 − 0.1250，数值大小约为DSMC对应参数的一到两倍。研究发现，热弛豫系数A的变化导致转动温度和转动热流稍有变化，而对密度、速度和应力偏量几乎没有影响。
接着，保持总热导率不变，且选择Atr=−5/(6Z),Art=−1/(3Z)，研究热导率分量的改变对流动的影响。图10展示了马赫数为4时不同平动Eucken因子下吴模型的数值结果。在实际中，较小数值的平动Eucken因子ftr是可能存在的，特别是在极性气体中，真实的平动Eucken因子可能远小于 2.5。例如，水（H2O）的平动Eucken因子为ftr=1.78，甲醇分子（CH3OH）[53]的平动Eucken因子为ftr=0.41。对比可见，所有的宏观量都随ftr改变而有明显的变化。首先，较大的ftr使得平动温度更快上升到其最大值，并在波后更快趋于下游温度，对于应力偏量和总热流也有相同的趋势。其次，ftr变化所导致的转动温度变化则集中于激波结构的中心处，较小的ftr导致更高的转动温度。最后，较大的ftr使得密度升高更快。
图10 总Eucken因子不变时，不同平动Eucken因子对马赫数为4的正激波结构的影响. 图中数据来源于文献[60]Fig. 10 Influence of the translational Eucken factor on the structure of normal shock wave with Mach number 4. The figures are from Ref. [60]
由此，在正激波问题中，由热弛豫系数A的变化引起的不确定性较小，而热导率分量的变化则会带来显著的不确定性。例如，ftr=1.5的转动热流大约是ftr=2.5的两倍，几乎与平动热流相当。
4.2 麦克斯韦妖驱动的微流动
使用与正激波算例中相同的系数矩阵A，使ftr和frot均保持不变，来考察麦克斯韦妖驱动的微流动中，热弛豫系数对流动速度和热流的影响。图11第一行图片数据为Kn=0.2的结果。可以看出，不同于正激波问题，此问题中不同热弛豫系数A对应的流动速度与热流有相当明显的差别。流动速度和平动热流的最大相对不确定性分别为1 6.7%和 1 7.6%。同时可以看出，流动的中间区域出现了较大的不确定性，而壁面附近的速度滑移与热流几乎不变。
为进一步研究平动Eucken因子的影响，考察Eucken因子ftr=1.5,2.0,2.5的算例，其他参数为Kn=0.2，fu=1.993，Atr=−5/6Z，Art=−1/3Z。如图11第二行图片数据所示，速度曲线与平动热流随着平动Eucken因子的变化而剧烈变化。ftr=2.5时的流动速度和热流比ftr=1.5时 增加了大约 68%。在第一行图中，在确定的ftr下，壁面速度滑移与热流不受平动Eucken因子的影响。而在第二行图中，速度与热流依赖于平动Eucken因子，例如壁面附近的速度与热流均与ftr呈正相关。这说明，在此问题中，平动Eucken因子ftr起主导作用。
图11 (第一行)热流弛豫速率对麦克斯韦妖驱动的热蠕动速度和热流分布的影响，以及吴模型中改变热弛豫速率的结果对比. 红色实线对应的 A数 值来自于DSMC，蓝色阴影部分来自吴模型计算结果，对应的取值范围为 A rt∈[−0.3124,0.0]、 A tr∈[−0.1250,0.0]，其他参数为Kn=0.2 , ftr=2.365 以 及 frot=1.435； (第二行)总Eucken因子 feu保 持不变，改变平动Eucken因子的结果对比，K n=0.2.图中数据来源于文献[60]Fig. 11 (First row) Influence of the thermal relaxation rates in the creep flow driven by the Maxwell demon, and comparison of the results of varying the thermal relaxation rate in the Wu model. Red solid lines are the results with A extracted from DSMC, blue shade region shows the results from the kinetic model (103), with A rt∈[−0.3124,0.0] and A tr∈[−0.1250,0.0]. Other parameters are K n=0.2,ftr=2.365 and frot=1.435. (Second row) Influence of the translational Eucken factor, while feu is fixed. The Knudsen number is K n=0.2.The figures are from Ref. [60]
平动Eucken因子在此问题中的重要性可得到如下解释：从方程（73）可以看出，热弛豫系数Atr和Art与平动能与转动能之间的能量交换有关。当Atr与Art为零时，平动热流与转动热流的弛豫过程互相解耦。例如，更多的计算结果显示，在Art=0的情况下，转动热流恒为零。这是因为麦克斯韦妖驱动的热蠕动流动中，外力项仅能直接影响平动能，而转动能和热流则需通过能量交换才能发生改变，具体取决于Atr、Art和Z。由于矩阵的非对角元素相对于其他元素为小量，这导致qrot≈0 ，而平动效应（qtr或者ftr）占主导。
5 确定热流弛豫速率的方法
上文的研究结果揭示了热流弛豫速率不确定性对稀薄气体流动的影响，因此有必要通过获取并使用正确的相应系数，来减少或消除模型的不确定性。我们认为可以尝试以下两种方法。
5.1 分子动力学模拟
从物理上来说，只要知道了分子间的相互作用势，根据牛顿第二运动定律，就可以通过分子动力学方法模拟气体的行为。由此，输运系数能够通过Green-Kubo公式获得[92-93]：
其中，V和T分别是模拟系综的体积和温度，τ为时间，尖括号表示系综平均；Jp=mcici−pVI，Jq=Eici，I是单位矩阵。对于单原子气体，Ei=mcici/2得到平动热导率。对于分子气体，取Ei=εi为 第i个分子的内部能量，则式（110）给出内部热导率。
然而，对于如何获得热流弛豫矩阵系数尚不清楚。一个可行的方法是仿照第2.3节中介绍的DSMC方法，即在分子动力学中设置带有初始热流的状态，模拟低密度气体的演化过程，然后通过拟合公式（73）获得弛豫系数。虽然目前尚未有相关工作的报道，但是原理上是可行的。
5.2 瑞利-布里渊散射
瑞利-布里渊散射指的是光在气体中传播时，被气体分子的自发密度涨落散射。由于散射光谱包含气体信息，因此可以用来无损测量声速、流速、温度和体积黏性[94-96]。例如，欧洲航天局于2018年发射的ADM-Aeolus卫星，测量从地球表面到30 km高度的全球风速，用于提高天气预报精度。
如图12所示，波矢为ki的入射光被散射，散射角为θ，则散射波数为k=|ki−ks|=2|ki|sin(θ/2)。散射光强度与光谱密度函数S(fs)相关，即密度扰动相关函数G(r,t) 的 时空傅里叶变换。假设散射波沿x2方向传播，相关函数和谱密度函数的解[97-98]为：
图12 自发瑞利-布里渊散射示意图[54]，在气体中传播的光被气体分子自发的密度涨落散射Fig. 12 Schematic of the spontaneous Rayleigh-Brillouin scattering[54], where the light is scattered by the spontaneous density fluctuations in the gas
因此，频谱为：
式中，i为虚数单位，n(x2,t)为 密度扰动，fs为散射频率。
对于短波长和高散射频率，应该使用气体动理论来描述密度扰动的演化过程。以单原子气体为例说明瑞利-布里渊散射的计算。密度扰动量与式（20）中分布函数扰动量 φ有关，初始密度脉动可定义为[83,97]：φ(t=0,x2,ξ)∝δ(x2)。 取特征长度为散射波长 ℓ=2π/k，通过时间上的拉普拉斯变换与空间上的傅里叶变换
求解线性化玻尔兹曼方程（21），可得：
其中方程右端的1来自于初始密度扰动对应的拉普拉斯变换，(v) 分 别对应 φ的时空坐标下的拉普拉斯和傅里叶变换[17]。散射频谱可通过以下公式计算：
其中 ℜ 表示复数的实部。
此例中，我们根据荷兰自由大学的实验数据[91]来提取N2的体积黏性和平动Eucken因子。其中，激光波长为 366.8 nm，散射角为90°，因此等效波长为ℓ=259.4 nm。实验条件下N2的转动自由度与振动自由度分别为dr=2和0。剪切黏性和热导率根据Sutherland公式计算。通过调整吴模型中非弹性碰撞数和平动热导率，使得实验数据和数值模拟得到的频谱偏差最小。表2给出了不同温度下气体物性参数和提取的碰撞数与热导率分量。
表2 自发瑞利-布里渊散射实验中，N2相关参数汇总以及根据吴模型从实验数据中提取的体积黏性与平动/转动Eucken因子(剪切黏性与总热导率均使用国际单位制)Table 2 Experimental conditions in the spontaneous Rayleigh-Brillouin scattering, and the extracted bulk viscosity,translational/rotational Eucken factors from the Wu model. The SI units are adopted for the shear viscosity and thermal conductivity
实验数据与吴模型最佳预测的比较如图13所示。实验[99]测得N2在标准温度和压力下的转动弛豫时间为 τr=7.4×10−10s。在 296.7 K下，测得的N2体积黏性为 µb=1.22×10−5kg·m−1·s−1。根据式（71），计算得到N2的体积黏性为：
图13 基于实验频谱(圆圈)和吴模型预测结果(线条)，提取N2的体积黏性与平动Eucken因子[54]：第一行的实线为计算光谱，通过计算给定 ftr 和 Z 范围内的最小误差；第二行给出了实验光谱与理论光谱之间的误差.实验与计算中所用到的参数汇总在表2中，频率 fs由 v m/ℓ归一化Fig. 13 Extraction of the bulk viscosity and translational Eucken factor of Nitrogen from the experimental spectra (circles)[54]. Lines in the first row show the spectra obtained from the Wu model, while those in the second row show the corresponding residuals between the experimental and theoretical spectra. Experimental conditions and extracted parameters are summarized in Table 2. The frequency fs is normalized by vm/ℓ
这一结果与实验值吻合较好。另外，提取的平动Eucken因子也在合理的范围。
6 总结与展望
本文回顾了国内外学者在稀薄效应动理学建模中所取得的进展，包括单原子气体与分子气体。由于高温真实气体效应，当分子内自由度增加时，分子间复杂的相互作用通常伴随着更加新颖、复杂的物理现象，这对稀薄分子气体动理学建模提出了很高的要求。我们指出，应力偏量、能量交换、热流的弛豫过程是输运现象的本质，而剪切黏性、体积黏性、平动/内部热导率等输运系数则是表象，即，随着克努森数的变化，这些输运系数对应的本构关系会失效，而相应的弛豫过程和速率不会改变。为了准确预测稀薄分子气体的非平衡行为，气体动理学的建模过程应尽可能准确地捕捉应力偏量、能量交换、热流的弛豫过程和速率。
本文还介绍了DSMC中提取热弛豫速率的方法，并讨论了DSMC方法恢复热导率的能力。虽然经典的Larsen-Borgnakke模型可以通过调整非弹性碰撞概率改变体积黏性，但是很难同时调整热导率及其分量。而这些分量的比例以及本质的热流弛豫时间对稀薄气体流动的影响巨大。因此，非常有必要对DSMC碰撞模型进行修改和优化。
非常有意思的是，基于王承书等的工作，分子气体的线性化理论和方法已经非常成熟。而非线性模型虽然五花八门，但都与线性化模型存在着本质联系，即所有的非线性模型都可以通过对线性化Gross-Jackson模型和Hanson-Morse模型的非线性化而来。本文首次给出一个非线性化方案。
对于临近空间高超声速飞行器而言，其周围的空气动力学非常复杂。如何建立描述多组分气体、化学反应的稀薄气体动理学模型是个极大的挑战。本文提出的弛豫过程本质论、输运系数表象论、从线性化模型中反推非线性模型的方法，以及从分子动力学和相关实验中提取弛豫速率的设想，将为解决此类问题提供极大助力。
致谢：感谢爱丁堡大学苏微博士的有益讨论。
