基于DSMC 方法的再入飛行器微燒蝕研究1)

甘 馳 陳 松,2) 張 俊

* (北京航空航天大學中法工程師學院/國際通用工程學院,北京 100191)

? (北京航空航天大學航空科學與工程學院,北京 100191)

引言

飛行器在以高超聲速再入大氣層時將會產(chǎn)生一道強激波,使得飛行器壁面周圍溫度急劇升高,導致空氣發(fā)生振動激發(fā)、離解甚至電離等一系列物理化學變化[1].針對這種極端的氣動熱環(huán)境,常常需要進行熱防護系統(tǒng)(TPS)設計.相比于可反復使用的熱防護系統(tǒng)而言,高超聲速條件下更需要燒蝕性的熱防護系統(tǒng).這種燒蝕性TPS 一般是多孔設計,通過材料質量的損失例如熱解、氧化和升華等一系列復雜的物理化學反應來為飛行器緩解極端的熱載荷.隔熱材料熱解后產(chǎn)生的氣體會沿著材料向外移動到邊界層區(qū)域,由溫度相對較低的熱解氣體帶走部分熱量.同時由于TPS 材料的多孔設計,高溫的邊界層氣體也可能會向內(nèi)流入材料微觀結構中,引起更深層次的化學反應,影響結構內(nèi)部的熱流.因此研究局部區(qū)域的微燒蝕問題對于整個熱防護系統(tǒng)的設計顯得十分重要.

當飛行器因強烈氣動加熱而發(fā)生燒蝕時,燒蝕外形的變化會進一步影響流場情況,因此在模擬燒蝕過程時需要進行耦合計算[2].為此,除了進行氣動熱環(huán)境分析外,還需要建立準確的燒蝕模型,以及如何在燒蝕過程中高效地處理燒蝕移動邊界的后退問題.Liang 等[3]針對“天宮一號”的再入過程建立了一維燒蝕模型和傳熱模型;周印佳等[4]采用Sutton-Graves 和Tauber-Sutton 理論計算駐點熱流以及Landau 坐標系實現(xiàn)燒蝕動網(wǎng)格,通過求解瞬態(tài)有限差分熱傳導方程實現(xiàn)了燒蝕熱防護一體化計算方法;李偉等[5]采用有限體積法建立了碳纖維增強復合材料的微觀燒蝕數(shù)值模型,并引入了固體相體積分數(shù)與界面重構技術;胥建宇[6]提出了基于比色測溫的燒蝕模型表面溫度測量方法,并建立了基于LSTM神經(jīng)網(wǎng)絡的燒蝕量預測模型;王曉婕[7]組合使用生死單元法、網(wǎng)格自適應技術和網(wǎng)格重劃分技術建立了氣動熱環(huán)境-固體導熱-燒蝕一體化耦合計算方法;楊凱威等[8]針對有限元軟件采用了自編熱流加載和燒蝕移動邊界用戶子程序方法,建立了高速空氣舵前緣三維燒蝕和溫度場耦合計算方法.

然而,在相當大的高度時,飛行器所處的氣體環(huán)境十分稀薄,克努森數(shù)Kn＞ 0.01,這將使得傳統(tǒng)的計算流體方法難以給出準確的模擬結果.在高度不斷下降的過程中,雖然飛行器外部周圍的流體可能整體上是連續(xù)的,但由于材料尺度較小,并且飛行器頭部存在各種端口、緊固件等微小粗糙元,局部的克努森數(shù)可能仍處于過渡區(qū)甚至稀薄區(qū).因此,在這種情況下需要采用適用于稀薄和近連續(xù)領域的具有高保真度的氣體動力學方法.

直接模擬蒙特卡洛方法(DSMC)于20 世紀60年代由學者Bird 提出[9],從統(tǒng)計學的角度求解玻爾茲曼方程,通過模擬大量微觀粒子的運動和碰撞來反映真實的流場情況.DSMC 方法適用于模擬稀薄氣體流動,并且發(fā)展至今已能夠應用于許多領域[10].現(xiàn)有的眾多程序都能夠很好地實現(xiàn)DSMC 方法,如Bird 開發(fā)的代碼DS2V[11]和DS3V[12],美國國家航空航天局NASA 開發(fā)的代碼DAC[13-14],Dietrich 等[15]開發(fā)的代碼MONACO,Scamlon 等[16-17]開發(fā)的并行代碼dsmcFoam 等.

