陳 曦，涂國(guó)華，萬(wàn)兵兵，*，袁先旭，陳堅(jiān)強(qiáng)，陳久芬

(1.空天飛行空氣動(dòng)力科學(xué)與技術(shù)全國(guó)重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000；2.中國(guó)空氣動(dòng)力研究與發(fā)展中心 超高速空氣動(dòng)力研究所，綿陽(yáng) 621000)

0 引 言

高超聲速邊界層轉(zhuǎn)捩會(huì)造成表面熱流和摩阻劇烈變化，成為飛行器設(shè)計(jì)需要考慮的重要問(wèn)題。針對(duì)諸如零迎角圓錐、平板外形等簡(jiǎn)單二維邊界層的轉(zhuǎn)捩現(xiàn)象，目前已有較充分的研究并取得了較大進(jìn)展[1-5]。相較而言，雖然三維邊界層轉(zhuǎn)捩更能體現(xiàn)真實(shí)飛行的流動(dòng)特征，但由于其流動(dòng)結(jié)構(gòu)復(fù)雜而研究難度大，相關(guān)研究尚處于起步階段。典型的三維邊界層轉(zhuǎn)捩研究模型包括有攻角圓錐[6-13]、美澳HiFIRE-5的橢錐[14-20]、美國(guó)空軍AFOSR的BOLT模型[21-22]以及中國(guó)空氣動(dòng)力研究與發(fā)展中心（CARDC）的升力體標(biāo)模HyTRV[23-24]等。這些模型的共同點(diǎn)是由于周向壓力梯度的存在，驅(qū)使流體從高壓區(qū)向低壓區(qū)匯聚，進(jìn)而形成3種特征流動(dòng)區(qū)域，即流向渦區(qū)（流線匯聚）、橫流區(qū)和接觸線區(qū)（流線發(fā)散）。下面重點(diǎn)介紹流向渦轉(zhuǎn)捩研究現(xiàn)狀。

流向渦轉(zhuǎn)捩研究最早在風(fēng)洞開(kāi)展[25-26]，風(fēng)洞實(shí)驗(yàn)主要關(guān)注轉(zhuǎn)捩位置的參數(shù)影響規(guī)律，少量實(shí)驗(yàn)開(kāi)展了擾動(dòng)波的定量測(cè)量。由于缺乏空間場(chǎng)測(cè)量數(shù)據(jù)，難以判斷流向渦的主導(dǎo)失穩(wěn)模態(tài)。隨著算法和計(jì)算機(jī)技術(shù)的發(fā)展，近年來(lái)涌現(xiàn)出不少涉及流向渦轉(zhuǎn)捩的數(shù)值模擬結(jié)果。數(shù)值模擬可依賴動(dòng)態(tài)模式分解（dynamic mode decomposition, DMD）技術(shù)[27]得到不同頻率的擾動(dòng)信息，通過(guò)傅里葉分析可得到流向渦轉(zhuǎn)捩中擾動(dòng)頻譜、幅值等物理量的定量演化過(guò)程[22]。相比風(fēng)洞實(shí)驗(yàn)和數(shù)值計(jì)算，穩(wěn)定性分析可以更快速、更直接獲得流向渦主導(dǎo)失穩(wěn)模態(tài)。早期的流動(dòng)穩(wěn)定性分析忽略了流向渦的周向流動(dòng)變化，而采用局部一維線性穩(wěn)定性理論（linear stability theory, LST）分析當(dāng)?shù)匚恢命c(diǎn)的流動(dòng)穩(wěn)定性特征。最早將流向渦的展向變化納入考慮的是Choudhari等的工作[15]，他們將橢錐短軸流向渦穩(wěn)定性問(wèn)題歸結(jié)為二維矩陣特征值問(wèn)題（BiGlobal），發(fā)現(xiàn)BiGlobal結(jié)果與一維分析問(wèn)題結(jié)果相差較大，于是他們得出結(jié)論：流向渦問(wèn)題需要考慮展向流場(chǎng)變化。Choudhari等求解的是時(shí)間模式BiGlobal問(wèn)題，不能描繪流向渦模態(tài)在空間上的失穩(wěn)特征。Paredes等[17]開(kāi)展了空間模式下的BiGlobal問(wèn)題研究，發(fā)現(xiàn)流向渦存在對(duì)稱和反對(duì)稱模態(tài)，在低雷諾數(shù)條件下兩者具有相似的性質(zhì)，而在高雷諾情形下對(duì)稱模態(tài)比反對(duì)稱模態(tài)更不穩(wěn)定。李曉虎等[28]根據(jù)模態(tài)形狀函數(shù)分布特征將對(duì)稱和反對(duì)稱模態(tài)進(jìn)一步細(xì)分為Y模態(tài)和Z模態(tài)，其中Y模態(tài)主要分布在流向渦肩部，而Z模態(tài)主要分布在流向渦柄部（stem region），此外Z模態(tài)頻率、相速度和增長(zhǎng)率等都較Y模態(tài)低。最近，Choudhari等[20]采用面推進(jìn)拋物化穩(wěn)定性方程（3D parabolized stability equations, PSE3D）研究了小迎角飛行工況下橢錐流向渦擾動(dòng)演化問(wèn)題，他們發(fā)現(xiàn)反對(duì)稱模態(tài)最不穩(wěn)定，其基于PSE3D（全局穩(wěn)定性理論）的eN方法得到的N值在飛行試驗(yàn)轉(zhuǎn)捩位置達(dá)到15左右。針對(duì)HyTRV升力體標(biāo)模，流向渦結(jié)構(gòu)主要存在于標(biāo)模下表面中心區(qū)和上表面肩部區(qū)域。下表面中心流向渦結(jié)構(gòu)類似于橢錐短軸上的流向渦，而上表面肩部流向渦是一種非對(duì)稱的三維渦結(jié)構(gòu)。Chen等[24]針對(duì)該非對(duì)稱流向渦開(kāi)展了穩(wěn)定性分析，發(fā)現(xiàn)在流向渦肩部外緣區(qū)的失穩(wěn)模態(tài)（稱為外模態(tài)）占主導(dǎo)。

