童福林，周桂宇，周浩，張培紅,*，李新亮2,

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

激波與湍流邊界層的相互作用是現(xiàn)代高速飛行器氣動設計中不可或缺的基礎問題，至今仍未被充分理解。激波/湍流干擾區(qū)內會出現(xiàn)流動分離/再附、局部強壓力脈動和熱流峰值，這會對飛行器氣動性能、防熱層以及結構疲勞等方面產生顯著影響。因此，進一步深入研究干擾區(qū)內復雜流動現(xiàn)象有助于加深對該問題的理解認識，為工程應用提供重要的理論參考依據(jù)。Dolling[1]和Gaitonde[2]從熱流預測、激波的非定常運動特性以及流動控制等方面對該問題進行了詳細的綜述。

依據(jù)激波的產生方式不同，可以將激波與邊界層干擾劃分為[3]：壓縮拐角、入射激波干擾、雙錐、后掠壓縮拐角、單楔、雙楔以及內流道問題。壓縮拐角和入射激波干擾是其中兩類最具代表性的流動構型。自20世紀40年代以來，國內外大量學者對這兩類構型進行了系統(tǒng)的風洞試驗和數(shù)值模擬研究。

在風洞試驗方面，Settles和Fitzpatrick[4]研究了不同激波強度下壓縮拐角干擾區(qū)內物面壓力、摩阻及平均速度的演化規(guī)律。Ardonceau[5]和Smits等[6]分析了激波干擾對湍流脈動的增強機制。結果表明，激波干擾湍流剪切應力的影響要明顯強于雷諾應力的其他分量，干擾區(qū)下游湍流邊界層內質量通量脈動強度顯著增強。分離激波的低頻振蕩現(xiàn)象及其物理機制一直以來都是風洞試驗的研究熱點。Andreopoulos和Muck[7]發(fā)現(xiàn)上游湍流邊界層的猝發(fā)現(xiàn)象與低頻振蕩運動密切相關。但隨后，Erengil和Dolling[8]的研究結果表明，上游壓力脈動才是激波非定常運動的主要來源。目前，對于該問題的實驗研究尚未達成共識。大量研究表明，激波運動機制與上游邊界層速度型剖面[9]、大尺度高低速條帶結構[10]等因素密切相關。

在數(shù)值模擬方面，隨著計算速度和數(shù)值格式的飛速發(fā)展，直接數(shù)值模擬(Direct Numerical Simulation，DNS)方法已逐漸成為激波/湍流邊界層干擾復雜流動機理方面的重要研究手段。相較于風洞試驗，DNS可以直接給出干擾區(qū)內試驗難以測量的流場信息。Adams[11]首次采用DNS方法研究了壓縮性對干擾區(qū)下游湍流結構的影響機制，發(fā)現(xiàn)邊界層猝發(fā)頻率與激波運動頻率較為接近。Ringuette等[12]數(shù)值研究了雷諾數(shù)對干擾區(qū)物面壓力脈動頻譜特性、分離區(qū)長度以及湍流脈動的影響規(guī)律。Priebe等[13-14]著重探討了激波運動的低頻振蕩現(xiàn)象。低通濾波后的DNS瞬時流場表明，分離激波的非定常運動與分離泡的膨脹/收縮存在較強關聯(lián)。此外，李新亮等[15]研究了激波低頻振蕩機制以及干擾區(qū)湍動能的輸運機制。結果證實了激波的低頻振蕩與上游邊界層擬序結構無關。近年來，童福林等[16-18]開展了大量的激波/湍流邊界層干擾直接數(shù)值模擬研究，探討了激波強度、壁面溫度、馬赫數(shù)等因素對干擾區(qū)內復雜流動現(xiàn)象的影響規(guī)律。

總體來看，國內外在激波與湍流邊界層相互作用問題上取得了長足的進步，并在一些復雜問題流動機理方面達成共識。但在干擾區(qū)內物面剪切應力統(tǒng)計特性方面，相關試驗及DNS的研究報道較為少見。Murthy和Rose[19]對馬赫數(shù)Ma=2.9下的入射激波湍流邊界層干擾問題進行了實驗研究，獲得了干擾區(qū)內物面剪切的平均量和脈動量。進一步深入開展干擾區(qū)內物面剪切應力統(tǒng)計特性的演化規(guī)律研究，有助于為改進現(xiàn)有湍流模型和亞格子模型提供理論支撐。

