李曉虎, 張紹龍, 劉建新, 黃章 張永明,3,*

(1.天津大學(xué) 力學(xué)系, 天津 300072; 2.中國工程物理研究院 流體物理研究所, 四川 綿陽 621900； 3.天津市現(xiàn)代工程力學(xué)重點(diǎn)實(shí)驗(yàn)室, 天津 300072)

0 引 言

高超聲速扁平形升力體的氣動(dòng)力和氣動(dòng)熱問題是當(dāng)前航空航天領(lǐng)域遇到的重要問題，其短軸處的轉(zhuǎn)捩位置和熱流明顯與兩側(cè)不同[1-2]。這是由于短軸處基本流存在流向渦[3-6]，與兩側(cè)的基本流有本質(zhì)的區(qū)別，導(dǎo)致短軸處穩(wěn)定性與兩側(cè)有本質(zhì)區(qū)別，最終導(dǎo)致轉(zhuǎn)捩位置和熱流不同。為了預(yù)測短軸處的轉(zhuǎn)捩和熱流，可以橢圓錐作為模型，研究其短軸處的穩(wěn)定性。該處基本流沿展向和法向的變化同量級(jí)，而沿流向的變化要小得多，所以是二維全局穩(wěn)定性問題。

人們認(rèn)識(shí)到短軸處是二維全局穩(wěn)定性問題，經(jīng)歷了一個(gè)過程。20世紀(jì)90年代中期，人們開始注意到短軸處的不穩(wěn)定性[7]，但一直用一維剪切層不穩(wěn)定性來考慮問題[8-9]，所得計(jì)算結(jié)果始終無法得到實(shí)驗(yàn)數(shù)據(jù)[10]的支持。直至2009年，Choudhari等[11]才明確指出該處的穩(wěn)定性與兩側(cè)完全不同，需要用全局穩(wěn)定性分析研究。直到最近，Paredes和Theofilis[12-13]才在長短軸比為2∶1橢圓錐，來流馬赫數(shù)7，單位雷諾數(shù)分別為8×105/m和1.89×106/m的短軸處找到了不穩(wěn)定的外模態(tài)，其特征是特征函數(shù)峰值位置遠(yuǎn)離壁面，靠近邊界層外緣，相速度較高。除了外模態(tài)，Zhang和Luo[14]在不可壓展向局部型自由湍流帶條紋邊界層中還找到過內(nèi)模態(tài)，其特征是特征函數(shù)峰值位置靠近壁面，相速度較低。

在Paredes和Theofilis[12-13]研究的基礎(chǔ)之上，有以下幾個(gè)更深入的具體問題需要回答：(1) 是否有不穩(wěn)定的內(nèi)模態(tài)；(2)有哪些不穩(wěn)定的外模態(tài)；(3)同一流向位置，不穩(wěn)定模態(tài)隨波數(shù)(頻率)的變化規(guī)律；(4)不穩(wěn)定模態(tài)沿流向的演化規(guī)律；(5)不穩(wěn)定模態(tài)的物理本質(zhì)，是否都由無粘的拐點(diǎn)不穩(wěn)定性造成的。本文正是針對(duì)這幾個(gè)問題進(jìn)行研究。

1 數(shù)值計(jì)算方法

本文使用Arnoldi方法和Rayleigh商迭代法計(jì)算橢圓錐短軸處二維全局穩(wěn)定性問題，Arnoldi方法能夠求解指定范圍的全部特征值，能夠很方便的用來尋找特征模態(tài)的種類[14]。Rayleigh商迭代法，對(duì)于追蹤某一特定模態(tài)隨著某一參量的演化更為方便。

1.1 基本流求解方法

橢圓錐基本流是通過直接數(shù)值模擬求解三維守恒型可壓縮N-S方程得到的，控制方程如下：

(1)