有迎角圓錐是工程常見(jiàn)的典型標(biāo)模，但其背風(fēng)流向渦轉(zhuǎn)捩研究還相對(duì)較少。陳曦等[29]針對(duì)高超聲速常規(guī)風(fēng)洞工況下有迎角圓錐背風(fēng)流向渦失穩(wěn)開(kāi)展了穩(wěn)定性分析和直接數(shù)值模擬工作，基本流采用轉(zhuǎn)捩場(chǎng)的平均流，研究發(fā)現(xiàn)反對(duì)稱模態(tài)最不穩(wěn)定，其N值在DNS轉(zhuǎn)捩位置可達(dá)10左右。Li等[30]針對(duì)相同工況，以層流場(chǎng)開(kāi)展了穩(wěn)定性分析，進(jìn)一步明確在上游反對(duì)稱的Z模態(tài)占據(jù)主導(dǎo)地位，可能導(dǎo)致轉(zhuǎn)捩提前發(fā)生。由于該模態(tài)分布主要集中在流向渦柄部的內(nèi)剪切層區(qū)域，因此也被稱為內(nèi)模態(tài)。他們還發(fā)現(xiàn)流場(chǎng)中存在被流向渦修正的Mack（第二）模態(tài)不穩(wěn)定性。隨后他們還考察了背風(fēng)中心線小突起粗糙元對(duì)流向渦穩(wěn)定性的影響，發(fā)現(xiàn)外模態(tài)被抑制，流向渦轉(zhuǎn)捩可能被延遲[31]。Paredes等[32]研究了迎角對(duì)鈍錐背風(fēng)面流向渦的影響，發(fā)現(xiàn)隨著迎角從0°增大到5°，背風(fēng)面中心線失穩(wěn)位置逐漸前移，但由于與風(fēng)洞實(shí)驗(yàn)轉(zhuǎn)捩位置對(duì)應(yīng)的N值都低于2，他們推測(cè)模態(tài)失穩(wěn)前的非模態(tài)增長(zhǎng)可能參與了轉(zhuǎn)捩過(guò)程。最近，Zhang等[33]利用BiGlobal和PSE3D分析了飛行工況下有迎角鈍錐流向渦失穩(wěn)特性，發(fā)現(xiàn)內(nèi)卷渦是其主導(dǎo)渦結(jié)構(gòu)，導(dǎo)致內(nèi)模態(tài)具有與外模態(tài)相當(dāng)?shù)氖Х€(wěn)頻率，與前人研究中外卷渦主導(dǎo)的結(jié)果不同；他們還發(fā)現(xiàn)由于冷壁效應(yīng)，該內(nèi)模態(tài)具有聲輻射性質(zhì)，而Mack模態(tài)最不穩(wěn)定。Wang等[34]利用直接數(shù)值模擬和BiGlobal分析考察了壁溫比對(duì)鈍錐背風(fēng)面流向渦轉(zhuǎn)捩的影響，發(fā)現(xiàn)增加壁溫比雖然可以抑制流向渦不穩(wěn)定性，但卻導(dǎo)致流向渦轉(zhuǎn)捩提前，他們認(rèn)為轉(zhuǎn)捩提前的原因在于高壁溫比促進(jìn)了側(cè)面橫流不穩(wěn)定性，后者影響了流向渦轉(zhuǎn)捩進(jìn)程。

綜上所述，有迎角圓錐背風(fēng)面流向渦失穩(wěn)研究還不充分，特別是缺乏理論與實(shí)驗(yàn)的細(xì)致對(duì)比研究。本文將針對(duì)不同工況下的有迎角圓錐背風(fēng)流向渦開(kāi)展全局穩(wěn)定性分析，分析流向渦不同模態(tài)的失穩(wěn)機(jī)制，并結(jié)合風(fēng)洞實(shí)驗(yàn)分析轉(zhuǎn)捩N值。

1 風(fēng)洞實(shí)驗(yàn)條件和結(jié)果

