高超聲速雙橢球氣動熱環(huán)境預測

朱志斌，尚 慶，潘宏祿，姜寶森

(中國航天空氣動力技術研究院， 北京 100074)

高超聲速飛行器表面氣動熱環(huán)境的準確預測一直是高超聲速流動數(shù)值模擬的重點和難點。尤其是超聲速飛行器存在的激波/邊界層干擾會在局部區(qū)域產(chǎn)生嚴重的氣動熱載荷，對飛行器防熱系統(tǒng)設計帶來嚴峻的考驗。高超聲速氣動熱是由流體粘性起主導作用的物理現(xiàn)象，在計算中受網(wǎng)格分布、數(shù)值格式、邊界條件算等因素的影響，其交錯影響決定了熱流計算的復雜性[1]。因此，合理有效的數(shù)值方法對于高超聲速飛行器熱環(huán)境的準確預測具有重要意義，進而能夠為飛行器熱防護系統(tǒng)的低冗余度設計提供有力支持。

雙橢球模型是典型的高超聲速飛行器外形，具有飛機機身與座艙組合的特征，其超聲速繞流流場中存在三維弓形頭激波、二次激波以及三維分離流動。國內外對此外形開展了大量的風洞試驗研究[2-6]，測量數(shù)據(jù)包括壓力、熱流密度、紋影照片、表面油流顯示等，這些結果相互補充，用來描述其復雜流場特征。李素循[5]針對該外形，選擇四座中等尺度的高超聲速風洞，分別開展了測力與測熱實驗，提供了詳細的測點數(shù)據(jù)，為驗證計算流體力學數(shù)值模擬方法提供了依據(jù)。

高超聲速雙橢球繞流流場中包含激波與邊界層的干擾、旋渦、分離和再附運動，具有復雜的干擾特征，精確預測其流動結構及其氣動熱環(huán)境存在較大難度。國內外已有文獻采用不同數(shù)值方法對該問題進行了模擬分析。Riedelbauch等[7]采用半隱式有限差分方法對高超聲速雙橢球繞流流場進行了層流模擬，捕捉到了激波/激波干擾、激波邊界層干擾以及流動分離等現(xiàn)象，與實驗圖像一致。劉昕等[8]通過層流計算，對比了不同格式對雙橢球熱流分布的影響，認為高階格式得到的熱流密度更準確可靠。耿云飛等[9]對比評估了不同湍流模型對與高超聲速雙橢球繞流的適用性，發(fā)現(xiàn)考慮可壓縮性修正得到的壁面熱流更準確。賀立新等[10]將DG方法應用于雙橢球高超聲速粘性繞流計算，得到的熱流曲線與實驗數(shù)據(jù)分布規(guī)律一致。以上數(shù)值模擬研究均捕捉到了基本流場結構特征，側重于對比驗證數(shù)值格式、湍流模型、離散方法等因素對雙橢球熱流分布計算結果的影響，而很少關注真實物理現(xiàn)象，并且沒有針對性地考慮流動分離誘導的轉捩現(xiàn)象對氣動熱環(huán)境的影響。

本文針對高超聲速雙橢球模型，采用層流、轉捩、全湍流不同數(shù)值模擬方法進行流場氣動熱環(huán)境模擬，通過與風洞試驗數(shù)據(jù)對比，分析計算網(wǎng)格、計算格式和湍流模型等氣動熱環(huán)境預測的影響；并采用大渦模擬，進一步對激波邊界層干擾誘導以及流動分離、轉捩、再附等物理現(xiàn)象進行更細致地模擬計算，分析認識高超聲速雙橢球模型的復雜流場結構，從而為該類外形高超聲速氣動熱環(huán)境的數(shù)值模擬預測提供借鑒和參考。

1 計算模型及來流條件

計算模型及來流條件與文獻[5]中FD-14A激波風洞的測熱試驗一致。雙橢球模型幾何外形如圖1所示。

圖1 雙橢球模型示意圖

模型幾何外形表面方程用以下公式描述(單位mm):