近年來,Plimpton 等[18]基于DSMC 方法開發(fā)出了開源并行程序SPARTA.SPARTA 采用多層笛卡爾網(wǎng)格對粒子進行分組,能夠高效地模擬粒子間的碰撞以及化學反應,同時允許在流體網(wǎng)格中加入固體壁面對象,并模擬粒子與壁面之間的碰撞與反射過程.程序采用了Marching Square 算法來生成以及實時更新燒蝕外形,并記錄粒子與壁面碰撞時傳遞給壁面的能量,模擬燒蝕面吸收的熱量來扣除網(wǎng)格的節(jié)點值,進而模擬燒蝕過程.該程序能夠同時考慮多組分化學反應,具備模擬高超聲速化學非平衡流動的能力.流場初始條件以及網(wǎng)格信息通過讀取外部輸入文件來設置.

隨著統(tǒng)一隨機粒子(USP)方法的提出,SPARTA還在不斷發(fā)展之中[19].針對于飛行器頭部局部粗糙區(qū)域的燒蝕問題,本文采用最新的標準SPARTA 源碼,并在此基礎上進行了改進.首先對粒子與壁面碰撞時動能的記錄規(guī)則進行了修改,避免了負能量出現(xiàn)的情況;其次添加了縮放系數(shù)使網(wǎng)格節(jié)點值與粒子動能值處于同一量級,減少一個燒蝕循環(huán)內(nèi)的溢出損失;最后對程序中燒蝕速率的計算方法進行了改進,以更準確地描述燒蝕物理過程.本文模擬了幾種典型幾何外形的燒蝕過程,并給出了最終的燒蝕外形以及流場圖像.在不考慮轉捩現(xiàn)象[20]以及氣體離解反應和電離反應的情況下,球錐燒蝕的模擬結果與參考文獻[21]的燒蝕結果在駐點線附近取得了較好的一致性,驗證了DSMC 方法以及SPARTA 程序對于燒蝕問題研究的可行性.同時,通過斜楔體的燒蝕計算本文對帶有微小粗糙元的飛行器表面的燒蝕問題進行了探討.如何將更加復雜的因素加入到程序中從而給出更加準確的燒蝕計算結果也是今后主要的研究方向.

1 研究方法

1.1 分子動理論與DSMC 方法

分子動理論是物理學中研究物質微觀粒子運動行為的理論.該理論認為分子的運動是隨機的,并遵循牛頓力學和統(tǒng)計物理學的原理.所有微觀分子的運動和相互作用共同影響了氣體的宏觀性質,例如密度、溫度和壓強等性質.當氣體密度很低時,分子平均自由程遠大于分子自身的尺寸,連續(xù)介質假設失效,此時的氣體被稱之為稀薄氣體.在稀薄氣體領域,分子間的相互作用可以忽略不計,分子個體的行為所產(chǎn)生的影響更加顯著.在這種情況下,分子被認為是非常小的硬球,其主要行為包括平動和碰撞.分子在平均自由程內(nèi)可視為作直線運動.

稀薄氣體動力學的核心是玻爾茲曼方程.該方程描述了分子在流場中的運動和碰撞.然而由于玻爾茲曼方程形式復雜,直接求解相當困難.因此,許多數(shù)值模擬方法被發(fā)展出來用以求解稀薄氣體問題.其中最常用的方法是DSMC 方法.

DSMC 方法是基于統(tǒng)計物理學的數(shù)值模擬方法,最初被用來模擬稀薄氣體的動力學行為,經(jīng)過數(shù)十年的發(fā)展已經(jīng)在燒蝕問題以及與流場的耦合計算等研究中具有重要作用.在該方法中,每個模擬粒子代表了大量相同的真實氣體分子.粒子的運動和碰撞在一個時間步內(nèi)可認為是解耦的.粒子在該步長下沿直線運動,并被隨機選取來模擬分子間的碰撞.常用的分子碰撞模型有VHS 模型、VSS 模型.此外,碰撞過程中的化學反應由熱化學模型來模擬,例如TCE 模型和QK 模型等.在燒蝕問題中,粒子會與固體表面發(fā)生碰撞并反射或被吸收從而反映氣-固間的相互作用過程.

1.2 SPARTA 程序

SPARTA 是基于DSMC 方法的開源并行程序,其運行效率高,用戶可以方便地添加新功能.該程序中的對象可分為3 類: 模擬粒子(simulated particles)、計算網(wǎng)格(grid cells) 以及流體中的壁面(surface elements).計算網(wǎng)格始終是結構化網(wǎng)格,可用于二維和三維流動的模擬.在二維流動中,網(wǎng)格為與兩條軸相互平行的矩形,壁面外形由數(shù)個點兩兩連接而成的線段組成;而在三維流動中,網(wǎng)格則為長方體形狀,壁面由眾多點三三連接形成的三角形組成.此外,網(wǎng)格可以根據(jù)需要進行多層加密,以及在不同區(qū)域進行不同程度的加密.

