基于DSMC方法的超低軌道衛(wèi)星阻力特性分析

杭曉晨 李彥斌 孔祥宏 陳 強(qiáng) 費(fèi)慶國(guó)

(1南京林業(yè)大學(xué)機(jī)電工程學(xué)院, 南京 210037)

(2東南大學(xué)空天機(jī)械動(dòng)力學(xué)研究所, 南京 210096)

(3上海衛(wèi)星工程研究所, 上海 201109)

超低軌道衛(wèi)星運(yùn)行于180～300 km高度地球軌道.相比于傳統(tǒng)低軌道,超低軌道衛(wèi)星偵查地面的光學(xué)成像分辨率成倍提高,雷達(dá)功耗降低,使得衛(wèi)星小型化成為可能,從而大大降低了衛(wèi)星的研制成本.超低軌道衛(wèi)星在空間信息對(duì)抗方面具有優(yōu)勢(shì),美國(guó)、歐洲、日本等國(guó)家相繼對(duì)此類衛(wèi)星投入了大量研究[1-3].衛(wèi)星在超低軌道上處于稀薄氣體大氣環(huán)境,氣動(dòng)阻力不可忽略,經(jīng)常需要離子發(fā)動(dòng)機(jī)進(jìn)行速度補(bǔ)償,以保證衛(wèi)星的軌道機(jī)動(dòng)和姿態(tài)控制.因此,研究超低軌道衛(wèi)星氣動(dòng)阻力特性和精確預(yù)測(cè)方法,對(duì)該類衛(wèi)星的合理結(jié)構(gòu)設(shè)計(jì)具有重要意義.

180～300 km軌道高度屬于稀薄氣體環(huán)境,氣體分子間距大,克努森數(shù)高,不滿足傳統(tǒng)航空領(lǐng)域的氣體連續(xù)性假設(shè),基于納維-斯托克斯方程的連續(xù)流體分析方法不再適用.目前,預(yù)測(cè)計(jì)算衛(wèi)星阻力的方法主要有工程近似算法和數(shù)值仿真方法.Titov等[4]提出了低軌道衛(wèi)星阻力計(jì)算的近似表達(dá)式,認(rèn)為阻力與大氣密度、迎風(fēng)面投影面積以及來(lái)流速度平方呈正比關(guān)系.該表達(dá)式能夠在衛(wèi)星結(jié)構(gòu)設(shè)計(jì)初期階段提供快速的阻力近似預(yù)測(cè).數(shù)值仿真方法包括面元積分法[5]、格子玻爾茲曼法[6]、直接模擬蒙特卡洛(direct simulation Monte Carlo, DSMC)方法[7]、試驗(yàn)粒子蒙特卡洛法[8]等.其中,DSMC方法是一種基于分子動(dòng)力學(xué)的改進(jìn)數(shù)值算法,利用概率統(tǒng)計(jì)判斷分子間是否發(fā)生碰撞,相比于傳統(tǒng)分子動(dòng)力學(xué)仿真可以極大降低計(jì)算資源消耗.

Moe等[9]采用DSMC方法研究了稀薄氣體環(huán)境分子-壁面碰撞模型,認(rèn)為能量調(diào)節(jié)系數(shù)兼顧了不同的氣動(dòng)作用,并估算了阻力系數(shù)的不確定性.Mehta等[10]利用DSMC方法仿真分析了漫反射模型和準(zhǔn)鏡面模型對(duì)衛(wèi)星阻力特性的影響,發(fā)現(xiàn)2種壁面碰撞模型在GRACE衛(wèi)星上的阻力差異約為15%.DSMC方法也被用于分析航天器再入大氣層時(shí)的稀薄大氣環(huán)境,如文獻(xiàn)[11]研究了航天器再入階段的氣動(dòng)特性和動(dòng)力學(xué)穩(wěn)定性,文獻(xiàn)[12]研究了氣動(dòng)加熱效應(yīng)等.李志輝等[13]將DSMC方法應(yīng)用于模擬阿波羅指令艙稀薄氣體動(dòng)力學(xué)特征中.周偉勇等[14]采用部分計(jì)算總體疊加的思路,分析了復(fù)雜外形衛(wèi)星的氣動(dòng)力,給出了衛(wèi)星減阻的設(shè)計(jì)建議.黃飛等[15]基于DSMC方法分析了GOCE衛(wèi)星的氣動(dòng)特性,討論了不同物面反射系數(shù)對(duì)阻力特性的影響.靳旭紅等[16-18]采用試驗(yàn)粒子蒙特卡洛方法,分析了大氣模型、飛行高度、軌道緯度等對(duì)超低軌道衛(wèi)星阻力特性的影響規(guī)律.

