李舟波，方蜀州

(北京理工大學(xué)宇航學(xué)院，北京 100081)

0 引 言

航天器的表面實(shí)際上存在很多不連續(xù)的情況，不同幾何形狀的縫隙經(jīng)常出現(xiàn)在各種氣動表面上，如航天器殼體、燃?xì)廨啓C(jī)通道、熱交換器和微電子芯片中的肋狀表面[1]。近年來，縫隙結(jié)構(gòu)更被廣泛應(yīng)用于超燃沖壓發(fā)動機(jī)，以延長超音速空氣和噴油在燃燒室中的停留時間[2]。但是，縫隙結(jié)構(gòu)會嚴(yán)重影響氣流的流動狀態(tài)與傳熱特性。這是因?yàn)樵诳p隙入口處由于邊界層分離與再附將產(chǎn)生局部高熱流區(qū)；在某些工況下縫隙入口處會產(chǎn)生激波；縫隙內(nèi)會卷起旋渦結(jié)構(gòu)等。同時，旋渦結(jié)構(gòu)會向縫隙中卷入大量熱量，使得傳熱特性由熱傳導(dǎo)變?yōu)闊醾鲗?dǎo)加對流傳熱，縫隙壁面熱流隨之大大增加。因此，高超聲速飛行器表面縫隙的流場結(jié)構(gòu)和熱流分布受到了廣泛的關(guān)注。近年來，國內(nèi)外許多研究者對縫隙進(jìn)行了各種數(shù)值和實(shí)驗(yàn)研究[3-15]。

Sudarshan等[3]研究了圓形和拋物形縫隙的傳熱變化和內(nèi)部激波特性，在馬赫數(shù)分別為 10.9和 10.2 的高超聲速風(fēng)洞中進(jìn)行了實(shí)驗(yàn)。結(jié)果發(fā)現(xiàn)，縫隙結(jié)構(gòu)使機(jī)頭表面的熱通量大幅減少，這對緩解氣動加熱是非常有效的。Dhankarghare等[4]將支板組合成一個支板縫隙，并在不同馬赫數(shù)和長徑比的情況下，將其與壁面縫隙進(jìn)行比較。結(jié)果表明，兩種縫隙在馬赫數(shù)和寬高比較低時都表現(xiàn)出內(nèi)部壓力振蕩，而在高馬赫數(shù)和寬高比時振蕩減小。Krishna等[5]通過實(shí)驗(yàn)研究了在超聲速條件下，矩形縫隙后壁面幾何形變的影響。研究表明，縫隙壁面坡道結(jié)構(gòu)使得縫內(nèi)整體聲壓和脈動壓力分別降低了10 dB和50%。

國內(nèi)方面，王世芬等[6]、唐貴明[7]分別在激波風(fēng)洞上開展了縫隙熱流測量實(shí)驗(yàn)，研究了平板上的矩形縫隙和單個橫縫內(nèi)的熱流分布。張昊元等[8]針對高超聲速飛行器前緣防熱瓦結(jié)構(gòu)，通過求解三維可壓縮Navier-Stokes方程進(jìn)行了數(shù)值模擬，研究了縫隙誘導(dǎo)形成的三維旋渦的空間分布特征和旋渦運(yùn)動對物面氣動加熱的影響規(guī)律。肖志祥等[9]通過數(shù)值求解雷諾平均的Navier-Stokes方程組，數(shù)值模擬導(dǎo)彈跨、超聲速流動，研究導(dǎo)彈環(huán)向凹槽等典型分離流動的被動減阻措施對阻力縮減的影響。通過分析計算結(jié)果發(fā)現(xiàn)環(huán)向凹槽有利于降低摩擦阻力，但底阻減小不明顯。陸海波等[10]對高超聲速飛行器鼻錐使用迎風(fēng)凹腔結(jié)構(gòu)作為熱防護(hù)系統(tǒng)時，凹腔結(jié)構(gòu)的防熱效能進(jìn)行了數(shù)值研究。結(jié)果表明迎風(fēng)凹腔結(jié)構(gòu)能夠有效對高超聲速飛行器的鼻錐尤其是駐點(diǎn)區(qū)域進(jìn)行冷卻，凹腔越深，其冷卻效果越好。鼻錐氣動加熱的最大熱流并不在尖銳唇緣的頂點(diǎn)，而是位于凹腔內(nèi)的側(cè)壁面上，凹腔的深度(L)變化對最大熱流的出現(xiàn)位置影響很小。除非凹腔很淺(L/D<0.5)，凹腔底面的熱流值都非常小，基本可以忽略。Wang等[11]通過風(fēng)洞實(shí)驗(yàn)和數(shù)值計算，揭示了圓形縫隙可壓縮流動的換熱特性。在縫隙后緣附近熱流比無縫隙模型顯著增加，同時隨著縫隙直徑(D) 與深度(H) 之比的減小，縫隙內(nèi)的升溫幅度減小。