SPARTA 的輸入文件包括了許多計算參數(shù)的設置,例如計算域大小、網(wǎng)格尺寸、來流條件、采樣設置以及需要計算的流體性質等.在流場中,粒子的性質也需要進行設置,如粒子的質量、轉動自由度、振動自由度等.SPARTA 使用VHS 或VSS 等碰撞模型來模擬粒子之間的相互作用.用于燒蝕的隱式壁面(implicit surface)信息由二進制文件存儲.

SPARTA 工作流程示意圖如圖1 所示.在開始計算過程之后,程序首先會將輸入文件中的命令全部讀入,之后再開始執(zhí)行相應的命令.程序會根據(jù)輸入的網(wǎng)格信息自動生成結構化網(wǎng)格并初始化.接著,程序會在指定的邊界處生成模擬粒子并為其賦予初始速度等信息,然后按照該速度來移動粒子.在每個時間步內(nèi),計算域外的粒子會被移除,而與壁面碰撞的粒子會發(fā)生反射或被吸收.此外,程序會隨機選取粒子對并模擬碰撞過程,碰撞后粒子的速度會發(fā)生改變.循環(huán)內(nèi)程序會每隔一定時間步進行采樣并在循環(huán)的最后輸出計算結果,如包括粒子和壁面在內(nèi)的流場圖像、每個網(wǎng)格的坐標以及網(wǎng)格內(nèi)流場的溫度、壓強、密度等信息.在輸入文件中會規(guī)定總共模擬的步數(shù)Ntotal,如果當前模擬步數(shù)N達到了總共的模擬步數(shù)Ntotal,程序就會停止運算.

圖1 SPARTA 基本工作流程圖Fig.1 The basic flowchart of SPARTA

1.3 燒蝕模型

目前,計算燒蝕速率的方法主要包括解析燒蝕外形方程[2]、求解燒蝕邊界熱平衡方程[3,22-23]、采用氣動熱化學燒蝕方法[24-26]等.SPARTA 采用類似于熱平衡邊界假設的燒蝕原理,認為氣動加熱量與燒蝕吸熱量達到平衡.

SPARTA 采用Marching Square 算法對燒蝕后的壁面外形進行實時更新.該算法與Marching cubes 方法類似,能夠根據(jù)輸入的二維標量場生成云圖或等值線.該算法一般用于作地形圖上的等高線或者氣象圖中的等壓線,而在SPARTA 中用于生成燒蝕外形.算法的輸入可以看作是一個二維矩陣,矩陣中的每一個元素是一個節(jié)點,4 個節(jié)點可以連接形成一個矩形,相當于一個計算網(wǎng)格.因此每個網(wǎng)格都具有4 個節(jié)點值.通常情況下,該算法會被給定一個閾值,大于該閾值的節(jié)點值被視為1,小于則視為0,因此每個網(wǎng)格的4 個節(jié)點分別可能為0 或者1,總共有16 種不同的情況.對應這16 種情況,算法給定一張查找表,不同的情況下網(wǎng)格會有不同的線段幾何外形,并且每個網(wǎng)格一定能在該索引表中尋找到相應的線段幾何外形,如圖2 所示.

圖2 Marching Square 算法的填充規(guī)則Fig.2 Look up table of Marching Square algorithm

Marching Square 算法的基本步驟如下:

(1) 遍歷每一個網(wǎng)格,獲取其4 個節(jié)點的數(shù)值;

(2) 根據(jù)節(jié)點數(shù)值的分布情況確定該網(wǎng)格對應的索引值,采用二進制編碼的方式,將4 個節(jié)點的數(shù)值映射到一個4 位二進制數(shù)上,例如0000,0001,0010 等;

(3) 利用預先建立的查找表,在查找表中找到對應的線段幾何外形,例如水平線段、垂直線段和斜線段等;

(4) 最終將所有網(wǎng)格的線段合并起來構成整體壁面外形.

在SPARTA 中,被線段所圍的區(qū)域被視為燒蝕體,能夠與模擬粒子發(fā)生碰撞并燒蝕.主程序根據(jù)輸入條件如計算域尺寸、網(wǎng)格尺寸和加密信息等自動生成結構化隱式網(wǎng)格.其中每個網(wǎng)格為正方形(二維計算)或立方體(三維計算),每個二維正方形網(wǎng)格具有4 個節(jié)點,三維立方體網(wǎng)格具有8 個節(jié)點.在給定閾值之后,根據(jù)每個網(wǎng)格的4 個節(jié)點值,使用Marching Square 算法高效地劃分出壁面外形.以二維計算為例,初始燒蝕外形的生成如圖3 所示.該算法能夠根據(jù)實時的節(jié)點數(shù)值信息來高效地更新燒蝕外形.

