朱志斌, 馮 峰, 沈 清

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

引 言

邊界層轉捩會對飛行器摩擦阻力、 熱流密度和流動分離再附特征等產(chǎn)生顯著影響, 是高超聲速飛行器研制中必須關注的氣動現(xiàn)象. 但是由于影響流動轉捩的因素眾多, 且轉捩過程復雜, 對邊界層轉捩特性和機制認識還十分有限, 嚴重制約了高超聲速氣動技術的進一步發(fā)展.

巡航、 滑翔類高超聲速飛行器普遍采用扁平、 后掠的高升阻比幾何構型, 受壓力梯度和后掠角的共同作用, 飛行器表面邊界層呈現(xiàn)顯著的三維效應, 且邊界層轉捩通常受橫流速度拐點導致的無黏失穩(wěn)機制主導[1]. 這顯著區(qū)別于傳統(tǒng)平板、 圓錐等簡單幾何外形誘導的高超聲速二維邊界層流動, 二維邊界層轉捩一般由速度剖面上廣義速度拐點引起的第1模態(tài)、 第2模態(tài)等高階模態(tài)流向行波不穩(wěn)定性引起[2]. 目前, 對高超聲速邊界層橫流轉捩過程和機理的認識不足, 導致對橫流轉捩的精確預測難度較大, 也難以準確解釋橫流失穩(wěn)影響因素及作用規(guī)律[3].

美國和澳大利亞聯(lián)合開展的高超聲速國際飛行研究實驗(Hypersonic International Flight Research Experimentation, HIFiRE)項目中設計的第5次飛行試驗(HIFiRE5)專門針對高超聲速邊界層橫流轉捩問題進行研究. HIFiRE5項目以橫縱比2∶1的橢圓錐構型為飛行載荷, 開展了大量的風洞實驗、 數(shù)值分析工作,并于2016年5月成功完成飛行試驗. 飛行試驗前, Holden等[4]在LENS1風洞中, 采用全尺寸模型, 通過溫敏漆測熱技術, 觀測到了起始于橢圓錐模型中心線及側方下游的轉捩形態(tài). Berger等[5]采用全局磷光熱成像技術對38.1%縮比模型表面熱流進行了測試, 對比分析了迎角、 側滑角及Reynolds數(shù)對熱流分布的影響. Juliano等[6-7]以相同縮比HIFiRE5橢圓錐模型為對象, 在普渡大學Mach數(shù)為6的靜音風洞中(Boeing/AFOSR Mach-6 Quiet Tunnel, BAM6QT), 采用溫敏漆技術對模型表面溫度進行整體測量, 研究了來流噪聲、 壁面粗糙度、 迎角和Reynolds數(shù)對邊界層轉捩的影響. Borg等[8-10]通過壓力傳感器、 溫敏漆和油流顯示技術, 觀測到HIFiRE5縮比模型在靜音條件下的定常橫流渦和行波擾動; 而在噪聲條件下, 只通過油流圖像觀測到定常橫流渦. 飛行試驗由熱流判定的轉捩位置與前期地面試驗結果整體定性一致, 如HIFiRE5b飛行器表面邊界層出現(xiàn)觀察到3個葉片狀的轉捩線形態(tài), 轉捩起始位置分別位于飛行器中心線、 側前緣以及二者之間的中部區(qū)域[11].

HIFiRE5轉捩預測的數(shù)值方法主要基于流動穩(wěn)定性分析. 線性穩(wěn)定性理論應用于三維外形的固有缺陷, Choudhari等[12]采用基于準平行性假設和曲率效應結合的穩(wěn)定性分析方法, 對典型飛行和試驗工況下HIRiRE5邊界層轉捩特征進行了預測. Gosse等[13]采用拋物化穩(wěn)定性分析方法, 對HIFiRE5飛行彈道點下對稱面的邊界層流動進行了轉捩預測分析. 研究發(fā)現(xiàn), 側緣平面邊界層可能會出現(xiàn)第2模態(tài)的增長, 而中心線處流動受兩側旋渦影響, 呈現(xiàn)出更加復雜的不穩(wěn)定性增長模式. Li等[14]采用經(jīng)典的穩(wěn)定性分析方法, 結合風洞試驗數(shù)據(jù), 對HIFiRE5縮比模型邊界層轉捩現(xiàn)象進行了分析, 并指出需要依靠對三維邊界層轉捩現(xiàn)象的直接數(shù)值模擬來改進經(jīng)典穩(wěn)定性理論在實際應用中的基礎假設.

