郭晨林，陳方，趙艷彬，尤超藍(lán)，錢(qián)勇，李文龍，楊麗麗

（1.上海交通大學(xué) 航空航天學(xué)院，上海 200240；2.上海衛(wèi)星工程研究所研發(fā)中心，上海 201109）

0 引言

超低軌飛行器由于其重訪周期短、對(duì)地成像觀測(cè)性能強(qiáng)的特點(diǎn)，在通信、對(duì)地探測(cè)等領(lǐng)域展現(xiàn)出了巨大的潛力［1］，但超低軌飛行器受到的氣動(dòng)攝動(dòng)大，造成軌道維持與姿態(tài)控制的困難。同時(shí)，由于其運(yùn)行高度在120～300 km 之間，屬于過(guò)渡流區(qū)與自由分子流區(qū)，對(duì)于氣動(dòng)力的估算缺少經(jīng)驗(yàn)理論與工具，造成了總體設(shè)計(jì)與動(dòng)力選型上的困難。因此，當(dāng)前急需有關(guān)超低軌衛(wèi)星的快速氣動(dòng)力計(jì)算工具，以實(shí)現(xiàn)合理的動(dòng)力選型，延長(zhǎng)衛(wèi)星工作壽命［2］與加速超低軌飛行器的初步設(shè)計(jì)迭代。

對(duì)于工作在自由分子流區(qū)中的飛行器的氣動(dòng)力特性求解，BIRD［3］提出的直接模擬蒙特卡洛法（Direct Simulation Monte Carlo，DSMC）求解分子之間的碰撞效應(yīng)高精度求解其氣動(dòng)力特性，DAVIS［4］提出的測(cè)試粒子蒙特卡洛法（Test Particle Monte Carlo，TPMC）省略了分子間的碰撞效應(yīng)，雖然降低求解精度，但是大幅減小了計(jì)算量，提升了計(jì)算效率。Maxwell 通過(guò)建立簡(jiǎn)單的面元與氣體粒子的碰撞模型，求解自由分子流區(qū)中的面元受力狀況，進(jìn)而積分求解飛行器受力［5］。該方法成為自由分子流區(qū)中快速氣動(dòng)力計(jì)算的基本方法，后續(xù)許多學(xué)者根據(jù)分子與面元碰撞時(shí)的動(dòng)量損失與漫反射效應(yīng)對(duì)Maxwell 提出的模型進(jìn)行優(yōu)化，擴(kuò)展了該模型的精度與適用范圍［6-7］。

在過(guò)渡流區(qū)中，氣體的運(yùn)動(dòng)規(guī)律更為復(fù)雜，難以抽象出簡(jiǎn)化的物理模型。當(dāng)前求解過(guò)渡流區(qū)中的流動(dòng)特性多依賴于DSMC 方法［8］。由于在該區(qū)域中大氣密度較高，意味著DSMC 模擬的分子數(shù)量多，所需要的計(jì)算量大，因此，求解過(guò)渡流區(qū)氣動(dòng)特性的方法為：首先對(duì)連續(xù)流區(qū)的氣動(dòng)特性與自由分子流區(qū)中的氣動(dòng)特性進(jìn)行求解，再使用橋函數(shù)加權(quán)來(lái)近似過(guò)渡流區(qū)的氣動(dòng)特性。常用的橋函數(shù)有sine-square 橋函數(shù)［9］與erf-log 橋函數(shù)［10］，通過(guò)假定多個(gè)待定系數(shù)函數(shù)擬合連續(xù)流區(qū)-過(guò)渡流區(qū)-自由分子流區(qū)的部分計(jì)算結(jié)果，實(shí)現(xiàn)對(duì)其他工況的結(jié)果進(jìn)行估計(jì)。以上橋函數(shù)的擬合只對(duì)特定的氣流來(lái)流條件與當(dāng)前幾何外形的飛行器適用，當(dāng)外部條件發(fā)生變化時(shí)，需要新的DSMC 或試驗(yàn)結(jié)果進(jìn)行標(biāo)定。

