凹腔對高超聲速邊界層穩(wěn)定性的影響

馮 爍，王伯福，涂國華，武健輝，陳堅強，楊 強，*

(1.上海大學 力學與工程科學學院，上海市應用數(shù)學和力學研究所，上海 200072；2.空天飛行空氣動力科學與技術全國重點實驗室，綿陽 621000；3.中國運載火箭技術研究院，北京 100076)

0 引言

層流到湍流的轉捩可導致飛行器表面摩阻熱流劇烈變化，在高超聲速時氣動力/熱上升3～5 倍，嚴重影響飛行安全，是制約高超聲速飛行器氣動設計的重要因素[1]。邊界層轉捩作為一個復雜的初邊值問題，對壁面粗糙度極為敏感。壁面凹陷和凸起是實際飛行器表面常見的粗糙元形式。典型粗糙元依據(jù)其形狀可以分成凹陷型和凸起型，其中凹陷型又稱為凹腔。對于矩形凹腔，寬深比d/h是影響凹腔內流動結構的重要因素之一[2-4]，依據(jù)不同寬深比，凹腔對應的流動狀態(tài)可分為三類：開式凹腔、過渡式凹腔和閉式凹腔。對于閉式凹腔，隨著凹腔深度的增加，邊界層主層對平均流的修正逐漸增強。當凹腔過渡為開式凹腔時，凹腔的核心區(qū)產生由主流邊界層和腔內回流區(qū)組成的混合層，這種混合層會在凹腔局部產生增長率更快的Kelvin-Helmholtz 不穩(wěn)定性。Guo 等[5]指出來流條件也是影響腔內流動結構的重要因素之一。Palmer 等[6]發(fā)現(xiàn)高超聲速邊界層在凹腔周圍和下游的流動形態(tài)與凹腔的幾何形狀和來流條件密切相關。本文研究了多個寬度、深度的凹腔模型，針對這種密集的尺寸，舊有的分類方法界限模糊，不利于得到定量的數(shù)據(jù)結論，因此本文研究討論還是直接使用了尺寸的大小分類。董昊等[7]從幾何形狀和分布形式兩方面對粗糙元的分類進行了較為完整的總結。

粗糙元尺寸對于邊界層穩(wěn)定性起著至關重要的作用。對于孤立表面凸起粗糙元，目前對其控制轉捩機理及效果的研究主要集中于粗糙元的幾何參數(shù)，其中粗糙元高度對轉捩效果起著至關重要的作用[8]，當高度改變時，粗糙元邊界層的感受性和失穩(wěn)機制也隨之改變。Sumariva 等[9]研究了不同長度、高度和形狀的凸起孤立粗糙元引起的轉捩提前現(xiàn)象。Brehm等[10]針對粗糙元形狀、高度和間距等因素對擾動的演化影響進行了研究，發(fā)現(xiàn)存在一個粗糙元高度閾值，低于該閾值的粗糙元對擾動增長的影響很小，與Gaster[11]的結論一致；同時發(fā)現(xiàn)形狀和間距對二維擾動波的增長具有顯著影響。

借助局部散射理論對擾動波控制效果開展定量分析，Wu 等[12]針對較小尺度粗糙元開展研究，發(fā)現(xiàn)當粗糙元尺度與不穩(wěn)定波波長相當時，粗糙元對T-S波起促進作用。李闖和董明[13]采用直接數(shù)值模擬研究了局部矩形凹腔對高超聲速平板邊界層中第二模態(tài)擾動演化的影響，發(fā)現(xiàn)對于較淺的凹腔，低頻第二模態(tài)擾動被促進，而高頻擾動則被抑制。通常，隨著凹腔深度的增加，凹腔對擾動的抑制作用增強，但當凹腔的深度超過某一臨界值時，凹腔內流動出現(xiàn)新的失穩(wěn)機制，且與凹腔內基本流的拓撲結構密切相關。對于淺凹腔，由于前后緣可以獨立地看作后向和前向臺階，前后緣分離泡相距較遠，相互作用較弱；而對于深凹腔，分離泡會占據(jù)整個凹腔，進而引起更多動能向回流區(qū)域傳遞[14]。