風(fēng)洞實(shí)驗(yàn)是在CARDC口徑1 m暫沖吹吸式高超聲速風(fēng)洞中進(jìn)行，可通過(guò)選定噴管和控制總溫總壓參數(shù)獲得所需來(lái)流馬赫數(shù)和單位雷諾數(shù)，并配備快速送機(jī)機(jī)構(gòu)，便于瞬時(shí)測(cè)熱實(shí)驗(yàn)，詳細(xì)介紹可參考文獻(xiàn)[35]。實(shí)驗(yàn)?zāi)Ｐ蜑榘脲F角7°的圓錐模型，總長(zhǎng)約800 mm，其中錐頂端采用可更換結(jié)構(gòu)，頭部半徑rn=0.05、5 mm，見(jiàn)圖1。圖1還展示了迎角α0、直角坐標(biāo)系（x為軸向，y指向背風(fēng)中心線，z指向垂直于x-y平面方向）及相應(yīng)的速度分量（U, V, W）、隨體坐標(biāo)系（ ξ,η,ζ）及相應(yīng)的速度分量（Uξ,Vη,Wζ）和半錐角（ θ0）。模型共加工兩套，一套是金屬模型，用于脈動(dòng)壓力測(cè)量，采用8個(gè)高頻脈動(dòng)壓力傳感器（PCB，采樣頻率3 MHz，安裝在同一條母線上，具體位置x=125、205、285、365、445、525、605、685 mm）；一套是非金屬模型，由金屬頭部和聚四氟乙烯尾部組成，金屬頭部理論長(zhǎng)度為165 mm，非金屬段長(zhǎng)度為635 mm，非金屬段可用于紅外熱圖測(cè)量，采用紅外熱像儀系統(tǒng)。紅外測(cè)溫技術(shù)可以在不破壞模型表面的情況下獲得模型表面溫度分布，然后根據(jù)層流到湍流的明顯溫升特點(diǎn)判斷轉(zhuǎn)捩位置或陣面。本文考察的風(fēng)洞工況見(jiàn)表1。

表1 風(fēng)洞實(shí)驗(yàn)工況Table 1 Wind tunnel experimental conditions

圖1 研究模型、來(lái)流條件、直角坐標(biāo)系和隨體坐標(biāo)系等示意圖Fig.1 Sketch of the model, freestream conditions, the Cartesian coordinate system ( x, y, z )and the body-oriented coordinate system ( ξ, η, ζ )

利用紅外熱像儀測(cè)得模型表面溫度隨時(shí)間的變化，采用階躍熱流法[36]求得表面熱流分布。在階躍熱流法中，假設(shè)模型表面熱流變化為階躍形式，據(jù)此可推得熱流值Q如下：

其中， ?Tw為壁溫在 ?t時(shí)間內(nèi)的變化值， ρ1、c1和k1分別為模型材料的密度、比熱和熱傳導(dǎo)系數(shù)。本文采用流動(dòng)穩(wěn)定段 ?t=1 s的溫升數(shù)據(jù)計(jì)算熱流值。圖2給出表1中3個(gè)工況下背風(fēng)面熱流分布。熱流數(shù)值根據(jù)以往實(shí)驗(yàn)?zāi)Ｐ偷奈镄詤?shù)標(biāo)定，存在一定誤差，但熱流變化對(duì)物性參數(shù)并不敏感。若參考背風(fēng)中心線熱流分布，以熱流開(kāi)始快速抬升位置作為轉(zhuǎn)捩起始位置，以熱流峰值位置作為轉(zhuǎn)捩終止位置，則3個(gè)工況下的轉(zhuǎn)捩發(fā)生區(qū)（圖2綠色帶狀區(qū)）分別為x∈(165mm,250mm) （ 工況1）、x∈(200mm,280mm)（工況2）和x∈(200mm,310mm)（工況3）。我們注意到轉(zhuǎn)捩區(qū)長(zhǎng)度大致在100 mm左右，遠(yuǎn)大于相近工況下的數(shù)值模擬結(jié)果[34]。這說(shuō)明我們的轉(zhuǎn)捩區(qū)劃定偏保守，擾動(dòng)的非線性效應(yīng)在轉(zhuǎn)捩區(qū)前期可能并不顯著。圖中還可看到，工況1和工況2的轉(zhuǎn)捩陣面與工況3的轉(zhuǎn)捩陣面明顯不同，可能是由于流場(chǎng)結(jié)構(gòu)的差異導(dǎo)致。從第3節(jié)中可見(jiàn)，工況1和工況2的背風(fēng)面中心線附近流場(chǎng)在轉(zhuǎn)捩位置附近還未發(fā)生明顯卷曲，而工況3的流場(chǎng)已出現(xiàn)明顯的流向渦結(jié)構(gòu)。此外，還注意到工況2和工況3的轉(zhuǎn)捩陣面并不對(duì)稱，推測(cè)可能的因素有兩個(gè)：一是模型安裝存在微小的偏航角；二是表面因加工誤差存在不均勻粗糙度?？紤]到中心線附近轉(zhuǎn)捩結(jié)果的局部對(duì)稱性較好，而本文關(guān)注的流動(dòng)正好集中在背風(fēng)面中心線附近，因此該實(shí)驗(yàn)結(jié)果對(duì)本文的理論分析研究仍具有參考意義。

圖2 紅外試驗(yàn)測(cè)得的圓錐背風(fēng)面轉(zhuǎn)捩陣面（左列，虛線表示中心線）和對(duì)應(yīng)的背風(fēng)中心線熱流值Q分布（右列）（a）工況1，（b）工況2，（c）工況3Fig.2 Infrared measurement of transition fronts (left) and axial variations of the wall heat flux in the leeward center lines (right), (a) Case 1，(b) Case 2，(c) Case 3