目前，過(guò)渡流區(qū)的計(jì)算高度依賴于DSMC，TPMC 等高精度計(jì)算工具，快速氣動(dòng)力計(jì)算工具十分有限，常用的氣動(dòng)估算方法存在一定的局限性。同時(shí)，面元法應(yīng)用的前景廣泛，通過(guò)對(duì)面元法進(jìn)行進(jìn)一步改進(jìn)，以上問(wèn)題能被較好的解決。本文旨在解決超低軌飛行器的快速氣動(dòng)計(jì)算需求，提出一套適用于超低軌衛(wèi)星從過(guò)渡流區(qū)到自由流區(qū)的快速氣動(dòng)估算方法。在面元法的基礎(chǔ)上，對(duì)迎風(fēng)單元篩選算法與橋函數(shù)進(jìn)行了設(shè)計(jì)，達(dá)到了提高精度，增加適用范圍的效果。

1 超低軌衛(wèi)星快速氣動(dòng)力計(jì)算方法

本文提出的超低軌衛(wèi)星快速氣動(dòng)力計(jì)算方法主要分為3 個(gè)部分，如圖1 所示。包括輸入數(shù)據(jù)、數(shù)據(jù)的預(yù)處理方式和氣動(dòng)力計(jì)算方法。其中，有2 個(gè)過(guò)程顯著影響計(jì)算精度和計(jì)算效率：①數(shù)據(jù)預(yù)處理與迎風(fēng)單元與背風(fēng)單元的劃分算法；②氣動(dòng)力的計(jì)算方法和橋函數(shù)模型。本文將對(duì)此進(jìn)行較詳細(xì)的闡述。

圖1 超低軌衛(wèi)星快速氣動(dòng)力計(jì)算方法流程Fig.1 Flow chart of the rapid aerodynamic calculation method for ultra-low orbit spacecrafts

1.1 數(shù)據(jù)預(yù)處理

在輸入的模型STL 格式文件中，包含著衛(wèi)星的幾何尺寸信息，首先需要對(duì)其表面單元進(jìn)行離散，得到其表面的三角形單元如圖2 所示。在建模的過(guò)程中，常以飛行器的體坐標(biāo)系進(jìn)行建模，為方便后續(xù)的迎風(fēng)背風(fēng)單元算法及氣動(dòng)力計(jì)算，需要對(duì)離散單元在風(fēng)軸系和體軸系下進(jìn)行自由轉(zhuǎn)換。通過(guò)以下旋轉(zhuǎn)矩陣可以實(shí)現(xiàn)轉(zhuǎn)換過(guò)程：

圖2 迎風(fēng)單元與背風(fēng)單元篩選Fig.2 Screening diagram of windward and leeward units

式中：Cbo為體坐標(biāo)系ObXbYbZb至風(fēng)軸坐標(biāo)系OwXwYwZw下的旋轉(zhuǎn)矩陣；α、β分別為當(dāng)前飛行器的迎角、側(cè)滑角。

進(jìn)行坐標(biāo)轉(zhuǎn)換后，需對(duì)表面的單元進(jìn)行分類(lèi)，分為迎風(fēng)單元與背風(fēng)單元。根據(jù)單元屬性的不同，在后續(xù)的氣動(dòng)力計(jì)算中使用不同的計(jì)算方式。本文設(shè)計(jì)了可用于非凸表面的迎風(fēng)單元篩選算法。通過(guò)驗(yàn)證發(fā)現(xiàn)，隨著網(wǎng)格密度達(dá)到一定要求時(shí)，其具有較高的計(jì)算精度，滿足對(duì)復(fù)雜外形的超低軌衛(wèi)星氣動(dòng)力估算的需求。迎風(fēng)單元和背風(fēng)單元篩選算法的輸入包含如下數(shù)據(jù)：

1）節(jié)點(diǎn)位置Pi=[xi，yi，zi]，為編號(hào)為i節(jié)點(diǎn)其對(duì)應(yīng)在風(fēng)軸系下的空間坐標(biāo)位置；

2）單元列表Un=[u，v，w]，n為單元編號(hào)，u，v，w為三角形單元的3 個(gè)頂點(diǎn)節(jié)點(diǎn)編號(hào)。

