歐吉輝, 趙 磊,2, 陳 杰,*

(1. 天津大學機械工程學院高速空氣動力學研究室, 天津市現代工程力學重點實驗室, 天津 300350；2. 中國空氣動力研究與發(fā)展中心 超高速空氣動力研究所, 綿陽 621000)

0 引 言

對于近空間高超聲速飛行器，其飛行高度一般在40～70km高空，飛行速度一般在5～20倍聲速范圍內，隨著飛行高度的增加，基于Navier-Stokes (N-S)方程的傳統(tǒng)CFD方法逐漸失效，需要研究新的空氣動力學問題[2]。近空間飛行的流場特點是，其邊界層中溫度很高，除迎風面壓力較高外，其它部位的分子自由程會顯著增大, 同時邊界層中的剪切很強，因而會出現文獻[1]中所考慮的那種稀薄氣體效應。

對于存在局部稀薄效應的近空間高超聲速飛行器流場模擬，現有的方法都存在一些局限，缺少模型可靠且計算量滿足工程實際需求的方法[3-4]。Bird[5]提出的直接模擬蒙特卡洛方法(DSMC)被認為是稀薄氣體流體動模擬方面最成功的方法之一，目前該方法已經可以用于模擬航天飛行器這類復雜三維外形的流場[6]，并且可以考慮分子內能激發(fā)、電離、離解/復合化學反應、熱輻射等復雜的物理化學變化[7-8]。但是，該方法將分子移動和碰撞解耦處理，要求計算網格小于分子的平均自由程、時間步長小于分子的平均碰撞時間，因此用于近連續(xù)流動時對計算資源提出了很高的要求。對于飛行高度在80km及以下的飛行器，計算量巨大，目前還無法應用于飛行器的設計上。Boltzmann方程雖然可以描述整個流域，但該方程是關于速度分布函數的七維積分微分方程，求解十分困難[9]。介于宏觀與微觀之間介觀層次的格子Boltzmann方法(lattice Boltzmann method, LBM)被廣泛用于模擬微尺度流動以及不可壓縮流動[10-11]。作為LBM的一個變體，離散玻爾茲曼方法(discrete Boltzmann method, DBM)已經可以用于模擬可壓縮化學反應流動[12]。不過這種方法在高超聲速流動方面的工程應用前景有待進一步研究。最近我國學者在跨流域統(tǒng)一算法方面取得了重要進展[13-14]，該算法將每個時間步長內宏觀量的更新和微觀氣體分布函數的更新緊密地耦合在一起。但對于高速流動，離散速度對內存的需求以及計算成本高，限制了其在工程中的應用。另外發(fā)展起來的NS-DSMC耦和模型需要對流場分區(qū)，不同區(qū)域采用不同的計算模型和算法。但分區(qū)的判據帶有一定的不確定性，分區(qū)計算界面處的數據交換也會限制計算效率的提升，已有文獻報道NS-DSMC相比全DSMC的計算速度提升可達3倍左右[15-16]，但這還不足以解決實際應用中DSMC和N-S在計算量上數量級的差別。

在CFD計算中，通過采用滑移邊界條件可以將N-S方程的應用范圍拓展到滑流區(qū)域[17-19]。但是隨著稀薄程度進一步增加，N-S方程的線性本構關系逐漸失效，僅使用滑移邊界條件不足以給出滿意的結果。一方面，Lockerby[20]等基于均勻剪切流的線化Boltzmann方程的解提出了一個壁面函數，從而對努森層內的黏性系數進行修正，來考慮壁面附近努森層的影響。Lofthouse[17]等直接將該壁面函數納入CFD滑移邊界條件，并用于高超聲速圓柱繞流計算，結果顯示該方法取得的效果還比較有限。另一方面，在微槽道流的計算中，通過考慮壁面對分子運動的限制效果從而得到等效黏性系數以對傳統(tǒng)的N-S本構關系進行修正[21-22]。然而，這種基于等效平均自由程的等效黏性系數的推導過程與幾何模型密切相關，并且在推導過程中沒有考慮空間密度的不均勻性，還不能直接用于高超聲速流動計算。

陳杰和趙磊在文獻[1]中，基于稀薄高超聲速邊界層流動特點，給出了一個判別氣體稀薄效應大小的參數Zh，并采用DSMC對簡單的平板剪切流進行了模擬分析。結果表明：由傳統(tǒng)連續(xù)介質模型所計算的剪應力的誤差隨著Zh參數單調增大，繼而由DSMC計算的剪應力在線性本構關系下所導出的等效黏性系數與傳統(tǒng)的黏性系數之比隨著Zh參數單調減小。據此，文獻[1]提出，在流場存在局部稀薄效應的情況下可以根據流場的每一點的Zh值，對連續(xù)介質模型中的黏性系數加以修正，應該是一個既考慮了流場的局部氣體稀薄效應，而又不顯著增加計算的復雜性的辦法。本文將以一個在70 km高空，以馬赫數15飛行的小迎角平板為例，來檢驗上述方法是否可行。