式中，(x,y,z)表示空間坐標(biāo)，t表示時(shí)間變量，U=(ρ,ρu,ρv,ρw,ρe)T為守恒變量，ρ為密度，e=cvT+(u2+v2+w2)/2，cv為比定容熱容，T為溫度；u、v、w分別表示x、y、z方向的速度分量，E、F、G分別為x、y、z方向的對(duì)流項(xiàng)通量，Ev、Fv、Gv分別為x、y、z方向的粘性項(xiàng)通量。

本文采用Sutherland公式計(jì)算粘性系數(shù)及熱傳導(dǎo)系數(shù)，即：

(3)

其中，Cκ、Cμ為參考溫度，本文均取為110.4 K。在實(shí)際情況中，因?yàn)闄E圓錐模型較為復(fù)雜，使得計(jì)算網(wǎng)格具有不規(guī)則性和復(fù)雜性，在一般曲線坐標(biāo)系計(jì)算更為方便。通過坐標(biāo)變換，得到一般曲線坐標(biāo)系下的控制方程：

(4)

式中：

(5)

1.2 穩(wěn)定性分析方法

二維全局穩(wěn)定性分析方法(Bi-Global)，是一維線性穩(wěn)定性理論的擴(kuò)展，是在當(dāng)?shù)厍蠼馓卣髦祮栴}。無量綱化后的非守恒型N-S方程如下：

(6)

將瞬時(shí)物理量寫成基本流和擾動(dòng)量之和，消去基本流滿足的方程，忽略二階及以上擾動(dòng)量，得到非守恒形式的線性擾動(dòng)方程。通過鏈?zhǔn)角髮?dǎo)法則，將直角坐標(biāo)系的導(dǎo)數(shù)轉(zhuǎn)化到一般曲線坐標(biāo)系下，得到一般曲線坐標(biāo)系下的線性擾動(dòng)方程:

(7)

橢圓錐短軸附近的流向渦，在展向變化劇烈，不能像LST一樣設(shè)成行進(jìn)波形式，但流向變化依然較為緩慢，因此流向引入局部平行假設(shè)，設(shè)成行進(jìn)波：

(8)

代入上式，即可以得到二維線性穩(wěn)定性方程：

(9)

其中:

(10)

2 計(jì)算模型及程序驗(yàn)證

2.1 計(jì)算模型

計(jì)算采用的是HIFiRE-5飛行試驗(yàn)的38.1%縮比模型，模型為長短軸比2∶1的橢圓錐，短軸所在平面頭半徑1 mm，半錐角7°?？紤]出口對(duì)基本流的影響，模型計(jì)算域長度延長到380 mm。模型的示意圖如圖1所示，基本流計(jì)算的來流條件，參考了Juliano實(shí)驗(yàn)[2]的工況,如表1所示。

模型零攻角工況下流場具有對(duì)稱性，為節(jié)省計(jì)算資源，只取1/4模型進(jìn)行計(jì)算，網(wǎng)格數(shù)601×251×201，為了得到一個(gè)更好的流場，計(jì)算網(wǎng)格在模型頭部、邊界層內(nèi)和激波處適當(dāng)加密了網(wǎng)格。特征長度取1 mm，數(shù)值計(jì)算中，時(shí)間推進(jìn)采用具有TVD性質(zhì)的三階Runge-Kutta格式，對(duì)流項(xiàng)采用AUSM+格式進(jìn)行通量分裂。流場中含有激波，采用2階NND格式捕捉激波。展向采用對(duì)稱邊界條件，極軸處存在奇點(diǎn)，采用極軸邊界條件，出口采用線性外推邊界條件。外邊界采用遠(yuǎn)場邊界條件，壁面采用無滑移等溫條件。具體格式以及各控制方程具體系數(shù)矩陣可以參見張紹龍博士論文[15]。

表1 基本流計(jì)算的來流條件Table 1 Flow condition for the basic flow calculation

