0 引言

我国的空气动力学自20世纪50年代中在钱学森先生回国后的全面规划建议下，至今已得到极大的发展^[1]。其中计算流体力学(Computational Fluid Dynamics, CFD)是非常突出的一个方面^[2]，以张涵信先生的团队为代表。他们的一个鲜明特点是无论是在发展算法，还是在分析数值结果时，都以清晰的物理分析为基础。一方面成功构建了多种有效的算法，从而在实际的复杂流动计算中提供了高保真数据^[3]。另一方面又能给出清晰的流动物理图像，为航空航天技术的发展做出了巨大贡献。

但在近年来的近空间飞行器的研发中发现，原有的空气动力学及相应的CFD不能完全应对遇到的情况，需要研究新的空气动力学问题^[1]。

近空间高超声速飞行器飞行高度一般在40-70 km的高空, 飞行速度一般在5-20倍声速的范围^[1]，其研发在中国及世界范围内已受到越来越多的关注。在此类飞行器的飞行过程中，流场局部的高温和低的压力等因素会导致局部区域内存在气体稀薄效应，使得连续介质模型失效，传统的CFD方法对气动力，尤其是摩擦阻力及热流的预测可能会有较大误差。

研究涉及稀薄气体效应的流动，首要的问题就是用什么参数判别稀薄效应的程度。钱学森先生最先根据稀薄程度将稀薄气体流分为三大流域^[4]。根据Knudsen数的大小，其定义为平均分子自由程λ与飞行器某一特征长度L的比值，三个流域依次为滑移流(0.01 < Kn < 0.1)、过渡流(0.1 < Kn < 10)和自由分子流(Kn > 10)。当时关注的主要是飞行器的压力分布，这样的划分基本上是可以的。后来，对于简单管道微尺度流动，L通常取管道的水力直径，当Kn数小于0.2时，连续介质模型结合二阶滑移边界条件可以获得与实验一致的结果^[5]。

但对于高速飞行器周围的流场，情况极为复杂。用一个参数Kn不足以刻画其流动状态。针对不同问题，要选取不同的物理特征长度。例如，在飞行器头部附近，温度和压力变化迅速，用飞行器整体的特征长度和来流的气体平均分子自由程不能反映当地的情况，即原先定义的Kn数并不适用，需要考虑局部的稀薄效应。因此，最直接的方法是使用局部的平均自由程和局部流动的特征尺度，如激波脱体距离，或一些宏观量梯度的标尺长度^[6]，建立局部的Kn数。最早Bird针对膨胀流的连续介质失效判定提出的P参数^[7]，它正比于Maλ(dρ/dx)/ρ。Boyd等认为在高超声速稀薄流动中应同时考虑局部的密度和温度梯度^[8]，提出了Kn_GLL=λ/(Q/|dQ/dl|)，Q为密度或温度，导数应沿梯度最大方向l。数值计算表明对于一维正激波和二维圆柱绕流的连续失效阀值为Kn_GLL-D=0.05。随后，Wang等将其写为Kn_Q=λ/(Q/|∇Q|)^[9]。Kn_GLL、Kn_Q是目前在NS-DSMC耦合方法^[10-11]中常用的对计算区域进行划分的参数之一。在一些动理论方法中，连续失效准则可以由某一类物理量(如剪应力、热流、速度分布函数等)与连续方法计算的值之间的相对误差的形式建立。Lockerby等提出了Kn_L=max(Kn_σ, Kn_q)^[12]，其中Kn_σ、Kn_q分别定义为R-13矩方法^[13]所求的剪应力、热流与在R-13流场结果上依据线性本构关系所求的剪应力、热流之间的相对误差。随后，Meng等直接从动理论的分布函数出发，定义了基于分布函数的失效准则E_c^NSF[14]，主要用于多尺度动理论方法^[15]中低阶模型与高阶模型的选择。另外，热力学不可逆过程熵增条件在非平衡效应判断和动理学模型构造中引起了关注^[16]。Camberos等基于只要是非平衡状态都有熵产生，将熵产率作为可以普遍预测非平衡的重要参数之一^[17]，并作为连续失效的参数。Xiao等基于从非平衡态到平衡态熵增函数数学模型收敛为Rayleigh-Onsager函数提出了新的近平衡失效参数$ \hat R$^[18]。但现有参数对于近空间飞行器的壁面摩擦阻力和热流问题并不适用，其壁面附近流动的特点为高温、低密度、强剪切流动。

