基于模態分解的軸對稱超聲速射流嘯聲產生位置數值分析1)

李虎 羅勇 韓帥斌 王益民 武從海 劉旭亮

(中國空氣動力研究與發展中心,空氣動力學國家重點實驗室,四川綿陽 621000)

(中國空氣動力研究與發展中心,計算空氣動力研究所,四川綿陽 621000)

引言

超聲速射流在航空航天飛行器的推進系統中廣泛存在和應用.由于噴管構型的限制,大多數超聲速射流是不完全膨脹的,其勢核中的激波柵格結構與剪切層中的湍流擬序結構的相互作用會產生強烈的激波噪聲.嘯聲是激波噪聲頻譜上的離散頻率聲模態,主要向上游傳播,聲壓級可高達160 dB.高強度的嘯聲不僅會造成嚴重的環境噪聲污染,還會使飛行器控制舵面、尾噴管結構部件產生聲疲勞、甚至斷裂,危及飛行安全,是發展軍用/民用超聲速飛行器急需解決的重要問題之一.

Powell[1]首次發現射流嘯聲,并將其產生機制解釋為非線性的自激聲反饋環機制.自此先驅性工作開始,普遍認為嘯聲反饋環包含4 種關鍵機制[2-4]:剪切層中的擬序結構的增長、激波柵格結構和擬序結構的相互作用、聲反饋和噴管唇口附近的感受性過程.此外,Powell[1]還發現嘯聲頻率在隨完全膨脹射流馬赫數(Mj)變化時存在模態跳躍行為.在圓射流中,Powell 確認了A、B、C、D 4 種嘯聲模態.Merle[5-6]發現A 模態可分為A1和A2兩個亞模態.Davies 和Oldfield[7]和Powell等[8]發現圓射流中的嘯聲模態與射流的內稟不穩定性直接相關,并且確認A1和A2模態是軸對稱模態,B 和D 模態是擺動模態,C 模態是螺旋模態.Ponton 和Seiner[9]發現擺動模態是由兩個同量級的反向螺旋模態疊加而成.Clem等[10]在完全膨脹射流馬赫數為1.0 ＜Mj＜1.7 的圓形收縮噴管實驗還捕捉到了嘯聲B′模態、E模態和F 模態.對于低馬赫數超聲速射流,軸對稱模態是主導嘯聲模態[11-12].Ponton等[13-14]在聲速噴管欠膨脹射流實驗中記錄了嘯聲A0模態,Li 和Gao[15]在軸對稱URANS 模擬中也捕捉到了該模態.Gao 和Li[16]根據他們提出的新的嘯聲波長(頻率)預測公式,提出Ponton等[13-14]在實驗中測量到的兩種A0模態可以歸類為A1和A2模態,給出了嘯聲A0模態的一種新解釋.

理解嘯聲的產生機制,需要對嘯聲反饋環的各個環節進行深入研究.作為環節之一,射流剪切層內擬序結構與勢核內激波柵格結構的相互作用涉及嘯聲聲波的產生過程和產生位置.為了解釋激波/擬序結構相互作用如何產生嘯聲聲波,基于模型問題,即剪切層與單個孤立的壓縮波-膨脹波結構單元的相互作用,Manning[17]提出了“激波泄漏(shock leakage)”機制.Suzuki 和Lele[18]進一步發現局部渦量對激波起到屏障作用,激波會在剪切層內旋渦之間的鞍點也就是局部渦量較弱的區域泄漏,也就是說,激波泄漏和關聯嘯聲的產生只能發生在射流剪切層大尺度湍流結構充分發展并且具有明晰的鞍點的區域.自此之后,激波泄漏機制受到了越來越多的關注,并逐漸被數值[19]和實驗[20-21]所證實.

