高超聲速有攻角錐裙直接數(shù)值模擬

賴江，范召林，王乾，董思衛(wèi)，童福林，*，袁先旭

1.中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000

2.中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

3.中國空氣動力研究與發(fā)展中心，綿陽 621000

高速飛行器內外流動中激波/邊界層干擾（Shock Wave/Boundary Layer Interaction， SWBLI）現(xiàn)象廣泛存在，對氣動力熱載荷具有巨大影響，尤其在高超聲速領域，干擾區(qū)流動分離/再附會導致局部極高熱流、摩阻，可能誘發(fā)熱力學和結構方面的問題，嚴重威脅飛行安全［1-4］。因此，激波/邊界層干擾的物理機制不僅是空氣動力學中亟待解決的關鍵基礎科學問題，在航空航天工程應用領域的重要性也不言而喻。

Gaitonde［3］、Gaitonde 和Adler［5］總 結 了 激波/邊界層干擾研究的典型標模，包括引入外部激波的反射干擾，自身誘發(fā)激波的壓縮拐角、軸對稱柱裙/錐裙，以及具有三維特征的后掠拐角和單楔/雙楔等構型。針對反射激波干擾和壓縮拐角這類二維構型已有大量公開發(fā)表的風洞試驗和數(shù)值模擬研究，對激波/邊界層干擾問題的認識已取得了較大的進展。例如，Pirozzoli 和Grasso［6］獲 取 了 馬 赫 數(shù)Ma∞=2.25 反 射 激 波 干擾中的流動分離、湍流增強、激波振蕩等現(xiàn)象并開展了相應機制分析，還提供了包含速度和壓力分布、湍動能輸運的精細流場數(shù)據(jù)庫［7］。Clem?ens 和Narayanaswamy［8］、范孝華等［9］對激波低頻振蕩研究進行綜述時指出，該現(xiàn)象可能由上游湍流邊界層的脈動驅動，也可能與下游分離流固有的不穩(wěn)定性有關。干擾區(qū)低頻非定常性的來源取決于流動分離區(qū)的尺寸，較大分離中由下游機制占主導，而較弱分離中則是上游和下游機制同時存在。Fang 等［10］對湍流增強機制進行了綜述，總結了早期包括激波非定常運動、激波/湍流直接干擾、聲波和熵波的產生、自由剪切層等各種可能的原因。隨著高速試驗技術和數(shù)值模擬方法的發(fā)展，自由剪切層被普遍認為是SWBLI 中湍流增強的主要因素。童福林等［11］通過分析比較3 個不同展向站位分離泡的非定常運動特性、分離微團幾何特征和相干結構等，定量給出了三維展向結構差異的影響規(guī)律。

隨著高超聲速飛行器的發(fā)展，對探明激波低頻振蕩原因、預測流動分離特征、發(fā)展減阻降熱技術等方面的需求越發(fā)強烈和迫切，但目前對上述部分問題認識尚不充分或還存在爭議。尤其是現(xiàn)有大多研究針對二維構型，必然會忽略三維特征，在數(shù)值模擬中還無法避免展向計算域和邊界條件等導致的誤差，從二維構型研究得到的結論在三維情況下的有效性和適用性也尚不明確，已不能涵蓋工程實際中面臨的問題和困難，因此需要面向更加接近工程實際的構型開展相關研究。

錐裙構型是SWBLI 的研究標模之一［3］，是高超聲速再入飛行器和導彈等的有效簡化，是本文所選取的計算構型。根據(jù)上游邊界層狀態(tài)，錐裙拐角附近流場可分為激波與層流邊界層干擾、激波與轉捩/湍流邊界層干擾。針對第1 類問題，大量試驗研究以25°～55°為典型幾何參數(shù)，著重刻畫干擾流場結構特征。Holden［12］歸納相關試驗研究，重點討論了Ma∞=10～12，雷諾數(shù)Re∞≈104～105條件下干擾區(qū)的熱流和壓力分布。Wright［13］、Nompelis［14］等刻畫了以激 波系統(tǒng) 為主的干擾流場，包含前錐上游的斜激波、分離激波，裙結構上的輸運激波、再附激波、弓形激波等，還給出了分離區(qū)、接觸面、超聲速噴流、聲速線、三叉點等典型結構。Candler［15］、Tumuklu［16］等采用蒙特卡洛方法獲取流場數(shù)據(jù)，詳細描述了拐角區(qū)域干擾流場的主分離泡和次分離泡的精細結構。針對第2 類問題，轉捩/湍流邊界層導致干擾流場更加復雜。Holden 等［17-18］針對6°前錐湍流錐裙模型，在高馬赫數(shù)Ma∞=11～16、高雷諾數(shù)Re∞≈106條件下開展了試驗研究，通過紋影、全息等一系列手段測量了壁面壓力、摩阻、熱流信息，并與理論預測值進行對比，但沒有詳細評估各壁面量之間的關系。岡敦殿［19］改變半錐角和前錐長度，測試了Ma∞=3 條件下多個錐裙的上游流態(tài)，建立了激波與湍流邊界層干擾流場，觀察到拐角處的分離流動和下游大尺度結構，并獲取5°～40°錐裙的時間相關納米粒子示蹤的平面激光散射圖像，發(fā)現(xiàn)軸對稱流場存在明顯的周向運動，但沒有對流動分離的非定常運動開展分析。Running等［20-21］對Ma∞=6 條件下7°前錐的錐裙結構開展風洞試驗，根據(jù)上游邊界層呈層流、轉捩、湍流3 種狀態(tài)，評估了頭部鈍度、裙角度、來流雷諾數(shù)等參數(shù)變化對邊界層分離、再附的影響；他們還采用紅外測熱技術定量分析了熱通量條紋的尺度和周向特征［22］，并利用高頻響應的壓敏漆技術成功捕捉到激波低頻振蕩現(xiàn)象［23］。Shiplyuk［24］、Bedarev等［25］對Ma∞=5.92 的7°～17°錐裙開展了試驗和計算研究，驗證了不同湍流模型對轉捩/湍流情況下的壓力分布和熱載荷預測的準確度，F(xiàn)rayssinet等［26］進一步獲取了該模型在1°攻角下的壓力和熱流分布，初步刻畫了分離點位置沿迎風區(qū)向背風區(qū)的演化，但并未針對不同周向子午面的壁面量關系和流動分離情況進行分析。

