高超聲速流存在局部稀薄效應(yīng)的一個判據(jù)及相應(yīng)的流動特性

陳 杰, 趙 磊(1. 天津大學(xué) 機(jī)械工程學(xué)院力學(xué)系, 天津 300350; 2. 中國空氣動力研究與發(fā)展中心 超高速空氣動力研究所, 四川 綿陽 621000)

0 引 言

我國的空氣動力學(xué)自20世紀(jì)50年代中在錢學(xué)森先生回國后的全面規(guī)劃建議下，至今已得到極大的發(fā)展[1]。其中計算流體力學(xué)(Computational Fluid Dynamics, CFD)是非常突出的一個方面[2]，以張涵信先生的團(tuán)隊為代表。他們的一個鮮明特點是無論是在發(fā)展算法，還是在分析數(shù)值結(jié)果時，都以清晰的物理分析為基礎(chǔ)。一方面成功構(gòu)建了多種有效的算法，從而在實際的復(fù)雜流動計算中提供了高保真數(shù)據(jù)[3]。另一方面又能給出清晰的流動物理圖像，為航空航天技術(shù)的發(fā)展做出了巨大貢獻(xiàn)。

但在近年來的近空間飛行器的研發(fā)中發(fā)現(xiàn)，原有的空氣動力學(xué)及相應(yīng)的CFD不能完全應(yīng)對遇到的情況，需要研究新的空氣動力學(xué)問題[1]。

近空間高超聲速飛行器飛行高度一般在40-70 km的高空, 飛行速度一般在5-20 倍聲速的范圍[1]，其研發(fā)在中國及世界范圍內(nèi)已受到越來越多的關(guān)注。在此類飛行器的飛行過程中，流場局部的高溫和低的壓力等因素會導(dǎo)致局部區(qū)域內(nèi)存在氣體稀薄效應(yīng)，使得連續(xù)介質(zhì)模型失效，傳統(tǒng)的CFD方法對氣動力，尤其是摩擦阻力及熱流的預(yù)測可能會有較大誤差。

研究涉及稀薄氣體效應(yīng)的流動，首要的問題就是用什么參數(shù)判別稀薄效應(yīng)的程度。錢學(xué)森先生最先根據(jù)稀薄程度將稀薄氣體流分為三大流域[4]。根據(jù)Knudsen數(shù)的大小，其定義為平均分子自由程λ與飛行器某一特征長度L的比值，三個流域依次為滑移流(0.01 10)。當(dāng)時關(guān)注的主要是飛行器的壓力分布，這樣的劃分基本上是可以的。后來，對于簡單管道微尺度流動，L通常取管道的水力直徑，當(dāng)Kn數(shù)小于0.2時，連續(xù)介質(zhì)模型結(jié)合二階滑移邊界條件可以獲得與實驗一致的結(jié)果[5]。

但對于高速飛行器周圍的流場，情況極為復(fù)雜。用一個參數(shù)Kn不足以刻畫其流動狀態(tài)。針對不同問題，要選取不同的物理特征長度。例如，在飛行器頭部附近，溫度和壓力變化迅速，用飛行器整體的特征長度和來流的氣體平均分子自由程不能反映當(dāng)?shù)氐那闆r，即原先定義的Kn數(shù)并不適用，需要考慮局部的稀薄效應(yīng)。因此，最直接的方法是使用局部的平均自由程和局部流動的特征尺度，如激波脫體距離，或一些宏觀量梯度的標(biāo)尺長度[6]，建立局部的Kn數(shù)。最早Bird針對膨脹流的連續(xù)介質(zhì)失效判定提出的P參數(shù)[7]，它正比于Maλ(dρ/dx)/ρ。Boyd等認(rèn)為在高超聲速稀薄流動中應(yīng)同時考慮局部的密度和溫度梯度[8]，提出了KnGLL=λ/(Q/|dQ/dl|)，Q為密度或溫度，導(dǎo)數(shù)應(yīng)沿梯度最大方向l。數(shù)值計算表明對于一維正激波和二維圓柱繞流的連續(xù)失效閥值為KnGLL-D=0.05。隨后，Wang等將其寫為KnQ=λ/(Q/|Q|)[9]。KnGLL、KnQ是目前在NS-DSMC耦合方法[10-11]中常用的對計算區(qū)域進(jìn)行劃分的參數(shù)之一。在一些動理論方法中，連續(xù)失效準(zhǔn)則可以由某一類物理量(如剪應(yīng)力、熱流、速度分布函數(shù)等)與連續(xù)方法計算的值之間的相對誤差的形式建立。 Lockerby等提出了KnL=max(Knσ,Knq)[12]，其中Knσ、Knq分別定義為R-13矩方法[13]所求的剪應(yīng)力、熱流與在R-13流場結(jié)果上依據(jù)線性本構(gòu)關(guān)系所求的剪應(yīng)力、熱流之間的相對誤差。隨后，Meng等直接從動理論的分布函數(shù)出發(fā)，定義了基于分布函數(shù)的失效準(zhǔn)則主要用于多尺度動理論方法[15]中低階模型與高階模型的選擇。另外，熱力學(xué)不可逆過程熵增條件在非平衡效應(yīng)判斷和動理學(xué)模型構(gòu)造中引起了關(guān)注[16]。Camberos等基于只要是非平衡狀態(tài)都有熵產(chǎn)生，將熵產(chǎn)率作為可以普遍預(yù)測非平衡的重要參數(shù)之一[17]，并作為連續(xù)失效的參數(shù)。Xiao等基于從非平衡態(tài)到平衡態(tài)熵增函數(shù)數(shù)學(xué)模型收斂為Rayleigh-Onsager函數(shù)提出了新的近平衡失效參數(shù)但現(xiàn)有參數(shù)對于近空間飛行器的壁面摩擦阻力和熱流問題并不適用，其壁面附近流動的特點為高溫、低密度、強(qiáng)剪切流動。