激波泄漏機制雖然指出了剪切層中的旋渦鞍點是影響射流嘯聲產生關鍵結構,但其還有不足:未能給出每種模態嘯聲的精確產生位置.為了定量的理解嘯聲反饋環,精準定位每種模態嘯聲的有效聲源位置至關重要.Powell等[8]基于紋影發現嘯聲A1、A2和B 模態的聲波是在噴管下游5 倍直徑處的區域輻射,而嘯聲C 模態的聲波輻射中心位于噴管下游六倍直徑的區域.Umeda 和Ishii[22]利用瞬時紋影發現螺旋C 模態的主控聲源位于第3 個激波柵格結構的尾緣.Edginton-Mitchell等[23]采用高分辨率平面PIV 也研究了螺旋C 模態,認為嘯聲聲波可能是在第2、第3 和第4 個激波反射點處,且擬序渦量經歷最大脈動時產生.Gao 和Li[16]基于數值分析認為前5 個激波柵格是嘯聲A1,A2,B 和C 模態的有效聲源區,并且主控聲波是在第2 和第4 個激波柵格之間產生.Panda[24]則通過嘯聲同相壓力脈動的近場映射發現嘯聲A2模態的聲波是從第3 和第4 個激波尖梢之間某個位置輻射.Shariff 和Manning[25]發現與A1模態相關的激波泄漏和噪聲輻射是在第3 和第4 個激波柵格的激波尖梢處.Edgington-Mitchell等[21]發現在第3 和第4 個激波柵格結構的激波尖梢處會產生與A1模態相關的噪聲輻射.基于近場聲測量和時間分辨的紋影成像,Mercier等[26-27]認為軸對稱嘯聲A1模態和A2模態都是在第4 個激波柵格的尾緣(激波尖梢)處產生,而B 模態是在第3 或第4 個激波尖梢處產生,具體取決于馬赫數.Li等[28-29]通過數值模擬研究發現嘯聲螺旋C 模態的有效聲源位于第4 個激波柵格,A1和A2模態的有效聲源位于第4 個激波柵格下游,B 模態的有效聲源位于第3 個激波柵格下游.

精確定位嘯聲的聲源位置是定量理解嘯聲反饋環和預測嘯聲頻率的一個關鍵所在[26].然而,截止目前,關于每種嘯聲模態究竟是單聲源還是多聲源,以及各模態嘯聲具體的聲源位置還存在一定的爭議.而且,同一類型嘯聲模態的產生位置還可能與具體的射流馬赫數相關.現有研究涉及的射流馬赫數范圍還比較有限,需要針對不同的射流馬赫數開展研究.本文采用高精度數值方法直接求解軸對稱可壓縮Navier-Stokes 方程,數值模擬了射流馬赫數為1.10 和1.15 的圓形聲速噴管欠膨脹超聲速冷射流,得到了A1和A2兩種軸對稱模態嘯聲.基于精細的非定常流場數據,通過傅里葉模態分解、本征模態分解和動態模態分解,研究了嘯聲關聯流動結構的空間演化,定位了兩種軸對稱模態嘯聲的有效聲源位置.

1 物理模型和數值模擬方法

1.1 控制方程

對于低射流馬赫數欠膨脹超聲速射流,嘯聲的主導模態是軸對稱模態,在產生嘯聲的流動范圍內,剪切層內的大尺度不穩定波或擬序結構也呈軸對稱形態,三維效應并不明顯.因而,控制方程采用軸對稱可壓縮Navier-Stokes (N-S)方程,該方程是從柱坐標系下的三維可壓縮N-S 方程退化而來.

無量綱的軸對稱可壓縮N-S 方程具體形式如下

是守恒變量,F和G是x和r方向的無黏通量

Fv和Gv是x和r方向的黏性通量,即

S是源項,即

E是單位體積總能,即

黏性應力項為

x和r方向的熱流率為

完全氣體狀態方程和聲速表達式如下

黏性系數根據Sutherland 公式計算

表1 欠膨脹超聲速射流的關鍵流動參數Table 1 The key flow parameters in underexpanded supersonic jet

1.2 數值方法