粗糙元位置對邊界層穩(wěn)定性也有重要影響。在高超聲速邊界層中，源于自由流的快、慢聲波在向下游演化過程中激發(fā)了邊界層中的快、慢模態(tài)（Mode F、Mode S），二者相速度相同時會出現(xiàn)同步[15]，進而激發(fā)Mack 模態(tài)[16-17]?？?、慢模態(tài)發(fā)生同步的位置稱為同步點。Duan 等[18]研究發(fā)現(xiàn)凸起粗糙元位于單頻擾動快、慢模態(tài)的同步點上游時對擾動起促進作用，而位于同步點下游時則對擾動起抑制作用。Chen、Lee[19]和Zhao 等[20]分別通過實驗和數(shù)值模擬確認了凹腔和快、慢模態(tài)同步點位置的關系，發(fā)現(xiàn)影響效果與凸起粗糙元類似。Park 等[21-22]采用拋物化穩(wěn)定性方程研究了不同形狀二維光滑凸起粗糙元對高超聲速平板邊界層線性不穩(wěn)定性的影響，也獲得了類似的結論。

綜上可見，針對凹腔對邊界層穩(wěn)定性影響的局部參數(shù)化研究已有較豐富的成果，還需進一步研究凹腔尺寸變化對Mack 模態(tài)的全面影響規(guī)律。為此，本文將借助線性穩(wěn)定性理論（linear stability theory,LST）、拋物化穩(wěn)定性方程（parabolic stability equation,PSE）、高精度直接數(shù)值模擬（direct numerical simulation,DNS）和全局穩(wěn)定性分解（global stability decomposition,GSD）系統(tǒng)分析凹腔對高超聲速和超聲速邊界層穩(wěn)定性的影響，以期給飛行器設計和轉捩控制提供參考。

1 物理模型及計算方法

1.1 物理模型

Fong等[23-24]發(fā)現(xiàn)了凸起粗糙元分別位于特定頻率（即同步頻率）的同步點前后時的擾動波幅值增長變化規(guī)律。Zhao、Dong、李斯特等[20,25-26]借助局部散射理論系統(tǒng)研究了粗糙元參數(shù)對擾動波演化的影響規(guī)律。Han 等[27]針對超聲速凹腔做了自激振蕩的模態(tài)分解。為了更好地與這些相關文獻對照補充，本文研究模型為帶凹腔的零攻角平板，來流馬赫數(shù)為Ma=5.92、1.38，與文獻[20,23-27]一致，計算工況詳見表1。本文重點主要在：1）超聲速和高超聲速來流條件下的對比；2）更廣泛的參數(shù)范圍；3）大尺度凹腔的全局不穩(wěn)定性研究。

表1 算例計算工況Table 1 Flow conditions for the simulation cases

為了方便表述，流體變量均基于特征參考量進行無量綱化[28]，對應的有量綱量用上標(●)*表示，參考長度選取為Lref=1.0 mm。計算域和邊界條件如圖1所示，坐標原點位于平板前緣。

圖1 計算域及邊界示意圖Fig.1 Schematic diagram of the computational domain and boundaries

為減少計算量，計算域流向區(qū)間為[x0,xl]=[100,400]，法向高度為Ly=100。左邊界為入口條件，給定可壓縮Blasius 相似解。右邊界為出口邊界，采用線性外推，計算擾動演化時，在出口附近采用網(wǎng)格拉伸處理，以消除虛假反射波的影響。上邊界采用遠場條件。下邊界為無滑移無穿透的固體壁面，對于高超聲速工況，給定壁溫Tw=334.9 K，而對于超聲速工況，給定絕熱壁面邊界，(?T/?y)|wall=0。已有研究發(fā)現(xiàn)，相對于絕熱壁面，等溫熱壁促進第一模態(tài)增長，而抑制第二模態(tài)增長，等溫冷壁則起相反的作用。在線性增長階段，不同壁溫條件對增長率的影響并不大[29-31]。

本文選取了不同尺寸的凹腔模型開展研究，參數(shù)如表2 所示。寬度d/δ99取值0.5、2.0、3.0、4.0、6.0，深度h/δ99取值0.5、1.0、1.5、2.0、2.5、3.0、4.0，其中δ99為凹腔中心位置處的當?shù)剡吔鐚雍穸?。設計算例可分為固定寬度變深度算例組、固定深度變寬度算例組。參照文獻[20]，凹腔中心點位置固定在x=185，對于Ma=5.92 和Ma=1.38 工況，對應的當?shù)剡吔鐚雍穸确謩e為δ99*=2.36 mm 和δ99*=2.45 mm。為方便描述，當凹腔尺寸小于當?shù)剡吔鐚雍穸然蚺c當?shù)剡吔鐚雍穸认喈敃r，稱為小尺度凹腔；當凹腔寬度和深度都已明顯大于當?shù)剡吔鐚雍穸葧r，稱為大尺度凹腔，在本文中，當凹腔尺寸滿足d/δ99＞6、h/δ99＞3 時，定義為大尺度凹腔。

