黃 飛，張 亮，程曉麗，沈 清

（中國航天空氣動力技術(shù)研究院，北京 100074）

0 引 言

傳統(tǒng)航天飛行器在80km以上的高空大氣環(huán)境中高速飛行，其流動表現(xiàn)出稀薄效應(yīng)，而對于未來采用升力體或乘波體構(gòu)型［1］的近空間（20km～100km）高超聲速飛行器，由于其外形特征多為扁平狀升力面，具有尖化前緣、復雜氣動外形、長時間飛行和高升阻比等氣動構(gòu)型和飛行軌道特征，飛行環(huán)境介于稠密大氣層和稀薄大氣層之間，其局部小尺寸部件在中等高空，甚至在較低的高空就會出現(xiàn)稀薄氣體效應(yīng)。在稀薄流區(qū)，要實現(xiàn)飛行器的飛行、機動和控制，必須了解它們的氣動力特性，同時為飛行器防、隔熱設(shè)計提供氣動熱數(shù)據(jù)。而隨著飛行高度的變化，摩擦力、熱流等氣體動力學特性較之連續(xù)流區(qū)發(fā)生很大的變化。這些變化對于傳統(tǒng)的鈍頭體再入飛行器軌道影響較小，因此以往很少細致考慮稀薄氣體效應(yīng)的影響，大余量的設(shè)計方法對于解決這類外型和軌道的飛行器已足夠，而對大升力面、尖化前緣復雜外形中低空飛行的稀薄流問題還有待研究。因此，中低Kn數(shù)的尖化前緣外形的氣動力、氣動熱的正確預測對近空間空域飛行器的設(shè)計和工程應(yīng)用都具有重要意義。

在低空問題研究中，由于大氣密度高，常采用傳統(tǒng)的CFD方法進行流場的模擬分析，而飛行器飛行過程中，隨著飛行高度的變化，其所處的外部流動環(huán)境是一個漸變的過程，那么連續(xù)流假設(shè)在什么條件下失效？連續(xù)流計算方法對近空域飛行器氣動力、熱預測的偏差有多大？就成為目前急需解決的問題。國外已開展了許多關(guān)于連續(xù)流失效性判據(jù)［2－6］的研究工作，而對近空域稀薄效應(yīng)對飛行器氣動特性變化的定量分析還不多見［7－9］。

本文首先采用了圓柱算例對四種滑移模型的有效性進行了對比分析，最后針對未來近空域飛行的高超聲速飛行器大升力面的構(gòu)型特征，采用典型外形，選用適當?shù)幕颇Ｐ秃蜔o滑移NS方法對近空間飛行器過渡流效應(yīng)引起的氣動特性變化進行了對比分析，定量地給出了稀薄氣體效應(yīng)對氣動特性的影響規(guī)律。

1 計算方法

1.1 N－S方程

在笛卡兒坐標系下，采用層流三維非定常可壓縮NS方程：

其中，Q為守恒變量；F、G、H為三個方向的無粘矢通量；Fv、Gv、Hv為三個方向的粘性矢通量。

1.2 滑移邊界類型

本文采用了四種滑移邊界模型，分別定義如下：Type1為經(jīng)典的Maxwell滑移邊界類型；Type2為Gokcen［10］滑移邊界；Type3為 Lockerby［11］提出的壁函數(shù)修正法；HS為分子自由程采用硬球模型簡化計算時的Maxwell滑移邊界。

Type2（對應(yīng)文獻［9］中CFD（2））：

Type3（對應(yīng)文獻［9］中CFD（3））：

HS：

σ動量調(diào)節(jié)系數(shù)，α能量調(diào)節(jié)系數(shù)。

上述結(jié)果表明，在非高斯干擾明顯的環(huán)境中，原有算法性能會急劇下降甚至失效，但基于Alpha穩(wěn)定分布的新算法可以有效完成雜波抑制，為接下來的動目標檢測奠定基礎(chǔ)。

1.3 無量綱參數(shù)定義

偏差分析量定義如下：

其中，qnoslip為無滑移NS方程熱流結(jié)果，qslip2為滑移模型Type2的熱流計算結(jié)果。

2 滑移模型的適應(yīng)性分析

本文滑移模型的算例驗證采用文獻［9］中圓柱的計算條件，來流氣體為氬氣，來流馬赫數(shù)分別為Ma＝10、25，來流溫度T∞＝300K，圓柱直徑D＝304.8mm，基于直徑D＝2R的努森數(shù)Kn為0.002、0.05、0.25。本算例分別采用了 Type1、Type2、Type3、HS、No＿slip等方法對圓柱擾流進行了計算分析，并與文獻［9］中的DSMC結(jié)果進行了對比。此算例中的粘性系數(shù)計算方法見文獻［9］，對應(yīng)于DSMC的VHS模型。