在稀薄气体数值模拟方面，Bird提出的直接模拟Monte Carlo方法(Direct Simulation Monte Carlo，DSMC)至今被认为是在模拟稀薄气体流动，尤其是在过渡流区域最成功的方法之一^[19-21]。该方法在涉及高速、高温、化学反应等的稀薄气体流动实际应用方面有很多成功例子^[22]。为了解决DSMC在低速流动中的统计涨落问题，我国沈青、樊菁发展了信息保存法(Information Preservation, IP)^[23]。另一方面，被认为可以描述整个流域的气体动理论基本方程——Boltzmann方程求解十分困难^[24]。最近我国学者在Boltzmann方程简化模型的数值求解方法，特别是跨流域统一格式方面，取得了重要进展^[25-28]，并引起了国外同行的普遍关注。但是对于三维近连续域气体动力学问题，其计算量相比传统CFD依旧很大，现阶段计算资源尚难满足近空间高超声速飞行器的工程计算。此外，大部分动理学格式建立在Boltzmann-BGK简化方程之上，此简化模型的理论基础为气体偏离平衡状态不太远。如本文下面所述，这一假定在实际的近空间飞行条件下与实际情况相去可能很远。另外一些学者从连续性方法出发，在小Kn数下从流动守恒和本构关系方面发展了一些基于连续模型的方法。Chapman和Enskog将分布函数在Maxwell平衡态附近展开为正比于Kn数幂次的级数^[29]，Burnnett以此为基础建立了Burnett方程等^[30]。Grad用Hermite多项式展开分布函数^[13]，得到了13个Boltzmann方程的矩方程，即Grad十三矩方程, 用于过渡流模拟^[31]。Eu从不可逆热力学过程出发，对Boltzmann方程模型化，建立了非线性耦合本构关系(Nonlinear Coupled Constitutive Relations，NCCR)理论^{[16, 32]}，提供了一种新的可能，该方法在真实三维稀薄流中的适用性仍有待验证^[33]。

综上分析，存在局部稀薄效应的近连续流的求解是目前最困难的问题之一，该问题正好对应于近空间高超声速飞行器的研发。传统连续介质模型中的线性牛顿黏性定律和傅里叶导热定律是否能给出足够精确的结果需要重新检验。而气体动理论方法由于其巨大的计算量难以应用于实际复杂外形。所以研究局部稀薄效应现象及其判别参数、认识传统的CFD适用范围及问题所在、考虑如何将局部稀薄效应包含于CFD计算中有非常重要的意义。

本文将选取具有代表性的高速剪切流为研究对象，对局部非平衡状态以及对剪应力计算等问题进行详细分析，进而找到一个可以反映剪切流局部非平衡效应的无量纲参数，并揭示此参数和剪切力与传统连续介质模型下剪切力偏移量的规律。

1 模型描述及计算方法 1.1 问题描述

周恒和张涵信在讨论与近空间飞行器对应的空气动力学面临的新问题^[1]时指出“气体分子运动论用于空气动力学时，仍然认为分子热运动速度的分布对Maxwell分布只有小的偏离。而这在流速很高的情况时是很可疑的。文献[1]中采用四种不同的数值方法模拟了长1 m，在70 km高空以马赫数15零迎角飞行的零厚度平板的绕流问题(见文献[1]中图 1的结果)，如暂时以图示计算结果为参考来估算，可知约在平板法向4 cm处(见文献[1]图 1(h))，气体密度最低，约为0.4，因而该处分子自由程约为5 mm(来流分子自由程为2 mm)。同时该处温度约为10(见文献[1]图 1(g))，对应于2000K(来流温度为200K)。该处速度梯度约为32 m·s^-1每毫米(见文献[1]图 1(e))，或在法向相差一个分子自由程的两层间速度差约为160 m·s^-1。如果假设分子热运动平均速度与温度的二分之一次方成正比，则该处的分子热运动平均速度约为1190 m·s^-1(在室温300K下，分子热运动平均速度约为450 m·s^-1)。粗略地讲，相隔一个分子自由程的两层间的分子碰撞是产生剪切力和热传导的主要机制。而在上述这个速度很高的例子下，两层分子的碰撞要使各自的平均速度改变160 m·s^-1，相对于其热运动平均速度，显然不是一个小量，它会使各自的分子热运动的速度分布显著偏离标准的Maxwell分布。”^[1]