本文采用直接數(shù)值模擬方法對入射激波/平板湍流邊界層相互作用問題進行數(shù)值研究。著重探討分離激波低頻振蕩運動對物面剪切應力功率譜密度的影響機制，研究分離泡內流向及展向剪切應力分量的概率密度分布規(guī)律。采用本征正交分解方法，分析比較了流向剪切應力脈動與物面壓力脈動的差異。為了便于比較和驗證結果，計算參數(shù)的選取與Bookey等[20]的實驗結果和Priebe等[13]的DNS結果相近。

1 計算參數(shù)

直接數(shù)值模擬的控制方程為三維可壓縮無量綱Navier-Stokes方程組：

(1)

式中：Q為守恒變量；F、G和H為3個方向上的無黏通量；Fv、Gv和Hv為3個方向對應的黏性通量，具體表達式參見文獻[21]。方程的無量綱化采用無窮遠處來流參數(shù)以及單位特征長度。計算時，為了抑制激波間斷區(qū)的數(shù)值振蕩同時保證對湍流邊界層內不同尺度流動結構的高分辨率，采用Martin等[22]優(yōu)化構造的WENO_SYMBO_LMT格式以及Steger-Warming流通量分裂方法計算無黏項。同時，采用八階中心差分格式對黏性項進行離散，時間推進采用三階Runge-Kutta方法計算。需要特別指出的是，本文DNS采用的高精度差分求解器OpenCFD-SC軟件已在多個激波/湍流邊界層干擾問題[15-18]中得到了成功的驗證和確認，可以保證DNS結果的準確和可靠。

圖1 計算模型示意圖Fig.1 Illustration of computation model

如圖1所示，計算模型為入射激波與平板湍流邊界層的相互作用問題。氣流方向為從左往右，來流馬赫數(shù)為2.9，基于單位長度的來流雷諾數(shù)為5 581.4 mm-1，來流靜溫為108.1 K，壁面溫度為307 K。計算域流向跨度為-363 mm

網格點數(shù)為3 200×200×140(流向×法向×展向)，計算網格采用代數(shù)解析方法生成，流向網格在激波與湍流邊界層的干擾區(qū)內均勻分布(如圖1所示)，法向網格在近壁區(qū)采用了雙曲正切函數(shù)的加密處理，展向網格均勻分布。這里以x=-60 mm處壁面量為參考量，干擾區(qū)內網格尺度分別為Δx+=4.5、 Δy+=0.5、 Δz+=5.0,與Priebe等[13]的DNS結果較為接近。如無特別說明，本文中上游湍流邊界層的統(tǒng)計變量均取自x=-60 mm處(位于充分發(fā)展湍流邊界層內)。表1分別給出了上游湍流邊界層的馬赫數(shù)Ma、邊界層厚度δ、位移厚度δ*、動量厚度θ和物面摩阻系數(shù)Cf。

表1 上游湍流邊界層參數(shù)Table 1 Parameters of incoming turbulent boundary layer

2 結果驗證

本節(jié)通過與以往數(shù)值模擬結果和風洞試驗數(shù)據(jù)的對比分析，進一步驗證計算結果的準確性，其中包括上游湍流邊界層的統(tǒng)計特性、干擾區(qū)內平均壓力和摩阻分布以及物面壓力脈動的功率譜等。

圖2 湍流邊界層統(tǒng)計特性Fig.2 Statistical characteristics of turbulent boundary layer

圖3分別給出了激波/湍流邊界層干擾區(qū)內物面壓力pw/p∞和摩阻系數(shù)沿流向的分布情況。為了便于比較說明，這里將壓力和摩阻分布的流向坐標均進行了平移和無量綱處理，其中x*為流向坐標平移后的值，xsep為平均分離點流向坐標。計算得到的壓力和摩阻分布與Priebe等[13]的數(shù)值結果均基本重合。與Bookey等[20]的試驗數(shù)據(jù)比較來看，干擾區(qū)下游的物面壓力值要明顯高于試驗值。在直接數(shù)值模擬時，展向取為周期性邊界條件，而實際風洞試驗時展向為真實固壁，因此洞壁干擾[13]是造成該差異的主要因素之一。

另一方面，圖4給出了計算得到的上游湍流邊界層和分離平均起始點物面壓力脈動的預乘譜(Pre-Multiplied Power Spectral Density,fPSD)，其中f為頻率，PSD為功率譜密度?？梢?，在無干擾區(qū)內，壓力脈動的無量綱峰值頻率出現(xiàn)在1.0U∞/δ附近,U∞為來流速度。由于干擾區(qū)內分離激波的低頻振蕩運動，分離點物面壓力脈動的低頻能量在(0.004～0.01)U∞/δ的范圍內急劇增強。計算值與以往激波湍流邊界層干擾直接數(shù)值模擬得到的低頻峰值頻率范圍(0.002～0.006)U∞/δ[13]較為接近，這也證實了本文DNS計算采用的數(shù)值方法和網格分辨率能夠準確捕捉到干擾區(qū)內的分離激波低頻振蕩現(xiàn)象。

