基于非线性耦合本构关系的改进边界条件

江中正，赵文文，袁震宇，陈伟芳

浙江大学 航空航天学院，杭州 310027

高超声速稀薄流动一直以来都是计算流体力学研究中的热点与难点问题。由于在远离局部平衡态流域，譬如激波间断处、靠近物面努森层内和尾流膨胀区，没有足够的分子碰撞去平衡气体输运行为，气体流动表现出极强的非平衡特点。采用一阶Maxwell-Smoluchowski (M-S) 滑移边界条件的线性牛顿黏性/傅里叶热传导本构模型(NSF)在这些流域中的模拟能力十分有限，无法准确地刻画非平衡流中所表现出的强非线性特点。为了准确模拟高超声速稀薄非平衡流动，国内外众多学者开展了大量研究[1]，提出了各种行之有效的理论和数值方法，包括粒子仿真方法——蒙特卡洛直接数值模拟(Direct Simulation Monte Carlo、DSMC)[2]、线性化玻尔兹曼方法[3]、模型方程方法(BGK、ES-BKG、Shakhov)[4]、离散速度法(Discrete Velocity Method, DVM)[5]、扩展动力学方法(Chapmann-Enskog)[6]、Burnett方程[7]、Grad矩方法[8]及其正则化理论[9]和全局双曲化理论[10]。近年兴起的统一气体动理论格式(Unified Gas-Kinetic Scheme, UGKS)[11]和气体运动论统一算法(Gas Kinetic Unified Algorithm, GKUA)[12]等气体动理学方法也为多尺度非平衡流动模拟提供了十分有效的手段。

为了弥补线性纳维-斯托克斯-傅里叶(Navier-Stokes-Fourier, NSF)本构模型的精度不足，同时克服粒子类方法和玻尔兹曼求解方法计算效率较低的缺陷，近年来一套非线性耦合本构关系模型(Nonlinear Coupled Constitutive Relations, NCCR)在由B.C. Eu提出的严格满足热力学第二定律的广义流体动力学方程(Generalized Hydrodynamic Equations, GHE)[13-14]的基础上被发展起来。该理论自提出之后得到了深入的研究和验证[15-19]，在非平衡流区域表现出优于Navier-Stokes方程的特点。然而，针对该理论的壁面边界条件研究为数不多，通常做法是直接采用一阶M-S边界条件开展计算[20-21]。在一阶M-S边界条件中，滑移速度或跳跃温度是速度或温度一阶梯度的线性比值，因此难以将NCCR模型的非线性特征表现出来，最终导致在壁面处人为降低了模型精度。

本文在Maxwell散射基础上开展了滑移边界理论研究，将努森层滑移跳跃效应通过与高阶非守恒量(应力和热流)表征的非平衡程度联系起来，提出一套与NCCR模型本身精度相适应的高阶边界条件，同时针对二维圆柱单原子气体绕流和平板绕流问题给出了不同稀薄程度下不同滑移边界条件NSF与NCCR模型的模拟结果，通过与DSMC的壁面结果进行比较，验证该边界条件模拟物面附近非平衡流动的能力。

1 非线性耦合本构理论

1.1 广义流体动力学方程

考虑外力影响的针对稀疏单组分单原子气体的玻尔兹曼方程表示为

(1)

式中：f表示包含时间、空间、速度七维相空间的速度分布函数；F为体积力；碰撞积分表示为

(2)

其中：g=|v-v1|表示粒子的相对速度；b为粒子碰撞参数，亦称为瞄准距离；φ为碰撞平面的方位角；v为分子运动速度。

为了推导广义流体动力学方程，在这里首先定义一个广义速度矩：

Φ(abc…l)=〈mCaCbCc…Clf〉

(3)

式中：角括号代表对分子速度的积分；C为分子热运动速度。对式(3)进行时间求导，得到广义矩的时间变化率，然后将玻尔兹曼方程式(1)代入，并定义碰撞耗散项为

Λ(Φ)(abc…kl)=〈mCaCbCc…ClC(f,f1)〉

(4)

便得到广义速度矩的演化方程为

(5)

若将关心的宏观物理量分为守恒量(ρ,u,E)和非守恒量(Π,Q)，均可以通过统一的微观式来表达，即

Φ(k)=〈h(k)f〉

(6)