在稀薄氣體數(shù)值模擬方面，Bird提出的直接模擬Monte Carlo方法(Direct Simulation Monte Carlo，DSMC)至今被認(rèn)為是在模擬稀薄氣體流動，尤其是在過渡流區(qū)域最成功的方法之一[19-21]。該方法在涉及高速、高溫、化學(xué)反應(yīng)等的稀薄氣體流動實際應(yīng)用方面有很多成功例子[22]。為了解決DSMC在低速流動中的統(tǒng)計漲落問題，我國沈青、樊菁發(fā)展了信息保存法(Information Preservation, IP)[23]。另一方面，被認(rèn)為可以描述整個流域的氣體動理論基本方程——Boltzmann方程求解十分困難[24]。最近我國學(xué)者在Boltzmann方程簡化模型的數(shù)值求解方法，特別是跨流域統(tǒng)一格式方面，取得了重要進(jìn)展[25-28]，并引起了國外同行的普遍關(guān)注。但是對于三維近連續(xù)域氣體動力學(xué)問題，其計算量相比傳統(tǒng)CFD依舊很大，現(xiàn)階段計算資源尚難滿足近空間高超聲速飛行器的工程計算。此外，大部分動理學(xué)格式建立在Boltzmann-BGK簡化方程之上，此簡化模型的理論基礎(chǔ)為氣體偏離平衡狀態(tài)不太遠(yuǎn)。如本文下面所述，這一假定在實際的近空間飛行條件下與實際情況相去可能很遠(yuǎn)。另外一些學(xué)者從連續(xù)性方法出發(fā)，在小Kn數(shù)下從流動守恒和本構(gòu)關(guān)系方面發(fā)展了一些基于連續(xù)模型的方法。Chapman和Enskog將分布函數(shù)在Maxwell平衡態(tài)附近展開為正比于Kn數(shù)冪次的級數(shù)[29]，Burnnett以此為基礎(chǔ)建立了Burnett方程等[30]。Grad用Hermite多項式展開分布函數(shù)[13]，得到了13個Boltzmann方程的矩方程，即Grad十三矩方程, 用于過渡流模擬[31]。Eu從不可逆熱力學(xué)過程出發(fā)，對Boltzmann方程模型化，建立了非線性耦合本構(gòu)關(guān)系(Nonlinear Coupled Constitutive Relations，NCCR)理論[16,32]，提供了一種新的可能，該方法在真實三維稀薄流中的適用性仍有待驗證[33]。

綜上分析，存在局部稀薄效應(yīng)的近連續(xù)流的求解是目前最困難的問題之一，該問題正好對應(yīng)于近空間高超聲速飛行器的研發(fā)。傳統(tǒng)連續(xù)介質(zhì)模型中的線性牛頓黏性定律和傅里葉導(dǎo)熱定律是否能給出足夠精確的結(jié)果需要重新檢驗。而氣體動理論方法由于其巨大的計算量難以應(yīng)用于實際復(fù)雜外形。所以研究局部稀薄效應(yīng)現(xiàn)象及其判別參數(shù)、認(rèn)識傳統(tǒng)的CFD適用范圍及問題所在、考慮如何將局部稀薄效應(yīng)包含于CFD計算中有非常重要的意義。

