尺度自适应的离散统一气体动理学格式及在可压缩流动中的应用
2020-05-20许丁,孙祥,刘欣
许 丁, 孙 祥, 刘 欣
(1. 西安交通大学 航天航空学院, 机械结构强度与振动国家重点实验室, 西安 710049; 2. 陕西省先进飞行器服役环境与控制重点实验室, 西安 710049)
0 引 言
近来在计算流体力学(Computational Fluid Dynamics, CFD)中，数值模拟宽范围努森数Kn的多尺度流动问题成为一个热点。随着Kn的变化，流动对应领域及流动特征尺度不尽相同，对传统CFD方法提出了挑战。通常根据Kn对气体流动领域进行了如下划分[1]：在Kn<0.001时属于连续流动领域，传统基于连续介质假设的N-S方程和黏性壁面处的无滑移边界条件是适用的，流动对应宏观动力学特征尺度，如绕流物体特征长度，这正是传统CFD方法蓬勃发展并得到广泛应用的领域；当0.001
尽管对于Kn较大的过渡流及自由分子流，Bird[3]提出的Direct Simulation Monte Carlo(DSMC)方法是一种有效方法，但是对于连续流动，DSMC所消耗的计算时间及计算资源将异常巨大。因此发展适用于宽范围Kn数变化、跨流域的统一算法对解决实际工程当中的多尺度流动是一个关键问题。
另一方面，直接求解分布函数满足的Boltzmann方程是气体动理学方法的研究思路，近来一些气体动理学格式也应运而生。如格子Boltzmann方法(Lattice Boltzmann Method, LBM)[4-6]、气体动理学统一算法(Gas Kinetic Unified Algorithms, GKUA)[7-11]、气体动理学格式(Gas Kinetic Scheme, GKS)[12-16]、统一气体动理学格式(Unified Gas Kinetic Scheme, UGKS)[17-22]、离散统一气体动理学格式(Discrete Unified Gas Kinetic Scheme, DUGKS)[23-25]等。这些方法与传统基于N-S方程的CFD方法相比，可以获得更多的流动信息，且已有大量成功的应用[26-33]。
本文从DUGKS格式出发，指出DUGKS在处理碰撞项时，采用显式和隐式的简单算术平均并未充分考虑不同流动领域碰撞项对应物理过程的不同，尤其没有充分考虑两个时间尺度即计算时间步长和碰撞松弛时间之间的相对大小对粒子碰撞过程宏观表现的影响。为此，本文将对碰撞项的离散采用显式和隐式的加权平均，且权系数紧密依赖于计算时间步长和当地碰撞松弛时间之比，即当地努森数。进一步通过物理上的分析与数学上的推演，表明该权系数可以根据当地流动特征尺度进行自适应调节，从而将原始DUGKS格式进一步发展成为具有尺度自适应特性的离散统一气体动理学格式(Scale Adaptive Discrete Unified Gas Kinetic Scheme, SADUGKS)。
1 Boltzmann-Shakhov模型方程及其沿特征线离散的一般形式
为了克服Boltzmann-BGK模型对应普朗特数Pr≡1与实际气体不相符，本文基于Boltzmann- Shakhov模型来构造计算格式。
在D维空间中，Boltzmann-Shakhov模型[24, 34]如下：
(1)
其中，f=f(xi,t,uj,ηk,ξl)为D维空间中t时刻、位置x=(x1,...,xD)处、速度为u=(u1,…,uD)的粒子分布函数，η=(η1,...,η3-D)为宏观速度恒为零的方向上的粒子运动速度，ξ=(ξ1,…,ξK)为粒子内自由度，K为粒子内自由度的维度，与粒子种类相关；τ为碰撞松弛时间，与动力学黏性系数μ及压力p相关，τ=μ/p；fs为Shakhov平衡态分布函数，可表示为Maxwell平衡态分布函数及Pr数修正两部分。
为节约内存开销和计算量，将f对内自由度和宏观速度恒为零的速度方向上积分得到两个约化分布函数：
(2)
(3)
宏观守恒量Q=(ρ,ρU,ρE)T、应力张量τ及热流q可由分布函数求矩得到：
(4)
(5)
(6)
其中c=u-U，geq为Maxwell平衡态分布：
(7)
对式(1)积分可得到约化分布函数g及h的控制方程：
(8)
(9)
其中，gs和hs分别为g和h的Shakhov平衡态分布函数，二者均可表示为Maxwell平衡态及Pr修正两部分：
(10)
(11)
(12)
heq=(K+3-D)RTgeq
(13)
(14)
根据粒子碰撞过程满足质量、动量与能量守恒定律，约化分布函数g及h满足相应相容关系：
(15)
式(8)与式(9)形式上一致，可统一表示为：
(16)
其中φ=g或h。对于式(16)，t+s时刻的分布函数可形式上表示为[12]：
φ(xi,t+s,uj)=
e-s/τφini(xi-uis,uj)
(17)
其中：t为初始时刻，s为时间发展长度，x′i=xi-ui(s-t′)为粒子运动轨迹；φini为初始t时刻的气体分布函数，即φini(xi,uj)=φ(xi,t,uj)。
尽管Boltzmann-Shakhov方程(16)基于微观观点建立，其所蕴含的流动尺度却不限于微观动理学尺度，这一点可以从式(17)清楚看出：式(17)右边第一项积分项反映粒子之间相互大量碰撞作用的积累效果，该项已被证明可以看作为宏观连续流动动力学模型[12, 35]，对应尺度为宏观动力学尺度，如绕流物体特征长度；式(17)右边第二项反映粒子的迁移运动，对应尺度为微观动理学尺度，即分子平均自由程和平均碰撞松弛时间。式(17)将宏观动力学尺度和微观动理学尺度有机地结合到一起，可以很好反映多尺度流动特性。但是式(17)并不能直接使用，因为右边第一项积分项所涉及的φs是未知的，需要满足相容关系式(15)。因此从这个角度上来说式(17)只是式(16)形式上的解。
另一方面，认识到Boltzmann-Shakhov方程(16)自身已将宏观动力学尺度和微观动理学尺度有机地结合到一起，因此可以直接从Boltzmann-Shakhov方程(16)出发来构造多尺度计算格式。这里将式(16)沿特征线进行离散，可得到：
(18)
其中：
φs,ini(xi,uj)=φs(xi,t,uj)
(19)
分别代表初始时刻t的平衡态分布函数和碰撞项。式(18)右端的碰撞项采用了显式和隐式加权平均处理。α为隐式部分对应的权系数，其为碰撞松弛时间τ和时间发展长度s的比值τ/s的函数，具体形式待定。式(18)可进一步整理为：
φ(xi,t+s,uj)=ϑ1φs(xi,t+s,uj)+
ϑ2φs,ini(xi-uis,uj)+ϑ3φini(xi-uis,uj)
(20)
其中，系数ϑ1、ϑ2、ϑ3为：
(21)
式(20)为Boltzmann-Shakhov方程沿特征线离散的一般形式，其右端前两项和平衡态分布函数相关反映粒子之间的碰撞作用，将和式(17)右端第一项积分项对应；式(20)右端第三项和初始非平衡态分布函数相关反映粒子的迁移效应，将和式(17)右端第二项相对应。通过对比式(17)和式(20)中的非平衡态分布函数相关项，并令ϑ3=e-s/τ，可确定权系数α为：
(22)
式(22)表明权系数α数学上直接依赖于时间比值τ/s，物理上的讨论将在后面第3节给出。确定了α后将其带入到式(21)，可得系数ϑ1、ϑ2和ϑ3具体形式：
(23)
为了进一步分析式(20)的数学物理特性，这里将从两种极限流动的角度来讨论，首先将式(20)整理为以下形式：
φ(xi,t+s,uj)=φs(xi,t+s,uj)+
β1[φini(xi-uis,uj)-φs,ini(xi-uis,uj)]+
β2[φs,ini(xi-uis,uj)-φs(xi,t+s,uj)]
(24)
其中：
β1=e-s/τ,β2=(1-e-s/τ)τ/s
(25)
对于Kn→0的宏观连续流动问题，物理上在一个时间步长内粒子经历的碰撞次数足够多，即τ≪s，数学上β1→0，β2→τ/s，式(24)可化为：
φ(xi,t+s,uj)≈φs(xi,t+s,uj)+
≈φs(xi,t+s,uj)-τDtφs(xi,t+s,uj)
(26)
其中Dt=∂t+uj∂xj。对式(8)进行Chapman-Enskog展开至O(τ)(对应为宏观N-S方程)同样可以得到式(26)[12]，因此说明在宏观连续流动区域式(20)与N-S方程一致, 此时流动尺度为宏观动力学尺度。
对于Kn≫0的稀薄流动或微尺度流动问题，物理上τ≫s，从而β1→1，β2→1，式(24)化为：
φ(xi,t+s,uj)→φini(xi-uis,uj)
(27)
上式其实就是自由分子流运动的解，此时流动对应微观动理学尺度。
由上述分析可知，Boltzmann-Shakhov方程沿特征线离散的一般形式(20)将宏观动力学尺度和微观动理学尺度有机地结合到一起，其中引入的系数ϑ1、ϑ2和ϑ3依赖于比值τ/s且可根据当地流动领域和流动尺度的不同进行自适应调节，比值τ/s在调节中起到了关键作用。
2 尺度自适应的离散统一气体动理学格式
尺度自适应的离散统一气体动理学格式(Scale Adaptive Discrete Unified Gas Kinetic Scheme, SADUGKS)可视为在原始DUGKS格式[24]的基础上进行改进，二者主要区别在于权系数α的形式不同。
基于有限体积法框架，将式(16)在网格单元Vj及时间步长(tn,tn+Δt)内积分，得：
(28)
(29)
对式(28)中通量的时间积分项采用中点公式处理：
(30)
对式(28)中碰撞项的时间积分采用加权梯形公式：
(31)
其中α(Δt)为式(22)给出的权系数:
(32)
需要指出的是在原始DUGKS格式中[24]，式(31)中的权系数取为αDUGKS≡0.5，即代表碰撞项的显式和隐式算术平均。
将式(30)和式(31)代入到式(28)中，可得：
(33)
从式(33)可以看出该式右端包含了通量和碰撞项的隐式部分，因此在使用该式进行迭代递推之前必须先解决通量和碰撞项的隐式离散问题，下面对其分别进行处理。
(34)
(35)
上式也等价于:
(36)
其中系数φ1和φ2为，
(37)
则式(34)可化为：
(38)
(39)
(40)
由式(4)、(15)和(35)可知:
(41)
结合式(38)位于界面xi,b处，tn+r时刻的宏观守恒量为:
Q(xi,b,tn+r)
(42)
(43)
由相容关系式(15)同时结合定义式(43)可知:
(44)
将式(43)代入到式(33)，整理得：
(45)
再次指出的是在原始DUGKS格式中，从式(34)至式(39)，式(43)和式(45)中出现的所有权系数α恒取为αDUGKS≡0.5。
需要说明的是上述在推导SADUGKS中速度空间是连续的，实际在使用SADUGKS时还需将速度空间离散，其离散方法与现有DUGKS[23-24]格式相同。关于边界条件的处理分两类，针对宏观连续流动问题，SADUGKS格式在固体壁面处采用了无滑移边界条件；而对较大Kn流动问题，则采用了完全漫反射边界条件，具体实现与DUGKS格式[23-24]亦相同，不再赘述。
3 尺度自适应的离散统一气体动理学格式的性质
首先，作为DUGKS格式的一种改进，SADUGKS格式同样具有渐近保持性[23-24]，这一点可以清楚地从式(26)看出。
其次，SADUGKS格式具有尺度自适应特性，能够根据当地流动特征尺度的不同进行自适应调节。SADUGKS格式中，时间步长按照CFL条件确定：
(46)
其中，Δ为空间网格尺度，CFL为CFL数；|u|max为粒子最大离散速度，通常与声速a同量级。同时，碰撞松弛时间τ可近似为：
(47)
τ/Δt～ls/Δ～Knlocal
(48)
其中Knlocal为当地努森数，反映当地的流动状态和流动尺度。在SADUGKS格式中，用来计算界面通量的式(40)和碰撞项加权离散的式(31)中皆出现权系数α，式(32)给出权系数α的定义式反映出α依赖于比值τ/Δt，即依赖于当地努森数，因此α取值的大小会随着当地流动状态和尺度的不同而进行自适应调节。
α→1,τ≪Δt
α→1/2,τ≫Δt
(49)
≈Ωφ(xi,t+Δt,uj)
(50)
可见数学结果与上述物理分析一致。
(51)
可见此时数学结果也与上述物理分析相一致。
在原始DUGKS格式中，本文所有出现权系数α的地方都取为常数αDUGKS≡1/2。通过上面的分析可以看出，这种常数取法主要是从数值离散的角度考虑，并未充分考虑在不同流动领域碰撞项对应物理过程的不同，尤其没有考虑两个时间尺度即计算时间步长和碰撞松弛时间之间的相对大小对粒子碰撞过程宏观表现的影响。
根据上述分析可知，式(32)给出的权系数α紧密依赖于时间比值τ/Δt，并可根据当地流动状态和流动尺度的不同进行自适应调节。当地努森数越小、越接近连续流动时α越接近于1；当地努森数越大、流动非平衡效应越强时α越接近于1/2。综合来说，SADUGKS是一种尺度自适应的格式，适用于从连续流动到滑移流动、过渡流动及自由分子流所有流动领域的问题，在第4节将通过若干典型算例对这一点进行检验。
4 尺度自适应的离散统一气体动理学格式在可压缩流动中的应用
4.1 一维激波结构问题
本问题中将用SADUGKS格式研究一维激波结构问题，并将结果与文献[36,37]的结果进行对比。气体为氩气，分子内自由度设置为K=0，普朗特数Pr=2/3，比热比γ=5/3，黏性系数与温度相关：μ∝Tw，其中w与分子模型相关。分子平均自由程λ与黏性系数直接相关[24]：
(52)
选取计算区域为-25λ1≤x≤25λ1，并采用100个均匀网格将其离散，每个网格单元大小Δx=0.5λ1，由于网格单元Δx和分子平均自由程λ1量级相当，计算将给出激波内部结构中物理量的变化。速度空间采用32×32半平面Gauss-Hermite积分离散[38]。初始时刻，x≤0处的流场采用激波上游条件及其Maxwell平衡态分布进行初始化，x>0处的流场采用激波下游条件及其Maxwell平衡态分布进行初始化。左右边界均采用外推边界条件,CFL数设置为0.5。
本文采用硬球模型w=0.5，图1展示了Ma=1.2数值模拟结果，并与文献[36,37]中的结果进行了对比。其中物理量的无量纲化方法如下：
(53)
另外，参照文献[36]给出固定激波所在位置的方法，选取ρ(x0)=(ρ1+ρ2)/2处为坐标原点，将坐标进行无量纲化：
(54)
从图1中可以看出，SADUGKS的结果与UGKS[37]及完整Boltzmann方程[36]都能很好地吻合。由于本问题是分辨激波内部结构，流动对应尺度为分子平均自由程尺度，按照前面关于权系数α的讨论，此时α应接近0.5，从图1中也可以清楚看出全场α范围在0.503至0.505之间。
图1 SADUGKS格式求解激波结构问题的密度、温度、热流、应力与权系数α分布(Ma=1.2)Fig.1 Distribution of density, temperature, heat flux, stress and α for shock structure problem (Ma=1.2)
4.2 从激波结构分辨格式到激波捕捉格式
在4.1小节问题中，由于所用网格单元Δx和分子平均自由程λ1量级相当，因此计算结果可以分辨激波内部结构。接下来本节将重点讨论当网格单元大小变化时对数值结果的影响，以验证SADUGKS格式的尺度自适应特性。
气体参数、激波上下游初值等参数与4.1小节一致，CFL数取为0.9。速度空间采用32×32半平面Gauss-Hermite积分离散[38]。选取计算区域为-50Δx≤x≤50Δx，区域总长为L=100Δx，采用100个均匀网格将其离散，网格单元Δx选取了四种不同大小，分别是Δx=0.5λ1，λ1，10λ1和100λ1。图2给出了无量纲密度、温度以及权系数α在不同Δx时的空间分布曲线，其中横坐标x*为采用激波上游的分子平均自由程λ1进行无量纲化后的空间位置，x*具体形式见式(54)，同时无黏Euler方程的精确解也画在图中作为对比。从图2可以看出，当Δx较小，满足Δx～λ1时，真实激波内部结构可以被分辨出来，此时SADUGKS表现为激波结构分辨格式(shock structure resolving scheme)；随着Δx增加，如当Δx≥10λ1时，SADUGKS给出的结果将逼近无黏Euler方程的精确解，说明对于网格单元远大于分子平均自由程时，激波被作为理想间断面来处理，SADUGKS通过自身格式黏性最后得到的是数值激波，此时SADUGKS表现为激波捕捉格式(shock capturing scheme)。
这个例子说明SADUGKS可以根据比值Δx/λ1的变化而在激波结构分辨格式和激波捕捉格式之间做出自适应调整，这对于实际复杂工程问题获得所用计算网格下的合理物理解是很有益的。因为一个问题在求解前进行网格划分时并不知道流场物理量的局部细节特性，所以划分的网格带有较强的人为因素，但是SADUGKS格式自身可以通过一套统一的算法自动获得流场在所用网格尺度下的合理物理解。
SADUGKS格式在激波结构分辨格式和激波捕捉格式之间的切换得益于权系数α的引入，从图2权系数α的变化曲线可以清楚看出：当Δx≫λ1时，α取值接近1，流动可认为是连续流动；当Δx～λ1时，α取值接近0.5，流动趋于稀薄流动，流动的非平衡效应增强。需要强调的是，通常所说的连续流动和稀薄流动都是一个相对的概念，和选用的网格单元尺度与分子平均自由程尺度相对大小有关。对同一个流动，网格单元尺度选用的不同，可获得流动的信息量不同，描述流动的方法也是有区别的。网格单元尺度越小，观察到的流动就越接近于一个个运动的微观粒子，获得的信息量就越丰富；而当网格单元尺度远大于分子平均自由程尺度时，可获得的流动信息是大量微观粒子运动的统计平均，此时就是连续介质框架下的描述方法。通过上述讨论说明SADUGKS是一种尺度自适应的气体动理学格式，在该格式中，计算网格发挥着积极作用，最终可得到计算网格所对应尺度下的流动物理信息。
作为对比，对于传统在连续介质框架下基于N-S方程建立的CFD方法来说，N-S方程建立后被作为一个数学偏微分方程来处理，网格进行加密的主要目的是减小N-S方程离散误差，网格的作用未被充分发挥出来。另外，从物理上考虑，N-S方程是基于连续介质框架得到的，N-S方程蕴含的流动物理尺度是远大于分子平均自由程的，因此即便将网格加密到分子平均自由程的量级，也只是获得N-S方程作为一个数学偏微分方程更为精确的数值解，但所获得的结果不可能超过N-S方程物理上成立的范围，即无法给出在分子平均自由程尺度上的流动物理信息。
图2 SADUGKS格式求解激波结构问题的密度、温度与权系数α随网格大小变化，CFL=0.9, Ma=1.2Fig.2 Density, temperature and α profiles of the shock structure with different cell sizes，CFL=0.9, Ma=1.2
4.3 从自由分子流领域到连续流领域
Sod问题为一维Riemann问题中的一个经典问题，其无黏解包含激波、膨胀波及接触间断多种波系结构。该问题被广泛用来检验数值格式捕捉间断的能力，其初始条件如下：
(55)
本文这里考虑在上述初值下，引入气体黏性后的流动。气体考虑为空气，分子内自由度设置为K=2，普朗特数Pr=0.72，比热比γ=1.4。气体黏性系数与温度满足μ=μ0(T/T0)w，其中w与分子模型相关，这里采用硬球模型w=0.5，μ0为参考黏性系数，T0为参考温度。参考分子平均自由程λ0与参考黏性系数直接相关：
(56)
其中ρ0为参考密度。本文中，选用左侧(x<0)初值ρ1、p1、T1分别作为参考密度、压力和温度。
μ0=10时的密度、温度、速度以及此时流场中式(32)给出的权系数α随空间分布如图3(a)所示。在此黏性系数下，左侧初值对应分子平均自由程λ0=12.77，比值λ0/Δx=1277，此时流动对应为自由分子流。从图2中可以看出，SADUGKS结果与UGKS结果基本完全重合，且这种远偏离平衡态的流动与无黏Euler解有着明显的区别：一方面，由于流动更为稀薄，分子平均自由程远大于网格单元，因此空间网格可以分辨包括激波在内的各种波系内部结构；另一方面，流动黏性很大，耗散作用增强，各种波系被显著光滑抹平，物理量在整个空间上连续变化。图中给出的权系数α在整个空间几乎是一条平直线，值为0.5，这一点和前面第3节中提到的对于大Kn数的流动α→1/2结论一致。当黏性系数μ0减小到1时，左侧初值对应分子平均自由程λ0=1.277，比值λ0/Δx=127.7，从图3(b)可以看出，相比μ0=10时物理量分布变化不大，同样此时权系数α在整个空间都接近0.5。
(a) μ0=10.0
(b) μ0=1.0
(c) μ0=0.1
(d) μ0=10-3
(e) μ0=10-5图3 SADUGKS格式求解Sod问题的密度、温度、速度和权系数α分布Fig.3 Distribution of density, temperature, velocity and α for Sod problem by SADUGKS
随着黏性系数μ0减小到0.1时，左侧初值对应分子平均自由程λ0=0.1277，比值λ0/Δx=12.77，相比前两个较大黏性系数的工况来说，图3(c)给出的物理量分布已出现一些变化，尤其速度峰值已明显增加，并且权系数α在整个空间也出现起伏，不再是平直直线了。
当黏性系数μ0减小到10-3时，左侧初值对应分子平均自由程λ0=1.277×10-3，比值λ0/Δx=0.1277，在此时的网格下已不能分辨激波内部结构了，从图3(d)中可见，SADUGKS结果逐渐接近无黏Euler解。权系数α在整个空间已出现显著变化，从最左侧的0.5293逐渐减小到右侧的0.5033。
当黏性系数μ0进一步减小到10-5时，左侧初值对应分子平均自由程λ0=1.277×10-5，比值λ0/Δx=1.277×10-3，此时流动属于连续流领域，并且由于黏性非常小，从图3(e)可以看出，SADUGKS的结果几乎与无黏Euler解完全重合。当前空间网格尺度远大于分子平均自由程尺度，对应网格不再能分辨激波内部流动细节，激波表现为宏观物理量的间断。此时，SADUGKS表现为激波捕捉格式。从权系数α分布图中可以看到α随空间显著变化，从左侧的0.9716逐渐减小到右侧的0.7660，这反映当地努森数和当地特征尺度随空间的变化。
本算例计算网格单元大小不变，改变μ0实现分子平均自由程λ0的变化，验证了SADUGKS格式的尺度自适应特性。对于实际工程问题，使用一种能够根据当地流动尺度进行自适应调节的统一算法显得尤为重要，而SADUGKS格式正是这样一种尺度自适应的算法，能够反映出当地流动物理特性。且根据当地流动尺度不同，SADUGKS格式可在激波结构分辨格式和激波捕捉格式之间进行自适应切换。
4.4 高速可压缩Couette流动问题
最后将SADUGKS格式应用于高速可压缩Couette流动问题，流动参数与采用IP方法[39]进行计算的文献一致。流动示意图如图4所示，两无限大平行平板之间距离为1，两板之间充满氩气，分子内自由度K=0，普朗特数Pr=2/3，比热比γ=5/3，黏性系数与温度相关μ=μ0(T/T0)w，采用变径硬球模型w=0.81。两板温度固定为273 K，下板静止，上板在自身平面内匀速运动，速度为U2=300 m/s，对应马赫数Ma约为0.97。流体与壁面之间存在剪切作用，同时黏性耗散作用会使得流场内部温度升高。通过改变初始氩气密度以改变流动努森数Kn，本文共选取五组Kn=0.01，0.1，1.0，10和100，涵盖了从连续流动到自由分子流的流动领域。不同Kn下流动将表现出不同的温度及速度分布。
两板间布置数目为31的均匀网格，速度空间采用28×28半平面Gauss-Hermite积分离散[38]。入口与出口采用周期性边界条件，上下壁面采用完全漫反射边界条件[19]。CFL数设置为0.5。
图4 高速可压缩Couette流动示意图Fig.4 Schematic diagram for the compressible Couette flow
不同Kn下的速度剖面如图5所示，从图中可以看出，当Kn=0.01时，流动可认为仍满足连续介质假设，但是上下壁面已出现轻微的速度滑移现象。随着Kn不断增加，非平衡效应变强，壁面处的滑移效应逐渐变得显著，尤其当Kn=100时，板间气体速度分布趋于均匀，约为150 m/s，和上、下板速度远不相同。
图5 SADUGKS计算所得不同Kn下的速度剖面图Fig.5 Velocity profiles for the compressible Couette flow at different Knudsen numbers
图6所示为不同努森数下的温度剖面。当Kn=0.01时，流动在壁面已出现轻微温度跳跃现象。由于上板克服壁面黏性剪切做功最终转化为热能，并逐渐向流道中间集中，造成温度峰值出现在中心对称线上。随着Kn不断增加，非平衡效应变强，壁面处的温度跳跃现象更为明显，板间气体的温度梯度逐渐减小。尤其当Kn=100时，板间流动趋于自由分子流，板间气体温度几乎相同，约为308 K，和上、下板温度远不相同。
为了反映气体Pr数对结果的影响，这里考虑了Pr=1时Kn=0.01和Kn=0.1两种工况，并和前面Pr=2/3的对应结果进行了对比，见图7给出的温度分布曲线，可以清楚看出Pr=1时的结果与Pr=2/3及文献IP方法[39]的结果出现了较大的偏差。从式(10～14)可以看出，Boltzmann-Shakhov模型中取Pr=1时退化为Boltzmann-BGK模型，这说明为了正确反映气体Pr数的影响从而获得准确温度场，基于Boltzmann-Shakhov模型构造计算格式是适当的。
无论对于速度还是温度，在所有Kn情况下，SADUGKS格式所得结果均能与文献IP方法[39]的结果很好地吻合，表明了SADUGKS格式对于宽Kn流动计算的有效性和尺度自适应特性。
图6 SADUGKS计算所得不同Kn下的温度剖面图Fig.6 Temperature profiles for the compressible Couette flow at different Knudsen numbers
图7 Pr=1和Pr=2/3的结果对比(Kn=0.01和Kn=0.1)Fig.7 Temperature profiles for the compressible Couette flow at different Knudsen numbers with Pr=1 and Pr=2/3
5 结 论
文中提出了Boltzmann方程沿特征线离散的一般形式，式中碰撞项的离散采用显式和隐式加权平均的方法，权系数α依赖于当地碰撞松弛时间和计算时间步长的比值。通过权系数α的引入，对文献现有离散统一气体动理学格式(DUGKS)进行了改进，提出了具有尺度自适应特性的离散统一气体动理学格式(SADUGKS)。对SADUGKS进行分析，明确权系数α与反映当地流动尺度的当地努森数直接相关，表明SADUGKS格式能够根据当地流动尺度进行自适应调节。SADUGKS格式的自适应特性和有效性通过典型可压缩流动进行了检验，取得了较好的计算结果，表明SADUGKS格式适用于宽范围Kn变化、跨流域的多尺度流动问题的数值模拟。后续将选用更为复杂的流动问题对SADUGKS做进一步的检验和验证。
