超聲速/高超聲速邊界層轉(zhuǎn)捩后期流場(chǎng)的模態(tài)分析

余 明, 黃偉希, 許春曉

(清華大學(xué) 航天航空學(xué)院, 北京 100084)

0 引 言

高速可壓縮邊界層流動(dòng)的轉(zhuǎn)捩存在多種不穩(wěn)定機(jī)制，而且非線性作用關(guān)系復(fù)雜[1-2]，人們開(kāi)展了一系列的工作對(duì)轉(zhuǎn)捩機(jī)理及預(yù)測(cè)方法進(jìn)行研究[3-4]。對(duì)流場(chǎng)進(jìn)行模態(tài)分解可使人們抓住流場(chǎng)中重要的物理機(jī)制和流動(dòng)結(jié)構(gòu)，是轉(zhuǎn)捩機(jī)理的研究的重要手段[5]。目前常用的模態(tài)分解方法主要有兩種：本征正交分解(Proper Orthogonal Decomposition，POD)和動(dòng)力模態(tài)分解(Dynamic Mode Decomposition，DMD)。POD是基于能量的模態(tài)分解方法，而DMD是基于頻率和增長(zhǎng)率的模態(tài)分解方法[6]。

DMD方法首先由Rowley和Schmid等提出[6-7]，并應(yīng)用于圓柱尾跡、方腔流動(dòng)、射流等流動(dòng)問(wèn)題，相繼又有很多學(xué)者提出了對(duì)該方法的修正和應(yīng)用，但是在可壓縮邊界層轉(zhuǎn)捩方面的應(yīng)用較少。He等[8]用實(shí)驗(yàn)方法分析了由圓柱尾跡誘導(dǎo)的不可壓縮平板邊界層轉(zhuǎn)捩，然后用DMD進(jìn)行數(shù)據(jù)分析，發(fā)現(xiàn)第一增長(zhǎng)階段的頻率與圓柱尾渦脫落頻率一致，第二增長(zhǎng)階段存在多個(gè)頻率的擾動(dòng)。Sayadi等[9]對(duì)Ma=0.3的可壓縮平板邊界層H-型和K-型轉(zhuǎn)捩的直接數(shù)值模擬數(shù)據(jù)進(jìn)行DMD和三重分解分析，獲得了相干結(jié)構(gòu)對(duì)雷諾切應(yīng)力的貢獻(xiàn)，發(fā)現(xiàn)與發(fā)卡渦腿對(duì)應(yīng)的低頻模態(tài)對(duì)雷諾切應(yīng)力貢獻(xiàn)最大，并用低頻模態(tài)準(zhǔn)確預(yù)測(cè)了整個(gè)轉(zhuǎn)捩過(guò)程中沿壁面法向雷諾應(yīng)力的變化；Subbareddy等[10]對(duì)Ma=6.0孤立粗糙元引起的平板邊界層轉(zhuǎn)捩進(jìn)行了動(dòng)力模態(tài)分解，并提取出了能量占優(yōu)模態(tài)的高頻流場(chǎng)結(jié)構(gòu)。

本文應(yīng)用DMD方法對(duì)Ma=2.25和Ma=6平板邊界層轉(zhuǎn)捩后期(非線性發(fā)展階段)流場(chǎng)進(jìn)行分析。通過(guò)各階模態(tài)我們可以清楚地看到流場(chǎng)結(jié)構(gòu)與壁面阻力和熱流之間的關(guān)系，以及超聲速和高超聲速邊界層轉(zhuǎn)捩中流動(dòng)結(jié)構(gòu)的差異。Ma=2.25的邊界層轉(zhuǎn)捩中的流場(chǎng)結(jié)構(gòu)與亞聲速和不可壓流動(dòng)相似，而Ma=6.0的邊界層轉(zhuǎn)捩流場(chǎng)中存在多種不穩(wěn)定的模態(tài)，對(duì)轉(zhuǎn)捩流場(chǎng)結(jié)構(gòu)和轉(zhuǎn)捩機(jī)制都有很大影響。

1 物理模型和模態(tài)分析方法

1.1 物理模型和計(jì)算參數(shù)