HIFiRE-5模型由于在長軸子午面和短軸子午面之間存在壓力梯度，導(dǎo)致邊界層內(nèi)流線的彎曲，基本流在短軸附近存在一個(gè)“蘑菇”狀的流向渦，該渦結(jié)構(gòu)沿流向變化緩慢，且緩慢抬升，沿展向和法向變化較為劇烈。圖2給出了短軸附近流向渦的穩(wěn)定性分析計(jì)算域選取范圍示意圖以及流向渦沿流向演化數(shù)值計(jì)算結(jié)果。

(a)

(b)

2.2 線性穩(wěn)定性分析程序驗(yàn)證

為了驗(yàn)證二維全局穩(wěn)定性分析的程序，與一維線性穩(wěn)定性方法Muller法(LST)和直接數(shù)值模擬計(jì)算結(jié)果進(jìn)行了對(duì)比。

2.2.1 與Muller法對(duì)比

二維全局穩(wěn)定性分析與一維線性穩(wěn)定性分析的最主要區(qū)別體現(xiàn)在考慮了基本流沿展向變化，因此擾動(dòng)在展向不能設(shè)成波動(dòng)形式。但當(dāng)基本流展向均勻的時(shí)候，二維全局穩(wěn)定性分析仍有效，且與一維線性穩(wěn)定性計(jì)算結(jié)果應(yīng)該是一致的。選取平板邊界層，具體參數(shù)見表2，全局穩(wěn)定性分析與用Muller法求解一維線性穩(wěn)定性分析的特征值結(jié)果對(duì)比如表3所示，所得特征函數(shù)如圖3所示，可見兩種方法所得結(jié)果一致，程序是可靠的。需要指出，這里我們計(jì)算了時(shí)間模式和空間模式兩種情況。

表2 驗(yàn)證Bi-Global程序所用基本流的工況Table 2 Basic flow freestream condition forverifying of Bi-Global program

(a) 時(shí)間模式

(b) 空間模式

αrαiωrωiMuller(temporal)1.601.46570.0466Bi-Global(temporal)1.601.46650.0462Muller(spatial)1.5917-0.06311.460Bi-Global(spatial)1.5907-0.06241.460

2.2.2 與直接數(shù)值模擬對(duì)比

為了驗(yàn)證基本流沿展向變化較大時(shí),全局穩(wěn)定性程序的計(jì)算是否可靠，與直接數(shù)值模擬做了對(duì)比。基于全局穩(wěn)定性分析基本流，以X*=192mm處的流場剖面，沿著流向均勻延拓，構(gòu)造了計(jì)算域?yàn)?0倍波長的平行流。將全局穩(wěn)定性分析結(jié)果作為直接數(shù)值模擬的入口條件，向下游進(jìn)行計(jì)算。圖4是用密度、速度、溫度和壓力不同參量得到的增長率，可見，經(jīng)過入口附近的調(diào)整，增長率最終都與全局穩(wěn)定性分析(空間模式)吻合。DNS所得特征值見表4，特征函數(shù)見圖5，都與全局穩(wěn)定性分析的結(jié)果吻合。

圖4 全局穩(wěn)定性與直接數(shù)值模擬得到的擾動(dòng)增長率Fig.4 Perturbation growth rate obtained by globalstability and direct numerical simulation

ωαr-αiDNS0.30.50500.0222Bi-Global0.30.50510.0219

圖5全局穩(wěn)定性分析與直接數(shù)值模擬對(duì)比：特征函數(shù)(等值線為全局穩(wěn)定性分析結(jié)果，彩色云圖為直接數(shù)值模擬結(jié)果)
Fig.5Comparisonbetweenglobalstabilityanalysisanddirectnumericalsimulation:eigenfunction

(Thecontouristheglobalstabilityanalysisresult,thecolorcloudimageisthedirectnumericalsimulationresult)

3 計(jì)算結(jié)果及討論

3.1 是否有不穩(wěn)定的內(nèi)模態(tài)