1 計算模型

1.1 CFD模型

在本文中計算采用的方程為二維守恒型可壓縮N-S方程，本構關系從形式上是線性關系，但在計算過程中，需要根據DSMC的結果對黏性系數進行修正。下面將這種修正N-S方程黏性系數的方法稱為NS-VC (Navier-Stokes-Viscosity Correction) 方法，而黏性系數沒有修正的則稱NS方法。二維無量綱可壓縮N-S方程為：

剪應力與熱流計算式為：

其中μ、κ分別為黏性系數與熱傳導系數，計算式為：

其中R*為氣體常數，對于空氣R*=287 J/(kg·K)。

陳杰等在文獻[1]中通過DSMC對簡單平板間均勻剪切流分析得到，由DSMC結果獲得的等效黏性系數與傳統(tǒng)連續(xù)介質模型中的黏性系數在連續(xù)流基本一致，但隨著稀薄效應參數Zh的增大，二者之比會單調減小。兩者之比，也是式(3)中黏性修正系數AZh，其與Zh參數的關系如圖1所示。圖中小黑點為文獻[1]中所給結果，實線為由三次樣條曲線擬合所得，以用于CFD計算。圖1中DSMC的數據包含了不同壁面馬赫數的工況，而不同工況槽道中心點的溫度差別很大，因此圖1的規(guī)律與溫度無關。從氣體的流變性質看，圖1的規(guī)律體現出氣體剪切變稀的性質，對于偏離平衡態(tài)比較遠的強剪切流動，氣體的等效黏性系數不再只與溫度有關，而是隨著剪切率的增加而逐漸減小。另外，雖然圖1的結果只顯示出了修正系數對Zh的依賴關系，但是Couette流同時也存在很大的溫度梯度與熱流，實際上圖1的結果也包含著熱流對黏性系數的影響。因此作為初步嘗試，可以將圖1的規(guī)律用于CFD計算看是否可以改進傳統(tǒng)CFD的結果。

圖1 黏性修正系數AZh隨Zh數的變化 Fig.1 The relationship between AZh and Zh

1.2 計算模型與工況

計算模型為鈍平板，模型示意圖如圖2所示，其中x、y分別為沿軸向以及垂直于軸向，坐標原點選在鈍頭圓心處。點A所示位置為鈍平板前緣，ξ為從前緣沿壁面的貼體弧長，AOA為迎角，Rn為有量綱圓頭半徑。在下文中，迎風面以及下平板指的是流場或平板對應y<0的區(qū)域，背風面以及上平板指的是流場或平板對應y>0的區(qū)域。本文計算了不同迎角、不同球頭半徑下馬赫數15、壁面溫度為2000K的工況，具體參數在表1列出。來流為70 km高空氣體，其物理參數為：

圖2 計算模型示意圖Fig.2 Sketch map of computation model

球頭半徑R?n/mm迎角AOA /(°)Case 145Case 247.5Case 3410Case 4205Case 52005

1.3 數值方法

由于NS-VC方法只對N-S方程的黏性系數進行了修正，其方程形式與N-S方程并無本質區(qū)別，因此傳統(tǒng)的CFD求解方法均有效?？紤]到二階NND[24](Non-oscillatory, No free parameter and Dissipative)格式具有無振蕩捕捉激波、魯棒性很好、計算效率高等優(yōu)點，在本文計算中，對流項離散采用二階NND格式，黏性項離散采用二階中心差分，時間推進采用二階Runge-Kutta法。為突出稀薄效應對黏性系數的影響，壁面沒有采用滑移邊界條件，而仍用無滑移等溫條件。遠場為來流條件，出口為外推邊界條件。

對于高超聲速稀薄流動，對黏性項的刻畫至關重要。由此，在圖3中給出了Case 1工況二階NND格式與五階WENO[25](Weighed Essentially Non-Oscillatory)格式(黏性項采用六階中心差分)計算得到的壁面摩擦系數與速度剖面比較的結果。其中ξ為沿壁面弧長，壁面摩擦系數計算式為:

由圖3可以看到，對于NND的計算結果，其背風面壁面摩擦系數和x*=1 m處速度剖面與五階WENO的結果都吻合得很好。因此，在本文的計算工況中采用二階NND格式可以很好的刻畫黏性項，后面的計算都采用二階NND格式進行計算。

(a) 背風面壁面摩擦系數沿壁面的分布

(b) 背風面x*=1 m處速度剖面

在CFD中實施黏性修正的具體計算過程為：先用NS方法計算一段時間(不需要等到收斂)，然后切換到NS-VC方法，在每一時間步，根據當時的流場，計算每一點的Zh值，得到相應的AZh，據此修正每一點的黏性系數，沿時間推進直至收斂。

實際上，NS方法和NS-VC方法的差別僅在于后者在每一步迭代后，需對所得流場的每一網格點計算Zh值，并據之對該處的黏性系數做修正，因此每一迭代步所需時間要增加約17%。此外，在啟動計算時，由于整個流場是從靜止態(tài)開始，而邊界條件則是給定值，所以在局部地方的剪切率會非常大，從而相應地該處的Zh值也很大，超出了可知的范圍而無法運行。因此要先用NS方法計算一段時間后，才能切換到NS-VC方法上來。這樣，總的迭代步數也會增加一些。對我們的算例，所需總的計算時間比單純用NS方法計算時約多50%左右。

2 結果及分析

2.1 黏性修正的影響

圖4給出的是用NS方法計算得到的Zh值在空間分布的等值線圖?？梢钥吹?，Zh值比較大的地方主要集中在激波、頭部和背風面壁面附近，說明這些位置是流動中稀薄效應最強的區(qū)域。圖5(a)給出的是NS方法計算得到的Zh值在背風面和迎風面沿壁面的分布。可以看到Zh值都是先增大后減小，在平板圓頭與板身相接的位置取得最大值。圖5(b)給出的是NS-VC計算收斂時背風面與迎風面壁面修正系數沿流向的分布。與圖5(a)中Zh值變化趨勢相反，黏性修正系數先減小后增大，在平板圓頭與板身相接處最小。這是由于Zh值越大，表征稀薄效應越強，因此黏性修正得越多。對于ξ>100，背風面黏性修正系數介于0.7與0.8之間，即修正后的等效黏性系數為原來黏性系數的70%～80%。迎風面的Zh值比背風面的小很多，因此對黏性系數的修正量也很小。

圖4 NS: Zh值空間分布的等值線圖Fig.4 Contour map of Zh value in spatial distribution

(a) NS：背風面與迎風面壁面Zh值沿壁面的分布

(b) NS-VC：壁面修正系數沿壁面的分布

圖6給出了有黏性修正與無修正情況下壁面壓力、壁面摩擦系數以及相對變化率(σDev)沿壁面的變化，其中相對變化率是(NS-VC結果-NS結果)/(NS結果)。從圖中可以看到黏性修正使得壁面壓力與壁面摩擦系數均減小，其中對背風面的影響比迎風面大很多。由圖6(a)可以看到，黏性修正使得背風面的壁面壓力相對于無修正來說在大部分地區(qū)均約減小6%左右。而圖6(b)可看出，黏性修正使得背風面的壁面摩擦系數相對于無修正來說減小最多為57%，而在ξ>100以后，則從11%逐步降至7%左右。即稀薄效應對背風面的壁面壓力以及摩擦力均會造成比較大的影響。但對迎風面來說，則稀薄效應的影響很小。

圖7給出了背風面與迎風面在x*=1 m處速度沿壁面法向分布的剖面。由圖7可以看到，黏性修正使得邊界層內的速度梯度增大，在背風面的修正效果比迎風面強很多。

2.2 迎角的影響

(a) 壁面壓力

(b) 壁面摩擦系數

(a) 背風面

(b) 迎風面

圖8給出的是迎角分別為5°、7.5°、10°情況下背風面壁面Zh值隨壁面弧長的變化，右上角為頭部的局部放大圖，圓頭與平板相切的位置大約在ξ=1.57 處?？梢钥吹讲煌堑谋诿鎆h值變化規(guī)律相似，都是先增大后減小，在頭部附近急劇變化，并在圓頭與板身相切的位置附近取得最大值。對比三個迎角的結果可以發(fā)現，迎角越大，其背風面壁面Zh值越大，表示其稀薄程度也越大。

圖8 不同迎角背風面壁面Zh值沿壁面的變化Fig.8 Variation of Zh value along the wall surface of the leeward with different angles of attack

表2給出的是不同迎角下有黏性修正和無修正情況下頭部的壓力和摩擦力產生的升力和阻力，表3給出的是整體的升力、阻力以及升阻比，表4給出的是板身部分上下板的壓力以及摩擦力分別對升力、阻力的貢獻，同時表2、表3、表4中都給出了黏性修正后的結果相對于無修正的變化率(σDev)。