圖3 SPARTA 燒蝕外形劃分示意圖Fig.3 The schematic of how SPARTA generates ablation geometry

在整個燒蝕模擬的過程中,模擬粒子不斷地在來流邊界生成并以指定的速度進行運動.在每一個時間步內(nèi),粒子進行直線運動或與其他粒子發(fā)生碰撞.同時部分粒子可能會與燒蝕壁面發(fā)生碰撞并依據(jù)壁面的溫度以一定的速度發(fā)生漫反射或鏡面反射.因此在每一時間步內(nèi),燒蝕外形上的不同部分會經(jīng)受不同程度的粒子碰撞,對應于不同程度的氣動熱環(huán)境.SPARTA 會記錄下每個燒蝕循環(huán)內(nèi)每一時間步內(nèi)壁面上發(fā)生碰撞的信息,例如與粒子發(fā)生碰撞的次數(shù),粒子碰撞前后的速度、動量等,并根據(jù)記錄下來的物理量對壁面進行燒蝕,例如計算粒子傳遞給壁面的總動能

此外,SPARTA 的燒蝕外形變化是在每次燒蝕循環(huán)后進行更新的,而每一個燒蝕循環(huán)由若干時間步組成.SPARTA 燒蝕循環(huán)流程圖如圖4 所示.圖4中N為當前模擬步數(shù),Nac為一個燒蝕循環(huán)所包含的模擬步數(shù),k為正整數(shù).以使用粒子動能變化作為燒蝕判據(jù)為例,每一個時間步內(nèi)粒子傳遞給壁面的動能都會被記錄下來,直到該燒蝕循環(huán)結束.燒蝕循環(huán)結束時,每一個壁面單元所吸收的能量將會被累加起來形成一個“燒蝕量”.程序會讀取該壁面單元所在網(wǎng)格的4 個節(jié)點值,從節(jié)點值中扣除該“燒蝕量”.在該燒蝕循環(huán)的最后,程序會根據(jù)更新后的節(jié)點值來生成新的燒蝕外形.例如當燒蝕循環(huán)步數(shù)Nac為100 時,每100 個時間步內(nèi)壁面外形的“耐久值”都在不斷減小,最后在第100 步時壁面外形會被更新.SPARTA 會將動能值轉化成與節(jié)點值相當量級的數(shù)值并從對應的網(wǎng)格節(jié)點值中扣除,直至該節(jié)點值為0,然后從該網(wǎng)格的其他節(jié)點值中繼續(xù)進行扣除.當一個網(wǎng)格的4 個節(jié)點值均為0 后,該網(wǎng)格內(nèi)的壁面單元即為“完全燒蝕”,此時該網(wǎng)格中不再存在壁面單元,粒子則可以進入該網(wǎng)格與“更里面”的壁面再次發(fā)生碰撞,并將該過程循環(huán)下去直至計算結束.因此在每一個燒蝕循環(huán)內(nèi),計算網(wǎng)格中節(jié)點值的分布會不斷發(fā)生變動,在循環(huán)結束時對壁面的外形進行更新.圖5 為SPARTA 模擬燒蝕過程的示意圖.

圖4 SPARTA 燒蝕循環(huán)示意圖Fig.4 Flowchart of ablation cycle in SPARTA

圖5 SPARTA 模擬燒蝕過程示意圖Fig.5 Simulated ablation process of SPARTA

圖5 中,左側為初始的壁面外形,其內(nèi)部的節(jié)點值均大于閾值,程序根據(jù)Marching Square 算法生成了該近似于正方形的外形.經(jīng)過一個燒蝕循環(huán)后,每個壁面單元經(jīng)受了不同程度的粒子碰撞,受到了不同程度的氣動加熱,因此節(jié)點值分布發(fā)生了變化,如圖5 右側所示.程序會依照新的網(wǎng)格節(jié)點值更新燒蝕后的壁面外形,并將該循環(huán)繼續(xù)進行下去.在SPARTA 的燒蝕循環(huán)中,最外層的壁面單元總是會最先被燒蝕.并且粒子碰撞最頻繁、碰撞粒子動能最高的區(qū)域燒蝕量最大,以此來反映真實、物理的燒蝕過程.

在SPARTA 中,粒子與壁面碰撞后反射的速度由基于壁面溫度的分布函數(shù)給出.在這種條件下,溫度較高的區(qū)域粒子反射的速度較大,傳遞給壁面的能量較小,因此反而燒蝕的速率較慢,與真實的情況相背;甚至當壁面溫度較高時,可能會出現(xiàn)碰撞后速度大于碰撞前的情況,該能量值成為負數(shù).因此本文對此項進行了優(yōu)化,將式(1)中碰撞后的動能項舍棄,只考慮入射的速度,從而保證溫度越高的區(qū)域燒蝕速率越快.