現(xiàn)有高超聲速錐裙研究主要采用風洞試驗測量和雷諾平均模擬、蒙特卡洛方法等數(shù)值手段，但采用高精度數(shù)值模擬的研究極為少見。筆者前期已針對軸對稱錐裙開展了國內首個激波/湍流邊界層干擾直接數(shù)值模擬（Direct Numerical Simula?tion， DNS）研究［27］，但計算域僅局限于有限的周向范圍并采用周期性邊界條件。因此，基于周向平均結果的分析不可避免地忽視了三維橫流效應的影響。與此不同的是，本文對錐裙整個360°周向進行直接數(shù)值模擬，研究了攻角對干擾區(qū)流動的影響規(guī)律，重點關注0°、90°、180°這3 個周向子午面的流場結構和壁面量，采用功率譜分析、本征正交分解方法研究干擾區(qū)分離泡非定常運動特性，并比對三維橫流效應下再附邊界層演化過程的差異。

1 直接數(shù)值模擬方法

1.1 計算模型和網(wǎng)格

本 文 計 算 模 型 與Running 等［20］在AFRL 的馬赫6 路德維希管風洞試驗中所采用的錐裙構型一致，由7°半錐角的前錐和34°裙結構組成，其長度分別為609 mm 和76 mm。圖1 給出了計算域示意圖，定義z為軸向，r為徑向，φ為周向。將原點設置在前錐頂點處，記為z=0 mm。類似于Sivasubramanian 和Fasel［28］、Huang 等［29］對 尖 錐邊界層轉捩進行DNS 研究的處理，模擬區(qū)域不包含尖錐頭部，而是將計算域入口設置在尖錐頂點下游z=156 mm 處。計算域出口位于裙結構的尾部。計算域流向長度為517 mm，法向高度為40 mm，周向角度為360°。為了使入口層流發(fā)生轉捩，在壁面za=180～185 mm 范圍設置定常吹吸擾動，參考Pirozzoli 等［30］的設置，隨機法向速度擾動的形式表示為Vbs=Afbs（z）gbs（θ），其中擾動幅值A= 0.4U∞，空間函數(shù)分別定義為

圖1 計算域示意Fig.1 Sketch of computational domain

式中：Φl為0～1 的隨機數(shù)；Lφ為周向計算域寬度；θ為周向角度。

圖2 給出了計算網(wǎng)格的分布。定義x為壁面流向，y為壁面法向。圖2（a）中流向-法向截面中網(wǎng)格間隔10 個點顯示，沿壁面流向采用2 757 個網(wǎng)格點進行離散，其中上游轉捩區(qū)156 mm＜z＜576 mm 逐漸加密地分布了1 301 個網(wǎng)格點，核心干擾區(qū)576 mm＜z＜660 mm 均勻分布了1 407個網(wǎng)格點，在緩沖區(qū)660 mm＜z＜675 mm 對出口附近網(wǎng)格進行拉伸，布置了逐漸稀疏的51 個網(wǎng)格點抑制激波反射干擾。沿壁面法向采用指數(shù)加密技術使y方向網(wǎng)格逐漸向壁面聚集，設置319 個網(wǎng)格點。所有周向子午面的網(wǎng)格均與此一致。圖2（b）中法向-周向截面中網(wǎng)格間隔5 個點顯示，為了在周向采用周期條件處理，周向1 201 個網(wǎng)格點分布在θ=?90°～269°的范圍。在θ=?22.5°～205.5°的范圍內均勻分布了1 000 個網(wǎng)格點，網(wǎng)格向周向邊界逐漸變疏作為緩沖，分析取θ=0°～180°范圍結果進行后處理。在核心干擾區(qū)，均勻的流向網(wǎng)格尺度為Δx=0.06 mm，周向網(wǎng)格對應角度為θ=0.225°，但周向網(wǎng)格尺度隨流向錐裙半徑增大而增大，例如在上游參考點zref=580 mm 的尺度為Δrθ=0.280 mm，拐角zc=609 mm 的尺度為Δrθ=0.294 mm，下游zd=660 mm 的尺度為Δrθ=0.427 mm。壁面法向首層網(wǎng)格厚度為Δyw=0.01 mm。表1 列出了θ=0°，90°，180°這3 個典型周向子午面的網(wǎng)格分辨率，上標“+”表示按照當?shù)乇诿鎱⒖剂喀苔?、νw進行處理。本文中“平均”表示時間平均，對于變量η，可分解為η=ηˉ+η′和η=η?+η′′，其 中ηˉ為η的 雷 諾 平 均 值，為η的密度加權平均值，對應的脈動值分別表示為η′、η′。

圖2 網(wǎng)格分布Fig.2 Grid distribution

1.2 DNS 設置