圖3 物面平均壓力及摩阻系數(shù)分布Fig.3 Distribution of wall average pressure and skin friction coefficient

圖4 物面壓力脈動預乘譜Fig.4 Pre-multiplied power spectral density of wall pressure fluctuations

3 流場結構

圖5分別給出了入射激波與平板湍流邊界層干擾區(qū)內的無量綱瞬態(tài)密度梯度場和時間平均密度場。圖中紅色和藍色曲線為Ma=1和u=0的瞬態(tài)等值線，這里u為流向無量綱速度。可以看到，在入射激波和分離激波的相互作用下，干擾區(qū)內存在著強逆壓梯度，邊界層內出現(xiàn)了大范圍的流動分離，同時在下游邊界層外緣還存在較強的壓縮波系。此外，從兩者的定性比較來看，分離泡內存在強烈的間歇性和非定常特征。

研究結果表明，物面剪切流向分量要比其展向分量大了約一個量級，這表明干擾區(qū)內物面剪切以流向剪切為主。圖6分別給出了物面流向剪切應力τx的瞬態(tài)和時均分布云圖。為了便于下文比較說明，沿流向選取了5個典型特征位置，其中E1～E3分別位于上游無干擾邊界層內、分離泡和下游再附區(qū)內，S和R分別位于平均分離和再附點。從整體分布規(guī)律來看，上游無干擾邊界層內物面剪切以條帶結構特征為主，這與邊界層近壁區(qū)的高低速條帶緊密相關。隨后，該條帶結構在分離區(qū)內被破壞并消失，分離區(qū)內流向剪切表征為強間歇特性。在干擾區(qū)下游，流向剪切的量值呈逐漸增大趨勢，但其展向分布規(guī)律與上游存在明顯差異。從時均結果也可以看到，干擾區(qū)下游的流向摩阻沿展向表現(xiàn)為強烈的非均勻性。

圖5 瞬態(tài)密度梯度和時均密度流場Fig.5 Instantaneous density gradient and mean density flow fields

圖6 物面流向剪切應力瞬態(tài)和時均分布云圖Fig.6 Instantaneous and mean contours of streamwise component of wall shear stress

4 物面剪切應力統(tǒng)計特性

4.1 預乘譜分析

分析物面剪切應力信號的功率譜密度有助于理解干擾區(qū)內各剪切分量脈動特征的演化規(guī)律。此外，分離激波的低頻振蕩運動對物面剪切脈動的影響同樣值得關注。

圖8給出了干擾區(qū)內各特征位置物面剪切應力預乘譜分布的比較情況。可以看到，對于流向分量，相較于上游E1處，分離區(qū)內和再附區(qū)下游邊界層內峰值頻率略有升高，低頻區(qū)脈動能量也有一定的增強，但從整體分布規(guī)律來看，脈動能量仍然以高頻特征為主。干擾區(qū)內展向分量的演化規(guī)律與流向分量基本類似，只是在干擾區(qū)下游S和E3處的峰值頻率略有降低，但高頻脈動仍占主導。研究結果也進一步表明，分離激波的低頻振蕩運動不會對干擾區(qū)內物面剪切應力各分量的脈動產生實質影響。

圖7 湍流邊界層內物面剪切應力信號預乘譜Fig.7 Pre-multiplied power spectral density of wall shear stress in incoming turbulent boundary layer

圖8 干擾區(qū)內物面剪切應力信號預乘譜Fig.8 Pre-multiplied power spectral density of wall shear stress in interaction region

如圖4所示，分離激波的大尺度低頻振蕩運動對物面壓力脈動影響顯著，尤其是對平均分離點附近的低頻能量。而從圖8的結果來看，激波的非定常運動對流向和展向分量的影響則要弱得多。盡管在低頻部分，有一個數(shù)量級的升高，但脈動仍以高頻能量為主。從定性分析來看，認為造成這種差異的原因很可能有以下兩方面。首先，分離激波非定常運動的影響作用主要體現(xiàn)在流場中零階分量，如壓力、質量通量[28]等，而物面剪切應力是流向或展向速度的法向梯度，為流場參數(shù)的一階分量?？梢杂^察到，盡管干擾區(qū)內剪切脈動的低頻能量有一定的升高，但其增長速率要明顯低于壓力脈動。另一方面，從圖5中還可以看到，分離激波在邊界層外層逐步弱化為弱壓縮波系，因而其非定常運動對近壁區(qū)內流動參數(shù)的影響程度會急劇減弱。從圖4和圖8的比較來看，邊界層內激波強度的弱化對物面剪切的影響程度要遠強于物面壓力脈動，具體更為詳細的定量作用機制有待下一步深入研究。

