李志輝，吳俊林，彭傲平，唐歌實（.中國空氣動力研究與發(fā)展中心超高速所，綿陽62000；2.國家計算流體力學(xué)實驗室，北京009；.北京航天飛行控制中心航天飛行動力學(xué)技術(shù)重點實驗室，北京00094）

·工程技術(shù)·

天宮飛行器低軌控空氣動力特性一體化建模與計算研究

李志輝1，2，吳俊林1，彭傲平1，2，唐歌實3
（1.中國空氣動力研究與發(fā)展中心超高速所，綿陽621000；2.國家計算流體力學(xué)實驗室，北京100191；3.北京航天飛行控制中心航天飛行動力學(xué)技術(shù)重點實驗室，北京100094）

對非規(guī)則板艙組合體天宮飛行器300～200 km低軌道飛行過程空氣動力特性一體化計算建模，提出考慮復(fù)雜構(gòu)型物面遮蓋效應(yīng)面元解析法與經(jīng)修正的Boettcher/Legge非對稱橋函數(shù)，發(fā)展基于三角形面元逼近復(fù)雜外形通用處理方法，建立適于天宮飛行器復(fù)雜物形處理與面元氣動力系數(shù)計算規(guī)則；將DSMC方法與求解Boltzmann模型方程氣體運(yùn)動論統(tǒng)一算法應(yīng)用于天宮飛行器簡化外形，進(jìn)行氣動力當(dāng)?shù)鼗P(guān)聯(lián)參數(shù)計算修正，建立針對大型復(fù)雜結(jié)構(gòu)天宮飛行器低軌道飛行控制過程空氣動力特性一體化快速算法與程序軟件。對大尺度圓柱體外形與天宮飛行器300～200 km不同高度變軌飛行過程不同迎角與側(cè)滑角及帆板平面與本體主軸不同夾角復(fù)雜構(gòu)型氣動力特性計算分析驗證，表明天宮飛行器在200 km以上低軌道飛行控制過程中所受空氣動力系數(shù)隨飛行高度發(fā)生顯著變化（8%～50%），證實長期在軌運(yùn)行的大型航天器若采用統(tǒng)一固定的氣動力系數(shù)，誤差累積巨大，需要采取防護(hù)措施，低軌道飛控大氣阻力仍是制約航天器定軌預(yù)報精度最關(guān)鍵因素。

天宮飛行器；低軌控空氣動力特性；當(dāng)?shù)鼗瘶蚝瘮?shù)；關(guān)聯(lián)參數(shù)確定；工程計算；DSMC方法；統(tǒng)一算法

1　引言

航天器從離地面數(shù)百公里外層空間往地球低軌道再入飛行過程，是一個自數(shù)百到數(shù)十量級Knudsen（克努森）數(shù)高稀薄自由分子流動狀態(tài)多物理場復(fù)雜構(gòu)型極高超聲速流動問題［1］。低地球軌道環(huán)境，是對地觀測衛(wèi)星、氣象衛(wèi)星、通訊遙感衛(wèi)星、載人飛船、空間站等大型復(fù)雜結(jié)構(gòu)航天器運(yùn)行區(qū)域。距離地面200～1000 km的低地球軌道環(huán)境大氣壓力很低、氣體非常稀薄，分子間斷粒子效應(yīng)與稀疏氣體分子分布變化通過對航天器產(chǎn)生力/熱沖擊作用，直接影響航天器運(yùn)行軌道、姿態(tài)和壽命。為了研究和控制航天器低地球軌道環(huán)境飛行過程，在過去計算機(jī)條件不夠先進(jìn)情況下，傳統(tǒng)做法是近似處理，不同流動區(qū)域發(fā)展相應(yīng)的計算方法。對于Knudsen數(shù)很大的稀薄氣體自由分子流區(qū)，氣體流動由無碰撞或近于無碰撞的自由分子流理論控制。為了進(jìn)行航天器從340～200 km低軌道飛行試驗，通常認(rèn)為飛行器在200 km以上，氣體分子平均自由程達(dá)兩百多米以上，其Knudsen數(shù)數(shù)十以上，飛行繞流已屬于完全自由分子流，阻力、升力等空氣動力不再隨飛行高度發(fā)生變化，而采用固定的氣動力系數(shù)參與飛行力學(xué)軌道計算［2］。于是大型航天器開展低軌道飛行試驗采用固定阻力系數(shù)2.2，發(fā)現(xiàn)飛行偏差有時越來越大，以致使用其上RCS姿控系統(tǒng)強(qiáng)制控回標(biāo)稱軌道，造成燃料消耗，甚至有時難以準(zhǔn)確有效控制［3，4］。為此，提出這樣的問題：我國天宮飛行器低軌道飛行所受大氣阻尼系數(shù)是否會變化？在軌天宮飛行器本體直徑數(shù)米、長度達(dá)十余米、其攜帶的太陽電池翼帆板展開寬度數(shù)十米，對如此復(fù)雜巨大不規(guī)則板艙組合體飛行器實施低軌道飛行控制、軌道參數(shù)計算時，能否采用統(tǒng)一、固定的空氣動力系數(shù)？如何開展天宮飛行器低軌道飛行過程空氣動力特性一體化建模與計算研究［4］，剖析天宮飛行器氣動特性隨飛行高度變化規(guī)律？為此，本文在研究建立航天空氣動力學(xué)統(tǒng)一計算理論與各流動區(qū)域多物理場復(fù)雜流動機(jī)理模擬方法研究基礎(chǔ)上，嘗試開展上述問題計算研究。