控制方程的空間對流項離散采用五階精度的有限差分加權本質無振蕩(WENO-JS-5)格式[30],迎風處理采用Van Leer 流通量分裂方法[31];此外,還采用了改善WENO 格式收斂能力的技術[32-33]來消除激波的波后振蕩.空間黏性項離散采用6 階精度的有限差分中心格式,時間導數項離散采用3 階精度的3 步TVD Runge-Kutta 格式[30].

5階精度的有限差分WENO-JS-5 格式形式如下:

5點總模板上的數值通量可以表示為3 個3 點子模板上3 階數值通量的凸組合,即

子模板上的數值通量為

非線性權為

其中,ε=10-6和p=2 .線性權為

光滑因子為

控制方程的對流項和黏性項經過空間離散后,得到一組常微分方程

其中,L(Q) 表示空間離散算子.3 階精度的3 步TVD Runge-Kutta 格式形式如下

1.3 計算域和初邊值條件

計算域和邊界條件的設置如圖1 所示.在物理域 -2D≤x≤15D和0 ≤r≤6D內,計算網格均勻分布,Δx=Δr=0.02D.在物理域之外,網格開始拉伸以消除反射聲波.

圖1 計算域和邊界條件示意圖Fig.1 The schematic of the computational domain and boundary conditions

計算域的上邊界和左邊界采用Thompson 的特征遠場邊界條件[34-35];在下游邊界區域,采用無反射出流邊界條件[36].在射流軸線上,由控制方程中含“ 1/r”的項所導致的奇性軸問題采用極限法[37]轉換成偏導數進行處理.此外,對軸線上的偏導數進行有限差分離散時,計算域也需拓展到的r＜0 的區域,相應的變量通過與r＞0 的區域的鏡像關系求得.噴管壁面采用無滑移邊界條件.在噴管出口處,入流平面內移6 個網格以便激發出不穩定波,減小對反饋環的數值限制和影響.噴嘴出口處的入流平面為聲速平面,無量綱形式的流動變量設置如下[15]

其中,比熱比 γ=1.4 ;Tr為儲氣罐溫度,根據冷射流假設,Tr=1 .

需要強調的是,在整個計算過程中,入流條件處并沒有添加任何形式的人工擾動.因此,嘯聲反饋環的產生和維持都不需要任何外部激勵.

初始時刻,除了噴管出口處,整個計算域內都設置為環境流動條件

上述參數表示射流由噴管射入周圍靜止的空氣.

為了方便對數值結果進行分析,計算選取固定的無量綱時間步長,具體取值為 Δt=0.002.為了確保時序信號傅里葉分析的精度,計算需要推進大量的時間步,總的無量綱積分時長為T=2500.本文所采用的數值方法、計算程序及針對軸對稱超聲速射流計算的網格無關性驗證見文獻[38].

2 數值模擬結果與分析

2.1 激波柵格結構

嘯聲聲波是由射流勢核中的激波柵格結構與剪切層中的大尺度擬序結構相互作用產生,因此,準確地模擬激波柵格結構的位置對于定位嘯聲的有效聲源位置至關重要.

圖2 是完全膨脹射流馬赫數為1.10 和1.15 的欠膨脹超聲速射流軸線上的時間平均壓力分布.為了驗證當前軸對稱數值模擬的準確性,將其結果與André等[39-40]的實驗結果進行了對比.圖中顯示,對于前3 個激波柵格,兩種狀態的數值解都與實驗測量值吻合得很好,無論是激波位置還是幅值.對于第4 個激波柵格,數值解的激波幅值與實驗測量值有所偏差,但激波位置吻合.在第5 個激波柵格結構之后的區域,數值解與實驗測量值差別較大,其原因可能是在此區域軸對稱控制方程物理上不再成立.由前文可知,嘯聲一般是在第5 個激波柵格的上游區域產生,所以這些偏差并不影響后續關于軸對稱嘯聲模態產生位置的分析.

圖2 射流軸線上的時間平均壓力分布及其與實驗結果[39-40]的比較Fig.2 The time-averaged pressure along the jet axis and its comparison with the experimental results[39-40] of screeching jet