采用高精度有限差分代碼OpenCFD-SC，直接求解一般曲線坐標系下的完全氣體三維守恒型可壓縮Navier-Stokes 方程。該求解器已廣泛應用于可壓縮湍流計算，大量反射激波干擾和壓縮拐角研究證實其在激波/邊界層干擾問題計算精度和可靠性，對升力體、鈍錐等復雜構型的計算也能夠得到令人滿意的結果。因此其模擬能力滿足本文對錐裙結構的計算需求。采用來流參數(shù)對控制方程進行無量綱處理，參考長度為L=1 mm。黏性系數(shù)由Sutherland 公式計算，壓力和溫度關系滿足無量綱理想氣體狀態(tài)方程，Prandtl 數(shù)、比熱比分別為Pr=0.7、γ=1.4。計算條件如下：馬赫數(shù)為Ma∞=6.0，單位雷諾數(shù)為Re∞=10 800 mm?1，靜溫為T∞=65 K。除非另做說明，下標“∞”和“w”分別表示來流、壁面量。由于大攻角會導致前錐背風區(qū)流動分離，而小攻角沒有明顯三維橫流影響，為滿足對錐裙拐角區(qū)域SWBLI 的研究需求，這里將攻角設為α=5°。

數(shù)值方法方面，空間對流項離散采用由Martín 和Taylor［31］、Wu 等［32］提 出 的 帶 寬 優(yōu) 化WENO-SYMBO 格式和限制器，黏性項由八階中心差分格式處理。時間推進采用三階總變差衰減的Runge-Kutta 方法，時間步長為ΔtU∞/L=0.005，U∞為來流速度。對于邊界條件，計算域入口為層流邊界層，層流邊界層剖面是通過相同來流條件的三維層流尖錐模擬得到，圖3 給出了層流入口的密度、流向速度、溫度的分布云圖。出口為超聲速條件。計算域底部是模型壁面，為無滑移等溫壁條件，壁面溫度為Tw=296.4 K，與恢復溫度Tr的比值為Tw/Tr= 0.62，視為冷壁面。圖1 所示吹吸擾動區(qū)域設置定常吹吸擾動，以較大幅值ψ= 0.4U∞加速尖錐轉捩過程。此外，上邊界設置無反射條件［33］以抑制偽擾動的反射。周向邊界設置為周期性條件。

圖3 層流入口Fig.3 Laminar inlet

2 典型流場結構

圖4 給出了瞬時速度和時均溫度分布云圖，圖中-T為無量綱平均溫度。圖4（a）、圖4（b）通過三維視圖展示了前錐上游流場從有序到雜亂的過程，表明入口層流在吹吸擾動作用下發(fā)生了轉捩。圖4（c）、圖4（d）為沿錐裙中心軸從前向后的6 個典型流向截面二維視圖，第1 截面為入口層流狀態(tài)；第2 截面在背風中心線附近存在“蘑菇狀”結構，沿周向存在向兩邊分別發(fā)展了橫流渦結構，即出現(xiàn)轉捩過程；隨著吹吸擾動的持續(xù)作用和邊界層的進一步發(fā)展，第3、4 截面中流動結構逐漸破碎；經(jīng)過錐裙拐角干擾區(qū)后，第5、6 截面中流動已完全混亂。注意到在錐裙拐角后出現(xiàn)了激波/邊界層干擾典型現(xiàn)象，如圖4（c）中存在大范圍速度為負的流動分離區(qū)域，圖4（d）中出現(xiàn)了靠近壁面的溫度升高現(xiàn)象。

圖4 瞬態(tài)和時均流場Fig.4 Instantaneous and time-averaged flowfield

圖5 典型周向子午面瞬時密度梯度場Fig.5 Instantaneous density gradient flowfield at typi?cal cross sections

存在攻角的情況下，圖2（a）網(wǎng)格分布中所示的θ= 90°周向子午面中無黏流線與當?shù)厮俣戎g存在夾角，產生三維橫流，其附近周向區(qū)域在后文分析中被稱作橫流區(qū)。幾何上迎風中心線（θ= 180°）、背風中心線（θ= 0°）均與當?shù)亓骶€重合，附近周向區(qū)域分別被稱作迎風區(qū)、背風區(qū)。通常利用99%來流速度位置判斷邊界層厚度，但這里存在壓力梯度，因此以來流總焓的1.003 倍位置定義邊界層厚度。經(jīng)估算，迎風區(qū)、橫流區(qū)、背風區(qū)在上游參考點處邊界層厚度依次約為δref= 2.84， 3.63， 13.31 mm，表明邊界層厚度從迎風區(qū)到背風區(qū)經(jīng)歷了逐漸增厚的發(fā)展過程。對比0°攻角錐裙相應位置邊界層厚度（δref=5.88 mm），5°攻角下邊界層的發(fā)展在迎風區(qū)、橫流區(qū)被抑制，在背風區(qū)卻大幅增厚。，其 中?ρ為 密度梯度，下標“max”和“min”分別表示其最大、最小值，圖中顏色越深的區(qū)域表明密度梯度越大。相較于Tong 等［27］對軸對稱錐裙的研究，有攻角情況下激波強度變弱，并形成了更加復雜的相互作用波系。值得注意的是，背風中心線截面中難以辨識出規(guī)則的斜激波結構，這將導致背風區(qū)與其他周向區(qū)域的差異，其拐角處的流動分離可能主要由構型的變化導致，但未形成分離激波和再附激波。

3 平均壁面量分布