本文將選取具有代表性的高速剪切流為研究對象，對局部非平衡狀態(tài)以及對剪應(yīng)力計算等問題進(jìn)行詳細(xì)分析，進(jìn)而找到一個可以反映剪切流局部非平衡效應(yīng)的無量綱參數(shù)，并揭示此參數(shù)和剪切力與傳統(tǒng)連續(xù)介質(zhì)模型下剪切力偏移量的規(guī)律。

1 模型描述及計算方法

1.1 問題描述

周恒和張涵信在討論與近空間飛行器對應(yīng)的空氣動力學(xué)面臨的新問題[1]時指出“氣體分子運動論用于空氣動力學(xué)時，仍然認(rèn)為分子熱運動速度的分布對Maxwell分布只有小的偏離。而這在流速很高的情況時是很可疑的。文獻(xiàn)[1]中采用四種不同的數(shù)值方法模擬了長1 m，在70 km 高空以馬赫數(shù)15零迎角飛行的零厚度平板的繞流問題(見文獻(xiàn)[1]中圖1的結(jié)果)，如暫時以圖示計算結(jié)果為參考來估算，可知約在平板法向4 cm處(見文獻(xiàn)[1]圖1(h))，氣體密度最低，約為0.4，因而該處分子自由程約為5 mm(來流分子自由程為2 mm)。同時該處溫度約為10(見文獻(xiàn)[1]圖1(g))，對應(yīng)于2000 K(來流溫度為200 K)。該處速度梯度約為32 m·s-1每毫米(見文獻(xiàn)[1]圖1(e))，或在法向相差一個分子自由程的兩層間速度差約為160 m·s-1。如果假設(shè)分子熱運動平均速度與溫度的二分之一次方成正比，則該處的分子熱運動平均速度約為1190 m·s-1(在室溫300 K下，分子熱運動平均速度約為450 m·s-1)。粗略地講，相隔一個分子自由程的兩層間的分子碰撞是產(chǎn)生剪切力和熱傳導(dǎo)的主要機(jī)制。而在上述這個速度很高的例子下，兩層分子的碰撞要使各自的平均速度改變160 m·s-1，相對于其熱運動平均速度，顯然不是一個小量，它會使各自的分子熱運動的速度分布顯著偏離標(biāo)準(zhǔn)的Maxwell分布?！盵1]

從上面所引的周恒和張涵信的分析中可以看出，他們認(rèn)為在有氣體稀薄效應(yīng)而剪切又很強(qiáng)的情況下，分子熱運動速度的分布對Maxwell分布只有小的偏離這一假設(shè)是很可疑的。而起決定性作用的參數(shù)是相隔為一個分子自由程的兩層流體的速度差和分子自由運動速度之比這一無量綱參數(shù)，即

1.2 Couette流

本文重點研究剪切力受氣體稀薄效應(yīng)的影響，因此選取最簡單的剪切流，平面Couette流為分析模型，即由兩個相對運動的無限大平行平板表面剪切力驅(qū)動的流動。針對中等稀薄程度(即傳統(tǒng)的Kn數(shù)為0.005～0.5之間)，但存在強(qiáng)剪切的情況進(jìn)行詳細(xì)的研究。圖1描述了Couette流示意圖，壁面溫度為Tw=273.15 K。以氬氣作為模擬氣體，初始溫度為273.15 K，此時黏性系數(shù)為μ=2.12×10-7kg·m-1·s-1。當(dāng)壓力為標(biāo)準(zhǔn)大氣壓p0，初始溫度為T0=Tw時，分子平均自由程為λ0=62.50 nm，相應(yīng)的分子溫度的最可幾速度為Vmp=337 m·s-1。在計算中，壁面速度U分別取Vmp、3Vmp、5Vmp、8Vmp，兩板間距為L=200λ0，模擬中通過改變初始分子數(shù)密度而改變氣體的稀薄程度。

圖1 Couette流示意圖Fig.1 Couette flow

1.3 數(shù)值模擬方法

