邢好運 劉 卓 汪 球 趙 偉 高亮杰 劉中臣 錢戰(zhàn)森

* (中國科學(xué)院力學(xué)研究所高溫氣體動力學(xué)國家重點實驗室,北京 100190)

? (中國科學(xué)院大學(xué)工程科學(xué)學(xué)院,北京 100049)

** (中國航空工業(yè)空氣動力研究院,沈陽 110034)

引言

由于火星的大氣條件和地球比較相似,火星探測一直是國際深空探測發(fā)展的熱點.截止目前,全球各國一共向火星進行了48 次發(fā)射,但在所有的發(fā)射中,實現(xiàn)真正意義上軟著陸的國家只有中國與美國.在所有著陸嘗試中,有10 次完全成功的火星探測器軟著陸,成功的概率僅約為50%,因此火星也被稱為“探測器墳場”[1-3].特殊且缺乏足夠數(shù)據(jù)庫支撐的火星大氣環(huán)境及由此引起的進入飛行氣動特性不確定性是造成探測器飛行安全隱患的主要原因之一,亟待氣動力/熱的精確預(yù)測來提升飛行器氣動布局和熱防護系統(tǒng)設(shè)計.

火星大氣與地球大氣除氣體組分、密度等不同外,還有一個特別之處在于火星大氣中存在微小的塵埃顆粒,主要由鐵氧化物和硅氧化物組成.在火星平均每3～4 年就會發(fā)生一次全球級的沙塵暴,灰塵顆粒甚至能達(dá)到60 km 的高空,根據(jù)Viking 著陸器的觀測以及對火星大氣中垂直風(fēng)強度的估計,大尺寸顆粒(直徑約為5～10 μm)在一次大型沙暴開始后的20～50 d 內(nèi)依舊會懸浮在大氣上空[4],而在高超聲速來流中,即使顆粒的質(zhì)量濃度較低,顆粒的存在也會導(dǎo)致探測器表面熱流顯著增加[5].因此,開展含顆?；鹦谴髿猸h(huán)境對飛行器氣動特性的影響研究十分有必要,它有助于降低防熱結(jié)構(gòu)設(shè)計的冗余量或提升飛行器性能.

塵埃顆粒對高超聲速飛行器火星進入過程中的影響研究最早是由Papadopoulos 等[6]通過數(shù)值模擬進行評估,其結(jié)果表明火星沙塵暴會對飛行器熱防護系統(tǒng)(thermal protection system,TPS)表面造成嚴(yán)重侵蝕;然而,Palmer 等[7]對此重新評估,認(rèn)為其高估了由于顆粒撞擊造成的附加質(zhì)量損失,兩者的主要差異源于他們使用了不同的塵埃顆粒模型.Majid等[8]模擬了(mars sample return orbiter,MSRO)在塵埃顆粒環(huán)境下的進入過程,結(jié)果表明顆粒撞擊表面引起的熱流密度與對流熱相比可以忽略不計;Vasilevskii等[5]發(fā)現(xiàn)在無擾動流中添加直徑為0.15 μm 的微小顆粒,即使在濃度較低(約為1%)時,模型臨界點區(qū)域的熱流通量也會顯著增加.Palmer 等[9]最近也對這個問題進行了全面的概述,并模擬了ExoMars Schiaparelli 飛行器在3 種不同的沙塵暴條件下的進入過程,其結(jié)果表明,在全球性沙塵暴發(fā)生時,駐點處的隔熱層侵蝕量約為2.1 mm,占TPS 厚度的17%,而在區(qū)域性沙塵暴發(fā)生時侵蝕量僅為0.35 mm,在沙塵靜止條件下的侵蝕量可以忽略不計.另外,Ching等[10]發(fā)現(xiàn)熱通量預(yù)測對阻力系數(shù)非常敏感;且Ching 等[11-12]最近的工作也表明,使用不同的阻力模型和自由流顆粒分布時,進入過程中飛行器上某些位置的表面衰退出現(xiàn)顯著差異.雖然目前學(xué)者已經(jīng)提出了模擬該問題的一些方法,但是微尺寸顆粒相關(guān)的理論模型還不完備,不同的研究人員的結(jié)論之間也存在分歧,有必要進一步研究含顆粒環(huán)境下高超聲速飛行器火星進入過程的氣動問題.