式中：h(k)为两组物理量的分子微观表达式；宏观物理量(密度、动量、内能、应力张量和热流)分别表示为

(7)

其对应的分子微观表达式为

(8)

其中：m为分子质量；数学符号[A](2)表示二阶无迹对称张量，即

(9)

由于守恒量是碰撞不变量，不存在碰撞耗散项，因此由广义速度矩的演化方程式(5)在不考虑体积力影响下可以得到3大守恒方程为

(10)

(11)

(12)

为了推导非守恒量的演化方程以及封闭非守恒系统的分子碰撞耗散项式(4)，Eu[13, 22]构建了非平衡态分布函数的形态定义，在数学和物理的层面保证熵和Calortropy变化的非负属性。这个非平衡态分布函数充当了连接熵增源项和宏观非守恒量耗散演化过程之间的桥梁。该分布函数的单组分单原子气体分子的形态定义为

(13)

式中：μ为标准化因子，用以平衡气体分子数密度n与该模型分布函数的关系，如

其中：Xk为与宏观量有关的系数，其作用跟Grad矩方程中的系数相似；T、kB分别为温度、Boltzmann常数。利用Eu分布函数形态定义的形式计算熵增项，可以得到熵增项的表达式为

(14)

Eu认为能量耗散是系统的分子间碰撞导致的，微观上的分子碰撞又会导致非守恒量耗散演化的宏观过程，而熵增特性正是对能量耗散的一种直接体现。式(14)建立了耗散项、宏观非守恒量和熵增项的联系，这正是Eu修正矩方法的基石，也是和Grad矩方法有所区别的关键。通过对熵增项式(14)进行累积量展开[23]，并截取一阶近似，可得

(15)

当趋于平衡态流动时，κ趋于在近平衡态发展的Rayleigh-Onsager耗散函数，即

(16)

式中：d、p、η、λ分别为分子直径、压力、气体黏性系数和热传导系数。由式(16)可见Eu的熵增项累积量式(15)满足熵非负属性。在此基础上通过式(14)求得碰撞耗散项，最终可以得到一组针对单组分单原子气体分子的非守恒量演化方程为

(17)

式中：ψ(4)和ψ(5)为更高阶矩；cp为定压比热。式(17)结合3大守恒方程(式(10)～式(12))，即为Eu的广义流体动力学方程。

1.2 非线性耦合本构关系

Eu的非守恒量演化方程仍旧是开放的包含了更高阶矩的偏微分方程系统。为了封闭该方程组，Eu和Alghoul[24]提出了另一种不同于Grad封闭的基本封闭原则：

ψ(4)=ψ(5)=0

(18)

(19)

为方便对比，这里给出线性的牛顿黏性和傅里叶传热导定律：

(20)

2 边界条件

在高超声速稀薄流动模拟中，气固表面分子碰撞与作用机理是一个十分重要而且复杂的影响因素。入射气体分子在物体表面附近与物面反射分子相互碰撞聚集，这个聚集层通常被称为努森层，其厚度与平均分子自由程l是同一个量级。当来流并不稀薄，即来流平均分子自由程l→0时，努森层厚度趋于零，在微观分子层面表征为入射分子被物体表面吸附，在宏观层面即为最常用的附着壁和等温壁条件。然而，当来流稀薄程度增大，物面努森层效应逐渐显著，附着壁和等温壁边界条件不再成立时，气体流动与壁面之间会呈现不连续的现象，即引起速度滑移和温度跳跃现象。气固表面分子作用机理在准确刻画努森层流动乃至整个宏观稀薄流动中具有十分重要的作用。为此，有研究提出两种常见的气固分子碰撞模型：Maxwell散射理论和Langmuir表面吸附理论。作者前期研究工作[19]表明，尽管Langmuir吸附理论为气固碰撞提供了一种吸附等温物理模型，但由此衍生出的Langmuir滑移边界[28]并不能很好地预测稀薄流壁面处的滑移速度和跳跃温度。因此，本文主要基于Maxwell散射理论开展边界条件研究。该理论中Maxwell首先引入两个简单假设模型：镜面反射和完全漫反射。在来流反射分子中σ部分假定为完全漫反射，(1-σ)部分为镜面反射，σ为介于0～1之间的比例系数。