本文采用DSMC方法對(duì)一典型超低軌道衛(wèi)星進(jìn)行了三維分子動(dòng)力學(xué)仿真,預(yù)測(cè)了衛(wèi)星阻力特性,研究了衛(wèi)星不同幾何參數(shù)對(duì)氣動(dòng)特性的影響規(guī)律.

1 DSMC方法

圖1 VHS二元碰撞模型示意圖

(1)

σT=πd2

(2)

式中,碰撞參數(shù)b為質(zhì)心參考系中未擾動(dòng)軌跡的最近距離,可由b=1/2(d1+d2)sinθA計(jì)算得到,其中θA為相對(duì)速度與球心連線的夾角;d為可變硬球分子的直徑參數(shù),是氣體分子碰撞對(duì)中相對(duì)速度cr的逆冪律形式函數(shù).

稀薄流氣體分子與物面的相互作用直接影響DSMC方法對(duì)阻力特性的計(jì)算精度.針對(duì)衛(wèi)星結(jié)構(gòu)中與來(lái)流平行或小夾角的壁面,采用鏡面反射模型,即假設(shè)氣體分子與物面產(chǎn)生彈性碰撞.針對(duì)衛(wèi)星結(jié)構(gòu)中與來(lái)流垂直或大夾角的壁面,采用漫反射模型,即假設(shè)氣體分子與物面發(fā)生非彈性碰撞,且反射后的氣體分子向各個(gè)方向散射,散射時(shí)的分子速度滿足平衡的Maxwell分布.采用動(dòng)量協(xié)調(diào)系數(shù)來(lái)描述反射分子的動(dòng)量特性改變,法向動(dòng)量協(xié)調(diào)系數(shù)σn表示分子與物面碰撞過(guò)程的法向動(dòng)量改變,切向動(dòng)量協(xié)調(diào)系數(shù)στ表示切向動(dòng)量改變,計(jì)算公式分別為

(3)

(4)

式中,pi、τi分別為入射分子的法向壓力和切向壓力;pr、τr分別為反射分子的法向壓力和切向壓力;pw為壁面法向壓力.發(fā)生鏡面碰撞時(shí),σn=στ=0;發(fā)生漫反射碰撞時(shí),σn=στ=1.

DSMC方法將抽樣的氣體分子以概率計(jì)數(shù)的方式來(lái)表示真實(shí)環(huán)境中大量的真實(shí)分子,從而降低計(jì)算量.本文采用非時(shí)間計(jì)數(shù)法(no time counter,NTC)[7],在一個(gè)DSMC網(wǎng)格內(nèi),選取潛在可能碰撞對(duì),當(dāng)其碰撞概率P大于生成的隨機(jī)數(shù)R(0

2 低軌衛(wèi)星阻力特性

以一典型超低軌道衛(wèi)星為研究對(duì)象,軌道高度為268 km,環(huán)境大氣密度為55.09 ng/m3,環(huán)境溫度為951.5 K,飛行速度為7 740 m/s.基于DSMC方法研究衛(wèi)星的阻力特性,分析結(jié)構(gòu)構(gòu)型對(duì)阻力特性的影響.

衛(wèi)星幾何模型如圖2所示.衛(wèi)星本體柱段需容納一定體積的載荷,兩側(cè)太陽(yáng)翼總面積約為10.2 m2.以柱體為長(zhǎng)方體、橫截面為矩形、柱體長(zhǎng)3 m為基準(zhǔn)模型,基于DSMC方法分析衛(wèi)星在軌阻力特性,并研究網(wǎng)格密度、頭錐形狀、長(zhǎng)細(xì)比等參數(shù)對(duì)阻力特性的影響規(guī)律.

圖2 典型超低軌道衛(wèi)星幾何模型(單位:mm)

2.1 網(wǎng)格密度對(duì)DSMC結(jié)果收斂性的影響

采用DSMC方法仿真稀薄氣體時(shí),衛(wèi)星結(jié)構(gòu)的網(wǎng)格過(guò)疏可能會(huì)導(dǎo)致對(duì)物面的曲面模擬失真,計(jì)算結(jié)果精度較低;過(guò)密的網(wǎng)格則會(huì)因?yàn)槟M大量分子壁面碰撞而加大計(jì)算資源消耗.以衛(wèi)星的圓形平板頭錐和半球形頭錐為例,在200 km軌道大氣環(huán)境下,研究網(wǎng)格密度對(duì)DSMC方法計(jì)算氣動(dòng)阻力結(jié)果的影響.

