超聲速膨脹角入射激波/湍流邊界層干擾直接數值模擬

童福林，孫東，袁先旭，李新亮

1. 中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，綿陽 621000 2. 中國科學院 力學研究所 高溫氣體動力學重點實驗室，北京 100190 3. 中國空氣動力研究與發展中心 計算空氣動力研究所，綿陽 621000 4. 中國科學院大學 工程科學學院，北京 100049

激波/湍流邊界層干擾區內流動參數變化劇烈，涉及了湍流非平衡效應、可壓縮效應、邊界層分離和再附、激波/激波干擾、分離激波的非定常振蕩等多種復雜流動現象，流動機理極其復雜。Zheltovodov[1]將激波-湍流邊界層干擾的特殊物理現象歸納總結為以下6大類：① 分離激波非定常運動對邊界層內湍流的放大機制；② 非定常激波對外流湍流的放大機制；③ 膨脹波對湍流的抑制作用；④ 再附區邊界層特征；⑤ Taylor-G?rtler渦的形成；⑥ 分離區內湍流的層流化。

經過半個多世紀的深入探索，國內外學者在激波與湍流邊界層的相互作用機制方面已取得長足的進展，例如快速畸變近似[2]、初始分離準則[3]、干擾模式與相似律[4]、分離激波的低頻振蕩及其形成機制[5-6]等。當前，風洞試驗以及高精度數值模擬研究主要針對壓縮拐角和入射激波/平板兩類簡單構型。在風洞試驗方面，自20世紀50年代以來，研究人員采用熱線、激光風速儀、三維粒子成像測速技術、納米粒子平面激光散射技術等先進測量手段準確獲得了干擾區內高分辨率瞬態結構及脈動信息。研究結果表明，強干擾下激波的非定常運動是造成快速畸變理論分析結果嚴重偏離試驗數據的重要因素[7]；在激波作用后的下游湍流邊界內質量通量脈動強度顯著增強，湍流剪切應力的變化要明顯強于雷諾應力其他分量[8]。近些年，Humble和Scarano[9]通過試驗獲得了激波干擾問題中的三維瞬態流場結構，發現上游擬序結構對激波低頻振蕩特性影響顯著。與此同時，在直接數值模擬(Direct Numerical Simulation，DNS)研究方面，國內外也都取得了較大的突破。美國Martin等[10-13]對馬赫數Ma=2.9、雷諾數Reθ=2 300下激波與湍流邊界層的相互作用進行了一系列直接數值模擬研究，其計算問題為24°壓縮拐角和12°入射激波/平板干擾，計算結果也得到了Bookey等[14]低雷諾數試驗結果的驗證和確認。在國內，文獻[15-18]基于其DNS數據，系統開展了拐角角度、壁溫、拐角導圓、可壓縮性以及轉捩等因素對分離泡、邊界層演化特性以及激波振蕩特性的影響規律和作用機制。

超聲速膨脹角是流體力學基礎問題之一，普遍存在于各類高速飛行器的內外表面，尤其是在超燃沖壓發動機進氣道中。工程實踐表明，氣流在進氣道唇口附近急劇壓縮后形成較強的入射激波，與前體壓縮面以及后部的膨脹區邊界層將產生相互干擾作用。情況惡劣時，將會嚴重影響進氣道工作性能，甚至導致發動機不啟動。與以往傳統的激波/湍流邊界層干擾問題不同的是，此時，膨脹角內存在較強的順壓梯度和膨脹波系，這會對激波運動及湍流脈動特性均帶來顯著影響，進而膨脹區內相互干擾也將呈現不同的演化機制。因此，深入開展膨脹角激波/湍流邊界層干擾區內復雜流動現象的機理研究，將有助于進一步加深對該問題的理解和認識，具有重要的工程優化設計和應用背景。