本文采用基于Bird提出的標(biāo)準(zhǔn)DSMC算法建立的自編程程序。在計算中，碰撞對的選取采用NTC(No Time Counter)方法，碰撞模型采用硬球模型，即黏性系數(shù)與溫度的根號成正比。在壁面處采用完全漫反射條件。在所有DSMC計算中，計算網(wǎng)格均小于當(dāng)?shù)氐钠骄杂沙痰娜种?、時間步長小于平均碰撞時間的五分之一，使用了大量的模擬分子，初始每個網(wǎng)格內(nèi)分子數(shù)為100個，整體樣本統(tǒng)計量在108的量級，以確保得到較好的粒子速度分布函數(shù)。同時，我們驗證了網(wǎng)格尺寸對粒子速度分布函數(shù)的影響，網(wǎng)格尺寸在滿足基本條件基礎(chǔ)上繼續(xù)減小，對分布函數(shù)產(chǎn)生的影響非常小。

在本文2.3節(jié)中，我們用傳統(tǒng)的CFD模擬了高超聲速小鈍度鈍楔模型的繞流。其飛行參數(shù)對應(yīng)于近空間飛行器的參數(shù)范圍。目的是考察實際飛行條件下，是否真的需要考慮氣體的局部稀薄效應(yīng)。模擬采用無滑移等溫壁面，遠(yuǎn)場激波外為來流條件，出口為外推邊界條件。空間離散分別采用五階WENO格式和六階中心格式。因是定常流計算，時間方向采用2階Runge-Kutta方法。

2 結(jié)果分析

本文針對高速剪切流主要研究兩個問題：1) 在弱稀薄但強(qiáng)剪切的情況下，氣體偏離平衡狀態(tài)的程度？即氣體分子熱運動速度分布相對于Maxwell分布的偏離程度。如的確處于不平衡狀態(tài)，參數(shù)Zh是否是表征非平衡程度的一個合適的參數(shù)？對這一問題，基于DSMC方法的微觀特性，統(tǒng)計獲得粒子的概率密度分布函數(shù)與Maxwell平衡態(tài)比較是最直觀的方法。2.1節(jié)將給出這方面的結(jié)果。2) 在判斷了氣體是否處于平衡狀態(tài)之后，首要關(guān)心的問題就是對實際剪應(yīng)力的影響?；谶B續(xù)介質(zhì)的牛頓黏性假設(shè)在什么范圍可以給出正確的剪應(yīng)力計算？如果不正確，偏差如何？2.2節(jié)將詳細(xì)討論這一問題。

2.1 概率密度分布函數(shù)

我們計算了不同密度、不同壁面速度下30個算例，畫出了中心點的概率密度分布函數(shù)，計算了當(dāng)?shù)氐腪h值，我們發(fā)現(xiàn)，可以通過Zh的值將分布函數(shù)的形狀分為三種。

圖2、3、4分別給出了三種典型的概率密度分布函數(shù)以及相應(yīng)的當(dāng)?shù)豘h參數(shù)的數(shù)值。中心點的宏觀速度接近于0， 粒子速度分布以速度v=0對稱。在稀薄程度較小、剪切強(qiáng)度較小時，速度分布函數(shù)與Maxwell平衡態(tài)給出的一致，如圖2(a)所示，此時中心點的Zh為0.02。隨著氣體稀薄程度、剪切強(qiáng)度的增加，速度分布函數(shù)逐漸偏離Maxwell分布且偏離的程度越來越大。以中心點當(dāng)?shù)氐腪h值可以將分布函數(shù)的狀態(tài)大致分為三種如圖2所示，即單峰的近似正態(tài)分布(Zh<1.0)，三峰分布(1.01.5)。

圖2 概率密度函數(shù)(Zh<1.0)Fig.2 Probability density function(Zh<1.0)