由表2與表3對比可以看到，頭部的升力相對于總的升力來說很小，幾乎可以忽略，頭部的阻力和總的阻力相比不能忽略，黏性修正使得頭部阻力的減小大約在3%左右。由表4可以看到，當迎角由5°變到10°，黏性修正使得上、下板壓力與摩擦力均減小，其中上板摩擦力減小8.42%～13.13%，上板壓力減小4.31%～6.43%，下板摩擦力減小2.92%～1.60%，下板壓力減小0.85%～0.18%。因此，黏性修正主要對上板產生影響，而對下板影響比較小，并且迎角越大，上板的稀薄效應越強，黏性修正使得上板壓力、摩擦力減小越多。同時可以看到，升力主要由上、下板的壓力貢獻，而阻力由上、下板的摩擦力以及下板壓力貢獻。綜合頭部與板身的效果，由表3可以看到，黏性修正使得最終的升力增加0.5%左右，阻力減小3.4%～1.71%，升阻比增加4.13%～2.23%。因此，稀薄效應存在時， NS方法預測的升力會比實際的偏小，而阻力比實際的偏大，導致其預測得到的升阻比偏小。迎角越大，上板的稀薄效應越大，但從表3看到黏性修正后升阻比的增加反而越小。這是由于隨著迎角變大，下板壓力迅速增大，使得下板壓力對升力以及阻力的貢獻都迅速增大(見表4)，由此使得上板對升力阻力貢獻在總的升力阻力里面占的比重都越來越小，而黏性修正對下板影響很小，最終使得迎角越大升阻比變化越小。

表2 不同迎角下升力、阻力表(積分從x*=-4 mm到x*=0 mm，僅頭部)Table 2 Lift and drag with different angles of attack(Integral from x*=-4 mm to x*=0 mm, only the head)

2.3 球頭半徑的影響

圖9給出了在迎角為5°時頭部半徑分別為4 mm、20 mm、200 mm ( Case 1、Case 4、Case 5 ) 壁面Zh值隨壁面弧長的變化?？梢钥吹?，頭部半徑越大，Zh值越小，表明稀薄效應越小，其越趨于連續(xù)流。

圖9 不同頭部半徑壁面Zh值沿壁面的變化Fig.9 Variation of Zh value along wall surface with different head radii

表4 不同迎角上下板的壓力、摩擦力對升力、阻力貢獻表(單位：N/m)(積分從x*=0 mm到x*=3000 mm，僅平板部分)Table 4 The contribution of pressure and friction of upper and lower plate to lift and drag(Integral from x*=0 mm to x*=3000 mm, only the plate)

表5 不同圓頭半徑升力、阻力、升阻比表(積分從前緣到x*=3000 mm)Table 5 Lift, drag and lift-drag ratio with different head radii (Integral from leading edge to x*=3000 mm)

3 結論與討論

本文的主要目的是看文獻[1]中提出的通過等效黏性考慮局部稀薄效應的結果是否能方便地應用于空氣動力學計算中，也要看應用后得到的結果是否更符合實際。這兩點都已得到正面的結果。第一，的確可以方便地應用于空氣動力學技術中，而且所增加的時間非常有限，計算量的增加在50%以下。第二，所得結果和傳統(tǒng)的計算方法所得結果相比，壁面摩擦系數減小，升阻比增加，其改進的方向與現有飛行試驗結果定性相符。本文還進一步分析了迎角及球頭半徑的影響。

但是需要說明的是目前所做只是初步的嘗試，真正解決實際問題，還有以下問題需要進一步考慮：

1) 本文所用的黏性系數的修正系數與參數Zh的關系是針對單元子氣體氬氣獲得的，對真實的空氣來說，即使定性上相似，定量上也會有一定程度的不同。

2) 由于在文獻[1]中的DMSC計算中，黏性與溫度的依賴關系采用了硬球模型的根號關系，與傳統(tǒng)CFD中用的公式有一定的差別。在用于實際問題時，需要重新考慮。

3) 本文沒有考慮稀薄效應對CFD中壁面條件的影響。

4) 高超聲速流往往伴隨有高溫及熱傳導問題。這里稀薄效應對熱傳導的影響考慮的還不全面。因此，稀薄效應對熱傳導問題的影響，熱傳導對黏性的影響以及傳熱和黏性在稀薄效應下的相互影響也是今后需要考慮的問題。

致謝:感謝周恒院士給予的直接指導幫助。