表2 不同類型凹腔尺寸參數(shù)Table 2 Geometry parameters for different types of cavities

1.2 計算方法

本文采用基于LST 和PSE 的線性穩(wěn)定性分析、DNS 和GSD 三類方法開展凹腔流動分析。流動的控制方程為可壓縮Navier-Stokes（N-S）方程，計算采用理想氣體狀態(tài)模型，黏性系數(shù)、熱傳導系數(shù)由Sutherland 公式計算給出。對流項采用Steger-Warming通量分裂，在激波附近采用五階加權緊致非線性格式（weighted compact nonlinear schemes,WCNS）[32-33]，激波以外區(qū)域采用六階迎風格式，黏性項采用六階中心差分，時間項離散采用具有TVD 性質的3 步三階Runge-Kutta 法。

采用當?shù)鼗木€性穩(wěn)定性分析方法LST 進行穩(wěn)定性分析，模態(tài)解表示形式如式（1）所示：

對于大尺度凹腔，流動在流向也具有較強的非平行性，為此采用Tu 等[36]提出的基于無矩陣的GSD方法開展流動穩(wěn)定性分析。該方法的思想是通過反復調用CFD 求解器構造原特征矩陣A（m×m維）的子空間S（n?m）實現(xiàn)降維處理，進而找到工程上最關心的數(shù)個主特征值。具體思路如下：

選擇一系列線性無關的初始擾動：

通過高精度CFD 工具，讓擾動場演化 ?t時間后得到：

其中矩陣B=exp(A?t)?？梢酝队疤幚順嬙煨【仃嚕?/p>

其中，Q、R對 應Um的QR分解矩陣，滿足Um=QR，且可以證明矩陣S與矩陣A的特征值和特征向量滿足：特別地，當采用Arnoldi 方法構造子空間時，GSD 可退化為傳統(tǒng)的動態(tài)模分解DMD。

1.3 網(wǎng)格收斂性

圖2 展示了d2h1H 凹腔附近的網(wǎng)格分布情況，法向每隔4 個網(wǎng)格點顯示，流向每隔2 個網(wǎng)格點顯示。網(wǎng)格沿法向進行拉伸，邊界層內網(wǎng)格數(shù)保持在150 以上，以保證近壁區(qū)足夠的網(wǎng)格分辨率，沿流向對凹腔的前后緣進行加密以準確捕捉此處的波系結構。

圖2 凹腔局部區(qū)域網(wǎng)格分布Fig.2 Grid distribution around the cavity

圖3 為d6h1.5H 算例加密網(wǎng)格（流向和法向網(wǎng)格數(shù)為腔內501×201，腔外 2 001 × 1 001 ）與原始網(wǎng)格（腔內251×101，腔外 1 001 ×501）計算結果對比圖。隨著凹腔寬度和深度的逐漸增大，腔內網(wǎng)格點數(shù)也逐漸增多。為驗證計算結果的網(wǎng)格無關性，選取d6h1.5H 算例進行研究，可看到壁面壓力系數(shù)和壓力脈動幅值隨流向的演化曲線在流向和法向網(wǎng)格加密前后完全吻合，說明了本文計算結果具有網(wǎng)格無關性。

圖3 網(wǎng)格收斂性驗證Fig.3 Grid convergence verification

2 基本流特征

入口（x0=100）采用Blasius 剖面以獲得計算域內的邊界層流場。圖4 展示了Ma=5.92 下d2h1H、d2h3H、d6h1H 和d6h3H 算例的基本流特征。

圖4 不同凹腔壓力云圖及流線圖Fig.4 Pressure contours and streamlines around cavities