風洞試驗試驗段馬赫數(shù)Ma=8.04，迎角0°，來流條件如表1所示。

表1 試驗來流條件

以標準球頭(R=15 mm)模型零攻角下的駐點熱流密度值作為熱流參考值Qref=568.4 kW/m2。實驗數(shù)據(jù)包括上下表面中心線測點，以及距模型前端距離78 mm和120 mm的兩個剖面。

2 數(shù)值方法

1) 層流計算方法

層流計算不考慮轉捩和湍流影響，直接求解守恒形式的三維可壓縮Navier-Stokes方程，對流通量離散格式分別采用Roe格式[11]以及前期發(fā)展的滿足幾何守恒律的WENO格式[12]。Roe格式是典型的線性化Riemann解，不僅具有很好的激波等間斷分辨率，同時具有較高的粘性分辨率。滿足幾何守恒律的WENO格式采用守恒型網(wǎng)格導數(shù)計算方法，將標準的WENO格式分解為中心差分部分和數(shù)值耗散部分，能夠消除網(wǎng)格導數(shù)計算誤差，保證來流保持性，并具有較高的精度和分辨率。

粘性項離散選取二階中心差分逼近粘性二階導數(shù)平方項和交叉項。時間推進采用LU-SGS隱式方法。

2) 雷諾平均方法

RANS(Reynolds Averaged Navier-Stokes)方法求解雷諾平均N-S方程，可獲得湍流流動的時間平均信息，并給出工程問題關心的平均流場、雷諾應力和基本氣動力熱等信息，是工程中使用最為廣泛的湍流流動計算方法。RANS方法需要引入湍流模型來封閉方程，常用的湍流模型有一方程SA模型[13]以及兩方程SST模型[14]。SA模型從經(jīng)驗與量綱出發(fā)，只需求解一個渦粘性系數(shù)滿足的輸運方程，其構造簡單、魯棒性好，并且對壁面網(wǎng)格質量依賴小。SST模型是k-ε模型與k-ω模型的混合模型，由于該模型不需要顯式的壁面衰減函數(shù)，因而適用性好。

3) 轉捩模式方法

為模擬流動轉捩現(xiàn)象，采用基于SSTk-ω模型的γ-Reθ轉捩模式[15]，該模型把經(jīng)驗關聯(lián)方法和間歇因子方法有機地結合了起來，最終目的是求解間歇函數(shù)γ(即空間某點的流態(tài)是湍流的概率)，并通過它與湍流模型的聯(lián)合來控制轉捩的發(fā)生。該模型不僅建立了間歇因子的方程，而且建立了動量厚度雷諾數(shù)的方程，并且利用已有的一些試驗數(shù)據(jù)構建了它們彼此之間的關聯(lián)還能夠反映來流湍流度的影響，已在高超聲速鈍雙楔繞流流動中進行了驗證[16]。

4) 大渦模擬方法

大渦模擬(Large Eddy Simulation，LES)求解空間Favre濾波形式的三維非定常可壓縮Navier-Stokes方程，對大尺度運動通過求解運動微分方程直接計算，小尺度運動對大尺度運動的影響通過亞格子應力來模擬。本文采用隱式大渦模擬方法[17,18]，以數(shù)值格式耗散代替亞格子模型作用。大渦模擬方法要求采用高精度高分辨率的計算格式，以精細分辨流場小尺度結構。本文采用特征通量限制型緊致格式[19]，在光滑捕捉流場間斷的同時，精細刻畫小尺度流場結構。時間推進采用二階顯式Runger-Kutta方法。

5) 計算網(wǎng)格及邊界條件

本文不考慮側滑角影響，針對半模模型進行計算，壁面采用無滑移等溫條件，壁面溫度設為300 K。

