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

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

0 引 言

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

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

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

1 計(jì)算方法

1.1 N－S方程

在笛卡兒坐標(biāo)系下，采用層流三維非定?？蓧嚎sNS方程：

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

1.2 滑移邊界類型

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

Type2（對(duì)應(yīng)文獻(xiàn)［9］中CFD（2））：

Type3（對(duì)應(yīng)文獻(xiàn)［9］中CFD（3））：

HS：

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

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

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

偏差分析量定義如下：

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

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

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

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

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

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

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

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

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

圖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 梯形翼近空間氣動(dòng)特性的計(jì)算分析

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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

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

［1］李怡勇，沈懷榮.發(fā)展近空間飛行器系統(tǒng)的關(guān)鍵技術(shù)［J］.裝備指揮技術(shù)學(xué)院學(xué)報(bào)，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.