圖3給出工況1不同流向站位處壓力傳感器測(cè)得的擾動(dòng)功率譜。在x= 125 mm出現(xiàn)200 kHz左右高頻擾動(dòng)峰值，在x= 205 mm出現(xiàn)150 kHz左右的高頻峰值，之后下游功率譜變?yōu)閷捵V，沒(méi)有明顯的擾動(dòng)峰值，表明轉(zhuǎn)捩已經(jīng)發(fā)生。

圖3 工況1的功率譜Fig.3 PSD of Case 1

2 計(jì)算方法與分析理論

2.1 數(shù)值模擬方法

本文采用高精度數(shù)值模擬方法求解笛卡爾直角坐標(biāo)系(x,y,z)下的Navier-Stokes方程，方程中空間黏性項(xiàng)采用四階中心差分格式，對(duì)流項(xiàng)在經(jīng)過(guò)Lax-Friedrichs通量分裂后采用三階WENO差分格式；時(shí)間項(xiàng)采用三階Runge-Kutta離散格式進(jìn)行推進(jìn)，推進(jìn)的時(shí)間步長(zhǎng)采用當(dāng)?shù)鼗疌FL參數(shù)，這里設(shè)定為0.8。邊界條件為：壁面采用無(wú)滑移、無(wú)滲透的等溫條件，上邊界采用來(lái)流條件，出口采用一階外推法。

計(jì)算網(wǎng)格在流向、周向和法向分別布置501、301和301個(gè)點(diǎn)，并在近壁和頭部附近加密。為了準(zhǔn)確捕捉流向渦的周向流場(chǎng)變化，網(wǎng)格在背風(fēng)流向渦附近也做了加密處理。網(wǎng)格分布特征如圖4所示。為了驗(yàn)證網(wǎng)格無(wú)關(guān)性，針對(duì)工況2，補(bǔ)充計(jì)算一套不同網(wǎng)格，即法向和展向網(wǎng)格均增加到351個(gè)點(diǎn)。圖5給出基于這兩套網(wǎng)格的基本流流向渦結(jié)構(gòu)，不同顏色曲線代表不同網(wǎng)格下的流向速度等值線，對(duì)比了x= 550 mm和x= 950 mm的流向渦結(jié)構(gòu)?？梢钥闯龆呋緹o(wú)差異，說(shuō)明現(xiàn)有網(wǎng)格分辨率足夠計(jì)算流向渦結(jié)構(gòu)。

圖4 網(wǎng)格分布特征（流向間隔5個(gè)點(diǎn)繪制，法向和展向網(wǎng)格均間隔2個(gè)點(diǎn)繪制）Fig.4 Computational mesh(For clarity, every fifth points in the axial direction and every second points in other two directions are plotted)

圖5 工況2下兩套網(wǎng)格的流向渦結(jié)構(gòu)對(duì)比（流向速度等值線，間隔為0.1）Fig.5 Grid independent test of the laminar solution (axial velocity) for Case 2

2.2 全局穩(wěn)定性理論及eN方法

全局穩(wěn)定性理論有兩種分析方法，分別是BiGlobal和PSE3D。BiGlobal是求解當(dāng)?shù)囟S特征值問(wèn)題的理論方法，相比LST方法，增加考慮了展向流動(dòng)變化，可預(yù)示更復(fù)雜的渦流動(dòng)失穩(wěn)機(jī)制。PSE3D是面推進(jìn)的三維拋物化穩(wěn)定性方程方法，解決的是擾動(dòng)空間演化問(wèn)題，比直接數(shù)值模擬方法效率高且不失精度?；谌址€(wěn)定性的eN方法的基本理論框架是通過(guò)BiGlobal獲得當(dāng)?shù)財(cái)_動(dòng)模態(tài)；然后將其作為入口條件，利用PSE3D沿流向推進(jìn)獲得不同流向站位的擾動(dòng)增長(zhǎng)率；最后通過(guò)eN方法積分增長(zhǎng)率，獲得N值分布。下面具體介紹理論方法。

BiGlobal和PSE3D均要求基本流在流向是緩變的，這點(diǎn)在隨體坐標(biāo)系下是成立的，即基本流沿ξ方向（母線方向）變化緩慢。于是，下面均基于隨體坐標(biāo)系構(gòu)建BiGlobal和PSE3D的方法體系。

對(duì)于BiGlobal方程，將流場(chǎng)分解成層流場(chǎng)加擾動(dòng)場(chǎng)，其中擾動(dòng)場(chǎng)基于平行流假設(shè)，其三維擾動(dòng)形式為：

其中，q≡(ρ,u,v,w,T)T，代表時(shí)均場(chǎng)，為以溫度脈動(dòng)峰值歸一化的擾動(dòng)形狀函數(shù)；ε表征擾動(dòng)幅值大小，一般假設(shè)ε?1；ω為擾動(dòng)頻率；待求的α實(shí)部為擾動(dòng)波數(shù)，虛部負(fù)值為增長(zhǎng)率；c.c.代表共軛。本文中，增長(zhǎng)率和位置信息一般用有量綱來(lái)表示，而其他物理量均采用對(duì)應(yīng)的來(lái)流值無(wú)量綱化。將上式代入(ξ,η,ζ)坐標(biāo)系下的Navier-Stokes方程，通過(guò)線性處理（只保留一階小量）可得空間BiGlobal方程如下：