連接微觀分子動力學(xué)與宏觀流體力學(xué)的介觀Boltzmann（玻爾茲曼）速度分布函數(shù)方程［5］本身可描述各個流域氣體分子輸運(yùn)現(xiàn)象，該方程自1872年提出以來，一直是國際學(xué)術(shù)界追蹤研究的問題，它作為一個高度復(fù)雜高維積分、微分、多相空間多尺度非線性剛性問題［5，6］，精確求解描述各流域氣體流動特征的Boltzmann方程至今未成現(xiàn)實。為此，近五十年來世界航天發(fā)展了基于微觀分子動力學(xué)將分子運(yùn)動與碰撞解耦進(jìn)行隨機(jī)統(tǒng)計模擬的DSMC（直接模擬Monte Carlo）方法，該方法自1963年Bird［7］將其發(fā)展用于稀薄氣體流動模擬至今，已在稀薄氣體動力學(xué)領(lǐng)域獲得了廣泛的應(yīng)用與檢驗，國內(nèi)外學(xué)者將DSMC方法用于模擬包括內(nèi)能、化學(xué)反應(yīng)、電離、輻射等多種作用物理模型的稀薄流區(qū)氣體流動問題研究，已取得較好的應(yīng)用發(fā)展［8-17］，并利用巨型計算機(jī)研究并行算法，引入中樞網(wǎng)格系統(tǒng)開展并行計算［18，19］。在DSMC方法對較高Knudsen數(shù)稀薄氣體流動仿真取得巨大成功同時，該方法因受網(wǎng)格劃分、時間步長和模擬分子數(shù)等方法本身模擬準(zhǔn)則所限，難以對低Knudsen數(shù)復(fù)雜大尺度飛行器氣動問題數(shù)值仿真［14，20，21］。

為了探索跨流域氣體流動問題一體化模擬方法，通過跟蹤國際上關(guān)于Boltzmann方程研究現(xiàn)狀與發(fā)展趨勢［22-27］，結(jié)合從事DSMC方法與計算流體力學(xué)有限差分法研究基礎(chǔ)，過去十余年本文第一作者等從研究Boltzmann方程碰撞松弛間隔理論出發(fā)，確立描述各流域微觀分子輸運(yùn)現(xiàn)象統(tǒng)一的Boltzmann模型方程，提出并應(yīng)用離散速度坐標(biāo)法對分子速度分布函數(shù)數(shù)值離散，發(fā)展可用于速度空間宏觀流動取矩的離散速度數(shù)值積分法。將計算流體力學(xué)有限差分方法推廣拓展到基于時間、位置空間與速度空間的Boltzmann模型方程數(shù)值求解，先后建立起模擬各流域一維、二維、三維繞流問題氣體運(yùn)動論統(tǒng)一算法（GKUA）理論與系列計算技術(shù)［28-33］，并開展了飛船返回艙、近空間飛行器、箔條云等跨流域氣動問題研究［33-37］，解決了系列工程研制需求問題。由于上述求解Boltzmann模型方程統(tǒng)一算法將氣體分子速度分布函數(shù)基于位置空間與速度空間的信息均保存起來數(shù)值求解，需要在位置空間與速度空間組成的六維多相空間進(jìn)行計算，由此占用大量計算機(jī)內(nèi)存資源。如果進(jìn)行三維高超聲速氣體流動問題數(shù)值計算，如天宮飛行器跨流域繞流問題研究，必須開展大規(guī)模并行計算。

另一方面，航天器軌道計算可靠性直接依賴于軌道動力學(xué)方程空氣動力項求解的準(zhǔn)確性，需要與飛行力學(xué)數(shù)值推進(jìn)方法耦合，計算確定航天器任意時刻空間位置與速度，需要發(fā)展跨流域空氣動力特性快速計算方法。為了能在有限時間，建立天宮飛行器低軌道飛行過程空氣動力一體化計算平臺，擬以稀薄氣體高超聲速繞流當(dāng)?shù)鼗瘶蚝瘮?shù)理論為基礎(chǔ)，開展天宮飛行器低軌道飛行控制過程300～200 km空氣動力特性當(dāng)?shù)鼗瘶蚝瘮?shù)工程計算建模，研究使用DSMC方法、求解Boltzmann模型方程的氣體運(yùn)動論統(tǒng)一算法對天宮飛行器簡化外形典型繞流狀態(tài)計算，修正適于大尺度復(fù)雜結(jié)構(gòu)航天器當(dāng)?shù)鼗P(guān)聯(lián)參數(shù)，建立適于天宮飛行器低地球軌道環(huán)境不同高度、馬赫數(shù)與不同迎角、側(cè)滑角空氣動力特性快速計算分析應(yīng)用研究框架。

2　適于天宮飛行器低軌控不同高度空氣動力特性當(dāng)?shù)鼗焖偎惴?/h2>

