李錦，江定武，黎昊旻，萬釗，王新光，毛枚良，陳堅強

（中國空氣動力研究與發(fā)展中心計算空氣動力研究所，綿陽 621000）

1 引言

立方星（CubeSat）是指尺寸為10cm×10cm×10cm，質(zhì)量一般在1kg 左右的微納衛(wèi)星，簡稱1U，其概念由加州理工大學(xué)和斯坦福大學(xué)于1999 年首先提出[1]。僅就單顆衛(wèi)星而言，過小的體積意味著其缺少推進系統(tǒng)，功能有限?？紤]到成本很低，研制周期短，發(fā)射靈活，因而主要被大學(xué)用于教學(xué)目的。但這種衛(wèi)星高度模塊化，可以根據(jù)需要進行組裝，形成2U～6U，甚至幾十U 的星座組合形式，功能大為拓展，且可避免由于單顆衛(wèi)星失效而導(dǎo)致的任務(wù)失敗，能夠執(zhí)行更為復(fù)雜的科學(xué)探測、全球通信等任務(wù)，具有極強的應(yīng)用價值。這是一項很有前景的顛覆性創(chuàng)新技術(shù)，近年來成為全球各科研院校和機構(gòu)的研究熱點之一。

然而，立方星的大量發(fā)射加劇了近地軌道中的空間碎片問題?？臻g碎片將給在軌運行的航天器帶來風險。阻力帆提供了解決這一問題的有效方案——利用氣動阻力使得衛(wèi)星被動離軌并進入大氣層燒毀[2]。與使用推進器主動離軌的技術(shù)不同，阻力帆技術(shù)不需要衛(wèi)星維持功能部件及儲存推進劑[1]。

離軌周期的計算需要對稀薄氣動特性做出精確預(yù)測。王立武等[3]采用DSMC開源軟件SPARTA數(shù)值模擬了400～800km 高度六邊形阻力帆離軌裝置的氣動特性，發(fā)現(xiàn)阻力系數(shù)隨著攻角的增大而減小，對軌道高度不敏感。Alexeenko 等[4]利用SPARTA 軟件初步分析了Cubesat 衛(wèi)星在185km、125km 和100km 高度的氣動特性。文中將轉(zhuǎn)動和振動松弛碰撞數(shù)指定為常數(shù)，并指出為了增加結(jié)果的保真度，應(yīng)該將它們設(shè)置為溫度的函數(shù)。在Alexeenko 等工作的基礎(chǔ)上，本文采用自研的NNW-DSMC軟件[5]，根據(jù)單元溫度計算轉(zhuǎn)動和振動松弛碰撞數(shù)，并評估不同處理方式對結(jié)果的影響。本文計算也是對NNW-DSMC軟件計算精度的進一步驗證[6]。

2 問題描述

美國普渡大學(xué)開展了一項立方星氣動離軌試驗（Aerodynamic Deorbit Experiment CubeSat）。立方星的形式為1U，由宇宙神V火箭送至地球軌道并釋放。初始軌道為地球同步轉(zhuǎn)移軌道，軌道近地點高度為185km，近地點速度約為10.239km/s。一旦立方星布置到這一軌道，就立即展開阻力帆，啟動離軌試驗。

圖1 給出了計算外形示意圖。立方星后接一個金字塔型的阻力帆，阻力帆底面為正方形，邊長為1160mm。計算條件見表1。來流攻角只考慮0°，不考慮側(cè)滑角。

表1 計算條件Table 1 Calculation conditions

圖1 計算外形示意圖Fig. 1 Calculation of outline diagram

3 DSMC方法及軟件

3.1 DSMC方法

DSMC方法由G A Bird提出[7]，是模擬稀薄氣體流動應(yīng)用最廣泛的方法之一。在過去的半個多世紀，該方法被廣泛用于不同的領(lǐng)域并經(jīng)受了實驗結(jié)果的驗證。隨著稀薄程度的增加，NS方程解算器的可信度降低。適合描述所有稀薄程度條件下流體行為的方程是Boltzmann方程。

式中：nf代表數(shù)密度n與速度分布函數(shù)f的乘積，c是分子速度，cr是相對分子速度，F(xiàn)是外力，σ是碰撞截面積，Ω是立體角。角標*代表碰撞后的值，f和1f分別是兩個不同分子的分布函數(shù)。