可以看到，在臨近空間(60～100 km)中的高超聲速縫隙結(jié)構(gòu)的研究仍然很少，這可能是高超聲速條件下可靠的實(shí)驗(yàn)數(shù)據(jù)測量困難和臨近空間稀薄氣體效應(yīng)的原因。但隨著計算機(jī)技術(shù)的發(fā)展，尤其是直接模擬蒙特卡洛(Direct simulation Monte Carlo,DSMC) 方法的計算效率增強(qiáng)，對稀薄縫隙繞流的研究也越來越多。Jin等[12]采用DSMC方法對一組具有圓形分離/再附角的二維矩形縫隙模型進(jìn)行了高超聲速稀薄流動數(shù)值模擬。結(jié)果表明，對于圓形分離角，外部氣體分子更容易進(jìn)入縫隙，從而促進(jìn)了與縫隙內(nèi)氣體分子的動量和能量交換。而圓形再附角則會使得氣體分子從縫隙內(nèi)逸出，削弱動量和能量的交換。Guo等[13]采用直接模擬蒙特卡洛法,研究了帶后掠角縫隙的流動特性，模擬條件為海拔高度20～60 km,馬赫數(shù) 6～20?？梢杂^察到，增加高度，縫隙內(nèi)二次回流區(qū)的旋轉(zhuǎn)方向會發(fā)生逆轉(zhuǎn)。增加馬赫數(shù)，縫隙內(nèi)的速度、壓力和溫度均增加。Palharini等[14]采用DSMC方法研究了80 km高度稀薄過渡流區(qū)高超聲速縫隙繞流問題，探索縫隙幾何尺寸(長深比)對縫隙內(nèi)部流場結(jié)構(gòu)、壁面熱流和壓力的影響規(guī)律。結(jié)果表明，在80 km高度稀薄過渡流區(qū)，幾何尺寸對高超聲速縫隙繞流的影響與連續(xù)流區(qū)存在不同，主要表現(xiàn)為，當(dāng)長深比等于4時縫隙內(nèi)部就出現(xiàn)了雙渦結(jié)構(gòu)和剪切層在縫隙底面再附，而連續(xù)流區(qū)產(chǎn)生類似現(xiàn)象對應(yīng)的長深比為14。靳旭紅等[15]采用DSMC方法研究了70 km,75 km,80 km和90 km等4個飛行高度下稀薄流區(qū)高超聲速縫隙流動問題。結(jié)果表明：上述飛行高度下,外部流動的分離和再附在縫隙內(nèi)部形成一個充滿腔體的單渦結(jié)構(gòu)；稀薄氣體效應(yīng)對縫隙內(nèi)部流動結(jié)構(gòu)和壁面熱流影響明顯，隨著高度的增加，主渦渦心上移,其形狀逐漸變得“扁長”，右上角逐漸變尖，熱流越來越集中分布于縫隙下游側(cè)面的頂部區(qū)域。

從目前的研究成果來看，上述文獻(xiàn)的研究模型主要為不發(fā)生形變的矩形結(jié)構(gòu)，并且研究內(nèi)容主要集中在流場特性方面，而對表面氣動參數(shù)的研究較少。但是在實(shí)際飛行中，高超聲速流場中縫隙內(nèi)表面會受到嚴(yán)重的氣動加熱，縫隙形狀結(jié)構(gòu)尤其是兩側(cè)壁面會被強(qiáng)烈的表面熱流燒蝕而發(fā)生變化。因此，對形變縫隙進(jìn)行流場結(jié)構(gòu)和氣動加熱分析顯得十分重要。

本文采用直接蒙特卡洛方法進(jìn)行仿真計算，將縫隙模型簡化為壁面可變傾斜角的梯形模型。研究了在稀薄高超聲速環(huán)境中壁面傾角改變對流場結(jié)構(gòu)和氣動熱特性分布的影響規(guī)律，對飛行器表面優(yōu)化設(shè)計提供進(jìn)一步的理論幫助。

1 數(shù)值方法