近期, 國外已有部分研究者采用非定常精細數(shù)值模擬方法對HIFiRE5風洞試驗模型邊界層轉捩現(xiàn)象開展研究. Dinzl等[15]對風洞試驗模型的定常橫流渦進行了直接數(shù)值模擬, 采用分布式粗糙度來激發(fā)表面最不穩(wěn)定擾動波, 并發(fā)現(xiàn)高熱流條帶是由流向速度擾動沖擊壁面引起. Tufts等[16]采用高保真非定常精細模擬方法研究了HIFiRE5風洞試驗模型邊界層轉捩現(xiàn)象, 利用平面聲波線性疊加方法引入擾動, 通過調整擾動幅值復現(xiàn)了靜音及噪聲來流條件下的模型熱流分布特征, 并獲得了壓力脈動能譜. 上述研究展現(xiàn)了精細數(shù)值模擬方法在高超聲速邊界層橫流轉捩問題研究中的良好應用潛力.

本文采用大渦模擬方法對橢圓錐邊界層橫流轉捩現(xiàn)象進行了細致模擬分析, 深入認識高超聲速邊界層橫流失穩(wěn)轉捩特征及機制, 為高超聲速飛行器邊界層轉捩準確預測及氣動設計奠定了基礎.

1 計算模型及工況

以普渡大學靜音風洞開展的橢圓錐風洞試驗[6-7]為數(shù)值研究的參考對象, 即模型為38.1%縮比的HIFiRE5外形, 模型長度328 mm, 見圖1. 風洞試驗工況參數(shù)如表1所示, 本文只考慮0°迎角. 試驗通過調節(jié)泄流槽的打開與閉合實現(xiàn)不同來流噪聲條件, 在靜來流條件湍流度低于0.05%, 而在噪聲條件下達到約3%.

圖1 橢圓錐試驗模型[6-7]Fig. 1 Experimental model of the elliptic cone[6-7]

表1 風洞試驗工況

2 流場模擬方法

2.1 控制方程及數(shù)值格式

數(shù)值計算采用了隱式大渦模擬方法[17-18], 即針對Favre濾波三維可壓縮Navier-Stokes方程[19], 基于亞格子模型對不可解尺度流動建模求解假設, 依靠數(shù)值格式的耗散特征來抑制亞格子湍流動能積聚.

為精細刻畫流場小尺度湍流結構, 同時無振蕩捕捉激波間斷, 流通量離散采用通量限制型偏心緊致格式[20], 并通過與TVD格式結合, 保證算法對激波間斷的無振蕩捕捉能力. 黏性項通過4階中心差分計算. 時間推進采用2階顯式Runge-Kutta方法. 該大渦模擬方法已多次在高超聲速復雜非定常流動中得到驗證、 應用[21-22].

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

基于橢圓錐模型橫、 縱向對稱性, 采用1/4橢圓錐模型為對象進行計算建模, 以降低計算花費. 使用多塊結構化網(wǎng)格劃分計算域, 見圖2, 流向計算域為82 mm≤x≤328 mm, 以便于在模型適當遠的位置x=82 mm施加人工入口擾動, 相同的方法已在文獻[15-16]中采用. 空間網(wǎng)格分布與脫體激波位置相匹配呈扁平型. 法向網(wǎng)格在近壁區(qū)進行加密, 壁面法向第1層網(wǎng)格尺度為5×10-4mm. 展向網(wǎng)格均勻分布, 流向網(wǎng)格尺度則逐漸增大. 流向、 法向和展向網(wǎng)格數(shù)分別為906×161×301, 對應的網(wǎng)格單元總量為4.344×107.

圖2 計算網(wǎng)格示意圖Fig. 2 Schematic of computational grids