圖3給出了不同網(wǎng)格密度設(shè)計(jì)方案.圓形平板直徑為1 m,厚度為10 mm,低、中、高密度網(wǎng)格分別包含524、960、1 960個(gè)三角形網(wǎng)格,對(duì)應(yīng)的網(wǎng)格特征長(zhǎng)度分別為100、50、30 mm,相對(duì)尺度分別為10%、5%、3%.半球形頭錐底面外接圓直徑為1 m,低、中、高密度網(wǎng)格分別包含500、1 000、2 000個(gè)三角形網(wǎng)格,對(duì)應(yīng)的網(wǎng)格特征長(zhǎng)度分別為110、75、55 mm.不同網(wǎng)格密度模型受到的氣動(dòng)阻力見(jiàn)表1.

(a) 低密度

(b) 中密度

(c) 高密度

表1 不同網(wǎng)格密度模型受到的氣動(dòng)阻力

圖4給出了半球形頭錐的DSMC分析結(jié)果.由圖可知,迎風(fēng)面中心處壓強(qiáng)最大,越靠外緣則壓強(qiáng)越小.半球附近存在高壓區(qū),這是因?yàn)闅怏w分子在頭錐表面撞擊后獲得反向動(dòng)量,導(dǎo)致區(qū)域內(nèi)氣體分子密度增加.由表1可知,低密度網(wǎng)格與中密度網(wǎng)格的氣動(dòng)阻力差在0.3 mN以內(nèi),中密度網(wǎng)格與高密度網(wǎng)格的計(jì)算結(jié)果幾乎無(wú)差別,故后續(xù)數(shù)值分析均采用中等以上的網(wǎng)格密度方案,即在對(duì)衛(wèi)星表面劃分網(wǎng)格時(shí),采取網(wǎng)格相對(duì)尺度在3%～5%范圍內(nèi)值取.

(a) 表面壓強(qiáng)場(chǎng)

(b) 全域壓強(qiáng)場(chǎng)

2.2 頭錐幾何參數(shù)對(duì)阻力特性的影響

本節(jié)研究不同的衛(wèi)星頭錐形狀對(duì)受到氣動(dòng)阻力的影響.考慮錐度為60°、90°、120°的四棱錐、六棱錐、八棱錐、圓錐以及半球形頭錐和圓形平板頭錐總共14個(gè)分析模型(見(jiàn)圖5).實(shí)際工程中結(jié)構(gòu)設(shè)計(jì)必須考慮星體的載荷尺寸約束,將所有頭錐模型的尺寸均設(shè)為底面投影外接直徑為1 m的圓.

(a) 錐角60°

(b) 錐角90°

(c) 錐角120°

不同的頭錐模型在平飛(θ=0°)、豎飛(θ=90°)兩種飛行姿態(tài)下的氣動(dòng)阻力計(jì)算結(jié)果見(jiàn)表2.由表可知,圓形平板頭錐、半球形頭錐、圓錐在平飛姿態(tài)下具有完全相同的迎風(fēng)面投影面積0.785 m2,所受氣動(dòng)阻力相差不大,均為18.50～18.95 mN,故迎風(fēng)面投影面積是影響氣動(dòng)阻力的主要因素.四棱錐、六棱錐、八棱錐等其他模型底面投影外接圓,故迎風(fēng)面積均小于0.785 m2,受到的阻力幅值也較小.圖6給出了頭錐模型的迎風(fēng)面與阻力散點(diǎn)圖.由圖可知,頭錐受到的阻力與迎風(fēng)面積具有線性相關(guān)性.

圖6 不同頭錐模型的阻力特性