圖6 給出了時間平均壓力、摩阻和熱流的壁面分布三維視圖。沿流向觀察到干擾導致的各壁面量顯著增強現(xiàn)象。圖6（a）中前錐的壁面壓力水平遠小于裙結構上的值，在錐裙拐角前已出現(xiàn)壓力上升。在整個模型的干擾前后，壓力從迎風區(qū)經(jīng)橫流區(qū)至背風區(qū)逐漸降低，在拐角前的壓力劇增出現(xiàn)的位置逐漸后移。圖6（b）中壁面摩阻在前錐上游部分保持在較低水平，其值升高表示轉捩發(fā)生。經(jīng)過拐角后進一步增大，但僅保持一小段流向距離隨后降低。在攻角影響下，轉捩位置在迎風面相對靠前（上游），經(jīng)橫流區(qū)至背風區(qū)向后（下游）推遲，但在背風中心線附近再次前移。整體來看，壁面摩阻在前錐上保持在較低水平，進入錐裙拐角干擾區(qū)后摩阻驟降為負，表示該位置發(fā)生流動分離，此后由于流動再附再次導致摩阻劇增，并持續(xù)到模型尾部。需要注意的是，圖6（c）中壁面熱流的流向峰值區(qū)域在周向存在明顯的區(qū)別，迎風區(qū)的峰值更高。此外，背風中心線附近的分布與其他周向位置存在明顯區(qū)別，可能與該區(qū)域的不對稱小渦有關。

圖6 時均壁面量分布Fig.6 Distribution of time-averaged wall quantities

圖7 為3 個周向子午面中時均壁面壓力Pw、摩阻Cf、熱流Ch的流向分布。各壁面量在流動進入干擾區(qū)后都發(fā)生了一定程度的增加，但由于激波強度和邊界層狀態(tài)的不同，變化過程和定量程度在不同子午面上存在明顯的差異。圖7（a）中無量綱壁面壓力在錐裙拐角之前緩慢上升，在下游再附區(qū)中隨著周向位置變化存在不同變化趨勢，其中迎風區(qū)和橫流區(qū)的壓力達到峰值（約上游的10 倍），此后略微降低并分別保持，但背風區(qū)的壓力始終處于上升過程。整體來看，迎風區(qū)受到干擾最強，背風區(qū)受到的影響相對較小，橫流區(qū)變化過程與0°攻角情況最為接近，無量綱壁面壓力在上游約為0.037，經(jīng)過干擾最終保持在0.49附近。

圖7 典型周向子午面平均壁面量流向分布Fig.7 Streamwise distribution of time-averaged wall quantities at typical spanwise cross sections

圖7（b）中壁面摩阻在干擾區(qū)降低，在分離區(qū)為負值，此后由于再附而劇增并達到峰值。以zref處的壁面量為參考，經(jīng)過干擾區(qū)在背風、橫流、迎風區(qū)再附邊界層中摩阻最高增至原來的2.60、2.08、 1.72 倍（遠低于0°攻角的4.72 倍左右）。從迎風區(qū)到橫流區(qū)，以Cf＜ 0 定義的流向分離尺度（顯著大于0°攻角的4 mm 左右）逐漸增大，從23.7 mm 擴大至26.3 mm，但背風中心線處流向分離區(qū)域僅11.8 mm，這是由于背風區(qū)近壁分離泡可能并不與從迎風面發(fā)展的分離泡相連。圖7（c）中壁面熱流在流動分離處略低，然后緩慢增加，在拐角處劇增并在再附點附近達到峰值（約為上游的5 倍），此后逐漸降低但始終高于上游，各子午面上相對大小與上游一致。該現(xiàn)象不同于0°攻角情況，熱流在軸對稱錐裙拐角處劇增后經(jīng)歷了非常緩慢的下降過程。

以上壁面量流向分布特征表明，熱流的變化與壓力和摩阻均存在一定的關系，與流動分離/再附、拐角構型變化相關。因此，圖8、圖9 分別給出了QP85 關系式［34］和雷諾比擬因子［35］在各子午面上沿流向的分布。圖中，S 為分離點，C 為錐裙拐點，R 為再附點，下標0、90、180 表示圖2（b）所示的周向位置對應角度。

圖8 QP85 流向分布Fig.8 Streamwise distribution of QP85

圖9 雷諾比擬因子流向分布Fig.9 Streamwise distribution of RAF

平均熱流與壓力關系定義為

式中：下標“0”表示上游參考值。

如圖8 所示，各周向子午面中上游基本滿足QP85 ≈ 1.0，與Murray 等［36］的 試 驗 結 果 一 致。流動分離和錐裙拐角構型的變化誘發(fā)了2 次QP85 上升，在0.5～1.5 之間變化；隨后流動再附導致QP85 迅速降低，圖8（a）中迎風區(qū)和橫流區(qū)下游保持在0.3 左右。以上變化趨勢在圖8（b）的背風區(qū)中涉及的流場范圍更大，直到下游區(qū)域仍處于下降狀態(tài)。再附區(qū)的單調遞減趨勢與Priebe和Martín［37］在高超聲速壓縮拐角中的現(xiàn)象一致，與圖7（a）中壓力的持續(xù)上升有關。

平均熱流與摩阻的雷諾比擬關系定義為

4 分離泡非定常特性

由圖9 可知，RAF 在分離區(qū)完全失效，這是由分離區(qū)存在極小負值的摩阻所致。而圖9（a）中，RAF 在迎風區(qū)、橫流區(qū)的無干擾上游、下游恢復區(qū)均符合較好，其中上游RAF ≈ 1.07，與0°攻角 錐 裙 的 結 果 一 致，也 符 合Roy 和Blottner［35］在高超聲速零壓梯度邊界層中給出的0.9～1.3 范圍。下游則穩(wěn)定在RAF≈1.20，在流向尾部出現(xiàn)小幅振蕩，與0°攻角下逐漸增大的變化趨勢略有不同。RAF 在圖9（b）的背風區(qū)上游、下游都高于其他位置情況，上升、下降趨勢所涉及的流向范圍更廣，說明其干擾存在明顯不同。由此可知，在三維橫流影響下，迎風區(qū)、橫流區(qū)的無干擾上游、下游恢復區(qū)仍然可用QP85、RAF 對壓力、摩阻、熱流進行預測和評估。