大渦模擬以常規(guī)2階Roe格式求解得到的全模型定常層流解為初始流場進行非定常推進計算. 計算域流向入口從全模型層流流場中截取得到. 為模擬真實情況下來流擾動和壁面粗糙度等對流場的干擾作用, 在入口處引入三維速度小擾動, 其形式與三維最不穩(wěn)定T-S波相似, 曾在超聲速噴流混合流場模擬中得到應用驗證[23], 具體為

u′,v′,w′=AG(η)[cos(βθ±ωt)+ cos(0.5βθ±0.5ωt)]

其中,A為擾動幅值,G(η)為壁面距離η的函數(shù), 形式為

擾動參數(shù)ω=2πf, 擾動頻率f由當?shù)剡吔鐚油饩壦俣萓e和邊界層厚度δ確定, 即

f=Ue/(10δ)

θ為無量綱展向位置參數(shù), 范圍為0～1, 展向擾動取20個周期, 因此頻率參數(shù)β為20×2π.

法向遠場邊界設置為來流條件. 物面采用等溫無滑移邊界條件, 壁溫設為300 K. 展向對稱面采用對稱邊界條件. 出口邊界采用外插方法設置, 并對網(wǎng)格作緩沖層處理. 大渦模擬非定常計算無量綱時間步長為5×10-7, 特征長度L=1 m, 時間步長特征參考量為L/U∞, 待入口擾動波流出出口邊界2倍周期的時間即流場充分發(fā)展后作流場統(tǒng)計, 統(tǒng)計步數(shù)為1×105.

2.3 數(shù)值結果驗證

基于試驗來流噪聲條件, 入口擾動幅值分別設為0.05%和3%. 在核心計算區(qū)域內, 無量綱網(wǎng)格尺度分別為y+≈0.2,Δx+≈10,Δz+≈5, 已接近滿足DNS網(wǎng)格分辨率要求[24], 即y+<1, 10≤Δx+≤20, 5≤Δz+≤10, 表明網(wǎng)格尺度達到大渦模擬要求.

圖3給出了不同擾動幅值下數(shù)值結果與文獻[6-7]中風洞試驗結果的橢圓錐時均熱流云圖對比. 當小擾動幅值A=0.05%時, 見圖3(a), 側前緣和中心線間的中部區(qū)域出現(xiàn)兩組細長條高熱流條帶結構. 在湍流度0.05%的試驗靜來流條件下, 見圖3(c), 模型表面出現(xiàn)熱流條帶結構, 且中部區(qū)域高熱流條帶形態(tài)及位置均與數(shù)值模擬小擾動幅值A=0.05%的結果相符. 當擾動幅值為A=3%時, 見圖3(b), 時均熱流分布反映了邊界層流動的橫流轉捩形態(tài)特征, 整體模型表面出現(xiàn)多個峰狀轉捩陣面, 中心線和側前緣間的中部區(qū)域出現(xiàn)兩個轉捩陣面. 對應試驗噪聲來流條件下, 見圖3(d), 橢圓錐模型邊界層表面出現(xiàn)雙肺葉狀的轉捩形態(tài), 與計算結果對比可發(fā)現(xiàn), 大擾動幅值計算得到熱流分布反映了與試驗觀測數(shù)據(jù)相似的轉捩形態(tài)特征. 此外, 風洞試驗和時均計算結果在模型側前緣區(qū)域均沒有出現(xiàn)轉捩形態(tài)特征. 整體而言, 本文大渦模擬能夠有效預測橢圓錐邊界層失穩(wěn)轉捩現(xiàn)象, 且獲得了豐富的流場信息, 有助于獲得高超聲速邊界層流動橫流失穩(wěn)特征的機理性認識.

(a) Numerical simulation: A=0.05%

(b) Numerical simulation: A=3%

(c) Experiment: quiet inflow

(d) Experiment: noisy inflow圖3 數(shù)值結果與風洞試驗熱流數(shù)據(jù)對比Fig. 3 Heat-flux comparison between numerical results and wind tunnel experimental data

3 計算結果分析

3.1 流場結構特征

不加入口擾動計算得到的空間流場和壁面熱流分布如圖4所示. 可以看到在高超聲速來流條件下, 橢圓錐模型不同周向角位置激波強度不同, 側緣區(qū)域壓強顯著大于中心線區(qū)域, 形成的壓力梯度使得流動由再附線向中心線流動. 在模型側緣和中心線間的中部區(qū)域, 出現(xiàn)多條平行于來流方向的高熱流條帶.