算法的輸出包含以下幾個(gè)部分：

1）迎風(fēng)單元列表Wn=[u'，v'，w']，u'，v'，w'為三角形單元的3 個(gè)頂點(diǎn)節(jié)點(diǎn)編號(hào)；

2）背風(fēng)單元列表Ln=[u'，v'，w']，u'，v'，w'為三角形單元的3 個(gè)頂點(diǎn)節(jié)點(diǎn)的編號(hào)；

3）迎風(fēng)節(jié)點(diǎn)列表PWi=[xi，yi，zi]，xi，yi，zi為節(jié)點(diǎn)在風(fēng)軸系下的坐標(biāo)位置。

具體算法見(jiàn)表1。

表1 迎風(fēng)單元篩選算法流程Tab.1 Flow chart of the screening algorithm for windward units

從迎風(fēng)方向依次選取節(jié)點(diǎn)，判斷其是否在迎風(fēng)單元列表的投影內(nèi)，若不在任意投影內(nèi)，則將該節(jié)點(diǎn)加入迎風(fēng)節(jié)點(diǎn)列表中，并將與其連接的單元加入迎風(fēng)單元列表。最后將所有總單元列表Un減去篩選出的迎風(fēng)單元列表，則可得到背風(fēng)單元列表。具體篩選過(guò)程如圖2 所示，在對(duì)表面進(jìn)行離散后的衛(wèi)星網(wǎng)格中，黑色單元為篩選出的迎風(fēng)單元，灰色的是待篩選的單元列表。紅色節(jié)點(diǎn)與藍(lán)色節(jié)點(diǎn)為兩個(gè)待篩選節(jié)點(diǎn)。在圖示的沿著來(lái)流方向的投影平面上，紅色節(jié)點(diǎn)在黑色單元投影外，所以其為新的迎風(fēng)單元；而藍(lán)色節(jié)點(diǎn)在當(dāng)黑色單元的投影內(nèi)，其為背風(fēng)節(jié)點(diǎn)。

為驗(yàn)證該算法的準(zhǔn)確性，以及提出在后續(xù)計(jì)算中劃分網(wǎng)格的精細(xì)程度，本文設(shè)計(jì)了如下算例，如圖3（a）所示，該圖所示的梨狀外形，為一個(gè)非凸的幾何結(jié)構(gòu)，其在軸向投影方向的最大投影半徑R=0.5 以此作為模型的特征尺度，因此可以得到其準(zhǔn)確的投影面積Sproj=πR2=0.25π。本文依次增加特征尺度與三角形網(wǎng)格大小R/δx之間的比值，觀察通過(guò)迎風(fēng)單元與背風(fēng)單元算法計(jì)算出的投影面積與實(shí)際投影面積的差距，驗(yàn)證算法精度及求解合適的特征尺度與網(wǎng)格大小的比值關(guān)系R/δx。計(jì)算的結(jié)果如圖3（b）所示，發(fā)現(xiàn)在比值為10 時(shí)，算法計(jì)算出的投影面積與實(shí)際的投影面積十分接近，進(jìn)一步縮小網(wǎng)格尺度，增加特征尺度與網(wǎng)格大小的比值，精度不再改善。因此，本文認(rèn)為該算法的準(zhǔn)確性良好，在對(duì)后續(xù)的計(jì)算模型進(jìn)行表面單元的離散過(guò)程中，將選取特征尺度與網(wǎng)格大小的比值R/δx=10。

圖3 迎風(fēng)單元與背風(fēng)單元篩選算例模型Fig.3 Example model for screening windward and leeward units

1.2 氣動(dòng)力計(jì)算

在氣動(dòng)力的計(jì)算過(guò)程中，需要從大氣模型中根據(jù)低軌衛(wèi)星的飛行高度和經(jīng)緯度獲得當(dāng)前的流動(dòng)條件，獲得準(zhǔn)確的流動(dòng)條件對(duì)算法的精度與可靠性也有較大的影響。因此，選取合適的大氣模型也至關(guān)重要。當(dāng)前常用的大氣模型有以下幾種，US Standard 1976、MSISe-90、MSISe-00 和Jacchia-71等［11-13］。本文選取NRLMSISE-00 作為計(jì)算使用的大氣模型，主要基于其擁有以下良好性質(zhì)［14］：①使用廣泛，有豐富的程序接口可調(diào)用；②能夠?qū)崿F(xiàn)努森數(shù)（Knudsen Number，Kn）與海拔高度的轉(zhuǎn)換；③其數(shù)據(jù)豐富，滿足本文的計(jì)算需求。

本文選取的模型計(jì)算條件見(jiàn)表2。

表2 大氣模型參數(shù)Tab.2 Parameters of the atmospheric model

為方便后續(xù)計(jì)算，本文固定部分大氣模型參數(shù)（見(jiàn)表2），只保留高度? 作為變量。以上數(shù)據(jù)取太陽(yáng)豐度為平均值時(shí)的大氣系數(shù)，能反映當(dāng)前的大氣較為普遍的狀態(tài)。選取以上參數(shù)后，本文假設(shè)大氣模型為函數(shù)Fa，其輸入變量為高度?，則其可輸出Kn、T、m、ρ分別為當(dāng)前高度下的努森數(shù)、大氣溫度、大氣的摩爾質(zhì)量與大氣密度，即：

根據(jù)不同高度對(duì)應(yīng)的流區(qū)的不同，本文將計(jì)算分為3 個(gè)部分，見(jiàn)表3。

表3 對(duì)應(yīng)不同流區(qū)適用的不同氣動(dòng)力計(jì)算方式Tab.3 Aerodynamic calculation methods for different flow regions

1.2.1 自由分子流模型

在自由分子流中，分子之間的碰撞概率很小，經(jīng)碰撞后反射分子對(duì)流動(dòng)的影響可以忽略不計(jì)［6］。因此，主要關(guān)注的過(guò)程是氣體分子與衛(wèi)星表面發(fā)生碰撞（見(jiàn)圖4），分子與衛(wèi)星表面發(fā)生能量與動(dòng)量的交換情況。因此，自由分子流動(dòng)時(shí)衛(wèi)星的受力狀況是由氣體-表面相互作用性質(zhì)決定的。在假設(shè)氣體僅與衛(wèi)星表面發(fā)生動(dòng)量交換的前提下，作用在衛(wèi)星表面上的力等于氣體動(dòng)量的變化率。對(duì)于一個(gè)外凸形狀的衛(wèi)星，其單位面積受力大小與氣體的入射動(dòng)量與反射動(dòng)量的大小有關(guān)。以一個(gè)衛(wèi)星的表面單元為例，其受到分子的作用力可表示為

圖4 氣體分子與單元表面碰撞模型Fig.4 Model for the collision of gas molecule and unit surface（a）Maxwell 模型（b）牛頓理論

式中：f為分子的作用力；A為衛(wèi)星表面的單位面積；p為法向的動(dòng)量通量；τ為切向的動(dòng)量通量。下標(biāo)a、b 分別對(duì)應(yīng)著入射通量和反射通量，法向和切向入射動(dòng)量流pa、τa取決于入射速度（V）和質(zhì)量流（dQ）；其反射動(dòng)量通量為-pbn+τbt；n為在物面表面的法向分量；t為在物面表面的切向分量。

在自由分子流中，可以認(rèn)為

式中：dQ=ρVcosδ，δ為粒子入射方向與平面之間的夾角，ρ為分子流密度；V為分子流速度大小。

在流動(dòng)中氣體分子的速度服從Maxwell 速度分布函數(shù)［15］。則入射動(dòng)量通量可以表示為

式中：m為單個(gè)氣體分子的質(zhì)量；F(u)為麥克斯韋分布函數(shù)。

式中：s為分子速度比；V為在速度分布函數(shù)下最可能的速度大小；Rc為氣體常數(shù)；T∞為當(dāng)前高度下的來(lái)流溫度。

目前，如何求解確定反射通量是個(gè)難點(diǎn)。有許多采用簡(jiǎn)化的氣體-表面作用模型來(lái)對(duì)反射動(dòng)量通量進(jìn)行建模，如Maxwell［15］、Schamberg、Schaaf&Chambre［6］等模型。本 文選取Sentman 模型［16］進(jìn)行建模，假設(shè)一個(gè)粒子在撞擊后會(huì)依據(jù)一定的概率分布向不同方向進(jìn)行漫反射，定義分子與平面發(fā)生碰撞之后，平面的法向力系數(shù)與切向力系數(shù)分別為

式中：Tk，b、Tk，a為

式中：αacc為能量調(diào)節(jié)系數(shù)，根據(jù)飛行器的高度確定［16］，隨著飛行高度的增加，能量調(diào)節(jié)系數(shù)的取值不斷減?。籘w為飛行器的表面壁溫。

1.2.2 連續(xù)流區(qū)模型

牛頓理論（牛頓正弦平方律）闡述了一種推導(dǎo)平面與流體相互作用的理論模型，其假設(shè)氣體分子為質(zhì)點(diǎn)且互相孤立，氣體與平面撞擊后，其動(dòng)量在物面的法向上發(fā)生完全的動(dòng)量交換，分子沿著物面的切向動(dòng)量完全保留。其理論在低速流動(dòng)條件下偏差較大，但物體在高速連續(xù)流區(qū)運(yùn)動(dòng)中，其對(duì)氣動(dòng)力的估算效果良好。物面受到的法向力系數(shù)CP為

牛頓理論提出的模型僅僅與氣流方向與物面的夾角有關(guān)，在馬赫數(shù)不夠大時(shí)，會(huì)出現(xiàn)明顯的偏差。為了對(duì)該偏差進(jìn)行修正，Lester Lee 提出了新的計(jì)算方式：

其中：

式中：Ma∞為自由來(lái)流的馬赫數(shù)。

在牛頓理論中，背風(fēng)單元不受氣動(dòng)力的作用，但在連續(xù)流區(qū)中，忽略背風(fēng)單元的氣動(dòng)力會(huì)導(dǎo)致計(jì)算出現(xiàn)明顯的誤差。因此，本文使用以下普朗特-邁耶膨脹波理論估算背風(fēng)單元的氣動(dòng)系數(shù)［17］

式中：γ為絕熱氣體指數(shù)，在量熱完全氣體的假設(shè)下γ為1.4。

在實(shí)際估算中，以述方法能夠較好地估計(jì)連續(xù)流區(qū)中的氣動(dòng)力系數(shù)值，但由于缺少黏性力的估算，依然會(huì)有一定的偏差，為做進(jìn)一步修正，本文使用可壓縮流的摩擦系數(shù)計(jì)算公式對(duì)摩擦力系數(shù)進(jìn)行估算。

1.2.3 過(guò)渡流區(qū)模型

在過(guò)渡流區(qū)的氣動(dòng)力計(jì)算中，本文使用橋函數(shù)的方式進(jìn)行計(jì)算，即在計(jì)算出連續(xù)流區(qū)與自由分子流區(qū)的氣動(dòng)特性參數(shù)后，根據(jù)在過(guò)渡流區(qū)中的高度h由大氣模型推導(dǎo)出當(dāng)前高度下的努森數(shù)Kn對(duì)連續(xù)流區(qū)與自由分子流區(qū)的氣動(dòng)特性參數(shù)進(jìn)行權(quán)重分配，得到當(dāng)前位置的氣動(dòng)特性參數(shù)。

式中：Cfm為在自由分子流中求得的氣動(dòng)力參數(shù)；Ccont為在連續(xù)流區(qū)中求得的氣動(dòng)力參數(shù)；Pb為橋函數(shù)。當(dāng)前常用的橋函數(shù)有sine-squared 橋函數(shù)與log-erf 橋函數(shù)。

在設(shè)定合理的參數(shù)條件下，2 種函數(shù)均表現(xiàn)出良好的擬合性能。以上2 種函數(shù)都需要根據(jù)求解物體的幾何外形的變化與來(lái)流條件的改變對(duì)相關(guān)參數(shù)進(jìn)行修正。在實(shí)際使用過(guò)程中，需要使用部分實(shí)驗(yàn)數(shù)據(jù)或高精度的求解結(jié)果對(duì)參數(shù)進(jìn)行調(diào)整，這增加了工作量。

為了對(duì)上述問(wèn)題進(jìn)行改善，本文提出了以logistic 函數(shù)為基礎(chǔ)的橋函數(shù)，并在函數(shù)中加入了形狀修正因子，拓展了橋函數(shù)的適用范圍。本文設(shè)計(jì)的橋函數(shù)如下：

式中：q為形狀因子，其與氣體來(lái)流方向與物面間的夾角有關(guān)；橋函數(shù)為與當(dāng)?shù)豄n數(shù)和形狀有關(guān)的函數(shù)。