从上面所引的周恒和张涵信的分析中可以看出，他们认为在有气体稀薄效应而剪切又很强的情况下，分子热运动速度的分布对Maxwell分布只有小的偏离这一假设是很可疑的。而起决定性作用的参数是相隔为一个分子自由程的两层流体的速度差和分子自由运动速度之比这一无量纲参数，即

$ \lambda \left| {\frac{{{\rm{d}}u}}{{{\rm{d}}y}}} \right|/{V_{{\rm{mp}}}} $

其中，λ是分子自由程，$\left| {\frac{{{\rm{d}}u}}{{{\rm{d}}y}}} \right| $是当地宏观速度u在其垂直方向的梯度，V_mp是当地分子最可几热运动速度。以下用Zh表示这一参数。他们认为，当这个参数足够大时，气体分子自由运动速度的分布函数会显著偏离Maxwell分布。

1.2 Couette流

本文重点研究剪切力受气体稀薄效应的影响，因此选取最简单的剪切流，平面Couette流为分析模型，即由两个相对运动的无限大平行平板表面剪切力驱动的流动。针对中等稀薄程度(即传统的Kn数为0.005~0.5之间)，但存在强剪切的情况进行详细的研究。图 1描述了Couette流示意图，壁面温度为T_w= 273.15 K。以氩气作为模拟气体，初始温度为273.15K，此时黏性系数为μ =2.12×10^-7kg·m^-1·s^-1。当压力为标准大气压p₀，初始温度为T₀=T_w时，分子平均自由程为λ₀=62.50 nm，相应的分子温度的最可几速度为V_mp=337 m·s^-1。在计算中，壁面速度U分别取V_mp、3V_mp、5V_mp、8V_mp，两板间距为L=200λ₀，模拟中通过改变初始分子数密度而改变气体的稀薄程度。

1.3 数值模拟方法

本文采用基于Bird提出的标准DSMC算法建立的自编程程序。在计算中，碰撞对的选取采用NTC(No Time Counter)方法，碰撞模型采用硬球模型，即黏性系数与温度的根号成正比。在壁面处采用完全漫反射条件。在所有DSMC计算中，计算网格均小于当地的平均自由程的三分之一、时间步长小于平均碰撞时间的五分之一，使用了大量的模拟分子，初始每个网格内分子数为100个，整体样本统计量在10⁸的量级，以确保得到较好的粒子速度分布函数。同时，我们验证了网格尺寸对粒子速度分布函数的影响，网格尺寸在满足基本条件基础上继续减小，对分布函数产生的影响非常小。

在本文2.3节中，我们用传统的CFD模拟了高超声速小钝度钝楔模型的绕流。其飞行参数对应于近空间飞行器的参数范围。目的是考察实际飞行条件下，是否真的需要考虑气体的局部稀薄效应。模拟采用无滑移等温壁面，远场激波外为来流条件，出口为外推边界条件。空间离散分别采用五阶WENO格式和六阶中心格式。因是定常流计算，时间方向采用2阶Runge-Kutta方法。

2 结果分析

本文针对高速剪切流主要研究两个问题：1)在弱稀薄但强剪切的情况下，气体偏离平衡状态的程度？即气体分子热运动速度分布相对于Maxwell分布的偏离程度。如的确处于不平衡状态，参数Zh是否是表征非平衡程度的一个合适的参数？对这一问题，基于DSMC方法的微观特性，统计获得粒子的概率密度分布函数与Maxwell平衡态比较是最直观的方法。2.1节将给出这方面的结果。2)在判断了气体是否处于平衡状态之后，首要关心的问题就是对实际剪应力的影响。基于连续介质的牛顿黏性假设在什么范围可以给出正确的剪应力计算？如果不正确，偏差如何？2.2节将详细讨论这一问题。

2.1 概率密度分布函数

我们计算了不同密度、不同壁面速度下30个算例，画出了中心点的概率密度分布函数，计算了当地的Zh值，我们发现，可以通过Zh的值将分布函数的形状分为三种。

图 2、3、4分别给出了三种典型的概率密度分布函数以及相应的当地Zh参数的数值。中心点的宏观速度接近于0，粒子速度分布以速度v=0对称。在稀薄程度较小、剪切强度较小时，速度分布函数与Maxwell平衡态给出的一致，如图 2(a)所示，此时中心点的Zh为0.02。随着气体稀薄程度、剪切强度的增加，速度分布函数逐渐偏离Maxwell分布且偏离的程度越来越大。以中心点当地的Zh值可以将分布函数的状态大致分为三种如图 2所示，即单峰的近似正态分布(Zh < 1.0)，三峰分布(1.0 < Zh < 1.5)及双峰分布(Zh>1.5)。