圖4 無擾動流場Fig. 4 Flowfield without inlet disturbance

圖5展示了x=200 mm處流向截面流場. 可以看到流動從側緣向中心線匯聚, 使得邊界層從再附線到中心線逐漸增厚, 并在中心線附近形成蘑菇狀形態(tài)的反向旋轉旋渦結構. 在中心線兩側區(qū)域, 邊界層內、 外部展向速度和流向渦量方向相反, 表明邊界層外緣存在橫向流動剪切作用.

(a) Density

(b) X component velocity

(c) Z component velocity

(d) X component vorticity

(e) Temperature圖5 x=200 mm處流向截面流場Fig. 5 Flowfield of streamwise section at x=200 mm

圖6給出了中心線附近不同展向位置處的速度分量對比. 從中可發(fā)現(xiàn), 在0≤z≤8 mm范圍內, 流向速度分量呈S分布, 展向速度分量呈反S分布. 在z=4 mm處速度剖面曲線扭曲最為嚴重. 而當z>8 mm時, 流向速度剖面恢復至常規(guī)邊界層速度型, 展向速度剖面呈鉤狀分布, 出現(xiàn)速度拐點. 如圖7所示, 中心線區(qū)域外側的速度剖面顯示, 模型中部區(qū)域展向速度剖面仍存在速度拐點, 而模型側前緣處當?shù)卣瓜蛩俣葎t接近為0, 流向速度均保持類似二維邊界層分布形態(tài). 根據(jù)速度剖面呈現(xiàn)的流動及潛在不穩(wěn)定性特征, 可將橢圓錐邊界層流場大致分為3個不同的區(qū)域, 即中心線區(qū)域、 側緣區(qū)域以及中心線和側緣間的中部區(qū)域. 在中心線區(qū)域, 流動匯聚產(chǎn)生流向渦, 流向、 展向速度剖面均呈明顯扭曲形態(tài), 流動非常容易失穩(wěn). 在中心線與側緣之間的中部區(qū)域, 流向速度體現(xiàn)為黏性不穩(wěn)定性, 但展向速度剖面仍存在拐點, 潛在表現(xiàn)為橫流不穩(wěn)定性. 在側緣區(qū)域, 流動只表現(xiàn)為流向速度黏性不穩(wěn)定性.

(a) X component velocity

(b) Z component velocity圖6 中心線附近速度剖面Fig. 6 Velocity profile near the centerline

(a) z=25 mm

(b) y=0 mm圖7 中心線外側速度剖面Fig. 7 Velocity profile outside the centerline

3.2 擾動幅值影響

為研究橫流轉捩對擾動幅值的敏感性, 開展了不同擾動幅值下橢圓錐邊界層轉捩現(xiàn)象計算分析.圖8展示了大渦模擬得到的瞬態(tài)渦系結構, 以速度梯度第2不變量等值面展示(Q2=0.01).由圖可見, 不同擾動幅值下流場渦系結構具有明顯差異. 當擾動幅值非常小(A≤0.05%)時, 中心線區(qū)域出現(xiàn)大尺度的流向渦結構, 在中心線和側前緣間的中部區(qū)域出現(xiàn)若干條平行于來流方向的流向渦. 當擾動幅值大于0.3%后, 流場中出現(xiàn)顯著的橫流失穩(wěn)轉捩現(xiàn)象, 且擾動幅值越大, 流場渦系結構愈加豐富. 從瞬時渦系結構的分布形態(tài)看, 中心線區(qū)域渦結構尺度最大, 呈現(xiàn)出明顯的失穩(wěn)、 轉捩及湍流發(fā)展演化特征. 側前緣和中心線間的中部區(qū)域按轉捩形態(tài), 可細分為靠近中心線側和靠近外緣側兩個子區(qū)域. 在較小的擾動幅值下(A=0.3%), 兩個子區(qū)域易于區(qū)分, 靠近中心線側的渦系結構尺度相對較大、 影響區(qū)域較小, 而靠近外緣側區(qū)域的渦結構尺度小、 影響區(qū)域較大. 較大的擾動幅值下(A=3%), 中心線和中部區(qū)域間的渦系結構在展向逐步融合.