本文中選取的物理模型為超聲速和高超聲速平板邊界層，采用直接數(shù)值模擬(Direct Numerical Simulation, DNS)方法求解完全氣體可壓縮守恒型N-S方程，以獲取后文分析使用的流場(chǎng)數(shù)據(jù)樣本。DNS的求解采用中科院力學(xué)所李新亮研究員開(kāi)發(fā)的Opencfd軟件[11-12]，并參考文獻(xiàn)[13]和文獻(xiàn)[14]選取相應(yīng)的數(shù)值方法和計(jì)算參數(shù)。計(jì)算域如圖1所示，初始流場(chǎng)為層流，入口為二維層流剖面，壁面為無(wú)滑移等溫條件，壁溫取為流場(chǎng)的恢復(fù)溫度，上邊界和出口為無(wú)反射邊界條件，出口附加網(wǎng)格間距增大以衰減湍流脈動(dòng)，避免出口對(duì)流場(chǎng)的影響；展向邊界為周期性條件。采用有限差分方法求解方程，無(wú)粘通量采用七階WENO格式，粘性通量采用六階中心差分格式，時(shí)間推進(jìn)采用三步TVD RK格式，具體算法說(shuō)明參見(jiàn)文獻(xiàn)[11]。

求解無(wú)量綱控制方程，各物理量用來(lái)流密度、速度和溫度進(jìn)行無(wú)量綱化，特征長(zhǎng)度取為1 mm。本文主要選取了兩組參數(shù)進(jìn)行計(jì)算和分析：一組是超聲速邊界層，參考文獻(xiàn)[13]中的實(shí)驗(yàn)選取計(jì)算參數(shù)，Ma=2.25，Re=25000/mm，來(lái)流溫度T∞=169.4K，壁面相對(duì)溫度Tw/T∞=1.9，流向(x)、法向(y)和展向(z)計(jì)算域范圍分別為x=101.6～406.4,y=0～5.0,z=0～9.4，網(wǎng)格數(shù)分別為2193、72和256。在入口下游給定壁面吹吸的擾動(dòng)，形式為[15]：

圖1 計(jì)算域和網(wǎng)格示意圖Fig.1 Sketch of computational domain and meshes

vb=Aaf(x)g(z)h(t)

(1)

其中Aa為擾動(dòng)幅值，其他各項(xiàng)為：

(2)

(4)

其中kz為擾動(dòng)的展向波數(shù)，zmax為展向計(jì)算域長(zhǎng)度，βm為擾動(dòng)的頻率，xb～xe為擾動(dòng)沿流向的范圍，φl(shuí)和φm為擾動(dòng)的相位。在此算例中，取各參數(shù)如下：

xb=114.3,xe=127.0,lmax=3,mmax=3；

β1=0.155,kz1=0,φ1=φ1=0,

Aa=1,A1=0.004；

β2=β3=0.309,kz2=kz3=1,φ2=φ3=0,

φ2=0,φ3=π/2,A2=A3=0.004。

吹吸擾動(dòng)的強(qiáng)度較大，轉(zhuǎn)捩過(guò)程為從非線性階段開(kāi)始的bypass轉(zhuǎn)捩。如圖2所示為壁面摩擦阻力系數(shù)(Cf)沿流向的分布，轉(zhuǎn)捩起始的位置為x=165，轉(zhuǎn)捩后湍流部分Cf的分布與文獻(xiàn)[15]中給出擬合曲線相吻合。圖3中展示了從擾動(dòng)起始到轉(zhuǎn)捩結(jié)束的流向范圍(x=114～178)的渦結(jié)構(gòu)(由速度梯度張量的第二不變量Q的等值面表示)，可以看到在轉(zhuǎn)捩開(kāi)始位置的上游為規(guī)整的流向渦，逐漸抬升并發(fā)展出碎小的渦結(jié)構(gòu)。各個(gè)主渦誘導(dǎo)的小渦在展向的范圍逐漸增加，在下游形成完全紊亂的流動(dòng)。

圖2 壁面阻力系數(shù)曲線(Ma=2.25)Fig.2 Profile of skin friction coefficient (Ma=2.25)

圖3 轉(zhuǎn)捩流場(chǎng)渦結(jié)構(gòu)(Q=0.02)Fig.3 Vortex structure on transitional flow (Q=0.02)

另一組是高超聲速邊界層，Ma=6，Re=10000/mm，來(lái)流溫度T∞=79 K，壁溫Tw/T∞=3.72。流向(x)、法向(y)和展向(z)計(jì)算域范圍分別為x=240～1600,y=0～36,z=0～30，網(wǎng)格數(shù)分別為4482、140和256。在入口下游x=260～280壁面加隨機(jī)吹吸擾動(dòng)，即式(1)中g(shù)(z)和h(t)替換為隨機(jī)函數(shù)，隨機(jī)數(shù)的幅值取為0.04。