DSMC方法的基本思想是：用有限個仿真分子代替真實(shí)氣體分子，并在計算機(jī)中存儲仿真分子的位置坐標(biāo)、速度分量以及內(nèi)能等信息，并在模擬分子運(yùn)動、碰撞以及與邊界相互作用的過程中不斷改變這些信息的具體數(shù)值，最后通過統(tǒng)計網(wǎng)格內(nèi)仿真分子的運(yùn)動狀態(tài)實(shí)現(xiàn)對真實(shí)氣體流動問題的模擬，流場的宏觀參數(shù)可以由統(tǒng)計平均求得。該方法的關(guān)鍵是在時間步長內(nèi)將分子的運(yùn)動與碰撞解耦。在時間步長內(nèi)讓所有分子運(yùn)動一定的距離并考慮在邊界的反射，然后計算此時間步長內(nèi)具有代表性的分子間的碰撞。

本文數(shù)值計算工具選擇的是美國桑迪亞國家實(shí)驗(yàn)室的開源SPARTA程序，并且使用北京并行云科學(xué)技術(shù)計算中心的工作站來進(jìn)行數(shù)值模擬計算。模擬分子采用可變徑軟球模型，并考慮到了空氣分子的轉(zhuǎn)動能和平動能之間的能量傳遞過程。分子內(nèi)能和平動能之間的能量傳遞采用標(biāo)準(zhǔn)Larsen-Borgnakke模型。另外，本文選擇的時間步長為當(dāng)?shù)刈杂煞肿悠骄鲎矔r間的三分之一。模擬分子數(shù)在每個網(wǎng)格內(nèi)最少為10個。

同時，由于采用DSMC方法進(jìn)行數(shù)值計算時，模擬分子數(shù)和網(wǎng)格劃分引起的計算量非常大，需要開展并行計算。一般并行算法采用區(qū)域分解，根據(jù)計算的處理器數(shù)(或CPU核心數(shù))的多少將計算區(qū)域劃分為等量的子區(qū)域，每個處理器在其分配的子區(qū)域內(nèi)部獨(dú)立地計算模擬分子的碰撞和遷移，離開子區(qū)域的模擬分子把攜帶的信息傳遞給對應(yīng)子區(qū)域的處理器。但DSMC方法的一個顯著特征是模擬分子在流場中的分布是不斷變化的，計算時間依賴于子區(qū)域內(nèi)的分子個數(shù)、分子的碰撞次數(shù)以及分子與物面的碰撞次數(shù)。因?yàn)槟M前不可能預(yù)測到不同子區(qū)域內(nèi)的分子個數(shù)，因此在區(qū)域分解時將各子區(qū)域的計算負(fù)載劃分得比較均衡就特別困難。如果簡單地將流場均勻劃分成若干同等大小的子區(qū)域，就有可能造成處理器之間的負(fù)載分布不均勻，負(fù)載最少的處理器就要花費(fèi)較長的時間等待負(fù)載最大的處理器，造成同步等待的時間顯著增加。這種方式并不適合DSMC并行模擬。

本文發(fā)展了一套靜態(tài)隨機(jī)負(fù)載平衡方法用于解決不同處理器之間的計算時間同步問題，該方法基于概率近似原理，將計算區(qū)域的全部網(wǎng)格平均分配給指定的所有處理器。由于采用相同的概率隨機(jī)選取流場網(wǎng)格，按照均分后的數(shù)量分配給每個處理器，因此當(dāng)計算區(qū)域的網(wǎng)格數(shù)量較多時，每個處理器包含近似相等的物體邊界、高密度流動區(qū)域和稀薄氣體區(qū)域的網(wǎng)格數(shù)，這樣每個處理器的計算負(fù)載也非常接近。

2 仿真計算

為了驗(yàn)證SPARTA程序的可行性，首先對二維的平板模型進(jìn)行了模擬，平板的長度為100 mm，厚度為5 mm，并將模擬計算得到的結(jié)果與風(fēng)洞實(shí)驗(yàn)[16]測量得到的結(jié)果進(jìn)行比較。模擬仿真的來流條件設(shè)置為風(fēng)洞的實(shí)驗(yàn)條件：來流氣體為氮?dú)猓杂蓙砹鞯乃俣葹?504 m/s，平均自由程為1.6 mm，溫度為13.32 K，氣體分子數(shù)密度3.716×1020m-3。