Parades和Theofilis[12-13]在短軸處找到了不穩(wěn)定的外模態(tài)。而在不可壓二維全局穩(wěn)定性分析中，既發(fā)現(xiàn)了外模態(tài)，也發(fā)現(xiàn)了有不穩(wěn)定的內(nèi)模態(tài)，且內(nèi)模態(tài)的增長率較大[14,17]。因此，本文認(rèn)為有必要在短軸處觀察是否有內(nèi)模態(tài)。

在X*=192mm處做全局穩(wěn)定分析，所得特征值分布如圖6所示。

圖6 全局穩(wěn)定性分析結(jié)果：特征值分布Fig.6 Globalstability analysis result: the distrbution of eigenfuction

圖中有兩個(gè)分支，都是外模態(tài)，沒有內(nèi)模態(tài)的結(jié)果。圖6中兩個(gè)分支上的特征函數(shù)如圖7所示，都是典型的外模態(tài)。其它流向位置結(jié)果與之類似，不再贅述。

3.2 不穩(wěn)定的外模態(tài)種類

在X*=192 mm處做全局穩(wěn)定分析，圖6中有兩個(gè)分支，都是外模態(tài)，分別為Y模態(tài)和Z模態(tài)(Y、Z模態(tài)見文獻(xiàn)[15-16])，其特征函數(shù)如圖7所示。Y模態(tài)特征函數(shù)在蘑菇狀慢速區(qū)的“兩肩”，沿展向伸展成條帶狀，表明該模態(tài)主要由法向剪切造成。Z模態(tài)特征函數(shù)在蘑菇狀慢速區(qū)的“腰部”，為沿法向伸展成條帶狀，表明該模態(tài)主要由展向剪切造成。Y模態(tài)又分為兩種，分別為對(duì)稱和反對(duì)稱模態(tài)，特征函數(shù)的實(shí)部和虛部如圖8所示。Z模態(tài)也分為對(duì)稱和反對(duì)稱模態(tài)兩種，特征函數(shù)如圖9所示。

3.3 在流向同一個(gè)位置，不穩(wěn)定模態(tài)隨波數(shù)的變化規(guī)律

在X*=192 mm處做全局穩(wěn)定分析，改變波數(shù)α，觀察不穩(wěn)定模態(tài)隨其變化，結(jié)果如圖10所示，圖中畫出了三個(gè)最不穩(wěn)定的Y模態(tài)和三個(gè)最不穩(wěn)定的Z模態(tài)。可見，Y模態(tài)的最大增長率高于Z模態(tài)，但是在低波數(shù)(低頻率)段，Z模態(tài)增長率略高于Y模態(tài)。我們在X*=175mm處也做了計(jì)算，結(jié)果與之類似，如圖11所示。

圖7 全局穩(wěn)定性分析結(jié)果：特征函數(shù)(黑實(shí)線表示邊界層外緣)Fig.7 Globalstability analysis result:eigenfunction (the black solid lines refer to the outer edge of the boundary layer)

圖8 對(duì)稱和反對(duì)稱的Y模態(tài)特征函數(shù)：左中(右)Fig.8 Eigenfunction of symmetric and anti-symmetric

圖9 對(duì)稱和反對(duì)稱的Z模態(tài)特征函數(shù)：左中(右)Fig.9 Eigenfunction of symmetric and anti-symmetric

圖10 X*=192mm處增長率隨波數(shù)α和頻率ωr的變化Fig.10 Growth rate varies with wave-number and frequency at X*=192mm

3.4 在固定波數(shù)下，不穩(wěn)定模態(tài)沿流向的演化規(guī)律