圖1 滑移與非滑移下Ar的流場結(jié)構(gòu)（Ma＝10）Fig.1 Flow field with different methods（Ma＝10）

圖1（a）為Kn＝0.002時的流場壓力云圖，從中可以看出，滑移方法與非滑移方法在近連續(xù)流區(qū)所得到的流場結(jié)構(gòu)非常一致，滑移流方法能夠準確捕捉激波結(jié)構(gòu)和尾渦結(jié)構(gòu)。圖1（b）為Kn＝0.05時兩種方法所得的流場壓力云圖和速度型，從圖中可以發(fā)現(xiàn)由于在此努森數(shù)下存在較弱的稀薄效應(yīng)，兩種方法所得結(jié)果存在以下三點不同：滑移流方法所得的激波厚度比非滑移流方法稍厚；在物面處，滑移流方法所得的物面速度出現(xiàn)明顯的滑移速度，使得物面處速度型不同于非滑移流所得結(jié)果；在尾部區(qū)域，兩種方法所得的壓力等值線明顯存在不同。稀薄效應(yīng)引起的以上三點不同將很有可能影響物面的熱流、壓力等分布。為了進行定量分析下面給出了不同滑移模型、非滑移方法、DSMC方法的物面氣動特性分布。

圖2分別給出了馬赫數(shù)等于10時，Kn＝0.002、0.05、0.25時不同滑移模型與DSMC方法所得的熱流系數(shù)分布曲線，其中圖2（a）為本文計算結(jié)果，圖2（b）為文獻［9］計算結(jié)果。從圖中可以看出，在較弱稀薄效應(yīng)下，即努森數(shù)Kn＝0.002時，不同滑移模型所得熱流系數(shù)分布與DSMC結(jié)果吻合較好，當努森數(shù)Kn增加到0.05、0.25時不同模型所得結(jié)果之間的差距開始明顯顯現(xiàn)；隨著努森數(shù)的增加，無滑移NS方法所得熱流分布都一致的高于滑移方法及DSMC方法所得結(jié)果，并且從圖中可以發(fā)現(xiàn)，類型2的滑移模型即使在較大努森數(shù)Kn＝0.05、0.25時也表現(xiàn)出了與DSMC相一致的結(jié)果，其它模型所得結(jié)果在較大努森數(shù)下都一致高于DSMC方法所得結(jié)果，由此可見，類型2的滑移模型對較大努森數(shù)下具有更好的適應(yīng)性，其它模型的結(jié)果都將高估熱流而使得防熱設(shè)計趨于保守。通過與文獻結(jié)果對比可以發(fā)現(xiàn)，本文與文獻結(jié)果熱流分布吻合很好。

圖3分別給出了Kn＝0.002、0.05、0.25時本文與文獻［9］采用不同滑移模型與DSMC方法所得的壓力系數(shù)分布曲線，從圖中可以看出，當努森數(shù)Kn＝0.002、0.05時，不同滑移模型所得表面壓力系數(shù)分布與DSMC所得結(jié)果非常一致，當努森數(shù)增加至0.25時，各種模型所得壓力分布都一致高于DSMC方法所得結(jié)果，這主要是因為強烈的稀薄效應(yīng)使得NS方法所得激波結(jié)構(gòu)明顯薄于DSMC方法所得激波結(jié)構(gòu)，NS方法產(chǎn)生的激波壓縮性更強，使得波后壓力較高。從圖2熱流系數(shù)分布和圖3壓力系數(shù)分布對比結(jié)果可以看出，相對于熱流而言，壓力對稀薄效應(yīng)敏感性較弱。通過與文獻結(jié)果對比可以發(fā)現(xiàn)，本文與文獻結(jié)果壓力分布同樣吻合很好。

圖4、圖5分別給出了馬赫數(shù)等于25時，Kn＝0.002、0.05、0.25時本文采用不同滑移模型所得的熱流系數(shù)與壓力系數(shù)分布曲線與文獻［9］DSMC方法的結(jié)果比較，從圖中可以看出，馬赫數(shù)等于25時的熱流、壓力分布規(guī)律與馬赫數(shù)等于10時的分布規(guī)律不同之處僅在于量值上的不同，其它規(guī)律與馬赫數(shù)等于10時的結(jié)果相似，類型2的滑移模型在較大馬赫數(shù)下也同樣具有較好的適應(yīng)性。