稀薄氣體流動領(lǐng)域氣動力系數(shù)計算所用的當(dāng)?shù)鼗こ逃嬎惴椒ǎ腔诎虢?jīng)驗橋函數(shù)理論的當(dāng)?shù)貥蚬椒ǎ?8，39］。自由分子流和牛頓連續(xù)流兩個極限流態(tài)都可以用當(dāng)?shù)貧鈩恿Φ慕馕鰜肀硎?，?dāng)?shù)乇砻婷嬖系臍鈩恿ο禂?shù)只依賴于來流和當(dāng)?shù)匦再|(zhì)，如當(dāng)?shù)赜?、表面作用等。?dāng)?shù)剡B續(xù)流和自由分子流系數(shù)之間的間隙可以用半經(jīng)驗的橋公式光滑搭接起來。這是一種加權(quán)函數(shù)，依賴于獨立的關(guān)聯(lián)或相似參數(shù)（如克努森數(shù)），通常由實驗結(jié)果與數(shù)值計算分析確定。對連續(xù)流區(qū)域的壓力系數(shù)計算，采用修正的牛頓無粘流理論以及背風(fēng)真空效應(yīng)法進(jìn)行估算；對自由分子流區(qū)，采用基于不同模型材料修正的Nocilla壁面反射模型進(jìn)行壓力系數(shù)和摩擦力系數(shù)計算；對過渡流區(qū)域，基于修正Boettcher/Legge非對稱橋函數(shù)理論［40］，使用可分段描述的非對稱壓力與摩擦力系數(shù)關(guān)聯(lián)橋函數(shù)。通常針對軸對稱鈍體飛行器設(shè)計，能預(yù)測估算鈍柱形高超聲速飛行器的稀薄氣動力特性，但本文研究大型復(fù)雜結(jié)構(gòu)航天器如天宮飛行器，其空氣動力特性計算不能沿用常規(guī)旋成體當(dāng)?shù)鼗嬎惴椒ǎ?1］?？裳芯渴褂肈SMC方法、求解Boltzmann模型方程從稀薄流到連續(xù)流統(tǒng)一算法［28-37］，基于準(zhǔn)確計算進(jìn)行當(dāng)?shù)鼗瘶蜿P(guān)聯(lián)系數(shù)的選擇調(diào)試與檢驗修正，確定各關(guān)聯(lián)參數(shù)的最佳組合［35，42］，發(fā)展針對縱橫數(shù)十米大尺度復(fù)雜結(jié)構(gòu)航天器，跨越高稀薄自由分子流、過渡流到大氣層內(nèi)低高度各流區(qū)、高低不同馬赫數(shù)、不同迎角與側(cè)滑角氣動力系數(shù)快速計算方法。

2.1高超聲速空氣動力特性當(dāng)?shù)鼗嬎惴椒?/p>

當(dāng)?shù)鼗瘶蚬椒ㄊ紫瓤紤]的是在兩個邊界流域連續(xù)流和自由分子流區(qū)域，氣動力特性可用兩個不同的數(shù)學(xué)函數(shù)稱為連續(xù)流函數(shù)FCont和自由分子流函數(shù)FFM來描述，而在過渡流區(qū)域需要給出一個新的函數(shù)Fb，亦即橋函數(shù)［4，35，38，39］，使其在連續(xù)流或自由分子流分別趨近FCont和FFM。對給定入射角θ，當(dāng)?shù)孛嬖辖?jīng)歸一化壓力和摩擦力系數(shù)表達(dá)為式（1）：

式中，F(xiàn)b，p和Fb，τ分別是壓力和摩擦力橋函數(shù)，依賴于獨立參數(shù)Kn、TW/T0和θ。

對連續(xù)流區(qū)壓力系數(shù)計算，可采用諸如修正牛頓無粘流理論及背風(fēng)真空效應(yīng)法計算；對自由分子流區(qū)，可使用基于不同模型材料修正的Nocilla壁面反射模型進(jìn)行壓力與摩擦力系數(shù)計算；對過渡流區(qū)，研究修正的Boettcher/Legge非對稱橋函數(shù)理論［39］，發(fā)展可分段描述的非對稱壓力與摩擦力系數(shù)關(guān)聯(lián)橋函數(shù)，其中，非對稱壓力橋函數(shù)［38，39，42］可表示為式（2）：

并且Knm，τ、ΔKnτ1、ΔKnτ2為可調(diào)參數(shù)［40］。

當(dāng)?shù)鼗椒P(guān)于飛行器物形處理［38，43，44］有兩種方法，一種是用解析式分別描述飛行器的各部分，這種方法適用于外形簡單且易于解析表達(dá)的軸對稱旋成體飛行器［38，41，44］，但大多數(shù)飛行器外形復(fù)雜且很不規(guī)則，難于進(jìn)行解析表示。于是當(dāng)?shù)鼗嬎惴椒ㄋ蕾嚨奈镄翁幚沓３Ｊ褂昧硪环N較為通用而近似的面元處理方法［4，39，43］，該方法將飛行器物面劃分為若干個面元，利用四邊形或三角形來逼近復(fù)雜外形，計算出每個面元上中心點坐標(biāo)、中心點處的法向及面元面積。

飛行器外形曲面分成若干塊小的曲面，對于每一個小曲面，選用一個小的平面四邊形來代替。這樣，計算小曲面上的氣動力就轉(zhuǎn)換成計算面元上的氣動力，將這些面元上的氣動力加起來就可以得到整個飛行器的氣動力。當(dāng)面元分得很細(xì)小時，由此引起的誤差就更?。?5，44］。

2.2天宮飛行器計算模型構(gòu)建及表面非結(jié)構(gòu)網(wǎng)格物形表征技術(shù)

由于天宮飛行器與其左右寬度達(dá)數(shù)十米的太陽電池翼帆板及連接裝置構(gòu)成一個相當(dāng)復(fù)雜而極不規(guī)則巨大板艙組合體。為了解決天宮飛行器及附件復(fù)雜物形表征的困難，研究引入分區(qū)網(wǎng)格與非結(jié)構(gòu)網(wǎng)格生成技術(shù)，以此獲得能準(zhǔn)確表征復(fù)雜物面邊界的非結(jié)構(gòu)網(wǎng)格，以便進(jìn)行復(fù)雜物形表征處理。模型表面采用非結(jié)構(gòu)貼體網(wǎng)格，通過網(wǎng)格處理得到計算需要的表面網(wǎng)格，具體過程見文獻(xiàn)［4］，其中基于三角非結(jié)構(gòu)網(wǎng)格表征天宮飛行器物形，圖1繪出天宮飛行器計算模型表面三角非結(jié)構(gòu)網(wǎng)格布局。

圖1　天宮飛行器計算模型表面非結(jié)構(gòu)網(wǎng)格布局Fig.1 Com putational m odel and sur face grid layout of the TG spacecraft

2.3天宮飛行器跨流域氣動特性快速算法

