基于非線性耦合本構關系的改進邊界條件

江中正，趙文文，袁震宇，陳偉芳

浙江大學 航空航天學院，杭州 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粒子仿真結果。