圖2 射流軸線上的時間平均壓力分布及其與實驗結果[39-40]的比較 (續)Fig.2 The time-averaged pressure along the jet axis and its comparison with the experimental results[39-40] of screeching jet (continued)

2.2 嘯聲的頻率與聲壓級

在噴管唇口壁面上設置壓力監測點[0.0D,0.642D],記錄該位置處的時序壓力信號,采用快速傅里葉變換(單邊模式) 進行處理,得到近場噪聲頻譜信息.圖3是Mj=1.10,1.15 這兩種狀態欠膨脹超聲速射流的近場噪聲頻譜分析結果.圖中的橫坐標表示頻率,單位是 Hz,頻帶寬度是75 Hz;縱坐標表示聲壓級,單位是dB.對于Mj=1.10 的欠膨脹超聲速射流,有兩個尖銳的聲壓級凸起,分別是軸對稱A1模態及其諧頻.對于Mj=1.15 的欠膨脹超聲速射流,可以識別出3 個有意義的聲壓級凸起,分別是軸對稱A2模態及其第1、第2 諧頻.

圖3 射流近場噪聲頻譜信息:噴管唇口壁面監測點[0.0D,0.642D]Fig.3 Spectral information of jet near-field noise:monitor [0.0D,0.642D] located on the nozzle exit lip wall

表2 具體給出了上述離散頻率模態的頻率和聲壓級,為了進行對比,同時列出了Ponton等[13-14]的實驗測量值和Loh等[41]的數值模擬結果.對于Mj=1.10 的射流,當前數值結果與Loh等[41]的數值結果都捕捉到了軸對稱A1模態及其諧頻,且頻率和聲壓級吻合得很好.對于Mj=1.15 的射流,不論是模態類型,還是模態頻率和聲壓級,當前數值結果都與實驗結果[13-14]要更為吻合.

表2 嘯聲模態及其諧頻的頻率和聲壓級Table 2 The frequencies and amplitudes of screech modes and their harmonics

2.3 軸對稱A1 模態嘯聲的聲源位置

圖4 是在完全膨脹射流馬赫數Mj=1.10 的欠膨脹射流中與A1模態相關的瞬時流場結構和聲場.從圖中可以看到3 道向上游傳播的嘯聲聲波.采用Edgington-Mitchell等[21]的方法,在圖4 中疊加一個與嘯聲聲波重合的圓(藍色虛線),通過圓心來定位聲源位置.圖4 顯示,圓心位于流向位置x=2.29D處,該位置恰好是第4 個激波柵格結構的尾緣,也稱激波反射點或激波尖梢;根據紋影和渦量還可以確定,在此位置處,剪切層內同時出現了一個旋渦鞍點.激波泄漏機制[17-18]也指出激波會在剪切層鞍點處向外泄漏,形成嘯聲聲波.因此,可以確定第4 個激波柵格結構的尾緣是軸對稱嘯聲A1模態的有效聲源位置.

圖4 Mj=1.10 的欠膨脹超聲速射流中與嘯聲A1 模態相關的瞬時流動結構和聲場Fig.4 The instantaneous flow structures and acoustic field associated with screech A1 mode in the underexpanded supersonic jet with Mj=1.10

2.4 軸對稱A2 模態嘯聲的聲源位置

通過快速傅里葉變換,對射流壓力場進行傅里葉模態分解,即將時域空間的脈動壓力場時間序列變換到頻域空間,提取嘯聲頻率成分,逐點計算聲壓級和相位,得到嘯聲頻率上聲壓級和相位的空間分布.聲壓級最大的位置可認為是嘯聲模態的產生位置.本節的壓力場傅里葉模態分解共選取了8192 個時刻的數據,起始時刻為t=200,無量綱時間間隔為δt=0.1.

