0 引　言

经典流体力学与空气动力学建立在连续介质假设之上，即宏观微元中需包含足够多的微观分子。与之相反，稀薄气体动力学描述的则是在此假设不成立时气体的流动与热力学等性质。克努森数 ${Kn}$ 是稀薄气体研究中一个重要的无量纲参数，用来表征气体的稀薄程度，通常定义为气体分子平均自由程与流动特征长度的比值。因而，当分子平均自由程较大（即气体密度与压强较小）时，或在流动特征长度极小（即微纳米尺度下），气体流动均可呈现出显著的稀薄效应。

由此，不难发现，稀薄气体流动广泛存在于航空航天、工业技术和能源环境等重要的国防、民生领域，如返回舱再入大气层以及临近空间飞行器高空域飞行时的绕流^[1-5]、微机电系统中的气体流动^[6-9]以及页岩气开采中甲烷在复杂多孔介质中的流动^[10]等。相关科学技术的发展对我国国防、国民经济和社会发展意义重大，如下所述：

1）临近空间飞行器既可以避免绝大部分的地面攻击，同时也能够有效实施对地攻击和对航天器的打击，因此临近空间是进行空中军事活动的理想区域，发展潜力极大，我国和不少其他国家已经将其纳入信息化武器装备体系^[11]。在70 km以上的飞行空域，复杂构型的高超声速飞行器表面存在极为复杂的流动，因而需要建立新的流动模型和计算方法，来预测复杂的稀薄气体效应、高温真实气体效应以及解决激波/膨胀波与黏性作用的耦合干扰等问题。深入研究上述实际应用中的稀薄气体流动规律是优化飞行器设计的关键所在，具有重要意义。

2）机电系统微小集成化是提高能源利用效率的重要手段，因为微小化可大幅度减少传热传质阻力从而显著提高工质的相关输运能力。近年来，微流动系统在电子、化学、生物等领域已经取得了巨大的发展并得到了广泛的应用，但是微尺度下的空气动力学理论和模拟方法仍不完善，还存在许多亟待解决的问题^[12]：如在电容式膜片压力计中，当被测气体的温度与压力计的温度不同时，稀薄气体自动从低温区向高温区蠕动，从而导致压力测量偏差，需要考虑稀薄气体效应进行修正^[13]；工作在低压风洞的测量设备也需要考虑相应的稀薄气体效应修正以提高测量精度。

3）随着自然资源逐渐枯竭和环境破坏愈发严重，能源转型和能源利用效率提升成为各国科技发展的重要目标。相较于传统化石能源，以甲烷为主要成分的页岩气，因其具有弥补常规能源储量及产量、减少温室气体排放的巨大潜力，在世界能源供给中变得愈发重要。值得一提的是，我国页岩气储备位列全球第一，但复杂的地质条件要求对其输运机制有更清晰深入的理解，才能实现经济有效的开采；虽然页岩气埋藏在地表深处，气体压力大，但是页岩孔隙直径在纳米尺度，因而同样需要考虑稀薄气体效应^[14]。

图1列举了典型稀薄气体流动中纳维-斯托克斯方程的数值预测结果与实验测量存在明显差距的例子。由此可以看出，基于连续性假设得到的流体力学方程不再适用于此类问题，其根本原因在于，在稀薄效应较强的区域内，由宏观量的一阶导数构成的线性本构关系不再准确。另一方面，直接模拟所有气体分子运动与相互作用的分子动力学方法，提供了一个更为准确地模拟稀薄气体流动的途径。然而，因其庞大的计算量，该方法仅能用于极小时空尺度下的气体输运问题。因此，对于更为广泛存在的介观尺度下稀薄气体流动，如何建立准确的动理论模型方程并有效求解，则是目前最为关切与亟待解决的两个问题。

（1）研究现状。

一般来说，传统的纳维-斯托克斯方程只能描述连续流区（ $Kn < 0.001$ ）的气体流动。随着克努森数的增加，气体流动将分别经历滑移流区（ $ 0.001< {{Kn}}<0.1 $ ）、过渡流区（ $ 0.1<{{Kn}}<10 $ ）和自由分子流区（ $Kn > 10$ ）。对于处于滑移区的流动，气体的稀薄效应集中于边界处，因此采用适当的速度滑移和温度跳跃边界条件来求解传统宏观方程即可得到较为准确的结果。过渡流区内的稀薄气体效应由边界扩展到流动内部；自由分子流区中的气体相互碰撞极少，直接模拟蒙特卡罗方法（direct simulation Monte Carlo，DSMC）是一种常用的求解方法。然而，在许多实际问题中，气体的稀薄程度可能存在跨尺度的现象。例如，中国空气动力研究与发展中心的江定武等发现^[19]，返回舱载入过程中会经历从自由分子流到近地面连续流的多流动尺度的变化，同时，当地克努森数在飞行器迎风面和背风面有若干数量级的差别，因而同时存在混合多尺度的流动^[20]。因此，一个能够在全流域下统一描述稀薄气体流动的模型方程与求解方法十分必要。历史上，曾尝试通过对速度分布函数及宏观量的渐近展开以及不同阶次取矩等方法，得到了十余个宏观方程来描述稀薄气体行为^[17]。相较于纳维-斯托克斯方程，这些方程的数量和复杂性大幅增加，而稳定性大幅降低。例如，Burnett方程和矩方程都曾被提出用来描述稀薄气体流动^[21-23]，但由于方程本身的不稳定性和边界条件的难确定性，未能有效应用于全流域稀薄流模拟。相反，基于气体动理论方法的建模和数值方法则在全流域稀薄气体模拟中取得了一定的成功。

作为最为重要的气体动理论方程，玻尔兹曼方程描述了单原子气体的速度分布函数在时间与空间中的演化规律。玻尔兹曼方程的求解，在航空航天、微机电系统和页岩气开发等领域中都发挥着重要作用。其中，基于统计随机方法的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]，而实际气体大多是分子气体，即一个气体分子包含两个或两个以上的原子，例如占据地球大气绝大部分比例的N₂和O₂，占据火星大气绝大部分比例的CO₂等。相对于仅具有平动自由度的单原子气体，分子气体中还存在转动与振动的内能形式，并且各类型自由度的能量可以相互转化。对于常见气体，其转动能在室温下已充分激发，而振动能也在上千度高温下逐渐激发。随着温度的继续升高，分子气体的离解与电离也将逐渐发生，与此同时，化学反应也通常存在于此。此时，高温真实气体效应显著。而单原子玻尔兹曼方程无法包含这些过程，因此需要更复杂的模型方程来给出相应的模拟。目前，大部分以此为目标的建模工作是针对连续流提出的，例如对高超声速飞行器周围气体的模拟。由于分子气体的热流由平动能量和内部能量的输运与相互转换共同作用导致，因此需要通过额外的能量方程来描述高温气体能量输运，并且保证能量交换速率和总热导率与实验数值一样，此类连续流下的模拟方法就可以较为准确地模拟高温真实气体效应，国内外学者在这方面已经做了大量工作方面已经开展了大量工作^[51]。

然而，对于临近空间飞行器等类似问题而言，稀薄气体效应也非常显著^[19]，因此有必要将稀薄气体效应和高温真实气体效应结合起来考虑。对此，王承书和乌伦贝克建立的WCU方程^[52]给出了相对完备的描述，该方程以经典力学方式描述平动能量，而以量子力学方式描述内部（转动与振动）能量，并对每个能级下的速度分布函数建立方程。显然，该方法使得方程组自由度过于庞大，致使其求解变得相当困难从而难以应用于实际问题。因此建立相对简化的模型方程尤为必要。与此同时，由于各类自由度的特征碰撞时间不同，导致稀薄气体效应的程度也不同，因此动理论建模难度也大为增加。

对于分子气体的动理论建模，如何在模型方程中准确反应气体的各个输运系数是首要问题。相对于单原子气体，由于更多内部自由度的存在，分子气体中出现了更多的输运过程及相应参数。其一为体积黏性，其表征了气体在膨胀压缩的做功过程中能量与内部自由度转换的阻力，这一阻力来源于有限时间的能量弛豫过程；其二为转动、振动热导率，即内部自由度同时参与了热流的传递，与热流的弛豫速率紧密相关。实际上，在分子气体动理论的建模中，正确的体积黏性是容易实现的，只需保证正确的内部能量（转动、振动等）和平动能量的弛豫速率即可。但是，实现正确的平动、转动和振动热导率却非常困难。本文以仅有转动自由度的分子气体为例说明问题所在。历史上，Eucken首先将气体总热导率 $ \kappa $ 与剪切黏性 $ \mu $ 的关系表示为：

$ \begin{split} \frac{\kappa m}{\mu k_B}& = \frac{\kappa_{tr} m}{\mu k_B}+\frac{\kappa_{rot} m}{\mu k_B}\\ &\equiv\frac{3}{2}f_{tr}+\frac{d_r}{2}f_{rot}\\ &\equiv\frac{3+d_r}{2}f_{eu}\end{split}$

(1)

其中： $ m $ 是分子质量； $ k_B $ 是玻尔兹曼常数； $ d_r $ 是分子的转动自由度的数量，如N₂在常温下具有两个转动自由度，则 $ d_r $ = 2； $ \kappa_{tr} $ 和 $ \kappa_{rot} $ 分别是气体的平动和转动热导率； $ f_{tr} $ 、 $ f_{rot} $ 和 $ f_{eu} $ 分别是平动、转动和总Eucken因子。对于单原子气体， $ f_{tr} $ 非常接近2.5，但对于分子气体， $ f_{tr} $ 和 $ f_{rot} $ 的数值很难确定。最初Eucken从简单分子动理论出发，认为：

$ f_{tr} = \frac{5}{2}, \quad f_{rot} = 1 $

(2)

后来，人们意识到气体内部热导率与分子扩散有关，将内部Eucken因子改写为：

$ f_{tr} = \frac{5}{2}, \quad f_{rot} = \frac{\rho{D}}{\mu}\equiv\frac{1}{{Sc}} $

(3)

其中 $ \rho $ 为气体密度， $ 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能给出正确的总热导率，也还需要检查 $ f_{tr} $ 和 $ f_{rot} $ 的数值是否与实际气体一致。

在近连续流条件下，能量在平动自由度与内部自由度的分配趋于一致，因而平动温度与转动和振动温度的差别极小，此时保证总体热导率的正确性，即可得到相应正确的结果。因而，在稀薄效应较强的区域，保证各个热导率分量的正确性就显得尤为关键。例如，爱因斯坦等很早就发现，对于热蠕动等稀薄气体流动，温度梯度引起的气体流动只与 $ f_{tr} $ 有关^[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,{\boldsymbol{x}},{\boldsymbol{v}}) $ 表示，是时间 $ t $ 、空间坐标 $ {\boldsymbol{x}} = (x_1,x_2,x_3) $ 和分子速度 $ {\boldsymbol{v}} = (v_1,v_2,v_3) $ 的函数。定义 $f(t,{\boldsymbol{x}},{\boldsymbol{v}}){\rm{d}}{\boldsymbol{x}}{\rm{d}}{\boldsymbol{v}}$ 是体积为 ${\rm{d}}{\boldsymbol{x}}{\rm{d}}{\boldsymbol{v}}$ 的相空间上的分子数，则气体的数密度 $ n $ 、宏观速度 $ {\boldsymbol{u}} $ 、温度 $ T $ 、压力张量 $ p_{ij} $ 和热流 $ {\boldsymbol{q}} $ 可以分别通过对速度分布函数求矩得到。

$ \begin{split} &n(t,{\boldsymbol{x}}) = \int{}f(t,{\boldsymbol{x}},{\boldsymbol{v}}){\rm{d}}{\boldsymbol{v}} \\ &{\boldsymbol{u}}(t,{\boldsymbol{x}}) = \frac{1}{n(t,{\boldsymbol{x}})}\int{\boldsymbol{v}}f(t,{\boldsymbol{x}},{\boldsymbol{v}}){\rm{d}}{\boldsymbol{v}}\\ &T(t,{\boldsymbol{x}}) = \frac{m}{3k_Bn(t,{\boldsymbol{x}})}\int{}c^2f(t,{\boldsymbol{x}},{\boldsymbol{v}}){\rm{d}}{\boldsymbol{v}} \\ &p_{ij}(t,{\boldsymbol{x}}) = {}m\int{}c_ic_jf(t,{\boldsymbol{x}},{\boldsymbol{v}}){\rm{d}}{\boldsymbol{v}} \\ &{\boldsymbol{q}}(t,{\boldsymbol{x}}) = \frac{m}{2}\int{}c^2{\boldsymbol{c}}f(t,{\boldsymbol{x}},{\boldsymbol{v}}){\rm{d}}{\boldsymbol{v}} \end{split} $

(4)

其中， $ {\boldsymbol{c}} = {\boldsymbol{v}}-{\boldsymbol{u}} $ 是气体分子热运动速度，即气体分子速度与当地宏观速度的矢量差，而 $ c $ 是热运动速率。注意本文其他矢量也如此表示：黑体表示矢量，斜体表示矢量的模。定义应力偏量为：

$ \sigma_{ij} = p_{ij}-p\delta_{ij} $

(5)

其中， $ p = nk_BT $ ， $ \delta_{ij} $ 为克罗内克函数。

对于稀疏气体（分子间距远远大于分子直径），在外部施加的加速度 $ {\boldsymbol{a}} = (a_1,a_2,a_3) $ 的作用下，速度分布函数的演化由玻尔兹曼方程描述：

$ \frac{\partial f}{\partial t}+{\boldsymbol{v}}\cdot\frac{\partial f}{\partial {\boldsymbol{x}}}+{\boldsymbol{a}}\cdot\frac{\partial f}{\partial {\boldsymbol{v}}} = {Q(f)}$

(6)

方程（6）左边的三项分别表示速度分布函数在时间上的变化、在速度作用下在物理空间的变化以及在外力作用下在速度空间的变化，右边表示使得分布函数趋于平衡态的气体分子的碰撞过程。在玻尔兹曼方程中，两体碰撞的形式为：

$ Q(f) = \iint{}B(\theta,{v}_r) \Big [f({\boldsymbol{v}}'_*)f({\boldsymbol{v}}')-f({\boldsymbol{v}}_*)f({\boldsymbol{v}})\Big]{\rm{d}}{{\boldsymbol{\varOmega}}}{\rm{d}}{\boldsymbol{v}}_* $

(7)

其中， $ {\boldsymbol{v}} $ 和 $ {\boldsymbol{v}}_\ast $ 分别是碰撞前两个分子的速度，而 $ {\boldsymbol{v}}' $ 和 $ {\boldsymbol{v}}'_\ast $ 是碰撞后的速度。由于碰撞前后两分子的距离足够远，因此其相互作用可以忽略不计。根据动量和能量守恒定律，碰撞前后速度关系如下：

$ \begin{split} {\boldsymbol{v}}' =& \frac{{\boldsymbol{v}}+{\boldsymbol{v}}_\ast}{2}+\frac{|{\boldsymbol{v}}-{\boldsymbol{v}}_\ast|}{2}{\boldsymbol{\varOmega}}\\ {\boldsymbol{v}}'_\ast = &\frac{{\boldsymbol{v}}+{\boldsymbol{v}}_\ast}{2}-\frac{|{\boldsymbol{v}}-{\boldsymbol{v}}_\ast|}{2}{\boldsymbol{\varOmega}}\end{split} $

(8)