2 模型驗(yàn)證

為驗(yàn)證提出的快速算法的有效性及正確性，本文結(jié)合DSMC 的計(jì)算結(jié)果，對(duì)以上方法進(jìn)行驗(yàn)證。吳子牛等［10］曾對(duì)圓柱體與鈍頭雙楔體進(jìn)行DSMC與橋函數(shù)的研究，其研究對(duì)象如圖5 所示，2 個(gè)模型展向上的長(zhǎng)度均為1 000 mm。

圖5 幾何外形Fig.5 Geometric shape

計(jì)算條件見(jiàn)表4，給定模型壁溫分別500 K、300 K；自由來(lái)流溫度條件給定為300 K。選取圓柱模型的速度工況為馬赫數(shù)4 與馬赫數(shù)16，鈍頭雙楔體模型速度工況為馬赫數(shù)4 與馬赫數(shù)8；選取計(jì)算努森數(shù)Kn為10-3～102，在這個(gè)范圍內(nèi)大氣特性變化顯著。其余輸入條件由Nrlmsise-00 大氣模型給出。

表4 模型驗(yàn)證的計(jì)算條件Tab.4 Calculation conditions for model verification

為了進(jìn)行對(duì)比，本文對(duì)DSMC 結(jié)果及使用本文提出的logistic-log 函數(shù)與sine-squared 函數(shù)進(jìn)行比較。圖6 表示了在本節(jié)所研究?jī)煞N不同模型中，本文所提出的logistic-log 橋函數(shù)與sine-squared 橋函數(shù)在2 種不同模型中不同條件下的其努森數(shù)（Kn）與阻力系數(shù)（CD）之間的關(guān)系如圖6 所示。在圖6（a）中，紅色線代表的是Ma=16 結(jié)果，藍(lán)色線代表的是Ma=4 結(jié)果。其中，粗虛線為logistic-log 橋函數(shù)的計(jì)算結(jié)果，而細(xì)虛線為sine-squared 橋函數(shù)的計(jì)算結(jié)果，三角符號(hào)代表著使用自由分子流算法（Free Molecular Fuction，F(xiàn)MF）的計(jì)算結(jié)果，其余圖像中的散點(diǎn)代表的是DSMC 計(jì)算數(shù)據(jù)。由圖6 可知，對(duì)于圓柱的計(jì)算工況，在Ma=16 的高速工況下，logisticlog 橋函數(shù)與sine-squared 橋函數(shù)與DSMC 結(jié)果吻合都較好。在Ma=4 的工況下，sine-squared 橋函數(shù)的計(jì)算結(jié)果在Kn數(shù)較高時(shí)與實(shí)驗(yàn)數(shù)據(jù)產(chǎn)生了較為明顯的差異。在圖6（b）中，紅色線代表的是Ma=8結(jié)果，藍(lán)色線代表的是Ma=4 結(jié)果。其中，粗虛線為logistic-log 橋函數(shù)的計(jì)算結(jié)果，而細(xì)虛線為sinesquared 橋函數(shù)的計(jì)算結(jié)果，三角符號(hào)代表著使用自由分子流算法（Free Molecular Fuction，F(xiàn)MF）的計(jì)算結(jié)果，其余圖像中的散點(diǎn)代表的是DSMC 計(jì)算數(shù)據(jù)。同樣可以發(fā)現(xiàn)，對(duì)于鈍頭雙楔體的計(jì)算工況，在速度較高的情況下（Ma=8），2 種橋函數(shù)與實(shí)驗(yàn)結(jié)果基本吻合。當(dāng)Ma=4 時(shí)，使用sine-squared 橋函數(shù)在Kn=12 時(shí)出現(xiàn)了最大值，與DSMC 數(shù)據(jù)相差較大，也違背了隨著Kn數(shù)的增加，氣動(dòng)阻力系數(shù)也隨之增加的基本規(guī)律。以上結(jié)果說(shuō)明了在固定系數(shù)的情況下，sine-squared 橋函數(shù)無(wú)法滿足在不同外形，多種工況下均保持良好的估算精度；而logistic-log 橋函數(shù)可以實(shí)現(xiàn)在不改變系數(shù)的情況下滿足對(duì)多種工況的氣動(dòng)系數(shù)估算。