圖5是Mj=1.15 的欠膨脹超聲速射流中嘯聲A2模態頻率上的聲壓級空間分布和勢核區的時間平均壓力空間分布.在圖5(a)聲壓級云圖中可以看到沿射流剪切層外緣分布的嘯聲A2模態關聯駐波結構.駐波結構是向下游傳播的剪切層不穩定波和向上游傳播的聲波相互作用所產生[24].圖5(b)中的時間平均壓力分布展現的是射流勢核內的激波柵格結構.結合圖5(a)和圖5(b)可知,聲壓級最大的位置是在流向位置x=2.14D處,即第3 個激波柵格結構尾緣.因此,可以確定嘯聲A2模態的聲源位于第3 個激波柵格結構的尾緣.

圖5 Mj=1.15 的欠膨脹超聲速射流Fig.5 The underexpanded supersonic jet with Mj=1.15

圖6是Mj=1.15 的欠膨脹超聲速射流中嘯聲A2模態頻率上的相位空間分布,從圖中可以清晰地看到聲波傳播的波振面.為了進一步確定A2模態的聲源位置,圖7 提取了直線R/D=1.5 上的相位分布.相位斜率的符號表示聲波的傳播方向,而相位斜率符號變化的位置也可指示聲源位置[26,28-29].圖中顯示,相位斜率的符號在流向位置x=2.14D附近發生了明顯變化,在其上游,斜率為正值,聲波向上游傳播,而在其下游,斜率為負值,聲波則向下游傳播;該位置同時也是圖5(a)中聲壓級最大的位置,表明嘯聲A2模態的聲源位置的確是在第3 個激波柵格結構的尾緣.

圖6 Mj=115 的欠膨脹超聲速射流中嘯聲A2 模態頻率上的相位空間分布Fig.6 The spatial distribution of phase at the frequency of screech A2 mode in underexpanded supersonic jet with Mj=1.15

圖7 Mj=1.15 的欠膨脹超聲速射流中嘯聲A2 模態頻率上的相位沿直線R/D=1.5 的分布Fig.7 The distribution of phase at the frequency of screech A2 mode along the line of R/D=1.5 in underexpanded supersonic jet with Mj=1.15

2.5 基于本征模態分解的聲源位置分析

本征正交分解(proper orthogonal decomposition,POD)是提取、分析湍流中擬序流動結構的重要工具[42].POD 方法的實質是在最小二乘意義下提供能夠代表已知數據的一組正交基.POD 快照法[43]的目標是尋找脈動速度場數據快照u′的最小二乘表達式,即殘差ures的模最小;相應的數學表達式[23,44]如下

其中,bi是時間POD 模態,ψi是空間POD 模態.M是模態數,M≤N,N是數據快照序列的快照總數.

為了得到POD 模態,首先構造脈動速度場的快照矩陣V,隨后計算脈動速度場的自相關(對稱正定)矩陣R=VTV.

針對自相關矩陣進行特征分解Rai=λiai,求得自相關矩陣的特征值 λi和特征向量ai

空間POD 模態可以表示為脈動速度場快照與時間POD 模態乘積的線性組合

式中特征值 λi表征了每個空間POD 模態兩倍的不可壓脈動動能.因此,POD 方法能夠提取最富能的流動結構.特別地,為了表征每個POD 模態的能量占比,需要對特征值 λi進行歸一化處理,并按大小進行排序.需要注意的是,POD 模態通常是成對出現的.

本節的POD 分析選取了共2048 個時刻的速度場數據快照,無量綱時間間隔 δt=0.1 .

圖8是Mj=1.10 的欠膨脹超聲速射流速度場的時間POD 模態,包括波形、相圖和特征值.POD 模態都是成對出現的,圖8(a)和圖8(b)是1 階和2 階POD 時間模態[a1(tk) 和a2(tk),k是快照時刻]及其相圖.相圖能夠表征各階POD 模態間的時間相關性,表示兩個模態的平均幅值.1 階和2 階模態的波長是相匹配的;相圖中的點大致呈圓形散布,表明前兩階模態是時間強相關的.時間相關性以及匹配的波長說明1 階模態和2 階模態是耦合的,兩者共同表征了1 個振蕩過程或擬序結構.圖8(c)是各階POD 模態的能量占比,1 階和2 階模態含能最高,兩者之和接近45%,大于其它各階模態.