根據(jù)稀薄過渡流流動特征，當(dāng)?shù)鼗焖偎惴▽⑾”∵^渡流區(qū)壓力橋函數(shù)（2）式與摩擦力橋函數(shù)（3）式用非對稱誤差函數(shù)形式關(guān)聯(lián)，兩式中Knm，p、ΔKnp1、ΔKnp2及Knm，τ、ΔKnτ1、ΔKnτ2應(yīng)根據(jù)風(fēng)洞試驗數(shù)據(jù)或數(shù)值計算結(jié)果進(jìn)行調(diào)試確定，由于天宮飛行器的大尺度復(fù)雜結(jié)構(gòu)與在軌飛行高稀薄自由分子流區(qū)繞流特點，風(fēng)洞實驗條件并不現(xiàn)實，可依靠高性能計算機(jī)開展DSMC方法［4，13，17］或求解Boltzmann模型方程統(tǒng)一算法［30-34］典型飛行狀態(tài)數(shù)值計算，調(diào)試選擇與計算修正上述六個關(guān)聯(lián)參數(shù)，使當(dāng)?shù)鼗こ逃嬎憬Y(jié)果最大程度地滿足實際計算精度。由于本文旨在有限時間搭建天宮飛行器低軌道飛行過程300～200 km空氣動力特性初步計算平臺，作為算法驗證，橋函數(shù)中當(dāng)?shù)鼗烧{(diào)關(guān)聯(lián)參數(shù)的確定采用250 km高度與天宮飛行器具有同樣尺寸的圓柱體繞流數(shù)值計算結(jié)果進(jìn)行修正確定。首先利用數(shù)值計算得到該圓柱體繞流典型狀態(tài)所受到的氣動阻力系數(shù)為參考對象；改變上述六個關(guān)聯(lián)參數(shù)，使用當(dāng)?shù)鼗焖偎惴ㄟM(jìn)行調(diào)試計算修正；最終獲得一套針對天宮飛行器計算適用的六個關(guān)聯(lián)參數(shù)的最佳組合，使得經(jīng)關(guān)聯(lián)參數(shù)修正后的氣動特性當(dāng)?shù)鼗嬎惴椒軡M足本文關(guān)于天宮飛行器在不同高度、不同馬赫數(shù)、不同迎角與側(cè)滑角的氣動力系數(shù)計算要求。圖2給出了250 km高度0°迎角上述大尺度圓柱體繞流狀態(tài)壓力等值線云圖分布，可看出圓柱體前端表面壓力值遠(yuǎn)大于后表面，氣流完全附著物面強(qiáng)擾動，流場呈現(xiàn)出高稀薄自由分子流動特征。

圖2　 H=250 km，α=0°大尺度圓柱體繞流壓力等值線分布Fig.2 Pressure distribution contours around largescale cylinder body w ith H=250 km，α=0°

數(shù)值計算得到0°迎角時圓柱體模型阻力系數(shù)為3.169025，以此為基礎(chǔ)得到當(dāng)?shù)鼗焖偎惴ㄖ袠蚝瘮?shù)的可調(diào)關(guān)聯(lián)參數(shù)如式（5）～（6）：

以上述六個當(dāng)?shù)鼗P(guān)聯(lián)參數(shù)為基礎(chǔ)，發(fā)展可用于確定大尺度飛行器在不同高度、馬赫數(shù)、迎角與側(cè)滑角的氣動力系數(shù)快速算法。由此計算得到上述圓柱體繞流阻力系數(shù)3.169013，與數(shù)值模擬結(jié)果相當(dāng)一致，彼此偏差僅0.000379%。表1列出了該狀態(tài)數(shù)值計算與當(dāng)?shù)鼗焖偎惴ǚ謩e得到4°迎角上述圓柱體模型氣動力系數(shù)比較情況，可看出二者吻合很好，偏差范圍0.063%～2.85%，驗證了本文提出發(fā)展低軌道飛行過程大尺度飛行器300 km到200 km空氣動力特性快速計算平臺可行性。

表1　當(dāng)?shù)鼗焖偎惴ㄓ嬎鉎=250 km，α=4°圓柱體繞流氣動力系數(shù)驗證比較Table 1 Validation and comparison of aerodynam ic coefficients around the large-scale cylinder w ith H=250 km，α=4°

3　天宮飛行器簡化外形200 km至300 km空氣動力特性計算分析

3.1算法驗證

為了進(jìn)一步驗證本文發(fā)展的當(dāng)?shù)鼗焖偎惴Ｐ驼_性以及200～300 km高度飛行的氣動力系數(shù)變化規(guī)律，分別采用數(shù)值模擬方法和當(dāng)?shù)鼗焖偎惴ㄓ嬎懔?00 km和300 km上述大尺度圓柱體0°迎角氣動力系數(shù)，表2給出了兩類方法計算0°迎角200 km和300 km阻力系數(shù)Cd比較情況。可看出，兩種計算方法得到的氣動力系數(shù)完全吻合，最大偏差不超過0.48%，證實本文發(fā)展的當(dāng)?shù)鼗焖偎惴ㄕ_可靠；在300 km飛行高度阻力系數(shù)比200 km阻力系數(shù)均高出至少7.43%，說明近地軌道飛行過程200～300 km不同高度飛行器阻力系數(shù)是隨飛行高度增加而增大，這反映出傳統(tǒng)工程處理手段認(rèn)為200～300 km氣體分子平均自由程λ∞達(dá)235.2～2595.37m，遠(yuǎn)大于飛行器特征尺寸，來流氣體處于自由分子流態(tài)，飛行器繞流使用固定不變阻力系數(shù)統(tǒng)一處理方式對類天宮飛行器兩艙結(jié)構(gòu)大尺度圓柱體繞流會產(chǎn)生至少7.43%偏差。

3.2天宮飛行器200～300 km氣動力特性計算分析

