立方星阻力帆氣動(dòng)特性的DSMC 模擬

李錦，江定武，黎昊旻，萬釗，王新光，毛枚良，陳堅(jiān)強(qiáng)

（中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所，綿陽 621000）

1 引言

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

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

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

2 問題描述

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

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

表1 計(jì)算條件Table 1 Calculation conditions

圖1 計(jì)算外形示意圖Fig. 1 Calculation of outline diagram

3 DSMC方法及軟件

3.1 DSMC方法

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

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

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

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

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

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

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

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

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

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

3.2 NNW-DSMC軟件

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

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

4 結(jié)果與分析

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

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

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

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

振動(dòng)松弛碰撞數(shù)的計(jì)算公式為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

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

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

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