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

(1.東南大學(xué)國家預(yù)應(yīng)力工程技術(shù)研究中心,南京 211189;2.北京空間機(jī)電研究所,北京 100094;3.中國航天科技集團(tuán)有限公司航天進(jìn)入、減速與著陸技術(shù)實(shí)驗(yàn)室,北京 100094)

1 引言

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

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

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

本文首先詳細(xì)介紹了DSMC數(shù)值模擬方法的相關(guān)理論,然后對DSMC方法進(jìn)行了標(biāo)模驗(yàn)證;最后將DSMC方法應(yīng)用于増阻離軌裝置稀薄流場下的氣動特性預(yù)測。

2 DSMC數(shù)值模擬方法

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

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

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

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

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

2.2 分子碰撞模型

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

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

2.3 分子碰撞對選取

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

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

2.4 仿真程序簡介

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

DSMC方法在SPARTA程序中實(shí)現(xiàn)流場如圖1所示。SPARTA程序采用的是笛卡爾網(wǎng)格,程序開始后需要初始化網(wǎng)格信息,包括粒子數(shù)目、網(wǎng)格節(jié)點(diǎn)數(shù)目及邊界條件等;然后計算粒子在Δt時間步內(nèi)的運(yùn)動及粒子與物面的相互作用;更新粒子索引信息,將分子重新排序和編號;按照規(guī)定碰撞模型進(jìn)行粒子間碰撞計算;流場取樣;判斷是否達(dá)到計算時間;最后計算平均流場信息:非定常流場對各個重復(fù)過程進(jìn)行系綜平均,定常流動在定?；笄笞銐虼髸r間的時間平均。

圖1 DSMC實(shí)現(xiàn)流程框圖Fig.1 Flow chart of DSMC implementation

3 SPARTA程序標(biāo)模驗(yàn)證

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

本文采用SPARTA程序?qū)η蛩憷M(jìn)行了自由分子流域內(nèi)的氣動性能進(jìn)行了模擬仿真,并與DAC軟件結(jié)果進(jìn)行了對比。

數(shù)值模擬工況如表1所示。

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

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

圖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中可以看出,粒子速度在遠(yuǎn)場為來流速,隨著粒子靠近球面,粒子與球面發(fā)生碰撞完成動量交換后,粒子速度明顯降低;粒子經(jīng)過球面,在球面背風(fēng)區(qū)粒子數(shù)目相對較少,接近真空環(huán)境。

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

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

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

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

本文將SPARTA程序應(yīng)用于增阻離軌裝置在稀薄流域內(nèi)的氣動特性數(shù)值預(yù)測。仿真工況如下:

表2 計算工況Table 2 Calculation conditions

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

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

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

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

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

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

5 結(jié)論

本文首先通過SPARTA程序?qū)η驑?biāo)模進(jìn)行了模擬并與DAC軟件結(jié)果進(jìn)行了對比,兩者誤差較小;然后將SPARTA程序應(yīng)用于增阻離軌裝置在不同攻角工況氣動特性仿真,結(jié)果表明在0°～30°攻角下,增阻離軌裝置隨著攻角增加阻力系數(shù)下降較緩,而當(dāng)攻角在30°～90°時,阻力系數(shù)下降明顯,軌道高度對氣動阻力系數(shù)的影響較小。