從壓力云圖可以看出不同尺度的凹腔均呈現(xiàn)類似的波系結構，在凹腔前緣出現(xiàn)膨脹波系，而在后緣出現(xiàn)壓縮波系，且二者均隨凹腔尺度的增加而增強。從流線圖可以清晰地看出凹腔內的渦結構分布：對于d2h1H 算例（圖4（a）），凹腔內存在單個穩(wěn)定的順時針旋轉主渦結構，隨著深度的增加（d2h3H 算例，圖4（b）），凹腔下底面左右側開始形成小的角渦結構；隨著寬度的增加（d6h1H 算例，圖4（c）），主渦結構被拉長，后緣渦強度高于前緣渦強度，但前后緣不對稱性減弱。對于較大尺度凹腔d6h3H，上述渦結構特征變化更為明顯（圖4（d）），凹腔內呈現(xiàn)復雜的多渦系結構。

圖5 展示了Ma=1.38 時d6h3S 算例的非定常流動狀態(tài)。定常情況下Ma=1.38 的波系結構與Ma=5.92的高超聲速情況類似，當固定寬度增加深度或者固定深度增加寬度時，凹腔前緣膨脹波均增強，后緣開始出現(xiàn)明顯的壓縮波，并且凹腔內部的角渦范圍明顯增大，內部回流區(qū)流動變得復雜。從圖5 可見凹腔內發(fā)生自激震蕩并在下游形成規(guī)律性周期振蕩，表明大尺度凹腔結構自身引起的絕對不穩(wěn)定性變得重要。

圖5 d6h3S 算例凹腔不同時刻的|?ρ|Fig.5 Density gradient |?ρ| contours around d6h3S at different time instances

圖6 和圖7 展示了凹腔寬度固定時，變深度組凹腔不同站位的速度型和壁面處速度剪切沿流向的變化。從速度型剖面可以看出，凹腔前站位速度型與Blasius 相似解吻合較好，離開凹腔很快又恢復到相似解（圖6）。Ma=1.38 超聲速情況下（圖6（b）），凹腔附近速度型與高超聲速情況類似，但凹腔前緣膨脹波較弱，前緣壁面速度型吻合度更高。相比于凸起粗糙元，凹腔對下游邊界層的基本流修正要弱得多[20]。由于凹腔內部分離泡結構的存在，速度型在凹腔內存在較大差異。隨著凹腔深度的增加，流向速度在凹腔頂部剪切層附近（y=0）單調增加，但是由于腔內回流區(qū)的存在，流向速度先減小到負值再回到凹腔底部壁面處的零值。

圖6 變深度組凹腔流向速度剖面對比Fig.6 Streamwise velocity profiles for cavities with different depths

圖7 變深度組凹腔壁面處速度剪切特征對比Fig.7 Wall shear rates for cavities with different depths

對于壁面處速度剪切（見圖7），在接近凹腔前后緣時剪切突增，表明前后緣附近流速突然增加，后緣附近明顯大于前緣。小深度的回流區(qū)中心靠近右側壁面，并且隨著深度增加逐漸向前緣移動，但是當深度增加到2 倍邊界層厚度（圖7（a））時，再增加凹腔深度不再有明顯的剪切極值。相比于高超聲速情形，超聲速來流下（圖7（b））凹腔內部壁面處剪切剖面僅存在單峰結構且峰值更加明顯，說明回流區(qū)更加集中于凹腔后緣。

圖8（a、b）展示了Ma=5.92 下固定寬度組和固定深度組凹腔內壁面壓力系數(shù)分布。如圖8（a）所示，在凹腔前部，隨著深度從0.5δ99（d2h0.5H 算例）開始逐漸增大到4δ99（d2h4H 算例），負壓區(qū)大小變化很?。欢诎记粌炔?，淺凹腔（d2h0.5H 算例）產生了較強的逆壓梯度，且隨著深度增加逆壓梯度逐漸減弱。如圖8（b）所示，隨著寬度的增加（由d0.5h1H 算例到d6h1H 算例），前后緣局部壓力極值增大，凹腔前部區(qū)域擁有更廣的負壓區(qū)，這說明凹腔寬度增加導致了前后緣的壓力修正進一步增強。圖8（c、d）展示了Ma=1.38 的情況，值得注意的是，超聲速情況下腔內逆壓梯度明顯減弱。

圖8 變深度和變寬度凹腔壁面壓力系數(shù)Cp 對比Fig.8 Comparison of wall pressure coefficient Cp for cavities with varying depths and widths

經過以上DNS 分析，獲得了不同凹腔對基本流的影響規(guī)律，下文將進一步通過穩(wěn)定性分析手段研究凹腔對邊界層穩(wěn)定性特征的影響。

3 非定常擾動演化

3.1 入口擾動選取