如圖4所示為壁面阻力系數(shù)Cf沿流向的變化。在擾動(dòng)下游x=260～800的范圍內(nèi)Cf曲線與層流解一致，沒(méi)有明顯的變化，轉(zhuǎn)捩開(kāi)始的位置為x=850。圖5中展示了z=15平面上的密度和速度分布，可見(jiàn)在x<600的區(qū)域流場(chǎng)與層流沒(méi)有明顯的差別，x=600～800范圍內(nèi)可以觀察到小幅的擾動(dòng)波。速度的擾動(dòng)波只在邊界層內(nèi)比較明顯，而密度的擾動(dòng)波明顯可以分為兩個(gè)區(qū)域：一是邊界層內(nèi)貼近壁面的部分，與速度的擾動(dòng)波相似；二是邊界層外緣的部分，對(duì)應(yīng)著邊界層外向外轉(zhuǎn)播的Mach波[16]。Duan等的研究表明，在充分發(fā)展的湍流邊界層中，Ma越高邊界層外緣的湍流間歇性和邊界層內(nèi)外的膨脹和壓縮波越強(qiáng)[17]，這與密度的脈動(dòng)是密不可分的，因此后文中我們重點(diǎn)討論各階模態(tài)中流向速度和密度的擾動(dòng)形式。

圖4 壁面阻力系數(shù)曲線(Ma=6)Fig.4 Profile of skin friction (Ma=6)

圖5 密度和流向速度的瞬時(shí)場(chǎng)(z=15.0)Fig.5 Instant flow field of density andstreamwise velocity (z=15.0)

1.2 DMD原理和算法簡(jiǎn)介

動(dòng)力模態(tài)分解是通過(guò)對(duì)流場(chǎng)的數(shù)據(jù)樣本進(jìn)行分析，提取流場(chǎng)中具有不同頻率和增長(zhǎng)率的流動(dòng)結(jié)構(gòu)的方法。若已知流場(chǎng)中的(m+1)個(gè)不同時(shí)刻的數(shù)據(jù)樣本x0,x1,…,xm，為了提取流場(chǎng)中隨時(shí)間演化的主要結(jié)構(gòu)，可將數(shù)據(jù)樣本寫(xiě)成矩陣的形式：

X=[x0,x1,…,xm-1]，X′=[x1,x2,…,xm]

(5)

若樣本的時(shí)間間隔Δt很小，則可以認(rèn)為從X到X′可以經(jīng)過(guò)一個(gè)線性映射得到，即

X′=CX

(6)

求解矩陣C的特征值問(wèn)題，便可以得到不同特征值μi(i=1,2,…,m)及其對(duì)應(yīng)的特征模態(tài)。特征值μi表示模態(tài)隨時(shí)間的演化，|μi|>1表示模態(tài)是不穩(wěn)定的，|μi|=1是中性穩(wěn)定的，|μi|<1是穩(wěn)定的。取:

ωi=lg(μi)/Δt

(7)

其實(shí)部表示模態(tài)的增長(zhǎng)率，虛部表示模態(tài)的頻率。

在實(shí)際求解過(guò)程中，矩陣C是未知的，由式(5)中矩陣求偽逆得到的C通常是病態(tài)的，這給求解特征值問(wèn)題帶來(lái)很大麻煩。因此在計(jì)算時(shí)，首先對(duì)樣本矩陣進(jìn)行預(yù)處理。目前有兩種預(yù)處理方法，一是對(duì)樣本矩陣進(jìn)行QR分解[6]，二是對(duì)樣本矩陣進(jìn)行奇異值分解[18]。本文采用第二種方法，具體算法如下：

1) 將數(shù)據(jù)矩陣X作奇異值分解：X=UΣVH，其中U和V的列向量相互正交，上標(biāo)H表示矩陣的共軛轉(zhuǎn)置，Σ為以奇異值為對(duì)角線元素的對(duì)角陣；

2) 計(jì)算C矩陣在U上的投影為C′=UHX′VΣ-1；

3) 求解C′的特征值和特征向量C′W=WΛ，Λ為特征值矩陣，W為DMD模態(tài)在U上的投影，進(jìn)一步計(jì)算DMD特征模態(tài)Φ=X′VΣ-1W。

4) 按需要將特征值和特征模態(tài)進(jìn)行排序，本文采用文獻(xiàn)[19]中的Sparsity-Promoting DMD(SPDMD)的方法對(duì)模態(tài)進(jìn)行選擇，并按模態(tài)的能量大小進(jìn)行排序。SPDMD是在一般的DMD方法基礎(chǔ)上，引入了“促進(jìn)稀疏”的方法，通過(guò)罰函數(shù)來(lái)選擇“不重要的”模態(tài)：衰減很快的模態(tài)和能量較小的中性穩(wěn)定模態(tài)會(huì)被篩選并去除，保留下來(lái)的模態(tài)一般是能量較高的中性穩(wěn)定模態(tài)和增長(zhǎng)的模態(tài)。

2 Ma=2.25邊界層轉(zhuǎn)捩的模態(tài)分析