根據(jù)在圖4 觀察到的流動分離現(xiàn)象，圖10 進一步提取θ=0°， 90°， 180°這3 個周向子午面內錐裙拐角附近的干擾流場。在圖10（a）～圖10（c）給出的瞬態(tài)速度場（以來流速度U∞對瞬時速度u進行無量綱處理）中，將貼體流向速度小于0 的區(qū)域看作分離泡，用粉色點劃線表示時均分離泡邊緣，發(fā)現(xiàn)背風子午面的流動分離區(qū)域相對較小，但速度變化的法向范圍相對最大。

圖10 典型周向子午面的流動分離Fig.10 Flow separation at typical spanwise cross sections

圖10（d）～圖10（f）的時均溫度場（以來流溫度T∞對時均速度-T進行無量綱處理）中疊加繪制了黑色流線。各子午面均出現(xiàn)流向較長、法向較短的分離泡結構，且溫度較高的區(qū)域與回流區(qū)范圍一致。然而，從迎風區(qū)經(jīng)橫流區(qū)到背風區(qū)，溫度變化的法向范圍逐漸變大，尤其是在背風子午面中該范圍遠超回流區(qū)，推測這與從迎風區(qū)到背風區(qū)逐漸增大的邊界層厚度有關。注意到流線表征的回流區(qū)中，背風子午面的分離泡較為分散，而橫流區(qū)、迎風區(qū)表現(xiàn)為一個較為整體的分離泡。這可能是由于在攻角影響下，背風區(qū)結構不完全對稱，且背風中心線附近的近壁渦與其他周向位置不同的緣故。

圖7（b）給出的壁面摩阻流向分布、圖10（d）～圖10（f）所示平均流線趨勢明確揭示了錐裙拐角區(qū)域干擾流場的流動分離現(xiàn)象。通過提取干擾區(qū)中貼體流向速度小于0 區(qū)域，對分離泡面積進行統(tǒng)計，圖11 給出了以其最大值進行無量綱處理的面積脈動A′（t）隨時間的變化歷程，觀察到分離泡在不斷進行膨脹和收縮運動。圖10（a）～圖10（c）所示瞬態(tài)速度場中，貼體流向速度的負值范圍在平均分離泡邊緣之外，表明該時刻的分離泡正處于膨脹運動。圖11 繪制了以最大值進行無量綱處理的平均分離點S 位置的壁面壓力脈動隨時間t的變化歷程。

圖11 分離泡面積脈動平均分離點壓力脈動時間序列Fig.11 Time series of separation bubble area fluctuation and pressure fluctuation at mean separation point

為了進一步揭示干擾區(qū)流動分離與激波之間的關系，圖12 給出了分離泡面積脈動A′（t）與平均分離點S 的壁面壓力脈動p′w( )t之間的相關性系數(shù)CPA結果，其定義為

圖12 分離泡面積和平均分離點壓力的互相關系數(shù)Fig.12 Cross-correlation coefficient between separation bubble area and wall pressure at mean separation point

式中：Δt為時間間隔，橫流區(qū)子午面在Δt= 0 處達到峰值CPA≈ 0.45，即干擾區(qū)的激波低頻振蕩與分離泡的非定常運動密切相關。迎風區(qū)子午面中兩者也具有強相關性，但觀察到一定的時間延遲。而背風區(qū)子午面中兩者無關，這是由于其激波結構不明顯，因而流態(tài)區(qū)別較大。

圖13 給出了θ= 0°， 90°， 180° 3 個周向子午面上分離泡面積序列的預乘功率譜，進一步對干擾區(qū)頻域特征開展分析，其中，ω為角頻率，Φ（ω）為功率譜。這里采用了Welch 方法和Hamming窗，將數(shù)據(jù)分為4 段，重疊率設置為50%。Strou?hal 數(shù)定義為St=fL/U∞，其中f為頻率。顯然，各子午面上分離泡面積脈動的主頻均為St=0.009 7，背風區(qū)還在St= 0.029 存在局部峰值，可能與圖10（d）中觀察到回流區(qū)的分散的多泡結構相關。

圖13 分離泡面積脈動預乘功率譜Fig.13 Pre-multiplied spectra of separation bubble area fluctuation

根據(jù)圖11 中給出的各周向子午面上平均分離點位置的壓力時間序列，圖14 給出其預乘功率譜，并以脈動均方根的平方進行處理，同時給出zref的結果作為參考。各子午面的上游壓力脈動主頻均在高頻區(qū)，即背風區(qū)、橫流區(qū)、迎風區(qū)中zref站位的壓力脈動主頻分別為St=0.22，0.45，0.71。從干擾區(qū)站位S 結果來看，各子午面預乘譜與上游的分布趨勢相似，全局主頻略微向低頻區(qū)移動，分別降至St=0.14，0.25，0.35。重要的是，在迎風區(qū)、橫流區(qū)，分離點S 的預乘功率譜均在低頻區(qū)St=0.009 7 處存在局部峰值，該頻率與圖13 中分離泡面積脈動的主頻一致，相應時間尺度分別為上游的73.2、46.4 倍。與以往SWBLI 研究中典型的時間尺度在上游的10～100 倍之間的分離激波低頻振蕩現(xiàn)象相符。

圖14 平均分離點壁面壓力脈動預乘功率譜Fig.14 Pre-multiplied spectra of wall pressure fluctua?tion at mean separation point