早期研究主要以風洞試驗為主，計算模型也多采用入射激波與平板/膨脹角構型。Chew[19]對Ma∞=1.8、2.5下的6°膨脹角入射激波/湍流邊界層干擾問題進行了試驗研究，發現當入射激波再入點位于膨脹角角點附近(約3～4倍邊界層厚度范圍內)時，膨脹效應對上游流動的影響較為明顯。Chung和Lu[20]分析研究了Ma∞=8下入射激波再入點位置對膨脹區內物面壓力平均量及脈動量的影響規律，但由于入射激波強度較弱，干擾區內并沒有發生流動分離現象。隨后，White和Ault[21]進一步研究發現在強激波干擾下膨脹效應將會導致干擾區內分離泡尺度的減小。最近，Sathianarayanan和Verma[22]試驗研究了考慮側壁效應的入射激波與膨脹角湍流邊界層相互作用，著重探究了入射激波強度、膨脹角角度以及入射點位置等因素對分離區三維形態的影響機制。相較于風洞試驗，國內外在膨脹角激波/湍流邊界層干擾的高精度數值模擬研究方面，相關工作還較為少見。Konopka等[23]采用大渦模擬對Ma∞=1.76的入射激波/膨脹角湍流邊界層干擾問題進行了初步研究，激波入射點位于膨脹角的下游，重點關注了激波干擾對膨脹區內湍流統計特性的影響規律。

本文采用直接數值模擬方法對超聲速膨脹角入射激波與湍流邊界層相互作用問題進行系統研究。著重探討強激波干擾下入射激波入射點位置的改變對干擾區內復雜流動結構的影響規律，如分離泡、物面壓力脈動特性及激波非定常運動、湍流邊界層統計特征等。采用本征正交分解方法，分析比較不同入射位置下膨脹區內動力學過程的差異。為了便于比較和驗證結果，來流參數的選取與Bookey等[14]的試驗和Priebe等[13]的DNS相近。

1 計算參數

控制方程為曲線坐標系(ξ,η,ζ)下三維可壓縮無量綱Navier-Stokes方程組，流場變量和長度變量分別采用無窮遠處來流參數和單位毫米進行無量綱化，具體形式為

?tU+?ξ(F+Fv)+?η(G+Gv)+

?ζ(H+Hv)=0

(1)

式中：U=J-1[ρ,ρu,ρv,ρw,ρE] 為守恒變量;F和Fv分別為ξ方向上的無黏和黏性通量;(G,Gv)和(H,Hv)分別對應于η和ζ方向;J為坐標系轉換時的雅可比矩陣。各項的具體表達式參見文獻[10]。

采用高精度有限差分解算器OpenCFD-SC進行DNS計算。在本文作者前期研究中[15-18,24]，采用該軟件對大量壓縮拐角和入射激波/平板湍流邊界層干擾問題進行了直接數值模擬研究，DNS結果可靠性和準確性都得到了充分的驗證與確認。DNS計算時，我們采用了Martin等[25]優化構造的WENO_SYMBO_LMT格式以及Steger-Warming流通量分裂方法計算無黏通量。該數值格式在激波/湍流邊界層干擾問題中具有較好的魯棒性和計算精度，可以保證在精確捕捉湍流邊界層內不同尺度流動結構的同時，又能較好地抑制強激波間斷處的數值振蕩。另外，采用八階中心差分格式對黏性項進行計算，時間離散采用的是三階精度的Runge-Kutta方法。

如圖1所示，計算模型為入射激波與平板/膨脹角湍流邊界層的相互作用問題。氣流方向為從左往右，坐標系原點取為膨脹角角點處，膨脹角為10°，入射激波的激波角為30°。Xin和Xup分別對應為入射激波在物面上的名義入射點和在上邊界的入射點位置。通過改變計算域上邊界中Xup的流向位置使得入射激波打在膨脹角壁面上不同的流向位置Xin。計算域流向長度為Lx=470 mm，流向跨度為-363 mm

圖1 計算域及網格示意圖Fig.1 Illustration of computational domain and grid

具體來流條件如下：來流馬赫數為Ma∞=2.9，基于單位長度的來流雷諾數為Re∞=5 581.4 mm-1，來流靜溫為T∞=108.1 K，壁面溫度取為Tw=307 K。計算的DNS工況分別為Case1～Case4，其中Case1為無膨脹情況下入射激波/平板湍流邊界層干擾，Case2～Case4為膨脹角入射激波/湍流邊界層干擾，其壁面上名義入射點位置Xin分別對應為膨脹角點上游、角點處以及角點下游3種工況。為了更好地揭示膨脹效應的影響，這里通過調整Xup將Case1和Case3的壁面入射點位置Xin均取為0，Case2和Case4的入射點Xin則分別取在角點上下游約δ的位置處 (δ對應為上游參考點ref處的湍流邊界層厚度，見圖1)。各工況中其他DNS參數設置均完全相同。