圖6 不同橋函數(shù)對(duì)應(yīng)的氣動(dòng)特性估算結(jié)果Fig.6 Results of the aerodynamic characteristics corresponding to different bridge functions

為驗(yàn)證本文提出的形狀因子對(duì)其結(jié)果的影響，本 文對(duì)q=1+0.632 6cosδ+6.124 9與q=1 進(jìn)行對(duì)比，計(jì)算結(jié)果如圖7 所示。

在圖7（a）中，紅色線圖代表的是Ma=16 結(jié)果，藍(lán)色線代表的是Ma=4 結(jié)果。粗虛線為形狀因子q=1+0.632 6cosδ+6.124 9 的計(jì)算結(jié)果，而細(xì)虛線為形狀因子q=1 的計(jì)算結(jié)果，圖像中的散點(diǎn)代表的是DSMC 計(jì)算數(shù)據(jù)。由圖7 可知，對(duì)于圓柱的計(jì)算工況，在Ma=16 的情況下，2 種形狀因子與實(shí)驗(yàn)結(jié)果的吻合情況較好。在Ma=4 的計(jì)算條件下，不使用q=1+0.632 6cosδ+6.124 9 的形狀因子則會(huì)出現(xiàn)明顯的誤差。對(duì)于鈍頭雙楔體的計(jì)算工況，發(fā)現(xiàn)形狀因子在不同工況下的影響不大。因?yàn)樵阝g頭雙楔體中，表面外形的傾角相對(duì)固定，形狀因子對(duì)其影響則不明顯；而在圓柱工況中，圓柱各個(gè)位置的傾角均有變化，則形狀因子在速度較低的情況下起到了明顯的修正作用。在速度較低時(shí)，形狀因子起到的修正效果明顯，在高速情況下形狀因子的修正作用不顯著。為了進(jìn)一步量化計(jì)算結(jié)果與DSMC 數(shù)據(jù)的吻合程度，本文引入了均方誤差對(duì)計(jì)算結(jié)果進(jìn)行了比較，其計(jì)算公式如下：

