火星大氣環(huán)境模擬裝置設(shè)計及仿真分析研究

劉高同，孫宇，張磊

北京衛(wèi)星環(huán)境工程研究所，北京 100094

?

火星大氣環(huán)境模擬裝置設(shè)計及仿真分析研究

劉高同*，孫宇，張磊

北京衛(wèi)星環(huán)境工程研究所，北京 100094

對火星表面大氣環(huán)境特性進行了研究，通過選取合適的計算方法并結(jié)合FLUENT流體有限元計算軟件對火星表面稀薄氣體內(nèi)部環(huán)流進行了模擬仿真分析，提出了以動量源模擬風扇段內(nèi)流的仿真方法，并進行了可行性討論。進一步實現(xiàn)了針對圓柱形模擬裝置多工況下的內(nèi)部氣體流場穩(wěn)態(tài)和非穩(wěn)態(tài)計算仿真，并對計算結(jié)果進行了分析討論，為火星大氣環(huán)境模擬裝置的設(shè)計提供了技術(shù)支持和參考。

火星大氣；稀薄氣體；流體仿真；環(huán)境模擬裝置；計算流體動力學

火星探測計劃以及火星采樣返回、載人登陸和火星基地建設(shè)等可能的后續(xù)任務(wù)，對準確預測探測器在火星大氣環(huán)境中的氣動特性提出了嚴格要求。陌生的大氣環(huán)境、技術(shù)手段的局限和基礎(chǔ)數(shù)據(jù)的缺乏對氣動特性預測帶來了極大的困難?；鹦潜砻嫦”〈髿獬煞峙c地球大氣顯著不同，密度也只有地球大氣的1%，因此飛行器進入火星大氣與再入地球大氣的過程及遭遇的氣動環(huán)境非常不同[1]。目前,針對火星探測器氣動特性的研究主要集中在火星軌道飛行器氣動減速[2]、飛行器進入大氣層高速降落和飛行器傘降等方面[3-5]，針對火星近地表面探測器低速情況下的氣動特性研究較少，尤其是試驗研究國內(nèi)尚無人開展，因此，設(shè)計建造火星大氣環(huán)境模擬裝置，針對火星表面稀薄大氣流態(tài)進行模擬，進而研究火星表面大氣環(huán)境對探測器周圍流場結(jié)構(gòu)和氣動特性的影響，對于掌握火星探測進入階段的飛行器氣動特性及探測器外形設(shè)計，具有重要科學意義與工程價值。本文結(jié)合火星大氣特性選取適當?shù)挠嬎隳Ｐ歪槍鹦潜砻娲髿猸h(huán)境模擬裝置進行了內(nèi)部環(huán)流仿真計算，對裝置內(nèi)部流場進行了計算仿真，討論了不同構(gòu)型裝置在設(shè)計工況下的性能。

1　火星表面稀薄氣體流態(tài)仿真分　 析方法

由于火星表面大氣遠比地球大氣稀薄，因此在模擬計算方法選取的過程中需考慮火星大氣的稀薄特性。稀薄氣體動力學理論中一般采用努森數(shù)來表征某個氣象流的稀薄程度[6]，努森數(shù)即氣體分子平均自由程λ與流動特征長度L的比值，其具體表達式為：

(1)

根據(jù)給定的溫度和壓力范圍計算得到火星表面稀薄氣體分子平均自由程最大值為8.357×10-5m，最小值為2.235×10-6m，取氣體流動宏觀量梯度的標尺長度為1 m，努森數(shù)范圍為2.235×10-6～8.357×10-5，根據(jù)錢學森關(guān)于氣體流動區(qū)域的劃分[7]，根據(jù)稀薄流體空氣動力學，火星大氣裝置內(nèi)部流場屬于連續(xù)流場范圍，N-S方程在其所處的氣體條件范圍下是適用的。

2　火星大氣環(huán)境模擬裝置仿真分　 析模型

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ū)域，用于模擬涵道風扇，紅色區(qū)域為試驗段，計算過程中對試驗段風速進行監(jiān)控，保證該區(qū)域平均風速達到計算工況的要求。

圖2　圓形裝置流體計算網(wǎng)格劃分Fig.2　Fluid calculation grid of cylindrical environment simulation device

2.2內(nèi)部環(huán)流流體物性參數(shù)設(shè)定