在DSMC 方法中，模擬分子代表大量的氣體分子，可以相互碰撞并且運動。計算域被劃分為網(wǎng)格并且在其中選擇碰撞對。最常用的碰撞對選擇格式是非時間計數(shù)算法（No Time Counter,NTC）。在NTC方法中，待測試的最大碰撞對數(shù)目Nc可如下計算：

式中：N是單元內(nèi)的模擬分子數(shù)，NF是一個模擬分子所代表的真實分子數(shù)，代表碰撞截面積與碰撞相對速度乘積的最大值，tΔ 是時間步長，cV是單元體積。

在選擇好碰撞對之后，就可以根據(jù)不同的碰撞模型進行碰撞計算。碰撞模型分為彈性碰撞和非彈性碰撞。彈性碰撞模型有硬球（Hard Sphere,HS）模型、變徑硬球（Variable Hard Sphere, VHS）模型、變徑軟球（Variable Soft Sphere, VSS）模型、廣義硬球模型（Generalized Hard Sphere, GHS）和廣義軟件模型（Generalized Soft Sphere, GSS）等。其中VHS 模型是應(yīng)用最廣泛的?；赩HS 模型，氣體的平均自由程λ可如下計算：

式中：ω是黏性系數(shù)-溫度指數(shù)，m是氣體質(zhì)量，k是Boltzmann 常數(shù)，T是溫度，μ是黏性系數(shù)，ρ是密度。

非彈性碰撞通常采用Larsen-Borgnakke 模型來完成內(nèi)能的松弛計算。

模擬分子不僅相互之間碰撞，還與固體壁面發(fā)生碰撞。壁面反射模型采用完全漫反射模型。反射模型分子的速度根據(jù)Maxwell分布給定。

最終氣體的宏觀性質(zhì)，比如壓力、溫度、密度和速度等，通過對模擬分子的微觀信息抽樣而得到。

DSMC 方法的詳細執(zhí)行可參見Bird 的專著[7]。具體執(zhí)行中，要得到可信結(jié)果，需要對一些參數(shù)如網(wǎng)格大小、時間步長、每單元模擬分子數(shù)、抽樣數(shù)等進行合理的選取。通常網(wǎng)格大小要小于分子平均自由程的1/3，時間步長應(yīng)小于平均碰撞時間，模擬分子在一個時間步長內(nèi)的運動距離不超過網(wǎng)格大小。每單元模擬分子數(shù)應(yīng)至少在10以上，抽樣步數(shù)應(yīng)使得統(tǒng)計誤差收斂至可接受水平。

3.2 NNW-DSMC軟件

國家數(shù)值風洞（National Numerical Wind tunnel，NNW）工程是國家“十三五”重大示范項目之一。工程建設(shè)目標是自主研發(fā)功能先進、種類齊全的計算流體力學(xué)軟件系統(tǒng)?；贒SMC的稀薄流數(shù)值模擬子系統(tǒng)（NNW-DSMC）是CFD軟件及集成框架系統(tǒng)下屬專用軟件子系統(tǒng)的建設(shè)內(nèi)容之一。其建設(shè)目標是研發(fā)面向工程實際、適用于大規(guī)模并行計算的DSMC軟件研究及應(yīng)用平臺，達到與國際主流DSMC計算軟件相當水平，以滿足航天航空等領(lǐng)域的需要。NNW-DSMC軟件由中國空氣動力研究與發(fā)展中心計算所與中國科學(xué)院力學(xué)研究所[8]聯(lián)合開發(fā)，軟件功能需求和概要設(shè)計由氣動中心提出，詳細設(shè)計及代碼開發(fā)在氣動中心的組織下主要由力學(xué)所承擔，軟件測試由雙方共同完成。