本文基于Euler-Lagrange 框架,采用單向耦合方法,模擬ExoMars Schiaparelli 進入艙在塵埃顆粒環(huán)境下的火星進入過程,考慮高溫相變模型對不同粒徑顆粒運動軌跡的影響,模擬不同尺寸顆粒的運動軌跡以及計算對應(yīng)的撞擊分?jǐn)?shù),并選取模態(tài)半徑為0.35 μm 的火星塵埃顆粒分布模型來計算撞擊能量分?jǐn)?shù),相關(guān)研究有助于理解塵埃顆粒的運動規(guī)律,并可幫助建立塵埃顆粒對飛行器表面的侵蝕模型.

1 顆粒運動模型

當(dāng)考慮火星大氣含灰時,數(shù)值模擬即變成一個兩相流問題,需要考慮顆粒與流場之間的耦合方式,一般通過分散相的體積分?jǐn)?shù)αp來選擇適當(dāng)?shù)鸟詈戏绞?包括單向耦合、雙向耦合以及四向耦合[13].火星大氣塵埃顆粒的體積分?jǐn)?shù)αp約為10-6量級[14],含量較低,因此本文的計算中選用單向耦合方法來進行計算,即僅考慮流場對顆粒的作用而不考慮顆粒對流場的作用,單向耦合也是國內(nèi)外學(xué)者研究這個問題最常用耦合方式[6,9,15-16].另外,本文所做的主要假設(shè)還包括:

(1) 火星塵埃為離散的固體球形顆粒,密度為2940 kg/m3[6];

(2) 只考慮顆粒受流場對流傳熱影響,不考慮其受到的輻射加熱;

(3) 顆粒自身不存在溫度梯度,且顆粒在流場中一直為球形,在模擬過程中,相變(蒸發(fā)、熔化或升華)均勻,顆粒不發(fā)生碎裂;

(4) 假定質(zhì)點運動僅為平動,不考慮Magnus力、附加質(zhì)量力、Basset 力、Saffman 升力和其他外力(重力、電磁力);

(5) 不考慮顆粒間的碰撞.

1.1 顆粒動力學(xué)模型

顆粒運動采用拉格朗日方法進行分析,可求得顆粒軌跡上的具體信息.對顆粒運動使用牛頓第二定律進行計算,僅考慮流場作用在顆粒上的阻力

其中,ΔV為顆粒所在位置流場速度Vg與顆粒自身速度Vp的差值,即ΔV=Vg-Vp,mp為顆粒質(zhì)量,Cd為顆粒阻力系數(shù),rp為顆粒半徑,ρg為顆粒所在位置流體的密度,進一步由式(1)可得

其中,ρp為顆粒密度,dp為顆粒直徑.本文阻力系數(shù)模型采用Henderson[17]提出的阻力模型,該模型適用于連續(xù)、過渡和自由分子流狀態(tài)的球體阻力系數(shù)的計算,同時還考慮了相對馬赫數(shù)Mrel、相對雷諾數(shù)Rerel以及顆粒自身溫度Tp與流場溫度Tg的影響,Henderson 模型與Bailey 等[18]在球體上的實驗數(shù)據(jù)以及20 世紀(jì)初的理論推導(dǎo)保持一致,Henderson阻力模型計算如下.

對于Mrel＜ 1

最后通過Lagrange 方法來計算顆粒的位置,計算方法如下

由于本文開展的是二維模擬,因此式中Xp=[xp,yp],Vp=[up,νp].Xp代表顆粒所在位置,Vp代表顆粒速度,xp和yp為顆粒在x,y方向的位置分量,up和νp為顆粒在x,y方向的速度分量.

1.2 顆粒熱力學(xué)模型

在飛行器進入火星大氣的過程中,激波后流場溫度可達(dá)上千度,溫度可能會高于顆粒的融化溫度,使顆粒發(fā)生相變.此外,顆粒的溫度還會影響阻力系數(shù)的計算.因此,本文將考慮流場向顆粒傳熱導(dǎo)致其溫度升高及相變發(fā)生,具體的計算公式如下

式中,Ch為對流傳熱系數(shù),通過式(8),可以得到顆粒表面溫度的計算式為

要計算顆粒表面溫度Tp的變化,需要確定顆粒的對流傳熱系數(shù)Ch,其與Nusselt 數(shù)(Nu) 相關(guān),Nu表示在邊界處對流換熱與熱傳導(dǎo)的比值