計算域?qū)挾仍O(shè)置為180 mm，水平方向設(shè)置360個網(wǎng)格；高度設(shè)置為205 mm，垂直方向設(shè)置410個網(wǎng)格。每個網(wǎng)格的長寬都是0.5 mm。計算域與平板示意圖如圖1所示。平板壁面溫度為290 K。模擬粒子的仿真比例FNUM設(shè)置為4.56×1011，即一個模擬粒子代表4.56×1011個真實(shí)分子。時間步長設(shè)置為3.102×10-7s。同時，采用自適應(yīng)結(jié)構(gòu)網(wǎng)格，分子碰撞采用VHS模型，NTC方法進(jìn)行取樣，能量交換采用Larsen-Borgnakke模型，壁面恒溫，采用漫反射模型。不考慮氣體化學(xué)反應(yīng)。

圖1 平板測試算例示意圖Fig.1 Calculation example of plate tests

平板上平面的壓力仿真結(jié)果與文獻(xiàn)[17]中DAC方法的仿真結(jié)果以及風(fēng)洞實(shí)驗(yàn)結(jié)果對比，如圖2所示。SPARTA方法得到的結(jié)果與DAC的結(jié)果比較吻合，但是都比實(shí)驗(yàn)結(jié)果略大，但壓力分布的趨勢與實(shí)驗(yàn)結(jié)果大致相同。計算結(jié)果表明本文選擇的數(shù)值方法和計算程序?qū)τ谶@類高超聲速流場的計算模擬具有較好的可信度。

圖2 仿真結(jié)果比較Fig.2 Comparison of simulation results

3 模型和流場參數(shù)

3.1 幾何定義

臨近空間航天器上的表面不連續(xù)結(jié)構(gòu)或縫隙遠(yuǎn)遠(yuǎn)小于航天器本身的尺寸，因此將航天器表面上的縫隙結(jié)構(gòu)簡化為無限長的平面上的矩形縫隙。本文僅對二維情況的縫隙模型進(jìn)行分析，模型的幾何示意圖如圖3所示。

圖3 縫隙模型外形示意圖Fig.3 Schematic diagram of the gap model

傾角θ表示前壁面和底面的夾角，傾角θ分別取30°,45°,60°,75°,90°(即矩形縫隙)。定義縫隙的長度為L=9 mm，深度為D=3 mm?？p隙上游的平板長度為3D，下游平板的長度為D。

3.2 來流條件

本文選取的高度為70 km，來流的速度、壓強(qiáng)、溫度、密度、粒子數(shù)密度以及平均自由程均采用美國標(biāo)準(zhǔn)大氣(1976) 確定，計算工況如表1所示。氣體成分為78%的N2和22%的O2。根據(jù)70 km自由來流的平均自由程和縫隙模型的特征尺寸，得到全局克努森數(shù)Kn大約為0.3，該縫隙模型的繞流問題屬于稀薄流范疇。本文設(shè)定縫隙模型的壁面溫度Tw為220 K。

表1 自由流條件Table 1 Free stream conditions

4 流場和壁面參數(shù)計算公式

4.1 流場特性分析

在本文的研究中主要關(guān)注的流場特性參數(shù)有流場中的馬赫數(shù)、壓力和溫度。DSMC方法中流場的宏觀參數(shù)是由計算域內(nèi)每個網(wǎng)格單元中的模擬粒子的參數(shù)計算平均得到的。

每個網(wǎng)格單元的質(zhì)心速度矢量的計算方法為：

(1)

式中：N,mj和cj分別為每個網(wǎng)格單元中粒子的數(shù)量、質(zhì)量和速度矢量。粒子的熱運(yùn)動速度或擾動速度c′定義為

c′≡c-c0

(2)

網(wǎng)格中平動溫度TT的計算公式為：

(3)

轉(zhuǎn)動溫度TR與網(wǎng)格單元的轉(zhuǎn)動能量εR和轉(zhuǎn)動自由度ζR有關(guān)，計算公式為：

(4)

每個網(wǎng)格單元的壓力的計算方法為：

(5)

每個網(wǎng)格單元的密度的計算方法為：

(6)

4.2 壁面參數(shù)

本文在研究模型的壁面特性時主要關(guān)注表面壓力系數(shù)和表面熱流系數(shù)。

模型表面的壓力pw由每個時間步內(nèi)入射和反射的粒子在物面上的法向動量之和得出，表達(dá)式如下:

(7)

式中：A為物面面積；N為單位時間內(nèi)在單位面積上與物面碰撞的粒子的個數(shù)；v是j粒子在物面法線法向的速度分量；i,r分別表示入射和反射粒子。