對圖15 的預乘功率譜結果按照

圖15 平均分離點壁面壓力脈動預乘功率譜積分Fig.15 Integral of the pre-multiplied spectra of wall pressure fluctuation

進行積分并采用頻率最大值結果歸一化處理，得到各頻率能量分布。將St＜ 0.01 定義為低頻區(qū)，圖15（a）的背風區(qū)中，上游和分離點S 的低頻段能量占比分別為17.3%、11.0%。而圖15（b）中迎風區(qū)、橫流區(qū)則出現(xiàn)完全相反的結果，分離點S 的低頻能量占比劇增。具體而言，分別從參考點處的0.8%、2.4% 增至分離點S 處的13.5%、19.0%。因此，頻域分析結果進一步確定了迎風區(qū)、橫流區(qū)分離激波的低頻特征，但是背風區(qū)不存在明顯的激波結構（見圖5（a）），因此僅觀察到分離點S 處壓力脈動與上游一致的高頻特征。

為了進一步揭示分離泡非定常運動中的典型相干結構，采用本征正交分解（Proper Orthogonal Decomposition， POD）對 流 向574 mm ＜z＜625 mm、法向約0 mm＜r＜15 mm 采樣窗口區(qū)的瞬態(tài)流向速度場進行分析。采樣時間步長ΔtU∞/L= 0.4，共獲取N= 2 500 個POD 模態(tài)。圖16 給出了模態(tài)能量Ek分布和累積能量占比。圖16（a）表明Ek隨著模態(tài)階數(shù)k的增加而逐漸降低，例如第10、100 階模態(tài)的能量分別比第1 階模態(tài)下降約1、2 個數(shù)量級。圖16（b）給出的累積能量揭示第1 階模態(tài)能量占比最高，約20%左右，視作主能量模態(tài)。注意到各子午面的能量變化趨勢雖然沒有顯著差異，但圖16（a）表明背風區(qū)能量在低階模態(tài)中相對最低，而在高階模態(tài)中相對最高。圖16（b）證實了上述規(guī)律，背風區(qū)的累積能量在前100 階模態(tài)都是相對最低的。

圖16 POD 模態(tài)能量分布Fig.16 Energy distribution of POD modes

POD 模態(tài)時間系數(shù)ak（t）反映了采樣窗口區(qū)干擾流場的非定常特征。圖17（a）中ak（t）預乘功率譜云圖表明，占優(yōu)頻率隨著模態(tài)階數(shù)增加從低頻區(qū)迅速轉向高頻區(qū)，在k＜100 都低于上游壓力脈動主頻。具體而言，k＜5 低階模態(tài)的峰值頻率集中在低頻范圍，在分離泡面積脈動特征頻率附近，但橫流區(qū)主能量模態(tài)主頻略有升高。5 ＜k＜10 的迎風區(qū)峰值頻率也略有升高。圖17（b）為ak（t）和分離泡面積脈動的相關系數(shù)CMA隨模態(tài)階數(shù)的變化，定義為

圖17 模態(tài)時間系數(shù)的預乘功率譜和互相關系數(shù)Fig.17 Pre-multiplied spectra and cross-corre-lation coefficient of time coefficient of POD modes

整體來看，主能量模態(tài)均與分離泡具有較高的相關性，此后低階模態(tài)相關性隨階數(shù)增大而急劇衰減，高階模態(tài)的CMA維持在零附近。因此，分離泡低頻膨脹/收縮運動僅與主能量模態(tài)和第2階模態(tài)（迎風區(qū)）密切相關，這與在反射激波干擾問題中得到的POD 分析結果存在顯著區(qū)別，以往分析得到的分離泡呼吸與前10 階POD 模態(tài)都存在顯著的相關性。

圖18 給出了第1、2 階POD 模態(tài)在部分法向范圍的空間分布，并疊加顯示了平均分離線位置。各子午面中，主能量模態(tài)均對應靠近壁面的大尺度含能結構，出現(xiàn)在流動分離剪切層附近。第2 階模態(tài)的主要結構分布范圍基本沒有變化，但表現(xiàn)流向尺度減小并具有正負交替特征的多個結構。對比不同周向子午面來看，從迎風區(qū)到背風區(qū)，主要空間結構的法向范圍逐漸擴大，與逐漸增厚的邊界層變化相關。

圖18 POD 模態(tài)空間分布Fig.18 Spatial distribution of POD modes

選取前10 個低階POD 模態(tài)對干擾區(qū)速度場進行重構，圖19 給出了由分離泡面積最大值Amax進行無量綱處理的分離泡面積A時間序列。對比來看，POD 重構曲線與DNS 原始曲線吻合較好，重構的面積變化準確捕捉了分離泡膨脹、收縮運動的交替變化過程。評估重構精度均達到87%以上，證明低階模態(tài)在該區(qū)域占主導地位，與以往在反射激波干擾、壓縮拐角等研究中得到的結論類似。兩者之間存在的微小差異體現(xiàn)在局部高頻脈動部分，這是由于缺少高階模態(tài)的高頻區(qū)能量所導致的。

圖19 分離泡面積時間序列的重構Fig.19 Reconstruction of time series of separation bubble area

5 再附邊界層演化

為了更好地理解再附邊界層的恢復過程，進一步刻畫時均速度、溫度、湍動能、雷諾應力分量及其各向異性不變量等特征統(tǒng)計量的流向演化，并通過對比θ=0°，90°，180°這3 個周向子午面的異同評估三維橫流效應對干擾恢復過程的影響。所選取的流向站位分別記為R1～R5，相應流向坐標依次為z=610，617，626，635，650 mm。其中R1 位于分離泡內，R2 位于再附點后，R3～R5均位于下游再附邊界層中，同時給出了各子午面上游zref站位的結果作為參考。