本小節(jié)應(yīng)用DMD對(duì)Ma=2.25邊界層轉(zhuǎn)捩的流場(chǎng)結(jié)構(gòu)進(jìn)行分析。首先應(yīng)用上文中的DNS數(shù)據(jù)，取201個(gè)等時(shí)間間隔的數(shù)據(jù)樣本進(jìn)行分析。在DMD分析中我們僅取從擾動(dòng)結(jié)束到轉(zhuǎn)捩開(kāi)始的流向范圍進(jìn)行分析，即流向x=127～165的范圍和法向、展向的全計(jì)算域?？紤]到吹吸擾動(dòng)的頻率和高頻波的時(shí)間分辨率，取時(shí)間間隔Δt=1.27，總時(shí)間跨度ΔT=254。模態(tài)分析采用的數(shù)據(jù)為每個(gè)網(wǎng)格點(diǎn)上三個(gè)方向的速度(u,v,w)、密度ρ和溫度T。

如圖6(a)所示為DMD特征值，由SPDMD選出的主要模態(tài)由○圈出，散點(diǎn)的顏色表示模態(tài)初始能量的大小。計(jì)算結(jié)果表明，主要模態(tài)的特征值都分布在單位圓上，表明這些模態(tài)都是中性穩(wěn)定的。進(jìn)一步在圖6(b)中展示了DMD的頻譜，橫坐標(biāo)為模態(tài)的頻率，縱坐標(biāo)為模態(tài)的初始幅值A(chǔ)(由SPDMD算法給出)，黑色實(shí)心圓為SPDMD選出的模態(tài)，這些模態(tài)在頻率上分布范圍較大，但是能量較高的模態(tài)都集中在低頻范圍，高頻模態(tài)的能量較小。我們進(jìn)一步考察低階模態(tài)擾動(dòng)能量之和，定義擾動(dòng)動(dòng)能Ek、內(nèi)能Ee和總能Et如下：

(8)

(9)

Et(x)=Ek(x)+Ee(x)

(10)

圖7 模態(tài)重構(gòu)擾動(dòng)能量(Ma=2.25)Fig.7 Mode reconstruction of disturbanceenergy (Ma=2.25)

下面我們選取了兩個(gè)典型模態(tài)進(jìn)行分析，分別為低頻模態(tài)和頻率相對(duì)較高的模態(tài)，在圖6(b)中用○圈出。為了方便表述，我們稱這兩個(gè)低頻和高頻的模態(tài)分別為Mode1和Mode2。

如圖8(a)所示為Mode1的渦結(jié)構(gòu)，可見(jiàn)這些渦結(jié)構(gòu)都是沿流向拉長(zhǎng)的，圖8(b)和8(c)中分別為該模態(tài)引起的壁面摩擦阻力系數(shù)和熱流的分布，可以看到這些渦結(jié)構(gòu)與摩擦阻力系數(shù)和熱流升高的區(qū)域存在明顯的相關(guān)性。渦結(jié)構(gòu)引起的流向速度剖面變化是壁面阻力系數(shù)變化的主要原因，因此渦結(jié)構(gòu)強(qiáng)的區(qū)域壁面阻力系數(shù)變化也大，壁面熱流與壁面阻力系數(shù)類似。由于此時(shí)流體質(zhì)點(diǎn)本身的溫度脈動(dòng)較小，對(duì)熱流的作用可以忽略不計(jì)，渦結(jié)構(gòu)對(duì)流體質(zhì)點(diǎn)物理量的輸運(yùn)是引起熱流變化的主要原因，溫度脈動(dòng)和速度脈動(dòng)具有很強(qiáng)的相關(guān)性，因此壁面阻力系數(shù)和熱流都與渦結(jié)構(gòu)有明顯的相關(guān)性。我們進(jìn)一步對(duì)其他低頻模態(tài)進(jìn)行了考察(本文中未展示)，結(jié)果表明低頻高能量的模態(tài)與Mode1類似，隨著頻率的增加，渦結(jié)構(gòu)的流向尺度逐漸變小。這些模態(tài)中結(jié)構(gòu)對(duì)壁面阻力系數(shù)和熱流在x=130到轉(zhuǎn)捩x=165都有很重要的影響。而對(duì)于如Mode2的高頻低能量的模態(tài)(如圖9所示)，渦結(jié)構(gòu)尺度小且密集，主要集中在x=155～165的范圍內(nèi)。這些渦結(jié)構(gòu)是轉(zhuǎn)捩過(guò)程中非線性效應(yīng)的衍生物，它們?cè)谵D(zhuǎn)捩前的法向位置高于低頻高能量的模態(tài)，直到轉(zhuǎn)捩開(kāi)始時(shí)對(duì)壁面阻力和熱流才有比較重要的影響。關(guān)于速度和溫度脈動(dòng)的相關(guān)性，可以在重構(gòu)流場(chǎng)的能量中得以體現(xiàn)。如圖10所示進(jìn)一步展示了流向x=150位置處速度和溫度脈動(dòng)強(qiáng)度沿法向的變化，從二者脈動(dòng)強(qiáng)度的分布上來(lái)看，速度和溫度具有明顯的相關(guān)性。