如表1所示，各DNS工況的網格點數均為3 200×200×140(流向Nx×法向Ny×展向Nz)，計算網格采用代數解析方法生成。為了精確捕捉流動信息，流向網格點在膨脹角激波與湍流邊界層干擾區內均勻密集分布，法向網格在近壁區采用了雙曲正切函數的加密處理，以保證在整個邊界層內有120個網格點，展向網格點均勻分布。以x=-60 mm處壁面量為參考(見圖1中參考點ref)，膨脹角干擾區內流向和展向網格尺度分別為Δx+=5.0和Δz+=7.1,壁面及邊界層外緣的法向網格尺度分別Δyw+=0.7和Δye+=11.2，與Priebe等[13]的DNS較為接近。表2 還分別給出了上游參考點ref處湍流邊界層的邊界層厚度δ、位移厚度δ*、動量厚度θ、形狀因子H和物面摩阻系數Cf。可以看到，本文結果與DNS數據[13]及試驗結果[14]較為接近。

表1 DNS工況參數Table 1 Parameters for DNS cases

表2 參考點湍流邊界層參數

首先，為了驗證本文DNS計算展向寬度的合理性，圖2 給出了膨脹角干擾區內速度展向關聯函數Rαα(rz)在邊界層內不同法向位置處的分布情況，這里rz為展向間距，Ruu、Rvv和Rww分別為流向、法向和展向速度關聯函數，具體表達式參見文獻[26],yn為物面法向距離。可以看到，3個方向上脈動速度的關聯函數在展向距離rz大于半個計算寬度時，均衰減到0附近，這表明本文DNS計算的展向寬度是合理的，能夠有效捕捉干擾區大尺度結構。需要特別說明的是，本文將在后續結果分析和討論中，通過與DNS結果[13]以及試驗數據[14]的比對，進一步驗證和確認本文計算結果的可靠性。

圖2 兩點展向相關函數分布Fig.2 Distribution of two-point correlation as a function of spanwise spacing

DNS計算時，在經過兩個無量綱時間(Lx/U∞)后流場達到統計定常狀態，隨后對瞬態場進行統計采樣。總無量綱采樣時間跨度約為500δ/U∞，物面壓力脈動信號的采樣間隔為0.06δ/U∞，共獲得400個三維瞬態流場樣本。如無特別說明，本文所指的平均定義為時間和展向平均。

2 流場結構

圖3給出了膨脹角干擾區瞬態密度梯度場，這里粉色曲線代表Ma∞=1的瞬態等值線，圖中Xsep為平均分離點位置。為了更好地顯示流場結構，采用文獻[10]中定義的流場變量NS：

(2)

圖3 瞬態密度梯度場Fig.3 Instantaneous gradient of density

Simpson[27]依據瞬態分離的統計概率將邊界層分離劃分為以下3類：初始分離(Incipient Detachment, ID)、間歇性瞬變分離(Intermittent Transitory Detachment, ITD)以及瞬變分離(Transitory Detachment, TD)。前兩者分別對應回流的統計概率為1%和20%，后者則為50%，對應為平均意義上的流動分離。圖4給出了膨脹角內物面上流動分離的統計概率(γ)分布情況。這里的流動分離定義為?Us/?yn<0，其中Us為沿物面的流向速度。為了便于比較，流向坐標采用分離點坐標Xsep和上游湍流邊界層厚度δ進行歸一化。整體來看，Case1～Case4工況物面流動分離統計概率的分布規律基本類似，均以雙峰結構為主，呈現V字型分布。從定量分布來看，所有工況下峰值概率均超過0.5，這說明膨脹角內流動存在平均意義上的流動分離。可以看到，膨脹角上游區域的統計概率分布較為一致，而下游區域內流動分離特征則變化顯著，特別是Case4工況，物面上存在較大范圍的初始分離。

圖4 物面上流動分離的統計概率分布Fig.4 Statistical probability distribution of flow separation at wall