對于顆粒來說,特征長度L為顆粒直徑,κg為顆粒周圍流體的熱傳導(dǎo)系數(shù).Nusselt 數(shù)表達(dá)式通常用顆粒馬赫數(shù)和雷諾數(shù)表示,因此也是流體速度、黏性和流體導(dǎo)熱系數(shù)的函數(shù).本文中使用Fox 等[19]建立的Nusselt 數(shù)計算模型,該模型考慮了可壓縮和非連續(xù)流動影響,且該式適用于亞聲速和超聲速顆粒馬赫數(shù)

顆粒進入激波層之后會吸熱,表面溫度將升高,一旦超過顆粒自身的融化溫度將會導(dǎo)致顆粒直徑的減小,火星大氣中顆粒的成分一般為SiO2,熔化溫度約為1990 K,在單向耦合的假設(shè)下,僅考慮顆粒自身直徑的變化而不考慮由顆粒相變導(dǎo)致的流場質(zhì)量增加,具體的計算模型如下

式中,ζ為相變潛熱,SiO2的相變潛熱為8.6 MJ/kg;Tvapor為汽化溫度,對于球形顆粒,式(12)可以轉(zhuǎn)化為

需要注意的是,顆粒溫度首先達(dá)到的是熔化溫度,即1990 K,之后達(dá)到Tvapor才會蒸發(fā),而在熔化溫度與蒸發(fā)溫度之間時,假設(shè)液態(tài)的部分附著在固態(tài)表面,顆粒質(zhì)量不變;顆粒的汽化溫度Tvapor隨表面壓力的增加而升高.本文采用Schaefer 等[20]提出的汽化溫度線性擬合模型,該模型與參考文獻[21]中給出的實驗數(shù)據(jù)相吻合,在該模型中,氣體壓力p的單位為bar (1 bar=1.0×105Pa)

2 顆粒定位及流場信息插值方法

2.1 顆粒定位方法

使用歐拉-拉格朗日方法,我們需要確定顆粒在每個拉格朗日時間步長時所處的網(wǎng)格位置.無論是使用單向耦合還是雙向耦合,拉格朗日計算都需要這些信息,這是由于計算需要在每個拉格朗日時間步長將流場的性質(zhì)插值到當(dāng)前顆粒位置.本文使用文獻[8]中的顆粒定位方法,為了確定一個顆粒是否在一個給定的網(wǎng)格內(nèi),從逆時針方向?qū)W(wǎng)格的各節(jié)點進行編號,從顆粒到每個網(wǎng)格節(jié)點定義向量,如圖1 所示.

圖1 顆粒定位示意圖Fig.1 Particle locating algorithm description

如果顆粒滿足下列條件,則認(rèn)為該顆粒位于圖1 所示網(wǎng)格內(nèi)

顆粒只有在計算程序運行初期時需要通過對全流場網(wǎng)格進行遍歷定位,一旦顆粒初始所在網(wǎng)格位置確定之后,在隨后的過程中,通過時間步長的限制,將顆粒的移動距離限制在相鄰網(wǎng)格中.在本文的工作中,定義了顆粒CFL 數(shù)用于限制顆粒在單個時間步長內(nèi)的位移距離

式中,dt為時間步長,Δx為顆粒所在網(wǎng)格最小尺寸,Vp為顆粒速度,CFLp為顆粒CFL 數(shù),通過顆粒CFL 數(shù)的限制,顆粒在單個時間步長內(nèi)最多只能位移至相鄰網(wǎng)格,一旦顆粒初始位置確定后,顆粒定位只需要對顆粒所在網(wǎng)格以及相鄰的網(wǎng)格進行定位即可.由于過激波后溫度梯度和速度梯度較大,在激波附近將設(shè)定更小的顆粒CFL 數(shù)以限制顆粒在單個時間步長內(nèi)的位移距離,在本文計算中,當(dāng)顆粒運動到激波所在的網(wǎng)格內(nèi)(本文根據(jù)密度梯度定義,當(dāng)顆粒所在網(wǎng)格密度梯度大于5 kg/m4時)時將其CFL數(shù)設(shè)為0.01.

2.2 流場信息的插值方法

