增阻離軌裝置氣動(dòng)特性的DSMC研究

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

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

1 引言

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

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

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

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

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

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

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

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

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

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

2.2 分子碰撞模型

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

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

2.3 分子碰撞對(duì)選取

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

式中,σT是分子碰撞截面;cr是兩個(gè)分子的相對(duì)速度。

2.4 仿真程序簡介

本文采用由Sandia實(shí)驗(yàn)室開發(fā)的稀薄流計(jì)算開源代碼SPARTA進(jìn)行分析計(jì)算。SPARTA基于DSMC算法,采用笛卡爾網(wǎng)格對(duì)粒子進(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ù)目及邊界條件等;然后計(jì)算粒子在Δt時(shí)間步內(nèi)的運(yùn)動(dòng)及粒子與物面的相互作用;更新粒子索引信息,將分子重新排序和編號(hào);按照規(guī)定碰撞模型進(jìn)行粒子間碰撞計(jì)算;流場取樣;判斷是否達(dá)到計(jì)算時(shí)間;最后計(jì)算平均流場信息:非定常流場對(duì)各個(gè)重復(fù)過程進(jìn)行系綜平均,定常流動(dòng)在定?；笄笞銐虼髸r(shí)間的時(shí)間平均。

圖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),該軟件對(duì)于航天器在自由分子流區(qū)的氣動(dòng)分析具有很高的計(jì)算精度。

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

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

表1 球算例計(jì)算工況Table 1 Calculation conditions for ball calculation examples

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

圖2 球表面時(shí)均壓力云圖 (a)與對(duì)稱面時(shí)均速度云圖 (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ā)生碰撞完成動(dòng)量交換后,粒子速度明顯降低;粒子經(jīng)過球面,在球面背風(fēng)區(qū)粒子數(shù)目相對(duì)較少,接近真空環(huán)境。

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

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

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

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

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

表2 計(jì)算工況Table 2 Calculation conditions

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

圖4 增阻離軌裝置物面時(shí)均壓力云圖 (a)與對(duì)稱面時(shí)均速度云圖 (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ù)曲線,由圖中可以看出攻角的變化對(duì)漫反射碰撞模型的氣動(dòng)阻力系數(shù)影響較大,且阻力系數(shù)隨攻角增大有顯著減小的變化趨勢(shì);當(dāng)攻角在0°～30°之間變化時(shí),攻角變化對(duì)鏡面反射碰撞模型的氣動(dòng)系數(shù)的影響很小,而當(dāng)攻角在30°～90°之間變化時(shí),氣動(dòng)阻力系數(shù)對(duì)攻角的變化較為敏感,會(huì)隨攻角的增大有顯著快速衰減的變化趨勢(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°攻角條件下,不同軌高處氣動(dòng)阻力系數(shù)計(jì)算結(jié)果。

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

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

5 結(jié)論

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