郭晨林，陳方，趙艷彬，尤超藍，錢勇，李文龍，楊麗麗

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

0 引言

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

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

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

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

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

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

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

1.1 數(shù)據預處理

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

圖2 迎風單元與背風單元篩選Fig.2 Screening diagram of windward and leeward units

式中：Cbo為體坐標系ObXbYbZb至風軸坐標系OwXwYwZw下的旋轉矩陣；α、β分別為當前飛行器的迎角、側滑角。

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

1）節(jié)點位置Pi=[xi，yi，zi]，為編號為i節(jié)點其對應在風軸系下的空間坐標位置；

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

算法的輸出包含以下幾個部分：

1）迎風單元列表Wn=[u'，v'，w']，u'，v'，w'為三角形單元的3 個頂點節(jié)點編號；

2）背風單元列表Ln=[u'，v'，w']，u'，v'，w'為三角形單元的3 個頂點節(jié)點的編號；

3）迎風節(jié)點列表PWi=[xi，yi，zi]，xi，yi，zi為節(jié)點在風軸系下的坐標位置。

具體算法見表1。

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

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

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

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

1.2 氣動力計算

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

本文選取的模型計算條件見表2。

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

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

根據不同高度對應的流區(qū)的不同，本文將計算分為3 個部分，見表3。

表3 對應不同流區(qū)適用的不同氣動力計算方式Tab.3 Aerodynamic calculation methods for different flow regions

1.2.1 自由分子流模型

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

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

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

在自由分子流中，可以認為

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

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

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

式中：s為分子速度比；V為在速度分布函數(shù)下最可能的速度大?。籖c為氣體常數(shù)；T∞為當前高度下的來流溫度。

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

式中：Tk，b、Tk，a為

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

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

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

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

其中：

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

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

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

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

1.2.3 過渡流區(qū)模型

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

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

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

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

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

2 模型驗證

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

圖5 幾何外形Fig.5 Geometric shape

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

表4 模型驗證的計算條件Tab.4 Calculation conditions for model verification

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

圖6 不同橋函數(shù)對應的氣動特性估算結果Fig.6 Results of the aerodynamic characteristics corresponding to different bridge functions

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

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

圖7 形狀因子對計算結果的影響Fig.7 Effects of the shape factor on the calculation results

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

通過對均方誤差進行比較，計算出不同方法的誤差，見表5。可見本文提出的計算方式誤差最小。

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

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

3 典型超低軌飛行器氣動力計算

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

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

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

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

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

4 結束語

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

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

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