本文考慮的火星顆粒尺寸均在微米量級,而火星探測器的直徑一般為米量級,在劃分網(wǎng)格時,塵埃顆粒的尺寸一般遠(yuǎn)小于網(wǎng)格尺寸.為了得到顆粒所在位置處的流場信息,在Euler-Lagrange 框架下,無論哪種耦合方式,都需要通過流場信息向顆粒所在位置進行準(zhǔn)確插值,插值方法有很多種,如Newton插值法和Lagrange 插值法,為了適應(yīng)不同的計算域網(wǎng)格類型,本文采用反距離加權(quán)(IDW)插值方法[22].

如圖2 所示,Z0為待求插值點,Z1～Z5為已知點信息,計算未知點到已知點的距離分別記為d1～d5.

圖2 IDW 插值模型示意圖Fig.2 Schematic of IDW interpolation model

對Z0點的插值公式如下

IDW 插值方法也適用于非結(jié)構(gòu)網(wǎng)格,一旦顆粒位置確定后,通過顆粒附近網(wǎng)格點上的信息對顆粒所在位置進行流場信息插值,從而進行顆粒的后續(xù)計算.

3 計算模型與程序驗證

3.1 計算模型及流場程序驗證

本文選用2016 年10 月進入到火星大氣的ExoMars Schiaparelli 進入艙作為計算模型,該艙防熱罩為球錐外形[23],直徑為2.4 m (Rb=1.2 m),頭部半徑Rn=0.5Rb,肩部半徑Rs=0.05Rb,其外形如圖3 所示.Schiaparelli 太空艙的任務(wù)目標(biāo)之一是在沙塵暴豐富的環(huán)境中進行火星大氣風(fēng)速、濕度、壓力、大氣顆粒特性以及火星表面溫度的測量,雖然在著陸的過程中由于導(dǎo)航數(shù)據(jù)的計算錯誤,進而使降落傘過早釋放導(dǎo)致任務(wù)失敗,但其傳送回的遙測數(shù)據(jù)詳細(xì)記錄了飛行軌跡相關(guān)信息,因此選用Schiaparelli 提供的自由來流條件作為輸入來計算流場并進行顆粒侵蝕分析.

圖3 Schiaparelli 飛行器與計算域Fig.3 Configuration of Schiaparelli and the calculation region

Gülhan 等[24]給出了Schiaparelli 飛行軌跡過程中的自由來流條件,本文選取30 km 高度的自由來流條件進行后續(xù)計算分析,即來流密度ρ=1.322 ×10-3kg/m3,速度u=2913.7 m/s,溫度T=190.1 K.

Ching 等[12]開展了20°攻角和0°攻角進入條件下的三維模擬,發(fā)現(xiàn)不同攻角造成表面侵蝕差異可以忽略不計,因此本文采用二維軸對稱熱化學(xué)非平衡數(shù)值模型,對流項采用AUSMPW+格式進行離散[25],時間推進使用LU-SGS 格式[26],黏性項使用中心差分.化學(xué)反應(yīng)模型采用Johnston 于2014 年提出的5 組分(CO2,CO,O2,O,C)二氧化碳模型,熱力學(xué)非平衡效應(yīng)采用雙溫模型.振動能松弛采用Landau-Teller 模型,CO2的振動松弛時間由Camac[27]給出,其余組分松弛時間由Millikan 等[28]的關(guān)系式給出.壁面邊界條件為非催化等溫壁面,壁面溫度Tw=300 K.同時,計算程序基于MPI 并行處理以加快求解速度.

本文首先通過與程序LAURA[29]在兩個相同工況下的計算結(jié)果對比來驗證流場計算程序,具體工況參數(shù)如表1 所示,兩狀態(tài)均為中高焓工況,駐點線存在較大占比的熱化學(xué)非平衡區(qū)域,由于熱化學(xué)非平衡會顯著影響激波脫體距離,因此通過MSL 進入器的無量綱激波脫體距離 Δ/R對程序的熱化學(xué)非平衡流動模擬準(zhǔn)度進行驗證,兩組工況下的結(jié)果如圖2所示,本文計算結(jié)果與文獻數(shù)據(jù)吻合較好.

表1 程序驗證的來流條件Table 1 Freestream conditions for program test