(a) A=0%

(b) A=0.05%

(c) A=0.3%

(d) A=3%圖8 不同擾動幅值下瞬態(tài)渦系結構(Q2=0.01)Fig. 8 Instantaneous vortex structures under different disturbance amplitudes (Q2=0.01)

通過不同擾動幅值下瞬時和時均熱流云圖, 可進一步分析邊界層流動失穩(wěn)轉捩特征. 見圖9, 當擾動幅值非常小時(A≤0.05%), 微小的擾動也會迅速引起中心線區(qū)域旋渦結構的發(fā)展演化, 使得壁面產(chǎn)生對應的局部高熱流分布形態(tài). 而側前緣和中心線間的中部區(qū)域出現(xiàn)兩組細長條高熱流條帶結構, 對應展向坐標范圍分別為20～25 mm和40～60 mm. 這兩類熱流條帶瞬態(tài)值和時均值沒有明顯差異, 表明對應的流向渦結構是定常的. 當擾動幅值大于0.3%時, 熱流結果反映出流場中存在顯著的非定常流動結構, 在中心線區(qū)域流場中存在大尺度渦結構, 中心線和側緣間的中部區(qū)域出現(xiàn)流向條帶分布. 時均熱流分布與小擾動幅值下定常流向渦結構位置相對應. 模型表面中心線處轉捩起始位置最為靠前. 中心線和側前緣間的中部區(qū)域出現(xiàn)兩個轉捩陣面, 其中靠近側前緣一側轉捩形態(tài)隨擾動幅值變化相對較小, 而靠近中心線一側的轉捩起始位置隨擾動幅值增大迅速前移.

瞬態(tài)熱流分布清晰地顯示了流動失穩(wěn)特征, 從中可發(fā)現(xiàn), 在時均轉捩陣面之前, 流場中存在斜向的條紋結構, 條紋結構方向與壁面摩擦力線一致, 表明對應的流動結構存在于邊界層底層, 且尺度較小. 轉捩后的熱流條帶整體形態(tài)平行于來流方向, 表明邊界層流動由外緣處的大尺度流動結構主導. 中心線和側緣之間的中部區(qū)域靠近中心線側的流動受到中心線大尺度渦和橫流擾動的共同影響, 發(fā)生失穩(wěn)轉捩位置相對更靠前, 產(chǎn)生的流動結構尺度也更大. 在側前緣區(qū)域, 不同擾動幅值下的熱流分布均沒有出現(xiàn)轉捩形態(tài)特征.

(a) Instantaneous: A=0%

(b) Time-averaged: A=0%

(c) Instantaneous: A=0.05%

(d) Time-averaged: A=0.05%

(e) Instantaneous: A=0.3%

(f) Time-averaged: A=0.3%

(g) Instantaneous: A=3%

(h) Time-averaged: A=3%

3.3 橫流轉捩特征

為認識靜來流及噪聲來流條件下高超聲速邊界層橫流轉捩特征,圖10顯示了0.05%和3%擾動幅值下流向截面的瞬時密度分布. 可以發(fā)現(xiàn)小幅值擾動下, 隨流動向下游發(fā)展, 中心線處匯聚卷起的旋渦逐漸發(fā)生變形演化, 側緣和中心線間的中部區(qū)域在靠近壁面處出現(xiàn)兩組獨立的渦結構. 靠近中心線一側的渦數(shù)量少、 尺度相對較大, 并且其形成與主渦存在關聯(lián), 靠近側緣的渦數(shù)量相對較多, 并且在靠近側緣肩部位置逐漸減弱. 大擾動幅值下, 中心線處主渦結構失穩(wěn)產(chǎn)生大尺度的流動結構, 中心線和側緣間的中部區(qū)域非定常流動結構尺度相對較小, 且更貼近壁面.

(a) A=0.05%

(b) A=3%圖10 流向截面瞬態(tài)密度Fig. 10 Instantaneous density at streamwise sections

