何青松, 王廣興, 王奇, 張 章, 王立武, 賈賀

(1.北京空間機(jī)電研究所, 北京 100094;2.中國航天科技集團(tuán)有限公司 航天進(jìn)入、減速與著陸技術(shù)實(shí)驗(yàn)室, 北京 100094)

充氣式返回艙是一種新型的氣動減速裝置,再入大氣時利用充氣形成氣動外形,并由表面耐高溫柔性材料提供熱防護(hù),最后由自身的充氣結(jié)構(gòu)實(shí)現(xiàn)著陸緩沖。充氣式返回艙集合了傳統(tǒng)飛行器的熱防護(hù)系統(tǒng)、降落傘減速系統(tǒng)以及著陸/著水緩沖系統(tǒng),簡化了航天器回收系統(tǒng)的設(shè)計(jì)[1]。充氣式返回艙拓展了傳統(tǒng)航天器回收著陸技術(shù),具有廣泛的應(yīng)用前景,各航天大國均開展了廣泛的研究[2-8],我國也在2020年5月用長征五號B運(yùn)載火箭進(jìn)行了首次柔性充氣式貨物返回艙搭載試驗(yàn)。

不同于傳統(tǒng)的剛性結(jié)構(gòu)返回艙,充氣式返回艙是由柔性材料制成,在高速氣流的沖擊下,柔性結(jié)構(gòu)會發(fā)生變形甚至破壞。因此,充氣式返回艙的結(jié)構(gòu)和防熱的設(shè)計(jì)十分重要,而獲得準(zhǔn)確的流場特性參數(shù)是進(jìn)行充氣式返回艙結(jié)構(gòu)和防熱設(shè)計(jì)的重要前提。相比于實(shí)驗(yàn)研究,數(shù)值模擬能夠快速獲得充氣式返回艙的流場特性參數(shù),縮短研制周期,降低研制成本。為獲得準(zhǔn)確可靠的充氣式返回艙流場特性參數(shù),建立合理的數(shù)值仿真模型,選擇合理的物性參數(shù)十分重要。目前充氣式返回艙大氣再入氣動特性仿真大多基于理想氣體假設(shè),黏性系數(shù)和熱導(dǎo)率等物性參數(shù)采用Sutherland公式[9-13]。通常情況,常溫下空氣的黏性和熱導(dǎo)率可采用Sutherland公式計(jì)算,但當(dāng)溫度進(jìn)一步升高時,空氣中的氧氣和氮?dú)鈱l(fā)生解離反應(yīng),更高溫度時將發(fā)生電離反應(yīng),這一系列的化學(xué)反應(yīng)使空氣物性發(fā)生復(fù)雜變化,Sutherland公式將不再適用[14]。充氣式返回艙高速再入時流場溫度非常高[15],若采用Sutherland公式計(jì)算空氣物性參數(shù)會帶來巨大誤差。因此,為了獲得準(zhǔn)確的氣動特性參數(shù),首先需要獲得高溫下空氣的物性參數(shù)。

為解決高溫條件下空氣物性參數(shù)的計(jì)算問題,學(xué)者們已經(jīng)進(jìn)行了深入研究[16-19],目前應(yīng)用較為廣泛的是Gupta的高溫空氣擬合關(guān)系式[17],該擬合關(guān)系式已經(jīng)被廣泛應(yīng)用于超聲速等離子體風(fēng)洞、電弧加熱器、圓柱繞流以及平板邊界層流動等領(lǐng)域的研究[14,20-22]。

本文分別采用Gupta的高溫空氣擬合關(guān)系式和Sutherland公式計(jì)算空氣的物性參數(shù),對充氣式返回艙高速再入流場進(jìn)行數(shù)值模擬研究,比較采用這兩種空氣物性參數(shù)計(jì)算方法獲得的流場特性參數(shù)差異,進(jìn)而為充氣式返回艙的工程設(shè)計(jì)提供一定的理論參考。

1 計(jì)算模型

本文采用二維數(shù)值模擬對充氣式返回艙氣動特性進(jìn)行研究,求解可壓縮N-S方程,在笛卡爾坐標(biāo)系下控制方程可以寫成如下形式

式中:ρ為流體密度；xi為笛卡爾坐標(biāo);ui為速度分量;p為壓強(qiáng);E為內(nèi)能;H為焓值;τij和qi分別為黏性應(yīng)力項(xiàng)和熱通量項(xiàng)。

湍流模型采用k-ω模型,模型輸運(yùn)方程如(4)～(5)式所示

式中,β′=0.09,α=1.4,β=0.075,σk=2,σω=2。