另外,由于SPARTA 中燒蝕外形的更新是在每次燒蝕循環(huán)結束之后進行,因此每個循環(huán)內(nèi)只有最外層的網(wǎng)格承受熱量,而內(nèi)層的網(wǎng)格并不會被燒蝕.在某些情況下,粒子傳遞給外層網(wǎng)格的能量可能足以在循環(huán)結束前燒蝕掉外層網(wǎng)格.然而,由于外形沒有及時更新,最外層的網(wǎng)格仍會存在并繼續(xù)與粒子碰撞.在這種情況下,本應燒蝕內(nèi)層網(wǎng)格的能量在燒蝕實際并不存在的外層網(wǎng)格.為了解決這個問題,本文在記錄能量時添加了一個縮放系數(shù) α,以確保在一個燒蝕循環(huán)內(nèi)無法完全燒蝕掉一個完整的網(wǎng)格,從而減少能量溢出的損失[27]

式中,fnum為真實分子與模擬粒子的比值,即一個模擬粒子代表多少真實的分子;Acell為網(wǎng)格的尺寸,如二維方形網(wǎng)格的邊長;dt為時間步長;Emax為模擬粒子最高的動能.

此外,當前模擬中計算網(wǎng)格的尺寸會影響燒蝕速率.fnum中真實分子的數(shù)量與氣體的數(shù)密度有關,如下式所示

為了保證計算效率,每個網(wǎng)格單元中模擬粒子的數(shù)目npercell通常是個定值.如果其值過大,模擬過程將耗費大量時間;而如果太小,DSMC 方法的統(tǒng)計特性將會引起較大的誤差.另外,在同一算例中,體積數(shù)密度 ρnum始終保持不變.在二維條件下,沿z方向的網(wǎng)格單元長度lz也是常數(shù),可以忽略.因此,fnum與x和y方向的網(wǎng)格單元大小lx和ly的乘積有關.通過代入和簡化,并假設用k表示其中的常數(shù)項,則式(1)中的系數(shù) β 變?yōu)?/p>

因此,β 正比于lx和ly的乘積與壁面元素長度a的比值.而a的大小與網(wǎng)格尺寸非常接近,因此lx和a可以同時消去,還剩下ly.這表明燒蝕速率實際上受到網(wǎng)格尺寸的影響.較大的網(wǎng)格意味著更快的燒蝕速率,而較小的網(wǎng)格則會導致較慢的燒蝕速率.燒蝕過程本身的速率受到了模擬參數(shù)的影響,將會產(chǎn)生誤差.

為了解決這一問題,對于二維模擬,系數(shù) β 的分母被修改為了a2.修改后保證了燒蝕速率不再受網(wǎng)格單元大小所影響,而是取決于燒蝕過程本身.

由于SPARTA 程序受到文件數(shù)據(jù)格式的限制,此前節(jié)點最大值被設定為255.為了考慮材料的傳熱特性,本文采用一維熱傳導方程以及燒蝕表面的熱平衡方程計算燒蝕速率[3].在燒蝕表面的熱交換過程包括高超聲速來流的氣動加熱熱流qaero、材料結構吸熱熱流qcond以及燒蝕表面的輻射熱流qrad

在DSMC 方法中,氣動加熱是由粒子與壁面碰撞前后傳遞的動能所實現(xiàn)的,因此氣動加熱熱流為粒子傳遞的動能流量之和

燒蝕表面的輻射熱流與壁面的溫度相關

其中 σ 為Stephan-Boltzmann 常數(shù),其值為 5.67×10-8W/(m2·K4);ε 為發(fā)射率,介于0～1 之間.當壁面的溫度到達最高點后,結構通過熱傳導方式所吸收的熱量可以認為完全用于結構的燒蝕,即

式中,V為燒蝕表面的后退速率,Q為壁面單位質量的材料燒蝕時所吸收的熱量.因此通過燒蝕表面的熱平衡方程可以求出燒蝕進行的速率V為

在以Marching Square 算法實現(xiàn)的燒蝕問題中,網(wǎng)格的節(jié)點值G將根據(jù)該燒蝕速率來設置

以此來實現(xiàn)更加物理的燒蝕過程.

2 典型外形燒蝕模擬

2.1 二維方柱燒蝕模擬

基于SPARTA 的燒蝕模型,本文就一些簡單的外形進行了燒蝕模擬.表1 展示了來流條件以及SPARTA 中的網(wǎng)格尺寸、時間步長等條件.