圖5 物面平均摩阻系數分布Fig.5 Distribution of mean wall skin-friction coefficients

圖6 分離泡高度沿流向分布Fig.6 Streamwise variations of separation bubble height

3 物面壓力脈動

干擾區物面壓力及其脈動特征一直以來都是激波與湍流邊界層相互作用問題的研究重點及熱點。深入分析膨脹效應對物面壓力脈動的影響機制，將有助于理解膨脹角分離激波的非定常運動規律及其物理機制，特別是不同入射位置下膨脹效應對分離激波低頻振蕩運動的影響規律。本節將通過與激波/平板湍流邊界層干擾結果的比對分析，進一步深入揭示膨脹效應對物面壓力脈動強度、概率密度函數、時空關聯特性以及功率譜密度等方面的影響。

圖7分別給出了各工況干擾區內平均物面壓力Pw/P∞及脈動強度Prms/Pw的分布情況,這里下標w和∞分別代表物面和無窮遠來流，下標rms為脈動均方根。圖中符號S表示分離區起始點的位置。需要說明的是，為了便于比較，將流向坐標進行了平移及歸一化，使得Case1計算結果中壓力升高點與Priebe等[13]的DNS結果相互重合，如圖7(a)中橫坐標X*所示。可以看到，兩者在干擾區內分布規律非常吻合，這也充分證實了本文DNS數據的準確性。從Case2～Case4的分布規律來看，膨脹效應使得干擾區下游物面壓力降低，總壓差減小，但入射激波流向位置的改變對膨脹區內總壓差沒有實質影響，不同工況下壓力分布曲線在膨脹角下游吻合良好。流向位置變化對分離區內物面壓力流向分布規律的影響則更為直接。Case4時壓力的變化則相對較為緩慢，這主要是此時膨脹角角點的強膨脹波系對入射激波存在較強抑制作用，如圖3(d)所示。從圖7(b)中還可以看到，物面壓力脈動強度的極值點主要出現在分離激波點S處，這主要是由于分離激波的存在。同時，隨著入射激波位置往下游移動，脈動峰值也逐步下降，這與Chung和Lu[20]的試驗結果是一致的。

圖7 平均物面壓力及脈動強度分布Fig.7 Distribution of mean and fluctuation intensity of wall pressure

圖8分別給出了分離點及膨脹角再附區物面壓力脈動的概率密度函數(Probability Density Function, PDF)，這里采用當地脈動壓力均方根σp對脈動量進行歸一化處理。與高斯分布的比較結果表明，各工況下分離點壓力脈動的PDF分布均呈現非對稱特征，大概率事件出現在-1

圖8 物面壓力脈動概率密度函數Fig.8 Probability density function of fluctuating wall pressure

為了更好地考察膨脹效應對物面壓力脈動p′時空關聯特性的影響規律，這里采用流向時空關聯系數R(Δξ, Δζ, Δt)，具體定義為[32]

R(Δξ,Δζ,Δt)=

(3)

式中：(x0,z0)為參考點坐標；Δξ、Δζ和Δt分別為流向間距、展向間距和延遲時間。在本文中沿展向中心線從-60 mm

圖9分別給出了分離點和下游膨脹區內物面壓力的自相關系數R(0,0,Δt)。可以看到，分離點和下游膨脹區自相關系數分布的差異較為明顯。在圖9(a)中，Case1～Case4分離點自相關系數要明顯高于其上游湍流邊界層(TBL)，而且隨著延遲時間的增大，其自相關系數仍維持在較高的量值，特別是Case1。Muck等[33]研究表明，分離激波大尺度振蕩運動產生的強間歇性是該現象的主要誘導因素。從本文計算結果來看，Case2～Case4分離點的自相關系數分布仍基本符合這一規律，但是其量值急劇下降，峰值系數對應的延遲時間也相對減小。這表明膨脹效應使得分離激波的非定常運動間歇性減弱和時間尺度減小。另外，在膨脹區x/δ=4.64處，從圖9(b)中仍可以看到自相關系數隨著延遲時間呈現急劇下降的趨勢，壓力脈動時間尺度也呈現進一步減小趨勢，特別是在ΔtU∞/δ>1.0時，脈動壓力自相關性非常弱。入射點流向位置的改變對下游再附邊界層物面壓力自相關系數分布規律的影響可以忽略不計。