圖20 為壁面貼體流向速度和溫度的法向分布，圖中yn為距離壁面的距離。經(jīng)過干擾，速度沿流向先減小然后逐漸增大恢復，而溫度沿流向先升高再逐漸降低恢復。上述過程在迎風區(qū)、橫流區(qū)子午面比在背風區(qū)子午面更迅速。具體而言，圖20（a）～圖20（c）中背風區(qū)的上游速度剖面發(fā)生變形，與迎風區(qū)、橫流區(qū)較為飽滿的速度剖面存在顯著差異。R1 出現(xiàn)負速度，這與其存在于回流區(qū)的情況是對應的，R2 位于再附點下游不遠處，近壁區(qū)受到的干擾較大，R3、R4 在非常近壁區(qū)迅速接近于上游情況。但在相對外層區(qū)域，橫流區(qū)的R3～R5 存在2 段速度分別在0.80U∞、0.85U∞的穩(wěn)定區(qū)域。各子午面中R1～R5 的法向分布都沒有與zref達到一致，而在0°攻角下邊界層外緣基本可以恢復到與來流一致的水平。圖20（d）～圖20（f）中溫度變化趨勢與圖20（a）～圖20（c）是類似的，R1、R2 受到干擾影響溫度出現(xiàn)明顯的升高，在內層大幅度偏離，在外層逐漸靠近上游參考值。雖然R3～R5 站位在橫流區(qū)、迎風區(qū)的近壁區(qū)已開始恢復，且邊界層外緣的分布基本與上游剖面重合，但背風區(qū)的整體恢復相對較慢，直到R5 站位還保持著顯著的溫度升高。

圖20 流向速度和溫度法向分布的演化Fig.20 Evolution of normal distribution of time-averaged streamwise velocity and temperature

將湍動能k進行無量綱處理為

式中：密度ρ、摩擦速度uτ均采用當?shù)貐?shù)。

圖21 給出了各子午面內的法向分布。由于5°攻角帶來三維橫流和周向梯度變化，各子午面在zref站位存在明顯不同。整體來看，迎風區(qū)、橫流區(qū)、背風區(qū)的全局峰值分別為k+=6.43，7.31， 13.43，高于0°攻角下的k+= 5.2。其法向位置出現(xiàn)在內層y+n≈ 20～30，略高于Pirozzoli 和Grasso［6］在零壓力梯度可壓縮湍流邊界層中給出的y+n≈ 14。外層y+

圖21 湍動能法向分布Fig.21 Normal distribution of TKE

n≈ 200 附近還出現(xiàn)了局部峰值。湍動能在背風區(qū)外層仍然較高，而迎風區(qū)、橫流區(qū)法向分布規(guī)律與以往結果一致。激波干擾對下游區(qū)域的影響主要體現(xiàn)在2 個方面。首先是湍動能被增強，R2～R5 站位的湍動能法向分布曲線整體抬升，在R2 外層峰值偏離上游最多，在R3、R4 中降低并形成平臺區(qū)域，而在下游R5 的近壁區(qū)已出現(xiàn)與上游重合的部分，外層區(qū)域還存在較大差異。其次，干擾導致峰值的法向位置發(fā)生變化。各子午面的R2 峰值位置均移向外層，攻角影響下R3～R5 的背風區(qū)峰值始終在外層，而迎風區(qū)和橫流區(qū)峰值已回到內層。根據(jù)圖5 給出的瞬時密度梯度場可知，y+n～1 000 的小尖峰是激波導致的。背風區(qū)的激波影響明顯較弱，而迎風區(qū)R5 站位的小尖峰不止一個，體現(xiàn)了復雜波系中存在多個壓縮波相互作用。

圖22 給出了雷諾應力流向分量u′′u′′、法向分量v′′v′′的法向分布，圖23 給出了雷諾應力展向分量w′′w′′、剪切分量u′′v′′的法向分布，各分量均采用當?shù)啬Σ了俣萿τ進行處理。由于R1 位于分離點內，這里主要關注再附邊界層中R2～R5 站位的變化特征，圖中還給出了上游參考位置zref的結果以作對比。注意y+n～1 000 附近出現(xiàn)的尖峰現(xiàn)象表明外層區(qū)域明顯受到錐裙拐角干擾區(qū)產生的激波結果影響。

圖22 雷諾應力流向、法向分量的法向分布Fig.22 Normal distribution of streamwise and normal components of Reynolds stress

圖23 雷諾應力展向、剪切分量的法向分布Fig.23 Normal distribution of spanwise and shear components of Reynolds stress

顯而易見，經(jīng)過干擾后各分量均被增強，這與Fang 等［10］在Ma∞=2.25 的入射激波平板中和Loginov 等［38］在Ma∞=2.95 的 壓 縮 拐 角 構 型 上的現(xiàn)象一致。相對于上游參考值定義雷諾應力分量的放大因子，攻角的存在導致了再附點附近放大因子劇增。具體而言，0°攻角再附點附近（R1）4 個分量（u″u″、v″v″、w″w″、u″v″）的放大因子分別為3.0、16.0、9.0、17.0，而5°攻角下再附點附近（R2）在迎風子午面分別為3.95、22.19、12.79、8.48。對 比 圖23（a）～圖23（c）、圖23（d）～圖23（f），法向分量增大的倍數(shù)大于流向分量。與Schreyer 等［39］在高超聲速激波邊界層干擾中發(fā)現(xiàn)的法向脈動比流向脈動的放大更為顯著的結論類似。另一個重要現(xiàn)象是峰值位置顯著外移。以流向分量為例，在迎風區(qū)上游zref站位，沿壁面法向逐漸增大，在y+n=18.05 達到峰值，然后逐漸減小，該趨勢與以往文獻結果一致，但峰值出現(xiàn)的位置相較于Pirozzoli 等［40］在可壓縮湍流邊界層中得到的峰值位置y+n=13 略遠離壁面。經(jīng)過干擾，再附點附近R2 站位的全局峰值出現(xiàn)在顯著遠離壁面的位置，達到～100 量級，而原峰值位置在～10 量級保持并形成局部峰值。