一個好的計算平臺應(yīng)具備以下特征[9]：軟件系統(tǒng)架構(gòu)科學(xué)合理，易于拓展及維護；功能豐富、計算準確、覆蓋核心需求；高效計算，快速提供可靠解決方案；界面友好，方便用戶學(xué)習使用。瞄準上述目標，我們通過合理利用C++繼承、派生等語言特性，最大程度增加代碼復(fù)用性。在構(gòu)建軟件框架時，對主要過程（如粒子碰撞、粒子移動）進行合理的接口設(shè)計，利用C++多態(tài)機制，實現(xiàn)“面向接口”的軟件開發(fā)過程，保證軟件具有良好的可閱讀性、可維護性和二次開發(fā)能力。幾何模型支持二維、軸對稱以及三維情形。軟件功能不僅覆蓋DSMC 常規(guī)算法，還包括若干與DSMC 相關(guān)的衍生算法，如可降低低速流動統(tǒng)計噪聲的信息保存法（IP-DSMC）、用于跨流域與多尺度流動的大網(wǎng)格DSMC（LC-DSMC）以及NS-DSMC耦合算法等。軟件集成豐富的氣體組分以及化學(xué)反應(yīng)數(shù)據(jù)庫，可應(yīng)用于地球/火星等再入流動計算場景。大規(guī)模并行計算依托第三方開源庫METIS 完成。此外，軟件還提供完善的輸入和輸出接口，方便用戶快速學(xué)習、使用。軟件主要性能已實現(xiàn)與國際同類型主流軟件相當，對標軟件主要為開源軟件SPARTA[10,11]。

4 結(jié)果與分析

本文計算中，碰撞抽樣采用NTC 算法，彈性碰撞采用VHS 模型，非彈性碰撞采用Larsen-Borgnakke模型，轉(zhuǎn)動松弛碰撞數(shù)和振動松弛碰撞數(shù)為常數(shù)（分別為5 和50），具體參數(shù)見表2。化學(xué)反應(yīng)采用5 組分19 反應(yīng)模型[12]。壁面采用等溫完全漫反射模型，不考慮壁面催化反應(yīng)。并行分區(qū)由METIS開源軟件完成。

表2 VHS模型參數(shù)Table 2 VHS model parameters

為了得到保真度更高的數(shù)據(jù)，轉(zhuǎn)動和振動松弛碰撞數(shù)應(yīng)該設(shè)置為溫度的函數(shù)。

轉(zhuǎn)動松弛碰撞數(shù)的計算公式為

振動松弛碰撞數(shù)的計算公式為

式中：Ttr是平動溫度，T是總溫度。C2為固定參數(shù)，見表3。

表3 轉(zhuǎn)動和振動松弛碰撞數(shù)參數(shù)Table 3 Rotation and vibration relaxation collision parameters

圖2 給出氮氣和氧氣的轉(zhuǎn)動和振動松弛碰撞數(shù)隨溫度的變化情況，可以看出，兩種氣體的轉(zhuǎn)動松弛碰撞數(shù)都隨溫度的升高而增大，且氮氣的結(jié)果要高于氧氣的結(jié)果。當溫度超過1000k以后，轉(zhuǎn)動松弛碰撞數(shù)就超過常規(guī)設(shè)置的常數(shù)值（5）。因此按照常數(shù)設(shè)置，轉(zhuǎn)動松弛碰撞數(shù)將低于實際情況。兩種氣體的振動松弛碰撞數(shù)則隨溫度的升高而減小，氮氣的結(jié)果同樣要高于氧氣的結(jié)果。只有當溫度超過約30000k以后，振動松弛碰撞數(shù)才能夠與常數(shù)值（50）相當。因此按照常數(shù)設(shè)置，振動松弛碰撞數(shù)也將低于實際情況。

圖2 轉(zhuǎn)動和振動松弛碰撞數(shù)隨溫度的變化Fig. 2 The variation of rotational and vibrational relaxation collision numbers with temperature

圖3 給出不同高度下對稱面流場密度和溫度分布，可以看出，隨著高度的增加，密度越來越低，激波的影響范圍越來越大，呈現(xiàn)出明顯的激波內(nèi)部結(jié)構(gòu)，與連續(xù)流時的間斷分布完全不同。

圖3 不同高度下對稱面流場密度和溫度分布Fig. 3 Density and temperature distribution of symmetric plane flow field at different heights

圖4 給出100km 高度時駐點線上的溫度分布，可以看出存在明顯的熱力學(xué)非平衡效應(yīng)。在這個高度，分子間的碰撞已經(jīng)不是非常劇烈，以致于分子無法充分地將它們的能量向內(nèi)模態(tài)進行松弛。平動溫度要遠高于轉(zhuǎn)動和振動溫度。

圖4 駐點線上溫度分布（H=100km）Fig. 4 Temperature distribution on the stagnation line（H=100km）