已有研究表明，凹腔安裝位置相對同步點位置而言對邊界層穩(wěn)定性具有重要影響[14,19,22,37-38]。圖9給出了Ma=5.92、x=185 站位的快慢模態(tài)相速度隨頻率的變化，可以看出快慢模態(tài)在ω=0.98（對應的頻率時相速度相等，發(fā)生同步，演化曲線與文獻[20]的結果一致。

圖9 Ma=5.92、x=185 站位的快慢模態(tài)相速度隨頻率變化曲線Fig.9 Phase velocity variations with frequency for fast and slow modes at Ma=5.92 and x=185

圖10（a）給出了本文研究的Ma=5.92 平板邊界層的中性曲線（圖中藍色空心圓圈代表快慢模態(tài)同步點站位），可見第一模態(tài)的上支界與第二模態(tài)的下支界重合，平板前緣很快進入穩(wěn)定性理論分析得到的不穩(wěn)定區(qū)間。圖10（b）給出了Ma=1.38 工況的中性曲線，由于來流馬赫數(shù)較低，僅出現(xiàn)第一模態(tài)，未出現(xiàn)穩(wěn)定性理論預測的高階Mack 模態(tài)。

圖10 平板邊界層中性曲線Fig.10 Neutral curves for flat-plate boundary layers

針對Ma=5.92 時的凹腔尺寸，選取了7 個典型的不穩(wěn)定頻率，即ω=0.30、0.53、0.76、0.98、1.07、1.14 和1.45，對應的有量綱頻率、同步點流向站位及無量綱擾動波穩(wěn)定性特征參數(shù)見表3?？梢?，不穩(wěn)定波頻率越高，對應的同步點越靠近上游，其中最低頻的兩個擾動（40 kHz 和70 kHz）對應的同步點流向站位已超出計算域。

表3 Ma=5.92 時入口擾動波參數(shù)Table 3 Parameters for the inlet disturbance wave at Ma=5.92

對于Ma=1.38，由于不存在同步點，并且其不穩(wěn)定區(qū)間較窄，因此選取凹腔中心站位處最不穩(wěn)定頻率ω=0.095 附近的三個典型頻率開展分析，具體參數(shù)如表4 所示。計算域入口不穩(wěn)定波形狀函數(shù)由LST 求解得到。

表4 Ma=1.38 時入口擾動波參數(shù)Table 4 Parameters for the inlet disturbance wave at Ma=1.38

3.2 凹腔對邊界層線性穩(wěn)定性的影響

圖11 展示了Ma=5.92 來流條件下平板、基準凹腔（d2h1H 算例）、深凹腔（d2h3H 算例）和寬凹腔（d6h1H 算例）對應的擾動增長率云圖。結果表明，凹腔對于增長率的影響較為局部，且大尺度凹腔的局部影響更為顯著，但增長率在凹腔下游均重新恢復到指數(shù)增長形式。圖12展示了Ma=1.38 來流條件下的基準凹腔（d2h1S 算例）和寬凹腔（d6h1S 算例）對應的擾動增長率云圖，可見增長率在凹腔附近被局部促進，并且隨著凹腔寬度增大促進效果也增強。

圖11 Ma=5.92 來流下不同尺寸凹腔的邊界層擾動增長率云圖Fig.11 Growth rate contours for boundary layers with different sizes of cavities at Ma=5.92

圖12 Ma=1.38 來流下不同尺寸凹腔的邊界層擾動增長率云圖Fig.12 Growth rate contours for boundary layers with different sizes of cavities at Ma=1.38

由于凹腔附近局部非平行性較強，采用傳統(tǒng)的基于當?shù)鼗僭O的一維線性穩(wěn)定性理論（LST）進行穩(wěn)定性預測誤差較大，接下來我們將采用DNS 和GSD 方法研究凹腔內部流動的穩(wěn)定性（3.2.1～3.3 節(jié)內容）。但是此處我們關注的重點是凹腔的上游和下游，在稍遠離凹腔的位置流動非平行性影響比較小，圖11 和圖12 主要展示的是凹腔對不同頻率擾動在不同流向站位的影響的整體畫面，也可以給讀者帶來凹腔附近增長率畸變的宏觀感受。

3.2.1 凹腔位置影響