同時(shí)表2給出了不同頭錐模型平飛姿態(tài)時(shí)的阻力系數(shù)CD,圓形平板頭錐的阻力系數(shù)最大為2.13,60°圓錐阻力系數(shù)最小為2.08.在相同的迎風(fēng)面投影面積情況下,平飛姿態(tài)下頭錐的錐角越小,頭錐越尖,受到阻力越小,但錐角不同導(dǎo)致的平飛阻力差異在1 mN范圍內(nèi).在豎飛姿態(tài)下,錐角越小,迎風(fēng)面越大,因此豎飛時(shí)阻力水平越高.例如,60°錐角的圓錐平飛阻力為18.49 mN,為同一迎風(fēng)面積下的最低阻力值,但其豎飛時(shí)阻力為10.77 mN,阻力值較大.究其原因在于,在氣體粒子與物面的碰撞中,頭錐越尖,壁面與來(lái)流的夾角越小;氣體粒子撞擊壁面的過(guò)程中,法向動(dòng)量改變?cè)叫?切向動(dòng)量改變?cè)酱?即壓差阻力占比較小,摩擦阻力占比較大.

表2 不同頭錐阻力特性結(jié)果

為進(jìn)一步揭示壓差阻力和摩擦阻力的關(guān)系,研究了不同高度梯形圓錐模型的氣動(dòng)阻力變化(見(jiàn)圖7),與基準(zhǔn)的60°錐角圓錐模型進(jìn)行對(duì)比,結(jié)果見(jiàn)表3.顯然,梯形圓錐高度越低,平飛姿態(tài)下頭錐將受到越多來(lái)流方向氣體分子的撞擊,摩擦分量占比越小,總阻力略微升高.

圖7 梯形圓錐高度影響分析模型

表3 梯形圓錐阻力特性 mN

2.3 衛(wèi)星星體長(zhǎng)細(xì)比對(duì)阻力特性的影響

衛(wèi)星星體長(zhǎng)細(xì)比也是初期結(jié)構(gòu)設(shè)計(jì)時(shí)必須要考慮的重要因素.本節(jié)研究了不同結(jié)構(gòu)星體的阻力特性隨長(zhǎng)細(xì)比的變化規(guī)律.考慮載荷需要,衛(wèi)星模型星體橫截面均設(shè)計(jì)為外接直徑為1 m的圓,星體形狀分別為四棱柱、六棱柱、八棱柱與基準(zhǔn)圓柱,利用DSMC方法分析得到阻力特性,結(jié)果見(jiàn)表4.同時(shí)也給出了圓柱星體帶圓錐頭錐結(jié)構(gòu)所受的氣動(dòng)阻力值,對(duì)比分析了同投影面積下頭錐對(duì)氣動(dòng)阻力特性的影響.

表4 不同長(zhǎng)細(xì)比星體結(jié)構(gòu)阻力特性 mN

由表4可知,對(duì)于同樣的星體橫截面,在迎風(fēng)面投影面積相同的情況下,長(zhǎng)細(xì)比越大,星體所受阻力也越大.究其原因在于,DSMC仿真中,大量的氣體分子與星體側(cè)面發(fā)生小入射角的摩擦型碰撞,星體長(zhǎng)細(xì)比越大,氣體分子與星體側(cè)面接觸面積越大,引起的摩擦阻力分量也越大.針對(duì)本文采用的衛(wèi)星模型,長(zhǎng)細(xì)比每增加1,阻力增加約5%～8%.

2.4 基于DSMC方法的整星氣動(dòng)特性分析

利用歐洲GOCE衛(wèi)星算例,驗(yàn)證本文提出的DSMC方法.GOCE衛(wèi)星是世界上第1顆重力梯度測(cè)量衛(wèi)星,主體由衛(wèi)星本體、太陽(yáng)翼和穩(wěn)定翼構(gòu)成.星體為八棱柱,太陽(yáng)翼面積較大,穩(wěn)定翼面積較小,其分析模型如圖8所示.

分析域?yàn)?0 m×10 m×10 m的立方體.根據(jù)網(wǎng)格收斂性結(jié)論,將GOCE衛(wèi)星表面劃分為6 980個(gè)網(wǎng)格,網(wǎng)格相對(duì)尺度為3%.分析結(jié)果如圖9所示,GOCE衛(wèi)星的六邊形頭部壓強(qiáng)最大,太陽(yáng)翼、尾翼與來(lái)流接觸面呈現(xiàn)較大壓強(qiáng).由全場(chǎng)壓強(qiáng)云圖可以看出,氣體分子密度在衛(wèi)星頭部附近最高;相對(duì)于遠(yuǎn)處流場(chǎng),氣體分子在尾翼附近、星體側(cè)面更為聚集.該衛(wèi)星在268 km軌道大氣環(huán)境下,平飛時(shí)所受氣動(dòng)阻力為8.8 mN.將本文方法獲得的GOCE衛(wèi)星在典型軌道200～268 km的氣動(dòng)阻力結(jié)果與文獻(xiàn)[17]中GOCE衛(wèi)星在不同大氣模型下的阻力曲線進(jìn)行對(duì)比,對(duì)比結(jié)果見(jiàn)圖10.由圖可知,本文方法結(jié)果與文獻(xiàn)[17]中采用多種大氣模型的計(jì)算結(jié)果吻合良好.

