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

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

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

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

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

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

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

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

1 計算參數(shù)

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

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

?ζ(H+Hv)=0

(1)

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

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

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

圖1 計算域及網(wǎng)格示意圖Fig.1 Illustration of computational domain and grid

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

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

表1 DNS工況參數(shù)Table 1 Parameters for DNS cases

表2 參考點湍流邊界層參數(shù)

首先，為了驗證本文DNS計算展向?qū)挾鹊暮侠硇?，圖2 給出了膨脹角干擾區(qū)內(nèi)速度展向關(guān)聯(lián)函數(shù)Rαα(rz)在邊界層內(nèi)不同法向位置處的分布情況，這里rz為展向間距，Ruu、Rvv和Rww分別為流向、法向和展向速度關(guān)聯(lián)函數(shù)，具體表達式參見文獻[26],yn為物面法向距離?？梢钥吹剑?個方向上脈動速度的關(guān)聯(lián)函數(shù)在展向距離rz大于半個計算寬度時，均衰減到0附近，這表明本文DNS計算的展向?qū)挾仁呛侠淼?，能夠有效捕捉干擾區(qū)大尺度結(jié)構(gòu)。需要特別說明的是，本文將在后續(xù)結(jié)果分析和討論中，通過與DNS結(jié)果[13]以及試驗數(shù)據(jù)[14]的比對，進一步驗證和確認本文計算結(jié)果的可靠性。

圖2 兩點展向相關(guān)函數(shù)分布Fig.2 Distribution of two-point correlation as a function of spanwise spacing

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

2 流場結(jié)構(gòu)

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

(2)

圖3 瞬態(tài)密度梯度場Fig.3 Instantaneous gradient of density

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

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

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

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

3 物面壓力脈動

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

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

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

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

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

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

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

(3)

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

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

圖9 物面壓力脈動自相關(guān)Fig.9 Autocorrelation of wall pressure fluctuations

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

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

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

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

圖11 分離區(qū)內(nèi)物面壓力脈動時空關(guān)聯(lián)系數(shù)Fig.11 Time-space correlation coefficients of fluctua- ting wall pressure in separation region

圖12 下游膨脹區(qū)內(nèi)物面壓力脈動時空關(guān)聯(lián)系數(shù)Fig.12 Time-space correlation coefficients of fluctua- ting wall pressure in downstream expansion region

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

(4)

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

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

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

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

4 膨脹區(qū)湍流邊界層

本節(jié)將討論入射激波位置對膨脹區(qū)湍流邊界層統(tǒng)計特性的影響規(guī)律，如平均速度剖面、雷諾應(yīng)力及其各向異性張量不變量和湍動能輸運機制等。結(jié)果均取自于x/δ=4.64處，上標+表示采用當(dāng)?shù)乇诿媪窟M行的無量綱化。

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

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

為了進一步定量評估內(nèi)外層雷諾應(yīng)力演化機制的差異，圖16(c)還給出了結(jié)構(gòu)參數(shù)-τ12/τii在邊界層內(nèi)分布情況。Klebanoff[39]和Grilli等[40]研究發(fā)現(xiàn)，結(jié)構(gòu)參數(shù)在邊界層0.1

圖16 雷諾應(yīng)力分量及結(jié)構(gòu)參數(shù)Fig.16 Reynolds stress components and structure parameter

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

(5)

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

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

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

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

(6)

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

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

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

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

5 速度場本征正交分解

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

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

(7)

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

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

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

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

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

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

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

(8)

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

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

6 結(jié) 論

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

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

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

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

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

致 謝

感謝國家超級計算天津中心、國家超級計算長沙中心、國家超級計算廣州中心以及中國科學(xué)院網(wǎng)絡(luò)中心超級計算中心提供計算機時。