由以上分析可知，低頻高能量的模態(tài)的流動(dòng)結(jié)構(gòu)是引起轉(zhuǎn)捩的主要原因，而小尺度的結(jié)構(gòu)則是這些模態(tài)衍生的結(jié)構(gòu)。在非線性發(fā)展階段，流場(chǎng)中的主要渦結(jié)構(gòu)是流向拉長(zhǎng)的，壁面阻力和熱流的變化也與這些結(jié)構(gòu)存在直接的關(guān)系。這一結(jié)論與Sayadi等對(duì)Ma=0.3邊界層轉(zhuǎn)捩后期的分析結(jié)論是一致的[10]，可見(jiàn)在超聲速邊界層中，轉(zhuǎn)捩流場(chǎng)的結(jié)構(gòu)和物理機(jī)制與亞聲速和不可壓的流動(dòng)是類似的。

圖8 Mode1的渦結(jié)構(gòu)俯視圖、壁面阻力系數(shù)和熱流在壁面上的分布。(a)中渦結(jié)構(gòu)采用Q=0.0001的等值面顯示，采用流向速度u染色Fig.8 (a) Top view of vortex structure, (b)skin frictioncoefficient and (c)heat flux on the wall of Mode1, vortex in (a) is shown with isosurface of Q=0.0001, colored by streamwise velocity u

圖9 Mode2的渦結(jié)構(gòu)俯視圖、壁面阻力系數(shù)和熱流在壁面上的分布。(a)中渦結(jié)構(gòu)采用Q=0.0001的等值面顯示，采用流向速度u染色Fig.9 (a) Top view of vortex structure, (b)skin frictioncoefficient and (c)heat flux on the wall of Mode2, vortex in (a) is shown with isosurface of Q=0.0001, colored by streamwise velocity u

(a) 速度

(b) 溫度

3 Ma=6 邊界層轉(zhuǎn)捩的模態(tài)分析

本小節(jié)對(duì)Ma=6邊界層轉(zhuǎn)捩的流場(chǎng)進(jìn)行分析，選用的變量與Ma=2.25的DMD模態(tài)分析相同，分析的區(qū)域取為流向x=600～800的范圍，法向和展向的全計(jì)算域。Ma=6邊界層的吹吸擾動(dòng)是隨機(jī)的，因此數(shù)據(jù)樣本的時(shí)間間隔選取不能直接由某個(gè)特征頻率得到。高超聲速邊界層中多種不穩(wěn)定模態(tài)的頻率差別很大[20]，高低頻模態(tài)的頻率量級(jí)之比在103左右，受計(jì)算條件的限制，在DMD分析中很難完全地分辨出所有頻率的模態(tài)。我們采用的做法是在保證主要高頻模態(tài)分辨率的同時(shí)，盡量分辨低頻模態(tài)。經(jīng)過(guò)幾次嘗試，我們選取時(shí)間間隔Δt=1.5，總時(shí)間跨度ΔT=300，總時(shí)間樣本數(shù)為201。

如圖11(a)所示為DMD的特征值，由SPDMD選出的模態(tài)用○圈出，散點(diǎn)的顏色表示模態(tài)能量的大小。這些模態(tài)都是分布在單位圓上的，即時(shí)間上是中性穩(wěn)定的。圖11(b)中進(jìn)一步展示了DMD的頻譜。可見(jiàn)SPDMD選出的主要模態(tài)在頻率上可以分為三個(gè)分支，一是低頻的分支，集中在Im(w)=0附近，二是中等頻率的分支，集中在Im(w)=0.6附近，三是高頻的分支，集中在Im(w)=1.4附近。我們按能量大小對(duì)模態(tài)進(jìn)行排序，并取2～4階模態(tài)對(duì)流場(chǎng)的擾動(dòng)能量進(jìn)行重構(gòu)，如圖12所示。計(jì)算結(jié)果表明，在模態(tài)分析的區(qū)域內(nèi)，擾動(dòng)內(nèi)能大于擾動(dòng)動(dòng)能，用3～4階模態(tài)就可以較好地預(yù)測(cè)擾動(dòng)能量隨流向的變化。