壓力系數(shù)Cp的定義如下：

(8)

熱流qw是通過計算碰撞壁面的模擬粒子的凈能量通量來計算的，與入射、反射粒子的平動能、轉(zhuǎn)動能和振動能有關(guān)，表達(dá)式如下：

(9)

(10)

式中：eRj和eVj分別代表粒子的轉(zhuǎn)動能和振動能。

傳熱系數(shù)Ch的定義如下：

(11)

5 流場和壁面參數(shù)計算公式

5.1 流場特性分析

5.1.1流線分析

為了清晰地描述縫隙的流動特征，縫隙內(nèi)回流區(qū)的變化是研究的重點(diǎn)。流線圖能較好反應(yīng)出其變化。圖4顯示了不同側(cè)壁傾角下的流線分布。

圖4 縫隙內(nèi)流線分布圖Fig.4 Distribution of streamline traces inside the gaps

從圖中可以看出，當(dāng)θ=90°時，縫隙內(nèi)部存在連在一起的兩個旋渦，此時仍然認(rèn)為其是一個回流區(qū)且填滿整個縫隙，上方流線直接跨過縫隙。當(dāng)θ=75°時，隨著傾角的減小，左、右旋渦強(qiáng)度都開始減小，在縫隙上方的外部流體逐漸進(jìn)入到縫隙內(nèi)部，縫隙中間的旋渦受到剪切層的擠壓作用而往縫隙底部凹陷。在這兩種工況下，旋渦仍然是連在一起，是上一種情況的延續(xù)，在縫隙內(nèi)表現(xiàn)為單個旋渦的特性。但兩個渦心的距離開始增大。當(dāng)θ=60°時，在縫隙前、后側(cè)出現(xiàn)了明顯分開的兩個旋渦，此時剪切層也深入縫隙底部。縫隙右側(cè)旋渦尺寸變化不大，但縫隙的左側(cè)旋渦明顯變小。θ=45°時，左側(cè)角落出現(xiàn)非常小的旋渦。對于右側(cè)旋渦區(qū)域，它在每種工況下，其姿態(tài)和流向長度變化不大。當(dāng)傾角繼續(xù)減小，θ=30°時，縫隙的左側(cè)角落幾乎被剪切層充滿，左旋渦徹底消失，在縫隙左側(cè)壁面上邊界層沒有流動分離。

綜合比較圖4，當(dāng)傾角從90°變化到30°時，縫隙內(nèi)部的回流區(qū)發(fā)生了明顯的變化，并以60°和 30°作為分界點(diǎn)。當(dāng)θ≤60°時，縫隙內(nèi)部出現(xiàn)了兩個回流區(qū)，且上方的流體能夠進(jìn)入到縫隙的底面；θ≤30°，左側(cè)的旋渦區(qū)域完全消失。同時，可以看出，旋渦區(qū)域的旋轉(zhuǎn)方向均為順時針。縫隙右側(cè)的角部渦流與前壁的傾角無關(guān)，都保持大致相同的尺寸?？梢酝茢?，改變θ的大小對縫隙后壁面區(qū)域的影響較小，這會在下文表面參數(shù)分析中進(jìn)行驗(yàn)證。

最后，通過傾斜縫隙后壁面也可以消除右側(cè)的旋渦區(qū)域。為此，以30°的角度傾斜縫隙前后壁面，然后使用DSMC方法計算。圖5顯示了兩側(cè)壁面傾角都為30°時的流線分布圖。通過與圖4中的結(jié)果進(jìn)行比較，發(fā)現(xiàn)當(dāng)后壁面傾角為30°時，縫隙右側(cè)旋渦區(qū)域也完全消失。因此，一個重要的結(jié)論是，通過以適當(dāng)?shù)慕嵌葍A斜縫隙前后壁面，可以完全消除縫隙內(nèi)的旋渦區(qū)域。

圖5 流線分布圖Fig.5 Distribution of streamline traces

5.1.2馬赫數(shù)分析

圖6展示了不同前壁傾角θ下縫隙模型的流場馬赫數(shù)等值線圖。