其中，A0、A1和A2是包含層流場(chǎng)信息的算子，具體形式見(jiàn)文獻(xiàn)[24]。

針對(duì)該方程，采用四階有限差分離散。展向計(jì)算域根據(jù)流場(chǎng)特征選取為 0≤?≤23 ，其中 0代表背風(fēng)中心線，而另一個(gè)展向邊界（ ?=23）遠(yuǎn)離流向渦結(jié)構(gòu)；法向外邊界距離壁面17 mm左右，同樣遠(yuǎn)離流向渦結(jié)構(gòu)。經(jīng)過(guò)計(jì)算表明，進(jìn)一步擴(kuò)大計(jì)算域并不影響分析結(jié)果。在背風(fēng)中心線處采用對(duì)稱或反對(duì)稱邊界條件，在其他邊界處均為在計(jì)算過(guò)程中，分別采用對(duì)稱和反對(duì)稱條件來(lái)得到對(duì)稱模態(tài)和反對(duì)稱模態(tài)。對(duì)于工況1和工況2，我們發(fā)現(xiàn)反對(duì)稱模態(tài)和對(duì)稱模態(tài)性質(zhì)相近（主頻、增長(zhǎng)率、相速度和形狀函數(shù)等性質(zhì)），但對(duì)稱模態(tài)更不穩(wěn)定，因此本文展示的結(jié)果均為采用對(duì)稱條件得到的對(duì)稱模態(tài)結(jié)果；對(duì)于工況3，外模態(tài)對(duì)中心線處邊界條件并不敏感，而內(nèi)模態(tài)則在反對(duì)稱條件下更不穩(wěn)定，因此展示了對(duì)稱外模態(tài)和反對(duì)稱內(nèi)模態(tài)的結(jié)果。根據(jù)流場(chǎng)結(jié)構(gòu)和不穩(wěn)定性性質(zhì)的不同分別采用從100到201和從100到300不等的網(wǎng)格點(diǎn)。方程采用隱式Arnoldi方法[37]進(jìn)行迭代求解，當(dāng)找到一個(gè)物理不穩(wěn)定模態(tài)時(shí)，可進(jìn)一步采用牛頓迭代法計(jì)算該模態(tài)在頻率或流向位置連續(xù)變化時(shí)的特征值。牛頓迭代方法[33, 38]可寫作如下形式：

其中下標(biāo)j代表迭代步數(shù)，Ξ、F和JF可表示為：

f1是離散的全局穩(wěn)定性方程（與式（3）一致），可寫作如下矩陣形式：

f2是歸一化條件：

當(dāng)F的模（離散形式下有）小于1×10-7時(shí)，迭代停止，即可得到特征值。

對(duì)于PSE3D方程，考慮非平行效應(yīng)影響，將流場(chǎng)作如下分解：

同樣，代入Navier-Stokes方程后，基于層流場(chǎng)和擾動(dòng)形狀函數(shù)均沿流向ξ緩變，僅保留其一階流向?qū)?shù)項(xiàng)，并忽略非線性項(xiàng)O(ε2)，可得PSE3D方程，形式如下：

式中L和M為包含層流場(chǎng)信息的微分算子，具體形式見(jiàn)文獻(xiàn)[24]。

該方程的離散方式與BiGlobal方程一致，求解方法采用一階隱式歐拉格式的流向推進(jìn)方式，用來(lái)推進(jìn)的入口條件由BiGlobal給出。為保證形狀函數(shù)沿流向ξ緩變，需在每步迭代α使得：

其中：

當(dāng)前后兩步的波數(shù)差別的絕對(duì)值小于1×10-5時(shí)迭代收斂。得到α后，即可積分其虛部（負(fù)值為增長(zhǎng)率）得到N值：

其中，ξ0代表擾動(dòng)中性失穩(wěn)（即BiGlobal預(yù)測(cè)的擾動(dòng)增長(zhǎng)率為零）位置。

3 結(jié)果與分析

圖6給出3個(gè)工況下圓錐背風(fēng)區(qū)流向速度剖面云圖，可以看出圓錐背風(fēng)中心線附近都存在流向渦結(jié)構(gòu)，且其強(qiáng)度隨雷諾數(shù)或迎角的增大而增強(qiáng)。下面分別針對(duì)這3個(gè)工況的流向渦開(kāi)展全局穩(wěn)定性分析。

圖6 3個(gè)工況的背風(fēng)面流向渦結(jié)構(gòu)沿流向的演化（云圖代表流向速度，虛線為背風(fēng)中心線位置）Fig.6 Axial evolutions of leeward vortices in Case 1, Case 2 and Case 3（Dashed lines represent the leeward centerlines）

3.1 工況1

工況1為小迎角尖錐流動(dòng)，其背風(fēng)流向渦結(jié)構(gòu)并不顯著，如圖7所示。背風(fēng)對(duì)稱面附近的基本流在約401 mm之前只能看到微微隆起，推測(cè)上游主要由修正的二維邊界層不穩(wěn)定性主導(dǎo)。為了驗(yàn)證這一推測(cè)，對(duì)x= 401 mm之前的流向渦進(jìn)行全局穩(wěn)定性分析。