本文研究的流體介質(zhì)為空氣,分別采用Sutherland公式和Gupta的擬合公式對方程中涉及到的物性參數(shù)進(jìn)行計(jì)算。公式(6)和(7)為Sutherland公式表示的黏性系數(shù)和熱導(dǎo)率。Gupta擬合公式基于空氣11組分模型(O2,N2,O,N,NO,O+,N+,NO+,O2+,N2+,e)計(jì)算高溫空氣物性,考慮了空氣因高溫導(dǎo)致的空氣化學(xué)反應(yīng),溫度范圍最高達(dá)30 000 K[19]。公式(8)和(9)為Gupta擬合公式的黏性系數(shù)和熱導(dǎo)率。

(6)

(7)

(8)

(9)

式中:T,μ,κ分別為溫度、黏性系數(shù)和熱導(dǎo)率;μ0=1.716 1×10-5Pa·s,κ0=2.415×10-2W/(m·K);To=273.16 K,Tc=110.4 K,T0為攝氏零度,Tc為參考溫度,χμ,χκ,Aμ,Aκ等參數(shù)可參考文獻(xiàn)[19]。

本文以日本研制的充氣式返回艙實(shí)驗(yàn)?zāi)Ｐ虷WT-MAAC為研究對象,該模型為飛行演示驗(yàn)證試驗(yàn)中采用的充氣式返回艙SMAAC的縮比模型,其結(jié)構(gòu)示意如圖1所示,尺寸參數(shù)可參見文獻(xiàn)[23]。

圖1 HWT-MAAC模型結(jié)構(gòu)示意圖

為考察不同溫度范圍2種物性參數(shù)計(jì)算方法帶來的結(jié)果差異,本文在文獻(xiàn)[23]實(shí)驗(yàn)工況的基礎(chǔ)上,通過改變來流溫度和馬赫數(shù)來進(jìn)行計(jì)算研究。實(shí)驗(yàn)工況為來流溫度45 K,來流馬赫數(shù)10,來流壓強(qiáng)均為60 Pa,表1列出了本文的5個計(jì)算工況,所有計(jì)算工況的來流攻角均為0°。

表1 計(jì)算工況

圖2給出了本文的計(jì)算域,計(jì)算域?yàn)榘雸A形,直徑約為返回艙直徑的10倍,返回艙前端距計(jì)算域前端約為2倍返回艙的直徑。本文的計(jì)算工況為超聲速流動,流場的擾動不會向上游傳播,計(jì)算域能夠捕捉完整的激波結(jié)構(gòu),則計(jì)算域的選取對計(jì)算結(jié)果的影響較小。計(jì)算域入口給定來流馬赫數(shù)、來流壓強(qiáng)以及來流溫度,出口參數(shù)外推,壁面采用無滑移邊界條件,壁面溫度為300 K。

圖2 計(jì)算域以及網(wǎng)格劃分情況

2 計(jì)算結(jié)果

首先本文進(jìn)行了網(wǎng)格無關(guān)性驗(yàn)證,分別采用網(wǎng)格數(shù)為30 000,60 000以及120 000的網(wǎng)格數(shù)量對工況3進(jìn)行了計(jì)算,同時對近壁網(wǎng)格進(jìn)行加密保證y+小于1,比較了沿壁面的熱流密度,比較結(jié)果如圖3所示。從圖中可以看出網(wǎng)格數(shù)量為60 000和120 000的計(jì)算結(jié)果十分接近,而網(wǎng)格數(shù)量為30 000的計(jì)算結(jié)果與另外2個相差較遠(yuǎn)。因此,為了能夠保證計(jì)算的準(zhǔn)確性,同時也能節(jié)省計(jì)算資源,本文采用網(wǎng)格數(shù)量為60 000,圖4給出了網(wǎng)格劃分情況。

圖3 網(wǎng)格無關(guān)性驗(yàn)證

圖4 網(wǎng)格劃分(網(wǎng)格量60 000)

為了驗(yàn)證模型的可靠性,本文將工況3計(jì)算得到的充氣式返回艙迎風(fēng)面的熱流密度與文獻(xiàn)[23]的實(shí)驗(yàn)結(jié)果和三維仿真結(jié)果進(jìn)行了比較,如圖5所示。實(shí)驗(yàn)中采用高超聲速風(fēng)洞對充氣式返回艙實(shí)驗(yàn)?zāi)Ｐ虷WT-MAAC進(jìn)行了研究,采用輻射平衡假設(shè)對表面熱流密度進(jìn)行測量。從圖中看出本文采用Gupta擬合公式和Sutherland公式的計(jì)算結(jié)果與實(shí)驗(yàn)值均符合較好,在阻力面處,本文的計(jì)算結(jié)果比文獻(xiàn)[23]的仿真結(jié)果略高。

圖5 返回艙壁面熱流密度與文獻(xiàn)[23]結(jié)果比較