碰撞示意图见图2。其中，碰撞前两分子的相对速度为 $ {\boldsymbol{v}}_r = {\boldsymbol{v}}-{\boldsymbol{v}}_\ast $ ，碰撞后的相对速度为 $ {\boldsymbol{v}}'-{\boldsymbol{v}}'_\ast $ 。 ${\boldsymbol{\varOmega}}$ 为定义在单位球空间的矢量，其与碰撞后的相对速度同方向。于是相对速度的偏转角 $ \theta $ 与碰撞前相对速度满足如下关系：

$ \cos\theta ={\boldsymbol{ \varOmega}}\cdot\frac{{\boldsymbol{v}}_r}{{v}_r} $

(9)

最后，碰撞核 $ B(\theta,v_r) $ 是碰撞偏转角度和相对速度的函数，具体形式取决于分子间的作用力。

1.1 偏转角和微分散射截面

假设分子间通过中心力场作用，作用势 $ \phi(r) $ 已知，其中 $ r $ 为分子间距，则偏转角可以通过经典力学和量子力学两种方式求解。若气体温度不低（如氦气的温度高于100 K），两种方式得到的输运系数（如黏性和热导率）相同^[62]。这里仅介绍经典力学的计算结果，如图2所示，偏转角可以表示为：

$ \theta(b,{v}_r) = {\text {π}}-2\int_0^{W_1}\left[1-W^2-\frac{4\phi(r)}{m{v}_r^2}\right]^{-\frac{1}{2}}{\rm{d}}W $

(10)

其中， $ W = b/r $ 为瞄准距离 $ b $ 与分子间距 $ r $ 的比值，积分上限 $ W_1 $ 对应于最短分子间距，即式（10）中括号内表达式等于零的方程的正根。

在气体动理论中，经常考虑如下形式的逆幂律分子作用势：

$ \phi(r) = \frac{K}{\eta-1}r^{1-\eta} $

(11)

其中，K表征分子间相互作用的强度。从式（10）可知，偏转角只与变量 $ s $ 有关，即 $\theta = \varTheta(s)$ ：

$ s = \left[\frac{m(\eta-1)}{4K}\right]^{\frac{1}{\eta-1}}bv^{\frac{2}{\eta-1}}_r $

(12)

微分散射截面定义为：

$ \sigma_{_D} = \frac{b\;\big|{\rm{d}}b\big|}{\sin\theta\;\big|{\rm{d}}\theta\big|} = \left[\frac{m(\eta-1)}{4K}\right]^{\frac{2}{1-\eta}}v^{\frac{4}{1-\eta}}_r \frac{s{\rm{d}}s}{\sin\varTheta{{\rm{d}}\varTheta}} $

(13)

而碰撞核为：

$ B(\theta,v_r) = v_r \sigma_{_D} $

(14)

当式（11）中 $ \eta = 5 $ 时，即为麦克斯韦分子，此时碰撞核与相对碰撞速度无关，记为 $B(\theta,v_r) = \sqrt{2K/m}B(\theta)$ 。对于硬球分子模型，可看作式（11）中 $ \eta = \infty $ 。从图2可以看出，偏转角可通过如下公式确定： $ b = \sigma\cos\left({\theta}/{2}\right) $ ，其中 $ \sigma $ 为硬球直径。因此微分散射截面为 $\sigma_{_D} = {{\sigma}^2}/{4}$ ，而碰撞核为 $ B(\theta,{v}_r) = {\sigma}^2v_r/4 $ 。对于一般的气体，其行为介于麦克斯韦分子和硬球分子之间。

1.2 熵增原理

对于任意关于分子速度的函数 $\varPsi({\boldsymbol{v}})$ ，对玻尔兹曼碰撞项（式（7））在速度空间积分，具有如下对称性：

$ \begin{split} \int\varPsi({\boldsymbol{v}})Q({\boldsymbol{v}}){\rm{d}}{\boldsymbol{v}} = &\frac{1}{4}\iiint {\rm{d}}{\boldsymbol{\varOmega}} {\rm{d}}{\boldsymbol{v}}_\ast{{\rm{d}}{\boldsymbol{v}}} \Delta[\varPsi] B(\theta,v_r) \cdot \\ & \left[f({\boldsymbol{v}}'_{\ast})f({\boldsymbol{v}}')-f({\boldsymbol{v}}_{\ast})f({\boldsymbol{v}})\right] \end{split} $

(15)

其中 $\Delta[\varPsi] = \varPsi({\boldsymbol{v}}_{\ast})+\varPsi({\boldsymbol{v}}) -\varPsi({\boldsymbol{v}}'_{\ast})-\varPsi({\boldsymbol{v}}')$ 。若 $\varPsi$ 满足 $\displaystyle\int\varPsi({\boldsymbol{v}}) $ $ Q{\rm{d}}{\boldsymbol{v}} = 0$ ，则称其为碰撞不变量。根据质量、动量和能量守恒，可知 $\varPsi = 1、{\boldsymbol{v}}、v^2$ 为碰撞不变量，而碰撞不变量的线性组合也是碰撞不变量。

定义 $ H $ 函数为：

$ H = -\iint {f\ln}f{\rm{d}}{\boldsymbol{v}}{\rm{d}}{\boldsymbol{x}} $

(16)

H是与气体系统的熵相关的标量。在无外力的情况下，H的时间导数可写为：

$ \begin{split} \frac{\partial H}{\partial t} = & -\iint(1+{\ln}f)\frac{\partial{f}}{\partial{t}}{\rm{d}}{\boldsymbol{v}}{\rm{d}}{\boldsymbol{x}}\\ = &\iint{{\boldsymbol{v}}\cdot\frac{\partial{f}\ln{f}}{\partial{\boldsymbol{x}}}}{\rm{d}}{\boldsymbol{v}}{\rm{d}}{\boldsymbol{x}}-\iint(1+\ln{f}){Q} {\rm{d}}{\boldsymbol{v}}{\rm{d}}{\boldsymbol{x}}\\ = &\oint \int{{\boldsymbol{v}}\cdot{\boldsymbol{n}} {f}\ln{f}}{\rm{d}}{\boldsymbol{v}}{\rm{d}}S-\iint(1+\ln{f}){Q} {\rm{d}}{\boldsymbol{v}}{\rm{d}}{\boldsymbol{x}}\end{split} $

式中， $ {\boldsymbol{n}} $ 为系统表面微元 ${\rm{d}}S$ 的外法线方向。对于孤立系统，上式右端第一项为零；利用方程（15），右端项第二项可写为：

$ \begin{split} \frac{1}{4}\iiint{}&{\rm{d}}{\boldsymbol{\varOmega}} {\rm{d}}{\boldsymbol{v}}_\ast{{\rm{d}}{\boldsymbol{v}}} B \left[f({\boldsymbol{v}}'_{\ast})f({\boldsymbol{v}}')-f({\boldsymbol{v}}_{\ast})f({\boldsymbol{v}})\right]\\ &\times \Big\{ \ln[f({\boldsymbol{v}}'_{\ast})f({\boldsymbol{v}}')] -\ln[f({\boldsymbol{v}}_{\ast})f({\boldsymbol{v}})] \Big\} \end{split} $

因为碰撞核 $ B $ 非负，且对于任意的两个正整数 $ a $ 和 $ b $ ，不等式 $ (a-b)(\ln{a}-\ln{b})\geqslant 0 $ 恒成立，所以有：

$ \frac{\partial H}{\partial t}\geqslant 0 $

(17)

式（17）表明孤立系统的 $ H $ 函数不会减小，这就是著名的玻尔兹曼的熵增原理。 $ H $ 随时间单调增加，但存在有限的上界，上界对应式（17）中等号成立时的平衡态，即 $ \ln{f}({\boldsymbol{v}}_{\ast})+\ln{f}({\boldsymbol{v}}) = \ln{f}({\boldsymbol{v}}'_{\ast})+\ln{f}({\boldsymbol{v}}') $ 。因此， $ \ln{f} $ 也是碰撞不变量，可以表示为五个基本碰撞不变量的线性组合： $ \ln{f} = {\alpha}_1+{\boldsymbol{{\alpha}}}_2\cdot{\boldsymbol{v}}+{\alpha}_3v^2 $ ；给定系统的密度、速度和温度，参数 $ {\alpha}_1 $ 、 $ {\alpha}_3 $ 和 $ {\alpha}_2 $ 可以唯一确定。此时，麦克斯韦平衡态速度分布函数为：

$ F_{eq}(T) = n\left(\frac{m}{2{\text{π}} k_BT}\right)^{\frac{3}{2}}\exp\left(-\frac{mc^2}{2k_BT}\right) $

(18)

1.3 线性化玻尔兹曼方程

线性化玻尔兹曼方程在气体动理论中具有重要地位。第一，线性化玻尔兹曼碰撞项的本征值与本征函数，不仅在渐近展开推导纳维-斯托克斯方程的过程中至关重要，而且是发展简化动理学模型方程的源头理论。第二，在许多微机电系统中，气体的压力梯度与温度梯度非常小，因此使用线性化方程可以高效准确地模拟微流动。第三，虽然玻尔兹曼方程是单粒子速度分布函数确定性演化的平均方程，但是在某些问题中（例如：第5.2节中的瑞利-布里渊散射）可以通过线化玻尔兹曼方程研究粒子涨落带来的影响。

通常将速度分布函数在全局平衡态

$ f_{eq} = n_0\left(\frac{m}{2{\text{π}} k_BT_0}\right)^{\frac{3}{2}}\exp\left(-\frac{mv^2}{2k_BT_0}\right) $

(19)

下展开为：

$ f(t,{\boldsymbol{x}},{\boldsymbol{v}}) = f_{eq}({\boldsymbol{v}})\left[1+\phi(t,{\boldsymbol{x}},{\boldsymbol{v}})\right] $

(20)

其中扰动量 $ \phi $ 满足约束条件 $ |\phi|\ll1 $ 。不考虑外力项，只保留 $ \phi $ 的一次项，玻尔兹曼方程（6）可线性化为如下形式：

$ \begin{split} \frac{\partial\phi}{\partial t}+{\boldsymbol{\xi}}\cdot\frac{\partial\phi}{\partial {\boldsymbol{x}}} & = \frac{\ell}{v_m}{J}(\phi)\\ &\equiv -\frac{\ell}{v_m}\iint B \Delta[\phi]f_{eq}({\boldsymbol{v}}_{\ast}){\rm{d}}{\boldsymbol{\varOmega}} {\rm{d}}{\boldsymbol{v}}_\ast \end{split} $

(21)

这里的速度和空间坐标分别通过最概然速率 $ v_m $ 和特征尺寸 $ \ell $ 进行无量纲化，时间用 $ \ell/v_m $ 归一化；无量纲分子速度为：

$ {\boldsymbol{\xi}} = \frac{{\boldsymbol{c}}}{v_m(T_0)} $

(22)

其中：

$ v_m(T) = \sqrt{\frac{2k_BT}{m}} $

(23)

为气体分子的最概然速率，即与麦克斯韦速率分布的极大值对应的速率。

对于麦克斯韦分子，线性化玻尔兹曼碰撞项的本征值与本征函数为^[52]：

$\begin{split} J({{\varPhi} _{rlm}}) =& n{\lambda _{rl}}{{\varPhi} _{rlm}}\\ {{\varPhi} _{rlm}} =& {g_{rl}}({\xi} )Y_l^m(\hat {\xi} )\\ =& \sqrt {\frac{{2{{\text{π}} ^{\frac{3}{2}}}r!}}{{\left( {r + l + \dfrac{1}{2}} \right)!}}} S_{l + \frac{1}{2}}^{(r)}({{\xi} ^2}){{\xi} ^l}Y_l^m(\hat {\xi} )\\ {\lambda _{rl}} =& 2{\text{π}} \sqrt {\frac{{2K}}{m}} \int_0^{\text{π}} {\rm{d}} {\theta} \sin {\theta} B({\theta} )\times\\& \Bigg [ - 1 + {\cos ^{2r + l}}\frac{{\theta} }{2}{P_l}\left( {\cos \frac{{\theta} }{2}} \right)-\\& {\delta _{r0}}{\delta _{l0}} + {\sin ^{2r + l}}\frac{{\theta} }{2}{P_l}\left( {\sin \frac{{\theta} }{2}} \right)\Bigg] \end{split} $

(24)

式中， $ P_l(x) $ 是勒让德多项式， $ S^{(r)}_{l+\frac{1}{2}}(\xi) $ 是索南多项式， $ Y^m_{l}(\hat{\xi}) $ 为 $ {\boldsymbol{\xi}} $ 方向的球谐函数。

前三项本征值为 $ \lambda_{00} = \lambda_{01} = \lambda_{10} = 0 $ ，分别代表三大守恒律。另外两个非常重要的本征值为：

$ \begin{split} & {{\lambda} _{02}} = \left( { - \frac{3}{2}} \right){\text{π}} \sqrt {\frac{{2K}}{m}} \int_0^{\text{π}} {{\rm{d}}{\theta} } {\sin ^3}{\theta} B({\theta} )\\& {{\lambda} _{11}} = \frac{2}{3}{{\lambda} _{02}} \end{split} $

(25)

这两个本征值决定了剪切黏性系数和热导率的大小。若归一化的速度 $ {\boldsymbol{\xi}} $ 在 $ z $ 轴的投影为 $ \xi_z = \xi\cos\theta $ （其中 $ \theta $ 为极距角），与上述本征值对应的本征函数为：

$ \begin{split}& \varPhi_{000} = 1 ,\; \varPhi_{010} = \sqrt{2}\xi_z \\& \varPhi_{100} = \sqrt{\dfrac{2}{3}}\left(\dfrac{3}{2}-\xi^2\right) \\& \varPhi_{020} = \dfrac{1}{3\sqrt{3}}\left(\xi^2_z-\dfrac{\xi^2}{3}\right) \\& \varPhi_{110} = \dfrac{2}{\sqrt{5}}\left(\dfrac{5}{2}-\xi^2\right)\xi_z \end{split} $

其与气体的分子数密度、速度、温度、压力和热流的扰动密切相关：

$ \begin{split} &[{n},{\boldsymbol{u}}, T] = \int\left[1,{\boldsymbol{\xi}},\frac{2}{3}\left(\xi^2-\frac{3}{2}\right)\right] {f_{eq}\phi}{\rm{d}}{\boldsymbol{\xi}} \\ &[\sigma_{ij},{\boldsymbol{q}}] = \int \left[2\xi_{\langle i}\xi_{j\rangle}, \left(\xi^2 -\frac{5}{2}\right){\boldsymbol{\xi}}\right] f_{eq}\phi{}{\rm{d}}{\boldsymbol{\xi}}\end{split} $

(26)

式（26）中，下标中的尖括号表示该张量是无迹张量，即 $ \xi_{\langle i}\xi_{j\rangle} = \xi_i\xi_j-\xi^2\delta_{ij}/3 $ 。为避免使用过多符号，在线性化问题中，本文使用与原物理量相同的符号表示归一化的扰动物理量，如这里的扰动数密度 $ n $ 应理解为偏离参考数密度 $ n_0 $ 的那部分，再除以参考数密度；扰动温度 $ T $ 为偏离参考温度 $ T_0 $ 的部分，再除以参考温度；在全局平衡态下，参考速度、应力偏量和热流全部为零，因此上式中它们分别以 $ T_0 $ 下的最概然速率 $ v_m $ 、 $ n_0k_BT_0 $ 和 $ n_0k_BT_0v_m $ 归一化。