表1 SPARTA 計算參數(shù)Table 1 Simulation parameters in SPARTA

燒蝕材料采用典型的碳/碳復合材料的性能參數(shù),其密度為1800 kg/m3,表面發(fā)射率為0.8.材料在高空稀薄大氣環(huán)境下發(fā)生化學燒蝕,大氣的組成為77% 的N2和23% 的O2,反應體系中包括碳的燃燒、碳氮反應以及碳的升華等反應,在該計算條件下Q=36.95 MJ/kg.來流的氣動加熱熱流qaero由SPARTA 計算得到,在該條件下駐點處的氣動加熱熱流為6.586 MW/m2,除去熱輻射所耗散熱流后,剩余的熱量完全用于材料結構的燒蝕,由此可以計算出燒蝕速率約為0.098 3 mm/s,以及網(wǎng)格的節(jié)點數(shù)值,最后對SPARTA 的燒蝕過程進行“標定”.

SPARTA 輸出了整個燒蝕過程中燒蝕外形的變化圖6 和某一時刻流場的速度分布圖7.高超聲速來流會在方柱前方形成脫體弓形激波,激波正后方溫度急劇上升.

圖6 二維方柱燒蝕過程Fig.6 The ablation process of the rectangular cylinder

圖7 二維方柱繞流速度場云圖Fig.7 The velocity contour of the rectangular cylinder

在速度分布圖中,方柱前緣有較大一塊區(qū)域內(nèi)粒子速度明顯較低,而模擬粒子在方柱的棱角處的平均速度要高于前緣中心附近.因此棱角處表現(xiàn)出的平均溫度更高,燒蝕后退位移多于前緣中心,前緣將會逐漸燒蝕成弧形.與方柱側邊和后緣碰撞的粒子數(shù)量遠小于前緣與棱角,并且碰撞時的入射速度較小,所以靠近后緣處的燒蝕量較少,基本維持著原來的外形.根據(jù)燒蝕外形的變化圖6 可以發(fā)現(xiàn),燒蝕過程與速度場云圖所展示的情況基本符合,方柱截面逐漸變成了“子彈”形狀,前緣形成了圓弧形狀,反映了SPARTA 燒蝕模型的有效性.

2.2 二維圓柱燒蝕模擬

本文同時模擬了二維圓柱的燒蝕,其計算參數(shù)與方柱燒蝕相同.在圓柱燒蝕模擬中激波的形狀更加貼近前緣圓弧形,因此圓柱表面的溫度分布更加均勻.圓柱前緣駐點處承受的溫度最高,并且沿著兩側弧形溫度逐漸下降.肩部的溫度要低于頭部駐點,所以在圓柱表面上頭部駐點處的燒蝕量最大,離駐點越遠的地方燒蝕量越小,整體的燒蝕外形表現(xiàn)出弧形逐漸“變平”的特征.圖8 所展示的白色區(qū)域為圓柱截面的外形變化.圖9 展示了燒蝕過程中某一時刻的溫度場分布情況.

圖8 二維圓柱燒蝕過程Fig.8 The ablation of the two-dimensional cylinder

圖9 二維圓柱繞流溫度場云圖Fig.9 The temperature contour of the two-dimensional cylinder

在不考慮化學反應以及電離、離解反應的情況下,激波后駐點線上的最高溫度達到了25 000 K.前緣表面上的溫度遠高于側后方,因此在整個過程中后退位移要遠多于側后方,使得圓柱的前緣表面呈現(xiàn)明顯的扁平形狀,而后緣表面仍基本保持圓弧形狀.圖8 的燒蝕結果和溫度場的結果基本吻合.

2.3 二維球錐燒蝕模擬

本文就典型的球錐氣動外形[21,28]的高超聲速再入過程進行了模擬.該球錐頭的截面外形的頭部半徑為0.032 m,錐角9.8°,其示意圖如圖10 所示.因為球錐的外形為三維軸對稱外形,因此在模擬過程中選取其中一個對稱截面進行計算.球錐頭從70 km 高度的軌道以5800 m/s 的速度再入大氣層.在該高度下,大氣的密度為8.3×10-5kg/m3,氣體十分稀薄,DSMC 方法能夠很好地適用.表2 給出了模擬再入時的稀薄氣體環(huán)境條件以及SPARTA 模擬時的計算網(wǎng)格設置.

表2 來流條件以及SPARTA 設置Table 2 Free stream conditions and SPARTA parameters of the simulated blunt cone

圖10 球錐截面氣動外形Fig.10 Aerodynamics geometry of the blunt cone section