4.2 概率密度分布

為了研究干擾區(qū)內物面剪切應力演化的統(tǒng)計特性，圖9給出了E1～E3處流向剪切分量及其脈動的概率密度函數(shù)(PDF)。在本節(jié)中剪切應力分量采用當?shù)氐臅r空平均值τx,av和τz,av進行歸一化處理，分別表示為

(2)

從圖9可以看到，在上游干擾區(qū)E1和下游再附區(qū)E3，流向分量的PDF曲線都近似地呈現(xiàn)對數(shù)正態(tài)分布。盡管流動在這兩個特征位置均以附著流為主，但仍存在一定的概率出現(xiàn)負的流向剪切。以往零壓力梯度平板的試驗結果表明[29-30]，邊界層內存在較小概率的回流現(xiàn)象(Backflow)，上述研究成果也進一步證實了該流動現(xiàn)象的存在。另外，在分離泡內E2處，從統(tǒng)計意義來看，流動以負剪切為主要特征，但其函數(shù)分布在正剪切范圍內仍存在較高的可能性，這主要由于分離泡的強間歇性的緣故。對于流向剪切的脈動量，由于采用當?shù)鼐礁M行了歸一化處理，可以清楚看到，干擾區(qū)上游和下游的函數(shù)分布近似重合，且與Carlos等[26]的不可壓平板數(shù)值結果吻合較好，兩者只是在脈動量變化劇烈的區(qū)域存在較明顯的差別，這表明在這個區(qū)域內，局部應力脈動的統(tǒng)計特征是相似的。此外，在分離泡內的脈動量變化范圍要明顯大于前兩個區(qū)域。

圖10給出了干擾區(qū)內展向剪切及其脈動值的概率密度分布曲線。函數(shù)與流向分量的分布規(guī)律則完全不同。如圖所示，干擾區(qū)內的展向剪切均近似呈現(xiàn)正態(tài)分布。偏斜因子的計算表明，3個特征位置處的偏斜因子絕對值小于0.1，這說明分布函數(shù)以對稱特征為主。但在變化較劇烈的區(qū)域，分離泡內的發(fā)生概率要明顯高于其他兩個區(qū)域。與此同時，干擾區(qū)內歸一化后的展向脈動值均與Carlos等[26]的統(tǒng)計結果吻合。綜上所述，在激波湍流邊界層干擾區(qū)內，研究表明，相較于上游充分發(fā)展湍流邊界層，分離泡內的流向剪切統(tǒng)計特性變化劇烈，而展向剪切統(tǒng)計特性的變化則可忽略不計。

圖9 干擾區(qū)內流向剪切應力概率密度函數(shù)Fig.9 Probability density functions of streamwise component of wall shear stress in interaction region

為了更好地考察流向剪切與展向剪切之間的相互關系，定義兩者的夾角ψτ(t)為

ψτ(t)=arctan(τz(t)/τx(t))

(3)

圖10 干擾區(qū)內展向剪切應力概率密度函數(shù)Fig.10 Probability density functions of spanwise component of wall shear stress in interaction region

圖11給出了干擾區(qū)內夾角的概率密度函數(shù)分布。總體來看，函數(shù)呈對稱分布，在變化劇烈的區(qū)域要明顯高于高斯正態(tài)分布。對于上游湍流邊界層，夾角主要集中出現(xiàn)在[-45°,45°]范圍內，這與Jeon等[31]的研究結果是一致的。在分離泡內E2處，可以看到，小夾角事件概率略有降低。值得注意的是，E2處夾角大于25°的可能性則要明顯高于其他兩個區(qū)域，這主要是由于分離泡內流向剪切急劇降低，而展向剪切變化較小，兩者較為接近，導致大夾角事件的發(fā)生概率也相應升高。在干擾區(qū)下游再附區(qū)E3處，流向剪切逐漸恢復到初始值，因而大角度事件的發(fā)生概率又隨之降低。

圖11 物面剪切應力夾角概率密度函數(shù)Fig.11 Probability density functions of angle between wall shear stress components

圖12 物面剪切應力矢量夾角和幅值聯(lián)合概率密度函數(shù)Fig.12 Joint probability density functions of angle and magnitude of wall shear stress vector