在与周恒院士的讨论中，他用一个更简单的物理模型定性地说明了上述现象。他以相距各为一个分子自由程，平均速度差为u的三层流体作为分析对象(如图所示5)，分析中间一层流体的分子运动速度沿流向分量的分布函数的可能变化如下：

设初始时各层流体分子运动速度分布函数都是Maxwell分布，其中心速度就是该层流体的平均速度。当上层流体分子来到中间层，与该层分子碰撞后，其总的效果就是使得速度大于平均速度的分子数增加，即原先Maxwell分布的正向某处(略小于u值)附近的分子数增加，相应地原先对应于峰值处的分子数减小。而下层流体的分子的作用类似，即使得速度小于平均速度的分子数增加，使得原先Maxwell分布的负向某处(略大于-u值)附近的分子数增加，相应地原先对应于峰值处的分子数减小。中间层的分子也可以相互碰撞，但并不会将上述的正负效应相抵消，因为多出来的分子运动速度大于平均速度的分子正好和相应的分子运动速度小于平均速度的分子相碰的概率并不大。因此，总的效果是原先Maxwell分布的峰值处会降低，而两侧则会增高。

当Zh小于1，对应于u值较小的情况，因此其总的效果只是原先的Maxwell分布峰的两侧近处略微增厚，峰略微下降，因而总体仍可保持一个对称的单峰分布如图 2所示。当Zh值大到一定程度，u的值增加，原先Maxwell分布的峰的两侧较远处的分子数增加，就会出现两个峰而不仅是原先峰的两侧近处的增厚。但这时原先的峰还没有降得足够多，因此速度分布函数就会出现三个峰的情况，如图 3所示。当Zh继续增大，则上下层分子的影响下形成的两侧的峰相隔更远，且粒子速度偏向u和-u的分子数随Zh增大而增加，使得主峰降得更多，从而出现两个峰情况。此时粒子速度分布函数已经远远偏离Maxwell分布了。

因此，通过上述的数值模拟的结果和分析，从物理机理出发提出的参数Zh可以很好地表征剪切流中的局部非平衡效应。在Zh比较大的情况下，粒子速度分布函数已经严重偏离Maxwell分布。而目前被广泛采用的Boltzmann-BGK方程建立的依据是粒子分布函数偏离Maxwell分布不太远，并用简单的向Maxwell分布趋近的弛豫过程代替分子间的碰撞作用。因此，虽然基于Boltzmann-BGK方程的各种动理学方法简化了Boltzmann方程的计算，有效降低了计算成本，但用以模拟高速稀薄剪切流时则未必能给出正确结果。

2.2 剪应力

从连续介质观点出发，由于黏性作用，一层流体对相对运动的另一层流体产生阻力，这种阻力称为黏性切应力，并且建立了黏性切应力和切变率之间的线性关系，即牛顿黏性假设。而从微观角度看，气体剪应力是分子动量交换的宏观表现，快速层中的动量较大的分子跃入慢层，给慢层增加了动量，慢速层中的动量较小的分子跃入快速层，减少了快速层中的动量。在DSMC方法中，我们可以通过统计每一个网格上进出分子的动量改变量，因而获得相应的黏性剪应力。另一方面，利用在DSMC计算中所得的宏观速度剖面，可以依据牛顿黏性公式计算出另一个剪应力。我们预期，当气体稀薄效应较低时，两者应基本相符，而随着气体稀薄效应增强，二者的差别将扩大。

以DSMC统计所得的剪应力为准，定义二者的相对误差为Δτ=(μdu/dx-τ_DSMC)/τ_DSMC，则Δτ和Zh的关系如图 6所示。可以看出，由牛顿黏性假设计算的误差随Zh的增加而单调增加。图中不同类型的标记对应于不同的壁面速度。从图中可以看出，当Zh小于0.2时，偏差在10%以内，随着Zh的增加偏差急剧增大，当Zh约为0.67时，误差已达到100%。

值得注意的是，由于内摩擦作用，流场中的温度并不均匀，中心点温度最高，沿壁面法向成抛物线分布。图 6中不同Zh情况下中心点的温度并不相同，也不是随Zh单调增高或降低。例如，对应于Zh=0.05的点，它的温度是7.6T₀，但剪应力偏差仅为3.8%。而对应于Zh=0.51的点，其温度为1.5T₀，但剪应力的偏差为48%。所以图 6中所得的结果对温度不具有依赖性。

