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

(1. 天津大學(xué)機(jī)械工程學(xué)院高速空氣動(dòng)力學(xué)研究室, 天津市現(xiàn)代工程力學(xué)重點(diǎn)實(shí)驗(yàn)室, 天津 300350；2. 中國空氣動(dòng)力研究與發(fā)展中心 超高速空氣動(dòng)力研究所, 綿陽 621000)

0 引 言

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

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

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

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

1 計(jì)算模型

1.1 CFD模型

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

剪應(yīng)力與熱流計(jì)算式為：

其中μ、κ分別為黏性系數(shù)與熱傳導(dǎo)系數(shù)，計(jì)算式為：

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

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

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

1.2 計(jì)算模型與工況

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

圖2 計(jì)算模型示意圖Fig.2 Sketch map of computation model

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

1.3 數(shù)值方法

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

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

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

(a) 背風(fēng)面壁面摩擦系數(shù)沿壁面的分布

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

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

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

2 結(jié)果及分析

2.1 黏性修正的影響

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

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

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

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

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

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

2.2 迎角的影響

(a) 壁面壓力

(b) 壁面摩擦系數(shù)

(a) 背風(fēng)面

(b) 迎風(fēng)面

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

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

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

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

表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°時(shí)頭部半徑分別為4 mm、20 mm、200 mm ( Case 1、Case 4、Case 5 ) 壁面Zh值隨壁面弧長的變化。可以看到，頭部半徑越大，Zh值越小，表明稀薄效應(yīng)越小，其越趨于連續(xù)流。

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

表4 不同迎角上下板的壓力、摩擦力對升力、阻力貢獻(xiàn)表(單位：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 結(jié)論與討論

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

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

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

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

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

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

致謝:感謝周恒院士給予的直接指導(dǎo)幫助。