圖7 形狀因子對(duì)計(jì)算結(jié)果的影響Fig.7 Effects of the shape factor on the calculation results

式中：ε為誤差值為DSMC 的采樣數(shù)據(jù)為使用不同橋函數(shù)得到的數(shù)值結(jié)果數(shù)據(jù)；k為DSMC 的數(shù)據(jù)維度。

通過(guò)對(duì)均方誤差進(jìn)行比較，計(jì)算出不同方法的誤差，見(jiàn)表5。可見(jiàn)本文提出的計(jì)算方式誤差最小。

表5 不同橋函數(shù)計(jì)算結(jié)果均方誤差（ε）Tab.5 Mean square errors of the calculation results of different bridge functions

綜上所述，本文提出logistic-log 橋函數(shù)與當(dāng)前常用的橋函數(shù)相比，具有適用范圍廣，估算精度高的特點(diǎn)。本文提出的形狀因子q能夠根據(jù)計(jì)算外形對(duì)計(jì)算結(jié)果進(jìn)行修正。

3 典型超低軌飛行器氣動(dòng)力計(jì)算

歐洲空間局（European Space Agency，ESA）的地球重力與海洋環(huán)流探測(cè)衛(wèi)星GOCE 為典型的超低軌飛行器［18］，有關(guān)學(xué)者也使用了DSMC 或TPMC方法對(duì)其氣動(dòng)特性進(jìn)行了研究［19-20］，本文第2 章提出的方法對(duì)其在工作高度范圍內(nèi)的阻力系數(shù)進(jìn)行測(cè)算。