(a) 表面壓強(qiáng)場(chǎng)

(b) 全域壓強(qiáng)場(chǎng)

圖10 GOCE衛(wèi)星阻力特性與文獻(xiàn)值對(duì)比

針對(duì)典型超低軌道衛(wèi)星結(jié)構(gòu)(見(jiàn)圖2和圖11),研究整機(jī)在不同姿態(tài)和攻角飛行時(shí)所受的阻力,并分析平飛姿態(tài)下全動(dòng)尾翼轉(zhuǎn)動(dòng)對(duì)整體氣動(dòng)特性的影響.全動(dòng)的水平尾翼為梯形,安裝于星體尾部高于整星中軸線40 mm處,面積為0.2 m2.分析結(jié)果見(jiàn)表5.

圖11 衛(wèi)星尾翼位置(單位:mm)

由表5可知,與200 km軌道的頭錐、星體算例結(jié)果相比,GOCE衛(wèi)星每單位投影面積所受到的空氣阻力明顯減小.究其原因在于,268 km軌道的空氣更為稀薄,每立方米空氣分子數(shù)為1.144×1015,遠(yuǎn)小于200 km軌道中每立方米空氣分子數(shù)(5.8×1015),相應(yīng)的氣體分子碰撞次數(shù)和碰撞概率均降低,因此阻力較小.

表5 衛(wèi)星整機(jī)與尾翼多攻角阻力特性

全動(dòng)尾翼偏轉(zhuǎn)角會(huì)加大迎風(fēng)面積,攻角為10°時(shí)附加阻力為0.09 mN,升力為0.031 mN,升阻比為0.34;攻角為20°時(shí),附加阻力為0.35 mN,升力為0.041 mN,升阻比為0.12;攻角為30°時(shí),附加阻力為0.56 mN,升力為0.053 mN,升阻比為 0.10.因此,在稀薄氣體自由分子流情況下,升阻比隨攻角增加呈下降趨勢(shì);而在連續(xù)流環(huán)境下,隨攻角增加,升阻比先增大后減小.由此可見(jiàn),在升力特性方面,超低軌道的稀薄氣體環(huán)境與連續(xù)流環(huán)境具有較大差異.尾翼偏轉(zhuǎn)角由0°增加到30°時(shí),升力逐漸增大,升阻比在0.10～0.34范圍內(nèi)先增大后減小.

3 結(jié)論

1) 基于DSMC方法分析了典型超低軌道衛(wèi)星的阻力特性,采用小攻角鏡面反射、大攻角Maxwell漫反射聯(lián)合模型,研究了衛(wèi)星結(jié)構(gòu)幾何構(gòu)型對(duì)阻力特性的影響規(guī)律.結(jié)果表明,迎風(fēng)面投影面積是影響衛(wèi)星阻力的主要因素.在相同的迎風(fēng)面投影面積下,較尖的頭錐結(jié)構(gòu)能略微降低平飛阻力.

2) 星體長(zhǎng)細(xì)比對(duì)衛(wèi)星阻力有較大影響.相同的投影面積情況下,星體越長(zhǎng),衛(wèi)星側(cè)面與氣體分子的接觸面積越大,產(chǎn)生切向的摩擦阻力越大.對(duì)于本研究中的衛(wèi)星,長(zhǎng)細(xì)比每增加1,阻力增加約5%～8%.

3) 設(shè)計(jì)的衛(wèi)星結(jié)構(gòu)在268 km典型超低軌道下,尾翼偏轉(zhuǎn)角由0°增加到30°時(shí),升力逐漸增大,升阻比在0.10～0.34范圍內(nèi)先增大后減小.

4) 對(duì)于GOCE衛(wèi)星,本文分析得到的阻力特性結(jié)果與文獻(xiàn)[17]中多種大氣模型結(jié)果吻合良好,驗(yàn)證了DSMC方法在計(jì)算超低軌道衛(wèi)星阻力特性時(shí)的準(zhǔn)確性.