基于DSMC方法的超低軌道衛星阻力特性分析

杭曉晨 李彥斌 孔祥宏 陳 強 費慶國

(1南京林業大學機電工程學院, 南京 210037)

(2東南大學空天機械動力學研究所, 南京 210096)

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

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

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

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

本文采用DSMC方法對一典型超低軌道衛星進行了三維分子動力學仿真,預測了衛星阻力特性,研究了衛星不同幾何參數對氣動特性的影響規律.

1 DSMC方法

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

(1)

σT=πd2

(2)

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

稀薄流氣體分子與物面的相互作用直接影響DSMC方法對阻力特性的計算精度.針對衛星結構中與來流平行或小夾角的壁面,采用鏡面反射模型,即假設氣體分子與物面產生彈性碰撞.針對衛星結構中與來流垂直或大夾角的壁面,采用漫反射模型,即假設氣體分子與物面發生非彈性碰撞,且反射后的氣體分子向各個方向散射,散射時的分子速度滿足平衡的Maxwell分布.采用動量協調系數來描述反射分子的動量特性改變,法向動量協調系數σn表示分子與物面碰撞過程的法向動量改變,切向動量協調系數στ表示切向動量改變,計算公式分別為

(3)

(4)

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

DSMC方法將抽樣的氣體分子以概率計數的方式來表示真實環境中大量的真實分子,從而降低計算量.本文采用非時間計數法(no time counter,NTC)[7],在一個DSMC網格內,選取潛在可能碰撞對,當其碰撞概率P大于生成的隨機數R(0

2 低軌衛星阻力特性

以一典型超低軌道衛星為研究對象,軌道高度為268 km,環境大氣密度為55.09 ng/m3,環境溫度為951.5 K,飛行速度為7 740 m/s.基于DSMC方法研究衛星的阻力特性,分析結構構型對阻力特性的影響.

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

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

2.1 網格密度對DSMC結果收斂性的影響

采用DSMC方法仿真稀薄氣體時,衛星結構的網格過疏可能會導致對物面的曲面模擬失真,計算結果精度較低;過密的網格則會因為模擬大量分子壁面碰撞而加大計算資源消耗.以衛星的圓形平板頭錐和半球形頭錐為例,在200 km軌道大氣環境下,研究網格密度對DSMC方法計算氣動阻力結果的影響.

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

(a) 低密度

(b) 中密度

(c) 高密度

表1 不同網格密度模型受到的氣動阻力

圖4給出了半球形頭錐的DSMC分析結果.由圖可知,迎風面中心處壓強最大,越靠外緣則壓強越小.半球附近存在高壓區,這是因為氣體分子在頭錐表面撞擊后獲得反向動量,導致區域內氣體分子密度增加.由表1可知,低密度網格與中密度網格的氣動阻力差在0.3 mN以內,中密度網格與高密度網格的計算結果幾乎無差別,故后續數值分析均采用中等以上的網格密度方案,即在對衛星表面劃分網格時,采取網格相對尺度在3%～5%范圍內值取.

(a) 表面壓強場

(b) 全域壓強場

2.2 頭錐幾何參數對阻力特性的影響

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

(a) 錐角60°

(b) 錐角90°

(c) 錐角120°

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

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

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

表2 不同頭錐阻力特性結果

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

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

表3 梯形圓錐阻力特性 mN

2.3 衛星星體長細比對阻力特性的影響

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

表4 不同長細比星體結構阻力特性 mN

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

2.4 基于DSMC方法的整星氣動特性分析

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

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

(a) 表面壓強場

(b) 全域壓強場

圖10 GOCE衛星阻力特性與文獻值對比

針對典型超低軌道衛星結構(見圖2和圖11),研究整機在不同姿態和攻角飛行時所受的阻力,并分析平飛姿態下全動尾翼轉動對整體氣動特性的影響.全動的水平尾翼為梯形,安裝于星體尾部高于整星中軸線40 mm處,面積為0.2 m2.分析結果見表5.

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

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

表5 衛星整機與尾翼多攻角阻力特性

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

3 結論

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

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

3) 設計的衛星結構在268 km典型超低軌道下,尾翼偏轉角由0°增加到30°時,升力逐漸增大,升阻比在0.10～0.34范圍內先增大后減小.

4) 對于GOCE衛星,本文分析得到的阻力特性結果與文獻[17]中多種大氣模型結果吻合良好,驗證了DSMC方法在計算超低軌道衛星阻力特性時的準確性.