1.4 输运系数和弛豫率：表象与本质

玻尔兹曼方程和纳维-斯托克斯方程分别在介观和宏观尺度上描述气体系统的演化，它们之间的关系一直是数学家和力学家研究的重要课题^[63]。对玻尔兹曼方程求守恒量 $m{\text{、}} m{\boldsymbol{c}}{\text{、}}mc^2/2$ 的矩，可以得到如下宏观方程组：

$ \begin{array}{c} \dfrac{\partial \rho}{\partial t}+\dfrac{\partial(\rho{}u_j) }{\partial x_j} = 0 \\ \dfrac{\partial (\rho{u_i})}{\partial t}+\dfrac{\partial (\rho{}u_iu_j+p_{ij})}{\partial x_j} = \rho{}a_i\\ \dfrac{\partial \left(\rho{}E\right)}{\partial t}+\dfrac{\partial \left(\rho{}{E}u_j+u_i{p}_{ij}+q_j\right)}{\partial x_j} = \rho{}a_ju_j \end{array} $

(27)

其中 $ E = 3k_BT/2m+u^2/{2} $ 。但是该方程组并不封闭，因为压力张量和热流的表达式不能用密度、速度和温度表示。有三种主要方法尝试封闭宏观方程组，它们分别是希尔伯特展开法^[64]、Chapman-Enskog展开法^[21,65-66]和矩方法^[22,67]。这里介绍最常用的第二种方法。

在Chapman和Enskog^[65-66]方法中，首先将玻尔兹曼方程改写为：

$ \frac{\partial f}{\partial t}+{\boldsymbol{v}}\cdot\frac{\partial f}{\partial {\boldsymbol{x}}}+{\boldsymbol{a}}\cdot\frac{\partial f}{\partial {\boldsymbol{v}}} = \frac{Q(f)}{\varepsilon} $

(28)

然后将分布函数和玻尔兹曼碰撞项展开为 $ \varepsilon $ 的幂级数：

$ \begin{split}&\qquad\qquad\qquad f = \sum_{k = 0}^\infty {\varepsilon^k} f^{(k)}\\& Q = \sum_{k = 0}^\infty {\varepsilon^k} Q^{(k)} = Q(f^{(0)}) +\varepsilon{}{J}(f^{(1)})+{O}(\varepsilon^2)\end{split} $

(29)

值得注意的是， $ \varepsilon $ 是一个小的形式参数，主要用于监测展开的阶数。展开完成后，取 $ \varepsilon = 1 $ 。

同时，压力和热流也相应展开如下：

$ \begin{split} p_{ij} = \sum_{k = 0}^\infty {\varepsilon^k} p_{ij}^{(k)} \equiv \sum_{k = 0}^\infty {\varepsilon^k} \; m\int{}c_ic_jf^{(k)}{\rm{d}}{\boldsymbol{v}} \\ {\boldsymbol{q}} = \sum_{k = 0}^\infty {\varepsilon^k} {\boldsymbol{q}}^{(k)} \equiv \sum_{k = 0}^\infty {\varepsilon^k} \; \frac{m}{2}\int{}{\boldsymbol{c}}c^2f^{(k)}{\rm{d}}{\boldsymbol{v}} \end{split} $

(30)

但是，碰撞不变量对应的宏观量 $ C_M = \{\rho,{\boldsymbol{u}}, T\} $ 仅由分布函数的零阶展开决定，即：

$ \Big [{\rho},{\rho}{\boldsymbol{u}}, 3k_B{\rho}{}T \Big ] = m\int{} \Big [1,{\boldsymbol{v}},mc^2 \Big ]f^{(0)}{\rm{d}}{\boldsymbol{v}} $

(31)

而各个高阶展开对 $ C_M $ 的贡献均为零，这称为兼容性条件。

将式（30）代入式（27），可知时间导数亦可表示为 $ \varepsilon $ 级数^[67]，

$ \frac{\partial}{\partial{} t} = \sum_{k = 0}^\infty \varepsilon^k \frac{\partial}{\partial{} t_k} $

(32)

且对于零阶近似，可得到：

$ \begin{array}{c}\dfrac{\partial \rho}{\partial t_0}+\dfrac{\partial(\rho{}u_j) }{\partial x_j} = 0 \\ \dfrac{\partial (\rho{u_i})}{\partial t_0}+\dfrac{\partial (\rho{}u_iu_j+nk_BT\delta_{ij})}{\partial x_j} = \rho{}a_i\\ \dfrac{\partial \left(\rho{}E\right)}{\partial t_0}+\dfrac{\partial \left(\rho{}{E}u_j+u_ink_BT\delta_{ij}\right)}{\partial x_j} = \rho{}a_ju_j\end{array} $

(33)

从式（28）和式（29）易知 $ Q(f^{(0)}) = 0 $ 。因此，速度分布函数的零阶展开为：

$ f^{(0)} = F_{eq} $

(34)

而速度分布函数的一阶近似 $ f^{(1)} = F_{eq}\phi $ 满足以下方程^[67]：

$ \begin{split} {J}(\phi) =& \frac{1}{F_{eq}}\left[\frac{\partial F_{eq}}{\partial t_0}+{\boldsymbol{v}}\cdot\frac{\partial F_{eq}}{\partial {\boldsymbol{x}}}+{\boldsymbol{a}}\cdot\frac{\partial F_{eq}}{\partial {\boldsymbol{v}}} \right]\\ =& 2\xi_{\langle{i}} \xi_{j\rangle} \frac{\partial u_{\langle i}}{\partial x_{j\rangle}} +\left(\xi^2-\frac{5}{2}\right)\xi_i\sqrt{\frac{2k_BT}{m}}\frac{\partial \ln{T}}{\partial x_i}\end{split} $

(35)

其中，方程的推导过程中用到

$ \frac{\partial F_{eq}}{\partial t_0} = \frac{\partial F_{eq}}{\partial C_M} \frac{\partial C_M}{\partial t_0}, \quad \frac{\partial F_{eq}}{\partial {\boldsymbol{x}}} = \frac{\partial F_{eq}}{\partial C_M} \frac{\partial C_M}{\partial {\boldsymbol{x}}} $

以及式（33）。注意式（21）中的 $ f_{eq} $ 应改写为 $ F_{eq} $ 。

积分方程（35）的解可分为齐次部分与非齐次部分。齐次部分满足 $ {J}(\phi) = 0 $ ，这说明 $ \phi $ 必然是碰撞不变量的线性组合，而兼容性条件要求所有线性组合的系数都为零。非齐次部分满足如下形式：

$ \phi = -A(\xi)\xi_i\sqrt{\frac{2k_BT}{m}}\frac{\partial \ln{T}}{\partial x_i} -B(\xi)\xi_{\langle{i}} \xi_{j\rangle}\frac{\partial u_{\langle i}}{\partial x_{j\rangle}} $

(36)

其中， $ A(\xi) $ 和 $ B(\xi) $ 的解满足：

$ {J}(A\xi_i) = -\left(\xi^2-\frac{5}{2}\right)\xi_i, \;\; {J}(B\xi_{\langle{i}} \xi_{j\rangle}) = -2\xi_{\langle{i}} \xi_{j\rangle} $

(37)

且兼容性条件要求 $\int \xi^2A(\xi)F_{eq}{\rm{d}}{\boldsymbol{v}} = 0$ 。

一旦得到 $ A(\xi) $ 和 $ B(\xi) $ ，通过一阶展开就可以恢复牛顿黏性定理与傅里叶热传导定理：

$ \sigma^{(1)}_{ij} = -2\mu\frac{\partial u_{\langle i}}{\partial x_{j\rangle}},\quad q^{(1)}_i = -\kappa \frac{\partial T}{\partial x_i} $

(38)

且剪切黏性与热导率分别为：

$ \begin{split}& \mu = \frac{2p}{15{\text{π}}^{\frac3 2}}\int\exp(-\xi^2)B(\xi)\xi^4 {\rm{d}}{\boldsymbol{\xi}} \\& \kappa = \frac{2p}{3{\text{π}}^{\frac3 2}} \frac{k_B}{m} \int\exp(-\xi^2)A(\xi)\xi^4 {\rm{d}}{\boldsymbol{\xi}}\end{split} $

(39)

一般将 $ A(\xi) $ 和 $ B(\xi) $ 展开为索南多项式级数：

$ \begin{split} A(\xi)& = -\sum_{r = 1}^{n_a}a_rS^{(r)}_{\frac{3}{2}}(\xi^2) \\ B(\xi)& = \sum_{r = 0}^{n_b}b_rS^{(r)}_{\frac{5}{2}}(\xi^2)\end{split} $

(40)

然后通过式（37）和索南多项式的正交性求解输运系数。对于麦克斯韦分子，取第一项展开可得到准确的剪切黏性和热导率，即取 $ A(\xi) = a_1(\xi^2-5/2) $ 和 $ B(\xi) = b_0 $ ；而对于其他相互作用势，仅考虑第一项展开得到的输运系数的相对误差也不超过2%。定义等效黏度截面为：

$ \sigma_{\mu} = \int \frac{B(\theta,v_r)}{v_r}\sin^2\theta {{\rm{d}}\varOmega} $

(41)

则剪切黏性为：

$ \mu = \frac{5\sqrt{{\text{π}}{m}k_BT}}{8\left(\dfrac{m}{4k_BT}\right)^4\displaystyle\int_0^\infty {}v_r^7\sigma_{\mu}\exp\left(-\dfrac{mv^2_r}{4k_BT}\right){\rm{d}}v_r} $

(42)

而热导率为：

$ {\kappa} = \frac{15}{4}\frac{k_B}{m}\mu $

(43)

从而单原子气体的普朗特数为：

$ \begin{split} Pr & = {c_p}\frac{\mu}{\kappa} = \frac{5k_B}{2m}\frac{\mu}{\kappa} = \frac{2}{3}\end{split} $

(44)

对于逆幂律分子间相互作用势（式(11)），可以得出气体黏性和温度的关系：

$\mu(T)\propto{T^\omega}, \quad \omega=\frac{\eta+3}{2(\eta-1)}$

(45)

其中 $ \omega $ 是黏性温度幂指数。对于麦克斯韦分子和硬球分子， $ \omega $ 分别取1和0.5；对于其他气体， $ \omega $ 一般介于0.5和1之间。

将式（38）代回式（27），即可获得纳维-斯托克斯方程。若继续计算Chapman-Enskog展开的高阶近似，可以得到Burnett、super-Burnett等宏观方程组^[21]。但是可能由于应力和热流与宏观量 $ C_M $ 展开的不匹配，导致高阶方程组的线性不稳定性。同时，由于边界条件数量随着方程组阶数升高而增加，但对附加的边界条件提法缺少充分的研究，因此高阶方程组很少被实际应用。另一方面，即使这些高阶方程组在某些情况下能得到精确解，解的精度也不一定随着展开阶数的增加而增加^[17]。

值得一提的是，对于麦克斯韦分子，在空间均匀系统中，本征值 $ \lambda_{02} $ 和 $ \lambda_{11} $ 与应力偏量和热流的弛豫速率相关，而这些弛豫速率决定着剪切黏性系数和热传导的大小：

$ \begin{split}& \frac{\partial \sigma_{ij}}{\partial t} = n\lambda_{02} \sigma_{ij} = -\frac{p}{\mu}\sigma_{ij}\\& \frac{\partial q_{i}}{\partial t} = n\lambda_{11} q_{i} = -\frac{2}{3}\frac{p}{\mu}q_{i}\end{split} $

(46)

可以认为这些弛豫过程是本质，而黏性系数和热导率则为表象：在稀薄气体流动中，本构关系式（38）和等效输运系数会随着克努森数的改变而改变，但一旦分子作用势确定，弛豫系数就会确定不变。

利用弛豫时间和输运系数的关系，可以大大简化模型方程中输运系数的推导过程。即，如果应力偏量的弛豫时间为 $ \tau $ ，则剪切黏性为 $ p\tau $ 。如果热流的弛豫时间为 $ \tau/A $ ，则普朗特数为 $ A $ 。

1.5 模型方程

玻尔兹曼方程碰撞项的复杂性给数值求解带来极大挑战^[31]，因此，学者建议使用简化的碰撞模型来代替。在简化玻尔兹曼方程的碰撞项时，一般认为需要满足以下几点：

1）动理学模型必须保证质量、动量、能量守恒；

2）理想气体的局部平衡态速度分布函数为麦克斯韦分布；

3）从模型方程导出的黏性系数与热导率应与玻尔兹曼方程导出的结果保持一致；

4）满足熵增定理，当且仅当孤立系统处于平衡态时熵增为零，并达到其最大值。

其中，前两项是基本要求，第三项要求在连续流区域通过Chapman-Enskog展开恢复出纳维-斯托克斯方程。而第四项仅为从物理角度出发的考量，并不能保证模型方程的精度：如果一个动理学模型方程与玻尔兹曼方程完全一致，那么熵增速率也应该相同，然而这通常是不可能做到的。实际上，大多数时候，仅在线性化情况下满足熵增定理的Shakhov模型的精度优于完全满足熵增定理的ES-BGK模型。

由于玻尔兹曼碰撞项在形式上可记为 $ Q = Q^+ -f/\tau $ ，因此，动理学模型方程的碰撞项一般采取如下形式：

$ Q = \frac{f_{r}-f}{\tau} $

(47)

其中： $ f_{r} $ 为参考速度分布函数； $ \tau $ 为特征碰撞弛豫时间，表征速度分布函数趋向参考分布函数的快慢，为简单计，通常假设与分子速度无关。

1.5.1 BGK模型

BGK模型是最著名的气体动理学简化模型，由Bhatnagar、Gross和Krook提出^[32]，其参考速度分布函数是局部麦克斯韦分布：

$ f^{BGK}_r = F_{eq} $

(48)

易见此模型满足三大守恒律，且给出了在平衡态下的正确解 $ f = F_{eq} $ 。在空间均匀系统中，由于 $ \ln{F_{eq}} $ 是碰撞不变量的线性组合，有 $\int(F_{eq}-f)\ln{F_{eq}}{\rm{d}}{\boldsymbol{v}} = 0$ ，所以可以证明BGK模型满足熵增原理：

$ \begin{split}\frac{\partial H}{\partial t} = & -\int \frac{F_{eq}-f}{\tau}\ln{f} {\rm{d}}{\boldsymbol{v}}\\ = &-\int \frac{F_{eq}-f}{\tau}(\ln{}F_{eq}-\ln{f}) {\rm{d}}{\boldsymbol{v}}\\ \geqslant &\; 0 \end{split} $

(49)

空间均匀系统中应力偏量和热流的弛豫时间均为 $ \tau $ ，因此根据式（46）随后段落的描述，可知剪切黏性为 $ \mu = p\tau $ ，而普朗特数为1。然而，通过玻尔兹曼方程推导得到的普朗特数正确值为 $ 2/3 $ 。这说明BGK模型不可能同时得到正确的剪切黏性和热导率。对于黏性主导的流动，一般选取 $ \tau $ 以恢复剪切黏性；对于热传导主导的流动，则以恢复热导率为目标来选取 $ \tau $ 。

1.5.2 ES-BGK模型

Holway^[33]提出了ES-BGK模型以修正BGK模型中错误的普朗特数，其中参考速度分布函数是在给定质量、动量、能量和应力张量条件下，通过最大化熵函数式（16）而来：

