劉高同,孫宇,張磊
北京衛(wèi)星環(huán)境工程研究所,北京 100094
?
火星大氣環(huán)境模擬裝置設(shè)計及仿真分析研究
劉高同*,孫宇,張磊
北京衛(wèi)星環(huán)境工程研究所,北京 100094
對火星表面大氣環(huán)境特性進(jìn)行了研究,通過選取合適的計算方法并結(jié)合FLUENT流體有限元計算軟件對火星表面稀薄氣體內(nèi)部環(huán)流進(jìn)行了模擬仿真分析,提出了以動量源模擬風(fēng)扇段內(nèi)流的仿真方法,并進(jìn)行了可行性討論。進(jìn)一步實(shí)現(xiàn)了針對圓柱形模擬裝置多工況下的內(nèi)部氣體流場穩(wěn)態(tài)和非穩(wěn)態(tài)計算仿真,并對計算結(jié)果進(jìn)行了分析討論,為火星大氣環(huán)境模擬裝置的設(shè)計提供了技術(shù)支持和參考。
火星大氣;稀薄氣體;流體仿真;環(huán)境模擬裝置;計算流體動力學(xué)
火星探測計劃以及火星采樣返回、載人登陸和火星基地建設(shè)等可能的后續(xù)任務(wù),對準(zhǔn)確預(yù)測探測器在火星大氣環(huán)境中的氣動特性提出了嚴(yán)格要求。陌生的大氣環(huán)境、技術(shù)手段的局限和基礎(chǔ)數(shù)據(jù)的缺乏對氣動特性預(yù)測帶來了極大的困難?;鹦潜砻嫦”〈髿獬煞峙c地球大氣顯著不同,密度也只有地球大氣的1%,因此飛行器進(jìn)入火星大氣與再入地球大氣的過程及遭遇的氣動環(huán)境非常不同[1]。目前,針對火星探測器氣動特性的研究主要集中在火星軌道飛行器氣動減速[2]、飛行器進(jìn)入大氣層高速降落和飛行器傘降等方面[3-5],針對火星近地表面探測器低速情況下的氣動特性研究較少,尤其是試驗(yàn)研究國內(nèi)尚無人開展,因此,設(shè)計建造火星大氣環(huán)境模擬裝置,針對火星表面稀薄大氣流態(tài)進(jìn)行模擬,進(jìn)而研究火星表面大氣環(huán)境對探測器周圍流場結(jié)構(gòu)和氣動特性的影響,對于掌握火星探測進(jìn)入階段的飛行器氣動特性及探測器外形設(shè)計,具有重要科學(xué)意義與工程價值。本文結(jié)合火星大氣特性選取適當(dāng)?shù)挠嬎隳P歪槍鹦潜砻娲髿猸h(huán)境模擬裝置進(jìn)行了內(nèi)部環(huán)流仿真計算,對裝置內(nèi)部流場進(jìn)行了計算仿真,討論了不同構(gòu)型裝置在設(shè)計工況下的性能。
由于火星表面大氣遠(yuǎn)比地球大氣稀薄,因此在模擬計算方法選取的過程中需考慮火星大氣的稀薄特性。稀薄氣體動力學(xué)理論中一般采用努森數(shù)來表征某個氣象流的稀薄程度[6],努森數(shù)即氣體分子平均自由程λ與流動特征長度L的比值,其具體表達(dá)式為:
(1)
根據(jù)給定的溫度和壓力范圍計算得到火星表面稀薄氣體分子平均自由程最大值為8.357×10-5m,最小值為2.235×10-6m,取氣體流動宏觀量梯度的標(biāo)尺長度為1 m,努森數(shù)范圍為2.235×10-6~8.357×10-5,根據(jù)錢學(xué)森關(guān)于氣體流動區(qū)域的劃分[7],根據(jù)稀薄流體空氣動力學(xué),火星大氣裝置內(nèi)部流場屬于連續(xù)流場范圍,N-S方程在其所處的氣體條件范圍下是適用的。
2.1火星大氣環(huán)境模擬裝置
模擬裝置容器艙體仿真計算模型如圖1所示。其中,環(huán)境箱外殼為圓柱形結(jié)構(gòu)形式,具體幾何尺寸和艙內(nèi)流體環(huán)境如表1所示。
表1 環(huán)境模擬裝置參數(shù)
圖1 模擬裝置容器艙體仿真計算模型Fig.1 Simulation calculation model of environment simulation device
考慮到容器的對稱性,為了減小模型計算量,建立了環(huán)境模擬裝置的半個模型,并在對稱面設(shè)置對稱邊界條件,環(huán)境模擬裝置的流場計算網(wǎng)格劃分如圖2所示,其中黃色部分為動量源施加區(qū)域,用于模擬涵道風(fēng)扇,紅色區(qū)域?yàn)樵囼?yàn)段,計算過程中對試驗(yàn)段風(fēng)速進(jìn)行監(jiān)控,保證該區(qū)域平均風(fēng)速達(dá)到計算工況的要求。
圖2 圓形裝置流體計算網(wǎng)格劃分Fig.2 Fluid calculation grid of cylindrical environment simulation device
2.2內(nèi)部環(huán)流流體物性參數(shù)設(shè)定
火星表面稀薄大氣約含二氧化碳95.3%、氮?dú)?.7%和氬氣1.6%,還有微量的氧氣、一氧化碳、氖氣、氪氣和氙氣,在模塊流體參數(shù)設(shè)定階段,需根據(jù)這幾種氣體的質(zhì)量分?jǐn)?shù)以及流場壓強(qiáng)和溫度合成火星大氣的物性參數(shù),在現(xiàn)有氣體組分熱物性質(zhì)的基礎(chǔ)上,本文結(jié)合相關(guān)文獻(xiàn)[8],研究了多組分氣體熱物性參數(shù)的計算方法,并取了合適的計算模型,對混合氣體的密度、比熱、粘度、導(dǎo)熱系數(shù)可采用以下公式進(jìn)行計算。
(1)混合氣體密度
混合氣體密度主要根據(jù)各組分氣體的體積分?jǐn)?shù)及其在制定壓力和溫度條件下的密度計算得到。具體計算公式為:
(2)
式中:ρi為各組分氣體在指定氣壓和溫度下的密度;yi為各組分氣體在混合氣體中的濃度。
(2)混合氣體黏度
混合氣體的黏度可以通過各組分的純物質(zhì)黏度、相對分子質(zhì)量及濃度,根據(jù)一定的混合規(guī)則求得。嚴(yán)格的Chapaman-Enskog動力論可以推廣用于計算多元?dú)怏w混合物的黏度。經(jīng)簡化,略去二階的影響,則嚴(yán)格的數(shù)值解可近似地用級數(shù)表示:
(3)
式中:yi為組分氣體濃度;ηi為組分氣體的黏度;φij為組分氣體i和組分氣體j的結(jié)合因子。
(3)混合氣體導(dǎo)熱系數(shù)
文獻(xiàn)中提供了很多混合氣體導(dǎo)熱系數(shù)計算的方法,本文采用以下公式進(jìn)行計算:
(4)
式中:λ混合為混合氣體的導(dǎo)熱系數(shù);λi為組分氣體i的導(dǎo)熱系數(shù)。
(4)混合氣體的定壓比熱容
混合氣體密度主要根據(jù)各組分氣體的濃度及其在制定壓力和溫度條件下的定壓比熱容計算得到。具體計算公式為:
(5)
式中:cpi為各組分氣體在指定氣壓和溫度下的定壓比熱容。
2.3涵道風(fēng)扇模擬
涵道風(fēng)扇是引起模擬裝置內(nèi)氣體環(huán)流的驅(qū)動器,涵道風(fēng)扇模擬的準(zhǔn)確度直接影響到容器艙內(nèi)氣體環(huán)流的精度??紤]到目前涵道風(fēng)扇并未開始詳細(xì)設(shè)計,無法提供可用于CFD分析計算的的風(fēng)扇細(xì)節(jié)幾何模型,分析計算中將風(fēng)扇所在區(qū)域簡化為一個動量源[10],以動量源代替槳葉對流場的作用來模擬涵道風(fēng)扇的氣動特性。在連續(xù)流場中關(guān)于坐標(biāo)軸i的動量方程可形成如下形式:
(6)
式中:ρ為流體密度;ui為流體沿坐標(biāo)軸i方向的流動速度;p為流體壓力;Si表示風(fēng)扇沿坐標(biāo)軸i方向上的單位體積作用力,即為動量源項(xiàng)。
FLUENT軟件內(nèi)部只能針對指定區(qū)域添加等量的動量載荷,無法直接實(shí)現(xiàn)實(shí)際涵道風(fēng)扇槳盤平面上的動量載荷分布規(guī)律,為此需對FLUENT軟件的動量源施加模塊進(jìn)行二次開發(fā)[11],通過軟件提供的UDF接口編制相應(yīng)的用戶自定義函數(shù)[12],采用葉素理論對動量源施加區(qū)域內(nèi)的每一個單元進(jìn)行判定,根據(jù)單元中心點(diǎn)在風(fēng)扇槳盤中的位置施加相應(yīng)的動量載荷。
2.4多工況求解計算和結(jié)果處理
為了提高仿真分析的自動集成化,在求解階段基于Python語言編制了多工況自動迭代計算程序,實(shí)現(xiàn)了FLUENT求解器的外部程序調(diào)用和多工況的自動提交、結(jié)果分析以及反復(fù)迭代計算。工況計算過程中,Python程序首先生成用于設(shè)置相關(guān)工況參數(shù)的journal批處理文件,再通過batch命令調(diào)用FLUENT求解器提交journal文件進(jìn)行計算,同時在計算過程中Python程序持續(xù)對計算進(jìn)程進(jìn)行監(jiān)控,一旦計算完成便對輸出的結(jié)果進(jìn)行處理分析。
3.1穩(wěn)態(tài)多工況計算
考慮到箱內(nèi)氣壓和試驗(yàn)段風(fēng)速對環(huán)境模擬裝置氣體環(huán)流狀態(tài)影響最大,結(jié)合火星大氣環(huán)境模擬裝置的使用工況,選取了多種典型工況對環(huán)境模擬裝置內(nèi)部氣體環(huán)流進(jìn)行了計算仿真,具體工況參數(shù)如表2所示。
環(huán)境模擬裝置典型工況下的速度矢量分布和壓力分布如圖3和圖4所示,裝置內(nèi)部流場流線分布如圖5所示。由圖可以看出試驗(yàn)段基本都處在等速層流區(qū)域,容器艙體內(nèi)部氣體由艙體四周向風(fēng)扇一側(cè)回流,并在中間層流區(qū)的四周形成渦旋,風(fēng)扇槳盤進(jìn)風(fēng)口附近會形成一定的負(fù)壓區(qū)域,這是風(fēng)扇吸力作用導(dǎo)致的,同時在風(fēng)扇加速氣流與艙體內(nèi)壁接觸的地方會出現(xiàn)一定的正壓區(qū),這是由于氣流在該區(qū)域受到阻滯速度降低從而壓力升高造成的。其余部分,特別是試驗(yàn)段壓力維持在比較均勻的水平?;趧恿吭捶椒ǖ娘L(fēng)扇模擬方式可以有效地模擬風(fēng)扇細(xì)節(jié)流場特征,同時減少了網(wǎng)格數(shù)量,有利于實(shí)際工程設(shè)計應(yīng)用。
表2 計算工況參數(shù)表
圖3 環(huán)境模擬裝置典型工況下的速度矢量分布Fig.3 Velocity vector distribution of environment simulation devices under typical operating conditions
圖4 環(huán)境模擬裝置典型工況下的壓力分布Fig.4 Pressure distribution of environment simulation devices under typical operating conditions
圖5 環(huán)境模擬裝置回流氣流流線分布Fig.5 Streamline distribution of reflux airflow of the environment simulation devices
根據(jù)環(huán)境模擬器箱多工況穩(wěn)態(tài)分析結(jié)果可知,在所有工況下試驗(yàn)段基本都處在等速層流區(qū)域,箱內(nèi)氣體由箱體四周向風(fēng)扇一側(cè)回流,并在中間層流區(qū)的四周形成渦旋,由于試驗(yàn)段完全處于層流區(qū)域內(nèi),因此計算過程中的速度監(jiān)控曲線在達(dá)到一定的迭代步數(shù)之后趨于穩(wěn)定并收斂。
由各工況壓力分布圖可知,風(fēng)扇槳盤進(jìn)風(fēng)口附近會形成一定的負(fù)壓區(qū)域,是風(fēng)扇吸力作用導(dǎo)致的,同時在風(fēng)扇加速氣流與箱壁接觸的地方會出現(xiàn)一定的正壓區(qū),這是由于氣流在該區(qū)域受到阻滯速度降低從而壓力升高造成的。其余部分,特別是試驗(yàn)段壓力維持在比較均勻的水平。由于計算模型沒有添加熱源,因此整個環(huán)境箱內(nèi)部流場溫度都保持在初始溫度上。
圖6 不同工況對應(yīng)風(fēng)扇需用功率變化趨勢Fig.6 Demand power changing tendency of driving fan under different operating conditions
根據(jù)流態(tài)仿真計算得到的不同工況驅(qū)動風(fēng)扇需用功率如圖6所示,由于環(huán)境箱內(nèi)氣壓很低,因此風(fēng)扇需用功率較小,最大僅達(dá)到了3.5 W,從變化趨勢上可以看出,驅(qū)動風(fēng)扇的需用功率隨著環(huán)境模擬器內(nèi)部壓力和試驗(yàn)段指定風(fēng)速的增加而增加,同時呈現(xiàn)出顯著地非線性特征。
3.2非穩(wěn)態(tài)典型工況計算
實(shí)際試驗(yàn)過程中需要對一些電子設(shè)備進(jìn)行測試,這些電子設(shè)備在測試過程中會產(chǎn)生一定的熱量,從而引起自身溫度的升高。為了研究容器艙體內(nèi)部氣流對這些產(chǎn)熱設(shè)備散熱的影響,在穩(wěn)態(tài)分析的基礎(chǔ)上在模型中加入試驗(yàn)熱源,通過對火星表面大氣環(huán)境模擬裝置氣體環(huán)流的非穩(wěn)態(tài)計算仿真來得到一定發(fā)熱功率下熱源與流場溫度隨時間變化規(guī)律。
計算基于一種典型工況,具體參數(shù)如表3所示。
表3 典型工況參數(shù)
比熱與導(dǎo)熱系數(shù)參考金屬鋁的材料特性,試驗(yàn)熱源初始溫度與環(huán)境溫度一致,容器內(nèi)部流場循環(huán)在風(fēng)扇作用下保持穩(wěn)定。其中,試驗(yàn)熱源前端風(fēng)速保持在20 m/s左右,初始時刻試驗(yàn)熱源開始以100 W的功率產(chǎn)生熱量并在整個計算過程中保持不變。熱源產(chǎn)生的熱量一部分用于熱源本身升溫,另一部分通過容器內(nèi)部氣流和自身輻射與模擬裝置容器艙體壁面進(jìn)行熱量交換,其中模擬裝置容器艙體壁面在整個計算過程中保持恒溫。試驗(yàn)熱源位置如圖7所示。
圖7 環(huán)境模擬裝置內(nèi)部試驗(yàn)熱源位置Fig.7 Position of test heat source inside environment chamber
環(huán)境模擬裝置內(nèi)部流場速度矢量分布和壓力分布如圖8、圖9所示,可以看出,試驗(yàn)熱源前端氣流速度基本維持在20 m/s左右,達(dá)到了計算工況的要求,試驗(yàn)熱源后段形成了一個紊流區(qū),使得流場氣體回流點(diǎn)相對于第3.1節(jié)的穩(wěn)態(tài)計算結(jié)果向前移動,造成氣流循環(huán)路徑變短,進(jìn)而導(dǎo)致風(fēng)扇區(qū)氣流速度進(jìn)一步增加。由環(huán)境模擬裝置內(nèi)部壓力分布圖可以看出,由于試驗(yàn)熱源的阻滯作用,熱源前端形成了一個正壓區(qū),后端由于處于紊流狀態(tài),壓力分布較為均勻,由于風(fēng)扇區(qū)氣流速度增加,風(fēng)扇進(jìn)風(fēng)口處的負(fù)壓區(qū)壓力相對于穩(wěn)態(tài)計算結(jié)果進(jìn)一步降低,在風(fēng)扇槳尖處達(dá)到極值。
圖8 環(huán)境模擬裝置內(nèi)部流場速度矢量分布Fig.8 Interior flow velocity vector distribution of cylindrical andsquare environment simulation devices
圖9 環(huán)境模擬裝置內(nèi)部壓力分布Fig.9 Pressure distribution of environment simulation devices
環(huán)境模擬裝置內(nèi)部最終溫度分布和試驗(yàn)熱源溫度變化時間歷程曲線如圖10和圖11所示,可以看出,初始階段試驗(yàn)熱源升溫較為迅速,而隨著試驗(yàn)熱源和環(huán)境箱內(nèi)流體之間的溫差不斷增加,這兩者之間的換熱速率也不斷增加,從而造成熱源升溫速率逐漸下降,經(jīng)過2 000 s之后,試驗(yàn)熱源的溫度最終趨于穩(wěn)定,與環(huán)境箱內(nèi)流體的溫差保持在9℃左右,環(huán)境箱內(nèi)流體流過試驗(yàn)熱源表面之后會產(chǎn)生一定的升溫,因此在試驗(yàn)熱源后方沿流場流線方向會產(chǎn)生一定長度的高溫流體帶。
在非穩(wěn)態(tài)分析過程中,環(huán)境模擬裝置內(nèi)部試驗(yàn)熱源會對箱體內(nèi)部氣體環(huán)流產(chǎn)生顯著影響,在恒定發(fā)熱功率下,熱源升溫速率隨時間增加而不斷減小并最終達(dá)到一個穩(wěn)定狀態(tài)。
圖10 環(huán)境模擬裝置最終溫度分布Fig.10 Final temperature distribution of environment simulation devices
圖11 環(huán)境模擬裝置試驗(yàn)熱源平均溫度時間變化曲線Fig.11 Heat source temperature changing with time of environment simulation devices
本文采用三維CFD數(shù)值分析方法對不同狀態(tài)下的火星大氣模擬裝置內(nèi)部環(huán)流狀態(tài)進(jìn)行了穩(wěn)態(tài)和非穩(wěn)態(tài)模擬仿真,得到以下結(jié)論:
1)根據(jù)稀薄流體空氣動力學(xué),火星大氣裝置內(nèi)部流場屬于連續(xù)流場范圍,N-S方程在其所處的氣體條件范圍下是適用的。
2)基于動量源方法的風(fēng)扇模擬方式可以有效地模擬風(fēng)扇細(xì)節(jié)流場特征,同時減少了網(wǎng)格數(shù)量,有利于實(shí)際工程設(shè)計應(yīng)用。
3)通過對環(huán)境模擬裝置內(nèi)部環(huán)流狀態(tài)進(jìn)行穩(wěn)態(tài)模擬得到了不同工況下的驅(qū)動風(fēng)扇需用功率,對風(fēng)扇的后續(xù)設(shè)計提供了輸入?yún)?shù)。
4)在非穩(wěn)態(tài)分析過程中,環(huán)境模擬裝置內(nèi)部試驗(yàn)熱源會對箱體內(nèi)部氣體環(huán)流產(chǎn)生顯著影響,在恒定發(fā)熱功率下,熱源升溫速率隨時間增加而不斷減小并最終達(dá)到一個穩(wěn)定狀態(tài)。
References)
[1]達(dá)道安,楊亞天,涂建輝. 太陽系行星及行星際大氣環(huán)境特性研究[J]. 宇航學(xué)報,2006,27(6):1306-1313.
DA D A, YANG Y T, TU J H, Research on the environment characters of atmosphere of planetary and interplanetary space in the solar system[J]. Journal of Astronautics, 2006,27(6): 1306-1313(in Chinese).
[2]張文普,韓波,張成義. 大氣制動期間探測器的氣動特性和軌道計算[J]. 應(yīng)用數(shù)學(xué)和力學(xué), 2010, 31(9): 1016-1026.
ZHANG W P, HAN B, ZHANG C Y, Spacecraft aerodynamics and trajectory simulation during aerobraking[J]. Applied Mathematics and Mechanics, 2010, 31(9): 1016-1026(in Chinese).
[3]VOTTA R, SCHETTINO A, RANUZZI G, et al. Hypersonic low-density aerothermodynamics of Orion-like exploration vehicle[J]. Journal of Spacecraft and Rockets, 2009, 46(4): 781-787.
[4]WAY D, DESAI P, ENGELUNG W, et al. Design and analysis of the drop test vehicle for the Mars exploration rover parachute structural tests[J]. AIAA Paper, 2003, 2138: 2003.
[5]于瑩瀟,田佳林.火星探測器降落傘系統(tǒng)綜述[J]. 航天返回與遙感, 2007, 28(4): 12-16.
YUY X, TIAN J L. Mars explorer′s parachute system overview[J]. Spacecraft Recovery & Remote Sensing, 2007, 28(4): 12-16(in Chinese).
[6]李志輝,張涵信. 稀薄流到連續(xù)流的氣體運(yùn)動論統(tǒng)一算法研究[J]. 空氣動力學(xué)學(xué)報, 2003, 21(3):255-266.
LI Z H, ZHANG H X. Study on gas kinetic unified algorithm for flows from rarefied transition to continuum[J]. Acta Aerodynamica Sinica, 2003, 21(3): 255-266(in Chinese).
[7]沈青著.稀薄氣體動力學(xué)[M]. 北京:國防工業(yè)出版社, 2003:77-78.
SHEN Q. Rarefied gas dynamics[M]. Beijing: National Defense Industry Press, 2003:77-78(in Chinese).
[8]常勇強(qiáng),曹子棟,趙振興,等. 多組分氣體熱物性參數(shù)的計算方法[J]. 動力工程學(xué)報, 2010, 30(10): 772-776.
CHANG Y Q, CAO Z D,ZHAO Z X, et al. Calculation method for thermal properties of multi-component gas[J]. Chinese Journal of Power Engineering, 2010, 30(10): 772-776(in Chinese).
[9]TOLSON R H, DWYER A M, HANNA J L, et al. Application of accelerometer data to Mars Odyssey aerobraking and atmospheric modeling[J]. Journal of Spacecraft and Rockets, 2005, 42(3): 435-443.
[10]宋長紅,林永峰,陳文軒,等. 基于動量源方法的涵道尾槳CFD分析[J]. 直升機(jī)技術(shù), 2009(1):6-11.
SONG C H, LIN Y F,CHEN W X, et al. CFD analysis for the ducted tail rotor based on momentum-source method[J]. Helicopter Technique, 2009(1): 6-11(in Chinese).
[11]肖虹,高超,黨云卿,等. FLUENT軟件的二次開發(fā)及其在火箭氣動計算中的應(yīng)用[J]. 航空計算技術(shù),2009, 39(5): 55-57.
XIAO H, GAO C,DANG Y Q, et al. Secondary development of FLUENT and application in numerical simulation of aerodynamic characteristics for rockets[J]. Aeronautical Computing Technique, 2009,39(5): 55-57(in Chinese).
[12]Fluent Inc. Fluent manual. 2010[EB/OL].[2016-02-01].http:∥148.204.81.206/Ansys/150/ANSYS%20Fluent%20User%20Guide.pdf.
(編輯:高珍)
Analysis of design, simulation & calculation module for Martian atmosphere environment simulation device
LIU Gaotong*,SUN Yu,ZHANG Lei
Beijing Institute of Spacecraft Environment Engineering,Beijing 100094,China
The characteristics of Martian atmosphere were studied.A simulation calculation module of the Martian atmosphere simulation environment box was set up by choosing appropriate calculation method combined with FLUENT fluid finite element calculation software.The module realizes steady and unsteady calculation simulation of interior flow field under multiple working conditions of the cylindrical Martian atmosphere environment simulation devices.The calculation results was analyzed so as to provide technical support and reference for design of the Martian atmosphere environment simulation device.
Martian atmosphere; rarefied gas; fluid simulation; environment simulation device;computational fluid dynamics
10.16708/j.cnki.1000-758X.2016.0059
2016-03-03;
2016-04-27;錄用日期:2016-08-22;
時間:2016-09-2113:41:25
http:∥www.cnki.net/kcms/detail/11.1859.V.20160921.1341.007.html
劉高同(1986-),男,工程師,lgt_0903@163.com,研究方向?yàn)榭臻g環(huán)境模擬控制技術(shù)研究及航天器試驗(yàn)裝備研制
LIUGT,SUNY,ZHANGL.Analysisofdesign,simulation&calculationmoduleforMartianatmosphereenvironmentsimulationdevice[J].ChineseSpaceScienceandTechnology, 2016,36(5):65-71(inChinese).
V416.5
A
http:∥zgkj.cast.cn
引用格式:劉高同,孫宇,張磊.火星大氣環(huán)境模擬裝置設(shè)計及仿真分析研究[J].中國空間科學(xué)技術(shù), 2016,36(5):65-71.