圖3 概率密度函數(shù)(1.0

在與周恒院士的討論中，他用一個更簡單的物理模型定性地說明了上述現(xiàn)象。他以相距各為一個分子自由程，平均速度差為u的三層流體作為分析對象(如圖所示5)，分析中間一層流體的分子運動速度沿流向分量的分布函數(shù)的可能變化如下：

圖4 概率密度函數(shù)(Zh>1.5)Fig.4 Probability density function(Zh>1.5)

圖5 三層流體示意圖Fig.5 Schematic illustration of three-layers in fluid

設(shè)初始時各層流體分子運動速度分布函數(shù)都是Maxwell分布，其中心速度就是該層流體的平均速度。當(dāng)上層流體分子來到中間層，與該層分子碰撞后，其總的效果就是使得速度大于平均速度的分子數(shù)增加，即原先Maxwell分布的正向某處(略小于u值)附近的分子數(shù)增加，相應(yīng)地原先對應(yīng)于峰值處的分子數(shù)減小。而下層流體的分子的作用類似，即使得速度小于平均速度的分子數(shù)增加，使得原先Maxwell分布的負(fù)向某處(略大于-u值)附近的分子數(shù)增加，相應(yīng)地原先對應(yīng)于峰值處的分子數(shù)減小。中間層的分子也可以相互碰撞，但并不會將上述的正負(fù)效應(yīng)相抵消，因為多出來的分子運動速度大于平均速度的分子正好和相應(yīng)的分子運動速度小于平均速度的分子相碰的概率并不大。因此，總的效果是原先Maxwell分布的峰值處會降低，而兩側(cè)則會增高。

當(dāng)Zh小于1，對應(yīng)于u值較小的情況，因此其總的效果只是原先的Maxwell分布峰的兩側(cè)近處略微增厚，峰略微下降，因而總體仍可保持一個對稱的單峰分布如圖2所示。當(dāng)Zh值大到一定程度，u的值增加，原先Maxwell分布的峰的兩側(cè)較遠(yuǎn)處的分子數(shù)增加，就會出現(xiàn)兩個峰而不僅是原先峰的兩側(cè)近處的增厚。但這時原先的峰還沒有降得足夠多，因此速度分布函數(shù)就會出現(xiàn)三個峰的情況，如圖3所示。當(dāng)Zh繼續(xù)增大，則上下層分子的影響下形成的兩側(cè)的峰相隔更遠(yuǎn)，且粒子速度偏向u和-u的分子數(shù)隨Zh增大而增加，使得主峰降得更多，從而出現(xiàn)兩個峰情況。此時粒子速度分布函數(shù)已經(jīng)遠(yuǎn)遠(yuǎn)偏離Maxwell分布了。

因此，通過上述的數(shù)值模擬的結(jié)果和分析，從物理機(jī)理出發(fā)提出的參數(shù)Zh可以很好地表征剪切流中的局部非平衡效應(yīng)。在Zh比較大的情況下，粒子速度分布函數(shù)已經(jīng)嚴(yán)重偏離Maxwell分布。而目前被廣泛采用的Boltzmann-BGK方程建立的依據(jù)是粒子分布函數(shù)偏離Maxwell分布不太遠(yuǎn)，并用簡單的向Maxwell分布趨近的弛豫過程代替分子間的碰撞作用。因此，雖然基于Boltzmann-BGK方程的各種動理學(xué)方法簡化了Boltzmann方程的計算，有效降低了計算成本，但用以模擬高速稀薄剪切流時則未必能給出正確結(jié)果。

2.2 剪應(yīng)力

從連續(xù)介質(zhì)觀點出發(fā)，由于黏性作用，一層流體對相對運動的另一層流體產(chǎn)生阻力，這種阻力稱為黏性切應(yīng)力，并且建立了黏性切應(yīng)力和切變率之間的線性關(guān)系，即牛頓黏性假設(shè)。而從微觀角度看，氣體剪應(yīng)力是分子動量交換的宏觀表現(xiàn)，快速層中的動量較大的分子躍入慢層，給慢層增加了動量，慢速層中的動量較小的分子躍入快速層，減少了快速層中的動量。在DSMC方法中，我們可以通過統(tǒng)計每一個網(wǎng)格上進(jìn)出分子的動量改變量，因而獲得相應(yīng)的黏性剪應(yīng)力。另一方面，利用在DSMC計算中所得的宏觀速度剖面，可以依據(jù)牛頓黏性公式計算出另一個剪應(yīng)力。我們預(yù)期，當(dāng)氣體稀薄效應(yīng)較低時，兩者應(yīng)基本相符，而隨著氣體稀薄效應(yīng)增強(qiáng)，二者的差別將擴(kuò)大。

以DSMC統(tǒng)計所得的剪應(yīng)力為準(zhǔn)，定義二者的相對誤差為Δτ=(μdu/dx-τDSMC)/τDSMC，則Δτ和Zh的關(guān)系如圖6所示?？梢钥闯?，由牛頓黏性假設(shè)計算的誤差隨Zh的增加而單調(diào)增加。圖中不同類型的標(biāo)記對應(yīng)于不同的壁面速度。從圖中可以看出，當(dāng)Zh小于0.2時，偏差在10%以內(nèi)，隨著Zh的增加偏差急劇增大，當(dāng)Zh約為0.67時，誤差已達(dá)到100%。

值得注意的是，由于內(nèi)摩擦作用，流場中的溫度并不均勻，中心點溫度最高，沿壁面法向成拋物線分布。圖6中不同Zh情況下中心點的溫度并不相同，也不是隨Zh單調(diào)增高或降低。例如，對應(yīng)于Zh=0.05的點，它的溫度是7.6T0，但剪應(yīng)力偏差僅為3.8%。而對應(yīng)于Zh=0.51的點，其溫度為1.5T0，但剪應(yīng)力的偏差為48%。所以圖6中所得的結(jié)果對溫度不具有依賴性。

圖6 DSMC統(tǒng)計所得到剪應(yīng)力和牛頓假設(shè)計算的剪應(yīng)力Fig.6 Deviation of shear stress calculated by Newton hypothesis

也可以把圖6中的結(jié)果，表達(dá)成兩種剪應(yīng)力之比Ae=τDSMC/(μdu/dx)和Zh的關(guān)系，如圖7所示。

圖7 黏性修正系數(shù)隨Zh的變化曲線Fig.7 Correctional coefficients of viscosity vary with parameter Zh

由以上結(jié)果我們分析可知，牛頓黏性公式在Zh較大的時候不再成立，而且所給出的剪應(yīng)力之值總是大于由DSMC給出之值。如果可以認(rèn)為DSMC所給之值更可信，則在CFD中，可將剪應(yīng)力寫成如下形式：

其中μe為等效黏性系數(shù)，μ為氣體的原本黏性系數(shù)，此黏性系數(shù)仍與溫度有關(guān)。Ae為黏性修正系數(shù)，其值隨Zh的變化已示于圖7中。

2.3 高超聲速小鈍度鈍楔模型繞流

為了檢驗在飛行器飛行的實際流場中是否存在Zh比較大的情況。我們采用CFD模擬了高超聲速小鈍頭楔的繞流問題。在此算例中，楔頭半徑R0=1 mm、半楔角5°，鈍楔總長3 m。來流馬赫數(shù)為15，迎角為15°。氣體性質(zhì)取70 km高空空氣之值，黏性系數(shù)與溫度的變化采用蘇士蘭(Sutherland)公式。 壁面溫度固定為1000 K，采用無滑移邊界條件。

圖8給出了整個流場的Zh分布的云圖。由圖可以看出，在背風(fēng)面的絕大部分地方Zh值都超過0.2。在流場中沿壁面法向搜索Zh的最大值，結(jié)果發(fā)現(xiàn)都是在壁面處。圖9給出了Zh沿壁面流向的分布，可以看出，背風(fēng)面的Zh值明顯大于迎風(fēng)面的值。這是因為背風(fēng)面的壓力明顯低于迎風(fēng)面的壓力，相應(yīng)地背風(fēng)面的分子自由程要顯著大于迎風(fēng)面的分子自由程。無論是背風(fēng)面還是迎風(fēng)面，Zh的最大值都出現(xiàn)在楔的頭部附近。

圖8 背風(fēng)與迎風(fēng)面流場的Zh值分布Fig.8 Zh distribution on the leeward side and windward side

圖9 背風(fēng)與迎風(fēng)壁面處的Zh值沿流向的分布Fig.9 Zh distribution on the leeward side and windward side

Zh值在前緣處迅速減小，主要是由于在前緣處邊界層厚度急速增加，速度梯度急劇減小。在背風(fēng)面，Zh的值從x=133 mm時的0.5緩慢地降為x=3000 mm時的0.24。在這一段中，本文的結(jié)果應(yīng)該都有效。但在更上游部分，Zh沿流向的變化很快，本文的在近似于宏觀定常的情況下所得結(jié)果，未必可以直接應(yīng)用，需要另做研究。

3 討 論

本文采用DSMC方法分析了高速剪切流中局部的非平衡效應(yīng)，提出了一個新的描述參數(shù)Zh，即一個平均自由程上的速度差與當(dāng)?shù)氐姆肿舆\動最可幾速度的比值。通過數(shù)值分析我們得到：