圖13 展示了Ma=5.92 來流下擾動同步點附近兩個頻率擾動波（ω=0.53、1.07）經過凹腔d2h1H 的擾動演化瞬時壓力場。可以看到擾動波穿過凹腔前緣的膨脹波時發(fā)生了強烈的畸變，且高頻擾動能更快地恢復到壓縮波前的狀態(tài)。擾動波能直接傳入凹腔內部，但對于同一流向站位，凹腔內擾動與凹腔外擾動存在約180o相位差。

圖13 Ma=5.92 來流下不同頻率擾動波經過凹腔d2h1H 的瞬時壓力脈動場Fig.13 Instantaneous pressure fluctuations around d2h1H at Ma=5.92 for disturbance waves with different frequencies

圖14 展示了Ma=1.38 來流下最不穩(wěn)定頻率附近兩個擾動波（ω=0.07、0.12）經過凹腔d2h1S 的擾動演化瞬時壓力場。與Ma=5.92 情況類似，擾動在經過凹腔上方時向勢流區(qū)傳播，但相比于高超聲速情況，膨脹波和壓縮波波角更大，壓縮效果更弱。

圖14 Ma=1.38 來流下不同頻率擾動波經過凹腔d2h1S 的瞬時壓力脈動場Fig.14 Instantaneous pressure fluctuations around d2h1S at Ma=1.38 for disturbance waves with different frequencies

圖15 進一步展示了高超聲速算例d2h1H 中擾動在帶凹腔平板上的幅值演化曲線。圖中綠色實線為DNS 計算得到的壓力脈動包絡，粉色虛線為LPSE 預測結果，灰色陰影區(qū)對應凹腔安裝位置。由圖可見，線性拋物化方程（LPSE）在擾動的初始發(fā)展階段能較好地預測擾動幅值沿流向的增長。對DNS 獲得的壓力脈動采用Hilbert 函數(shù)提取幅值包絡[39]，可以發(fā)現(xiàn)，擾動在凹腔前呈指數(shù)增長，在凹腔處產生復雜的波動，在凹腔下游又迅速恢復至近似指數(shù)增長。這說明凹腔處局部突變對擾動波的影響范圍有限。為此，我們引入Wu 等[12]提出的局部散射理論中的透射系數(shù)?來定量描述凹腔對擾動波和轉捩過程的影響：

其中，xf是凹腔下游透射系數(shù) ?趨于常數(shù)定值時的流向站位，Ap(xf) 是xf站位處壓力脈動幅值的法向最大值，下標“0”對應光滑無凹腔情形。當 ?＜1時，表示凹腔對擾動波起抑制作用，且其值越接近于0 抑制作用越強；當 ?＞1時，表示凹腔對擾動波起促進作用，其值越大表明促進效果越明顯。

小尺度凹腔情形的透射系數(shù)與擾動的增長率云圖可以直接關聯(lián)。在凹腔尺寸比較大時，此時凹腔下游擾動增長率誤差較大，透射系數(shù)僅能依靠擾動演化的DNS 結果換算出來。

圖16 顯示了高超聲速工況下凹腔同步點附近三個不同頻率擾動波的透射系數(shù)隨凹腔尺寸的變化情況?？梢?，不同尺寸的凹腔在同步頻率ω=0.98 附近對擾動波的抑制效果最佳。凹腔對擾動波的透射系數(shù)與凹腔寬度和深度并非呈完全單調。觀察不同寬度凹腔的透射系數(shù)分布可以發(fā)現(xiàn)，隨著寬度的增加，凹腔對其上游同步頻率擾動的促進作用逐漸增強，而對凹腔附近同步頻率擾動的抑制作用也逐漸增強。并且對于大寬度、小深度的寬淺腔d4h0.5H 和d6h1.5H算例，其在同步點附近的抑制效果更為突出。圖17展示了Ma=1.38 的三個典型工況d2h1S、d4h3S 和d6h2S算例對應的透射系數(shù)隨著擾動波頻率的變化。在所考察的頻率范圍內，凹腔一直促進擾動波的增長，并且促進效果隨頻率升高逐漸增強。這一結果與高超聲速情況下凹腔對低頻擾動的作用效果類似（圖16）。

圖16 Ma=5.92 來流下凹腔同步點附近不同頻率擾動波的透射系數(shù)隨凹腔尺寸的變化Fig.16 Variations of correction coefficients with different sizes of cavities at Ma=5.92 for disturbance waves with different frequencies near the synchronization point of the cavity