常见的一阶Maxwell滑移条件有3种形式：Maxwell滑移条件、Gokcen滑移条件、Lockerby滑移条件。Maxwell[29]在1879年基于分子运动论和小努森数假设，推导了一阶泰勒展开的滑移条件为

(21)

式中：u为靠近物面处气体切向速度；uw为物面速度；σu为速度适应系数；y为壁面法向；x为壁面切向；右端第2项则是在非等温壁下的热蠕动效应。

对于温度跳跃效应，Smoluchowski[30]在1898年也推出了相类似的温度边界条件：

(22)

式中：T为靠近物面处气体温度；Tw为壁温；σT为温度适应系数；γ为比热比；Pr为普朗特数。式(21)和式(22)也被称为M-S边界条件。M-S边界条件结合Navier-Stokes方程，既有效地增强了Navier-Stokes方程模拟滑移稀薄流的能力又保证了较小的计算量，在工程上得到广泛应用。

然而，Gökçen和MacCormack[31]研究认为基于小努森数一阶展开得到的M-S边界条件模拟努森数范围十分有限，在大的局部努森数区域预测的壁面结果与粒子仿真方法的结果有很大差异。因此，他们通过研究壁面切向动量损失和引入努森层处的速度，提出另一种通用的滑移条件:

(23)

式中:a代表速度或者温度；下标l表示一个平均分子自由程处的物理量，w表示物面处的物理量。该条件被认为在小努森数时能够还原Maxwell条件，在较大努森数条件下能得到更优的壁面结果。但是等式右端依然用到线性的NSF本构关系(应力张量随应变张量线性变化)，实际上这与努森层内流动速度分布的非线性特点不符。因此Gökçen边界条件的可用范围值得商榷。

此外，Lockerby[32]针对切应力均匀分布平板壁面的Cercignani线性玻尔兹曼速度解进行曲线拟合，通过引入一个壁面函数式(24)来修正努森层黏性ηnew=ηΨ-1的方式来修正线性的速度梯度项。

(24)

同时Lockerby等的研究[33]回归了不考虑热蠕动Maxwell边界条件的一般形式，认为右端项采用应力形式比应变率形式更具有普适性。

(25)

(26)

在Maxwell一般形式(式(25)和式(26))中，采用由壁面函数修正的有效黏性ηnew=ηΨ-1以及线性的牛顿黏性定律和傅里叶热传导定律，便可得到Lockerby边界条件。Lockerby通过进一步研究认为式(25)前面仍需要乘上系数0.8。该边界条件中通过构造非线性函数去修正线性的应力应变关系的做法，为尝试捕捉努森层的非线性物理量分布特点提供了一种新的思路。

此外，精确描述速度滑移和温度跳跃现象的另一条途径是构造二阶乃至更高阶的边界条件，Cercignani[34]、Beskok[35]、Deissler[36]、Hsia[37]等在Maxwell一阶边界条件基础上做了富有成效的研究。尤其Deissler通过物理推导手段以及Hsia和Domoto通过实验手段，均指出了单纯通过对Maxwell条件进行泰勒展开的数学方法来构造高阶边界条件会存在二阶乃至偶数阶滑移项反向过度修正的问题。为了区别传统的二阶M-S边界条件，现将偶数项修正之后的Hsia-Domoto边界条件表示为

(27)

(28)

对比Maxwell边界条件一般形式(式(25)和式(26))和不考虑热蠕动的常见形式(式(21)和式(22))，其主要的差异在线性的NSF本构模型(式(20))。参考Lockerby利用有效黏性修正真实气体黏性的非线性做法，对Maxwell边界一般形式采取另一种非线性修正方法：在保留真实气体黏性系数前提下结合NCCR本构模型，利用NCCR的应力和热流修正NSF的应力与热流项。考虑到前期研究工作[19]定量指出二阶修正项在准确预测努森层滑移边界值的重要性，在提出的修正边界条件也引入了二阶修正项(式(27)和式(28)) 的影响，并忽略热蠕动项，具体表示为

(29)

(30)

式(29)和式(30)中的应力和热流通过在壁面求解NCCR方程得到。此外，为了防止滑移速度和跳跃温度变化剧烈而导致高阶边界条件在迭代计算过程中的不稳定，Lockerby[38]和包福兵[39]在壁面处采用边界松弛方法进行计算，并取得较好的计算稳定性。因此，本文也采取相同的边界松弛方式:

(31)

式中：us和Tj分别为滑移速度和跳跃温度；上标r、n+1、n分别表示松弛之后、新时刻、上一时刻的值；Rf为迭代松弛因子，根据具体算例选取不同的值。

3 数值算法

选取无穷远来流参数对3大守恒方程和非线性耦合本构模型(NCCR)进行无量纲化，即

(32)

式中:下标∞表示无穷远来流参数。为了简化表达，下文均忽略上标星号，最终得到无量纲化的流动控制方程为

(33)

式中：守恒量U、无黏通量Fc和黏性通量Fv分别为

(34)

针对单原子气体的非线性耦合本构关系，有

(35)

式中：E为比总能；下标0代表线性的NSF应力和热流，上标^表示参数有如下关系：

(36)

其中：γ代表参考值。

3.1 空间离散

在有限体积法框架下离散求解流动控制方程式(33)。为了克服AUSM+的数值振荡和AUSMD的“红玉”现象，Kim等[40]结合这两种格式的优势并且考虑到计算效率和鲁棒性，通过引入修正的压力权函数，提出了AUSMPW+通量计算格式。该格式计算效率较高，鲁棒性好，而且激波捕捉能力强，在高速流动中得到了广泛应用，因此在本文的无黏通量计算中也采用AUSMPW+格式。AUSMPW+格式的通量表达式为

(37)

(38)

式中：

(Δ-)i=qi-qi-1

(Δ+)i=qi+1-qi

其中重构主要针对原始变量q进行，ε是防止分母为零的一个小量，而黏性通量均采用中心差分离散。

3.2 时间推进

本文采用LU-SGS (Lower-Upper Symmetric Gauss-Seidel)隐式时间推进方法，该方法具有较好的计算稳定性和收敛性，是最常用的隐式时间推进方法之一。对控制方程式(33)进行时间一阶隐式离散后可得

[Ι+Δt(DξA+DηB+DζC)]ΔQn=-ΔtRHS

(39)

式中：Dξ、Dη、Dζ为微分算子；RHS为残差项；A、B和C为无黏雅克比矩阵。通过矩阵分裂和LU分解得到LU-SGS方法的表达式为

LD-1UΔQ=-RHS

(40)

式中：

(41)

其中：A±、B±和C±为近似雅克比矩阵；ρ(A)、ρ(B)和ρ(C)为无黏雅克比矩阵的谱半径。式(40)可通过向前扫描、向后扫描和更新3个步骤进行求解。一般在式(39)中对黏性通量进行显式时间离散，但在黏性影响较大的区域容易造成计算不稳定。参考Navier-Stokes方程求解方法，本文也通过修改近似雅克比矩阵特征值，采用谱半径代替黏性雅克比矩阵，对NCCR黏性项进行了隐式处理，以ξ方向为例:

(42)

3.3 本构模型求解

NCCR本构模型式(35)是应力和热流的隐式表达，在3个方向上各项相互耦合影响，并不能像NSF本构关系简单地通过守恒量梯度一步显式求解，需要额外求解步骤。为此，Myong[17-18]提出的分裂算法通过忽略一些耦合项，将流动分解成3个方向之间互不影响的一维问题，然后进行线性叠加。前期研究[20-21]表明，分裂算法容易造成计算不稳定，尤其在尾流的模拟中表现明显，同时分裂算法忽略的一些耦合项也人为降低了模型的计算精度。因此，本文采用完整耦合的求解思路，利用最速下降法对NCCR本构模型进行求解，有效保证模型的计算精度。最速下降法不过分依赖初值的选取，是求解单组分单原子NCCR模型的有效算法。对于NCCR模型，可以表示成如下非线性方程组的一般形式：

F(x)=0

(43)

(44)

(45)

如何选取合适的步长αk，使得

(46)

是最速下降法的关键。最直接的解析方法是构造有关步长α的单值函数，即

(47)

通过对α求导解一元方程的根来寻找函数极小值。但是这种做法计算成本太高，本文采用一维搜索近似得到每步较优步长来进行迭代。选取3个步长α1=0<α2=α3/2<α3，使得h(α3)

(48)

4 计算结果与分析

4.1 圆柱绕流