除激波脫體距離外,為了驗證模擬結(jié)果,本文對參數(shù)為表1 中HYPULSE749 狀態(tài)的自由來流進行壁面熱流計算,將數(shù)值結(jié)果與密歇根大學(xué)開發(fā)的LeMANS 求解器[30]以及NASA HYPULSE 膨脹管的氣動試驗數(shù)據(jù)[31]進行了對比,計算中CO2各個振動模態(tài)松弛速率統(tǒng)一采用彎曲模態(tài)速率,壁面邊界條件為300 K 等溫?zé)o滑移非催化壁面,計算外型為Mars Pathfinder 前體,其與MSL 探測器結(jié)構(gòu)類似.圖4 給出了本文程序與LeMANS 在HYPULSE749工況下壁面熱流的計算與試驗數(shù)據(jù)的對比結(jié)果,其中,試驗數(shù)據(jù)的測量不確定度包含試驗標(biāo)準(zhǔn)偏差以及試驗精度誤差,根據(jù)文獻[31],試驗數(shù)據(jù)的測量不確定度為 ±11%.結(jié)果表明,本文計算程序與LeMANS 求解器得到的壁面熱流值接近,且均在NASA HYPULSE 膨脹管的氣動試驗數(shù)據(jù)誤差范圍之內(nèi).通過對比可以認(rèn)為本文流場計算程序模擬高溫?zé)峄瘜W(xué)非平衡流動具有較好的準(zhǔn)確性,計算結(jié)果對比見表2.

表2 本文程序和LAURA 的無量綱激波脫體距離計算結(jié)果對比Table 2 Comparison of shock standoff distances of sphere-cone model calculated by different codes

圖4 壁面熱流計算結(jié)果對比Fig.4 Comparison of surface heat flux calculation

3.2 顆粒計算程序驗證

本文通過與Ching 等[10]在相同來流條件下不同初始位置顆粒在流場中的運動軌跡進行對比來驗證顆粒計算程序,計算模型為球體,其半徑為Rs=0.006 m,來流氣體組分為N2,來流馬赫數(shù)Ma∞=6.1,總壓Pt,∞=17.5 bar,總溫Tt,∞=570 K,壁面溫度Twall=300 K;顆粒密度ρd=2264 kg/m3,顆粒直徑dp=0.19 μm.兩顆粒初始位置與駐點線距離ds分別為0.1 mm 和1.8 mm;本文使用的顆粒阻力模型及Nusselt 數(shù)計算模型均與文獻[10]相同,圖5 為顆粒運動軌跡本文計算程序與文獻數(shù)據(jù)對比結(jié)果.

圖5 顆粒運動軌跡結(jié)果對比Fig.5 Comparison of particle trajectory

通過對比可以發(fā)現(xiàn)距駐點線1.8 mm 的顆粒均在到達(dá)壁面前受流場影響發(fā)生偏轉(zhuǎn),運動轉(zhuǎn)向并運動至計算域外,而初始位置ds=0.1 mm 的顆粒與壁面發(fā)生碰撞;在初始位置ds=1.8 mm 的顆粒軌跡后半段存在一定偏差,這是由于單向耦合下顆粒軌跡計算結(jié)果不僅與流場計算結(jié)果有關(guān),而且顆粒在流場中的插值方法和顆粒運動的時間步長選擇等對計算結(jié)果均有影響,而文獻[10]中并未包含完整的流場信息、選用的插值方法以及如何確定時間步長;為了對計算結(jié)果的偏差進行量化,將計算結(jié)果與文獻[10]參考數(shù)據(jù)的縱向差值與顆粒進入激波后的運動軌跡長度的比值設(shè)為偏差值Δe,通過計算得ds=1.8 mm,直徑dp=0.19 μm 的顆粒在運動軌跡上的最大偏差值Δe=0.73%,在顆粒運動軌跡規(guī)律一致的情況下認(rèn)為該偏差是可以接受的.

3.3 網(wǎng)格無關(guān)性與顆粒CFL 敏感性分析

首先,采用3 組不同的網(wǎng)格數(shù)目來對本文流場計算進行網(wǎng)格無關(guān)性說明,分別為250×180 (第一層壁面網(wǎng)格高度Δn=0.1357 mm),350×180 (Δn=0.0965 mm),450×180 (Δn=0.0744 mm),僅對垂直于壁面方向網(wǎng)格加密.針對本文選取的30 km 高度的自由來流條件,圖6 給出了3 種不同網(wǎng)格數(shù)目下駐點線上的x方向速度分布以及駐點線上的溫度分布,可以看出本文計算結(jié)果基本不受網(wǎng)格數(shù)量的影響.后續(xù)研究中采用350×180 網(wǎng)格量的結(jié)果,其中垂直于壁面方向為350 個網(wǎng)格.