采用貼體結構網(wǎng)格對計算域進行空間離散，在相貫線附近以及近壁區(qū)域均進行了適當?shù)募用?，基準網(wǎng)格物面法向間距為0.002 6 mm，對應網(wǎng)格雷諾數(shù)為對應網(wǎng)格雷諾數(shù)為30，網(wǎng)格單元總量1 998 000。為分析計算網(wǎng)格的影響，在基準網(wǎng)格基礎上，一方面，對物面法向網(wǎng)格進行進一步加密，第一層網(wǎng)格尺度設為0.001 mm，對應網(wǎng)格雷諾數(shù)為11；另一方面，采用拼接網(wǎng)格對分離區(qū)影響區(qū)域進行加密，拼接網(wǎng)格單元總量7 626 000。計算網(wǎng)格如圖2所示。

圖2 計算網(wǎng)格示意圖

大渦模擬計算中采用拼接網(wǎng)格在流動分離影響區(qū)域進行加密，以解析流動分離誘發(fā)轉捩等復雜流動現(xiàn)象，計算網(wǎng)格如圖3所示。網(wǎng)格單元總量 84 240 000，其中加密區(qū)域網(wǎng)格量為83 016 000。

3 結果分析

3.1 層流計算對比

首先通過層流計算，對比分析計算網(wǎng)格、數(shù)值格式對雙橢球外形高超聲速氣動熱環(huán)境計算結果影響。

1) 計算網(wǎng)格影響

在基準網(wǎng)格下采用Roe格式計算得到的流場如圖4所示，可以看到，高超聲速來流在下橢球頭部產(chǎn)生弓形激波，并與上橢球激波發(fā)生干涉，兩橢球相貫線附近出現(xiàn)三維流動分離。

圖4 高超聲速雙橢球繞流流場

采用不同網(wǎng)格得到的對稱面物面熱流同實驗數(shù)據(jù)對比如圖5所示，基準網(wǎng)格和加密網(wǎng)格得到的熱流值幾乎重合，表明基準網(wǎng)格雷諾數(shù)已達到計算要求；而采用拼接網(wǎng)格得到的計算在相貫線附近熱流曲線出現(xiàn)振蕩，這是由于拼接網(wǎng)格在分離區(qū)域空間網(wǎng)格尺度小，對小尺度流動結構分辨更為精細，捕捉到了流動非定常特性。

圖5 不同計算網(wǎng)格對稱面熱流對比

從對稱面回流區(qū)流線計算結果對比(圖6)可以看到，基準網(wǎng)格和加密網(wǎng)格計算得到的分離渦形態(tài)一致，而采用拼接網(wǎng)格計算得到的渦結構更為復雜，在一定程度上分辨出了分離渦的發(fā)展演化形態(tài)。

圖6 不同網(wǎng)格對稱面流線對比

2) 數(shù)值格式影響

分別采用Roe格式以及滿足幾何守恒律的WENO格式，采用基準網(wǎng)格進行層流計算，結果對比分析如下。

從對稱面和流向位置處的熱流曲線對比(圖7)，可看到，不同精度WENO計算格式得到的熱流值與Roe格式結果略有差別，熱流分布規(guī)律一致，采用高階格式得到的計算結果與實驗數(shù)據(jù)誤差并無明顯改進。

圖8給出了Roe格式與5階WENO格式流場計算結果，對比可發(fā)現(xiàn)，不同數(shù)值格式均光滑捕捉到了下橢球頭部弓形激波及其與上橢球激波的干擾，從壁面流線可以看出激波、邊界層干擾引起的一次和二次分離。高階格式對激波的捕捉更為精細。

圖7 不同格式對稱面熱流結果對比

圖8 不同格式等Ma線，壁面摩擦線及熱流分布計算結果對比

從對稱面回流區(qū)流線對比(圖9)可以看出，高精度計算格式得到的回流區(qū)渦形態(tài)愈加復雜，渦由一個主渦和二次渦變?yōu)槿齻€渦。這是由于高階格式具有較高的精度和分辨率，對流動小尺度流動結構的分辨能力強，從而捕捉到了流動分離渦的運動演化。

圖9 不同格式對稱面流線對比

3.2 轉捩模式及湍流模型計算

分別采用SA模型和SST模型進行全湍流計算，并采用轉捩模式對分離誘導產(chǎn)生的轉捩進行預測。