圖7 工況1在x = 101、401 mm流向站位處的流向速度分布Fig.7 Axial velocity profiles in y-z planes at x = 101 mm and x = 401 mm for Case 1

之前的工作以及研究[17]表明反對(duì)稱模態(tài)與對(duì)稱模態(tài)性質(zhì)相似，但后者更不穩(wěn)定。圖8為x= 101 mm、f= 190 kHz的對(duì)稱模態(tài)特征譜，可見(jiàn)特征譜上只有一個(gè)不穩(wěn)定模態(tài)。計(jì)算表明上游流場(chǎng)只存在這一個(gè)不穩(wěn)定模態(tài)，該模態(tài)隨流向站位和頻率的變化關(guān)系如圖9所示。圖9（a）顯示，模態(tài)的溫度脈動(dòng)形狀函數(shù)上下雙層結(jié)構(gòu)，且其主頻隨流向位置向下游發(fā)展而逐漸減小，屬于典型的Mack模態(tài)特征，因此可判斷該模態(tài)屬于Mack模態(tài)。圖9（a）的理論預(yù)示結(jié)果與圖3風(fēng)洞實(shí)驗(yàn)頻譜在轉(zhuǎn)捩之前隨流向站位的變化規(guī)律一致。此外，從圖9（b）可以看出，Mack模態(tài)的無(wú)量綱相速度集中在0.85～0.95，且隨頻率增加先減小后緩慢增加。這些結(jié)果與二維情形一致。在401 mm處，除了主Mack模態(tài)外（綠色五角星），還出現(xiàn)高頻不穩(wěn)定模態(tài)（藍(lán)色和淺藍(lán)色五角星）。從圖9給出的形狀函數(shù)分布看，高頻不穩(wěn)定模態(tài)同樣也屬于Mack模態(tài)。

圖8 工況1在x = 101 mm流向站位處頻率為190 kHz的特征譜Fig.8 Eigenvalue spectrum for frequency 190 kHz at x = 101 mm in Case 1

圖9 工況1在6個(gè)流向站位處的不穩(wěn)定模態(tài)分布，云圖表示x = 401 mm處的前兩階模態(tài)（最不穩(wěn)定頻率）的溫度形狀函數(shù)（實(shí)部）分布Fig.9 Variations of (a) growth rates and (b) phase velocities with frequency at six axial stations.Insets are the temperature disturbance(real part) distributions of the two most unstable modes in y-z planes

利用PSE3D方法對(duì)上述主導(dǎo)模態(tài)的增長(zhǎng)率進(jìn)行積分獲得N值分布，如圖10所示，圖中不同顏色曲線代表不同頻率的N值。可以看出，頻率越低，N值越大，在x= 500 mm處，最大N值為5左右。與圖2（a）的實(shí)驗(yàn)轉(zhuǎn)捩位置相比可知，轉(zhuǎn)捩區(qū)對(duì)應(yīng)的N值約在2～2.5之間。此外，高頻擾動(dòng)在x= 400 mm后開(kāi)始重新增長(zhǎng)，表明流向渦不穩(wěn)定性開(kāi)始出現(xiàn)。

圖10 工況1不同頻率流向渦失穩(wěn)模態(tài)的N值（綠色區(qū)域表示實(shí)驗(yàn)測(cè)量的轉(zhuǎn)捩位置范圍）Fig.10 Axial variations of N-factors of unstable modes for the leeward vortex with several frequencies in Case 1（The green zone represents the experimental measured transition zone）

3.2 工況2

相比工況1，工況2的迎角仍為2°，但模型鈍度和來(lái)流雷諾數(shù)均有增加。圖11展示了4個(gè)流向站位處的背風(fēng)中心線附近的流向速度云圖，可見(jiàn)在300 mm前流場(chǎng)都保持較好的二維特征，300 mm處流場(chǎng)在背風(fēng)中心線附近有近似均勻的隆起（即低速區(qū)），該低速區(qū)在350 mm處出現(xiàn)明顯的周向扭曲，并在400 mm處出現(xiàn)流向渦形態(tài)的卷曲。工況2的基本流在300 mm之前與工況1的基本流定性一致，之后則呈現(xiàn)更復(fù)雜的流動(dòng)特征，這是因?yàn)轭^部鈍度和來(lái)流雷諾數(shù)較工況1有顯著增加的緣故。

圖11 工況2在x = 100、300、350、400 mm流向站位處的流向速度層流場(chǎng)Fig.11 Axial velocity profiles in y-z planes in Case 2 at x = 100, 300, 350, 400 mm

基本流在x= 300 mm前后的差異直接導(dǎo)致流動(dòng)不穩(wěn)定性在300 mm前后呈現(xiàn)迥異的特征，如圖12所示。從圖12（a）可以看出，在x= 300 mm之前，Mack模態(tài)為唯一不穩(wěn)定性；隨著流向站位向下游變化，其主頻從350 kHz遞減至100 kHz，其最大增長(zhǎng)率則先增后減。這與工況1中的最大增長(zhǎng)率沿流向逐漸遞減不同。此外在上游工況2中Mack模態(tài)明顯比工況1中的穩(wěn)定，這些差異主要源自鈍度影響：更大的鈍度抑制了Mack模態(tài)的增長(zhǎng)。