圖8(d)是1 階和3 階POD 模態的相圖.圖中的點呈雙扭線形分布,表明1 階、3 階模態之間沒有類似1 階、2 階模態的時間相關性.圖8(e)和圖8(f)是5 階和6 階POD 時間模態[a5(tk)和a6(tk) ] 及其相圖,同樣表示兩者的平均幅值.與前兩階模態類似,5 階和6 階模態的波長相匹配,并且是時間強相關的,兩者也共同表征了1 種擬序結構.

圖8 脈動速度場的POD 時間模態和特征值,Mj=1.10Fig.8 Temporal POD modes and their eigenvalues of fluctuating velocity field,Mj=1.10

圖9是Mj=1.10 的欠膨脹超聲速射流速度場第1 對和第3 對POD 時間模態的傅里葉分析.一般來說,每個POD 模態包含多種頻率組分.圖9 顯示,這兩對POD 模態都有一個很強的主導頻率,其中,第1 對的主導頻率是4396 Hz,第3 對的主導頻率是8595 Hz.第3 對POD 模態的主導頻率與嘯聲A1模態的頻率(8560 Hz)十分一致,即第3 對POD 模態表征的是與嘯聲A1模態相關的擬序流動結構.第1 對POD 模態的主導頻率接近第3 對POD 模態的1/2,因此,其所表征的擬序流動結構是嘯聲A1模態關聯擬序流動結構的亞諧頻.

圖9 脈動速度場的POD 時間模態的傅里葉分析,Mj=1.10Fig.9 The FFT analysis of the first and the third pairs of temporal POD modes of fluctuating velocity field,Mj=1.10

圖10是Mj=1.10 的欠膨脹超聲速射流速度場的第1 對和第3 對POD 空間模態,包括u分量和v分量.圖中顯示,空間POD 模態的u分量存在著非常明顯的交界面,在交界面兩側符號相反,這一交界面是剪切層與射流主流(勢核)的交界面;此外,每1 對空間POD 模態之間都存在著1/4 波長的相位差.第1 對POD 模態所表征的擬序結構在流向位置x＞3.8D的區域達到飽和態;第3 對POD 模態所表征的嘯聲相關擬序結構發展的更快,在流向位置1.5D＜x＜3.8D的區域就達到飽和態,隨后開始衰減.前文2.3 節所述嘯聲A1模態的聲源位置(x=2.29D)就位于該飽和態區域內.在第3 對POD 模態中,還可以清楚地看到從飽和態區域發出的,規律性地向上游和下游輻射的條紋,而向上游傳播的條紋即為A1模態的嘯聲.根據主導頻率和空間演化,第1 對POD 模態所表征的擬序結構是由嘯聲A1模態關聯擬序結構在下游發生旋渦合并所產生.

圖10 脈動速度場的POD 空間模態,Mj=1.10Fig.10 The first and the third pairs of spatial POD modes of fluctuating velocity field,Mj=1.10

通過對不同流向位置處的脈動速度場切片進行POD 分析,可以研究擬序結構能量占比沿流向的發展變化.在每一流向站位,提取特定時間間隔內的

切片速度場時間序列組成局部數據快照矩陣,然后對局部數據快照矩陣進行POD 分析.

圖11是Mj=1.10 的欠膨脹超聲速射流速度場POD 模態特征值沿流向的變化曲線.圖中顯示,在流向位置x=3.1D處,第2 對POD 模態特征值之和的極大值點恰好對應第1 對POD 模態特征值之和在此處的極小值點;第1 對POD 模態所表征的擬序結構在此位置前經歷了較為強烈的能量釋放,而第2 對POD 模態所表征的擬序結構在此位置前一直在吸收能量,表明這兩種擬序結構之間存在著能量傳遞.第3 對POD 模態特征值之和沿流向的變化表明,嘯聲A1模態關聯擬序結構在其飽和態區域(1.5D＜x＜3.8D)內,先后經歷了能量吸收再釋放的過程,并且在能量釋放子區間內,其能量占比,即第3 對POD 模態的特征值之和,在流向位置x=2.29D處產生了一個局部極小值點,這個局部極小值點就是前文2.3 節中確定的嘯聲A1模態的產生位置.