圖13 聯(lián)合概率密度隨物面應力矢量幅值的變化Fig.13 Variation of joint PDFs with magnitude of wall shear stress vector

4.3 本征正交分解

為了進一步深入分析物面剪切脈動的能量結構，采用本征正交分解(POD)方法探究了非定常物面剪切脈動場的典型相干結構。通過POD方法可以對復雜高維度的流場進行低階近似，提取出非定常演化歷程中能量占優(yōu)的特征模態(tài)。假設非定常物面剪切場為T(x,z,t)，POD分析可以確定一族正交基函數(shù)φj(x,z)，j=1,2,…，具體分解過程如下[32]：

T(x,z,t)=T(x,z,t)

(4)

POD分析針對400個流向/展向平面內瞬態(tài)物面流向剪切脈動場進行操作。依據(jù)模態(tài)特征值對模態(tài)能量Ej進行排序，歸一化的模態(tài)能量定義為

(5)

圖14給出了物面流向剪切場POD模態(tài)能量的分布情況。可以看到，高能量模態(tài)主要集中在前20個模態(tài)，約占總能量的50%。

圖15分別給出了物面流向剪切脈動場的能量占優(yōu)主模態(tài)的空間結構。圖中云圖為歸一化后的正交基向量，虛線表征了分離激波非定常運動的流向范圍，通過計算壓力脈動信號的間歇因子[16]得到。為了便于比較差異，圖中還給出了物面壓力脈動場對應的能量模態(tài)結構，這里SS1、SS2、SS3和SS40分別代表流向剪切脈動場的第1、2、3和40個POD模態(tài)，而P1、P2、P3和P40分別為壓力脈動場對應的POD模態(tài)。如圖15(a) 和圖15(b)所示，剪切場的第1個模態(tài)空間結構與壓力場較為相似，沿展向近似呈現(xiàn)二維分布，約占總能量的16%。從分布規(guī)律來看，能量結構主要集中在平均分離點S及間歇區(qū)內，這也說明了分離激波沿流向的非定常低頻振蕩運動主導了物面壓力和流向剪切脈動場?？梢郧宄吹?，兩者不同之處在于，流向剪切場在平均再附點下游還有一定的能量結構。

圖14 歸一化POD模態(tài)能量分布及累積能量Fig.14 Distribution of normalized energy of POD modes and cumulative energy

圖15 物面壓力和流向剪切應力POD模態(tài)空間分布Fig.15 Spatial distribution of POD modes of wall pressure and streamwise shear stress

如圖15(d)和圖15(f)所示，物面剪切場第2和第3模態(tài)的空間結構則與第1模態(tài)完全不同，此時，空間結構主要以再附點下游沿展向正負交替大尺度結構為主，平均分離點附近結構強度則要弱得多。以往的動態(tài)模態(tài)分解(DMD)研究表明[33]，再附點下游的摩阻分布與G?rtler-like渦結構密切相關。本文的POD結果也進一步證實了該研究結果。此外，盡管壓力場的第2和第3模態(tài)結構也集中在再附點附近，但其空間分布沿展向仍呈二維結構，在流向表征為正負交替結構，這很可能與下游分離泡沿流向的膨脹/收縮運動有關。圖15(g)和圖15(h)給出了第40個模態(tài)的空間分布規(guī)律，該模態(tài)為低能量模態(tài)，僅占總量的0.5%?？梢?，剪切脈動場與壓力脈動場的空間分布差別較小，均以小尺度結構特征為主。

5 結 論

本文采用直接數(shù)值模擬方法研究了來流馬赫數(shù)2.9、12°激波角的入射激波與平板湍流邊界層相互作用問題，詳細地分析了干擾區(qū)內物面剪切應力場的典型統(tǒng)計特征，如預乘譜、概率密度分布和相干結構等，得到以下結論：

1) 數(shù)值模擬準確捕捉到了分離激波的非定常運動。與物面壓力脈動不同的是，分離激波低頻振蕩運動對物面剪切應力預乘譜沒有實質影響，干擾區(qū)內流向及展向切應力的脈動能量仍以高頻特征為主。

2) 干擾區(qū)內分離泡對物面剪切應力各分量統(tǒng)計特性的影響機制差異明顯。流向剪切概率密度函數(shù)變化劇烈，分離泡內不再滿足對數(shù)正態(tài)分布規(guī)律，而展向剪切概率密度函數(shù)的變化則較小，近似于正態(tài)分布。

3) 主能量模態(tài)的空間結構表明，分離點附近物面剪切脈動與分離激波的低頻振蕩運動密切相關，而下游再附區(qū)內則由大尺度G?rtler-like流向渦結構占主導。

致 謝

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