為使本文當(dāng)?shù)鼗焖偎惴ㄓ嬎?00～200 km期間天宮飛行器氣動特性結(jié)果能盡可能真實可靠，本文對任一電池翼帆板平面與天宮飛行器本體主軸夾角0°、22.5°、45°、67.5°、90°確定的計算模型，進(jìn)行300～200 km不同高度300 km、280 km、250 km、220 km、200 km，不同迎角α= -13.5°，0°，13.5°，22°，38°與側(cè)滑角β=5.5°，2.6°，0°，-2.6°，-4.1°的天宮飛行器氣動特性計算，提供了大量可供軌道控制彈道飛行力學(xué)所用空氣動力系數(shù)計算數(shù)據(jù)。下面僅列兩組典型狀態(tài)。

表2　 200 km和300 km飛行高度0°迎角圓柱體繞流阻力系數(shù)驗證比較Table 2 Com parison of drag coefficients around the large-scale cylinder w ith H=200、300 km，α=0°

3.2.1帆板平面與天宮飛行器主軸夾角為0°氣動特性計算

擬定天宮飛行器250 km、高度13.5°迎角、0°側(cè)滑角外形，根據(jù)數(shù)值計算調(diào)試確定的當(dāng)?shù)鼗瘶蚝瘮?shù)關(guān)聯(lián)參數(shù)為Knm，p=0.9、ΔKnp1=3.5、ΔKnp2=3.5，Knm，τ=0.85、ΔKnτ1=3.5、ΔKnτ2=3.0。圖3繪出天宮飛行器簡化外形展向水平面與縱向豎直平面內(nèi)壓力等值線分布，可看出繞流流場附著物面流動，僅在物面附近出現(xiàn)強(qiáng)擾動，逐漸擴(kuò)展到遠(yuǎn)前方，表現(xiàn)出高稀薄自由分子流動特征，等值線分布范圍相當(dāng)廣，達(dá)天宮飛行器本體尺寸2～3倍。

圖3　天宮飛行器簡化外形帆板與主軸夾角0°高度250 km迎角13.5°、側(cè)滑角0°軸對稱面壓力等值線分布Fig.3 Pressure distribution contours around the sim p lified TG-spacecraft w ithθ=0°，H=250 km，α= 13.5°，β=0°

表3列出了帆板平面與本體主軸夾角θ=0°時，考慮帆板與天宮飛行器兩艙本體三腳架連接（見圖1）天宮一號目標(biāo)飛行器計算模型，采用本文算法計算得到α=-13.5°、β=5.2°條件下，不同飛行高度氣動力系數(shù)，包括軸向力Ca、法向力Cy、側(cè)向力Cz、阻力Cd、升力Cl、側(cè)力Czl，其中“最大差異”指300 km的氣動力系數(shù)與200 km相應(yīng)氣動力系數(shù)值之間的相對偏差?？煽闯?，軸向力系數(shù)與法向力系數(shù)最大差異在5%～10%，而側(cè)向力系數(shù)最大差異在7%～13%。但是，由軸向力、法向力和側(cè)向力系數(shù)經(jīng)坐標(biāo)變換得到的風(fēng)軸系下阻力、升力和側(cè)力系數(shù)的最大差異分別達(dá)到了10%、30%和80%。由此可見，在200～300 km的大型復(fù)雜結(jié)構(gòu)航天器飛行過程，若采用固定統(tǒng)一的氣動力系數(shù)進(jìn)行飛行控制會產(chǎn)生較大誤差。

表3　天宮飛行器帆板與本體0°夾角、-13.5°迎角、5.2°側(cè)滑角200～300 km各高度氣動力系數(shù)Table 3 Aerodynam ic coefficients around the TG-spacecraft from 200～300 km withθ=0°，α=-13.5°，β=5.2°

3.2.2帆板平面與天宮飛行器主軸夾角為45°氣動特性計算

對θ=45°，天宮飛行器與帆板形成計算模型見圖4（a），根據(jù)數(shù)值計算確定Knm，p=0.9、ΔKnp1=3.2、ΔKnp2=3.5，Knm，τ=0.65，ΔKnτ1=3.1，ΔKnτ2=3.0。圖4（b）繪出天宮飛行器簡化外形繞流展向水平面與縱向豎直平面內(nèi)壓力等值線分布云圖，同樣看出繞流流場附著物面流動，僅在物面附近出現(xiàn)強(qiáng)擾動，并逐漸擴(kuò)展到遠(yuǎn)前方，表現(xiàn)出高稀薄自由分子流動特征。

表4介紹了本文方法計算天宮飛行器θ= 45°、α=0°、β=0°條件下，200～300 km間隔20 km各飛行高度氣動力系數(shù)。可看出，300 km與200 km的氣動力系數(shù)差別很明顯，隨飛行迎角狀態(tài)不同，這種差異也不同。對于某些飛行狀態(tài)，出現(xiàn)側(cè)力系數(shù)正負(fù)符號變化，法向力系數(shù)“最大差異”達(dá)25%以上。

圖4　天宮飛行器簡化外形帆板與主軸夾角θ= 45°，H=250 km，α=0°，β=0°模型與對稱面壓力分布Fig.4 Pressu re distribution contours aroundthe simp lified TG-spacecraft withθ= 45°，H=250 km，α=0°，β=0°

表4　θ=45°，α=0°，β=0°飛行過程200～300 km各高度氣動力系數(shù)Table 4 Aerodynam ic coefficients around the TG-spacecraft from 200～300 km w ithθ=45°，α=0°，β =0°

4　結(jié)論

針對在軌服役天宮飛行器板艙組合體大型復(fù)雜構(gòu)型特點，研究修正的Boettcher/Legge非對稱橋函數(shù)，發(fā)展了可靠模擬不同高度、不同迎角與側(cè)滑角及不同復(fù)雜外形高超聲速空氣動力特性快速計算方法；利用三角形逼近復(fù)雜物形，研制適于天宮飛行器復(fù)雜物面表征的非結(jié)構(gòu)網(wǎng)格生成技術(shù)；研究通用面元處理法，發(fā)展適于天宮飛行器復(fù)雜物形處理與面元氣動力系數(shù)計算規(guī)則；建立了適于天宮飛行器近地軌控200～300 km空氣動力特性一體化快速算法應(yīng)用研究平臺，揭示大尺度天宮飛行器空氣動力變化規(guī)律，并得到如下結(jié)論：