為進(jìn)一步驗(yàn)證本文模型的可靠性,本文選擇文獻(xiàn)[24]的實(shí)驗(yàn)結(jié)果進(jìn)行驗(yàn)證,實(shí)驗(yàn)?zāi)Ｐ蜑?0 mm圓柱體,采用高焓激波風(fēng)洞對圓柱體迎風(fēng)面的熱流密度和壓強(qiáng)進(jìn)行了測量,來流馬赫數(shù)為8.78,來流溫度694 K,來流壓強(qiáng)為687 Pa,比較結(jié)果如圖6所示。從圖中看出采用Gupta擬合公式獲得的表面熱流密度與實(shí)驗(yàn)結(jié)果符合得更好,采用這兩種物性計(jì)算方法得到的表面壓強(qiáng)與實(shí)驗(yàn)結(jié)果符合均較好。

圖6 圓柱體計(jì)算結(jié)果與文獻(xiàn)[24]結(jié)果比較

溫度對空氣物性參數(shù)有十分重要的影響,圖7給出了本文計(jì)算工況的溫度分布情況,從圖中看出采用Gupta擬合公式計(jì)算得到的激波位置更靠近返回艙表面,流場的溫度也顯著比Sutherland的計(jì)算結(jié)果低。另外,從圖中還看出這種差異隨著來流馬赫數(shù)和溫度的增加更加明顯,比如工況3的來流馬赫數(shù)和溫度相對來說較低,Gupta擬合公式和Sutherland公式計(jì)算結(jié)果相對來說也較為接近。工況3中流場最高溫度不超過1 000 K,Sutherland公式還能較好地描述空氣的物性,Sutherland公式和Gupta擬合公式計(jì)算得到的空氣物性差別不大。而隨著來流馬赫數(shù)和溫度的增加,流場中的溫度逐漸上升,Gupta擬合公式和Sutherland公式計(jì)算得到的空氣物性參數(shù)差別較為明顯,兩者計(jì)算得到的流場特性差別逐漸顯現(xiàn)。

圖7 流場中溫度分布情況

對于高超聲速激波形狀和距離的研究始于20世紀(jì)50年代,涉及理論和實(shí)驗(yàn)研究[25]。以大量實(shí)驗(yàn)數(shù)據(jù)為基礎(chǔ)得出的激波脫體距離工程算法計(jì)算量小,計(jì)算速度快,Inouye[26]建立了多種氣體環(huán)境下的繞半球體脫體激波距離關(guān)系式如(10)式所示

(10)

式中：Δ為激波脫體距離；R為半球半徑；ρ∞和ρ分別為波前和波后氣體的密度。從這個關(guān)系式中可以發(fā)現(xiàn),脫體激波的距離跟波前和波后氣體密度的比值成正比。為探究采用2種物性參數(shù)方法計(jì)算的激波位置產(chǎn)生差異的原因,圖8給出了沿流場中心線氣體的密度。從圖8中看出在波前采用2種物性參數(shù)計(jì)算方法得到的密度相差不大,而在波后Gupta擬合公式計(jì)算得到的氣體密度大于Sutherland公式計(jì)算得到的氣體密度,即Gupta擬合公式獲得的波前和波后的氣體密度的比值小于Sutherland公式,因此脫體激波的距離更接近物面。

圖8 沿中心線氣體的密度

不管是來流馬赫數(shù)增加還是來流溫度的增加,流場溫度的上升是導(dǎo)致兩者計(jì)算結(jié)果差異的重要因素,表2列出了5個工況流場中的最高溫度,TG和TS分別表示采用Gupta擬合公式和Sutherland公式計(jì)算得到的最高溫度,比較發(fā)現(xiàn)隨著來流馬赫數(shù)和溫度的增加,流場中最高溫度的差距增加十分明顯。產(chǎn)生流場溫度差異的主要原因是Gupta擬合公式考慮了空氣高溫下因振動、解離以及電離等化學(xué)反應(yīng)對物性參數(shù)的影響,空氣的這些化學(xué)反應(yīng)能夠吸收返回艙外部的部分能量,導(dǎo)致返回艙外部溫度降低。

表2 流場中最高溫度比較

在高速氣流沖擊下,充氣式返回艙迎風(fēng)面將承受劇烈的氣動力和氣動熱,準(zhǔn)確的氣動熱和氣動力特性對充氣式返回艙的結(jié)構(gòu)和防熱設(shè)計(jì)十分重要。圖9給出了返回艙迎風(fēng)面的熱流密度分布,從圖中看出熱流密度在頭錐和兩側(cè)圓環(huán)處均存在峰值,最大峰值在頭錐處。從圖中還可以看出返回艙表面的熱流密度隨著來流馬赫數(shù)和溫度增加上升十分顯著,以Gupta擬合公式計(jì)算結(jié)果為例,來流Ma=10,T=45 K,頭錐處峰值的熱流密度為79 kW/m2,來流Ma=10,T=135 K, 頭錐處峰值的熱流密度為451 kW/m2,來流Ma=18,T=45 K,頭錐處峰值的熱流密度為1 022 kW/m2。比較發(fā)現(xiàn),Sutherland公式計(jì)算的熱流密度總體上低于Gupta擬合公式的結(jié)果,Sutherland公式低估了返回艙表面的熱流密度,如工況5所示,來流Ma=18,T=45 K,Gupta擬合公式獲得的熱流密度峰值高出Sutherland公式獲得的熱流密度峰值29%。其中工況3由于兩者計(jì)算得到的流場溫度相差不大,兩者物性參數(shù)差別不大,Gupta公式的計(jì)算結(jié)果僅略微高于Sutherland公式的結(jié)果。