圖6 網(wǎng)格無關(guān)性研究Fig.6 Grid independence studies

不同的顆粒CFL 數(shù)會改變計算的時間步長,影響到顆粒運動軌跡的計算,本節(jié)對顆粒計算程序的顆粒CFL 數(shù)進行敏感性分析,在激波附近區(qū)域存在較大的速度梯度與溫度梯度,所以當(dāng)顆粒運動到激波附近時(本文根據(jù)密度梯度定義,當(dāng)顆粒所在網(wǎng)格密度梯度大于5 kg/m4時),將顆粒CFL 數(shù)調(diào)整為0.01,而遠(yuǎn)離激波時顆粒CFL 數(shù)分別設(shè)置為0.1,0.3 和0.5.不同的顆粒CFL 數(shù)對應(yīng)的顆粒軌跡如圖7 所示,當(dāng)顆粒CFL 數(shù)在0.1～0.5 范圍內(nèi)時,計算所得顆粒的運動軌跡差異約為流場尺度的0.1%,可以認(rèn)為顆粒的運動軌跡基本重合,在本文網(wǎng)格尺度下,顆粒CFL 數(shù)取為0.3 (激波附近為0.01).

圖7 顆粒CFL 數(shù)無關(guān)性檢驗Fig.7 Particle CFL number independence test

4 結(jié)果與討論

4.1 顆粒高溫相變模型的影響

單向耦合方法下,顆粒在計算域網(wǎng)格入口邊界均勻間隔的初始位置生成并添加至穩(wěn)態(tài)的流場中,通過顆粒運動方程進行迭代以及更新顆粒的屬性,直到顆粒撞擊在壁面或離開計算域停止計算.

顆粒在流場的運動中伴隨著熱量交換,不同尺寸的顆粒到達(dá)其融化溫度需要吸收的熱量不同,圖8中給出的是dp=0.5,1.5 和2.5 μm 的顆粒在不考慮和考慮相變情況下的運動軌跡.

圖8 相變模型對顆粒軌跡影響Fig.8 Effect of phase transition model on particle trajectories

根據(jù)阻力計算式(1),可得顆粒所受到的加速度為

對于dp=0.5 μm 的顆粒,其粒徑小,慣性力小,不考慮溫升模型時大部分顆粒均受流場影響發(fā)生偏轉(zhuǎn)后,隨著流線運動至肩部以外而不與壁面發(fā)生撞擊,如圖8(a)所示;考慮溫升模型后,顆粒吸收熱量達(dá)到熔點,半徑減小,根據(jù)式(20)可知,顆粒所受加速度與其半徑成反比,因此小尺寸的顆粒加速度更大.

在多相流體動力學(xué)中,定義斯托克斯數(shù)St來衡量粒子跟隨性[32],當(dāng)斯托克斯數(shù)很小時,粒子具有足夠時間去響應(yīng)流場的變化,即跟隨性越好.St=τp/τt,τt為流動特征時間,τp為顆粒松弛時間或響應(yīng)時間[33],定義如下

式中,ρp為顆粒密度;dp為顆粒直徑;μ為流體的黏性系數(shù);Knd=λ/dp為顆?？伺瓟?shù),λ為分子平均自由程.因此dp越小對應(yīng)的St越小,對應(yīng)著顆粒的跟隨性越好,本文的計算結(jié)果也表明小尺寸顆粒展現(xiàn)出相對更好的跟隨性;考慮顆粒的高溫相變造成小尺寸顆粒體積變化量相比于其原始體積占比更重,運動一定距離后會完全融化甚至蒸發(fā)消失(本文判定依據(jù)為dp＜ 1.0×10-12m),因此dp=0.5 μm 的顆?；静粫Ρ诿媲治g產(chǎn)生影響;dp=1.5 μm 的顆粒運動軌跡與另外兩組直徑顆粒受溫升影響下軌跡偏離程度相對更大,靠近駐點線的顆粒運動時間相對較短,在其軌跡未發(fā)生明顯偏轉(zhuǎn)時已撞擊在壁面,靠近肩部的顆粒發(fā)生相變后粒徑減小,加速度相對變大,顆粒響應(yīng)時間減小,St數(shù)逐漸減小,流場跟隨性增加,且由于吸收的熱量不足以使顆粒融化至消失,其與不考慮溫升模型的顆粒軌跡相比差異最大;dp=2.5 μm 的顆粒受溫升影響相對較小,主要原因是顆粒自身粒徑大,體積大,對應(yīng)的St數(shù)相對較大,跟隨性相對較差,且在流場運動中吸收的熱量不足以達(dá)到其融化溫度,因此其運動軌跡變化相對較小.