在不同側(cè)滑角下，不同攻角下隨著工作高度的變化GOCE 的阻力系數(shù)與阻力大小的變化趨勢(shì)如圖8 所示。隨著工作高度的增加，阻力系數(shù)的變化范圍越小，GOCE 在俯仰與偏航時(shí)所受的氣動(dòng)力變化相對(duì)較低。但是隨著運(yùn)行高度降低，大氣密度顯著增加，阻力呈幾何量級(jí)增大。在平飛狀態(tài)下，阻力的最大值為42.526 N，最小值為0.005 41 N。因此，超低軌衛(wèi)星在進(jìn)行工作時(shí)，需要及時(shí)進(jìn)行姿軌控制維持其高度。

圖8 不同高度下GOCE 的阻力與阻力系數(shù)大小與其姿態(tài)角度關(guān)系Fig.8 Relationships between the drag or drag coefficient of GOCE and its attitude angle at different altitudes

續(xù)圖8 不同高度下GOCE 的阻力與阻力系數(shù)大小與其姿態(tài)角度關(guān)系ContinuedFig.8 Relationships between the drag or drag coefficient of GOCE and its attitude angle at different altitudes

GOCE 形態(tài)細(xì)長(zhǎng)，姿態(tài)變化導(dǎo)致的氣動(dòng)力變化明顯，在運(yùn)行高度在100 km 時(shí)，其阻力系數(shù)變化范圍在3.2～9.8；在高度為260 km 時(shí)其阻力系數(shù)的變化范圍在4.1～9.5。其太陽(yáng)能帆板占主要浸濕面積，在偏航時(shí)，迎風(fēng)面積變化明顯，導(dǎo)致了其在側(cè)滑時(shí)的阻力系數(shù)與阻力變化明顯高于俯仰時(shí)的情況。