首先研究一组在不同稀薄来流程度下的高超声速单原子气体圆柱绕流问题，该算例的DSMC数据以及Gokcen和Lockerby边界条件的NSF数据均来自文献[41]。圆柱直径为12英寸，如图1所示。无穷远来流速度为U∞=2 624.1 m/s，马赫数为10。采用等温壁模型，壁温为500 K。不同稀薄来流下的气体密度由表1给出，其中Kn∞表示采用无穷远来流条件计算的克努森数。

图1 几何外形定义Fig.1 Definition of geometry

Kn∞Mass density/(kg·m-3)Number density/(particles·m-3)0.0021.408×10-42.124×10210.012.818×10-54.247×10200.055.636×10-68.494×10190.251.127×10-61.699×1019

单原子气体为氩气，比热比为1.67，气体常数为208.16，普朗特数为0.67，其黏性计算采用变径硬球(Variable Hard Sphere, VHS)模型：

(49)

(50)

式中：黏性指数s为0.734；参考分子直径dr为3.595×10-10m；参考温度Tr为1 000 K；kB取1.380 658×10-23J·K-1；m为氩气的分子质量。

本文研究重点是量化比较线性(NSF)与非线性(NCCR)两个宏观本构模型在几组不同边界条件影响下预测的圆柱物面压力、摩阻和热流系数与DSMC结果之间的差异。因此物面无量纲压力系数Cp、摩阻系数Cf和热流系数Ch定义为

(51)

图2～图5分别给出了Kn=0.002,0.01,0.05,0.25 这4组条件下，压力系数Cp、摩阻系数Cf和热流系数Ch沿圆柱壁面随着几何角度θ的变化情况，并结合一阶M-S、Gokcen、Lockerby边界的NSF预测结果以及结合一阶M-S和二阶非线性修正边界的NCCR结果，与DSMC粒子仿真结果(由MONACO软件计算得到)进行比较。来流Kn=0.002的流动处于连续流域范畴。

图2 物面系数曲线(Kn=0.002)Fig.2 Curves of surface coefficient (Kn=0.002)

由图2可以看出，线性NSF模型和非线性NCCR模型的预测壁面压力、摩阻、热流系数结果与DSMC的结果非常吻合。考虑到计算成本的因素，线性NSF模型仍然是连续流域模拟中既有效又高效的理论方法。而从另一方面看，本文基于NCCR模型提出的二阶修正边界条件在连续流域能有效回归NSF的预测结果，体现了其在近平衡流附近的模拟能力。

图3 物面系数曲线(Kn=0.01)Fig.3 Curves of surface coefficient (Kn=0.01)

由图3可看出，当Kn=0.01时，除了在θ=120°尾迹区附近的物面，结合一阶M-S边界NSF预测的摩阻系数稍微偏高之外，其余位置不同方法的计算结果依然吻合较好。Gokcen和Lockerby边界比一阶M-S边界有效提升了NSF的摩阻预测能力。相比NSF本构模型，结合一阶M-S边界和二阶非线性修正边界的NCCR都准确模拟3组物面系数。来流Kn=0.01可以被认为是无滑移流域的上限，是连续流和稀薄滑移流域的界限。尽管流动在物面附近开始出现连续性假设失效，预测的结果表明滑移边界条件仍能够有效拓展宏观本构模型的非平衡流动模拟能力。

如图4所示，随着来流努森数继续增加，不同方法壁面压力系数的预测依然较为一致，但摩阻和热流系数预测的结果出现明显的差异。尽管来流Kn=0.05的流动处于滑移流域，但尾迹壁面附近区域的非平衡效应更加强烈，连续性失效更加严重。对于结合一阶M-S边界和Lockerby边界的NSF结果，尾迹区的物面摩阻和热流系数预测偏高，而Gokcen边界相对更为准确。与所有NSF结果相比，NCCR预测结果更加趋近于DSMC仿真结果，这表明非线性耦合本构模型(NCCR)比线性的牛顿黏性和傅里叶热传导模型(NSF)对非平衡流动模拟精度更高，能为有效解决连续性失效流域的模拟问题提供更加可靠的理论手段。

