增阻離軌裝置氣動特性的DSMC研究

王立武,馮瑞,張九雙,王廣興,王奇,何青松

(1.東南大學國家預應力工程技術研究中心,南京 211189;2.北京空間機電研究所,北京 100094;3.中國航天科技集團有限公司航天進入、減速與著陸技術實驗室,北京 100094)

1 引言

隨著人類太空探索活動的日益頻繁,發射任務過程中產生的大量火箭拋棄物和廢棄航天器滯留在太空中[1]。其中,近地軌道 (LEO)區域尤為嚴重,容納了約90%的空間碎片,對當前和未來的航天器在軌運行造成了巨大的安全隱患[2]。控制空間碎片數量已經成為業界共識,并成為世界主要航天國家的優先事項之一。目前,最有效的方法是將壽命末期的航天器從軌道上及時移除[3,4],主要的離軌方式包括推力離軌、太陽帆離軌、電動力繩離軌以及增阻離軌等。其中,增阻離軌技術以大收攏比、質量輕、成本低、不需額外燃料等優勢,在LEO空間碎片減緩方面具有廣闊的應用前景。其原理是通過展開大面積的增阻面來有效提升壽命末期航天器的面質比,同時減小彈道系數,充分利用稀薄大氣阻力對航天器進行快速降軌并最終進入稠密大氣層燒毀[5-8]。因此,準確預測所選增阻離軌裝置稀薄流區的氣動特性對總體設計具有重要意義。

在稀薄流域,分子之間間距較大,經典連續介質假設不再成立,廣泛用于飛行器流場數值模擬的NS方程、Euler方程等不再適用。故需要從微觀分子模型角度出發,采用稀薄氣體動力學方法開展研究。在求解稀薄氣體動力學問題的諸多方法中,DSMC方法被認為是20世紀下半葉稀薄氣體動力學的代表性成果,是當前工程實踐中求解稀薄流域問題的主流方法。

DSMC方法由G.A.Bird于1963年提出,并首次應用于氣體分子的內能松弛問題[9]。而后用于二維、三維流動數值模擬,并且成功應用于多個型號飛行器的氣動特性計算。所得結果與實際飛行測量得到的數據相符度較高,如美國航天飛機升阻比的預估。后來眾多學者不斷從網格技術、分子碰撞模型、計算效率等方面對該方法予以發展完善,并多次用于工程實踐中三維復雜外形流場模擬,比如AFE飛船、三角翼、“哥倫比亞”號航天飛機事故分析、火星探測器減速裝置Ballute氣動加熱分析。國內沈青、樊菁等最早針對DSMC方法也開展了相關研究[10-12]。

本文首先詳細介紹了DSMC數值模擬方法的相關理論,然后對DSMC方法進行了標模驗證;最后將DSMC方法應用于増阻離軌裝置稀薄流場下的氣動特性預測。

2 DSMC數值模擬方法

DSMC方法是基于氣體動力學理論發展而來的數值仿真方法,DSMC方法被嚴格證明是與玻爾茲曼方程等價的[11]。DSMC方法基本原理是在計算機中用大量模擬分子代表真實氣體,模擬分子的空間坐標、速度、內能等存儲在計算機中,因分子運動、碰撞以及與邊界的相互作用而隨時間變化。該方法的核心思想是在一個時間步長內,將分子運動與碰撞解耦,所有分子首先按照各自的速度,運動一段距離,然后再按照與真實分子相同的碰撞頻率,根據一定的碰撞模型,從所有可能組合的分子對中選出碰撞對,以幾率的方式分配碰撞后的分子速度和內能、決定是否有化學反應發生等。

下面主要對DSMC中分子與物面相互作用模型、分子碰撞模型以及碰撞對選取等進行介紹。

2.1 分子與物面相互作用模型

稀薄流氣體分子與物面相互作用過程的模擬,直接影響DSMC方法對氣動力、熱的計算精度。本文采用鏡面反射和漫反射兩種相互作用模型。其中鏡面反射假設氣體分子與物面發生彈性碰撞;漫反射模型假設氣體分子與物面發生非彈性碰撞,且反射后的氣體分子向各個方向散射,散射后的分子速度滿足平衡的Maxwell分布,該分布只與反射后分子溫度有關,與入射分子的動量速度無關。入射分子的動量取決于自由來流的條件,反射分子的動量特性用動量協調系數來刻畫,不同的協調系數取值則代表了不同的碰撞模型。法向動量協調系數σN表示分子與面碰撞過程中法向動量的改變,切向動量協調系數σT則表示分子與面碰撞過程中切向動量的改變,計算公式如下:

式中,pi、 τi和pr、 τr分別為入射分子與反射分子的法向壓力和切向壓力;pw為壁面法向壓力。由此可知,發生鏡面反射碰撞時有σN=σT=0,發生漫反射碰撞時則有σN=σT=1。

2.2 分子碰撞模型

稀薄流場條件下的氣體分子滿足三體碰撞不重要條件d3/δ3?1,即分子之間發生三體或者以上碰撞的概率要遠小于雙體碰撞。因此在考慮分子碰撞模型時,僅考慮雙體碰撞。此時,利用動量守恒和能量守恒定律可以得到:

式中,b是分子間距;d是分子直徑。得到偏轉角以后,結合公式 (3),可求得碰撞后的分子速度,且速度方向在各個方向均勻散布。

2.3 分子碰撞對選取

分子碰撞對的選取常有的有兩種方法:時間計數器 (Time Counter,TC)法、無時間計數器 (No Time Counter,NTC)法。本文中的分子碰撞對選取采用常用的無時間計數器 (No Time Counter,NTC)法[13]。用假設的單個模擬分子來代替的一定數量的真實分子,則計算碰撞對發生概率P時,僅需考慮所有碰撞對的一部分并同時將P按相同的倍數加以放大。倘若該倍數可使得最大碰撞概率為1,則可以理解為計算效率最高。可以得到,需要做計算處理的碰撞概率P的表達式為:

式中,σT是分子碰撞截面;cr是兩個分子的相對速度。

2.4 仿真程序簡介

本文采用由Sandia實驗室開發的稀薄流計算開源代碼SPARTA進行分析計算。SPARTA基于DSMC算法,采用笛卡爾網格對粒子進行跟蹤和分組;其包含多種粒子碰撞模型,化學反應模型以及壁面邊界條件模型。近年來已經在很多工程實踐和學術研究中得到了充分的驗證和應用。

DSMC方法在SPARTA程序中實現流場如圖1所示。SPARTA程序采用的是笛卡爾網格,程序開始后需要初始化網格信息,包括粒子數目、網格節點數目及邊界條件等;然后計算粒子在Δt時間步內的運動及粒子與物面的相互作用;更新粒子索引信息,將分子重新排序和編號;按照規定碰撞模型進行粒子間碰撞計算;流場取樣;判斷是否達到計算時間;最后計算平均流場信息:非定常流場對各個重復過程進行系綜平均,定常流動在定常化后求足夠大時間的時間平均。

圖1 DSMC實現流程框圖Fig.1 Flow chart of DSMC implementation

3 SPARTA程序標模驗證

NASA基于DSMC開發出一套行業標準軟件工具DAC(DSMC Analysis Code),該軟件對于航天器在自由分子流區的氣動分析具有很高的計算精度。

本文采用SPARTA程序對球算例進行了自由分子流域內的氣動性能進行了模擬仿真,并與DAC軟件結果進行了對比。

數值模擬工況如表1所示。

表1 球算例計算工況Table 1 Calculation conditions for ball calculation examples

圖2分別給出了球面時均壓力與對稱面時均速度云圖。從圖中可以看出在球的迎風面存在一個高壓區,壓力量級約10-5,在背風區由于來流粒子與球面動量交換較少,因此壓力明顯較低。

圖2 球表面時均壓力云圖 (a)與對稱面時均速度云圖 (b)Fig.2 Hourly average pressure on the surface of the ball(a)and time-averaged velocity distribution cloud diagram of symmetry plane(b)

從圖2中可以看出,粒子速度在遠場為來流速,隨著粒子靠近球面,粒子與球面發生碰撞完成動量交換后,粒子速度明顯降低;粒子經過球面,在球面背風區粒子數目相對較少,接近真空環境。

圖3給出了SPARTA和DAC軟件不同工況阻力系數對比曲線。不同工況下由文中計算得出的阻力系數與DAC計算得出結果吻合程度很好,其中鏡面反射碰撞模型的最大計算誤差僅為0.53%。由此可知文中采用的計算方法和程序可對自由分子流區航天器增阻面的阻力系數給出可靠的計算結果。

圖3 SPARTA與DAC軟件球阻力系數對比Fig.3 Comparison of ball drag coefficient between SPARTA and DAC software

4 增阻離軌裝置稀薄流氣動特性預測

增阻離軌空間碎片清除技術是一項適用于近地軌道區域任務后航天器離軌的新技術,對于軌道高度小于850km的任務后航天器,通過增阻離軌系統形成很大迎風面,增加大氣阻力,從而逐漸降低并脫離運行軌道,在一定時間內隕落,再入大氣層燒毀。該技術可用于廢棄衛星、微小衛星、廢棄的運載火箭上面級的離軌,應用前景廣闊。

本文將SPARTA程序應用于增阻離軌裝置在稀薄流域內的氣動特性數值預測。仿真工況如下:

表2 計算工況Table 2 Calculation conditions

圖4給出了物面壓力和對稱面速度云圖分布,與球算例類似,迎風面存在高壓區;物面附近速度降低,物面后存在近似完全真空區域。

圖4 增阻離軌裝置物面時均壓力云圖 (a)與對稱面時均速度云圖 (b)Fig.4 Time-averaged pressure cloud diagram of the object surface of the drag-increasing de-orbit device(a)and time-averaged velocity cloud diagram of symmetry plane(b)

圖5為增阻離軌裝置在不同攻角下阻力系數曲線,由圖中可以看出攻角的變化對漫反射碰撞模型的氣動阻力系數影響較大,且阻力系數隨攻角增大有顯著減小的變化趨勢;當攻角在0°～30°之間變化時,攻角變化對鏡面反射碰撞模型的氣動系數的影響很小,而當攻角在30°～90°之間變化時,氣動阻力系數對攻角的變化較為敏感,會隨攻角的增大有顯著快速衰減的變化趨勢。

圖5 增阻離軌裝置在不同攻角下阻力系數曲線圖Fig.5 The drag efficient curve of the drag-enhancing de-orbit device at different angles of attack

不同軌高處的稀薄大氣在分子的種類、數密度及來流溫度等方面存在一定差異,表3給出了7.5km/s、1000K來流溫度、0°攻角條件下,不同軌高處氣動阻力系數計算結果。

表3 不同軌道阻力系數Table 3 Drag coefficient of different tracks

結果表明,阻力系數隨軌道高度的變化不明顯,但由于不同軌道高度處大氣密度差距較大,不同軌高處的氣動阻力差距也較大。

5 結論

本文首先通過SPARTA程序對球標模進行了模擬并與DAC軟件結果進行了對比,兩者誤差較小;然后將SPARTA程序應用于增阻離軌裝置在不同攻角工況氣動特性仿真,結果表明在0°～30°攻角下,增阻離軌裝置隨著攻角增加阻力系數下降較緩,而當攻角在30°～90°時,阻力系數下降明顯,軌道高度對氣動阻力系數的影響較小。