火星表面稀薄大氣約含二氧化碳95.3%、氮氣2.7%和氬氣1.6%，還有微量的氧氣、一氧化碳、氖氣、氪氣和氙氣，在模塊流體參數(shù)設(shè)定階段，需根據(jù)這幾種氣體的質(zhì)量分數(shù)以及流場壓強和溫度合成火星大氣的物性參數(shù)，在現(xiàn)有氣體組分熱物性質(zhì)的基礎(chǔ)上，本文結(jié)合相關(guān)文獻[8]，研究了多組分氣體熱物性參數(shù)的計算方法，并取了合適的計算模型，對混合氣體的密度、比熱、粘度、導熱系數(shù)可采用以下公式進行計算。

(1)混合氣體密度

混合氣體密度主要根據(jù)各組分氣體的體積分數(shù)及其在制定壓力和溫度條件下的密度計算得到。具體計算公式為：

(2)

式中：ρi為各組分氣體在指定氣壓和溫度下的密度；yi為各組分氣體在混合氣體中的濃度。

(2)混合氣體黏度

混合氣體的黏度可以通過各組分的純物質(zhì)黏度、相對分子質(zhì)量及濃度，根據(jù)一定的混合規(guī)則求得。嚴格的Chapaman-Enskog動力論可以推廣用于計算多元氣體混合物的黏度。經(jīng)簡化，略去二階的影響，則嚴格的數(shù)值解可近似地用級數(shù)表示：

(3)

式中：yi為組分氣體濃度；ηi為組分氣體的黏度；φij為組分氣體i和組分氣體j的結(jié)合因子。

(3)混合氣體導熱系數(shù)

文獻中提供了很多混合氣體導熱系數(shù)計算的方法，本文采用以下公式進行計算：

(4)

式中：λ混合為混合氣體的導熱系數(shù)；λi為組分氣體i的導熱系數(shù)。

(4)混合氣體的定壓比熱容

混合氣體密度主要根據(jù)各組分氣體的濃度及其在制定壓力和溫度條件下的定壓比熱容計算得到。具體計算公式為：

(5)

式中：cpi為各組分氣體在指定氣壓和溫度下的定壓比熱容。

2.3涵道風扇模擬

涵道風扇是引起模擬裝置內(nèi)氣體環(huán)流的驅(qū)動器，涵道風扇模擬的準確度直接影響到容器艙內(nèi)氣體環(huán)流的精度?？紤]到目前涵道風扇并未開始詳細設(shè)計，無法提供可用于CFD分析計算的的風扇細節(jié)幾何模型，分析計算中將風扇所在區(qū)域簡化為一個動量源[10]，以動量源代替槳葉對流場的作用來模擬涵道風扇的氣動特性。在連續(xù)流場中關(guān)于坐標軸i的動量方程可形成如下形式：

(6)

式中：ρ為流體密度；ui為流體沿坐標軸i方向的流動速度；p為流體壓力；Si表示風扇沿坐標軸i方向上的單位體積作用力，即為動量源項。

FLUENT軟件內(nèi)部只能針對指定區(qū)域添加等量的動量載荷，無法直接實現(xiàn)實際涵道風扇槳盤平面上的動量載荷分布規(guī)律，為此需對FLUENT軟件的動量源施加模塊進行二次開發(fā)[11]，通過軟件提供的UDF接口編制相應的用戶自定義函數(shù)[12]，采用葉素理論對動量源施加區(qū)域內(nèi)的每一個單元進行判定，根據(jù)單元中心點在風扇槳盤中的位置施加相應的動量載荷。

2.4多工況求解計算和結(jié)果處理

為了提高仿真分析的自動集成化，在求解階段基于Python語言編制了多工況自動迭代計算程序，實現(xiàn)了FLUENT求解器的外部程序調(diào)用和多工況的自動提交、結(jié)果分析以及反復迭代計算。工況計算過程中，Python程序首先生成用于設(shè)置相關(guān)工況參數(shù)的journal批處理文件，再通過batch命令調(diào)用FLUENT求解器提交journal文件進行計算，同時在計算過程中Python程序持續(xù)對計算進程進行監(jiān)控，一旦計算完成便對輸出的結(jié)果進行處理分析。

3　火星大氣環(huán)境模擬裝置氣體流　 態(tài)仿真分析

3.1穩(wěn)態(tài)多工況計算