也可以把图 6中的结果，表达成两种剪应力之比A_e=τ_DSMC/(μdu/dx)和Zh的关系，如图 7所示。

由以上结果我们分析可知，牛顿黏性公式在Zh较大的时候不再成立，而且所给出的剪应力之值总是大于由DSMC给出之值。如果可以认为DSMC所给之值更可信，则在CFD中，可将剪应力写成如下形式：

$ \tau = {\mu _e}\frac{{{\rm{d}}u}}{{{\rm{d}}x}} = \mu {A_e}\frac{{{\rm{d}}u}}{{{\rm{d}}x}} $

其中μ_e为等效黏性系数，μ为气体的原本黏性系数，此黏性系数仍与温度有关。A_e为黏性修正系数，其值随Zh的变化已示于图 7中。

2.3 高超声速小钝度钝楔模型绕流

为了检验在飞行器飞行的实际流场中是否存在Zh比较大的情况。我们采用CFD模拟了高超声速小钝头楔的绕流问题。在此算例中，楔头半径R₀=1 mm、半楔角5°，钝楔总长3 m。来流马赫数为15，迎角为15°。气体性质取70 km高空空气之值，黏性系数与温度的变化采用苏士兰(Sutherland)公式。壁面温度固定为1000K，采用无滑移边界条件。

图 8给出了整个流场的Zh分布的云图。由图可以看出，在背风面的绝大部分地方Zh值都超过0.2。在流场中沿壁面法向搜索Zh的最大值，结果发现都是在壁面处。图 9给出了Zh沿壁面流向的分布，可以看出，背风面的Zh值明显大于迎风面的值。这是因为背风面的压力明显低于迎风面的压力，相应地背风面的分子自由程要显著大于迎风面的分子自由程。无论是背风面还是迎风面，Zh的最大值都出现在楔的头部附近。

Zh值在前缘处迅速减小，主要是由于在前缘处边界层厚度急速增加，速度梯度急剧减小。在背风面，Zh的值从x=133 mm时的0.5缓慢地降为x=3000 mm时的0.24。在这一段中，本文的结果应该都有效。但在更上游部分，Zh沿流向的变化很快，本文的在近似于宏观定常的情况下所得结果，未必可以直接应用，需要另做研究。

3 讨论

本文采用DSMC方法分析了高速剪切流中局部的非平衡效应，提出了一个新的描述参数Zh，即一个平均自由程上的速度差与当地的分子运动最可几速度的比值。通过数值分析我们得到：

1) Zh参数能有效地刻画气体速度分布函数偏离Maxwell的程度，以Zh的值可以将分布归为三类：单峰的正态分布(Zh < 1.0)，三峰(1.0 < Zh < 1.5)及双峰分布(Zh>1.5)。

2) 由基于连续介质的牛顿黏性假设所计算的剪应力相对于由DSMC所得之值的差随着Zh单调增加。

3) 气体在一定稀薄程度、强剪切情况下呈现非牛顿现象，本文提出一个等效的黏性系数，并给出了修正系数。

4) 因为2.3节中的情况，其参数在实际的近空间飞行器的范围之内，因此本文的结果对设计飞行器是有实际意义的。

但本文终究只是研究了一种情况，即气体稀薄效应对剪切力的影响，而且是在流动演化是非常缓慢的条件下。不难想象，要解决真正的实际问题，还有不少问题需要研究。例如，气体稀薄效应对热流的影响，这时参数Zh是否仍是唯一需要考虑的参数?反过来，温度梯度如果很大，是否也会对剪切力和热流产生影响?决定其影响是否需要另一个无量纲参数？本文中的结论虽然和温度值无关，但如果温度或温度梯度很高，这个结论是否仍然正确？对应2.3节中接近于楔的头部处，Zh值很大，而且沿流向的变化非常快，本文在流动变化缓慢条件下所得结论显然不一定能直接应用，而必须考虑快速变化或非定常演化的因素，等等。我们将在今后的工作中陆续针对上述诸问题继续工作。希望最终可以得到能解决在高超声速飞行器设计中遇到的空气动力学新问题的结果。

致谢:

感谢周恒院士给予的直接指导帮助，以及张涵信院士在和周恒院士讨论中给出的间接帮助和鼓励。感谢天津大学青年人才自主科研基金支持。