在不同周向子午面內，雷諾應力分量的增強和峰值外移的基本特征沒有發(fā)生本質上的改變。但是，相對于迎風區(qū)，橫流區(qū)法向分量的放大因子有所降低，峰值位置也更加遠離壁面。背風區(qū)情況則非常不同，三維橫流效應對背風子午面各分量的外層影響最大，導致其峰值位置沿R1～R5 持續(xù)外移；尤其是對于包含法向速度脈動的法向分量和剪切分量，一方面，在橫流和迎風子午面中近壁區(qū)域存在全局或局部的峰值，而背風子午面內完全不存在內層峰值；另一方面，橫流、迎風子午面中放大因子沿R3～R5 呈逐漸降低的恢復趨勢，而背風子午面內則存在波動。最后，從干擾恢復過程來看，沿再附邊界層向下游發(fā)展，R3～R5 的峰值位置逐漸靠近壁面。R2～R5 的峰值放大因子逐漸降低，內層分布逐漸靠近上游參考值，表明再附剪切層逐漸衰減，近壁結構逐漸恢復。相較于0°攻角流向分量恢復最慢的特征，5°攻角下各子午面內流向分量的恢復卻是最為迅速，表明攻角導致的三維橫流對法向、周向速度脈動的影響強于流向速度脈動。內層恢復比外層快，直到下游R5 位置，雷諾應力各分量在所選取的周向子午面上均未完全恢復，尤其是背風區(qū)子午面中與其上游分布差異較大。

圖24 進一步給出下游R5 站位處雷諾應力張量各向異性不變量的法向分布情況，根據(jù)Lum?ley［41］給出的三角形輔助線以觀察湍流狀態(tài)變化，IIIb、IIb分別為雷諾應力各向異性張量第三、第二不變量。圖中的3 個頂點分別標注為一組分（One-Component， 1C）、兩組分軸對稱 （Two-Component Axisymmetric， 2CA）和 各 向 同 性（Isotropic， ISO）狀態(tài)，3 個邊界分別標注為軸對稱壓縮（Axisymmetric Compression， AC）、軸對稱膨脹（Axisymmetric Expansion， AE）和兩組分（One-Component， 2C）狀態(tài)。對比圖中3 個典型周向子午面情況，發(fā)現(xiàn)近壁區(qū)的湍流都呈2C 狀態(tài)，并很快達到各向異性峰值，與可壓縮湍流邊界層的近壁狀態(tài)一致，但是從迎風面到背風面各向異性的程度逐漸降低。隨著法向位置增大，湍流各向異性減弱，折返向AC 狀態(tài)靠近，但橫流區(qū)更接近AE 狀態(tài)，而迎風區(qū)、背風區(qū)更靠近2C 狀態(tài)。此外，隨著法向位置進一步增大，湍流狀態(tài)再次出現(xiàn)反轉，然后沿著AE 狀態(tài)發(fā)展，由于尚未完全從激波干擾中恢復，僅在橫流區(qū)子午面觀察到比較接近ISO 的狀態(tài)。

圖24 R5 站位湍流狀態(tài)法向演化Fig.24 Evolution of normal distribution of turbulent states at location R5

6 結 論

對高超聲速下有攻角錐裙的激波/邊界層干擾問題進行直接數(shù)值模擬研究，通過對比背風區(qū)、橫流區(qū)、迎風區(qū)3 個周向子午面的典型流場結構、分離泡非定常運動特性、下游再附邊界層演化過程，得到以下主要結論：

1） 激波/邊界層干擾導致壁面壓力、摩阻和熱流顯著升高，分別達到上游的10、2、5 倍左右。熱流與壓力在上游滿足QP85 ≈ 1.0 關系式，在干擾區(qū)由分離、拐角構型變化誘發(fā)2 次上升，再附后迅速降低。熱流與摩阻的RAF 關系式在分離區(qū)完全失效，在迎風區(qū)、橫流區(qū)的無干擾上游、恢復區(qū)下游均符合較好。

2） 各周向子午面內均出現(xiàn)流動分離現(xiàn)象，分離泡面積脈動主頻均為St≈ 0.009 7，采用低階POD 模態(tài)重構干擾流場，表明分離泡非定常運動的低頻特征。在三維橫流影響下，平均分離點壓力脈動在迎風區(qū)和橫流區(qū)的低頻段增強，對應時間尺度分別為上游的73.2、46.4 倍，表明激波低頻振蕩與分離泡膨脹/收縮運動密切相關，而在背風區(qū)兩者無關。

3） 再附邊界層演化研究結果表明，干擾后速度、溫度、湍動能在近壁區(qū)迅速恢復，但外層由于恢復較慢，攻角的存在導致再附點附近雷諾應力分量放大因子劇增，流向分量恢復最迅速。三維橫流效應影響主要體現(xiàn)在背風區(qū)，雷諾應力各分量峰值位置持續(xù)外移，下游近壁區(qū)的各向異性峰值減弱，湍流邊界層仍在恢復中。