1）天宮飛行器在300 km到200 km變軌飛行過程中，氣動力系數(shù)隨飛行高度變化幅度較大：體軸系中軸向力系數(shù)與風(fēng)軸系中阻力系數(shù)變化幅度5%～9%；法向力、側(cè)向力等其它氣動力系數(shù)變化幅度更大，有的甚至超過50%，某些狀態(tài)（如帆板平面與本體主軸夾角67.5°，以0°迎角飛行）法向力系數(shù)“最大差異”達(dá)48%，而一些飛行狀態(tài)升力和側(cè)力系數(shù)出現(xiàn)了正負(fù)符號與方向差別。對長期在軌飛行大型航天器，這種誤差累積后果難以控制，須采取防護(hù)措施；

2）對天宮飛行器這類由太陽電池翼與主體構(gòu)成大型復(fù)雜板艙結(jié)構(gòu)航天器低軌道飛行300～200 km氣動力系數(shù)隨飛行高度、姿態(tài)與太陽電池翼繞連接軸轉(zhuǎn)動外形變化而改變，變化幅度與飛行狀態(tài)相關(guān)；不同迎角、側(cè)滑角氣動力系數(shù)存在較大差異；同時，太陽能電池翼圍繞主體連接軸旋轉(zhuǎn)對天宮整器氣動力特性也產(chǎn)生較大影響；

3）天宮飛行器近地軌道飛行過程，氣動力系數(shù)不僅發(fā)生變化，且變化幅度與方向不同。如軸向力和阻力系數(shù)隨飛行高度300～200 km變化的“最大差異”在5%～10%，該類氣動力系數(shù)本身基數(shù)大，隨高度變化量絕對值較大（一般在0.18左右）；與之相應(yīng)法向力系數(shù)“最大差異”有時達(dá)48%，但其隨高度變化量絕對值為0.06；而某些狀態(tài)，如太陽電池翼平面與本體主軸夾角45°，以0°迎角、2.6°側(cè)滑角所受側(cè)力系數(shù)200 km為-0.0001，而300 km變?yōu)?.0032，出現(xiàn)正負(fù)側(cè)力方向的改變，它們之間“最大差異”就會出現(xiàn)異常嚴(yán)重，而該側(cè)力系數(shù)隨高度變化量絕對值僅0.0033。這種由太陽電池翼旋轉(zhuǎn)所致氣動力方向上的變化會直接危及航天器精細(xì)軌控。

（致謝：本文寫作過程中和李革非研究員、李勰工程師及課題組蔣新宇助理工程師等進(jìn)行了交流討論，部分計算在總參三部北方計算中心完成，特此感謝?。?/p>

［1］PilinskiM D，Argrow BM，Palo SE.Drag coefficientsofsatellites with concave geometries：comparingmodels and observations［J］.Journal of Spacecraft and Rockets，2011，48（2）：312-325.

［2］Moe K，Bowman B R.The effects of surface composition and treatment on drag coefficients of spherical satellites［C］// AAS/AIAA Astrodynamics Specialist Conference，Americal Astronautical Soc.，Paper 05-258，Springfield，VA，2005.

［3］Moe K，Moe M M，Wallace SD.Improved satellite drag coefficient calculations from orbitalmeasurements of energy accommodation［J］.Journal of spacecraft and rockets，1998，35（3）：266-272.

［4］李志輝，吳俊林，彭傲平.天宮飛行器低軌道飛行控制過程空氣動力建模與計算研究［R］.中國國防科技報告，綿陽：中國空氣動力研究與發(fā)展中心超高速研究所，2013. Li Z H，Wu JL，Peng A P.Aerodynamicsmodeling and calculation during low-orbit flying control of TG-1 target spacecraft［R］.Technical Report，Hypersonic Aerodynamics Institute，China Aerodynamic Research and Development Centre，2013.（in Chinese）

［5］Cercignani C.The Boltzmann Equation and its Applications［M］.Springer Verlag，Berlin，1988：55-58.

［6］應(yīng)純同.氣體輸運(yùn)理論及應(yīng)用［M］.北京：清華大學(xué)出版社，1990：123-136. Ying C T.Gas transport theory and application［M］.Tsinghua Univ.Press，Beijing，1990：123-136.（in Chinese）

［7］Bird G A.Approach to translational equilibrium in a rigid sphere gas［J］.Physics of Fluids（1958-1988），1963，6（10）：1518-1519.

［8］Pham-Van-Diep G，Erwin D，Muntz E P.Nonequilibrium molecular motion in a hypersonic shock wave［J］.Science，1989，245（4918）：624-626.

［9］Carlson A B，Hassan H A.Direct simulation of reentry flows with ionization［R］.AIAA 90-144，1990.

［10］Haas B L，Boyd ID.Models for direct Monte Carlo simulation of coupled vibration-dissociation［J］.Physics of Fluids A：Fluid Dynamics，1993，5（2）：478-489.

［11］Bird G A.Molecular Gas Dynamics and Direct Simulation of Gas Flows［M］.Oxford University Press，1994：78-189

［12］樊菁，沈青.過渡領(lǐng)域高超聲速圓柱繞流的直接模擬［J］.空氣動力學(xué)學(xué)報，1995，13（4）：405-413. Fan J，Shen Q.Direct simulation of hypersonic transitional flow around circular cylinder，Acta Aerodynamica Sinica，1995，13（4）：405-413.（in Chinese）