由于選取的球錐外形的尺寸較小,相應地,為保證計算的準確度,網(wǎng)格尺寸設置得足夠小,并且小于分子的平均自由程;時間步長的選取同樣小于分子平均自由時間.壁面采用恒溫條件,在燒蝕過程中壁面的平均溫度可達到4000 K 左右,因此認為在燒蝕過程中壁面溫度保持在該水平上.來流的氣動加熱熱流最高為126.009 MW/m2,此時沿駐點線的燒蝕速率約為1.72 mm/s.

燒蝕過程中某時刻的溫度云圖如圖11 所示.SPARTA 計算的最終燒蝕外形以及參考文獻[21]中的原始外形和燒蝕外形由圖12 給出.燒蝕過程中初步考慮了O2,N2,O,N,NO 五組分的化學反應.不過圖11 的結果中顯示的溫度最高達到了14 000 K以上,事實上在該溫度下氣體已經(jīng)發(fā)生了電離和離解反應,反應過程會吸收很大一部分熱量,因此真實情況下的最高溫度將低于14 000 K.

圖11 球錐截面燒蝕溫度場云圖Fig.11 The temperature contour of the ablated blunt cone section

圖12 SPARTA 預測燒蝕外形與文獻[21]結果對比Fig.12 The surface recession predicted by SPARTA and the comparison with the result in Ref.[21]

由圖12 可以看出,在頭部駐點線附近,SPARTA模擬得到的燒蝕外形與參考文獻結果相符,駐點后退量為9.50 mm,而文獻中的燒蝕后退量為9.79 mm,誤差在3%以內(nèi).在再入過程中,球錐前方形成了弓形脫體激波,激波后頭部承受的熱流最大,因此燒蝕量最大.肩部區(qū)域所承受的熱流低于頭部區(qū)域,因此燒蝕后退位移少于頭部.對比參考文獻中的燒蝕外形,SPARTA 計算的燒蝕外形整體呈現(xiàn)較為平滑的弧線型,而文獻中的結果在肩部附近呈現(xiàn)“削平”的狀態(tài).本文結果中肩部區(qū)域的燒蝕量明顯較少,這是由于參考文獻中考慮了轉捩現(xiàn)象的影響.肩部的轉捩點后的熱流會迅速升高,局部的燒蝕量較大,因此肩部區(qū)域會產(chǎn)生凹陷.并且由于外形的變化,轉捩點不斷沿肩部移動,使得整個肩部區(qū)域的燒蝕量要多于不考慮轉捩現(xiàn)象的情況,導致外形越來越尖.所以SPARTA 計算的外形仍呈圓弧形,而參考文獻為典型的“雙錐形”.SPARTA 計算中暫未考慮轉捩現(xiàn)象的影響,因此導致了部分誤差.如何將復雜的轉捩機制加入SPARTA 程序中也是今后的改進之處.

3 帶有微小粗糙元的二維斜楔體燒蝕模擬

本文模擬了二維條件下帶有微小的方塊狀粗糙元的斜楔體燒蝕[29],其截面幾何外形如圖13 所示.

圖13 斜楔體截面外形示意圖Fig.13 Geometry of the wedge with an obstacle

相對于整個斜楔體而言,方塊狀粗糙元的尺寸非常小,可以視作飛行器表面上的螺栓等小零件.來流條件參考NASA 于2008 年的一項實驗數(shù)據(jù)[29]給出.具體的計算參數(shù)如表3 所示,

表3 來流條件以及SPARTA 設置Table 3 Free stream conditions and SPARTA parameters of the simulated wedge

為了還原該幾何外形,保證微小粗糙元能夠以較好的分辨率重現(xiàn)在SPARTA 中,并綜合計算量的考慮,網(wǎng)格的尺寸設置為0.000 2 m.相應地,時間步長為1×10-8s.壁面的初始溫度設置為360 K.計算中考慮了高溫條件下大氣組分中N2,O2,N,O,NO五組分的化學反應.而化學燒蝕模型考慮碳的氧化、碳氮反應以及碳的升華反應.在42.5 km 的高度下,大氣的稀薄程度有所降低,此時單位質量的燒蝕材料完全熔化時所吸收的熱量Q取值為10.20 MJ/kg.駐點處的熱流為5.37 MW/m2,燒蝕速率約為0.081 mm/s.

通過初步的分析,斜楔體的頭部較為尖銳,來流在前方形成脫體弓形激波,頭部會承受較高的熱流.同時微小粗糙元附近也會形成局部的高熱流區(qū)域,因此燒蝕會首先發(fā)生在頭部以及粗糙元區(qū)域.SPARTA計算結果如圖14～圖16 所示,圖14 呈現(xiàn)了燒蝕過程中微小粗糙元的外形變化,圖15 和圖16 分別展示了流場中溫度和克努森數(shù)的分布情況.