圖11給出了0.05%和3%來流擾動幅值下不同展向位置截面瞬態(tài)密度分布.A=0.05%小幅值擾動下, 模型后部中心線區(qū)域(z=0 mm)渦外緣位置出現(xiàn)了微弱的流場擾動, 而在展向中部(z=25 mm)及靠近側緣處(z=50 mm), 沒有出現(xiàn)擾動形態(tài), 表明邊界層流動處于穩(wěn)定狀態(tài). 而A=3%擾動幅值下, 瞬態(tài)流場密度分布形態(tài), 橢圓錐邊界層流動向下游發(fā)展過程中出現(xiàn)了明顯的失穩(wěn)轉捩現(xiàn)象.

(a) A=0.05%: z=0 mm

(b) A=3%: z=0 mm

(c) A=0.05%: z=20 mm

(d) A=3%: z=20 mm

(e) A=0.05%: z=50 mm

(f) A=3%: z=50 mm圖11 展向截面瞬態(tài)密度Fig. 11 Instantaneous density at spanwise sections

圖12對比了不同擾動幅值下的密度脈動均方根和湍動能分布, 以展現(xiàn)二者湍流發(fā)展差異. 在較小的入口擾動幅值下, 只有在中心線區(qū)域流場出現(xiàn)較為微弱的流場脈動, 中心線以外流場可認為是定常的. 而在入口擾動幅值較大時, 流場非定常脈動特征顯著, 其中邊界層外緣處脈動幅值較大. 結合瞬態(tài)流場渦系結構特征, 可確定較小擾動幅值下, 中心線和側緣間的橫流渦結構是定常的, 來流擾動提高會促發(fā)定常橫流渦的二次失穩(wěn), 進而產(chǎn)生橫流失穩(wěn)轉捩現(xiàn)象.

(a) Root mean square of density fluctuation

(b) Turbulent kinetic energy圖12 不同擾動幅值密度脈動均方根、 湍動能對比Fig. 12 Comparison of root mean square of density fluctuation and turbulent kinetic energy between different disturbance amplitudes

3.4 非線性動力學分析

采用壓強頻譜分析、 相空間以及Lyapunov指數(shù)等方法對三維邊界層的失穩(wěn)轉捩非線性動力學特性進行分析.

非定常統(tǒng)計過程中在橢圓錐表面布置了12個探測點, 其位置如圖13所示. 探測點1～4位于中心對稱線, 探測點5～8和9～12分別位于不同展向位置處.

圖13 探測點位置示意圖Fig. 13 Schematic of probe positions

圖14為通過快速Fourier變換得到的壓力脈動頻譜曲線.圖14(a)顯示探測點1處脈動幅值較小, 對應于該點處于層流流動狀態(tài). 探測點2在低頻200～400 kHz范圍內存在較大脈動幅值. 而探測點3和4呈現(xiàn)出顯著的寬頻脈動特征, 頻率范圍約為200～1 000 kHz, 可推斷當?shù)亓鲃右堰_到湍流狀態(tài).圖14(b)和(c)脈動頻譜對比發(fā)現(xiàn), 在側緣處(探測點12), 壓強脈動幅值微小, 流場脈動微弱, 而在接近側緣位置(探測點8和11), 壓力開始出現(xiàn)較強震蕩, 表明擾動進入發(fā)展放大階段. 此外, 在轉捩發(fā)展的起始位置處(探測點6和7), 壓力脈動具有較低的特征頻率頻譜, 而在湍流區(qū)域(探測點5, 9, 10), 壓力在寬頻范圍內較連續(xù)光滑地過渡.

(a) Probe 1～4

(b) Probe 5～8

(c) Probe 9～12