圖9 返回艙迎風(fēng)面的熱流密度分布 圖10 返回艙迎風(fēng)面的壓強(qiáng)分布

由于返回艙表面的氣動力主要來源于表面壓強(qiáng)的作用,圖10給出了返回艙迎風(fēng)面處的壓強(qiáng)分布情況,從圖中看出壓強(qiáng)的分布規(guī)律與熱流密度的分布規(guī)律類似,在頭錐和兩側(cè)圓環(huán)處均存在峰值,最大峰值在頭錐處。從圖10a)中看出,來流溫度對迎風(fēng)面的壓強(qiáng)影響并不大,迎風(fēng)面壓強(qiáng)隨來流溫度增加僅略微增加。從圖10b)中看出,來流馬赫數(shù)對迎風(fēng)面壓強(qiáng)影響較大,迎風(fēng)面的壓強(qiáng)隨來流馬赫數(shù)的增加顯著上升。比較發(fā)現(xiàn),Gupta擬合公式計(jì)算得到的迎風(fēng)面壓強(qiáng)大于Sutherland公式的計(jì)算結(jié)果,但相對于表面熱流密度的差異,表面壓強(qiáng)的差距并不是很大。

充氣式返回艙的柔性熱防護(hù)系統(tǒng)通常由防熱層、隔熱層和阻力承力層構(gòu)成,最外層的防熱層主要用來阻擋熱流密內(nèi)部氣囊滲漏,保持充氣結(jié)構(gòu)的形狀[27]。準(zhǔn)確的氣動力和氣動熱特性是各層材料選擇以及材料鋪層數(shù)量選取的重要參考,Sutherland公式在高溫下不再適用,低估了充氣式返回艙表面的熱流密度和壓力,則有可能造成設(shè)計(jì)的充氣結(jié)構(gòu)防熱性能和強(qiáng)度不足,因此在高溫環(huán)境下采用Gupta擬合公式進(jìn)行研究更為合理。

熱導(dǎo)率是引起熱流密度差異的重要原因,以工況5為例,圖11給出了流場中的熱導(dǎo)率分布,從圖中看出Gupta擬合公式獲得的熱導(dǎo)率顯著高于Sutherland公式獲得熱導(dǎo)率,流場中最大熱導(dǎo)率分別為0.152和0.107 W/(m·K)。高溫下,空氣的熱導(dǎo)率由平動熱導(dǎo)率和反應(yīng)熱導(dǎo)率組成,Gupta擬合公式考慮了空氣的反應(yīng)熱導(dǎo)[19],因此采用Gupta擬合公式獲得的熱導(dǎo)率更大。由圖7可知,Sutherland公式計(jì)算得到的流場溫度更高,在近壁區(qū)域的溫度梯度也更高,但Sutherland獲得的近壁區(qū)域熱導(dǎo)率比Gupta擬合公式低,因此導(dǎo)致壁面的熱流密度更低。

圖11 流場中熱導(dǎo)率分布情況

3 結(jié) 論

本文分別采用Gupta擬合公式和Sutherland公式計(jì)算空氣物性參數(shù),比較了采用這兩種計(jì)算公式獲得的充氣式返回艙流場特性的差異,研究表明:①采用Gupta擬合公式獲得的激波位置比Sutherland公式的結(jié)果更靠近返回艙,這種差異在高來流馬赫數(shù)和溫度時更加明顯;②采用Gupta擬合公式獲得的流場溫度更低,且隨著來流馬赫數(shù)和溫度的增加,流場最高溫度的差異更加明顯;③Sutherland公式低估了返回艙迎風(fēng)面的熱流密度,尤其在來流溫度和馬赫數(shù)較高的工況,Gupta擬合公式獲得的熱流密度顯著高于Sutherland公式的結(jié)果, 工況5中Gupta擬合公式獲得的熱流密度峰值高出Sutherland公式獲得的熱流密度峰值29%;④Sutherland公式獲得的迎風(fēng)面壓強(qiáng)低于Gupta擬合公式的結(jié)果,但是壓強(qiáng)的差異相對于熱流密度的差異小很多。