圖9 物面壓力脈動自相關Fig.9 Autocorrelation of wall pressure fluctuations

本文還進一步給出膨脹效應對物面壓力脈動流向時空關聯系數的影響規律。圖10首先給出了上游湍流邊界層內物面壓力脈動的流向互相關，這里參考點(x0,z0)取為(-60 mm,7 mm)。如圖10(a)所示，隨著流向間距Δξ的增大，互相關系數的極值逐漸減小，延遲時間逐漸增大。通常，脈動壓力的對流速度可由壓力測點流向距離Δξ除以互相關曲線峰值對應的延遲時間τopt求得。本文由此計算得到的對流速度Uc約為0.6U∞～0.8U∞，這與Willmarth和Wooldridge[34]的試驗結果較為接近。圖10(b)還給出了無量綱流向間距Δξuτ/(Ucδ)與互相關系數極值的關系。這里的無量綱間距表征了兩個不同時間尺度的比值:Δξ/Uc和δ/uτ，其中前者對應大尺度渦特征對流時間，而后者則是含能渦時間尺度。可以看到，本文計算結果與試驗結果[33,35-37]吻合較好。

圖10 上游湍流邊界層物面壓力脈動流向時空 關聯Fig.10 Time-space correlations of fluctuating wall pressure in upstream turbulent boundary layer

圖11和圖12還分別給出了分離區以及下游膨脹區的流向時空關聯。相較于上游湍流邊界層，分離區及膨脹區流向互相關曲線分布規律主要有以下兩個方面值得重點關注。一方面，從圖11(a)中可以清楚看到，分離點物面壓力脈動的時空關聯曲線在負延遲時間區域仍存在較強關聯性。以往研究結果表明[31]，造成這一現象的物理機制與分離激波非定常運動密切相關。可以看到，隨著激波入射點位置往下游移動，負延遲時間區域內的相關性急劇下降，這表明膨脹效應對分離激波的非定常運動特性影響顯著。可以看到，特別在Case3和Case4工況，此時關聯曲線與上游湍流邊界層的差別較小。

另一方面，無論是在分離區還是膨脹區，在正延遲時間區域內，相關曲線均存在明顯的局部峰值。在分離區內，膨脹效應對關聯系數峰值大小的影響較為顯著，而延遲時間位置的變化則相對較小。如圖11(a)和圖11(b)所示，局部峰值對應的延遲時間ΔtU∞/δ主要位于0.46～0.78的范圍內，由此估算的對流速度Uc約為(0.41～0.69)U∞，略小于上游湍流邊界層。Muck等[33]在壓縮拐角分離泡也發現類似流動現象。在圖12中還可以看到，膨脹區內關聯曲線峰值及其延遲時間的分布則相對較為集中，峰值大小出現在0.5～0.6范圍內，延遲時間ΔtU∞/δ約為0.35，由此估算得到的對流速度約為0.91U∞，這說明膨脹區內物面壓力波的傳播存在一個加速過程，分離區與膨脹區內對流速度的差異很可能是由于湍流大尺度結構在分離區與膨脹區內不同的演化機制。在分離區內，盡管物面壓力測點位于分離泡底部，但其傳播特性主要由剪切層外部大尺度結構所決定。上游邊界層內的大尺度結構在穿過分離激波后，其對流速度將相對降低；而在膨脹區內，外層的大尺度結構又經歷一個急劇的加速膨脹過程，導致物面壓力波的傳播速度將增大。

圖11 分離區內物面壓力脈動時空關聯系數Fig.11 Time-space correlation coefficients of fluctua- ting wall pressure in separation region

圖12 下游膨脹區內物面壓力脈動時空關聯系數Fig.12 Time-space correlation coefficients of fluctua- ting wall pressure in downstream expansion region

為了定量化描述膨脹效應對分離激波非定常運動特性的影響規律，圖13還給出了物面壓力脈動加權功率譜密度云圖(Weighted Power Spectral Density, WPSD)，圖中符號S和R分別代表平均分離點和再附點的流向位置，黑色實線代表膨脹角角點位置(EC)。加權功率譜密度具體定義為[31]