1)Zh參數(shù)能有效地刻畫氣體速度分布函數(shù)偏離Maxwell的程度，以Zh的值可以將分布?xì)w為三類：單峰的正態(tài)分布(Zh<1.0)，三峰(1.01.5)。

2) 由基于連續(xù)介質(zhì)的牛頓黏性假設(shè)所計算的剪應(yīng)力相對于由DSMC所得之值的差隨著Zh單調(diào)增加。

3) 氣體在一定稀薄程度、強(qiáng)剪切情況下呈現(xiàn)非牛頓現(xiàn)象，本文提出一個等效的黏性系數(shù)，并給出了修正系數(shù)。

4) 因為2.3節(jié)中的情況，其參數(shù)在實際的近空間飛行器的范圍之內(nèi)，因此本文的結(jié)果對設(shè)計飛行器是有實際意義的。

但本文終究只是研究了一種情況，即氣體稀薄效應(yīng)對剪切力的影響，而且是在流動演化是非常緩慢的條件下。不難想象，要解決真正的實際問題，還有不少問題需要研究。例如，氣體稀薄效應(yīng)對熱流的影響，這時參數(shù)Zh是否仍是唯一需要考慮的參數(shù)?反過來，溫度梯度如果很大，是否也會對剪切力和熱流產(chǎn)生影響?決定其影響是否需要另一個無量綱參數(shù)？本文中的結(jié)論雖然和溫度值無關(guān)，但如果溫度或溫度梯度很高，這個結(jié)論是否仍然正確？對應(yīng)2.3節(jié)中接近于楔的頭部處，Zh值很大，而且沿流向的變化非常快，本文在流動變化緩慢條件下所得結(jié)論顯然不一定能直接應(yīng)用，而必須考慮快速變化或非定常演化的因素，等等。我們將在今后的工作中陸續(xù)針對上述諸問題繼續(xù)工作。希望最終可以得到能解決在高超聲速飛行器設(shè)計中遇到的空氣動力學(xué)新問題的結(jié)果。