［13］李志輝，吳振宇.阿波羅指令艙稀薄氣體動力學(xué)特征的蒙特卡羅數(shù)值模擬［J］.空氣動力學(xué)學(xué)報，1996，14（2）：230-233. Li Z H，Wu Z Y.DSMC simulation of hypersonic rarefied flow past Apollo-CM，Acta Aerodynamica Sinica，1996，14（2）：230-233.（in Chinese）

［14］Ivanov M S，Gimelshein S F.Computational hypersonic rarefied flows［J］.Annual Review of Fluid Mechanics，1998，30（1）：469-505.

［15］吳其芬，陳偉芳.高溫稀薄氣體熱化學(xué)非平衡流動的DSMC方法［M］.長沙：國防科技大學(xué)出版社，1999：165-176. Wu Q F.Chen W.F.DSMC method for high-temperature thermochemical nonequilibrium flow of rarefied gas［M］.National University of Defense Technology Press，ChangSha，1999：165-176.（in Chinese）

［16］Sun Q，F(xiàn)an J，Boyd ID.Improved sampling techniques for the direct simulation Monte Carlomethod［J］.Computers& Fluids，2009，38（2）：475-479.

［17］梁杰，閻超，李志輝，等.稀薄過渡流區(qū)橫向噴流干擾效應(yīng)數(shù)值模擬研究［J］.空氣動力學(xué)學(xué)報，2013，31（1）：27-33. Liang J，Yan C.，Li Z H，et al.Numerical investigation of lateral jet interaction effects in rarefied transition flow regime，Acta Aerodynamica Sinica，2013，31（1）：27-33.（in Chinese）

［18］Wong B C，Long L N.Direct simulation Monte Carlo（DSMC）on the connection machine［R］.AIAA 92-0564，1992.

［19］Scanlon T J，RoohiE，White C，etal.An open source，parallel DSMC code for rarefied gas flows in arbitrary geometries［J］.Computers&Fluids，2010，39（10）：2078-2089.

［20］李志輝，楊彥廣，梁杰，等.飛船返回艙高空稀薄氣動特性研究［C］.全國載人航天工程氣動工作總結(jié)暨關(guān)鍵問題研討會，6月，北京，2004. Li Z H.Yang Y G.，Liang J，etal.Study of rarefied aerodynamics for spacecraft capsule in high-altitude，Symposium of front key problems for national manned spaceflight，Beijing，June，2004.（in Chinese）

［21］Yoon S，Gnoffo P A，White JA，et al.Computational challenges in hypersonic flow simulations［R］.AIAA 2007-4265，2007.

［22］Bhatnagar P L，Gross E P，Krook M.A model for collision processes in gases.I.Small amplitude processes in charged and neutral one-component systems［J］.Physical review，1954，94（3）：511-525.

［23］Holway Jr L H.New statistical models for kinetic theory：methods of construction［J］.Physics of Fluids（1958-1988），1966，9（9）：1658-1673.

［24］Shakhov E M.Generalization of the Krook kinetic relaxation equation［J］.Fluid Dynamics，1968，3（5）：95-96.

［25］Yang JY，Huang JC.Rarefied flow computations using nonlinearmodel Boltzmann equations［J］.Journal of Computational Physics，1995，120（2）：323-339.

［26］Titarev V A，Shakhov E M.Heat transfer and evaporation from a p lane surface into a half-space upon a sudden increase in body temperature［J］.Fluid dynamics，2002，37（1）：126-137.

［27］Mieussens L.Discrete-velocitymodels and numerical schemes for the Boltzmann-BGK equation in plane and axisymmetric geometries［J］.Journal of Computational Physics，2000，162（2）：429-466.

［28］李志輝，張涵信.稀薄流到連續(xù)流的氣體運(yùn)動論統(tǒng)一數(shù)值算法初步研究［J］.空氣動力學(xué)學(xué)報，2000，18（3）：255-263. Li ZH.，Zhang H X.Study on gas-kinetic algorithm for flows from rarefied transition to continuum，Acta Aerodynamica Sinica，2000，18（3）：255-263.（in Chinese）

［29］Li Z，Zhang H.Study on gas kinetic algorithm for flows from rarefied transition to continuum［C］//Proc.of 22ndInternational Symposium on Rarefied Gas Dynamics，Sydney，Australia，Jul.9-14，2000：628-636.

［30］李志輝.從稀薄流到連續(xù)流的氣體運(yùn)動論統(tǒng)一數(shù)值算法研究［D］.綿陽：中國空氣動力研究與發(fā)展中心研究生部，2001. Li Z H.Study on gas kinetic unified algorithm for flows from rarefied transition to continuum，Ph.D thesis，China Aerodynamics Research and Development Center，Jan.，2001.（in Chinese）

［31］Li ZH，Zhang H X.Study on gas kinetic unified algorithm for flows from rarefied transition to continuum［J］.Journal of Computational Physics，2004，193（2）：708-738.

［32］李志輝，張涵信.基于Boltzmann模型方程的氣體運(yùn)動論統(tǒng)一算法研究［J］.力學(xué)進(jìn)展，2005，35（4）：559-576. Li Z H.，Zhang H X.Study on gas kinetic numerical algorithm using Boltzmann model equation.Advances in Mechanics 2005；35（4）：559-576.（in Chinese）

［33］Li Z，Zhang H.Gas-kinetic numericalmethod for solvingmesoscopic velocity distribution function equation［J］.Acta Mechanica Sinica，2007，23（2）：121-132.

［34］Li Z H，Zhang H X.Gas-kinetic numerical studies of threedimensional complex flows on spacecraft re-entry［J］.Journal of Computational Physics，2009，228（4）：1116-1138.