從圖中可以看出，當(dāng)θ=90°時，在縫隙內(nèi)的大部分區(qū)域，其馬赫數(shù)小于 1，說明在內(nèi)部回流區(qū)的氣流速度很小，且由于壁面的黏性作用，在縫隙中間的流體速度大于兩側(cè)的流體，縫隙底面中間位置受到的侵蝕最大。當(dāng)θ減小時，縫隙外剪切層從縫隙中部緩慢地深入，導(dǎo)致最大流速區(qū)往前移動，亞聲速區(qū)被縫隙外剪切層往下壓縮到縫隙底部，流場內(nèi)的整體速度也開始增大。這是由于縫隙前壁面往上游傾斜時，一方面使縫隙具有更寬的開口，一方面讓縫隙入口處流體的再附點(diǎn)從縫隙底面移動到縫隙前壁面，外部流動可以更深地作用縫隙內(nèi)部流動，因此可以看到縫隙內(nèi)流體速度開始增加。當(dāng)θ=30°時，亞聲速區(qū)只存在于縫隙右半?yún)^(qū)域，縫隙外超聲速流場進(jìn)入縫隙左側(cè)區(qū)域。

5.1.3密度云圖

圖7為不同前壁面傾角時縫隙內(nèi)部的無量綱密度(當(dāng)?shù)孛芏扰c來流密度之比)等值線圖。

從圖中可以看出，當(dāng)θ=90°時，縫隙左側(cè)完全被低密度區(qū)占據(jù)，其值小于來流密度。又由于氣流在后壁面被壓縮，縫隙右側(cè)流體密度明顯比左側(cè)高很多，密度最高處在縫隙右側(cè)角落。隨著傾角θ的減小，流場內(nèi)的整體密度變化不大，只是由于左側(cè)壁面向前傾斜，縫隙內(nèi)的低密度區(qū)范圍也隨之往上游延伸。而在縫隙右側(cè)，隨著傾角減小，密度只在縫隙右側(cè)角落處略微增大?？梢钥闯銮氨诿鎯A角的改變對流場密度的影響較小。

5.1.4壓力云圖

圖8為不同前壁面傾角下縫隙內(nèi)部的無量綱壓力(當(dāng)?shù)貕毫εc來流壓力之比)等值線圖。

從圖中可以看出，當(dāng)θ=90°時，在縫隙左側(cè)的壓強(qiáng)較低。這是因?yàn)楫?dāng)來流經(jīng)過縫隙的前壁面時，相當(dāng)于流經(jīng)一個后向臺階，流場的面積突然增加，自由來流在這里產(chǎn)生流動分離，因此在前壁面頂點(diǎn)處附近出現(xiàn)了局部壓強(qiáng)較低的區(qū)域。隨著前壁傾角的減小，縫隙左側(cè)的低壓區(qū)域被顯著壓縮。隨著縫隙前壁面的傾斜，縫外高壓氣流加速流入，壓力呈上升趨勢。

在縫隙右側(cè)，還可以觀察到在后壁面頂點(diǎn)的左上方流場處出現(xiàn)了一個局部高壓區(qū)域。從流線圖中可以看出，這是高速流體撞擊在后壁面頂點(diǎn)區(qū)域所致。隨著θ的減小，該高壓區(qū)域開始緩慢增大。這也是因?yàn)榍氨诿鎯A斜，進(jìn)入縫隙內(nèi)的流體速度開始增大。其最終會在右側(cè)壁面開始重新被壓縮，并且壓縮效應(yīng)隨速度的增大而增大，所以，前壁傾角的減小會增大縫隙右側(cè)的壓力。又因?yàn)橄”⌒?yīng)，減弱了縫隙內(nèi)的氣流速度對氣體壓縮的影響，所以，其改變幅度也并不顯著。

圖6 縫隙內(nèi)馬赫數(shù)分布Fig.6 Mach number distribution in the gaps

最后，還可以明顯發(fā)現(xiàn)前壁傾角的改變對流場壓強(qiáng)的影響比對流場密度的影響更大。

5.1.5溫度云圖

圖9展示了不同前壁傾角下縫隙內(nèi)部的無量綱溫度(當(dāng)?shù)販囟扰c來流溫度之比)等值線圖。

在縫隙左側(cè)，當(dāng)θ=90°時，低溫區(qū)占據(jù)縫隙左下方區(qū)域，也正好是左旋渦區(qū)域的范圍。當(dāng)傾角減小時，縫隙前壁面附近的低溫區(qū)開始減小。高溫區(qū)開始前移，縫隙內(nèi)的整體溫度也開始有所提高。當(dāng)θ=30°時，隨著左側(cè)旋渦區(qū)域消失，縫隙外的高溫氣流進(jìn)入縫隙的左側(cè)角落，區(qū)域溫度明顯升高，遠(yuǎn)大于縫隙右側(cè)角落。同樣在縫隙右側(cè)上方，由于壓縮作用形成了一個高溫區(qū)。當(dāng)傾角減小時，右側(cè)的溫度分布基本沒有變化，只是縫隙右側(cè)上方形成的高溫區(qū)范圍略微增大。這驗(yàn)證了前壁傾角的改變對縫隙后壁面附近的影響較小。