圖12 不同流向站位的不穩(wěn)定模態(tài)增長(zhǎng)率隨頻率的變化：（a） x = 100 ~300 mm； （b）x = 350 mm；（c） x = 400 mm.五角形代表Mack模態(tài)，圓圈代表流向渦外模態(tài)和內(nèi)模態(tài)；（d） x = 400 mm處最不穩(wěn)外模態(tài)形狀函數(shù)（溫度脈動(dòng)實(shí)部）分布； (e) x=400 mm處最不穩(wěn)內(nèi)模態(tài)的形狀函數(shù)（溫度脈動(dòng)實(shí)部）分布Fig.12 Variations of the growth rate with frequency at (a) x = 100 ~300 mm; (b) x = 350 mm; (c) x = 400 mm.(d) and (e) are, respectively,the temperature disturbance (real part) distributions of the most unstable inner and outer modes in y-z planes at x = 400 mm(Pentagon: Mack modes; Circle: inner and outer modes of leeward vortex)

在x= 350 mm處，Mack模態(tài)仍占據(jù)主導(dǎo)地位，但還出現(xiàn)了多個(gè)流向渦模態(tài)，如圖12（b）所示。其中，在Mack模態(tài)頻帶附近出現(xiàn)兩支外模態(tài)（擾動(dòng)集中在外剪切層），而在低頻區(qū)出現(xiàn)兩支內(nèi)模態(tài)（擾動(dòng)集中在內(nèi)剪切層）。外模態(tài)的失穩(wěn)頻帶非常寬，其中一支外模態(tài)（紅色符號(hào)）在100 kHz附近還呈現(xiàn)三維Mack模態(tài)特征。隨著基本流在下游進(jìn)一步卷曲，流向渦模態(tài)在x= 400 mm成為主導(dǎo)失穩(wěn)模態(tài)，見(jiàn)圖12（c）。圖12（d、e）展示了x= 400 mm處最不穩(wěn)定流向渦外模態(tài)和內(nèi)模態(tài)的形狀函數(shù)分布，可見(jiàn)外模態(tài)分布在流向渦肩部處的外剪切層，而內(nèi)模態(tài)位于流向渦柄部?jī)?nèi)剪切層，與之前研究結(jié)果[30]一致。

由于流向渦內(nèi)外模態(tài)遠(yuǎn)在實(shí)驗(yàn)觀察到的轉(zhuǎn)捩位置之后才逐漸成為主導(dǎo)失穩(wěn)模態(tài)，因此它們不大可能是轉(zhuǎn)捩的主導(dǎo)因素。在轉(zhuǎn)捩位置之前只有Mack模態(tài)，于是在圖13只給出了Mack模態(tài)的N值變化曲線。與工況1結(jié)果相似，擾動(dòng)頻率越低，失穩(wěn)區(qū)越靠后而N最大值越大。與實(shí)驗(yàn)轉(zhuǎn)捩位置（見(jiàn)圖2（b））對(duì)應(yīng)的最不穩(wěn)定頻率約為120～160 kHz，轉(zhuǎn)捩區(qū)對(duì)應(yīng)的N值約在0.4～2.5之間。

3.3 工況3

工況3比工況2的迎角增加2°，迎角的增加導(dǎo)致周向壓力梯度增大，進(jìn)而在圓錐背風(fēng)面產(chǎn)生更顯著的流向渦結(jié)構(gòu)，如圖14所示。流向渦結(jié)構(gòu)早在x= 180 mm處就已顯現(xiàn)，在x= 250 mm處已發(fā)展成熟。

圖14 工況3在x = 120、180、250、350 mm流向站位處的流向速度層流場(chǎng)Fig.14 Axial velocity profiles in y-z planes in Case 3 at x = 120, 180, 250, 350 mm

圖15展示了6個(gè)流向站位處的不穩(wěn)定模態(tài)擾動(dòng)增長(zhǎng)率分布。可以看出，上游Mack模態(tài)依然為主導(dǎo)不穩(wěn)定性（圖15（a）），其主頻是工況2中相同位置處Mack主頻的一半左右，這是因?yàn)楣r3迎角增加導(dǎo)致背風(fēng)面邊界層厚度增大的緣故。隨著邊界層的快速卷曲，低頻流向渦內(nèi)模態(tài)（藍(lán)色和綠色圓圈）率先出現(xiàn)（圖15（b））。然后在180 mm處成為主導(dǎo)不穩(wěn)定性（圖15（c）），而Mack模態(tài)的最大增長(zhǎng)率和主頻皆沿下游單調(diào)下降，直到在200 mm之后完全消失（圖15（d））。在200 mm處，失穩(wěn)頻率范圍達(dá)20～350 kHz的外模態(tài)（紅色圓圈）出現(xiàn)，并最終在下游成為主導(dǎo)模態(tài)（圖15（e、f））。

圖15 工況3不同流向站位處的不穩(wěn)定模態(tài)擾動(dòng)增長(zhǎng)率分布（五角形代表Mack模態(tài)，圓圈代表流向渦外模態(tài)和內(nèi)模態(tài)）Fig.15 Variations of growth rate with frequency for unstable modes in Case 3 at x= 120, 140, 180, 200, 220, 250 mm(Pentagon: Mack modes; Circle: inner and outer modes of leeward vortex)