致謝:感謝周恒院士給予的直接指導(dǎo)幫助，以及張涵信院士在和周恒院士討論中給出的間接幫助和鼓勵。感謝天津大學(xué)青年人才自主科研基金支持。

[1]Zhou H, Zhang H X. New problems of aerodynamics[J]. Sci Sin-Phys Mech Astron, 2015, 45: 104709.(in Chinese)周恒, 張涵信. 空氣動力學(xué)的新問題[J]. 中國科學(xué): 物理學(xué) 力學(xué) 天文學(xué), 2015, 45: 104709

[2]Zhang S H, Li Q, Zhang L P, et al. The history of CFD in China[J]. Acta Aerodynamica Sinica, 2016, 34(2): 157-174.(in Chinese)張樹海, 李沁, 張來平, 等. 中國CFD史(英文)[J]. 空氣動力學(xué)學(xué)報, 2016, 34(2): 157-174

[3]Zhang H X. Investigations on fidelity of high order accurate numerical simulation for computational fluid dynamics[J]. Acta Aerodynamica Sinica, 2016, 34(1): 1-4. 張涵信. 關(guān)于CFD高精度保真的數(shù)值模擬研究[J]. 空氣動力學(xué)學(xué)報, 2016, 34(1): 1-4.(in Chinese)[4]Tsien H S. Superarodynamics, mechanics of rarefied gas[J]. J Aerospace Sci, 1946, 13: 342

[5]Kandlikar S G, Garimella S, Li D, et al. Heat transfer and fluid flow in minichannels and microchannels[M]. Second Edition, Elsevier, 2013

[6]Shen Q. Rarefied gas dynamics[M]. Beijing: National Defense Industry Press, 2003.(in Chinese)沈青, 稀薄氣體動力學(xué)[M]. 北京: 國防工業(yè)出版社, 2003

[7]Bird G A. Breakdown of translational and rotational equilibrium in gaseous expansions[J]. AIAA J, 1970, 8(11): 1998-2003

[8]Boyd I D, Chen G, Candler G V. Predicting failure of the continuum fluid equations in transitional hypersonic flows[J]. Phys Fluids, 1995, 7(1): 210-219

[9]Wang W L, Boyd I D. Predicting continuum breakdown in hypersonic viscous flows[J]. Phys Fluids, 2003, 15(1): 91-100

[10]Sun Q H, Boyd I D, Candler G V. A hybrid continuum/particle approach for modeling subsonic, rarefied gas flows[J]. J Comput Phys, 2004, 194(1): 256-277

[11]Li Z H, Li Z H, Li H Y, et al. Research on CFD/DSMC hybrid numerical method in rarefied flows[J]. Acta Aerodynamica Sinica, 2015, 33(2): 266-271.(in Chinese)李中華, 李志輝, 李海燕, 等. 過渡流區(qū)N-S/DSMC耦合計算研究[J]. 空氣動力學(xué)學(xué)報, 2015, 33(2): 266-271