綜合以上計算分析可以得出以下結(jié)論：通過本文計算結(jié)果與文獻［9］所得結(jié)果的對比可以發(fā)現(xiàn)，本文結(jié)果與文獻結(jié)果吻合較好，本文所發(fā)展的代碼具有一定的可靠性；滑移模型2的適應(yīng)范圍較其它模型更廣，它能夠在較大努森數(shù)下取得相對較為滿意的熱流結(jié)果；在較大努森數(shù)下，滑移模型所得熱流、壓力結(jié)果都一致高于DSMC結(jié)果，使得防熱設(shè)計趨于保守；相對于熱流而言，壓力對稀薄效應(yīng)的敏感性更弱；經(jīng)過數(shù)值模擬還可以發(fā)現(xiàn)，模型2雖然適應(yīng)性廣，但相對于其它模型而言，計算時較為耗時，收斂性較差，對網(wǎng)格要求較高，如果在飛行器防熱設(shè)計條件允許下，可以考慮采用其它幾種滑移模型。

圖2 不同滑移模型在不同努森數(shù)下的物面熱流系數(shù)分布（Ma＝10）Fig.2 Heat transfer rate on the surface at different Knnumbers（Ma＝10）

圖3 不同滑移模型在不同努森數(shù)下的物面壓力系數(shù)分布（Ma＝10）Fig.3 Pressure coefficient on the surface at different Knnumbers（Ma＝10）

圖4 不同滑移模型在不同努森數(shù)下的物面熱流系數(shù)分布Fig.4 Heat transfer rate on the surface at different Knnumbers（Ma＝25）

圖5 不同滑移模型在不同努森數(shù)下的物面壓力系數(shù)分布Fig.5 Pressure coefficient on the surface at different Knnumbers（Ma＝25）

3 梯形翼近空間氣動特性的計算分析

由于滑移模型2在稀薄氣動特性預測方面的適應(yīng)范圍更廣，預測結(jié)果更為接近DSMC結(jié)果，故本文針對馬赫數(shù)15，飛行高度h＝50km～80km的梯形翼主要采用了模型2進行了稀薄氣動特性的計算分析。飛行攻角為10°，壁溫為500K，翼前緣直徑為30mm，翼根弦長約1.96m，翼尖弦長約0.63m，翼展長為0.5m，前緣后掠角約20°。

圖6為本文所采用的計算網(wǎng)格示意圖，圖7給出了梯形翼不同位置處的截面示意圖，其中Z＝5mm為梯形翼對稱面部位，Z＝300mm為翼中部位。為了進行定量分析，以下給出了截面位置處不同方法所得結(jié)果的差異性比較分析。

圖8、圖9分別為攻角10°下截面Z＝5mm、Z＝300mm處不同方法所得的表面熱流分布及偏差曲線結(jié)果，從上圖熱流分布結(jié)果可以發(fā)現(xiàn)，在相對較低的飛行高度，滑移方法與無滑移方法所得結(jié)果較為一致，隨著飛行高度的增加，兩種方法所得結(jié)果的差距逐漸增大。從圖9偏差曲線圖的定量結(jié)果可以明顯看出，由于迎風面氣體壓縮，背風面氣體膨脹，使得滑移的影響首先出現(xiàn)在背風面，故無滑移方法與滑移方法之間的差距在迎風面明顯低于背風面，并且隨著飛行高度的增加，稀薄效應(yīng)越明顯，滑移方法與非滑移方法所得的熱流結(jié)果之間的偏差逐漸增大。從圖中還可發(fā)現(xiàn)，此偏差最大值主要發(fā)生在翼前緣背風面膨脹區(qū)附近，當此種飛行器從50km至80km飛行時，在飛行器翼面大面積上的偏差約在5%至15%之間變化。

表1為不同高度下滑移方法和無滑移方法所得氣動力特性的結(jié)果比較，從中可以看出，隨著飛行高度的增加摩阻系數(shù)、軸向力系數(shù)、法向力系數(shù)都有所增加，其中摩阻增加較大，故而可以發(fā)現(xiàn)，升阻比隨著高度的增加逐漸減小。從不同方法的偏差結(jié)果可以看出，隨著飛行高度的增加，氣體稀薄度增加，稀薄效應(yīng)增強，傳統(tǒng)無滑移假設(shè)下的NS方法已不能正確描述壁面邊界條件，故而兩種方法所得結(jié)果的偏差量逐漸增大。