(4)

式中：f為頻率；Ψ(f)為功率譜密度。

如圖13(a)所示，對于入射激波/平板湍流邊界層干擾問題，上游湍流邊界層物面脈動壓力的峰值頻率約為fδ/U∞=1，而在分離點附近，由于分離激波的大尺度流向低頻振蕩，物面脈動壓力的低頻能量急劇增加，峰值頻率出現在fδ/U∞=0.01附近，隨后在下游再附邊界層內，峰值頻率又恢復到較上游略低的高頻區。本文Case1的計算結果與DNS結果[13]以及風洞試驗[14]均較為一致。從圖13(b)～圖13(d)中可以清楚看到，膨脹角分離點附近脈動壓力加權功率譜密度的低頻能量呈現急劇降低的趨勢，在fδ/U∞<0.1頻段內Case2～Case4工況均無明顯能量峰值，這表明膨脹效應極大地抑制了分離激波的大尺度低頻振蕩。在定量比較方面，圖14還分別給出各頻段能量在總脈動能量的占比(Fr)沿流向的分布情況，其中選取的低頻段(LF)為fδ/U∞<0.1，而中頻段(MF)為0.1

圖13 物面壓力脈動加權功率譜密度云圖Fig.13 Weighted power spectral density map of fluctuating wall pressure

圖14 物面壓力脈動各頻段占比分布Fig.14 Frequency-band percentage of fluctuating wall pressure

4 膨脹區湍流邊界層

本節將討論入射激波位置對膨脹區湍流邊界層統計特性的影響規律，如平均速度剖面、雷諾應力及其各向異性張量不變量和湍動能輸運機制等。結果均取自于x/δ=4.64處，上標+表示采用當地壁面量進行的無量綱化。

圖15 平均速度剖面分布Fig.15 Distribution of mean velocity profile

圖16(a)和圖16(b)分別給出了湍流邊界層雷諾正應力τ11和切應力τ12的分布情況。從量值大小來看，由于激波對湍流的增強作用，雷諾應力呈現顯著升高的趨勢。相較于上游TBL，Case1雷諾正應力和切應力分別增大了約4倍和10倍。以往研究表明[38]，強膨脹效應對湍流脈動具有抑制作用，會使得湍流邊界層出現層流化的趨勢。可以看到，Case2～Case4雷諾應力變化規律也進一步證實了該結論。此時，膨脹區內各分量量值急劇減小，與上游TBL已較為接近，但兩者在分布規律方面仍存在著較大差異。如圖所示，膨脹區內峰值雷諾應力出現在邊界層外層區域yn/δ>0.3，而上游充分發展湍流邊界層則集中在0.01

為了進一步定量評估內外層雷諾應力演化機制的差異，圖16(c)還給出了結構參數-τ12/τii在邊界層內分布情況。Klebanoff[39]和Grilli等[40]研究發現，結構參數在邊界層0.1

圖16 雷諾應力分量及結構參數Fig.16 Reynolds stress components and structure parameter

圖17給出了膨脹區內湍流邊界層雷諾應力各向異性張量不變量的分布情況,圖中橫坐標為第3不變量IIIb，縱坐標為第2不變量IIb，定義為[40]

(5)

式中：上標～代表Favre平均；K代表湍動能。

與以往研究結果[40]一致：在上游可壓縮TBL內，靠近壁面的區域，流動以兩組元湍流為主，隨后在近壁區則逐步逼近一組元湍流，而邊界層外緣流動表征為各向同性狀態。對于Case1，平板下游再附邊界層內湍流狀態在靠近壁面區域呈現兩組元軸對稱，近壁區內湍流則沿著軸對稱壓縮特征線變化，并在邊界層外緣附近趨近于軸對稱壓縮狀態，如圖17第2個圖(放大圖)所示，這與以往壓縮拐角激波/湍流邊界層干擾再附區的研究結果[40]是一致的。從Case3與Case1的結果比較來看(兩者激波入射位置完全相同)，膨脹效應使得下游再附邊界層近壁區內湍流沿著兩組元特征線逐步趨近一組元湍流狀態，而邊界層外緣呈現逼近各向同性狀態的態勢，這表明膨脹效應加速了下游再附邊界層的恢復過程。從Case2和Case4的分布規律來看，激波入射流向位置的改變對近壁區內湍流狀態的影響更為明顯，而外層的變化則相對要小得多。可以看到，當入射激波位于膨脹角下游，近壁區湍流呈現遠離一組元湍流的態勢，與入射激波位于角區上游的演化規律完全相反，這很可能是此時恢復區較短的緣故。