為研究模態(tài)沿流向的演化，取固定波數(shù)α=0.5，觀察不穩(wěn)定模態(tài)沿流向的變化，結(jié)果如圖12所示，圖中畫出了三個(gè)最不穩(wěn)定的Y模態(tài)和三個(gè)最不穩(wěn)定的Z模態(tài)?？梢?，幾個(gè)模態(tài)一直在相互競爭，沒有一個(gè)模態(tài)始終是最不穩(wěn)定的。從相速度隨流向演化的曲線來看，Y模態(tài)相速度較高，主要分布在0.9左右，Z模態(tài)相速度與Y模態(tài)有明顯差別，比其小很多。值得注意的是，在流向X*=210 mm附近，有很多不穩(wěn)定的模態(tài)在此站位交匯，很可能會(huì)引起模態(tài)間的轉(zhuǎn)換，具體轉(zhuǎn)換機(jī)理還得做深入的研究。我們又取波數(shù)α=0.7做了計(jì)算，結(jié)果與之類似，如圖13所示。

圖11 X*=175 mm處增長率隨波數(shù)α和頻率ωr的變化Fig.11 Growth rate varies with wave-numberand frequency at X*=175 mm

圖12 不穩(wěn)定模態(tài)沿流向的變化(α=0.5): 增長率ωi和相速度CFig.12 Variation of unstable modes along streamwise(α=0.5): the growth rateωi and phase velocity C

圖13 不穩(wěn)定模態(tài)沿流向的變化(α=0.7)： 增長率ωi和相速度CFig.13 Variation of unstable modes along streamwise (α=0.7): the growth rate ωi and phase velocity C

3.5 不穩(wěn)定模態(tài)的物理本質(zhì)

Y模態(tài)和Z模態(tài)的特征函數(shù)在形式上有所區(qū)別，但都是由于基本流的剪切造成。由一維基本流的穩(wěn)定性理論可知，無粘的拐點(diǎn)不穩(wěn)定性需要滿足Fjortoft拐點(diǎn)條件。橢圓錐短軸處基本流流向速度分布是二維的，在Y*-Z*平面內(nèi)既沿Y*方向有變化，又沿Z*方向有變化。而變化最劇烈的方向，即為基本流流向速度在Y*-Z*平面內(nèi)的梯度方向，基本流沿這個(gè)方向變化最劇烈，即剪切最強(qiáng)。所用Fjortoft拐點(diǎn)條件為：

(12)

其中?/?l為當(dāng)?shù)鼗玖髟赮*-Z*平面內(nèi)沿流向速度梯度方向的導(dǎo)數(shù)。對(duì)X*=192mm處的基本流，將Fjortoft拐點(diǎn)位置用綠色圓圈標(biāo)記在圖14中，相應(yīng)的Y模態(tài)和Z模態(tài)的流向速度特征函數(shù)也畫在圖中?？梢?，特征函數(shù)集中的位置與Fjortoft拐點(diǎn)吻合。這表明它們都具有相同的產(chǎn)生機(jī)理，即由基本流剪切的無粘拐點(diǎn)不穩(wěn)定性造成。

由圖14還可以看出，“蘑菇”所在區(qū)域沒有符合內(nèi)模態(tài)特性的Fjortoft拐點(diǎn)，這就解釋了為什么我們找不到內(nèi)模態(tài)。圖4中紅色云圖表示梯度方向二階導(dǎo)數(shù)大于零；藍(lán)色表示小于零；白色表示等于零；綠色圓圈表示Fjortoft拐點(diǎn)位置；黑色等值線為特征函數(shù)。

(a) Y-mode

(b) Z-mode

4 結(jié) 論

找到了不穩(wěn)定的外模態(tài)，尚未發(fā)現(xiàn)不穩(wěn)定的內(nèi)模態(tài)。外模態(tài)有兩個(gè)分支，即Y模態(tài)和Z模態(tài)，二者都存在對(duì)稱和反對(duì)稱兩種模態(tài)。在同一流向位置，Y模態(tài)最大增長率高于Z模態(tài)，但在低波數(shù)段，Z模態(tài)增長率略高于Y模態(tài)。觀察不穩(wěn)定模態(tài)沿流向的演化，發(fā)現(xiàn)幾個(gè)模態(tài)一直在相互競爭，沒有一個(gè)模態(tài)始終是最不穩(wěn)定的。所找到的不穩(wěn)定模態(tài)都出現(xiàn)在基本流沿其梯度方向的Fjortoft拐點(diǎn)處，這表明它們都具有相同的產(chǎn)生機(jī)理，即由基本流剪切的無粘拐點(diǎn)不穩(wěn)定性造成。因?yàn)闆]有符合內(nèi)模態(tài)特性的Fjortoft拐點(diǎn)，所以沒有內(nèi)模態(tài)。