圖17 Ma=1.38 來流下三個典型工況對應的透射系數(shù)隨擾動波頻率的變化Fig.17 Variations of correction coefficients with different disturbance wave frequencies at Ma=1.38 for three typical cases

考慮到寬度和深度變化時流動特征變化顯著，因此需要進一步進行尺寸參數(shù)化研究。

3.2.2 凹腔尺寸影響

針對Ma=5.92 來流條件，我們將凹腔中心點固定在頻率為ω=0.98 擾動波的同步點位置（x=185），進一步研究凹腔的尺寸變化對擾動演化的影響。圖18 展示了ω=0.98 的入口擾動波經過表2 中4 個典型凹腔的瞬時壓力脈動場?？梢钥吹?，隨著凹腔深度和寬度的增加，擾動波穿透進凹腔并形成復雜的波系結構。擾動波在凹腔前緣被拉伸，而在凹腔后緣則被壓縮，并沿凹腔底部法向形成反射波。當凹腔寬度顯著大于擾動波波長時（如圖18（c）），可以觀察到法向反射波基本垂直于腔體壁面，并與入射波形成90°波角。

圖18 Ma=5.92 來流下擾動波（ω=0.98）經過不同尺度凹腔的瞬時壓力脈動場Fig.18 Instantaneous pressure fluctuations for the disturbance wave with ω=0.98 passing through cavities at Ma=5.92

圖19 展示了來流Ma=1.38 下ω=0.095 的入口擾動波經過3 個典型凹腔的瞬時壓力脈動場。結果表明，隨著凹腔深度和寬度的改變，凹腔上方勢流區(qū)傳播的擾動逐漸增強，并且相對于深度的增加，寬度的改變對下游邊界層擾動傳播影響更大。

圖19 Ma=1.38 來流下擾動波（ω=0.095）經過不同尺度凹腔的瞬時壓力脈動場Fig.19 Instantaneous pressure fluctuations for the disturbance wave with ω=0.095 passing through cavities at Ma=1.38

由于凹腔位于同步點位置附近時抑制效果最佳，選擇展示Ma=5.92、ω=0.98 工況下不同尺寸凹腔的擾動系數(shù)云圖（圖20），圖中星號表示4 個典型凹腔算例：d2h1H、d2h3H、d6h1H 和d6h3H?？梢赃M一步將參數(shù)區(qū)間分成三個區(qū)域：左側區(qū)域（d/δ99＜2）的小寬度凹腔抑制效果最弱；右下角區(qū)域（d/δ99＞2，h/δ99＜1.5）的大寬度小深度凹腔抑制效果最好；右上角的大寬度大深度凹腔區(qū)域次之。圖中最優(yōu)凹腔尺寸為寬淺腔d6h0.5H 算例，達到了近40%的擾動幅值抑制效果?？傮w而言，凹腔深度不變，抑制擾動幅值的效果隨凹腔寬度單調增強；而固定寬度增加深度的抑制效果并非單調變化的，這與較小尺度下的凹腔深度變化規(guī)律類似[13]。當凹腔深度增加時，腔內分離泡充滿整個區(qū)域，繼續(xù)增加深度反而產生了額外的促進作用。

圖20 Ma=5.92、ω=0.98 工況下不同尺寸凹腔透射系數(shù)云圖Fig.20 Correction coefficient contours for cavities with different sizes at Ma=5.92 and ω=0.98

對于高超聲速來流，擾動波頻率ω=0.98 時，擾動透射系數(shù)隨凹腔寬度d變化規(guī)律如圖21（a）所示：隨著寬度逐漸增加，擾動波抑制效果逐漸增強，且在(d/δ99,h/δ99)=(6,1)時達到最佳擾動幅值抑制效果近35%。擾動波抑制效果隨凹腔深度h變化如圖21（b）所示：寬度較小時，隨著深度的增加，抑制效果變化不明顯；當寬度d達到4 倍當?shù)剡吔鐚雍穸葧r，隨著深度增加，抑制效果出現(xiàn)較大的非單調波動。對于超聲速來流，擾動波頻率ω=0.095 時，隨著凹腔寬度的增加，促進效果逐漸增強（圖21（c））；隨著凹腔深度的增加，促進效果變化不明顯，對于寬腔（d/δ99=4 和6）甚至出現(xiàn)了促進效果減弱的現(xiàn)象（圖21（d））。最終研究發(fā)現(xiàn)，大寬度凹腔對擾動波幅值的控制效果（抑制或促進）均優(yōu)于小寬度凹腔。