下面我們對(duì)三個(gè)分支分別取能量最高的脈動(dòng)模態(tài)(頻率成對(duì)出現(xiàn)的模態(tài))進(jìn)行分析，在圖11(b)中用○圈出，下文中分別稱這三個(gè)模態(tài)為Mode1(低頻)，Mode2(中等頻率)和Mode3(高頻)。

如圖13所示為低頻分支Mode1的流向速度(a)和密度(b)的分布。圖中z=13.3展向截面，速度的擾動(dòng)在邊界層內(nèi)外緣都較強(qiáng)，而密度的擾動(dòng)只在邊界層外緣較強(qiáng)。我們進(jìn)一步在邊界層內(nèi)外分別取y=1.3和y=3.0兩個(gè)法向截面，在這兩個(gè)截面上速度和密度的擾動(dòng)都是流向拉長(zhǎng)的條帶，流向速度在y=1.3的擾動(dòng)強(qiáng)度大于y=3.0，而密度在y=3.0上的擾動(dòng)強(qiáng)度大于y=1.3，這表現(xiàn)了低頻模態(tài)擾動(dòng)在法向位置分布的差異。從條帶形式來(lái)看，速度和密度的擾動(dòng)分布具有比較強(qiáng)的相關(guān)性，在轉(zhuǎn)捩后期體現(xiàn)得更為明顯。這些條帶結(jié)構(gòu)對(duì)壁面阻力和熱流的影響如圖14所示?？梢钥吹皆诘皖l模態(tài)中，轉(zhuǎn)捩初期也存在很明顯的二維擾動(dòng)波，在后期流向的條帶結(jié)構(gòu)占主導(dǎo)。二維擾動(dòng)波在熱流分布中流向延伸得更長(zhǎng)，表明熱流與速度的條帶不完全對(duì)應(yīng)，表明此時(shí)熱力學(xué)量的脈動(dòng)與速度脈動(dòng)是同量級(jí)的，熱流不僅有速度脈動(dòng)輸運(yùn)的影響，還有流體質(zhì)點(diǎn)溫度脈動(dòng)的影響。

(a) DMD特征值

(b) DMD頻譜

圖12 模態(tài)重構(gòu)擾動(dòng)能量(Ma=6)Fig.12 Mode reconstruction of disturbance energy (Ma=6)

圖13 Mode1的流向速度和密度擾動(dòng)分布Fig.13 Contours of (a) streamwise velocity and (b) density and temperature of Mode1

圖14 Mode1的壁面阻力系數(shù)和熱流Fig.14 (a) Skin friction coefficient and (b) Heat flux on the wall of Mode1

如圖15所示為中等頻率模態(tài)Mode2的速度和密度分布，所取的截面與圖14一致。在此模態(tài)中，這些擾動(dòng)波是二維的，流向波長(zhǎng)遠(yuǎn)比低頻模態(tài)小，速度的擾動(dòng)在邊界層內(nèi)部較強(qiáng)，而密度的擾動(dòng)在邊界層內(nèi)外都有較強(qiáng)的脈動(dòng)，脈動(dòng)延伸到邊界層外緣并向外輻射，這些結(jié)構(gòu)產(chǎn)生的壁面阻力和熱流的分布(如圖16所示)也是呈現(xiàn)為二維擾動(dòng)波的形式。高頻模態(tài)Mode3與Mode2類似(如圖17所示)，高頻模態(tài)與中等頻率模態(tài)除了在尺度上的差異外，其他性質(zhì)都是類似的。因此以下為了討論方便，將中等頻率模態(tài)和高頻模態(tài)統(tǒng)稱為高頻模態(tài)。

圖18中給出了壁面阻力和熱流沿流向的相關(guān)系數(shù)，它們由圖14和圖16中高低頻模態(tài)阻力和熱流展向平均得到。從圖中可以定量地看到，同不可壓、亞聲速和超聲速邊界層類似，在高超聲速邊界層中，低高頻模態(tài)的壁面阻力和熱流也具有明顯的相關(guān)性。

圖15 Mode2的流向速度和密度擾動(dòng)分布Fig.15 Contours of (a) streamwise velocity and (b) density and temperature of Mode2

圖16 Mode2的壁面阻力系數(shù)和熱流Fig.16 (a) Skin friction coefficient and (b) Heat flux on the wall of Mode2