$ f_r^{ES} = \frac{n}{\sqrt{{\rm{det}} [2\pi{\boldsymbol{\lambda}}_{ij}]}}\exp\left(-\frac{1}{2}{\boldsymbol{\lambda}}_{ij}^{-1}c_i{c_j}\right) $

(50)

其中：

$ {\boldsymbol{\lambda}}_{ij} = (1-b)\frac{k_BT}{m}\delta_{ij}+b\frac{p_{ij}}{nm} = \frac{p\delta_{ij}+b\sigma_{ij}}{nm} $

(51)

若 $ b = 0 $ ，二阶张量 ${\boldsymbol{\lambda}}_{ij}$ 中仅有主对角元素非零，模型恢复为BGK模型。

容易证明，应力偏量的弛豫时间为 $ \tau/(1-b) $ ，而热流的弛豫时间为 $ \tau $ 。因此，剪切黏性和普朗特数分别为：

$ \mu = \frac{1}{1-b}{p\tau}, \;\; Pr = \frac{1}{1-b} $

(52)

为确保弛豫时间为正，则有 $ b<1 $ 。另外， $ f_r^{ES} $ 有界的充要条件是 $b\geqslant -1/2$ 。因此，ES-BGK模型的普朗特数范围为 $ [2/3,+ \infty) $ 。为了恢复单原子气体的普朗特数 $ 2/3 $ ，取

$ b = -\frac{1}{2} $

(53)

可以看出，当碰撞项为零时，ES-BGK方程给出 $ f = f_{r}^{ES} $ ，当且仅当 $ \sigma_{ij} = 0 $ 时，速度分布函数退化为麦克斯韦平衡态分布。而这并不违背动理学建模的第二个要求，因为当局部平衡态为 $ f = f_{r}^{ES} $ 时，应力偏量的演化满足如下形式 $ \partial \sigma_{ij}/{\partial t} = -{p}\sigma_{ij}/\mu $ ，即应力偏量将随时间推进逐渐衰减至零。因此，ES-BGK模型中趋于平衡态的过程可以解释为：当速度分布函数朝着椭球分布（式（50））松弛，椭球分布也朝着麦克斯韦分布演化，最终实现 $ f = F_{eq} $ 。

Andries等^[34]证明了ES-BGK模型满足熵增定理，使之成为唯一满足熵增定理且能正确恢复各输运系数的动理学模型，因此，ES-BGK模型受到广泛关注。

1.5.3 Shakhov模型

不同于ES-BGK模型在参考速度分布函数中引入应力修正项，Shakhov模型通过在参考速度分布函数中引入热流修正项来恢复正确的输运系数：

$ f_r^S = F_{eq}\left[1+(1-{Pr})\frac{2m{\boldsymbol{q}}\cdot {\boldsymbol{c}}}{5n(k_BT)^2}\left(\frac{mc^2}{2k_BT}-\frac{5}{2}\right)\right] $

(54)

此模型同样满足三大守恒律。从偏应力和热流的弛豫时间可以看出，该模型的黏性为 $ \mu = p\tau $ ，且能实现正确的普朗特数。此外，由于热流弛豫过程满足 $ \partial {\boldsymbol{q}}/\partial t = -(2p/3\mu){\boldsymbol{q}} $ ，在速度分布函数向 $ f_r^S $ 靠近的过程中，参考速度分布函数也将逐渐演化为麦克斯韦分布，从而使Shakhov模型在平衡态时满足 $ f = F_{eq} $ 。

线性化的Shakhov模型满足熵增定理。而对于非线性流动，熵增定理未被证明但也无从证伪；另外，速度分布函数可能出现非物理的负值。尽管如此，Shakhov模型仍能在许多问题中给出准确结果，因此得到了广泛应用。

1.5.4 Gross-Jackson模型及其非线性化

针对麦克斯韦分子，根据线化玻尔兹曼碰撞项的本征值与本征函数，Gross和Jackson提出的以下形式的动理学模型去逼近线性化玻尔兹曼碰撞项 $ J(\phi) $ ：

$ \begin{split} {L}_{GJ} = &\lambda_{st}\phi+\sum_{nl}(\lambda_{nl}-\lambda_{st})\frac{2l+1}{4\pi} \cdot \\ &\int{}P_l(\cos\theta')g_{nl}(\xi)g_{nl}(\xi_1)\phi({\boldsymbol{\xi}}_1) {\rm{d}}{\boldsymbol{\xi}}_1 \end{split} $

(55)

式中， $ \lambda_{st} $ 是一个任意负数， $ \lambda_{nl} $ 是线性化玻尔兹曼项的本征值（式（24））， $ \theta' $ 为 $ {\boldsymbol{v}} $ 和 $ {\boldsymbol{v}}' $ 的夹角。该模型的物理意义在于保证模型方程和线性化玻尔兹曼方程各阶矩的弛豫率相同。

Gross和Jackson将 $ {L}(h) $ 的 $ N $ 阶近似表示为：

$ \begin{split} {L}_{GJ}^{(N)} = &\lambda_{0N}\phi -\lambda_{0N}\sum_{2n+l\leqslant{N}} \left(1-\frac{\lambda_{nl}}{\lambda_{0N}}\right)g_{nl}(\xi)\cdot \\ &\frac{2l+1}{4\pi} \int{}P_l(\cos\phi')g_{nl}(\xi_1)\phi({\boldsymbol{\xi}}_1) {\rm{d}}{\boldsymbol{\xi}}_1 \end{split} $

(56)

即保留了所有阶数小于或等于 $ N $ 的多项式，并且取 $ \lambda_{st} = \lambda_{0N} $ ，因此该截断模型的最大弛豫率为 $ \lambda_{0N} $ 。

在线性化情况下，即把速度分布函数在全局平衡态下根据式（20）展开，BGK、ES-BGK和Shakhov模型有如下形式：

$ \begin{split} &\frac{\partial \phi}{\partial{}t} +{{\boldsymbol{\xi}}}\cdot\frac{\partial \phi}{\partial {\boldsymbol{x}}} = \frac{\ell}{v_m}{L} \\ &{L}_{BGK} = \frac{p}{\mu}\left[{n}+2{\boldsymbol{u}}\cdot{\boldsymbol{\xi}}+T\left(\xi^2-\frac{3}{2}\right)-\phi\right] \\ &{L}_{ES} = \frac{2}{3}{L}_{BGK}-\frac{2p}{3\mu} \frac{\sigma_{ij}}{2}\xi_{\langle i}\xi_{j\rangle} \\ &{L}_{S} = {L}_{BGK} +\frac{p}{\mu}\frac{4(1-{Pr})}{5}{\boldsymbol{q}}\cdot{\boldsymbol{\xi}}\left(\xi^2-\frac{5}{2}\right) \end{split} $

(57)

其中扰动宏观物理量的定义可参考式（26）。它们都可以看作是Gross-Jackson模型的特例。当 $ N = 2 $ 时，从式（56）可以得到 $ {L}_{GJ}^{(2)} = {L}_{BGK} $ ；令 $ N = 3 $ ，根据式（24）有：

$ \lambda_{03} = \frac{3}{2}\lambda_{02} = \frac{3p}{2\mu} $

(58)

从而式（56）给出：

$ \begin{split} {L}_{GJ}^{(3)} = &\frac{3}{2}L_{BGK}+\frac{p}{\mu}\left[\frac{\sigma_{ij}}{2}\xi_{\langle i}\xi_{j\rangle} +\frac{2}{3}{\boldsymbol{q}}\cdot{\boldsymbol{v}}\left( \xi^2-\frac{5}{2}\right)\right]f_{eq} \\ = &-\frac{3}{2}{L}_{ES}+\frac{5}{2}{L}_{S}\end{split} $

虽然，这既不符合线性化ES-BGK模型，也不符合线性化Shakhov模型，但是可将其视为这两个模型的线性组合。另一方面，若给定约束 $ 2n+l\leqslant3 $ ，当 $ \lambda_{st} = -2p/3\mu $ 时，式（55）恰好展开为线性化ES-BGK模型；而当 $ \lambda_{st} = -p/\mu $ 时，式（55）展开为线性化Shakhov模型。

因此，通过Gross-Jackson模型和线性化BGK、ES-BGK、Shakhov模型的关系，可以反推出各种形式的非线性动理学模型。例如，陈松泽等将ES-BGK模型与Shakhov模型线性组合为一个复合模型，通过改变两模型所占的比例和普朗特数，复合模型能够分别恢复到BGK、ES-BGK和Shakhov模型式^[68]。类似地，吴雷在博士论文也提出一个带有自由参数 $ b $ 的复合模型^[69]：

$ \begin{split} f_r = &F_{eq} \left\{1+ \frac{\sigma_{ij}}{p} \frac{c_{\langle i}c_{j\rangle}}{2k_B\frac{T}{m}} +\right.\\ &\left.\Bigg[1-{Pr}(1-b)\Bigg]\frac{2m{\boldsymbol{q}}\cdot {\boldsymbol{c}}}{5n(k_BT)^2}\left(\frac{mc^2}{2k_BT}-\frac{5}{2}\right) \right\}\\ \tau = &(1-b)\frac{\mu}{p} \end{split} $

(59)

该模型实际上是在Gross-Jackson模型（55）中取 $ 2n+l\leqslant{N} = 3 $ 和 $ \lambda_{st} = 1/\tau $ ，并通过将如下的线性化碰撞项变为非线性碰撞项得来：

$ \begin{split} {n}+2{\boldsymbol{u}}\cdot{\boldsymbol{\xi}}+T\left(\xi^2-\frac{3}{2}\right)-\phi \rightarrow& {F_{eq}}-f\\ \sigma_{ij}\frac{\xi_{\langle i}\xi_{j\rangle}}{2} \rightarrow& \frac{\sigma_{ij}}{p}\frac{c_{\langle i}c_{j\rangle}}{2k_BT} F_{eq} \\ {\boldsymbol{q}}\cdot{\boldsymbol{\xi}}\left(\xi^2-\frac{5}{2}\right) \rightarrow& \frac{m{\boldsymbol{q}}\cdot {\boldsymbol{c}}}{2n(k_BT)^2}\left(\frac{mc^2}{2k_BT}-\frac{5}{2}\right)F_{eq} \end{split} $

(60)

需要注意到，含有应力偏量的项既可以被吸收到指数函数里使参考分布函数变成类似式（50）的形式，也可以放在外面（此时对 $ b $ 的限制仅为 $ b<1 $ ），两者对数值计算的结果影响不大^[69]。通过将Gross-Jackson模型中的其他阶本征函数非线性化为相应的形式，可以构造更加准确的动理学模型。

2 分子气体的特性

玻尔兹曼方程是仅针对单原子气体提出的，为模拟具有两个或多个原子的分子结构的气体的动力学行为，则需要在气体动理论层面重新建模。对于分子气体，两体碰撞保持动量守恒，但是能量却可以在各形式自由度和能级之间进行交换。因此，相比于仅具有剪切黏性的单原子气体，分子气体还存在体积黏性。另外，除了由分子平动引起的平动热导率外，还有由分子转动和振动引起的内部热导率。复杂的物理过程和更多的输运参数都给分子气体动理论建模带来挑战。

分子的内能可以表示为转动能量、振动能量以及电离能之和，分子气体内部自由度导致了更加复杂的碰撞项。王承书和乌伦贝克从量子力学角度出发建立了WCU方程，定义 $ f_i(t,{\boldsymbol{x}},{\boldsymbol{v}}) $ 为第 $ i $ 个能量为 $ e_i $ 的能级的速度分布函数，并对每个能级的分布函数都写出一个玻尔兹曼方程，其中迁移项保持不变，而碰撞项 $ Q_i $ 是玻尔兹曼碰撞项的唯象扩展，

$\begin{aligned}[b] Q_i=&\sum_{jkl}\iint P_{ij}^{kl}B(\theta,{v_r}) {\rm{d}}{\boldsymbol{\varOmega}} {\rm{d}}\boldsymbol{v}_\ast \cdot \\ &\left[\omega_{ij}^{kl}f_l(\boldsymbol{v}'_{\ast})f_k(\boldsymbol{v}')-f_j(\boldsymbol{v}_{\ast})f_i(\boldsymbol{v})\right] \end{aligned}$

(61)

其中， $\omega_{ij}^{kl} = g_kg_l/(g_ig_j)$ ， $ g_i $ 表示第 $ i $ 能级的简并度。此碰撞积分表征两个“粒子”之间的碰撞，第一个粒子在碰撞前处于 $ i $ 能级，平动速度为 $ {\boldsymbol{v}} $ ，在碰撞后跃迁到 $ k $ 能级，平动速度变为 $ {\boldsymbol{v}}' $ ；第二个粒子在碰撞前处于 $ j $ 能级，平动速度为 $ {\boldsymbol{v}}_\ast $ ，在碰撞后跃迁到 $ l $ 能级，平动速度变为 $ {\boldsymbol{v}}'_\ast $ 。两粒子分别从 $ (i,j) $ 能级跃迁至 $ (k,l) $ 能级的概率为 $ P_{ij}^{kl} $ ，由于两体碰撞的总能量守恒，上述跃迁只在以下条件发生：

$ |{\boldsymbol{v}}'_r|^2 = |{\boldsymbol{v}}_r|^2+\frac{4}{m}(e_{i}+e_{j}-e_{k}-e_{l})>0 $

(62)

碰撞后的分子速度为：

$ {\boldsymbol{v}}' = \frac{{\boldsymbol{v}}+{\boldsymbol{v}}_\ast}{2}+\frac{v'_r}{2}{\boldsymbol{\varOmega}}, \; \; {\boldsymbol{v}}'_\ast = \frac{{\boldsymbol{v}}+{\boldsymbol{v}}_\ast}{2}-\frac{v'_r}{2}{\boldsymbol{\varOmega}} $

(63)

碰撞核和跃迁概率的具体表达式是非常复杂的。这里以N₂为例，在只存在转动自由度的情况下，对WCU方程的计算复杂性进行说明。在Lennard-Jones作用势下，转动能可表达为：

$ e_{i} = 2.9{i}(i+1)\; {K} $

而转动能级的简并度为：

$ g_i = 2i+1, \quad i = (0,1,\cdots,\infty) $

对照分子动力学模拟数据^[70-71]，将跃迁概率拟合为：

$ \begin{split} P_{ij}^{kl} = P_0\omega_{ij}^{kl}&\left[\frac{2e_{tot}}{5e_{tr0}}\exp(-\varDelta_1-\varDelta_2-\varDelta_3-\varDelta_4)+\right.\\ &\left.\frac{5e_{tr0}}{2e_{tot}}\exp(-\varDelta_3-\varDelta_4)\right] \end{split} $

其中， $ P_0 $ 为一常数，

$ \begin{split}& \varDelta_1 = \dfrac{|e_{i}+e_{j}-e_{k}-e_{l}|}{e_{tr0}} \\& \varDelta_2 = 2\dfrac{|e_{i}+e_{l}-e_{k}-e_{j}|}{e_{tot}} \\& \varDelta_3 = 4\dfrac{|e_{i}-e_{k}|}{e_{tr0}+e_{i}} \\& \varDelta_4 = 2\dfrac{|e_{j}-e_{l}|}{e_{tr0}+e_{j}} \\& e_{tr0} = \dfrac{m}{4}v_r^2 \\& e_{tot} = \dfrac{m}{4}v_r^2+e_{i}+e_{j} \end{split}$

显然，如果总能级数为 $ N $ ，式（61）的计算量为单原子玻尔兹曼碰撞项（式（7））的 $ N^4 $ 倍，计算时间成本巨大。