[12]Lockerby D A, Reese J M, Struchtrup H. Switching criteria for hybrid rarefied gas flow solvers[J]. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science, 2009, 465(2105): 1581

[13]Grad H. On the kinetic theory of rarefied gases[J]. Communications on Pure and Applied Mathematics, 1949, 2(4): 331-407

[14]Meng J P, Dongari N, Reese J M, et al. Breakdown parameter for kinetic modeling of multiscale gas flows[J]. Phys Rev E, 2014, 89(6): 63306

[15]Meng J P, Zhang Y H, Shan X W. Multiscale lattice Boltzmann approach to modeling gas flows[J]. Phys Rev E, 2011, 83(4): 046701

[16]Eu B C. Kinetic theory and irreversible thermodynamics[M]. Canada: John Wiley & Sons, Inc., 1992

[17]Camberos J, Schrock C, McMullan R, et al. Development of continuum onset criteria with direct simulation Monte-Carlo using Boltzmann’s h-theorem: review and vision[C]//9th AIAA/ASME Joint Thermophysics and Heat Transfer Conference, AIAA, 2006: 2006-2942

[18]Xiao H, Myong R S, Singh S. A new near-equilibrium breakdown parameter based on the Rayleigh-Onsager dissipation function[J]. AIP Conference Proceedings, 2014, 1628: 527

[19]Bird G A. Moleculargas dynamics and the direct simulation of gas flows[M]. Oxford: Clarendon Press. 1994

[20]Bird G A. Recent advances and current challenges for DSMC[J]. Computer Math Applic, 1998, 35: 1-14

[21]Shen C. Rarefiedgas dynamics: Fundamentals, simulations and Micro-flows[M]. Springer, Berlin, 2005

[22]Fan J. Rarefied gas dynamics: Advances and applications[J]. Advances in Mechanics, 2013, 43(2): 185-201.(in Chinese)樊菁. 稀薄氣體動力學(xué)[J]. 進(jìn)展與應(yīng)用, 力學(xué)進(jìn)展, 2013, 43(2): 185-201

[23]Fan J, Shen C. Statistical simulation of low-speed rarefied gas flows[J]. J Comput Phys, 2001, 167: 393-412

[24]Wu L, White C, Scanlon T C, et al. Deterministic numerical solutions of the Boltzmann equation using the fast spectral method[J]. J Comput Phys, 2013, 250(1): 27-52

[25]Xu K, Huang J C. A unified gas-kinetic scheme for continuum and rarefied flows[J]. J Comput Phys, 2010, 229(20): 7747-7764

[26]Li Z H, Zhang H X. Study on gas kinetic unified algorithm for flows from rarefied transition to continuum[J]. J Comput Phys, 2004, 193(2): 708-738

[27]Li Z H, Zhang H X. Study on gas kinetic unified algorithm for flows from rarefied transition to continuum[J]. Acta Aerodynamica Sinica, 2003, 21(3): 255-266.(in Chinese)李志輝, 張涵信. 稀薄流到連續(xù)流的氣體運動論統(tǒng)一算法研究[J]. 空氣動力學(xué)學(xué)報, 2003, 21(3): 255-266

[28]Xu K. Direct modeling for computational fluid dynamics: construction and application of unified gas-kinetic schemes[M]. World Scientific Publishing Co, 2015

[29]Chapman S, Cowling T. The mathematical theory of non-uniform gases[M]. London: Cambridge Univ. Press, 1970

[30]Burnett D. The distribution of molecular velocities and the mean motion in a non-uniform gas[J]. Proceedings of the London Mathematical Society, 1936, 2(1): 382-435

[31]Gu X J, Emerson D R. A high-order moment approach for capturing non-equilibrium phenomena in the transition regime[J]. J Fluid Mech, 2009, 636: 177

[32]Myong R S. Thermodynamically consistent hydrodynamic computational models for high-Knudsen-number gas flow[J]. Physic of Fluids, 1999, 11(9): 2788-2802

[33]Xiao H, Shang Y H, Wu D, et al. Nonlinear coupled constitutive relations and its validation for rarefied gas flow[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(7): 2091-2104.(in Chinese)肖洪, 商雨禾, 吳迪, 等. 稀薄氣體動力學(xué)的非線性耦合本構(gòu)方程理論及驗證[J]. 航空學(xué)報, 2015, 36(7): 2091-2104.