最后，為了更深入了解其變化，研究其縫隙內(nèi)部的在稀薄氣體流動中的熱不平衡效應(yīng)，對轉(zhuǎn)動溫度也進(jìn)行研究。圖10展示了不同前壁傾角下縫隙內(nèi)部的無量綱轉(zhuǎn)動溫度(當(dāng)?shù)剞D(zhuǎn)動溫度與來流轉(zhuǎn)動溫度之比)等值線圖。

綜合比較圖9，可以看出，當(dāng)θ=90°時，在縫隙左側(cè)，由于整體溫度都較低，平動溫度TT與轉(zhuǎn)動溫度TR的差距并不明顯，該處的熱不平衡性較小。而在縫隙右側(cè)，由于氣流在后壁面處受到強(qiáng)壓縮，該區(qū)域的平動溫度TT遠(yuǎn)大于轉(zhuǎn)動溫度TR，存在強(qiáng)烈的熱不平衡效應(yīng)。而當(dāng)傾角減小時，平動溫度和轉(zhuǎn)動溫度都開始增加，但轉(zhuǎn)動溫度增大的幅度更大。當(dāng)θ=30°時，平動溫度和轉(zhuǎn)動溫度間的差距小于θ=90°時兩者間的差值。這說明傾角減小能減弱縫隙內(nèi)流場中的熱不平衡效應(yīng)。

圖7 縫隙內(nèi)密度分布Fig.7 Density ratio distribution in the gaps

5.2 表面參數(shù)分析

5.2.1表面熱流系數(shù)分析

圖11分別展示了不同前壁傾角下的縫隙的前壁面、底面、后壁面三個表面上的熱流系數(shù)分布。

從圖中可以看出，在前壁面，當(dāng)θ=90°時，熱流系數(shù)從表面頂點(diǎn)往下緩慢增大，在表面中部位置附近達(dá)到最大值。之后，熱流系數(shù)往下再緩慢下降直到表面底端。當(dāng)θ=75°和 60°時，在表面的上半部分熱流分布與θ=90°區(qū)別不大。在表面下半部分，熱流系數(shù)隨前壁傾角減小而增大。這是因?yàn)?，?dāng)前壁傾角減小時，其底部附近的流體速度也加快，傳熱速度也隨之增大。當(dāng)θ繼續(xù)減小到 45°和 30°時，前壁面上整體熱流系數(shù)有了大幅提高。熱流系數(shù)從表面頂點(diǎn)往下單調(diào)增大，并在底部位置達(dá)到最大值。

在縫隙底面上，對于不同前壁傾角，沿來流方向，表面熱流系數(shù)都是先增加后減少，因此存在一個熱流峰值，峰值對應(yīng)在底面上的x軸坐標(biāo)都是約為 0.8，該位置也是縫隙右側(cè)旋渦的起始位置。同時，可以看到，隨著前壁傾角的減小，在底面的前半部分，整體熱流系數(shù)也隨之增大。這是因?yàn)殡S著傾角的減小，流體越來越深入縫隙，造成熱流系數(shù)的整體增加。當(dāng)傾角越小，增加的速度也越大。在底面的后半部分，不同工況下的熱流系數(shù)開始趨于定值。

在后壁面上，各工況下的表面熱流系數(shù)均在頂點(diǎn)處達(dá)到最大值，且從頂點(diǎn)沿后壁面向下越來越小。同時，隨著傾角的減小，在后壁面的同一位置處，表面熱流系數(shù)雖然逐漸增大，但幅度很小。這說明改變前壁面傾角θ，并不會對右側(cè)壁面產(chǎn)生較大影響。前壁傾角改變對縫隙內(nèi)表面的影響主要集中在前壁面以及底面的前半部分。

圖8 縫隙內(nèi)壓強(qiáng)分布Fig.8 Pressure distribution in the gaps