從對稱面熱流曲線對比(圖10)可以看到，一方面，轉捩模式和全湍流計算得到的熱流峰值較層流明顯提高，其中轉捩模式與實驗數(shù)據(jù)吻合最好，全湍流SST模型次之；另一方面，全湍流計算的下橢球熱流值明顯高于層流和轉捩模式結果。

不同流向位置處熱流分布對比如圖11所示，從圖11中可發(fā)現(xiàn)，在上表面(0°～90°)轉捩模式計算結果介于層流和湍流之間，分布規(guī)律與實驗數(shù)據(jù)一致，在下表面(90°～180°)轉捩模式與層流熱流計算結果重合。

圖10 不同模型對稱面熱流結果對比

圖11 不同模型流向位置處熱流結果對比

圖12為轉捩模式計算得到的湍流區(qū)域示意圖(以湍流與層流粘性系數(shù)之比的等值面表示)，可以看到，流動在下表面保持層流，上表面流動在分離影響區(qū)域發(fā)生轉捩及湍流化；上橢球頂部附近由于順壓梯度，流動發(fā)生再層流化。而全湍流計算得到的湍流范圍較廣，得到的熱流密度值偏高。 從表面極限流線(圖13)和對稱面流線(圖14)對比可發(fā)現(xiàn)，層流計算結果流動在相貫線附近發(fā)生二次分離，而轉捩模式和全湍流計算結果只得到一個主渦，其中轉捩模式得到的渦的大小與層流接近，全湍流計算結果相對較小。

圖12 湍流區(qū)域示意圖

圖13 不同模型極限流線計算結果

圖14 不同模型對稱面流線

3.3 大渦模擬分析

為進一步認識高超聲速雙橢球繞流流場結構特征及熱流分布規(guī)律，采用大渦模擬方法進行模擬分析。

圖15顯示了流場瞬時渦系結構，從中可以看到，相貫線附近激波與邊界層干擾引起流動分離，誘導產(chǎn)生的湍流相干結構沿流動方向進一步發(fā)展和演化。圖16為統(tǒng)計平均熱流分布，上橢球局部區(qū)域熱流值較高；對比可發(fā)現(xiàn)，流場渦系結構與熱流分布直接相關。

圖15 瞬時渦系結構(Q準則，流向速度著色)

圖16 平均熱流分布

從熱流分布曲線對比(圖17和圖18)可以看到，大渦模擬計算得到的熱流曲線抖動明顯，表明流場非定常特征導致局部區(qū)域的熱流密度值不穩(wěn)定，熱流數(shù)值介于層流和全湍流結果之間，與實驗數(shù)據(jù)出現(xiàn)偏差的原因可能在于來流雷諾數(shù)高，采用的計算網(wǎng)格未能充分解析小尺度邊界層湍流結構。

圖17 對稱面熱流曲線對比

圖18 不同流向位置處熱流曲線對比

4 結論

1) 高超聲速雙橢球繞流流場中存在激波/邊界層干擾、流動分離、再附等非定常流動現(xiàn)象，流動分離旋渦剪切作用會誘發(fā)流動轉捩，從而改變邊界層流動狀態(tài)，并直接導致相貫線附近熱流急劇升高。

2) 層流計算無法模擬流動轉捩和湍流對氣動加熱影響，得到的熱流密度值偏低；全湍流計算得到熱流峰值與實驗數(shù)據(jù)吻合較好，但高估了湍流影響范圍；而采用關聯(lián)轉捩模型方法得到的計算結果與試驗數(shù)據(jù)符合最好，對流動分離引起的轉捩和湍流化現(xiàn)象的模擬更為準確。大渦模擬計算獲得了三維分離渦發(fā)展演化并誘發(fā)邊界層轉捩的時空發(fā)展過程，有利于揭示流動演化及其熱流分布的物理機理。

3) 基于轉捩模式的流場計算方法，能夠較為準確的模擬轉捩位置、湍流區(qū)域等流動現(xiàn)象，從而得到正確合理的氣動熱環(huán)境，同時計算量較小，在高超聲速飛行器氣動熱環(huán)境預測中具有廣泛的應用前景。