随着努森数进一步增加，流动离开滑移流域进入过渡流域。在来流Kn=0.25的流动中气体已经非常稀薄，由于缺乏足够的气体粒子碰撞，非平衡效应非常强烈。图5的物面系数分布表明，在过渡流域模拟中，滑移边界不再能有效地提升线性NSF本构模型的预测能力。强烈的非平衡效应不仅存在壁面努森层流动中，而且几乎扩展到整个流场范围。在过渡流域中，若远离物面流动仍然使用线性的NSF本构模型进行模拟，将严重影响到物面系数的预测精度。NSF的壁面压力和热流系数均比DSMC的预测结果偏高，除了采用Gokcen边界的NSF摩阻系数偏低外，其他滑移边界的摩阻系数也普遍偏高，壁面摩阻峰值偏高约12.5%。相比NSF的预测结果，NCCR在壁面系数的预测中表现更优，比如圆柱前部的压力系数和后部的热流系数更加趋近DSMC的结果。对于壁面摩阻系数的预测，一阶M-S和非线性修正边界的差异在过渡流域明显开始增大。随着气流加速流过圆柱，采用一阶M-S滑移边界的NCCR物面系数预测结果明显偏离DSMC，在θ=100°处(尾流膨胀区)偏离值达到最大，而结合非线性修正滑移边界的NCCR摩阻预测在迎风与背风区域与DSMC的结果基本吻合。可见，基于NCCR的非线性修正滑移边界在壁面处保留非线性优势，保证了NCCR理论模型的精度一致性，相比一阶M-S滑移边界极大提升了NCCR壁面系数的预测能力。

图4 物面系数曲线(Kn=0.05)Fig.4 Curves of surface coefficient (Kn=0.05)

图5 物面系数曲线(Kn=0.25)Fig.5 Curves of surface coefficient (Kn=0.25)

4.2 平板绕流

平板阻力预测是经典的流体力学问题，其外形简单，流动信息丰富，是组成复杂飞行器结构最基本的几何元素，因此在低速不可压、高速可压缩等领域得到广泛深入的研究[42-45]。对于高超声速平板绕流问题，激波与边界层干扰严重，黏性作用加剧，在平板前缘区域存在非常强烈的稀薄非平衡效应，这些叠加的复杂效应都会影响表面摩阻的准确预测，因此也是一个边界条件研究的经典验证算例。

本节以零厚度零迎角平板高超声速绕流进行数值模拟研究，进一步分析结合非线性修正边界条件的NCCR模型在复杂流动的预测能力。流域左边界和上边界是无穷远来流，来流温度为300 K，右边界为出口，下边界为X/L=0～1的零厚度平板，平板壁温为300 K，其余为对称边界。为了减少双原子气体在高速流动中转动自由度激发带来的影响，该算例依然采用单原子氩气，气体属性均采用4.1节的描述，重点研究了Ma=10、Kn=0.002和Ma=5、Kn=0.01两组平板绕流算例。壁面采用全漫反射条件。

图6和图7分别展示了两组来流条件下(Ma=10、Kn=0.002和Ma=5、Kn=0.01)由NSF本构模型结合一阶M-S滑移边界和NCCR模型结合非线性修正边界条件模拟的马赫数Ma、局部克努森数KnGLL、无量纲压力p/p∞和密度ρ/ρ∞云图情况。局部克努森数是由Boyd等[46]提出，用来具体刻画当地流域的连续性失效情况，即

(52)

该参数的方向导数采取当地梯度方向，Q可以根据具体的情况选取密度、温度或者压力。Boyd等[46]认为0.05可以作为当地连续性失效的判据。本文基于当地密度梯度来求局部克努森数。由图6和图7中两个算例的局部克努森数云图可以看出，整个流域除了平板前缘、激波和壁面附近外，失效区域和程度都不大，因此由NSF和NCCR模拟的流场大部分都相似，这也可以通过图中的马赫数、压力和密度云图的对比看出。当流动的连续性失效程度不大，稀薄非平衡效应不明显时，结合非线性修正滑移边界条件的NCCR能够回归传统的NSF解。从图6和图7中的马赫数云图还能够非常清晰地分辨在平板前缘之后的激波与边界层融合层以及之后紧邻的高超声速黏性干扰区，在Ma=5、Kn=0.01的算例还能观察到在尾部激波与边界层之间的无黏区，而在Ma=10、Kn=0.002的算例中，由于马赫数增大，激波角减少，边界层增大，其无黏区范围不明显。在融合层，激波尚未充分发展，而边界层相对较厚，此时激波与边界层相互融合，难以区分，这里也会存在比较显著的速度滑移和温度跳跃的稀薄现象；而之后的高超声速黏性干扰区，激波已经充分发展，尤其在Ma=10时非常靠近边界层外缘，边界层则迅速增长，黏性效应增强，与激波相互影响着波后的无黏区和平板表面的热流和摩阻分布。