圖19(a,b)中展示了流向x=700位置用2-4階模態(tài)重構(gòu)的密度和速度脈動(dòng)均方根沿壁面法向的變化和與DNS結(jié)果的比較。密度和速度脈動(dòng)沿法向的變化與上文中的定性分析結(jié)果相一致。圖19(c)中進(jìn)一步給出了溫度脈動(dòng)沿法向的變化，可以看到溫度的脈動(dòng)在壁面附近和邊界層外緣分別存在一個(gè)峰值，前者體現(xiàn)了溫度脈動(dòng)對(duì)熱流作用的影響，后者體現(xiàn)了邊界層外緣劇烈的熱量交換。值得注意的是，與上文中Ma=2.25的結(jié)果相差不同，速度和溫度脈動(dòng)沿法向的分布相關(guān)性很小，進(jìn)一步證實(shí)了溫度脈動(dòng)不完全是由速度脈動(dòng)輸運(yùn)引起的，還包括流體質(zhì)點(diǎn)溫度脈動(dòng)的影響。

圖17 Mode3的流向速度和密度擾動(dòng)分布Fig.17 Contours of (a) streamwise velocity and (b) density and temperature of Mode3

圖18 壁面阻力和熱流的相關(guān)系數(shù)Fig.18 Correlation coefficient betweenskin friction and heat flux

以上我們考察了Ma=6邊界層轉(zhuǎn)捩后期的高低頻模態(tài)，可以看到轉(zhuǎn)捩邊界層中的擾動(dòng)分布在兩個(gè)法向位置，分別在壁面和邊界層外緣。低頻模態(tài)和高頻模態(tài)在結(jié)構(gòu)上存在明顯的差異，低頻模態(tài)的結(jié)構(gòu)主要是流向的條帶，而高頻模態(tài)的結(jié)構(gòu)則是二維擾動(dòng)波。這一現(xiàn)象在Li 等[20]研究鈍錐邊界層轉(zhuǎn)捩中就給出了類似的結(jié)論。Li等采用了濾波的方法研究了流場(chǎng)中沿流向和法向分布的幾個(gè)點(diǎn)上的速度脈動(dòng)信號(hào)，大尺度和小尺度表現(xiàn)為不同頻率的正弦擾動(dòng)，在本文中利用DMD方法得到的模態(tài)空間分布上可以看到，高低頻模態(tài)的結(jié)構(gòu)在空間上也分別對(duì)應(yīng)于不同波長(zhǎng)的正弦形式的擾動(dòng)，與Li等的研究結(jié)果是一致的。值得注意的是，在計(jì)算的擾動(dòng)區(qū)域是隨機(jī)的吹吸擾動(dòng)，可以認(rèn)為是高頻的脈動(dòng)，即在邊界層發(fā)展的初期是不存在低頻成分的，隨著擾動(dòng)向下游的發(fā)展，低頻結(jié)構(gòu)逐漸出現(xiàn)，因此低頻模態(tài)是不穩(wěn)定高頻模態(tài)誘導(dǎo)產(chǎn)生的。

(a) 密度

(b) 速度

(c) 溫度脈沖均方根

4 結(jié) 論

超聲速和高超聲速邊界層轉(zhuǎn)捩中存在多種不同的物理機(jī)制，其相互作用復(fù)雜，模態(tài)分析是研究這一問(wèn)題的重要工具。本文對(duì)Ma=2.25和Ma=6邊界層流動(dòng)轉(zhuǎn)捩后期的流場(chǎng)進(jìn)行動(dòng)力模態(tài)分解，并對(duì)模態(tài)進(jìn)行分析討論，重點(diǎn)考察了模態(tài)中流場(chǎng)的結(jié)構(gòu)和壁面阻力和熱流的關(guān)系。計(jì)算和分析結(jié)果表明，超聲速邊界層轉(zhuǎn)捩后期的流動(dòng)與不可壓縮和亞聲速流動(dòng)相似，在轉(zhuǎn)捩后期以低頻的流向渦結(jié)構(gòu)為主，對(duì)壁面阻力和熱流的變化起主要作用；高頻小尺度結(jié)構(gòu)則是低頻流向渦的衍生物，在轉(zhuǎn)捩開(kāi)始后才對(duì)壁面阻力和熱流有比較重要的作用；高超聲速邊界層轉(zhuǎn)捩中存在多種不穩(wěn)定模態(tài)，低頻模態(tài)的流場(chǎng)結(jié)構(gòu)為流向條帶，高頻的模態(tài)流場(chǎng)結(jié)構(gòu)為二維擾動(dòng)波，低頻模態(tài)是由高頻模態(tài)誘導(dǎo)產(chǎn)生的。高低頻模態(tài)之間的相互作用可以通過(guò)條件統(tǒng)計(jì)和相關(guān)性分析等手段得出，這一方面的研究正在進(jìn)行。

致謝:感謝清華信息科學(xué)與技術(shù)國(guó)家實(shí)驗(yàn)室(籌)公共平臺(tái)與技術(shù)部提供的計(jì)算機(jī)時(shí)。