圖21 擾動波的透射系數(shù)隨凹腔尺寸變化Fig.21 Variations of correction coefficients with different sizes of cavities

3.3 凹腔自激不穩(wěn)定性

當凹腔尺度達到一定閾值后，如Ma=1.38 來流下的d6h3S 算例，凹腔內開始出現(xiàn)自激振蕩，給下游邊界層引入新的擾動源，將觸發(fā)新的轉捩機制?？紤]到凹腔內部的強非平行性，本節(jié)采用GSD[36]開展分析。采用Gavish 和Donoho[40]提出的基于漸近分析的最優(yōu)截斷秩，通過判斷奇異值最優(yōu)截斷，篩除掉非物理的模態(tài)。所獲得的特征譜如圖22 所示。圖23 展示了4 個典型的GSD 模態(tài)的空間分布（圖22 中M1、M2、M3 和M4），對應的頻率分別為ω=0.033、0.073、0.106 和0.140?？梢园l(fā)現(xiàn)不同模態(tài)擾動分布的集中區(qū)域不同，如剪切層區(qū)域主要由M1 和M2 模態(tài)主導，尾跡區(qū)主要由M1 模態(tài)主導且擾動分布類似Mack 模態(tài)。同時凹腔內部還出現(xiàn)了新的擾動形式，如M3 模態(tài)。M4 模態(tài)的擾動能量也主要集中在凹腔上方的剪切層區(qū)域，但能量幅值較弱。凹腔上方存在向下游傳播的聲波，近壁面處存在行進波，與亞聲速凹腔流動結果類似[41]。

圖22 凹腔流動d6h3S 的全局穩(wěn)定性特征譜Fig.22 Eigen spectrum for the global stability of the cavity flow d6h3S

圖23 典型主導模態(tài)形狀函數(shù)Fig.23 Modal shapes of typical dominant modes

圖24 展示了DNS 流場中凹腔尾跡區(qū)域和前緣近壁區(qū)域采樣點，即(x,y)=(196.6,0.7)和(x,y)=(179.3,-4.9)處的流向速度時間序列及對應的傅里葉頻譜。可見，凹腔的尾跡區(qū)中存在較強的單頻主導擾動，頻率為ω=0.035，與GSD 的M1 模態(tài)接近；而在凹腔的前緣近壁區(qū)除主導頻率ω=0.035 擾動外，還存在ω=0.071、0.107 和0.141 三種頻率的擾動，與GSD 模態(tài)M2、M3 和M4 接近，頻率相對誤差小于3%。這說明凹腔前緣近壁區(qū)不穩(wěn)定性較尾跡區(qū)更為復雜，由多種模態(tài)主導。

圖25 展示了通過前數(shù)階主導GSD 模態(tài)對流場快照累計積分重構的結果?？梢园l(fā)現(xiàn)重構流場與原始流場隨時間的演化過程曲線整體吻合較好，反映了凹腔流動的GSD 模態(tài)能較好地捕捉到大尺度凹腔自激不穩(wěn)定的動力學特征，可作為原始系統(tǒng)轉捩過程研究的低維降階模型。

圖25 積分重構的速度u 與原始采樣信號對比Fig.25 Reconstructed velocity u and the original sampled signal

4 結論

本文針對高超聲速和超聲速兩種典型工況下帶凹腔的平板邊界層開展了系統(tǒng)的穩(wěn)定性和擾動演化分析，探討了不同尺度凹腔對邊界層穩(wěn)定性的影響。研究發(fā)現(xiàn)，凹腔內存在多尺度渦結構，且隨凹腔尺度的增加其變得愈加復雜。

在高超聲速情況下，當凹腔位于同步點上游時會促進擾動發(fā)展，位于同步點附近或下游時則起到了抑制作用。當固定凹腔深度、增大寬度時，擾動波抑制效果逐漸增強；當固定凹腔寬度、增加深度時，擾動波抑制效果呈現(xiàn)非單調變化。超聲速情況下，凹腔在擾動波不穩(wěn)定區(qū)間內一直起促進效果，當固定凹腔深度、增大寬度時，凹腔對擾動波的促進效果逐漸增強；但當凹腔尺度增加到一定程度，腔體本身自激振蕩激發(fā)出非定常不穩(wěn)定波，給凹腔下游邊界層引入新的轉捩途徑，值得后續(xù)進一步研究。