圖6 計算網(wǎng)格示意圖Fig.6 Grid in the computation region

圖7 不同位置處的截面示意圖Fig.7 Slice at different stations

圖8 不同位置處的熱流分布結(jié)果Fig.8 Heat flux at different stations

圖9 偏差分析結(jié)果Fig.9 Error at different stations

表2為不同方法所得峰值熱流隨高度的變化，從中可明顯看出，隨著高度的增加，氣體密度降低，飛行器前緣的峰值熱流降低較快，而采用滑移與無滑移兩種方法所得的結(jié)果偏差則隨著高度的增加而增大，這主要是由于高度的增加使得稀薄度增大，傳統(tǒng)連續(xù)流假設(shè)下的無滑移方法不再有效。

表1 不同方法下的氣動力特性結(jié)果Table 1 Aerodynamics with different methods

表2 不同方法下的峰值熱流結(jié)果Table 2 Peak heat flux with different methods

4 結(jié) 論

高超聲速飛行器飛行過程中經(jīng)歷了不同大氣環(huán)境密度的變化，本文針對中低空域大升力面飛行的高超聲速飛行器的結(jié)構(gòu)特征，采用典型外形研究了連續(xù)流失效后不同滑移模型對表面氣動特性的預測能力，最后給出了大升力面梯形翼近空間氣動特性的變化規(guī)律。結(jié)果表明：

（1）本文所得規(guī)律與文獻結(jié)果吻合很好，所建立的模擬手段具有一定的可靠性。

（2）滑移模型2比其它模型適應(yīng)范圍更廣，能夠在較大努森數(shù)下取得相對較為滿意的結(jié)果。

（3）在較大努森數(shù)下，除了滑移模型2以外，其它模型所得熱流、壓力結(jié)果都一致高于DSMC結(jié)果，使得防熱設(shè)計趨于保守。

（4）相對于熱流而言，壓力對稀薄效應(yīng)的敏感性較弱。

（5）在本文所研究的尺寸下，當大升力面梯形翼的飛行高度從50km增加至80km時，滑移方法與無滑移方法在翼面大面積上所得的熱流偏差量約在5%至15%之間變化；而駐點處峰值熱流的偏差量約在1.6%至14.5%之間變化；相應(yīng)的升阻比偏差結(jié)果約在0.06%至14.38%之間變化。

［1］李怡勇，沈懷榮.發(fā)展近空間飛行器系統(tǒng)的關(guān)鍵技術(shù)［J］.裝備指揮技術(shù)學院學報，2006，17（5）.

［2］BIRD G A.Breakdown of translational and rotational equilibrium in gaseous expansions［J］.AIAAJournal，1970，8（11）：1998－2003.

［3］TIWARI S.Coupling of the boltzmann and Euler equations with automatic domain decomposition［J］.Journal ofComputationalPhysics，1998，144：710－726.

［4］CAMBEROS J A，SCHROCK C R，McMULLAN R J.Development of continuum onset criteria with direct simulation Monte－Carlo using boltzmann＇s H－theorem：review and vision［R］.Proceedings of the 9th AIAA／ASME Joint Thermophysics and Heat Transfer Conference，San Francisco，California，June 2006.

［5］BOYD I D，CHEN G，XANDLER G V.Predicting failure of the continuum fluid equations in transitional hypersonic flows［J］.PhysicsofFluids，1995，7（1）：210－219.

［6］WANG W L.A hybrid particle／continuum approach for nonequilibrium hypersonic flows［D］.［Ph.D.Thesis］.The University of Michigan，2004.

［7］BOYD I D，PADILLA J F.Simulation of sharp leading edge aerothermodynamics［R］.AIAA Paper 2003－7062.

［8］LOFTHOUSE A J，BOYD I D，etc.Effects of continuum breakdown on hypersonic aerothermodynamics［R］.AIAA Paper 2006－993，2006

［9］lOFTHOUSE A J，SCALABRIN L C，BOYD I D.Velocity slip and temperature jump in hypersonic aerothermodynamics［R］.AIAA Paper 2007－0208，2007.

［10］GOKCEN T，MACCORMACH R W.Nonequilibrium effects for hypersonic transitional flows using continuum approach［R］.AIAA Paper 1989－0461，1989.

［11］LOCKERBY D A，REESE J M，GALLIS M A.Capturing the knudsen layer in continuum－fluid models of nonequilibrium gas flows［J］.AIAAJournal，2005，43（6）：1391－1393.