立方星阻力帆氣動特性的DSMC 模擬

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

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

1 引言

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

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

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

2 問題描述

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

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

表1 計算條件Table 1 Calculation conditions

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

3 DSMC方法及軟件

3.1 DSMC方法

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

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

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

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

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

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

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

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

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

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

3.2 NNW-DSMC軟件

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

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

4 結果與分析

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

表2 VHS模型參數Table 2 VHS model parameters

為了得到保真度更高的數據，轉動和振動松弛碰撞數應該設置為溫度的函數。

轉動松弛碰撞數的計算公式為

振動松弛碰撞數的計算公式為

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

表3 轉動和振動松弛碰撞數參數Table 3 Rotation and vibration relaxation collision parameters

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

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

圖3 給出不同高度下對稱面流場密度和溫度分布，可以看出，隨著高度的增加，密度越來越低，激波的影響范圍越來越大，呈現出明顯的激波內部結構，與連續流時的間斷分布完全不同。

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

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

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

圖5 給出化學反應以及松弛碰撞數對駐點線溫度分布的影響。考慮化學反應后，由于吸熱的離解反應占據主導，流場中溫度會下降。而轉動和振動松弛碰撞數考慮為溫度的函數后，轉動和振動溫度降低，平動溫度則升高，流場總溫度也增加了約10%。這是因為隨著溫度的升高，轉動和振動松弛碰撞數增大，而轉動和振動松弛的概率降低。

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

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

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

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

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

5 結論

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

（1）流場存在顯著的熱力學非平衡現象，隨著高度的增加，激波的影響范圍越來越大，呈現出明顯的激波內部結構，與連續流時的間斷分布完全不同。

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

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

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

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

本文結果可為后續的立方星離軌動力學分析提供支撐。為得到更高的保真度，應將轉動和振動松弛碰撞數考慮為溫度的函數。