考慮到箱內(nèi)氣壓和試驗段風速對環(huán)境模擬裝置氣體環(huán)流狀態(tài)影響最大，結(jié)合火星大氣環(huán)境模擬裝置的使用工況，選取了多種典型工況對環(huán)境模擬裝置內(nèi)部氣體環(huán)流進行了計算仿真，具體工況參數(shù)如表2所示。

環(huán)境模擬裝置典型工況下的速度矢量分布和壓力分布如圖3和圖4所示，裝置內(nèi)部流場流線分布如圖5所示。由圖可以看出試驗段基本都處在等速層流區(qū)域，容器艙體內(nèi)部氣體由艙體四周向風扇一側(cè)回流，并在中間層流區(qū)的四周形成渦旋，風扇槳盤進風口附近會形成一定的負壓區(qū)域，這是風扇吸力作用導致的，同時在風扇加速氣流與艙體內(nèi)壁接觸的地方會出現(xiàn)一定的正壓區(qū)，這是由于氣流在該區(qū)域受到阻滯速度降低從而壓力升高造成的。其余部分，特別是試驗段壓力維持在比較均勻的水平?；趧恿吭捶椒ǖ娘L扇模擬方式可以有效地模擬風扇細節(jié)流場特征，同時減少了網(wǎng)格數(shù)量，有利于實際工程設(shè)計應用。

表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é)果可知，在所有工況下試驗段基本都處在等速層流區(qū)域，箱內(nèi)氣體由箱體四周向風扇一側(cè)回流，并在中間層流區(qū)的四周形成渦旋，由于試驗段完全處于層流區(qū)域內(nèi)，因此計算過程中的速度監(jiān)控曲線在達到一定的迭代步數(shù)之后趨于穩(wěn)定并收斂。

由各工況壓力分布圖可知，風扇槳盤進風口附近會形成一定的負壓區(qū)域，是風扇吸力作用導致的，同時在風扇加速氣流與箱壁接觸的地方會出現(xiàn)一定的正壓區(qū)，這是由于氣流在該區(qū)域受到阻滯速度降低從而壓力升高造成的。其余部分，特別是試驗段壓力維持在比較均勻的水平。由于計算模型沒有添加熱源，因此整個環(huán)境箱內(nèi)部流場溫度都保持在初始溫度上。

圖6　不同工況對應風扇需用功率變化趨勢Fig.6　Demand power changing tendency of driving fan under different operating conditions

根據(jù)流態(tài)仿真計算得到的不同工況驅(qū)動風扇需用功率如圖6所示，由于環(huán)境箱內(nèi)氣壓很低，因此風扇需用功率較小，最大僅達到了3.5 W，從變化趨勢上可以看出，驅(qū)動風扇的需用功率隨著環(huán)境模擬器內(nèi)部壓力和試驗段指定風速的增加而增加，同時呈現(xiàn)出顯著地非線性特征。

3.2非穩(wěn)態(tài)典型工況計算

實際試驗過程中需要對一些電子設(shè)備進行測試，這些電子設(shè)備在測試過程中會產(chǎn)生一定的熱量，從而引起自身溫度的升高。為了研究容器艙體內(nèi)部氣流對這些產(chǎn)熱設(shè)備散熱的影響，在穩(wěn)態(tài)分析的基礎(chǔ)上在模型中加入試驗熱源，通過對火星表面大氣環(huán)境模擬裝置氣體環(huán)流的非穩(wěn)態(tài)計算仿真來得到一定發(fā)熱功率下熱源與流場溫度隨時間變化規(guī)律。

計算基于一種典型工況，具體參數(shù)如表3所示。

表3　典型工況參數(shù)

比熱與導熱系數(shù)參考金屬鋁的材料特性，試驗熱源初始溫度與環(huán)境溫度一致，容器內(nèi)部流場循環(huán)在風扇作用下保持穩(wěn)定。其中，試驗熱源前端風速保持在20 m/s左右，初始時刻試驗熱源開始以100 W的功率產(chǎn)生熱量并在整個計算過程中保持不變。熱源產(chǎn)生的熱量一部分用于熱源本身升溫，另一部分通過容器內(nèi)部氣流和自身輻射與模擬裝置容器艙體壁面進行熱量交換，其中模擬裝置容器艙體壁面在整個計算過程中保持恒溫。試驗熱源位置如圖7所示。

圖7　環(huán)境模擬裝置內(nèi)部試驗熱源位置Fig.7　Position of test heat source inside environment chamber