綜合分析三個壁面上的熱流系數(shù)，熱流系數(shù)最大的區(qū)域是在縫隙的后壁面，這是因?yàn)榱黧w在縫隙后壁面產(chǎn)生壓縮。同時，隨著θ的減小，前壁面和底面的熱流系數(shù)大幅提高，這會導(dǎo)致前壁面和底面更容易燒蝕，并進(jìn)一步減小θ值。

5.2.2表面壓力系數(shù)分析

圖12分別展示了不同前壁傾角下的縫隙的前壁面、底面、后壁面三個表面上的壓力系數(shù)分布。

當(dāng)傾角θ從 90°減小到 45°時，在前壁面底部的壓力系數(shù)變化不大，但在底部上方的大部分區(qū)域內(nèi)，壓力系數(shù)值隨著前壁傾角的減小而減小。當(dāng)θ=30°時，在表面上方部分的壓力系數(shù)分布曲線和θ=45°時基本相同，但在表面底部，壓力系數(shù)遠(yuǎn)大于θ=45°時的分布。這是因?yàn)楫?dāng)θ=30°時，縫隙左側(cè)角落旋渦開始徹底消失，流體與表面完全貼合，消除了縫隙角落處旋渦的影響，使得壓力系數(shù)在底部突然增大。

在縫隙底面上，各工況的表面壓力系數(shù)都有一個沿來流方向持續(xù)增大的趨勢，其最大值所對應(yīng)的位置均在底面最右端取得。同時，底面上整體壓力系數(shù)隨著前壁傾角的減小而增大，不過增大的幅度很小。

在縫隙后壁面上，壓力系數(shù)的變化規(guī)律與熱流系數(shù)相同。對于不同的工況，表面壓力系數(shù)在后壁面的頂點(diǎn)處取得最大值，從頂點(diǎn)沿后壁面向下，表面壓力系數(shù)迅速減小。到表面中部區(qū)域時，減速放緩。

圖9 縫隙內(nèi)溫度分布Fig.9 Temperature distribution in the gaps

圖10 縫隙內(nèi)轉(zhuǎn)動溫度分布Fig.10 Rotational temperature distribution in the gaps

圖12 表面壓力系數(shù)分布Fig.12 Pressure coefficient distribution on the surface

并且在表面的底部區(qū)域，表面壓力系數(shù)基本穩(wěn)定。同時,在后壁面的同一位置處，傾角越小，對應(yīng)的表面壓力系數(shù)越大。這是因?yàn)?，隨著θ減小，縫隙內(nèi)氣流速度增大，對后壁面的壓縮效應(yīng)也隨之增強(qiáng)。

6 結(jié) 論

本文聚焦高超聲速稀薄流航天器表面的縫隙，使用直接模擬蒙特卡洛方法進(jìn)行數(shù)值模擬，對可變壁面傾角的梯形縫隙的流場情況及壁面參數(shù)進(jìn)行分析，結(jié)果表明：

1) 側(cè)壁傾斜式變形對縫隙內(nèi)部流場形成了高溫引流和加速導(dǎo)流兩種效應(yīng)。前者可以誘導(dǎo)更多高溫流體更深入縫隙內(nèi)部，流場參數(shù)整體有所提高。后者能夠加快縫內(nèi)氣體的切向流速，使其更快溢流。兩種效應(yīng)都將促使縫隙底部熱流密度升高，而這里是最接近飛行器內(nèi)部零件的位置，因此壁面傾斜形變將給飛行器縫隙內(nèi)的熱防護(hù)帶來極為嚴(yán)峻的挑戰(zhàn)。

2) 當(dāng)傾角改變時，縫隙內(nèi)部的回流區(qū)發(fā)生了明顯的變化，并以60°和30°作為分界點(diǎn)。當(dāng)θ≤60°時，縫隙內(nèi)部從一個大旋渦徹底分裂為兩個旋渦；當(dāng)θ≤30°，左側(cè)的旋渦區(qū)域完全消失，在前壁面消除了流動分離。同時，通過改變縫隙兩側(cè)壁面的角度大小，能完全消除縫隙內(nèi)的旋渦區(qū)域，改善其氣動特性。

3) 對于縫隙內(nèi)的表面而言，前壁面傾角的改變主要影響前壁面和底面前半部分。隨著傾角θ的減小，在前壁面，傳熱系數(shù)顯著增大，而壓強(qiáng)系數(shù)減小。而在底面上，傳熱系數(shù)隨θ減小而增大，壓力系數(shù)也小幅增大。在后壁面上，隨著θ減小，熱流系數(shù)和壓力系數(shù)都有增大，但幅度很小。