［35］李志輝，梁杰，李四新，等.箔條云跨流域整體氣動特性計算研究［J］.空氣動力學(xué)學(xué)報，2011，29（01）：59-67. Li Z H.，Liang J，Li SX，et al.Computing study on holistic aerodynamics of chaff cloud covering various flow regimes，Acta Aerodynamica Sinica，2011，29（1）：59-67.（in Chinese）

［36］Li Z H.Gas-kinetic unified algorithm for re-entering complex flows covering various flow regimes by solving boltzmann model equation［M］//Hall J.Advances in Spacecraft Technologies.Croatia：InTech，2011：273-332.

［37］李志輝，彭傲平，吳俊林，等.稀薄氣體自由分子流到連續(xù)流跨流域氣動力熱繞流統(tǒng)一算法研究［J］.載人航天，2013，19（2）：81-91. Li ZH.，Peng A P，Wu JL，etal.Gas-kinetic unified algorithm for aerothermodynamics from free molecular regime to continuum，Manned Spaceflight，2013，19（2）：81-91.（in Chinese）

［38］Gentry A E，Smyth D N，OliverW R.The Mark IV supersonic-hypersonic arbitrary-body program，Vol.11：Program Formulation［R］.AFFDL TR-73-159，1973.

［39］Kashkovskyt A V.Computational tools for rarefied aerodynamics［C］//Rarefied Gas Dynamics：Technical Papers from the Proceedings of the Eighteenth International Symposium on Rarefied Gas Dynamics，University of British Columbia，Vancouver，British Columbia，Canada，1992：115-126.

［40］Mosanov SV，F(xiàn)reedlender O G，Nikiforov A P，et al.Experimental determination ofmomentum transfer coefficients in hypersonic freemolecular flow and distribution function recovery of reflected molecules［C］//Proceedings of the XIII International Symposium Rarefied Gas Dynamics，Vol.1，Plenum，New York，1985：669-676.

［41］謝硯儒，李志輝.尖鈍錐誘餌再入阻力特性工程計算研究［R］.綿陽：中國空氣動力研究與發(fā)展中心，中國國防科技報告，1992. Xie Y R.，Li ZH.，Study ofengineering calculation of re-entry drag characteristics for sharp-cone body［R］.Technical Report，Hypersonic Aerodynamics Institute，China Aerodynamic Research and DevelopmentCentre，1992.（in Chinese）

［42］李志輝.箔條及箔條云整體氣動特性建模方法研究［R］.綿陽：中國空氣動力研究與發(fā)展中心，中國國防科技報告，2005. Li Z H.，Modeling study on holistic aerodynamic characteristic of foil and foil-cloud cluster［R］.Technical Report，Hypersonic Aerodynamics Institute，China Aerodynamic Research and Development Centre，2005.（in Chinese）

［43］謝硯儒，劉濤.彈頭稀薄氣動特性的工程計算方法研究［R］.綿陽：中國空氣動力研究與發(fā)展中心，中國國防科技報告，1990.Xie Y R.，Liu T.，Study ofengineering calculation method of rarefied aerodynamics for bullet［R］.Technical Report，Hypersonic Aerodynamics Institute，China Aerodynamic Research and Development Centre，1990.（in Chinese）

［44］粱杰，李志輝.目標(biāo)稀薄氣體特性工程計算［R］.綿陽：中國空氣動力研究與發(fā)展中心，國防科技報告，1999. Li Z H.，Liang JEngeering calculation of target rarefied aerodynamic characteristics［R］.Technical Report，Hypersonic Aerodynamics Institute，China Aerodynamic Research and Development Centre，1999.（in Chinese）

Unified Modeling and Calculation of Aerodynam ics Characteristics during Low-O rbit Flying Control of the TG Vehicle

LIZhihui1，2，WU Junlin1，PENG Aoping1，TANG Geshi3
（1.Hypervelocity Aerodynamics Institute，China Aerodynamics Research and Development Center，Mianyang 621000，China；2.National Laboratory for Computational Fluid Dynamics，Beijing 100191，China；3.Science&Technology on Aerospace Flight Dynamics Laboratory，Beijing Aerospace Control Center，Beijing100094，China）

The surface analyticalmethod，considering surface shielding effect of the complex structure and the modified Boettcher/Legge nonsymmetric bridge correction function，was proposed to computationally model aerodynamic characteristics of the large-scale Tiangong spacecraft of irregular plate-capsule assembly.The complex configuration processing and computing rules of surface element aerodynamic coefficientswere setup by developing a general triangle elementapproximation for complex shapes.A unified fast algorithm and computing software for aerodynamic characteristics of large-scale complex spacecraft structures were developed，by computing correction of local correlation parameters during Tiangong spacecraft’s low-earth orbit flight control，in which the DSMC method and the gas-kinetic unified algorithm solving the Boltzmann model equation were applied.Itis indicated that，by computing the aerodynamics of the Tiangong spacecraft in 300～200 km altitude-orbit flight process，with various flying heights，various angles of attack and various angles between the panel and the principal axis，the aerodynamic coefficients varied remarkably，with a range of 8%～50%，with the altitude change，and the atmospheric drag in low-earth orbit flight control was the key factor of the spacecraft orbit prediction accuracy.It is validated that，if the fixed and constant unity aerodynamic coefficient is used for long-term orbit flight of a large spacecraft，the error accumulation will be huge，and protectivemeasures are needed.

TG-1 target spacecraft；aerodynamic characteristic during low-orbit control；local bridge function；correlation parameter evaluation；engineering calculation；DSMCmethod；Gas-Kinetic U-nified Algorithm（GKUA）

V411.4

A

1674-5825（2015）02-0106-09

2014-09-12；

2015-02-26

國家973計劃（2014CB744100）；國家自然科學(xué)基金（11325212、91016027）；國防基礎(chǔ)科研基金（51313030104）作者簡介：李志輝（1968-），男，博士，研究員，研究方向為跨流域空氣動力學(xué)。E-mail：zhli0097@x263.net