圖11 脈動速度場軸向切片POD 分析中相關模態的特征值沿流向的變化,Mj=1.10Fig.11 The variations of eigenvalues along streamwise of the POD modes based on the POD analysis of velocity field’s slices in different axial position,Mj=1.10

2.6 基于動態模態分解的聲源位置分析

動態模態分解[45](dynamic mode decomposition,DMD)方法則是基于動力系統Koopman 分析發展而來的分析工具,其通過時間正交的特征模態來重構數據快照序列.POD 模態強制空間正交性(結構去相關),每個獨立的POD 模態都包含了多種頻率成分,而DMD 模態是時間正交的(單一頻率),一般空間不正交.如果數據序列足夠長,DMD 方法能夠識別主頻,提取相關的流動特征.對源自于飽和非線性過程的實驗或數值數據來說,相較于提取按能量排序的空間去相關結構,提取帶有特定頻率信息的時間去相關結構更為重要.

基于DMD 模態的數據快照重構表達式[46]為

式中,?k是DMD 模態,K是DMD 模態的數目,ωk=lnμk/Δt的虛部是DMD 模態的圓頻率,而實部則表示指數增長或衰減,具體取決于符號;bk(0) 是每個DMD 模態的初始幅值.初始時刻的數據快照x1=Φb,所以DMD 模態的初始幅值可通過b=Φ+x1求得.

本節針對壓力場的DMD 分析選取與2.5 節POD 分析中相同時刻的數據快照序列.具體求解過程如下.

其中,μj和yj分別是特征值和特征向量.

其中,U和W是酉矩陣,Σ 是奇異值對角矩陣.

最終可得DMD 模態

需要說明的是,本節后續分析中的DMD 模態都是包括了初始幅值的DMD 模態,即bj?j.

圖12是Mj=1.10 的欠膨脹超聲速射流的壓力場DMD 模態.各階DMD 模態的頻率分布見圖12(a),其中,幅值最高的前4 階DMD 模態的頻率依次是4383 Hz,8570 Hz,4252 Hz 和4186 Hz.對比2.5 節中Mj=1.10 的欠膨脹射流速度場POD 分析結果可知,壓力場的1 階、3 階和4 階DMD 模態的頻率都十分接近速度場第1 對POD 模態的主導頻率(4396 Hz),并且三者具有相似的空間形態,這里只給出了3 階DMD 模態,見圖12(b),其所表征的擬序結構也是在流向位置x＞3.8D達到飽和態.2 階DMD模態的頻率接近第3 對速度場POD 模態的主導頻率(8595 Hz),此頻率也是嘯聲A1模態的頻率(8560 Hz),因此,2 階DMD 模態表征的是嘯聲A1模態的關聯擬序結構,其飽和態區域為 1 .5D＜x＜3.8D,見圖12(c),從中可以觀察到從飽和態區域向外輻射的嘯聲A1模態聲波條紋.根據頻率,3 階DMD 模態表征的擬序結構是2 階DMD 模態所表征擬序結構的亞諧頻;而且在空間演化上,兩者的飽和態區域是接續的,因此,前者是由后者在下游發生旋渦合并后所產生.

圖12 Mj=1.10 的欠膨脹超聲速射流的脈動壓力場DMD 模態Fig.12 The DMD modes of fluctuating pressure field in the underexpanded supersonic jet with Mj=1.10

圖12(d)是Mj=1.10 的欠膨脹超聲速射流壓力場2 階DMD 模態的相位空間分布,清晰地表征了A1模態嘯聲聲波傳播的波振面.圖13 提取了2 階DMD 模態在直線R/D=2.1 上的相位分布.前文2.4 節中已經說明相位斜率符號改變的位置可以指示嘯聲的聲源位置.2 階DMD 模態相位斜率的符號是在流向位置x=2.29D處發生改變,該位置也是前文2.3 節和2.5 節中所確定的嘯聲A1模態的產生位置,這再次證實嘯聲A1模態的有效聲源位置是在第四個激波柵格結構的尾緣.