4 結(jié)束語(yǔ)

1）本文提出了一種超低軌飛行器的快速氣動(dòng)力計(jì)算工具方法。使用面元法進(jìn)行求解，通過(guò)將飛行器進(jìn)行便面單元離散為迎風(fēng)單元與背風(fēng)單元，分別求解不同種類(lèi)單元的氣動(dòng)力系數(shù)，通過(guò)對(duì)離散單元的數(shù)值結(jié)果進(jìn)行積分得到了飛行器的氣動(dòng)力特性。通過(guò)重新設(shè)計(jì)橋函數(shù)，降低計(jì)算方法的數(shù)據(jù)依賴程度，提高了應(yīng)用范圍。

2）使用2 種橋函數(shù)計(jì)算，并對(duì)比鈍頭雙楔體與圓柱體的DSMC 結(jié)果，說(shuō)明了算法的可靠性和準(zhǔn)確性，解釋了本文提出的在計(jì)算不同工況氣動(dòng)力時(shí)，形狀因子對(duì)計(jì)算結(jié)果有修正作用。

3）分析闡述了典型超低軌飛行器GOCE 在100～260 km 工作范圍內(nèi)隨著俯仰角與偏航角變化氣動(dòng)特性的變化趨勢(shì)。發(fā)現(xiàn)在發(fā)生姿態(tài)變化時(shí)的阻力系數(shù)變化顯著，其變化范圍在3.2～9.8，隨著高度的降低阻力系數(shù)變化的范圍越大。因此，飛行器在低軌工作時(shí)，需要為飛行器設(shè)計(jì)舵面或陀螺儀對(duì)其姿態(tài)進(jìn)行控制。在過(guò)渡流區(qū)內(nèi)，隨著高度增加，氣動(dòng)阻力減小，阻力最大值為42.526 N，最小值為0.005 41 N。由此可見(jiàn)，超低軌飛行器需要在寬域過(guò)渡流區(qū)內(nèi)運(yùn)行，需要額外的動(dòng)力裝置。