致謝:感謝國家重點(diǎn)研發(fā)計(jì)劃“大科學(xué)裝置前沿研究”重點(diǎn)專項(xiàng)“高超聲速邊界層轉(zhuǎn)捩機(jī)理、預(yù)測及控制方法研究”(2016YFA0401200)提供資助。

參考文獻(xiàn)：

[1]Wheaton M B, Juliano J T, Berridge C D, et al.Instability and transition measurements in the Mach-6 Quiet Tunnel[R].AIAA 2009-3559.

[2]Juliano J Thomas.Instability and transition on the HIFiRE-5 in a Mach-6 quiet tunnel[D].Purdue University.Indiana, US, 2010.

[3]Schmisseur J, Schneider S, Collicott S.Receptivity of the Mach-4 boundary-layer on an elliptic cone to laser-generated localized free stream perturbations[R].AIAA-98-0532, 1998.

[4]Schmisseur J, Schneider S, Collicott S.Response of the Mach-4 boundary layer on an elliptic cone to laser-generated free stream perturbations[R].AIAA-99-0410, 1999.

[5]Poggie J, Kimmel R.Traveling instabilities in elliptic cone boundary-layer transition atMa=8[R].AIAA-98-0435, 1998.

[6]Huntley M, Smits A.Transition studies on an elliptic cone in Mach 8 flow using filtered Rayleigh scattering[J].European Journal of Mechanics B-Fluids, 2000, 19(5): 695-706.

[7]Lyttle I, Reed H.Use of transition correlations for three-dimensional boundary layers within hypersonic flows[R].AIAA-95-2293, 1995.

[8]Kimmel R, Klein M, Schwoerke S.Three-dimensional hypersonic laminar boundary-layer computations for transition experiment design[J].Spacecraft Rockets, 1997, 34(4): 409-415.

[9]Gosse R, Kimmel R.CFD study of three-dimensional hypersonic laminar boundary layer transition on a Mach 8 ellipticcone[R].AIAA 2009-4053, 2009.

[10]Poggie J, Kimmel R, Schwoerke S.Traveling instability waves in a Mach 8 ow over an elliptic cone[J].AIAA J, 2000, 38(2): 251-258.

[11]Choudhari M, Chang C, Jentink T, et al.Transition analysis for the HIFiRE-5 vehicle[R].AIAA 2009-4056, 2009.

[12]Paredes P, Theofilis V.Spatial linear global instability analysis of the HIFiRE-5 elliptic cone model flow[R].AIAA 2013-2880.

[13]Paredes P, Theofilis V.Traveling global instabilities on the HIFiRE-5 elliptic cone model flow[R].AIAA 2014-0075, 2014.

[14]Zhang Y M, Luo J S.Application of Arnoldi method to boundary layer instability[J].Chinese Physics B, 2015, 24(12): 124701.

[15]Zhang S L.The instability and wave propagation in the hypersonic 2∶1 elliptic cone boundary layer[D].Tianjin University, 2016.(in Chinese)張紹龍.高超聲速2∶1橢圓錐邊界層的穩(wěn)定性特征及擾動(dòng)演化[D].天津大學(xué), 2016.

[16]Li F, Choudhari M M, Duan L, et al.Nonlinear development and secondary instability of traveling crossflow vortices[J].Phys Fluids, 2014, 26: 064104.

[17]Zhang Yongming, Zaki Tamer, Sherwin Spencer, et al.Nonlinear response of a laminar boundary layer to isotropic and spanwise localized free-stream turbulence.AIAA 2011-32[R].Reston: AIAA, 2011.