由上分析可知,高溫相變模型對于顆粒運動的影響和其初始位置以及粒徑有關(guān);從顆粒對壁面的撞擊規(guī)律來說,相變模型對于較小或較大的顆粒影響不大,較小顆粒和流動的跟隨性較好,較大顆粒的潛熱較大,從流場吸收的熱量對其運動軌跡影響不大.本文后續(xù)的討論都是考慮有高溫相變模型的工況.

4.2 顆粒運動分析

不同粒徑的顆粒在流場中運動,其運動軌跡及其對飛行器壁面的侵蝕差異較大,圖9 選取直徑為1 μm 和5 μm 顆粒作為參考,對不同初始位置的顆粒運動軌跡進行對比分析.

圖9 不同直徑顆粒軌跡示意圖Fig.9 Trajectories of particles with different diameters

對于dp=1 μm,靠近駐點線的顆粒受流場影響較小,由圖9 可知,氣體經(jīng)過激波后靠近駐點線處的流場在y方向速度分量νg較小,靠近肩部的νg較大,而顆粒經(jīng)過激波時速度不受影響,νp為0,此時靠近駐點線處流場與顆粒y方向速度差 Δv較小,

由式(1)可知,y方向速度差 Δv越小表明顆粒受到較小的y方向流場作用力,即顆粒在y方向的加速度越小;此外,在駐點線附近激波脫體距離比較小,氣體對顆粒作用時間短.因此,靠近駐點線的顆粒更容易與壁面發(fā)生撞擊;而靠近肩部的顆粒與流場y方向速度差Δν相對較大,作用時間也相對較長,受流場影響較大,更貼近流線運動;對于dp=5 μm的顆粒,其具有較大的質(zhì)量和慣性,在流場中受影響較小,基本都以x方向運動為主,最后撞擊到壁面.

圖10 為在同一初始位置下對不同粒徑顆粒軌跡的模擬.根據(jù)式(20)計算,ap與rp成反比關(guān)系,隨著粒徑的增大,顆粒受到的加速度變小,即流場對顆粒作用越晚體現(xiàn),dp=0.5 μm 的顆粒在進入激波后最早在流場作用下出現(xiàn)“轉(zhuǎn)向”現(xiàn)象,且在運動過程中由于吸熱完全蒸發(fā);其次是dp=1.0 μm 的顆粒,在流場作用下最終脫離計算域;而對1.5,2.0 和2.5 μm的顆粒流場作用較小,均穿過流場后撞擊壁面,顆粒越小,向上偏離的位置越遠(yuǎn).

圖10 不同粒徑軌跡示意圖Fig.10 Trajectory of different particle sizes

在飛行器進入的過程中,需要關(guān)注的是塵埃顆粒對TPS 的撞擊侵蝕,這里引入Connolly 等[34]提出的撞擊分?jǐn)?shù)η以及撞擊能量分?jǐn)?shù)χ的概念;顆粒撞擊分?jǐn)?shù)η被定義為撞擊防熱罩的顆粒數(shù)與基于防熱罩半徑Anose的截面積上的顆?？倲?shù)的比值.受影響的顆粒上游有一個固定的截面積Aimpact,可以利用截面積計算出各尺寸顆粒的撞擊分?jǐn)?shù)η

其中,rimpact是能夠發(fā)生撞擊的顆粒在上游徑向的最大半徑,rnose是防熱罩半徑,對應(yīng)圖3 中的Rb.

圖11 為不同尺寸直徑顆粒分別在考慮和不考慮相變模型下進行計算得到的撞擊分?jǐn)?shù)對比,由圖11可知,無論是否考慮相變模型,直徑1 μm 以下的顆粒撞擊分?jǐn)?shù)均低于0.1,主要原因是小尺寸顆粒的顆粒響應(yīng)時間較小,在到達(dá)壁面之前大部分均被流場“捕獲”運動至計算域以外,另外由于自身體積較小,相變造成顆粒體積變化量相比于其原始體積占比更重,甚至達(dá)到融化或完全蒸發(fā)的情況,因此并未對TPS 造成侵蝕;而顆粒直徑達(dá)到3 μm 以上時,其對應(yīng)的顆粒響應(yīng)時間比較大,顆粒沒有足夠的時間響應(yīng)流場的變化,且由于顆粒自身的體積較大,高溫相變模型對其半徑的影響較小,因此對應(yīng)較高的撞擊分?jǐn)?shù).