圖17 雷諾應力各向異性張量不變量Fig.17 Invariants of Reynolds stress anisotropy tensor

可壓縮湍動能的輸運方程為[41]

(6)

式中：C為對流項；T為湍流輸運項；P為生成項；V為黏性擴散項；П為壓力膨脹項；ε為黏性耗散項；M為可壓縮質量通量項。各項具體表達式可參見文獻[41]。以往研究表明[41]，湍動能輸運機制中起主要作用的是生成項P、輸運項T、擴散項V和耗散項ε。本節將重點研究膨脹效應對湍動能輸運方程各項的影響規律。

圖18 上游湍流邊界層湍動能輸運方程各項分布Fig.18 Budget of turbulent kinetic energy in upstream turbulent boundary layer

圖19分別給出了膨脹區內湍動能輸運方程各項的分布情況。總體來看，Case2～Case4膨脹區邊界層近壁區和外層區域的湍動能輸運機制均與上游湍流邊界層較為類似，膨脹效應的影響主要是體現在各項具體量值。為了反映真實量值的變化規律，這里各項均采用無窮遠處來流進行歸一化。值得特別關注的是圖19(a)中的生成項。如圖中黑色曲線所示，此時生成項在Case1的再附邊界層內呈現雙峰結構，外層峰值要明顯高于內層峰值，這主要是由于此時再附邊界層仍未恢復到平衡態，分離區剪切層產生的外層大尺度結構仍占主導。可以看到，膨脹效應使得內層生成項峰值升高，而外層生成項峰值降低。隨著入射激波位置往下游移動，這一演化趨勢有所緩解。從圖19(b)～圖19(d)中還可以清楚看到，膨脹效應顯著降低了近壁區黏性耗散項、湍流輸運項和湍流擴散項，特別是在內層yn/δ<0.01的范圍，但入射點流向位置改變對各項分布規律的影響可以忽略不計。

圖19 膨脹區湍動能輸運方程各項分布(x/δ=4.64)Fig.19 Budget of turbulent kinetic energy in expansion region(x/δ=4.64)

5 速度場本征正交分解

這里采用本征正交分解方法探究膨脹角入射激波與湍流邊界層干擾的非定常動力學過程。通過本征正交分解(Proper Orthogonal Decomposition, POD)方法可以對復雜高維度的流場進行低階近似，提取出非定常演化歷程中能量占優的特征模態。這里將著重分析入射激波在不同流向位置時膨脹效應對干擾區流場內非定常運動歷程的影響機制。具體的POD分析方法介紹可參見文獻[43]。

與Mustafa等[44]的研究相似，這里POD分析主要針對400個展向平均流向速度脈動場進行操作。假設非定常流向速度場Us(x,y,t)，POD分析可以確定一族正交基函數Φj(x,y)，j=1,2,…，具體分解過程為[43]

(7)

式中：〈Us(x,y,t)〉為時空平均場；aj(t)為第j個模態隨時間變化的模態系數；Nt為模態總數。

圖20給出了各模態能量分布。為了驗證POD分析結果的收斂性，圖中還給出了Case1工況Nt=600的結果。可以看到，模態總數Nt分別取400和600時，能量分布曲線基本重合。因此，后續的POD分析都基于Nt=400的樣本空間進行。從圖20(a)中整體分布趨勢來看，隨著模態階數增大，模態能量急劇降低。值得注意的是，Case1高階模態能量衰減符合j-11/9律，該研究結果與Mustafa等[44]在壓縮拐角中的發現相吻合。同時還可以看到，Case2～Case4工況在膨脹效應作用下其高階模態能量的衰減律與j-11/9律存在著較為明顯的差異。此外，隨著入射激波位置往下游移動，第1階模態單調下降，第2階模態能量變化較小，而高階模態能量(j>10)則呈現較為明顯的升高趨勢。對于第1階模態，其對總能量的貢獻最大，下文也稱為主能量模態。從圖20(b)中還可以清楚看到，相較于Case1、Case2～Case4工況模態累積能量曲線的梯度存在逐步升高的趨勢，這說明高階模態能量在總能量中的占比逐漸增強。以100～400階模態為例，Case1時其能量占比約為8%，而Case4時其能量貢獻增大到了約26%。