根据碰撞过程中平动动能是否守恒，WCU方程中的碰撞项可以分为弹性碰撞与非弹性碰撞。若处在 $ i $ 能级和 $ j $ 能级的两个分子发生碰撞，碰撞后仍然能分别处在 $ i $ 能级和 $ j $ 能级，则为弹性碰撞，否则为非弹性碰撞：

$ \begin{split} Q_i = &\sum_{j}\iint P_{ij}^{ij}B \left[f_j({\boldsymbol{v}}'_{\ast})f_i({\boldsymbol{v}}')-f_j({\boldsymbol{v}}_{\ast})f_i({\boldsymbol{v}})\right]{\rm{d}}{\boldsymbol{\varOmega}} {\rm{d}}{\boldsymbol{v}}_\ast + \\ &\sum_{jkl}(1-\delta_{jk}\delta_{il}) \iint P_{ij}^{kl}B \left[f_k({\boldsymbol{v}}'_{\ast})f_l({\boldsymbol{v}}')-f_j({\boldsymbol{v}}_{\ast})f_i({\boldsymbol{v}})\right]{\rm{d}}{\boldsymbol{\varOmega}} {\rm{d}}{\boldsymbol{v}}_\ast \end{split} $

对于弹性碰撞和非弹性碰撞，可以分别引入两个等效的弛豫时间^[72]：

$ \begin{split} \frac{1}{\tau_{el}} = &\sum_{j}\iint P_{ij}^{ij}B f_j({\boldsymbol{v}}_{\ast}){\rm{d}}{\boldsymbol{\varOmega}} {\rm{d}}{\boldsymbol{v}}_\ast \\ \frac{1}{\tau_{in}} = &\sum_{jkl}(1-\delta_{jk}\delta_{il})\iint P_{ij}^{kl}B f_j({\boldsymbol{v}}_{\ast}){\rm{d}}{\boldsymbol{\varOmega}} {\rm{d}}{\boldsymbol{v}}_\ast \end{split} $

而将WCU碰撞项形式地改写为：

$ \begin{split} Q_i = &\sum_{j}\iint P_{ij}^{ij}B f_j({\boldsymbol{v}}'_{\ast})f_i({\boldsymbol{v}}'){\rm{d}}{\boldsymbol{\varOmega}} {\rm{d}}{\boldsymbol{v}}_\ast + \\ &\sum_{jkl}(1-\delta_{jk}\delta_{il})\iint P_{ij}^{kl}B f_k({\boldsymbol{v}}'_{\ast})f_l({\boldsymbol{v}}'){\rm{d}}{\boldsymbol{\varOmega}} {\rm{d}}{\boldsymbol{v}}_\ast - \frac{f_i({\boldsymbol{v}})}{\tau_{el}}-\frac{f_i({\boldsymbol{v}})}{\tau_{in}} \end{split} $

这两个等效的碰撞时间，表征着能量在各种自由度或能级间的弛豫速度，因而与下面将讨论的分子气体输运系数、体积黏性、平动/转动热导率等，有着直接的关联。

2.1 体积黏性

体积黏性表征了气体无剪切膨胀或压缩时所受的阻力，对于稀疏气体，这种阻力来自于气体分子平动动能与内能的交换过程：在膨胀或者压缩过程中，压力做功会立刻改变平动动能，但在一定的延迟后才会通过非弹性碰撞来改变内能。为方便起见，在下述推导中，假设只存在一个转动能级。

引入非弹性碰撞所对应的平动能与转动能交换的弛豫时间 $ \tau_r $ ，在空间均匀系统中，转动温度 $ T_r $ 的弛豫过程由Jeans-Landau方程描述：

$ \frac{\partial{T_r}}{\partial t} = \frac{T_t-T_r}{\tau_r} $

(64)

式中， $ T_t $ 和 $ T_r $ 分别为平动和转动温度，而平衡态温度T可表达为：

$ T = \frac{3T_t+d_rT_r}{3+d_r} $

(65)

忽略剪切黏性和热传导的影响，通过膨胀压缩过程中能量守恒可得到：

$ \begin{split} \frac{3nk_B}{2}\frac{{\rm{D}}T_t}{{\rm{D}}t} = -p_t\nabla\cdot{\boldsymbol{u}}-\frac{d_r}{2}nk_B\frac{T_t-T_r}{\tau_r}\\ \frac{d_r}{2}nk_B\frac{{\rm{D}}T_r}{{\rm{D}}t} = \frac{d_r}{2}nk_B\frac{T_t-T_r}{\tau_r} \end{split} $

(66)

式中， ${\rm{D}}/{\rm{D}}t$ 为物质导数，推导得到：

$ \frac{{\rm{D}}}{{\rm{D}}t}(T_t-T_r) = -\frac{2}{3}T_t\nabla\cdot{\boldsymbol{u}}-\frac{3+d_r}{3}\frac{T_t-T_r}{\tau_r} $

(67)

将 $ T_t-T_r $ 基于小量 $ \tau_r $ 展开为 $T_t-T_r = T_t\tau_r+{\rm{O}}(\tau^2_r)$ ，可知式（67）左侧的时间导数项较右侧为更高阶小量，于是得到：

$ T_t-T_r\approx-\frac{2\tau_r}{3+d_r}{T_t\nabla\cdot{\boldsymbol{u}}}\approx-\frac{2\tau_r}{3+d_r}{T\nabla\cdot{\boldsymbol{u}}} $

(68)

这说明，内能弛豫的影响将压力从原始压力 $ p_0 = $ $ nk_BT $ 改变为平动压力 $ p_t = nk_BT_t $ 。结合式（65）可得：

$ p_t = nk_BT_t = p_0\left[1-\frac{2d_r}{(3+d_r)^2}\tau_r\nabla\cdot{\boldsymbol{u}}\right] $

(69)

式中速度散度前的系数即为体积黏性 $ {\mu_b} $ 。在此，可以引入无量纲转动碰撞数 $ Z $ ：

$ Z = \frac{3}{d_r+3}\frac{\tau_r}{\tau} $

(70)

其中由分子平动引起的平均碰撞时间定义为 $ \tau = {\mu}/{n_0k_BT_t} $ 。于是体积黏性与剪切黏性的比值有如下关系：

$ \frac{\mu_b}{\mu} = \frac{2d_rZ}{3(d_r+3)} $

(71)

从式（70）和式（71）可知，体积黏性正比于能量交换时间。对于单原子气体，若不考虑电离，不存在内部能量的交换，因此体积黏性为零。而对于一些分子气体，例如CO₂，体积黏性可以是剪切黏性的上千倍^[73]。但是，需要强调的是，式（68）成立的条件是 $ \tau_r $ 很小，即能量弛豫的时间远远小于系统的特征时间。如果系统特征频率 $ \varpi $ 过高，则等效的体积黏性一般与系统频率有如下近似关系^[74-77]：

$ {\mu}_b = 2n_0k_BT_t\dfrac{\dfrac{{d_r} {\tau_r}}{1+\varpi^2\tau^2_r}}{\Bigg(3+\dfrac{d_r}{1+\varpi^2\tau^2_r}\Bigg)^2} $

(72)

也就是说，在 $ \varpi\tau_r\rightarrow \infty $ 极限下，体积黏性趋向零而不是无穷大，这是因为能量交换时间过长以至于系统无法感知，此时可以认为内部自由度是冻结的。因此，与式（46）中提到的剪切黏性和热导率的弛豫速率是本质而剪切黏性和热导率是表象一样，这里Jeans-Landau方程里的能量交换也是本质，而体积黏性只是连续流下的一个表象。

2.2 热导率

分子气体的热流由两部分组成，一部分是由分子的平动引起的动能传递，一部分是由分子的扩散引起的内能的传递。由于气体分子迁移过程中存在动能和内能交换，这一非弹性碰撞改变了能量在不同自由度下的分布，从而间接改变了动能传递与内能随扩散传递的结果。因而，平动热流和转动热流的弛豫过程是相互耦合的。于是，在空间均匀系统中，类似于式（46），平动热流 $ {\boldsymbol{q}}_{t} $ 和转动热流 $ {\boldsymbol{q}}_{r} $ 的弛豫过程一般可以写成如下形式：

$ \begin{split}& \frac{\partial{\boldsymbol{q}}_{t}}{\partial{t}} = -\frac{p}{\mu}( A_{tt}{\boldsymbol{q}}_{t}+A_{tr}{\boldsymbol{q}}_{r}) \\& \frac{\partial{\boldsymbol{q}}_{r}}{\partial{t}} = -\frac{p}{\mu}( A_{rt}{\boldsymbol{q}}_{t}+A_{rr}{\boldsymbol{q}}_{r}) \end{split} $

(73)

而通过Chapman-Enskog展开得到的连续流下的平动热导率 $ \kappa_{tr} $ 和转动热导率 $ \kappa_{rot} $ 可表示为：

$ \begin{split}& A_{tt}\kappa_{tr}+A_{tr}\kappa_{rot} = \frac{5}{2}\frac{k_B\mu}{m} \\& A_{rt}\kappa_{tr}+A_{rr}\kappa_{rot} = \frac{d_r}{2}\frac{k_B\mu}{m} \end{split} $

(74)

将方程（74）代入方程（1），就可以通过热弛豫系数求解如下方程，得到对应的Eucken系数：

$ \begin{bmatrix} f_{tr} \\ f_{rot} \end{bmatrix} = {\begin{bmatrix} 3A_{tt} & d_rA_{tr} \\ 3A_{rt} & d_rA_{rr} \end{bmatrix}}^{-1} \begin{bmatrix} 5 \\ d_r \end{bmatrix} $

(75)

基于WCU方程（61），Manson和Manchick通过修正的Chapman-Enskog展开推导出四个无量纲热弛豫系数的具体形式：