圖13 2 階DMD 模態的相位在直線(R/D=2.1)上的分布,Mj=1.10Fig.13 The distribution of phase along the line R/D=2.1 for the second DMD mode,Mj=1.10

圖14是Mj=1.15 的欠膨脹超聲速射流的壓力場DMD 模態.圖14(a)是各DMD 模態的頻率分布,其中,幅值最高的前4 階DMD 模態的頻率依次是4579 Hz,4644 Hz,4710 Hz 和9289 Hz.前3 階DMD 模態的頻率近似相等,而且也具有相似的空間形態,因此,這里只給出了2 階模態,見圖14(b).4 階DMD 模態接近嘯聲A2模態的頻率(9250 Hz),因此,它表征的是嘯聲A2模態關聯擬序結構,圖14(c)顯示其飽和態區域為 1 .65 ＜x＜3.25,并且可以觀察到自此區域向外輻射的嘯聲A2模態聲波條紋.2 階DMD 模態是4 階DMD 模態的亞諧頻,其所表征的擬序結構在流向位置x＞3.25 處達到飽和態.同樣地,根據頻率和空間演化可知,前者也是由后者在下游發生旋渦合并所產生.

圖14 Mj=1.15 的欠膨脹超聲速射流的脈動壓力場DMD 模態Fig.14 The DMD modes of fluctuating pressure field in the underexpanded supersonic jet with Mj=1.15

圖14(d)是Mj=1.15 的欠膨脹超聲速射流壓力場4 階DMD 模態相位的空間分布.與圖12(d)類似,相位空間分布清晰地表征了A2模態嘯聲聲波傳播的波振面.圖15 提取了4 階DMD 模態在直線R/D=1.5 上的相位分布.沿流向,2 階DMD 模態有多個相位斜率符號發生改變的位置,但只有一個位于飽和態區域之內,即x=2.14D.圖14(c)也顯示嘯聲聲波是在飽和態區域內向外輻射,因此,飽和態區域外的相位斜率符號轉變位置與嘯聲的產生沒有直接關聯.此外,x=2.14D這一位置也與2.4 節中傅里葉模態分解得到的結果一致,進一步證實嘯聲A2模態的產生位置是在第3 個激波柵格結構的尾緣.

圖15 4 階DMD 模態的相位在直線(R/D=1.5)上的分布,Mj=1.15Fig.15 The distribution of phase along the line (R/D=1.5) for the fourth DMD mode,Mj=1.15

3 結論

為了研究軸對稱模態嘯聲的產生位置,本文采用高階精度有限差分方法直接求解軸對稱可壓縮Navier-Stokes 方程,數值模擬了完全膨脹射流馬赫數分別為1.10 和1.15 的圓形聲速噴管欠膨脹超聲速冷射流,得到了軸對稱嘯聲A1模態和A2模態.通過激波柵格結構位置和嘯聲模態頻譜與文獻實驗結果的對比,驗證了數值模擬的準確性.

基于傅里葉模態分解、本征模態分解和動態模態分解,分析了射流時序壓力場和速度場,研究了嘯聲關聯擬序流動結構的空間演化,精確定位了軸對稱模態嘯聲的有效聲源位置.

研究表明:嘯聲關聯擬序流動結構存在飽和態區域,并且嘯聲聲波是在其飽和態區域內產生并向外輻射.在本文所涉及的射流馬赫數范圍內,A1和A2兩種軸對稱模態嘯聲的有效聲源位置分別是在第4 和第3 個激波柵格結構的尾緣.因此,本文的數值模擬結果支持A1和A2兩種嘯聲模態是單聲源的觀點.

需要注意的是,本文只研究了低馬赫數超聲速射流中的軸對稱模態嘯聲,對于出現在高馬赫數射流中的擺動模態嘯聲和螺旋模態嘯聲沒有涉及,它們具體的聲源位置以及是否為單聲源,仍有待進一步研究.