圖20 POD模態能量分布Fig.20 Energy distributions of POD modes

圖21 POD模態系數aj(t)頻數分布直方圖Fig.21 Histogram of frequency distribution of time- varying coefficient aj(t) of POD modes

圖22分別給出了各工況下主能量模態和第100階模態的空間結構分布。為了便于比較說明，圖中黑色曲線和粉色點畫線分別代表平均流線和聲速線。如圖所示，各工況下主能量模態空間結構特征與高階模態完全不同，前者主要集中在分離激波附近及分離區內，特別是在分離泡剪切層角部、分離泡內及其下游再附區域。而高階模態結構則表征為正負交替的小尺度結構,這與湍流邊界層內小尺度的高頻脈動結構密切相關。可以看到，隨著入射激波位置往下游移動，主能量模態及高階模態結構的分布規律基本一致，只是特征結構分布區域有所減小。

此外，為了觀察主能量模態和高階模態動力學性質的差異，圖23分別基于POD主能量模態和高階模態對膨脹角內非定常流場進行了低階重構。這里定義如下變量Ms[45]：

圖22 POD模態1(左)和100(右)空間分布Fig.22 Spatial distributions of POD modes 1 (left) and 100 (right)

(8)

圖23 基于POD主能量模態低價重構得到的 分離泡流量隨時間變化Fig.23 Time series of separation bubble flux based on low-order reconstructions using the first POD mode

圖24 基于POD高階模態低價重構得到的分離 泡流量隨時間變化Fig.24 Time series of separation bubble flux based on low-order reconstructions using high-order POD mode

6 結 論

本文采用直接數值模擬方法研究了來流馬赫數2.9、30°激波角的入射激波與膨脹角湍流邊界層相互作用問題，詳細分析了入射激波分別位于膨脹角上游、膨脹角角點和膨脹角下游3種工況下膨脹角干擾區內復雜流動現象的一些基本問題，如分離泡、物面壓力脈動及激波的非定常運動特性，膨脹區湍流邊界層統計特性和相干結構動力過程等，得到以下結論：

1) 入射激波位置的改變對膨脹角分離泡影響顯著。隨著入射激波位置往下游移動，分離泡長度及高度急劇減小，特別是入射激波位于膨脹角角點和下游區域時，分離泡長度和高度分別約為無膨脹角工況的17%和10%。

2) 隨著入射激波位置往下游移動，物面壓力脈動急劇衰減，而概率密度函數的變化則可相對忽略不計。時空關聯分析表明，分離泡內物面壓力脈動對流速度均有所降低，而下游膨脹區對流速度將相對增大。研究發現，膨脹效應極大地抑制了分離激波的低頻振蕩運動。

3) 膨脹區再附湍流邊界層仍處于非平衡態。統計結果表明，膨脹效應對平均速度剖面對數區和尾跡區影響顯著。入射激波位置的改變將使得邊界層內層結構參數升高而外層降低，同時近壁區內湍流呈現遠離一組元湍流的態勢。湍動能輸運機制與無膨脹角工況基本類似，差異主要體現在輸運方程各項的量值。

4) 速度場本征正交分解結果表明，膨脹角主能量模態的能量占比急劇降低，高階模態能量占比逐漸增強。前者空間結構主要集中在分離激波以及剪切層根部附近，對應為分離泡的低頻膨脹/收縮過程；而后者則以小尺度正負交替脈動結構為主，表征為分離泡高頻脈動。入射激波位置的改變對膨脹區內動力學過程影響主要體現在分離泡脈動強度。

致 謝

感謝國家超級計算天津中心、國家超級計算長沙中心、國家超級計算廣州中心以及中國科學院網絡中心超級計算中心提供計算機時。