參考文獻(xiàn)：

[1]Xie S F, Yang W B, Shen Q.Review of progresses in hypersonic boundary layer transition mechanism and its applications[J].Acta Aeronautica et Astronautica Sinica, 2015, 36(3): 714- 723.(in Chinese)解少飛, 楊武兵, 沈清.高超聲速邊界層轉(zhuǎn)捩機(jī)理及應(yīng)用的若干進(jìn)展回顧[J].航空學(xué)報(bào), 2015, 36(3): 714-723.

[2]Luo J S.Transition and prediction for hypersonic boundary layers[J].Acta Aeronautica et Astronautica Sinica, 2015, 36(1): 357-372.(in Chinese)羅紀(jì)生.高超聲速邊界層的轉(zhuǎn)捩及預(yù)測(cè)[J].航空學(xué)報(bào), 2015, 36(01): 357-372.

[3]Chen J Q, Tu G H, Zhang Y F, et al.Hypersnonic boundary layer transition: what we know, where shall we go[J].Acta Aerodynamica Sinica, 2017, 35(3): 311-337.(in Chinese)陳堅(jiān)強(qiáng), 涂國(guó)華, 張毅鋒, 等.高超聲速邊界層轉(zhuǎn)捩研究現(xiàn)狀與發(fā)展趨勢(shì)[J].空氣動(dòng)力學(xué)學(xué)報(bào), 2017, 35(3): 311-337.

[4]Cao W.A study of the transition prediction of hypersonic boundary layer on plane and wedge flow[J].Acta Aerodynamica Sinica, 2009, 27(5): 516-523.(in Chinese)曹偉.高超聲速邊界層的轉(zhuǎn)捩問(wèn)題[J].空氣動(dòng)力學(xué)學(xué)報(bào), 2009, 27(5): 516-523.

[5]Rowley C W, Dawson S T M.Model reduction for flow analysis and control[J].Annual Review of Fluid Mechanics, 2017, 49: 387-417.

[6]Schmid P J.Dynamic Mode decomposition of numerical and experimental data[J].Journal of Fluid Mechanics, 2010, 656: 5-28.

[8]He G, Wang J, Pan C.Initial growth of a disturbance in a boundary layer influenced by a circular cylinder wake[J].Journal

of Fluid Mechanics, 2013, 718: 116-130.

[9]Subbareddy P K, Bartkowicz M D, Candler G V.Direct numerical simulation of high-speed transition due to an isolated roughness element[J].Journal of Fluid Mechanics, 2014, 748: 848-878.

[10]Sayadi T, Schmid P J, Nichols J W, et al.Reduced-order representation of near-wall structures in the late transitional boundary layer[J].Journal of Fluid Mechanics, 2014, 748: 278-301.

[11]傅德薰, 馬延文, 李新亮.可壓縮湍流直接數(shù)值模擬[M].北京: 科學(xué)出版社, 2010.

[12]Li X L.Direct numerical simulation techniques for hypersonic turbulent flows[J].Acta Aeronautica et Astronautica Sinica, 2015, 36(1): 147-158.(in Chinese)李新亮.高超聲速湍流直接數(shù)值模擬技術(shù)[J].航空學(xué)報(bào), 2015, 36(01): 147-158.

[13]Li X L, Fu D X, Ma Y W, et al.Acoustic calculation for supersonic turbulent boundary layer flow[J].Chinese Physics Letters, 2009, 26(9): 094701.

[14]Li X L, Fu D X, Ma Y W.Direct numerical simulation of a spatially evolving supersonic turbulent boundary layer atMa= 6[J].Chinese Physics Letters, 2006, 23(6): 1519.

[15]Pirozzoli S, Grasso F, Gatski T B.Direct numerical simulation and analysis of a spatially evolving supersonic turbulent boundary layer atM=2.25[J].Physics of Fluids, 2004, 16(3): 530-545.

[16]Williams J E F, Maidanik G.The Mach wave field radiated by supersonic turbulent shear flows[J].Journal of Fluid Mechanics, 1965, 21(4): 641-657.

[17]Duan L, Beekman I, Martin M P.Direct numerical simulation of hypersonic turbulent boundary layers.Part 3.Effect of Mach number[J].Journal of Fluid Mechanics, 2011, 672: 245-267.

[18]Tu J H, Rowley C W, Luchtenburg D M, et al.On dynamic Mode decomposition: theory and applications[J].arXiv preprint arXiv: 1312.0041, 2013.

[20]Li X, Fu D, Ma Y.Direct numerical simulation of hypersonic boundary layer transition over a blunt cone with a small angle of attack[J].Physics of Fluids (1994-present), 2010, 22(2): 025105.