图6 马赫数、局部克努森数、无量纲压力和密度云图对比(Ma=10、Kn=0.002)Fig.6 Comparison of detailed contour of Mach number, local Knudsen number, dimensionless pressure and density (Ma=10，Kn=0.002)

图7 马赫数、局部克努森数、无量纲压力和密度云图对比(Ma=5、Kn=0.01)Fig.7 Comparison of detailed contour of Mach number, local Knudsen number, dimensionless pressure and density (Ma=5,Kn=0.01)

图8 表面摩擦阻力系数曲线(Ma=10、Kn=0.002)Fig.8 Curves of surface friction coefficient (Ma=10，Kn=0.002)

图9 表面摩擦阻力系数曲线(Ma=5、Kn=0.01)Fig.9 Curves of surface friction coefficient (Ma=5,Kn=0.01)

准确预测超声速流动的摩阻一直是气动研究的重难点。为了定量分析基于NCCR的非线性修正边界的能力，本文具体比较了由DSMC、结合一阶M-S的NSF、结合一阶M-S和非线性修正边界的NCCR预测平板单侧表面摩阻系数分布的情况，如图8和图9所示。DSMC的结果采自文献[44]。在Ma=10、Kn=0.002的X/L≥0.05和Ma=5、Kn=0.01的X/L≥0.2平板范围，NSF、NCCR和DSMC三者预测的摩阻非常吻合，说明这些区域仍能满足连续性介质假设，稀薄效应不明显。而在这些范围之前平板前缘和融合层区域，三者预测的结果产生显著的差异，NSF即使结合一阶M-S边界也无法准确预测该处的摩阻，过高预估了该值。相比结合一阶M-S的NCCR解，NCCR结合非线性修正滑移边界条件对壁面摩阻分布和极值大小的预估精度有着比较明显的提升，更加吻合DSMC的结果。从图6和图7的局部克努森数云图可以发现，该处平板能够达到KnGLL=0.5以上，处于过渡流域，稀薄非平衡效应显著，因此结合滑移边界的NSF能力明显不足，无法胜任对于该处模拟的需求，需要精度更高的理论模型结合高精度的边界条件进行预测。另外，在该处的摩阻分布，还能观察到一个有趣的现象，表面摩擦阻力系数有个先增后降的过程，即在平板某处产生一个摩阻峰值。这个现象实际上在一定程度反映了平板前缘的非平衡效应。当地克努森数较大，壁面速度滑移明显，因此会引起壁面切向速度的法向梯度的减少，从而降低了表面摩阻系数，这就是造成在X/L=0处摩阻系数非常小的原因。

5 结 论

本文针对最新发展的非平衡气体动力学理论——非线性耦合本构关系模型(NCCR)开展了边界条件研究，提出了一套在壁面处与NCCR模型精度保持一致的非线性修正滑移边界条件。在有限体积的计算框架下，结合LU-SGS隐式时间技术和AUSMPW+现代通量计算格式，采用耦合求解NCCR模型的思路，对不同稀薄程度的高超声速单原子气体圆柱绕流和平板绕流问题进行数值模拟，重点比较了采用不同滑移边界的NSF和NCCR模型对物面压力、摩阻与热流系数的预测情况。

1) 结合滑移边界条件的NSF只能有效拓展到滑移流域的模拟，NCCR模型可准确模拟的非平衡流域范围比NSF更宽。

2) 结合一阶M-S边界的NCCR模型在非平衡效应强烈的努森层流动以及激波与边界层的平板融合干扰区的模拟中会弱化NCCR的非线性优势，摩阻预测明显偏离DSMC仿真结果。而本文的非线性修正边界条件有效增强了NCCR模型在过渡流域和激波/边界层融合干扰区的非线性模拟能力，在壁面流动中保证了与流动模型精度的一致性。

3) 在稀薄效应影响下，NSF预测的摩阻和热流结果偏高，而NCCR模拟结果更趋近于DSMC粒子仿真结果。