環(huán)境模擬裝置內(nèi)部流場速度矢量分布和壓力分布如圖8、圖9所示，可以看出，試驗熱源前端氣流速度基本維持在20 m/s左右，達到了計算工況的要求，試驗熱源后段形成了一個紊流區(qū)，使得流場氣體回流點相對于第3.1節(jié)的穩(wěn)態(tài)計算結(jié)果向前移動，造成氣流循環(huán)路徑變短，進而導致風扇區(qū)氣流速度進一步增加。由環(huán)境模擬裝置內(nèi)部壓力分布圖可以看出，由于試驗熱源的阻滯作用，熱源前端形成了一個正壓區(qū)，后端由于處于紊流狀態(tài)，壓力分布較為均勻，由于風扇區(qū)氣流速度增加，風扇進風口處的負壓區(qū)壓力相對于穩(wěn)態(tài)計算結(jié)果進一步降低，在風扇槳尖處達到極值。

圖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)部最終溫度分布和試驗熱源溫度變化時間歷程曲線如圖10和圖11所示，可以看出，初始階段試驗熱源升溫較為迅速，而隨著試驗熱源和環(huán)境箱內(nèi)流體之間的溫差不斷增加，這兩者之間的換熱速率也不斷增加，從而造成熱源升溫速率逐漸下降，經(jīng)過2 000 s之后，試驗熱源的溫度最終趨于穩(wěn)定，與環(huán)境箱內(nèi)流體的溫差保持在9℃左右，環(huán)境箱內(nèi)流體流過試驗熱源表面之后會產(chǎn)生一定的升溫，因此在試驗熱源后方沿流場流線方向會產(chǎn)生一定長度的高溫流體帶。

在非穩(wěn)態(tài)分析過程中，環(huán)境模擬裝置內(nèi)部試驗熱源會對箱體內(nèi)部氣體環(huán)流產(chǎn)生顯著影響，在恒定發(fā)熱功率下，熱源升溫速率隨時間增加而不斷減小并最終達到一個穩(wěn)定狀態(tài)。

圖10　環(huán)境模擬裝置最終溫度分布Fig.10　Final temperature distribution of environment simulation devices

圖11　環(huán)境模擬裝置試驗熱源平均溫度時間變化曲線Fig.11　Heat source temperature changing with time of environment simulation devices

4　結(jié)束語

本文采用三維CFD數(shù)值分析方法對不同狀態(tài)下的火星大氣模擬裝置內(nèi)部環(huán)流狀態(tài)進行了穩(wěn)態(tài)和非穩(wěn)態(tài)模擬仿真，得到以下結(jié)論：

1)根據(jù)稀薄流體空氣動力學，火星大氣裝置內(nèi)部流場屬于連續(xù)流場范圍，N-S方程在其所處的氣體條件范圍下是適用的。

2)基于動量源方法的風扇模擬方式可以有效地模擬風扇細節(jié)流場特征，同時減少了網(wǎng)格數(shù)量，有利于實際工程設(shè)計應用。

3)通過對環(huán)境模擬裝置內(nèi)部環(huán)流狀態(tài)進行穩(wěn)態(tài)模擬得到了不同工況下的驅(qū)動風扇需用功率，對風扇的后續(xù)設(shè)計提供了輸入?yún)?shù)。

4)在非穩(wěn)態(tài)分析過程中，環(huán)境模擬裝置內(nèi)部試驗熱源會對箱體內(nèi)部氣體環(huán)流產(chǎn)生顯著影響，在恒定發(fā)熱功率下，熱源升溫速率隨時間增加而不斷減小并最終達到一個穩(wěn)定狀態(tài)。

References)

[1]達道安,楊亞天,涂建輝. 太陽系行星及行星際大氣環(huán)境特性研究[J]. 宇航學報,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]. 應用數(shù)學和力學, 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ù)流的氣體運動論統(tǒng)一算法研究[J]. 空氣動力學學報, 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]沈青著.稀薄氣體動力學[M]. 北京：國防工業(yè)出版社, 2003:77-78.

SHEN Q. Rarefied gas dynamics[M]. Beijing: National Defense Industry Press, 2003:77-78(in Chinese).

[8]常勇強,曹子棟,趙振興,等. 多組分氣體熱物性參數(shù)的計算方法[J]. 動力工程學報, 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]. 直升機技術(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ā)及其在火箭氣動計算中的應用[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，研究方向為空間環(huán)境模擬控制技術(shù)研究及航天器試驗裝備研制

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].中國空間科學技術(shù), 2016,36(5)：65-71.