圖14 帶有局部粗糙元的斜楔體燒蝕過程Fig.14 The ablation process of the obstacle on the wedge

圖15 帶有局部粗糙元的斜楔體溫度場云圖Fig.15 The temperature contour of the wedge with an obstacle

圖16 帶有局部粗糙元的斜楔體克努森數(shù)云圖Fig.16 The Knudsen number contour of the wedge with an obstacle

圖14 中粗糙元區(qū)域的外形變化非常明顯,燒蝕開始后左上角迅速變成了圓弧形,并且粗糙元前方的角區(qū)開始向內(nèi)凹陷;然后粗糙元逐漸被削平,前方凹陷的部分也逐漸加深;最后整個粗糙元幾乎完全被燒蝕,和斜楔體表面共同后退.

通過溫度場云圖15 可以看出,斜楔體頭部和粗糙元前方溫度較高,最高溫度達到了3800 K.高溫區(qū)域主要分布在斜楔的頭部區(qū)域以及方塊狀粗糙元的前緣區(qū)域.溫度場的分布情況表明主要燒蝕的區(qū)域集中在頭部和粗糙元前緣區(qū)域.

圖16 為流場中克努森數(shù)的分布.通常情況下,當克努森數(shù)大于0.05 時,流體的連續(xù)介質假設將會失效,此時基于該假設的數(shù)值方法也不再適用[30].從圖16 可以看出,盡管云圖中有部分噪聲,但是絕大部分區(qū)域的克努森數(shù)均在0.1 以上.特別是在粗糙元區(qū)域后方以及靠近斜楔表面的區(qū)域,克努森數(shù)已經(jīng)遠超連續(xù)流體的范疇.由此可見,在42.5 km 的高度下,即使飛行器周圍所處的流體環(huán)境可能是整體上連續(xù)的,但靠近飛行器表面的流場仍處于過渡區(qū)甚至是稀薄區(qū).因此,DSMC 方法在該情況下依舊是重要的研究方法.

從該算例也可以反映出,飛行器頭部表面局部粗糙區(qū)域容易先于其他部分發(fā)生燒蝕.在飛行器頭部布置相關的部件導致形成類似于該粗糙元形狀的區(qū)域時,需要特別注意該部分的熱防護系統(tǒng)設計.

4 結論

本文采用基于分子動理論的DSMC 方法對高超聲速飛行器再入過程中的微燒蝕問題展開了研究,針對多種典型氣動外形進行了燒蝕計算與分析,并總結了飛行器表面微燒蝕的主要特征.

首先,本文定性分析了二維條件下方柱和圓柱的燒蝕.通常情況下,駐點處的溫度最高,因此燒蝕量也會最大.而模擬結果顯示,局部尖銳或粗糙的區(qū)域熱流會更加集中,具有更快的燒蝕速率.

其次,本文模擬了二維條件下典型軸對稱球錐的再入燒蝕.在目前設置條件下,計算得到的駐點附近燒蝕速率與文獻結果基本相符.不過由于未考慮轉捩現(xiàn)象以及復雜化學反應的影響,本文預測的燒蝕量產(chǎn)生了部分誤差,結果中激波后溫度也偏高.因此如何在程序中綜合考慮轉捩現(xiàn)象、電離以及離解反應的影響也是重要的改進方向.

最后,本文模擬了二維條件下帶有微小粗糙元的斜楔體的燒蝕過程.粗糙元附近存在有相當稀薄的流場區(qū)域,表明了DSMC 方法在微燒蝕問題中的重要意義.整個斜楔體的頭部以及粗糙元前方熱流較高,率先發(fā)生燒蝕,粗糙元區(qū)域會逐漸被削平并形成凹陷.因此熱防護系統(tǒng)的設計要考慮到局部的粗糙元,尤其是飛行器頭部較為重要的部件.

本文初步驗證了DSMC 方法對于飛行器再入過程中燒蝕問題研究的可行性.此外,Marching Square 算法能夠非常高效地更新燒蝕外形,與基于粒子模擬的DSMC 方法契合較好,在燒蝕耦合計算方面具有很大的潛力.在本文微燒蝕研究的基礎上,可以預測高超聲速飛行器再入過程中的燒蝕速率,為飛行器熱防護系統(tǒng)的設計提供理論基礎或者為在軌航天器退役后遺骸的再入燒蝕問題提供處理方案.不過本文目前的計算方法仍存在一定誤差,并且未考慮到更為復雜的化學反應體系.將以上屬性加入到燒蝕模擬中將成為本文接下來的主要工作方向.