$ \begin{split}& A_{tt} = \frac{2}{3}+\frac{5}{6Z}\frac{d_r}{3+d_r} \\& A_{rr} = \frac{\mu}{\rho{D'}}+\frac{3}{2(3+d_r)Z} \\& A_{tr} = -\frac{5}{2(3+d_r)Z} \\& A_{rt} = -\frac{d_r}{2(3+d_r)Z} \end{split} $

(76)

其中 $ D' $ 是有效扩散系数，如果极性分子发生强共振碰撞，该系数可能与自扩散系数 $ D $ 有明显差异，例如HCl的有效扩散系数满足 $ \;\rho{D}'/\mu = 0.89 $ ，而 $\; \rho{D}/\mu = 1.33 $ 。根据式（75），忽略 $ 1/Z $ 的高阶项，可将平动Eucken因子和转动Eucken因子写为^[53]：

$ \begin{split}& f_{tr} = \frac{5}{2}\left[1-\frac{5d_r}{4(d_r+3)Z}\left(1-\frac{2\rho{D'}}{5\mu}\right)\right] \\& f_{rot} = \frac{\rho{D'}}{\mu}\left[ 1+\frac{15}{4(d_r+3)Z}\left(1-\frac{2\rho{D'}}{5\mu}\right) \right] \end{split} $

(77)

注意到， $ f_{tr} $ 随着转动碰撞数 $ Z $ 的减小而减小，即内能交换越频繁，平动热导率降低。这是因为在分子气体的碰撞过程中，部分平动动能转化为内能，当 $ D' $ 增加， $ Z $ 减小时，能量转化效率增强。这在物理上是合理的。

Mason与Monchick的解析解依然存在的一个问题，即无法保证通过式（76）推导而来的无量纲热弛豫系数 $ {\boldsymbol{A}} $ 与实际气体参数保持一致。其原因在于，在黏性和扩散系数确定的条件下，转动碰撞数 $ Z $ 作为公式中唯一用来调节Eucken因子的参数，只能使得总热导率可能达到与实验数据相近的数值，但无法分别独立调整平动与内能Eucken因子。与此同时，由于 $ Z $ 也是决定分子气体体积黏性的唯一调节参数，参考式（71），则总热导率与体积黏性的准确性通常无法同时被满足。

2.3 DSMC中的热弛豫速率

对于单原子气体，在粒子数和空间网格足够多的情况下，DSMC已被证明收敛于玻尔兹曼方程的解^[25]。实际上，由于DSMC中两体碰撞后粒子速度的更新方式与玻尔兹曼方程一致，式（43）在可接受的误差范围内恒成立，即DSMC可确保平动Eucken因子 $ f_{tr} $ 非常接近 $ 2.5 $ 。但对于分子气体，DSMC不一定能够准确恢复热导率。在DSMC的碰撞模型中，可变软球模型的提出是为了获得准确的剪切黏性与自扩散系数；唯象的Larsen-Borgnakke模型^[55]可以通过恢复能量交换速率^[56-58]从而准确获取体积黏性。然而迄今为止，尚无研究证明DSMC能够在分子气体模拟中正确恢复热导率。

我们通过提取DSMC中的热弛豫速率 $ {\boldsymbol{A}} $ 来分析其恢复热导率的能力。具体而言，在Larsen-Borgnakke模型中选取不同的转动碰撞数，对非极性气体N₂与极性气体HCl进行研究。参照Bird书中的公式（4.50）和公式（5.67），DSMC中非弹性碰撞概率 $\varLambda$ 与本文中定义的转动碰撞数 $ Z $ 通过Jeans-Landau方程（64）联系在一起^[24]：

$ Z = \frac{1}{\varLambda}\frac{{\alpha}(5-2\omega)(7-2\omega)}{5({\alpha}+1)({\alpha}+2)} $

(78)

式中， $ \omega $ 为黏性温度幂指数，见式（45）； $ {\alpha} $ 是可变软球模型中与斯密特数相关的参数：

$ {Sc} = \frac{\mu}{\rho{D}} = \frac{5(2+{\alpha})}{3(7-2\omega){\alpha}} $

(79)

使用开源软件SPARTA在空间均匀系统中提取DSMC方法中的热弛豫系数矩阵 ${\boldsymbol{A}}$ 。将一百万个DSMC模拟粒子均匀地分布在边长为10 nm的立方体内，所有的边界设置为周期性边界条件。气体数密度为 $ {n_0 = 2.69\times10^{25}} $ $\text{m}^{-3}$ ，温度为 $ T_0 = 300 $ K。如图3（a）所示，在模拟初始时刻，沿 $ x $ 轴正方向运动的粒子的速度分布函数为200 K温度下的麦克斯韦分布，沿 $ x $ 轴负方向运动的粒子速度分布函数为 400 K温度下的麦克斯韦分布。类似地，沿 $x$ 轴正方向运动的粒子的转动能满足200 K温度下的麦克斯韦分布，而沿反方向运动的粒子转动能满足400 K温度下的麦克斯韦分布，见图3（b）。这样的设置是为了建立初始平动和转动热流分布。随后记录平动和转动热流随时间的演化，直至其在系统中按照熵增原理达到热平衡后消失。图3（c）展示的是通过100次系综平均得到的平滑计算结果。

图3（d）给出N₂的热流变化率的演化过程，可以看出平动热流的时间导数随时间先增加后减小，这可以解释为平动热流与转动热流之间存在强耦合效应，即在平动热流衰减的同时得到了转动热流的补偿，且根据式（73）可知， $ A_{tr} $ 必定小于零。

图4显示通过最小二乘法求解线性回归问题（式73）得到的热流弛豫速率A。可以看出，在相同的转动碰撞数Z和斯密特数下，N₂与HCl气体的热弛豫速率几乎相同。另外，当 $ Z\rightarrow \infty $ 时，有 $ A_{tt} = 2/3 $ ， $A_{rr} = {Sc}$ ， $ A_{tr} = A_{rt} = 0 $ ，可以得到与式（3）一致的结果。图中虚线为Mason和Manchick的结果（式（76）），由于忽略了 $ 1/Z $ 的高阶项，其结果在 $ 1/Z $ 较大时与DSMC的结果出现明显偏差。

图5给出根据式（74）和式（75）计算得到的Eucken式因子，并与常温下的实验结果对比。对于N₂，当施密特数取 $ 1/1.33 $ 时，通过选择特定的碰撞数 $1/\varLambda\approx3.5$ ，DSMC可以恢复实验测量得到的总热导率，但由于缺少相应实验数据，我们无法分辨其是否能正确恢复平动与转动热导率分量。若此温度下N₂的体积黏性要求 $1/\varLambda\not\approx3.5$ ，则DSMC给出的总热导率必然存在误差。对于极性气体HCl，实验中常温下测得的总Eucken因子约为 $ 1.7 $ ，但是从图5（a）中可见，如果取施密特数为 $ 1/1.33 $ ，从DSMC中导出的总Eucken因子最小只能达到 $ 1.8 $ 。因此，无论如何选择碰撞数，DSMC都不能恢复总热导率。

从式（77）可得，当气体种类确定时，转动碰撞数决定了体积黏性，因此可以通过调整DSMC中的扩散系数以恢复总热导率，见图5（b）。然而，即便在体积黏性与总热导率均与实验值相同的情况下，方程（77）的近似解与使用Larsen-Borgnakke模型的DSMC方法得到的平动热导率 $ f_{tr} $ 并不一致，且两者的准确性都尚不确定。对于HCl气体，极性分子的共振相互作用^[51]导致 $ \rho{}D'/\mu $ 从 $ 1.33 $ 降低至 $ 0.89 $ ，为了恢复实验测量的总热导率，从图5（b）可知，DSMC的非弹性碰撞概率约等于1，达到了平动能与内能交换的极限，因此并不是一个满足物理要求的值。若选择合理的非弹性碰撞概率，则将要求DSMC中的有效扩散系数远小于自扩散系数。这将带来另一个问题：对于多组分气体，扩散系数与热导率同等重要，而Larsen-Borgnakke模型难以同时正确恢复这两个输运系数，亟需改进。

3 分子气体的动理论模型及精度

WCU方程的每一个内能能级都对应一个速度分布函数，因而计算代价巨大。对于强激波问题，N₂的能级数可达到70，给数值求解带来了极大的困难^[71]。因此对WCU方程进行简化是非常必要的，而简化的基本原则与单原子气体相同，即需要满足必要的守恒律和平衡态分布，还需要尽量使各物理量的输运系数与WCU方程保持一致，如扩散系数、剪切黏性、体积黏性、热导率及其在各自由度下的分量等。更接近本质的说法是，各输运系数对应的物理量的弛豫速率与实际气体保持一致。学者们提出了许多描述分子气体稀薄效应的动理学模型方程，如Rykov模型和ES-BGK模型等。与在单原子气体中Gross-Jackson模型可以视为所有模型的源头类似，在分子气体中可以认为Hanson-Morse模型^[78]是所有动理论模型的源头。

3.1 Hanson-Morse线性化模型

此模型在线性化情况下逼近WCU方程（61）的碰撞项。这种构造方法类似于单原子气体的Gross-Jackson模型^[79]，即通过对玻尔兹曼碰撞项进行正交多项式展开，并保证各多项式对应的宏观量的弛豫速率与WCU方程一致而获得。对于分子气体，关于分子速度的多项式仍由麦克斯韦分子的本征函数（24）给出；关于内部能量的多项式的前两项取为1和 $ \varepsilon = e_i-\bar{e} $ ，其中 $ \bar{e} $ 为平均内部能量^[52]。

分子气体的平衡态速度分布函数可表示为 $ f_i = o_iF_{eq} $ ，其中

$ o_i = \frac{g_i\exp\Bigg(\dfrac{-e_i}{k_BT}\Bigg)}{\displaystyle\sum_j g_j\exp\Bigg(\dfrac{-e_j}{k_BT}\Bigg)} $

(80)

表征了内能为 $ e_i $ 的分子比例。在线性化问题中，通常考虑宏观速度为零的全局平衡态，并将速度分布函数写为如下形式：

$ f_i(t,{\boldsymbol{x}},{\boldsymbol{\xi}}) = o_if_{eq}({\boldsymbol{\xi}})\left[1+\phi_i(t,{\boldsymbol{x}},{\boldsymbol{\xi}})\right] $

(81)

其中每个内部能级对应的扰动分布函数 $ \phi_i $ 满足 $ |\phi_i|\ll1 $ 。注意到为与上文的单原子线性化模型保持一致，这里的速度和空间坐标分别通过最概然速率 $ v_m $ 和特征尺寸 $ \ell $ 进行无量纲化，时间通过 $ \ell/v_m $ 归一化，加速度通过 $ v_m^2/\ell $ 无量纲化，内能 $ e_i $ 通过 $ k_BT_0 $ 无量纲化。归一化后的扰动物理量可以通过对扰动分布函数求矩得到：

$ \begin{array}{c} {n}_i = \displaystyle\int{f_{eq}}\phi_i{\rm{d}}{\boldsymbol{\xi}} \\ \bar{e} = \displaystyle\sum_io_ie_i \\ [{\boldsymbol{u}},\sigma_{{\alpha}\beta}] = \displaystyle\sum_io_i\int{f_{eq}}\phi_i[{\boldsymbol{\xi}},\xi_{\langle {\alpha}}\xi_{\beta\rangle}]{\rm{d}}{\boldsymbol{\xi}} \quad\\ T_{t} = \displaystyle\sum_io_i\int{f_{eq}}\phi_i \left(\frac{2}{3}\xi^2-1\right){\rm{d}}{\boldsymbol{\xi}} \\ {\boldsymbol{q}}_{t} = \displaystyle\sum_io_i\int{f_{eq}}\phi_i \left(\xi^2-\frac{5}{2}\right){\boldsymbol{\xi}}{\rm{d}}{\boldsymbol{\xi}} \\ [T_{r},{\boldsymbol{q}}_{r}] = \displaystyle\sum_io_i\int{f_{eq}}\phi_i\varepsilon_i \left[\frac{2}{d_r},{\boldsymbol{\xi}}\right]{\rm{d}}{\boldsymbol{\xi}} \end{array} $

(82)

同时，总扰动密度为 ${n} = \displaystyle\sum_i{n}_io_i$ 。

若加速度为小量，WCU方程的线性化形式为：

$ \frac{\partial \phi_i}{\partial{}t} +{{\boldsymbol{\xi}}}\cdot\frac{\partial \phi_i}{\partial {\boldsymbol{x}}} -2{{\boldsymbol{a}}}\cdot{\boldsymbol{\xi}} = \frac{\ell}{v_m}{J}_i $

(83)

其中 $ {J}_i $ 的形式与式（21）类似，不再赘述。根据Hanson-Morse模型^[78]，该线性化碰撞项可近似为^[78]：

$ \begin{split} {J}_i = &-J_{030}(n_i+2{\boldsymbol{\xi}}\cdot{\boldsymbol{u}}-\phi_i) +{2\sigma_{{\alpha}\beta}}(J_{020}-J_{030}) \xi_{\langle {\alpha}}\xi_{\beta\rangle}+\\ & T_{t}\left[ (-J_{030}+J_{100})\left(\xi^2-\frac{3}{2}\right)-\frac{3}{d_r}J_{100}{\varepsilon} \right] +\\ & 2{\boldsymbol{\xi}}\cdot{\boldsymbol{q}}_{t}\left[ (J_{110}-J_{030})\frac{2}{5}\left(\xi^2-\frac{5}{2}\right) -J_{011}^{110}\sqrt{\frac{4}{5d_r}}{\varepsilon} \right]+\\ & T_{r}\left[ -J_{100}\left(\xi^2-\frac{3}{2}\right) +\left( \frac{3}{d_r}J_{100}-{J}_{030} \right){\varepsilon} \right]+\\& 2{\boldsymbol{\xi}}\cdot{\boldsymbol{q}}_{r}\left[ -{J}_{011}^{110}\sqrt{\frac{4}{5d_r}}\left(\xi^2-\frac{5}{2}\right) +2(J_{011}-J_{030})\frac{{\varepsilon}}{d_r} \right] \end{split} $

Hanson-Morse模型假设各内部能级的弛豫时间相同，则与式（76）对应的弛豫率为：

$ \begin{split}& \tau = \frac{\mu}{p_t} \\& J_{020} = -\tau^{-1} \\& J_{030} = -\frac{3}{2}\tau^{-1} \\& J_{100} = -\frac{1}{Z\tau}\frac{d_r}{3+d_r} \\& J_{110} = -\frac{2}{3\tau}-\frac{5}{6Z\tau}\frac{d_r}{3+d_r} \\& J_{011}^{110} = -\sqrt{\frac{5}{4d_r}}\frac{1}{Z\tau}\frac{d_r}{3+d_r} \\& J_{011} = -\frac{\mu}{\rho{D'}}-\frac{3}{2(3+d_r)Z}\end{split} $

(84)

在自发瑞利-布里渊散射中，如果模型中包含应力偏量 $ \sigma_{{\alpha}\beta} $ ，则该模型称为Tenti-S7模型^[80-81]，其缺点在于不能得到实验散射光谱。Tenti-S6模型通过移除应力偏量 $ \sigma_{{\alpha}\beta} $ ，即，令模型中的 $ \sigma_{{\alpha}\beta} = 0 $ ，获得与实验一致的散射光谱。实际上，如果不考虑内部能级，Tenti-S7和Tenti-S6模型分别退化为单原子ES-BGK模型和Shakhov模型，而在大多数情况下Shakhov模型的精度高于ES-BGK模型。

当在Hanson-Morse模型中移除应力偏量，为了复现应力偏量、能量弛豫、热流的弛豫过程，式（84）中的 $ J_{030} $ 应该改成 $ -1/\tau $ 。同时，将系数 $ J_{011} $ 修改为以下形式^[80,82]，从而保证正确恢复热导率。

$ J_{011} = -\dfrac{2d_r}{3(3+d_r)\tau} \dfrac{\dfrac{1}{5}(3+d_r)+\dfrac{6+d_r}{4Z}+\dfrac{9}{16Z^2}f_u} {\dfrac{2}{15}f_u\left(3+d_r\right)+\dfrac{d_r}{6Z}f_u-1} $

(85)

3.2 Hanson-Morse模型的约化及非线性化

为了减小计算量，引入如下两个约化速度分布函数：

$ g = \sum_i{o_i}f_{eq}{\phi_i}, \;\;r = \sum_i{o_i}f_{eq}{\phi_i}(e_i-\bar{e}) $

(86)

消去内能自由度，则Hanson-Morse模型可以简化为：

$ \begin{split}& \frac{\partial g}{\partial{}t} +{{\boldsymbol{\xi}}}\cdot\frac{\partial g}{\partial {\boldsymbol{x}}} -2{{\boldsymbol{a}}}\cdot{\boldsymbol{\xi}} = \frac{\ell}{v_m}{J}_g \\& \frac{\partial r}{\partial{}t} +{{\boldsymbol{\xi}}}\cdot\frac{\partial r}{\partial {\boldsymbol{x}}} -2{{\boldsymbol{a}}}\cdot{\boldsymbol{\xi}} = \frac{\ell}{v_m}{J}_r \end{split} $

(87)

其中：

$ \begin{split} {J}_g =& -J_{030}\left[ {n}+2{\boldsymbol{\xi}}\cdot{\boldsymbol{u}}+{\boldsymbol{\xi}}\cdot{\boldsymbol{q}}_{t} \frac{4}{15}\left(\xi^2-\frac{5}{2}\right)-g\right]-\\ &J_{030}T_{t}\left(\xi^2-\frac{3}{2}\right)+{2\sigma_{{\alpha}\beta}}(J_{020}-J_{030}) \xi_{\langle {\alpha}}\xi_{\beta\rangle} -\\ &(T_{r}-T_{t})J_{100}\left(\xi^2-\frac{3}{2}\right)-\\ &{\boldsymbol{\xi}}\cdot\left[{\boldsymbol{q}}_{t} \frac{2d_r}{3Z\tau(3+d_r)} +2{\boldsymbol{q}}_{r} {J}_{011}^{110}\sqrt{\frac{4}{5d_{r}}} \right]\left(\xi^2-\frac{5}{2}\right) \end{split} $

(88)

$ \begin{split} {J}_r = & -{J}_{030}\left(\frac{d_r}{2}T_{r}+2{\boldsymbol{\xi}}\cdot{\boldsymbol{q}}_{r}-r\right)-\\ &\frac{3}{2}J_{100}(T_{t}-T_{r})-\\ & 2{\boldsymbol{\xi}}\cdot{\boldsymbol{q}}_{t}J_{011}^{110}\sqrt{\frac{d_{r}}{5}} +2{\boldsymbol{\xi}}\cdot{\boldsymbol{q}}_{r} J_{011} \end{split} $

(89)

至此，注意到单原子气体中Gross-Jackson模型与BGK、ES-BGK、Shakhov等模型的关系，可以通过类似于式（60）中展示的非线性化程序构造出分子气体的动理学模型。即，首先写出扰动量 $ g $ 和 $ r $ 对应的速度分布函数，

$ \begin{split} &G(t,{\boldsymbol{x}},{\boldsymbol{v}}) = \sum_i{o_i}{f_i}\\ &R(t,{\boldsymbol{x}},{\boldsymbol{v}}) = \sum_i{o_i}{f_i}e_i \end{split} $

(90)

并定义宏观量如下：

$ \begin{split} &[n,{\boldsymbol{u}},T_t] = \int{}\left[1,\frac{{\boldsymbol{v}}}{n}, \frac{mc^2}{3nk_B}\right]G{\rm{d}}{\boldsymbol{v}} \\ &[p_{ij},{\boldsymbol{q}}_{t}] = \int{}\left[mc_i{c_j},\frac{{mc^2}{\boldsymbol{c}}}{2}\right]G{\rm{d}}{\boldsymbol{v}} \\ &[T_r,{\boldsymbol{q}}_{r}] = \int{}\left[\frac{2}{d_rnk_B},{\boldsymbol{c}}\right]R {\rm{d}}{\boldsymbol{v}}\end{split} $

(91)

总热流定义为 $ {\boldsymbol{q}} = {\boldsymbol{q}}_{t}+{\boldsymbol{q}}_{r} $ ，平动压力、转动压力和总压力分别定义为 $p_t = nkT_t {\text{、}} p_r = nkT_r$ 和 $ p = nkT $ 。注意此处 $ G $ 可以线性化为 $ G = f_{eq}(1+g) $ ，而 $ R $ 应线性化为 $ R = (d_rk_BT/2)f_{eq}(1+g+r) $ 。其次，把速度分布函数的演化过程形式地写为：

$ \begin{split} \frac{\partial G}{\partial t}+{\boldsymbol{v}}\cdot \frac{\partial G}{\partial {\boldsymbol{x}}} +{\boldsymbol{a}}\cdot \frac{\partial G}{\partial {\boldsymbol{v}}}& = \underbrace{\frac{G_t-G}{ \tau }}_{{\rm{弹性}}}+\underbrace{\frac{G_r-G_t}{ Z\tau}}_{{\rm{非弹性}}} \\ \frac{\partial R}{\partial t}+{\boldsymbol{v}}\cdot \frac{\partial R}{\partial {\boldsymbol{x}}} +{\boldsymbol{a}}\cdot \frac{\partial R}{\partial {\boldsymbol{v}}}& = \underbrace{\frac{R_t-R}{ \tau }}_{{\rm{弹性}}}+\underbrace{\frac{R_r-R_t}{ Z\tau}}_{{\rm{非弹性}}} \end{split} $

(92)

最后，根据式（88）和式（89）中的碰撞项，类比式（60），写出四个参考分布函数 $ G_t $ 、 $ G_r $ 、 $ R_t $ 和 $ R_r $ 的形式。

例如，考虑Hanson-Morse模型中不出现应力偏量。式（88）中的第一行和第二行即为线性化Shakhov模型的碰撞项，可以非线性化为：

$\begin{split}& \frac{p_t}{\mu} \left\{ F_{eq}(T_t)\left[1+\frac{2m{\boldsymbol{q}}_{t} \cdot {\boldsymbol{c}}}{15k_BT_tp_t}\left(\frac{mc^2}{2k_BT_t}-\frac{5}{2}\right)\right]-G \right\}\\&\qquad= \frac{p_t}{\mu} (G_t-G) \end{split} $

(93)

结合非线性模型（式92）中的非弹性碰撞项 $ (G_r-G_t)/Z\tau $ ，可以选取

$ \begin{split}& G_r = F_{eq}(T)\left[1+\frac{2m{\boldsymbol{q}}' \cdot {\boldsymbol{c}}} {15k_BTp}\left(\frac{mc^2}{2k_BT}-\frac{5}{2}\right)\right] \\& {\boldsymbol{q}}' = \left[1-\frac{5d_r}{2(3+d_r)}\right]{\boldsymbol{q}}_{t} +\frac{15}{2(3+d_r)}{\boldsymbol{q}}_{r} \end{split} $

(94)

以在线性化情况下退化为式（88）中的后两项，其中能量交换由 $ F_{eq}(T) $ 和 $ F_{eq}(T_t) $ 的差距而来，热流松弛由 $ {\boldsymbol{q}}' $ 决定。

3.3 非线性模型

然而，历史并没有选择非线性化路线，而是由许多学者各自提出了非线性动理学模型^{[33-34,50,61,83-87]}。下面先介绍主要的Rykov模型和ES-BGK模型，并揭示与Hanson-Morse模型的内在联系，最后介绍吴模型。

3.3.1 Rykov模型

原始的Rykov模型^[61]仅适用于无振动模态的双原子气体，但是构造思想可以直接拓展至多原子分子气体建模^[83]。在该模型中，四个参考速度分布函数定义为：

$ \begin{split}& G_t = F_{eq}(T_t)\left[1+\frac{2m{\boldsymbol{q}}_{t} \cdot {\boldsymbol{c}}}{15k_BT_tp_t}\left(\frac{mc^2}{2k_BT_t}-\frac{5}{2}\right)\right] \\& G_r = F_{eq}(T)\left[1+\omega_0\frac{2m{\boldsymbol{q}}_{t} \cdot {\boldsymbol{c}}} {15k_BTp}\left(\frac{mc^2}{2k_BT}-\frac{5}{2}\right)\right] \\& R_t = \frac{d_rk_BT_r}{2}G_t+F_{eq}(T_t)(1-{Sc})\frac{m{\boldsymbol{q}}_{r} \cdot {\boldsymbol{c}}}{p_t} \\& R_r = \frac{d_rk_BT}{2}G_r+F_{eq}(T)\omega_1(1-{Sc})\frac{m{\boldsymbol{q}}_{r} \cdot {\boldsymbol{c}}}{p} \end{split} $

(95)

选取弹性碰撞弛豫时间 $ \tau = {\mu}/{p_t} $ 以及适当的转动碰撞数 $ Z $ ，可保证恢复正确的剪切黏性和体积黏性。四个热流弛豫系数为：

$ \begin{split} &A_{tt} = \frac{2}{3}\left(1+\frac{1-\omega_0}{2Z}\right) \\ &A_{rr} = {Sc}+\frac{(1-{Sc})(1-\omega_1)}{Z} \\ &A_{tr} = A_{rt} = 0 \end{split} $

(96)

因此通过式（75）可以得到Eucken因子的表达式如下：

$ \begin{split} f_{tr}& = \frac{5}{2}\left(1+\frac{1-\omega_0}{2Z}\right)^{-1} \\ f_{rot}& = \left[{Sc}+\frac{(1-{Sc})(1-\omega_1)}{Z}\right]^{-1} \end{split} $

(97)

这说明，在Rykov模型中，可通过调整自由参数 $ \omega_0 $ 和 $ \omega_1 $ 修改Eucken系数。

当 $ Z\rightarrow \infty $ ，该模型退化为单原子气体的Shakhov模型。如果将其线性化，可以得到类似于第3.2节中的约化Hanson-Morse模型。不同之处在于空间均匀系统中，Rykov模型中平动和转动热流的弛豫不会相互干扰，这是因为式（74）中 $ A_{tr} = A_{rt} = 0 $ 。

3.3.2 ES-BGK模型

Holway通过最大熵原理提出了ES-BGK模型，该模型保留了标准BGK模型的数学简便性，并且能够恢复正确的普朗特数，同时满足熵增定理^[33-34]。速度分布函数的演化方程为：

$ \begin{split}& \frac{\partial G}{\partial t}+{\boldsymbol{v}}\cdot \frac{\partial G}{\partial x} +{\boldsymbol{a}}\cdot \frac{\partial G}{\partial {\boldsymbol{v}}} = \frac{1}{ \tau }\left(f_r^{ES}-G\right) \\& \frac{\partial R}{\partial t}+{\boldsymbol{v}}\cdot \frac{\partial R}{\partial x}+{\boldsymbol{a}}\cdot \frac{\partial R}{\partial {\boldsymbol{v}}} = \frac{1}{ \tau }\left(\frac{d_r}{2}k_BT_{rel}f_r^{ES}-R\right) \end{split} $

(98)

其中，弛豫温度 $ T_{rel} $ 定义为转动温度与总温度的加权平均值：

$ T_{rel} = \left(1-\frac{1}{Z}\right){T_r}+\frac{T}{Z} $

(99)

参考分布函数 $ f_r^{ES} $ 由式（50）给出，但是此处

$ \lambda_{ij} = \left(1-\frac{1}{Z}\right) \left[(1-b)\frac{k_BT_t}{m}\delta_{ij}+b\frac{p_{ij}}{nm}\right] +\frac{k_BT}{Zm}\delta_{ij} $

(100)

当 $ Z\rightarrow \infty $ ，转动和平动能量交换消失，因此模型方程退化为单原子气体的ES-BGK方程。

取弹性碰撞弛豫时间 $\tau = \mu/p_t Pr$ 及适当的转动碰撞数 $ Z $ ，则分子气体的剪切和体积黏性自动满足要求。其中普朗特数可以通过调整 $ b $ 的取值得到：

$ Pr = \frac{1}{1-b+\frac{b}{Z}} $

(101)

为使参考分布函数有界且在物理上有意义，要求 $Z\geqslant 1$ 和 $-1/2\leqslant {b} \leqslant 1$ 。因此，普朗特数范围为 $2/3 \leqslant Pr \leqslant +\infty$ 。对于大部分双原子分子来说，普朗特数大约为 $ 5/7 $ ，因此可以保证要求。

与Rykov模型一样，ES-BGK模型的平动热流和转动热流的弛豫过程互不干涉；从其弛豫率可知，Eucken因子可以表达为：

$ \begin{split} &f_{tr} = \frac{5}{3Pr} \\ &f_{rot} = \frac{1}{Pr} = \frac{3+d_r}{5+d_r}f_{eu} \end{split} $

(102)

注意到，对于多原子ES-BGK模型，平动Eucken因子与转动Eucken因子的比值始终为 $ 5/3 $ ，并且Eucken因子完全由普朗特数决定。

3.3.3 吴模型

上述气体动理学模型的弛豫时间 $ \tau $ 都与分子速度无关，与玻尔兹曼方程的弛豫时间不相符。因此，在正激波问题中，Rykov模型计算得到的结果可能会出现波前温升的误差，且随着马赫数的增加，这一偏离显著增长。由此，为了提高模型方程准确性，吴雷等首先将Rykov模型中的弹性碰撞项 $ ({G_t}-G)/\tau $ 替换为单原子玻尔兹曼方程的碰撞项 $ Q(G) $ ，以使弹性碰撞弛豫时间是分子速度的函数^[83]；然后，又考虑引入更为准确的热流弛豫过程^[60]，从而在相同参数条件下，得到更为符合DSMC的模拟结果。

吴模型的演化方程为：

$ \begin{split} &\frac{\partial G}{\partial t}+{\boldsymbol{v}}\cdot \frac{\partial G}{\partial x}+{\boldsymbol{a}}\cdot \frac{\partial G}{\partial {\boldsymbol{v}}} = Q(G)+\frac{G_r-G_t}{ Z\tau } \\& \frac{\partial R}{\partial t}+{\boldsymbol{v}}\cdot \frac{\partial R}{\partial x} +{\boldsymbol{a}}\cdot \frac{\partial R}{\partial {\boldsymbol{v}}} = \frac{R_t'-R}{ \tau }+\frac{R_r-R_t}{ Z\tau } \end{split} $

(103)

式中， $ R_t' = (d_r/2)k_BT_r(\tau Q+G) $ ，弛豫时间为 $ \tau = \mu/p_t $ ，对应的参考速度分布函数为：

$ \begin{split}& G_t = F_{eq}(T_t)\left[1+\frac{2m{{\boldsymbol{q}}_{t}} \cdot {\boldsymbol{c}}}{15k_BT_tp_t}\left(\frac{mc^2}{2k_BT_t}-\frac{5}{2}\right)\right] \\& G_r = F_{eq}(T)\left[1+\frac{2m{\boldsymbol{q}}' \cdot {\boldsymbol{c}}} {15k_BTp}\left(\frac{mc^2}{2k_BT}-\frac{5}{2}\right)\right] \\& R_t = \frac{d_rk_BT_r}{2}G_t \\& R_r = \frac{d_rk_BT}{2}G_r+F_{eq}(T)\frac{m{\boldsymbol{q}}'' \cdot {\boldsymbol{c}}}{p} \end{split} $

(104)

为了恢复热流弛豫过程（式（73）），选取 $ {{\boldsymbol{q}}'} $ 和 $ {{\boldsymbol{q}}''} $ 如下：

$ \begin{split}& {\boldsymbol{q}}' = \left[-3Z( A_{tt}-\frac{2}{3})+1\right]{\boldsymbol{q}}_{t}-3ZA_{tr}{\boldsymbol{q}}_{r} \\& {\boldsymbol{q}}'' = -Z\left[A_{rt} q_t+(A_{rr}-1){\boldsymbol{q}}_{r}\right] \end{split} $

(105)

当 $ Z\rightarrow \infty $ ，吴模型完全退化为单原子玻尔兹曼方程。另一方面，如果将玻尔兹曼碰撞项换回Shakhov碰撞模型，在线性化的情况下，吴模型则完全与第3.2节中的约化Hanson-Morse模型一致。

3.4 模型方程精度比较

本节以正激波结构、库特流、热蠕动、二维热驱动问题为例，比较不同模型描述稀薄气体流动的精度。模拟中，分子气体为N₂，转动自由度 $ d_r = 2 $ ，模型中的相关参数选为 $Z = 2.6671 {\text{、}} \omega = 0.74$ ，施密特数 ${Sc} = 1/1.33$ 。吴模型中的热弛豫系数矩阵 $ {\boldsymbol{A}} $ 来自第2.3节中的DSMC模拟，从而保证吴模型的Eucken因子为 $ f_{eu} = 1.993 $ 、 $ f_{tr} = 2.365 $ 和 $ f_{rot} = 1.436 $ ，与DSMC保持一致；此外，进一步假设这些参数不随温度变化。Rykov模型可以通过式（97）调整 $ \omega_0 $ 和 $ \omega_1 $ 恢复正确的Eucken系数。从式（102）可知，ES-BGK模型只有一个自由度来恢复总Eucken系数，即ES-BGK模型没有分辨平动Eucken系数和转动Eucken系数的机制。幸运的是，对于N₂，ES-BGK模型给出的平动和转动热导率与DSMC非常接近，见表 1中关于Eucken参数的汇总。再次强调，ES-BGK和Rykov模型中平动和转动热流的弛豫过程是相互独立的。

引入参考密度 $ n_0 $ 、参考温度 $ T_0 $ 、最概然速率 $ v_m = \sqrt{2k_BT_0/m} $ 、参考压力 $ n_0k_BT_0 $ 、参考热流 $ n_0k_BT_0v_m $ 。宏观量通过参考值进行无量纲化，且无量纲克努森数定义为：

$ {Kn} = \frac{\mu(T_0)}{n_0\ell}\sqrt{\frac{\text{π}}{2mk_BT_0}} $

(106)

用于表征流动的稀薄程度，其中 $ \ell $ 为流动特征长度。

在以下各研究问题中，气体在固体表面的反射均用完全漫反射作为边界条件，即气体反射速度分布为满足固体表面宏观量的麦克斯韦分布。在数值解法上，使用离散速度法求解动理学模型，快速谱方法求解玻尔兹曼碰撞项^[31]，开源软件SPARTA进行DSMC模拟。

3.4.1 正激波

图6给出马赫数为7的一维N₂正激波的结构。空间坐标用激波上游的分子平均自由程 $L = 16\mu/5n\sqrt{2{\text{π}}{m}k_BT}$ 对其进行无量纲化，密度与温度通过 $ (Q-Q_u)/(Q_d-Q_u) $ 归一化，其中 $Q = n {\text{、}} T_t {\text{、}} T_r$ ，下标 $ u $ 和 $ d $ 分别代表上游和下游的物理量。所有动理学模型的密度分布均与DSMC计算结果大体上相符，只有ES-BGK模型给出的密度在波前稍低一些。Shakhov模型和ES-BGK模型的平动温度在波前存在“过早温升”现象，其中ES-BGK的平动与转动温度分布偏离最大，这是因为这两个模型方程的碰撞时间与分子速度无关。而吴模型和DSMC中的碰撞时间都是分子速度的单调递增函数，因而得到最为接近的结果。图6（c、d）给出了压力与热流分布的计算结果，可见吴模型与DSMC符合得更好，而其他模型均存在显著误差，尤其是热流分布的对比。

3.4.2 库特流动

平板库特流动的上下两板温度相同，处在 $ x_2 = L/2 $ 的上板速度是 $ v_w = v_m $ ，处在 $ x_2 = -L/2 $ 的下板速度相反。图7给出各模型和DSMC的对比结果。当 ${Kn} = 0.1$ 时，所有动理学模型都给出了良好的密度分布。ES-BGK模型高估了平动温度和平动热流，却低估了转动热流。Rykov模型给出良好的温度和转动热流，却高估了平动热流。吴模型与DSMC的对比结果非常理想。当 ${Kn} = 1$ 时，吴模型同样能给出与DSMC较为一致的结果。Rykov模型方程的密度分布与DSMC的结果稍有差别，在两板中心区域略低。ES-BGK模型的转动温度也略低于其他模型。

3.4.3 麦克斯韦妖驱动的微流动

考察处在两平行静止的无限长平板间的分子气体，每个气体分子受到与其速度相关的水平外力的作用：

$a_{x_1}=a_0\left(\frac{v^2}{v^2_m}-\frac{3}{2}\right)$

(107)

在此例中，取 $ {2{a_0}{L}/{v^2_m} = 0.0718} $ ，以保证系统仅略微偏离平衡态。可以看出，当分子速率大于 $ \sqrt{3/2}v_m $ 时，其加速度为正，反之则为负。因此该加速度可以看出是一种特殊形式的麦克斯韦妖，见图8。

图8给出克努森数分别为0.1和1的模拟结果对比。吴模型与DSMC非常符合，速度剖面与热流分布的最大误差仅为 $ 2\% $ 。值得注意的是，ES-BGK模型和Rykov模型的转动热流值远大于吴模型和DSMC的结果。其原因在于，这两个模型方程中的平动热流弛豫与转动热流弛豫相解耦，即 $ A_{tr} = A_{rt} = 0 $ ，即便热导率相近，但热弛豫系数的差别仍可导致显著的宏观量的差异。这正说明了采用正确的热弛豫速率的重要性。

3.4.4 封闭方腔中的热蠕动

最后，考察二维封闭方腔内的热蠕动问题。方腔宽度与高度比值为5，固定左板和右板的温度分别为200 K和400 K，而上下两板的温度在此间线性分布。由于此问题以 $ x_2 = 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中，因碰撞模型的原因，固定体积黏性后，热弛豫系数也随之确定。在吴模型中，可以通过独立修改 $ {A_{ij}} $ 以量化热弛豫系数的影响。上文已验证吴模型的准确性，这里使用吴模型来量化这种不确定性。首先，在约束 $ {f_{tr}} $ 和 $ {f_{rot}} $ 不变的前提下修改 $ {A_{ij}} $ ，以研究热弛豫速率的影响。其次，在保持总热导率不变前提下，通过改变 $ {f_{tr}} $ 和 $ {f_{rot}} $ 来量化热导率分量的影响。

4.1 正激波

给定 $ {f_{tr}} $ 和 $ {f_{rot}} $ ，即体系中所有的热导率分量都被确定，若此时其余的输运系数也保持不变，则称此系统为宏观意义上的确定系统。然而，不同的 $ {A_{ij}} $ 的选择将可能导致结果的不确定性。

考察一维正激波问题， $ {A_{tr}} $ 、 $ {A_{rt}} $ 的取值范围分别为 $ [-5/(6Z),0] $ 和 $ [-1/(3Z),0] $ ，而 $ {A_{tt}} $ 和 $ {A_{rr}} $ 的值根据式（73）调整，以保持 $ {f_{tr}} $ 和 $ {f_{rot}} $ 不变。在DSMC中非对角线元素为 $A_{tr} = -0.201,\;A_{rt} = -0.059$ ，本例中，给定 $ {Z = 2.6671} $ ， $ {A_{rt}} $ 和 $ {A_{tr}} $ 的最小值分别为 $ -0.3124 $ 和 $ -0.1250 $ ，数值大小约为DSMC对应参数的一到两倍。研究发现，热弛豫系数A的变化导致转动温度和转动热流稍有变化，而对密度、速度和应力偏量几乎没有影响。

接着，保持总热导率不变，且选择 $A_{tr} = -5/(6Z),\; $ $ A_{rt} = -1/(3Z)$ ，研究热导率分量的改变对流动的影响。图10展示了马赫数为4时不同平动Eucken因子下吴模型的数值结果。在实际中，较小数值的平动Eucken因子 $ f_{tr} $ 是可能存在的，特别是在极性气体中，真实的平动Eucken因子可能远小于 $ 2.5 $ 。例如，水（H₂O）的平动Eucken因子为 $ f_{tr} = 1.78 $ ，甲醇分子（CH₃OH）^[53]的平动Eucken因子为 $ f_{tr} = 0.41 $ 。对比可见，所有的宏观量都随 $ f_{tr} $ 改变而有明显的变化。首先，较大的 $ {f_{tr}} $ 使得平动温度更快上升到其最大值，并在波后更快趋于下游温度，对于应力偏量和总热流也有相同的趋势。其次， $ {f_{tr}} $ 变化所导致的转动温度变化则集中于激波结构的中心处，较小的 $ {f_{tr}} $ 导致更高的转动温度。最后，较大的 $ {f_{tr}} $ 使得密度升高更快。

由此，在正激波问题中，由热弛豫系数 ${\boldsymbol{A}}$ 的变化引起的不确定性较小，而热导率分量的变化则会带来显著的不确定性。例如， $ {f_{tr} = 1.5} $ 的转动热流大约是 $ {f_{tr} = 2.5} $ 的两倍，几乎与平动热流相当。

4.2 麦克斯韦妖驱动的微流动

使用与正激波算例中相同的系数矩阵A，使 $ f_{tr} $ 和 $ f_{rot} $ 均保持不变，来考察麦克斯韦妖驱动的微流动中，热弛豫系数对流动速度和热流的影响。图11第一行图片数据为 ${Kn} = 0.2$ 的结果。可以看出，不同于正激波问题，此问题中不同热弛豫系数A对应的流动速度与热流有相当明显的差别。流动速度和平动热流的最大相对不确定性分别为 $ 16.7\% $ 和 $ 17.6\% $ 。同时可以看出，流动的中间区域出现了较大的不确定性，而壁面附近的速度滑移与热流几乎不变。

为进一步研究平动Eucken因子的影响，考察Eucken因子 $ {f_{tr} = 1.5, 2.0, 2.5} $ 的算例，其他参数为 ${Kn} = 0.2$ ， $ f_u = 1.993 $ ， $ A_{tr} = {-5/6Z} $ ， $ A_{rt} = {-1/3Z} $ 。如图11第二行图片数据所示，速度曲线与平动热流随着平动Eucken因子的变化而剧烈变化。 $ {f_{tr} = 2.5} $ 时的流动速度和热流比 $ {f_{tr} = 1.5} $ 时增加了大约 $ {68\%} $ 。在第一行图中，在确定的 $ f_{tr} $ 下，壁面速度滑移与热流不受平动Eucken因子的影响。而在第二行图中，速度与热流依赖于平动Eucken因子，例如壁面附近的速度与热流均与 $ {f_{tr}} $ 呈正相关。这说明，在此问题中，平动Eucken因子 $ {f_{tr}} $ 起主导作用。

平动Eucken因子在此问题中的重要性可得到如下解释：从方程（73）可以看出，热弛豫系数 $ {A_{tr}} $ 和 $ {A_{rt}} $ 与平动能与转动能之间的能量交换有关。当 $ {A_{tr}} $ 与 $ {A_{rt}} $ 为零时，平动热流与转动热流的弛豫过程互相解耦。例如，更多的计算结果显示，在 $ {A_{rt} = 0} $ 的情况下，转动热流恒为零。这是因为麦克斯韦妖驱动的热蠕动流动中，外力项仅能直接影响平动能，而转动能和热流则需通过能量交换才能发生改变，具体取决于 $ {A_{tr}} $ 、 $ {A_{rt}} $ 和 $ Z $ 。由于矩阵的非对角元素相对于其他元素为小量，这导致 $ q_{rot}\approx0 $ ，而平动效应（ $ q_{tr} $ 或者 $ f_{tr} $ ）占主导。

5 确定热流弛豫速率的方法

上文的研究结果揭示了热流弛豫速率不确定性对稀薄气体流动的影响，因此有必要通过获取并使用正确的相应系数，来减少或消除模型的不确定性。我们认为可以尝试以下两种方法。

5.1 分子动力学模拟

从物理上来说，只要知道了分子间的相互作用势，根据牛顿第二运动定律，就可以通过分子动力学方法模拟气体的行为。由此，输运系数能够通过Green-Kubo公式获得^[92-93]：

$ \mu = \frac{1}{V}\frac{1}{k_BT}\int_0^\infty \langle J_{xy}^p(0)J_{xy}^p(\tau) \rangle {\rm{d}}\tau $

(108)

$ \frac{4}{3}\mu+\mu_b = \frac{1}{V}\frac{1}{k_BT}\int_0^\infty \langle J_{xx}^p(0)J_{xx}^p(\tau) \rangle {\rm{d}}\tau $

(109)

$ \kappa = \frac{1}{V}\frac{1}{k_BT^2}\int_0^\infty \langle J_{x}^q(0)J_{x}^q(\tau) \rangle {\rm{d}}\tau $

(110)

其中， $ V $ 和 $ T $ 分别是模拟系综的体积和温度， $ \tau $ 为时间，尖括号表示系综平均； $ {\boldsymbol{J}}^p = m{\boldsymbol{c}}_i{\boldsymbol{c}}_i-pV{\boldsymbol{I}} $ ， $ {\boldsymbol{J}}^q = {E_i}{\boldsymbol{c}}_i $ ， $ {\boldsymbol{I}} $ 是单位矩阵。对于单原子气体， $ E_i = mc_ic_i/2 $ 得到平动热导率。对于分子气体，取 $ E_i = \varepsilon_i $ 为第 $ i $ 个分子的内部能量，则式（110）给出内部热导率。

然而，对于如何获得热流弛豫矩阵系数尚不清楚。一个可行的方法是仿照第2.3节中介绍的DSMC方法，即在分子动力学中设置带有初始热流的状态，模拟低密度气体的演化过程，然后通过拟合公式（73）获得弛豫系数。虽然目前尚未有相关工作的报道，但是原理上是可行的。

5.2 瑞利-布里渊散射

瑞利-布里渊散射指的是光在气体中传播时，被气体分子的自发密度涨落散射。由于散射光谱包含气体信息，因此可以用来无损测量声速、流速、温度和体积黏性^[94-96]。例如，欧洲航天局于2018年发射的ADM-Aeolus卫星，测量从地球表面到30 km高度的全球风速，用于提高天气预报精度。

如图12所示，波矢为 $ {\boldsymbol{k}}_i $ 的入射光被散射，散射角为 $ \theta $ ，则散射波数为 $ k = |{\boldsymbol{k}}_i-{\boldsymbol{k}}_s| = 2|{\boldsymbol{k}}_i|\sin\left({\theta}/{2}\right) $ 。散射光强度与光谱密度函数 $ S(f_s) $ 相关，即密度扰动相关函数 $ G({\boldsymbol{r}},t) $ 的时空傅里叶变换。假设散射波沿 $ x_2 $ 方向传播，相关函数和谱密度函数的解^[97-98]为：

$ G(x_2,t) = \int_{-\infty}^{\infty}{n}(l+x_2,t){n}(l,0){\rm{d}}l $

(111)

因此，频谱为：

$ S(f_s) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}G(x_2,t){\rm{e}}^{{\rm{i}}(kx_2-2\pi{f_s}t)}{\rm{d}}x_2{\rm{d}}t $

(112)

式中，i为虚数单位， $ {n}(x_2,t) $ 为密度扰动， $ f_s $ 为散射频率。

对于短波长和高散射频率，应该使用气体动理论来描述密度扰动的演化过程。以单原子气体为例说明瑞利-布里渊散射的计算。密度扰动量与式（20）中分布函数扰动量 $ \phi $ 有关，初始密度脉动可定义为^{[83, 97]}： $\phi(t = 0,{x_2},{\boldsymbol{\xi}})\propto{}\delta({x_2})$ 。取特征长度为散射波长 $\ell = 2{\text {π}}/k$ ，通过时间上的拉普拉斯变换与空间上的傅里叶变换求解线性化玻尔兹曼方程（21），可得：

$ 2{\text {π}}{{\rm{i}}}(f_s-v_2)\hat{\phi} = \frac{\ell}{v_m}J(\hat{\phi})+1 $

(113)

其中方程右端的1来自于初始密度扰动对应的拉普拉斯变换， $ \hat{\phi}({\boldsymbol{v}}) $ 分别对应 $ \phi $ 的时空坐标下的拉普拉斯和傅里叶变换^[17]。散射频谱可通过以下公式计算：

$ S(f_s) = \Re\left(\int \hat{\phi} f_{eq}({\boldsymbol{\xi}}){\rm{d}}{\boldsymbol{\xi}}\right) $

(114)

其中 $ \Re $ 表示复数的实部。

此例中，我们根据荷兰自由大学的实验数据^[91]来提取N₂的体积黏性和平动Eucken因子。其中，激光波长为 $ 366.8 $ nm，散射角为90°，因此等效波长为 $ \ell = 259.4 $ nm。实验条件下N₂的转动自由度与振动自由度分别为 $ d_r = 2 $ 和0。剪切黏性和热导率根据Sutherland公式计算。通过调整吴模型中非弹性碰撞数和平动热导率，使得实验数据和数值模拟得到的频谱偏差最小。表 2给出了不同温度下气体物性参数和提取的碰撞数与热导率分量。

实验数据与吴模型最佳预测的比较如图13所示。实验^[99]测得N₂在标准温度和压力下的转动弛豫时间为 $ \tau_r = 7.4\times10^{-10} $ s。在 $ 296.7 $ K下，测得的N₂体积黏性为 $\mu_b = 1.22\times10^{-5}\;{\rm{kg\cdot{m^{-1}}\cdot{s^{-1}}} }$ 。根据式（71），计算得到N₂的体积黏性为：

$ \mu_b = 2p\tau_r\frac{d_r}{(3+d_r)^2} = 1.18\times10^{-5}{\rm{kg\cdot{m^{-1}}\cdot{s^{-1}}}} $

这一结果与实验值吻合较好。另外，提取的平动Eucken因子也在合理的范围。

6 总结与展望

本文回顾了国内外学者在稀薄效应动理学建模中所取得的进展，包括单原子气体与分子气体。由于高温真实气体效应，当分子内自由度增加时，分子间复杂的相互作用通常伴随着更加新颖、复杂的物理现象，这对稀薄分子气体动理学建模提出了很高的要求。我们指出，应力偏量、能量交换、热流的弛豫过程是输运现象的本质，而剪切黏性、体积黏性、平动/内部热导率等输运系数则是表象，即，随着克努森数的变化，这些输运系数对应的本构关系会失效，而相应的弛豫过程和速率不会改变。为了准确预测稀薄分子气体的非平衡行为，气体动理学的建模过程应尽可能准确地捕捉应力偏量、能量交换、热流的弛豫过程和速率。

本文还介绍了DSMC中提取热弛豫速率的方法，并讨论了DSMC方法恢复热导率的能力。虽然经典的Larsen-Borgnakke模型可以通过调整非弹性碰撞概率改变体积黏性，但是很难同时调整热导率及其分量。而这些分量的比例以及本质的热流弛豫时间对稀薄气体流动的影响巨大。因此，非常有必要对DSMC碰撞模型进行修改和优化。

非常有意思的是，基于王承书等的工作，分子气体的线性化理论和方法已经非常成熟。而非线性模型虽然五花八门，但都与线性化模型存在着本质联系，即所有的非线性模型都可以通过对线性化Gross-Jackson模型和Hanson-Morse模型的非线性化而来。本文首次给出一个非线性化方案。

对于临近空间高超声速飞行器而言，其周围的空气动力学非常复杂。如何建立描述多组分气体、化学反应的稀薄气体动理学模型是个极大的挑战。本文提出的弛豫过程本质论、输运系数表象论、从线性化模型中反推非线性模型的方法，以及从分子动力学和相关实验中提取弛豫速率的设想，将为解决此类问题提供极大助力。

致谢：感谢爱丁堡大学苏微博士的有益讨论。