圖11 相變模型對撞擊分?jǐn)?shù)影響Fig.11 Effect of phase transition model on impact fraction

4.3 顆粒撞擊能量研究

對于本文關(guān)注的表面侵蝕來說,由于目前沒有完善的微尺寸顆粒侵蝕模型,一般認(rèn)為顆粒對壁面的侵蝕主要與其對壁面造成的撞擊能量有關(guān),而顆粒造成的撞擊能量主要與它的法向動能有關(guān),其法向動能的損失等于對壁面造成的撞擊能量,即

其中,KE⊥是法向動能,mp是顆粒的質(zhì)量,ν⊥是顆粒撞擊壁面時的法向速度.

n個直徑相同的顆粒造成的撞擊能量可以寫為

Palmer 對火星大氣塵埃顆粒的大小和數(shù)量密度進行了建模,本文選取其提出的模態(tài)半徑rm=0.35 μm對應(yīng)的顆粒分布模型,其質(zhì)量分?jǐn)?shù)分布如圖12 所示[9].

圖12 顆粒分布模型[9](rm=0.35 μm)Fig.12 Particle distribution model for rm=0.35 μm

所有尺寸的顆粒造成的總撞擊能量為

其中,md是直徑為d的顆粒質(zhì)量,將某特定尺寸直徑d的顆粒產(chǎn)生的能量對總撞擊能量的占比定義為撞擊能量分?jǐn)?shù)χd

圖13 為基于本文采用的顆粒分布模型進行模擬得到的撞擊能量分?jǐn)?shù)分布.由圖可知,平均直徑為5 μm 的顆粒造成的撞擊能量在該顆粒分布模型所造成的總撞擊能量中占比最高,約占13%;對撞擊能量主要貢獻的顆粒尺寸集中在3～10 μm,約占總撞擊能量的80%;由圖13 可知,小于3 μm 的顆粒對應(yīng)相對較小的撞擊分?jǐn)?shù),且由于自身質(zhì)量較小,撞擊在壁面產(chǎn)生的撞擊能量較小;直徑10 μm 以上的顆粒雖然具有較大的撞擊分?jǐn)?shù)以及單個顆粒造成較高的撞擊能量,但由于其較小的質(zhì)量分?jǐn)?shù),對應(yīng)較小的凈質(zhì)量濃度,因此直徑小于3 μm 或在10 μm 以上的顆粒在撞擊能量中占比并不高.

圖13 撞擊能量分?jǐn)?shù)分布Fig.13 Impact energy fraction distribution

5 結(jié)論

本文通過數(shù)值模擬方法研究了ExoMars Schiaparelli太空艙在火星大氣30 km 高度的來流條件下,大氣中微塵顆粒的運動軌跡;同時研究了在給定含塵顆粒質(zhì)量分?jǐn)?shù)模型下,不同尺寸顆粒對應(yīng)的撞擊分?jǐn)?shù)以及撞擊能量分?jǐn)?shù).本文得出結(jié)論如下.

(1) 引入高溫相變模型對不同尺寸的顆粒有不同的影響,直徑小于0.5 μm 的顆粒由于其體積小,易達(dá)到融化溫度直至完全蒸發(fā)消失;對于直徑為1.5 μm的顆粒,顆粒在吸熱后發(fā)生相變粒徑減小,高溫相變模型對其運動軌跡影響較大;直徑大于 2.5 μm 的顆粒由于體積較大,吸收的熱量不足以發(fā)生相變,高溫相變模型對運動軌跡影響非常小.

(2) 直徑1 μm 以下的顆粒跟隨性相對較強,對應(yīng)的撞擊分?jǐn)?shù)接近于0,即大部分顆粒不會與模型壁面發(fā)生碰撞;顆粒直徑3 μm 以上時,其撞擊分?jǐn)?shù)達(dá)95%以上,即幾乎所有初始位置的顆粒全部會撞擊到模型壁面.

(3) 考慮顆粒的撞擊能量主要與其法向動能有關(guān),直徑在3～10 μm 之間的顆粒是壁面侵蝕的主要來源,其造成的撞擊能量約占總撞擊能量的80%.