圖16給出了x= 250 mm處流向渦失穩(wěn)主導(dǎo)模態(tài)的形狀函數(shù)分布，圖16（a、b）表示圖15中前兩階最不穩(wěn)定的外模態(tài)（紅色圓圈和黑色圓圈），圖16（c、d）表示圖15中前兩階最不穩(wěn)定的內(nèi)模態(tài)（藍(lán)色圓圈和綠色圓圈）。可見(jiàn)外模態(tài)分布在流向渦肩部外緣，而內(nèi)模態(tài)分布在流向渦柄部。

圖16 工況3下x = 250 mm處前四階不穩(wěn)定模態(tài)的溫度（實(shí)部）形狀函數(shù)分布（云圖）：（a、b）前兩階外模態(tài)；（c、d）前兩階內(nèi)模態(tài).頻率為圖14（f）中各模態(tài)最不穩(wěn)定頻率，線圖為流向速度等值線圖，間隔為0.1Fig.16 Temperature disturbance (real part) distributions of the two most unstable (a, b) outer and (c, d) inner modes in y-z planes at x=250 mm overlaid by the axial velocity contour lines

圖17給出了不同頻率模態(tài)擾動(dòng)的N值演化曲線，可見(jiàn)Mack模態(tài)在上游率先失穩(wěn)，這與前兩個(gè)工況相似，但最大N值均難以超過(guò)2；而反對(duì)稱內(nèi)模態(tài)在200 mm后迅速失穩(wěn)并很快成為主導(dǎo)不穩(wěn)定性，與風(fēng)洞實(shí)驗(yàn)觀察到轉(zhuǎn)捩位置附近（圖2（c））對(duì)應(yīng)的最不穩(wěn)定頻率約為30 kHz，轉(zhuǎn)捩區(qū)對(duì)應(yīng)的N值約在1.8～6之間。

圖17 工況3不同頻率流向渦失穩(wěn)模態(tài)的N值（實(shí)線：反對(duì)稱內(nèi)模態(tài)；虛線：對(duì)稱外模態(tài)；點(diǎn)劃線：對(duì)稱Mack模態(tài).綠色區(qū)域表示實(shí)驗(yàn)測(cè)量的轉(zhuǎn)捩位置范圍）Fig.17 Axial variations of N-factors of unstable modes for leeward vortex with several frequencies in Case 3 (Solid lines:inner modes; dashed lines: outer modes; solid dotted lines: Mack modes.The green zone represents the experimental measured transition zone)

4 結(jié) 論

本文針對(duì)高超聲速有迎角圓錐的背風(fēng)流向渦，利用BiGlobal和PSE3D方法開(kāi)展了全局穩(wěn)定性與轉(zhuǎn)捩分析，并結(jié)合風(fēng)洞實(shí)驗(yàn)研究了流向渦轉(zhuǎn)捩N值，進(jìn)一步認(rèn)識(shí)了流向渦轉(zhuǎn)捩機(jī)制，主要結(jié)論如下：

1） 在流向渦較弱的階段，Mack模態(tài)為主導(dǎo)不穩(wěn)定性；當(dāng)流場(chǎng)出現(xiàn)明顯周向卷曲時(shí)，流向渦失穩(wěn)模態(tài)開(kāi)始出現(xiàn)，并隨著卷曲的加強(qiáng)而逐漸成為主導(dǎo)模態(tài)，同時(shí)Mack模態(tài)逐漸消失。流向渦模態(tài)可按擾動(dòng)分布分為外模態(tài)和內(nèi)模態(tài)，其中外模態(tài)具有較高主頻和相速度，內(nèi)模態(tài)主頻和相速度相對(duì)較低。

2） 結(jié)合風(fēng)洞實(shí)驗(yàn)測(cè)量的轉(zhuǎn)捩位置，對(duì)于2°小迎角圓錐，流向渦較弱，基于全局穩(wěn)定性得到的轉(zhuǎn)捩N值不到3；而對(duì)于4°迎角圓錐，流向渦相對(duì)較強(qiáng)，轉(zhuǎn)捩完成時(shí)的N值能達(dá)到6左右。

對(duì)于本文選擇的工況與研究模型，實(shí)驗(yàn)中轉(zhuǎn)捩發(fā)生位置對(duì)應(yīng)的模態(tài)失穩(wěn)的N值都比較小，說(shuō)明轉(zhuǎn)捩不全由模態(tài)失穩(wěn)貢獻(xiàn)完成的。模態(tài)失穩(wěn)擾動(dòng)需要通過(guò)風(fēng)洞中來(lái)流擾動(dòng)的激勵(lì)和演化而來(lái)，期間擾動(dòng)以非模態(tài)形式存在，可能經(jīng)歷非模態(tài)增長(zhǎng)。前人研究[39]表明非模態(tài)增長(zhǎng)是引起轉(zhuǎn)捩的重要機(jī)制之一，因此下一步將開(kāi)展非模態(tài)分析，以確定非模態(tài)增長(zhǎng)機(jī)制是否可能在上述工況的背風(fēng)區(qū)轉(zhuǎn)捩中扮演重要角色。

致謝：本文用到的風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)得到了中國(guó)空氣動(dòng)力研究與發(fā)展中心徐洋工程師等的幫助，在此表示感謝！