圖15為各探測點的壓強脈動相軌跡圖, 其中橫軸坐標為壓強, 縱軸為其對時間的導數(shù). 相軌跡曲線的拓撲結構可直觀展示邊界層轉捩的非線性演化特征. 在探測點1, 壓強脈動相圖的范圍較小, 可近似認為是一個不動點, 即該點流動接近于穩(wěn)定狀態(tài). 在探測點2處, 相圖軌跡在多個環(huán)形線圈間變換, 表明流動出現(xiàn)非線性演化及分叉, 呈現(xiàn)出向混沌狀態(tài)轉變趨勢. 探測點3和4形成內外側線圈纏繞結構的相圖軌跡, 顯示了流場中存在不同尺度、 多周期的非定常運動. 探測點6的相圖軌跡形態(tài)與探測點2相似, 但壓強導數(shù)幅值相對較小, 結合轉捩演化過程結果可發(fā)現(xiàn), 二者均處于轉捩的起始區(qū)域. 探測點5和7分別處于模型中部內外側橫流轉捩發(fā)展區(qū)域, 形成圍繞多個中心的線圈纏繞相圖軌跡形態(tài). 探測點8的相圖曲線圍繞單個空心的多重線圈形態(tài), 且相圖范圍較小, 可推斷流場中存在多周期的小尺度流動結構. 探測點9和10處于全湍流區(qū)域, 其相圖軌跡也與探測點3和4相似. 探測點11靠近模型側緣, 壓強導數(shù)相對較小, 說明對應的流動結構的尺度較小. 探測點12的相圖范圍較小, 且呈現(xiàn)規(guī)則的多重線圈, 表明此側緣處流動結構還未開始非線性演化.

(a) Probe 1

(b) Probe 2

(c) Probe 3

(d) Probe 4

(e) Probe 5

(f) Probe 6

(g) Probe 7

(h) Probe 8

(i) Probe 9

(j) Probe 10

(k) Probe 11

(l) Probe 12

Lyapunov指數(shù)可定量評價動力系統(tǒng)非線性特征, 表征了系統(tǒng)在相空間中相鄰軌道間收斂或發(fā)散的平均指數(shù)律. 最大Lyapunov指數(shù)大于零時, 意味著在系統(tǒng)相空間中, 初始兩條軌線間的差別會隨著時間的演化而呈指數(shù)律的增加以致達到動力學混沌狀態(tài)[25]. 采用Wolf等[26]提出的基于軌道跟蹤的Lyapunov指數(shù)計算方法, 從各探測點的壓強時間序列提取了最大Lyapunov指數(shù). 如圖16所示, 在中心線處(探測點1～4, 實線所示), 最大Lyapunov指數(shù)(Lyapunav exponent, LE1)由接近零急劇增大, 表明中心線處流動擾動經(jīng)歷強烈的非線性增長, 流動出現(xiàn)顯著的混沌特征. 最大Lyapunov指數(shù)沿展向的分布(虛線及點劃線)顯示, 模型中部區(qū)域壓強脈動的LE1值明顯低于中心線處. 靠近中心線的內側(探測點6和9)相對較大的Lyapunov指數(shù)值表明當?shù)亓鲃哟嬖谳^強非線性作用, 而在靠近側緣的外側區(qū)域, Lyapunov指數(shù)值較小, 表明當?shù)亓鲃臃蔷€性特征較弱.

圖16 壓強時間序列最大Lyapunov指數(shù)Fig. 16 Maximal Lyapunov exponent of pressure time series

4 研究結論

本文采用隱式大渦模擬方法, 通過在計算域入口引入三維速度擾動, 對不同擾動幅值下的HIFiRE5風洞試驗橢圓錐模型的橫流失穩(wěn)轉捩現(xiàn)象開展了細致數(shù)值研究, 得到如下結論:

(1)大渦模擬精細模擬了高超聲速邊界層橫流失穩(wěn)轉捩的非定常演化過程, 預測了靜音及噪聲來流條件下邊界層轉捩的形態(tài)特征, 獲得了與試驗數(shù)據(jù)接近的熱流分布形態(tài).

(2)橢圓錐流場中心線流動匯聚形成的流向渦結構非常容易失穩(wěn), 在中心線及側緣間的中部區(qū)域存在顯著的橫流不穩(wěn)定性, 兩種轉捩機制共同影響三維邊界層轉捩過程.

(3)擾動幅值對橢圓錐三維邊界層轉捩影響顯著. 在靜來流條件下, 橫流區(qū)域出現(xiàn)兩組獨立的定常橫流渦結構; 噪聲來流條件下, 中心線主渦和中部橫流渦發(fā)生失穩(wěn)轉捩, 在模型表面形成多峰狀的轉捩陣面.

(4)探測點壓強脈動的非線性動力學分析, 展示出三維邊界層發(fā)生失穩(wěn)轉捩的非線性演化特征, 深化解釋了模型表面邊界層內周期性小擾動發(fā)展至混沌的非線性動力學演化過程.