圖5 給出化學(xué)反應(yīng)以及松弛碰撞數(shù)對駐點線溫度分布的影響?？紤]化學(xué)反應(yīng)后，由于吸熱的離解反應(yīng)占據(jù)主導(dǎo)，流場中溫度會下降。而轉(zhuǎn)動和振動松弛碰撞數(shù)考慮為溫度的函數(shù)后，轉(zhuǎn)動和振動溫度降低，平動溫度則升高，流場總溫度也增加了約10%。這是因為隨著溫度的升高，轉(zhuǎn)動和振動松弛碰撞數(shù)增大，而轉(zhuǎn)動和振動松弛的概率降低。

圖5 駐點線上溫度分布（H=100km）Fig. 5 Temperature distribution on the stagnation line（H=100km）

表4 給出本文計算的阻力系數(shù)和最大熱流與參考文獻結(jié)果的對比（文獻中沒有給出185km 高度的計算結(jié)果）。從表中可以看出，NNW-DSMC計算的阻力系數(shù)與參考結(jié)果非常接近，相對誤差小于1%，最大熱流的相對誤差在6%以內(nèi)。這充分表明，NNW-DSMC 軟件具有良好的計算精度。盡管阻力系數(shù)隨著高度的降低而減小，但是阻力還是增加的。

表4 阻力系數(shù)和最大熱流的對比Table 4 Comparison of drag coefficient and maximum heat flux

圖6 給出立方星阻力帆的阻力隨高度的變化情況。本文計算結(jié)果與文獻及自由分子流理論計算結(jié)果吻合較好。隨著高度的降低，大氣密度逐漸升高，阻力隨之增大。特別地，當高度下降到120km 以下時，阻力急劇增大。100km 時的阻力約為46.5N，大約是120km 時（0.93N）的50 倍，185km時（0.05N）的930倍。需要指出的是，就工程實踐而言，立方星阻力帆在下降到約100km時，可能就已經(jīng)被大氣阻力撕裂而不再保持原始帆的形態(tài)。自研NNW-DSMC 軟件還不具備考慮阻力帆在下降過程中受“阻力”和“高溫”影響發(fā)生形變的能力，相關(guān)研究尚待未來進一步的工作。

圖6 阻力隨高度的變化（NNW-DSMC：本文計算結(jié)果，SPARTA：文獻[4]中結(jié)果，F(xiàn)M：自由分子流結(jié)果）Fig. 6 The variation of resistance with height

5 結(jié)論

本文采用自研NNW-DSMC 軟件，對立方星阻力帆離軌裝置在不同高度下的稀薄繞流流場進行數(shù)值模擬，分析了氣動特性的變化規(guī)律。結(jié)果表明：

（1）流場存在顯著的熱力學(xué)非平衡現(xiàn)象，隨著高度的增加，激波的影響范圍越來越大，呈現(xiàn)出明顯的激波內(nèi)部結(jié)構(gòu)，與連續(xù)流時的間斷分布完全不同。

（2）化學(xué)反應(yīng)對流場具有一定影響，100km時駐點線上的溫度峰值比不考慮化學(xué)反應(yīng)時要低3000k左右；對氣動特性的影響較小。

（3）將碰撞松弛數(shù)設(shè)置為溫度的函數(shù)后，轉(zhuǎn)動和振動溫度降低，平動溫度升高，流場總溫度也增加了約10%，模擬具有更高的保真度；對氣動特性的影響較小。

（4）阻力和最大熱流隨著高度下降而增大，阻力系數(shù)則隨著高度下降而減小；阻力系數(shù)與參考結(jié)果的相對誤差小于1%，最大熱流的相對誤差在6%以內(nèi)；NNW-DSMC 的計算結(jié)果與開源軟件SPARTA結(jié)果吻合良好。

（5）立方星阻力帆在下降到約100km 時，氣動阻力急劇增大，可能就已經(jīng)被大氣阻力撕裂而不再保持原始帆的形態(tài)。相關(guān)研究尚待未來進一步的工作。

本文結(jié)果可為后續(xù)的立方星離軌動力學(xué)分析提供支撐。為得到更高的保真度，應(yīng)將轉(zhuǎn)動和振動松弛碰撞數(shù)考慮為溫度的函數(shù)。