激波/湍流邊界層干擾分離泡直接數值模擬

童福林，董思衛，段俊亦，李新亮

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

數十年來，激波與湍流邊界層的相互作用一直是工程界和學術界都十分關注的熱點問題和基礎性問題，至今仍未得到徹底解決。在強激波的誘導作用下，干擾區內會出現大范圍的流動分離現象，并形成復雜的波系結構。大量研究表明，對于無分離或小分離邊界層流動，現有湍流模型往往有較好的預測精度，數值結果與試驗測量數據吻合較好。而對于存在大范圍分離流的激波/湍流邊界層干擾問題，工程上常用的湍流模型往往不能很好地重現干擾區內某些重要的流動現象，例如多數模型能夠較為準確的預測出再壓縮過程后的峰值熱流和壓力，但采用不同湍流模型往往給出了不同的分離區長度和分離過程中的壓升結果。因此，進一步深入研究干擾區內分離泡的典型特征非常有助于加深對該問題的理解認識，同時也可為湍流模型的改進提供重要的理論依據。

國內外學者針對壓縮拐角和平板入射激波這兩類經典構型開展了大量的風洞試驗和高精度數值模擬研究，特別是在初始分離準則和分離區長度的準確預測方面。Roshko等通過試驗測量數據發現，大尺寸分離區出現的臨界角隨雷諾數增長具有增大的趨勢，而Settles等的試驗結果進一步證實了，對于完全湍流情況，該臨界角與雷諾數無關。Zheltovodov等在總結歸納大量試驗數據和數值模擬結果的基礎上，給出了絕熱壁條件下考慮了雷諾數、馬赫數、拐角角度和邊界層厚度影響的分離區長度預測公式。隨后，Jaunet等通過試驗還發現壁面加熱對分離判據的影響較小。Zhu等通過直接數值模擬(Direct Numerical Simulation，DNS)研究了壁溫對分離泡大小的影響規律，發現隨著壁溫增加，分離泡變大。他們認為壁溫的影響主要體現在對近壁雷諾數的改變上，并由此建立了分離泡大小與壁溫的0.85 次方成正比的半理論模型。童福林等對入射激波與超聲速湍流邊界層干擾進行了一系列的DNS研究，分析了強膨脹作用下分離區長度的變化規律，給出了考慮膨脹效應的分離區長度歸一化公式。

近些年，分離泡空間結構特征及其非定常運動特性越來越受到重視。Grilli等采用動態模態分解方法對其大渦模擬數據進行模態分析，利用4個低頻模態成功重構了壓縮拐角內分離泡的膨脹和收縮。Priebe等數值研究了24°壓縮拐角分離激波低頻振蕩現象與分離泡的關聯性。低通濾波結果表明，低頻振蕩運動與下游分離泡的舒張和收縮運動密切相關。然而，以往研究大多是針對展向平均后的二維流場開展的，并沒有考慮分離泡展向三維結構的影響。實際上，現有風洞試驗和數值模擬結果都證實，即便是在壓縮拐角和平板入射激波這類準二維構型下，分離區內流動仍存在顯著的三維特征，如分離流線發散和匯聚、回流斑塊和類“Owl eye”結構等。

本文采用DNS方法對入射激波/平板湍流邊界層干擾區內分離泡進行數值研究。通過定量比較3個不同展向站位分離泡的典型特征，細致探討展向三維結構對分離泡非定常運動特性的影響規律，條件統計分析分離泡內瞬態分離微團的幾何形態。此外，采用本征正交分解(Proper Orthogonal Decomposition,POD)方法，比較分析各展向站位脈動速度場的差異，同時對分離泡的非定常運動過程進行低維重構。為了便于比較和驗證結果，計算參數的選取與Dupont等的試驗和Fang等的DNS相近。

1 計算設置

控制方程為三維曲線坐標系(,,)守恒型可壓縮Navier-Stokes方程組:

(1)

式中：守恒變量，無黏通量、和，黏性通量、和的具體表達式見文獻[18]。采用自由來流參數以及單位特征長度對方程進行無量綱化，黏性系數的計算采用Sutherland公式。DNS計算時，求解器為課題組自主開發的高精有限差分軟件OpenCFD-SC，該程序已在壓縮拐角、平板入射激波、超聲速膨脹角等復雜流動中得到廣泛應用，可以保證計算結果的準確性和可靠性。黏性項的計算采用WENO_SYMBO_LMT格式以及Steger-Warming流通量分裂方法，無黏項的計算采用八階中心差分格式，時間推進采用三階Runge-Kutta方法，詳細的數值方法介紹見文獻[18, 22]。

圖1給出計算模型示意圖。模型流向長度=137.6 mm，法向高度為=12.7 mm，展向寬度為=4.4 mm，來流方向為從左往右。如圖1所示，計算域入口處層流邊界層在壁面吹吸擾動(Blowing and suction)作用下轉捩生成非定常湍流邊界層，隨后在計算域上邊界數值生成一道入射斜激波，進而與湍流邊界層產生相互作用。這里入射斜激波波角取為33.2，上邊界層激波入射點取為=82.2 mm，為入射激波在壁面上的名義入射點。來流馬赫數為=2.25，來流靜溫為=169.44 K，基于單位米的來流雷諾數為=2.5×10。圖1中:為上游湍流邊界層參考點;為層流邊界層厚度。

圖1 計算模型示意圖Fig.1 Computational model

圖2 計算網格示意圖Fig.2 Computational grid

具體邊界條件如下：首先，在計算域入口取為相同來流條件下的層流解，層流邊界層厚度為=0.6 mm。計算域出口采用超聲速出口條件及緩沖區。計算域上邊界采用簡單無反射邊界條件，同時在激波入射點前后分別設置為自由來流參數和按照Rankine-Hugoniot關系式給出波后參數。計算域下邊界采用無滑移等溫壁條件，壁溫為=321.9 K，吹吸擾動帶的起始位置分別=7.62 mm和=20.32 mm，多頻擾動波的分布函數和擾動形式與Fang等的DNS完全一致。由于本文轉捩區相對較短，為了匹配上游參考點處湍流邊界層參數(見表1中邊界層厚度雷諾數、位移厚度雷諾數*和動量厚度雷諾數)，擾動幅值和擾動基頻取為=0.2和=0.628/。在本文分析討論中，為來流速度，為參考點處的邊界層厚度;、、分別對應流向、法向、展向上的速度。

表1 參考點xref處湍流邊界層參數

2 結果驗證

在流場達到統計定常后(約兩個無量綱時間/)，對三維瞬態流場開始進行統計取樣，共獲得600個流場樣本，統計時間為/=750(為采用時間)。

圖3給出了上游層流邊界層的轉捩過程，這里采用物面法向距離渲染的瞬態密度梯度等值面云圖顯示。在壁面吹吸擾動作用下，層流邊界層間歇性顯著增強，同時邊界層外層出現了大尺度湍流凸塊(Bulges)。圖4給出了參考點處的平均速度剖面和雷諾應力分布情況。本文中平均指的是時間和展向平均。可以看到，計算結果與Bookey等的試驗數據以及Fang等的DNS結果吻合良好。從圖4(b)中還可以看到，雷諾正應力峰值出現在近壁區=13處，這與Pirozzoli等的研究結果也較為一致。

圖5還給出了干擾區內無量綱化平均物面壓力=(-)(-)的分布情況，其中和分別為物面壓力和來流壓力;為入射激波波后壓力值。與Fang等相同，圖中流向坐標=(-)，其中為=(+)/2對應的流向位置。在強激波的干擾作用下，物面壓力急劇升高，隨后在干擾區下游逐步逼近于無黏解。計算得到的壓力分布與Dupont等的試驗結果和Fang等的DNS結果均符合較好，進一步證實了本文計算結果的可靠性。

圖3 上游邊界層轉捩過程瞬態密度云圖Fig.3 Instantaneous density gradient in upstream transitional boundary layer

圖4 參考點xref處湍流統計特征Fig.4 Turbulence statistics at reference station xref

圖5 物面壓力分布Fig.5 Wall pressure distribution

3 平均結構

圖6為干擾區內時間平均分離泡云圖，這里采用速度等值面=0進行顯示。圖中:和代表平均分離點，和代表平均再附點。圖7 為物面平均摩阻分布及不同展向站位截面內分離泡高度的比較情況。這里沿展向分別取了3個站位：=0.7，=1.1和=1.75，其中站位對應為展向中截面。各站位截面內分離泡高度通過速度等值線=0的法向位置來確定。

圖6 時間平均分離泡高度云圖Fig.6 Contour of mean separation bubble height

圖7 平均摩阻及不同展向截面內分離泡高度分布Fig.7 Distribution of mean skin-friction and separation bubble height in xOy plane at different spanwise locations

從總體構型來看，分離泡呈現復雜的三維結構特征，其流向尺度明顯大于其他兩個方向，特別是在法向，其峰值高度僅約為0.12,說明分離泡整體形態以扁平型為主。同時，以摩阻曲線的過零點來確定分離區流向范圍，可以看到，分離泡在流向上可以劃分為兩個區域：次分離泡(Second Bubble, SB)，位于-區間；主分離泡(First Bubble, FB)，位于-區間。可以看到，盡管次分離泡的展向范圍略大于主分離泡，但其流向長度和法向高度則明顯小于主分離泡，以站位為例，前者分別約為后者的42%和3%。從圖7中還可以看到，主分離泡法向高度在展向也存在較大的差異，站位和的峰值高度約為站位的25%和19%，這說明主分離泡沿展向呈現中間高兩邊低的山峰型結構。

4 非定常特性

圖8給出了站位截面內不同時刻的瞬態流向速度云圖，圖中黑色箭頭為速度矢量方向，白色實線為時間平均分離泡的高度。可見，在湍流的強間歇作用下，瞬態分離泡分布極為不規則，在下游再附區內仍存在一定的出現概率。值得特別關注的是，不同時刻下分離泡面積變化非常劇烈。以流向速度<0區域的面積和為表征(圖中深藍色區域)，在/=100時，瞬態分離泡遠大于平均分離泡，其峰值高度約為0.3；而在/=300時，其面積又急劇減小，此時分離泡以小尺度多塊結構隨機分布為主。這表明分離泡具有較強的非定常特性，對應為分離泡的膨脹和收縮過程，這與Priebe等在壓縮拐角流動中的研究結論較為類似。

為了定量描述分離泡的膨脹/收縮過程，這里對3個展向站位的瞬態分離泡進行了高頻采樣，總采樣時間跨度為/=442, 約為6個無量綱時間()。圖9分別給出了各站位分離泡面積的時間序列，其中黑色虛線表示時間平均值。采樣時，還進一步將瞬態分離泡按照圖6中的定義區分為主分離泡和次分離泡分別進行統計。

圖8 不同時刻展向中截面瞬態流向速度云圖Fig.8 Contours of instantaneous streamwise velocity in spanwise mid-plane

圖9 分離泡面積時間序列Fig.9 Time series of separation bubble areas

圖中藍色曲線對應為次分離泡，紅色曲線為主分離泡，黑色曲線則代表未區分時得到的結果。

從圖9中可以清楚看到，相較于主分離泡，次分離泡的瞬時面積非常小，其整體貢獻可以忽略不計。這主要是由于盡管次分離泡在流向上具有一定的范圍，但是其法向高度非常低，因而導致其占比非常小。而從主分離泡來看，各展向站位處主分離泡的瞬時面積差異較大，但各時間序列與未區分時的統計結果均吻合較好，這也證實了主分離泡在干擾區流動分離現象中占主導作用。因此，在本文后續的分析中主要以主分離泡為研究對象。

圖10給出了各站位分離泡面積預乘功率譜的比較情況。圖中:PSD為功率譜密度，為頻率，為Strouhal 數，定義為=/。以往研究表明，干擾區分離激波的低頻振蕩運動時間尺度約為10/～100/。從本文計算結果來看，各站位處分離泡的特征頻率也基本符合這一規律。可以看到，站位和的預乘功率譜較為接近，主要以雙峰結構為主，峰值頻率分別約為=0.02和=0.04，而站位處則以單峰特征為主，峰值頻率出現在=0.04處，這很可能是由于其分離區尺度相對較小的緣故。結果表明，干擾區分離泡的非定常特性表征為低頻膨脹/收縮過程，不同展向站位之間的峰值頻率略有差異。

圖10 分離泡面積預乘功率譜Fig.10 Pre-multiplied power spectral density for separation bubble areas

為了進一步給出各站位之間膨脹/收縮過程的相關程度，圖11還分別給出了站位分離泡面積與站位和分離泡面積的互相關系數。從圖11(a)中可以看到，站位與站位存在正的強相關性，兩者之間時間延遲約為τ=-11，而圖11(b)的結果表明，站位與站位之間則為負相關，延遲時間約為=-3/。顯然，站位分離泡的膨脹收縮過程明顯滯后于站位和。

圖11 分離泡面積互相關系數Fig.11 Cross-correlation coefficient between separation bubble areas

5 條件統計分析

Waindim等采用經驗模態分解(Empirical Mode Decomposition, EMD)對分離激波瞬時位置進行了條件分析，研究了分離泡收縮(激波往上游運動)和分離泡膨脹(激波往下游運動)的兩者情況下的流場結構特征。與Waindim等的不同之處在于，本節采用EMD方法直接對分離泡面積的時間序列進行分解，條件統計分析低頻膨脹/收縮過程中分離泡瞬態幾何特征。

首先，通過EMD方法，將原始分離泡脈動信號自適應分解為一系列具有不同特征時間尺度的本征模態函數(Intrinsic Mode Function, IMF)。隨后，選取特定的本征模態函數對分離泡的低頻膨脹/收縮過程進行了重構，并據此分別對膨脹和收縮兩個特定過程中的分離泡瞬態結構進行提取和統計。EMD方法的求解過程如下：

1) 獲得分離泡面積脈動原始時間序列()的所有局部極大值()和極小值()。

2) 采用三次樣條差值方法，依據局部極值的位置，分別構建原始信號的上包絡線()和下包絡落線()。

3) 從原始型號()中減去上下包絡線的均值()=05(()+())，判斷剩余部分()=()-()是否滿足本征模態函數的條件。如不滿足則將其代替()后再重復步驟1)～步驟3)，直至()滿足本征模態函數的條件，并將其記為分離泡面積脈動原始時間序列的第一個IMF分量(IMF)，依次類推，最后將()可分解為若干個本征模態函數和殘余項。更為詳細的經驗模態分解方法介紹可參見文獻[26]。

EMD方法具有良好的自適應頻率分辨率，無需人為給定濾波閥值。圖12分別給出了采用EMD方法對站位處分離泡面積脈動進行分解后得到9個本征模態函數(IMF～IMF)。可以清楚看到，不同的IMF都有著特定的頻段，且每個IMF包含的頻率成分和頻寬都是各不相同的，但整體來看，IMF～IMF對應為從高頻到低頻的排列，同時脈動幅值呈現增大的趨勢。

圖12 站位Z2分離泡面積本征模態函數Fig.12 IMFs of separation bubble areas at station Z2

圖13分別給出了IMF～IMF對應的預乘功率譜。可見，IMF～IMF的峰值頻率由≈4逐漸減小到≈0.01，與之前圖10中的研究結果較為符合，這也證實了EMD方法得到的本征模態函數是準確可靠的。從定量比較來看，IMF～IMF的主頻基本維持在0.1<<10范圍內，這里將其定義為高頻模態，而IMF～IMF的能量主要集中在<0.1區間，因而定義為低頻模態。圖14還給出了低頻模態和高頻模態重構信號的預乘功率譜，低頻重構信號為()=IMF()+IMF()+IMF()+IMF()+IMF()，高頻重構信號為()=()-()。如圖14所示，低頻重構信號的預乘功率譜在<0.1區間內基本與原始信號完全重合，而高頻重構信號在>0.3的范圍內則與原始信號吻合較好，這說明基于IMF～IMF的低頻重構信號能夠準確表征分離泡的低頻非定常運動過程。

圖13 本征模態函數預乘功率譜Fig.13 Pre-multiplied power spectral density for IMFs

圖14 分離泡面積脈動低頻和高頻成分預乘功率譜Fig.14 Pre-multiplied power spectral density for low-and high-frequency of separation bubble area fluctuations

在低頻重構出分離泡脈動信號后，采用式(12) 分別確定分離泡的膨脹(Dilation)和收縮(Contraction)運動：

(2)

圖15給出了站位處分離泡脈動的低頻重構信號以及提取得到的膨脹(紅色曲線)和收縮(藍色曲線)運動。可見，瞬態膨脹和收縮運動的時間尺度差異顯著。例如，在/=100，此時膨脹運動對應的時間跨度約14，而隨后的收縮運動則只持續了4，這說明此時分離泡的非定常運動以長時膨脹為主。而在/=420，收縮運動的時間尺度明顯大于前后時刻的膨脹運動，因而長時收縮運動占主導。總體來看，由于分離泡脈動處于統計定常態，因而膨脹和收縮運動在總時間上的占比基本一致。

圖15 站位Z2分離泡面積脈動條件取樣分析(紅：膨脹；藍：收縮)Fig.15 Conditional analysis of separation bubble area fluctuations at station Z2(red: dilation; blue: contraction)

圖16 分離微團瞬態云圖Fig.16 Contour of instantaneous separation micro-clusters

圖17 分離微團高度/長度比概率密度分布函數Fig.17 Probability density function of aspect-ratio for separation micro-clusters

圖18 分離微團面積和高度的散點圖(黑色：壓縮；紅色：膨脹)Fig.18 Scatter plots of separation micro-clusters area and height(black: contraction; red: dilation)

6 本征正交分解

為了進一步揭示分離泡非定常運動中的典型相干結構，采用POD方法對瞬態流向速度場進行了分析。在筆者前期的研究中，通過對展向平均后的非定常流場進行低階近似，得到了分離泡非定常演化過程中能量占優的特征模態。但由于POD分析針對的是展向平均場，其結果并不能精確反映出分離泡三維結構差異的影響規律。從本文之前的分析來看，此時分離泡沿展向存在變化劇烈的三維結構(見圖6)，因此這里將分別針對3個展向站位(～)的流向-法向剖面內非定常流向速度場開展POD分析。樣本總數為4 425，采樣區間為-3<(-)<06和0<<04。采樣時間取為0.1/，對應的可分辨范圍為0.009<<5。流向-法向剖面內非定常流向速度場(,,)的具體分解過程如下：

(3)

圖19給出了不同站位歸一化后的POD模態能量和累積能量的比較情況。這里第個模態的歸一化能量定義為：=∑，其中為第個模態的特征值。從圖19(a)中可以看到，模態能量隨著模態階數的增加而急劇降低，>100時，模態能量下降了約兩個數量級，其能量衰減率基本符合律，這與Mustafa等的研究結果是一致的。不同站位能量分布曲線也基本重合，低階模態能量略有差異。以第1階模態為例(能量貢獻最大，下文簡稱主能量模態)，站位處主能量模型的能量占比約為20.8%，略高于站位和站位處的15.6%和16.4%。從圖19(b)中累積能量分布還可以看到，站位～的前10階模態能量占比分別約為42.1%、44.9%和43.9%，這說明相較于高階模態，低階模態的貢獻明顯占優。

圖19 模態能量分布Fig.19 Mode energy distributions

為了考察POD模態的非定常特性，圖20分別給出了不同站位模態系數的預乘功率譜。可見，各站位的結果基本類似，隨著模態階數增加，占優頻率從低頻區往高頻區快速移動。對于<10的低階模態，其峰值頻率主要集中在0.01<<0.1范圍內，這與圖10中分離泡脈動的特征頻率較為接近；而對于>100的高階模態，其能譜則以>1的高頻特征為主，低頻區內則沒有明顯的脈動能量。圖21還給出了各模態時間系數()與分離泡面積()的相關系數。這里相關系數定義為

圖20 模態時間系數預乘譜云圖Fig.20 Contours of pre-multiplied spectra for time coefficient of POD modes

(4)

由式(4)可以看到，各站位主能量模態與分離泡面積脈動強相關，站位處為正相關=0.76，站位和處為負相關=-0.89和=-0.91。對于<10時，盡管相關系數急劇衰減，但兩者之間仍存在著弱相關。然而，高階模態與分離泡面積脈動的相關性系數則基本維持在零附近。研究結果表明，分離泡的低頻膨脹/收縮運動與低階模態(<10)密切相關，特別是主能量模態。

圖21 模態時間系數與分離泡面積相關系數Fig.21 Correlation between time coefficient and separation bubble area

圖22給出模態1和模態100的空間分布情況。圖中黑色虛線為時間平均聲速線，粉色實線為分離泡外邊界。從定性比較來看，盡管分離泡在展向存在復雜的三維結構，其流向和法向尺度差異顯著，但各站位主能量模態的空間結構基本一致，分離泡展向三維結構的影響可以忽略不計。如圖22所示，站位～主能量模態均對應為大尺度含能結構，模態能量集中在聲速線下方和分離泡上方的狹長區域，而在上游湍流邊界層和分離泡內并沒有出現明顯的大尺度含能結構。研究表明，主能量模態與分離泡上方剪切層密切相關，特別是其腳部區域，這與前人在平板入射激波中的研究結論是較為類似。造成這一現象的主要原因很可能是分離泡上方的剪切層仍以準二維結構特征為主，見圖22(a)～圖22(c)中聲速線的位置。高階模態的空間則與主能量模態完全不同，對應為含能較低的小尺度正負交替結構，表征了分離泡的高頻脈動過程。

圖22 POD模態空間分布Fig.22 Spatial distribution of POD modes

本文研究結果也進一步支持了激波湍流邊界層干擾低頻振蕩現象的下游機制。對于上游機制，大量學者認為低頻振蕩現象來源于上游湍流邊界層速度型、壓力脈動、擬序結構等因素，而下游機制則認為其誘因來源于下游干擾區內，如分離泡、剪切層等下游因素。從計算結果來看，低階模態與分離泡的低頻膨脹/收縮運動存在較強關聯(見圖21)，同時低階模態含能結構集中出現在分離泡上方的剪切層，而非上游湍流邊界層內。為了定量考察低階模態和分離泡低頻膨脹和收縮的內在關聯，這里采用前10個POD低階模態對各站位的非定常脈動速度場進行了重構，具體過程為

(5)

圖23分別給出了重構得到的站位～分離泡面積時間序列，為了便于比較，這里采用分離泡瞬時面積的最大值進行了歸一化處理。顯然，重構信號準確捕捉到了原始分離泡的低頻膨脹/收縮過程，兩者的差異主要體現在局部高頻脈動方面，這是由于前10階低頻模態峰值頻率主要集中在低頻區，因而其高頻區能量的貢獻則要的小得多。

圖23 前10階模態重構的分離泡面積時間序列Fig.23 Time series of separation bubble areas reconstructed from POD modes 1 to 10

7 結 論

本文采用DNS方法研究了來流馬赫數2.25，33.2°激波角的入射激波/平板湍流邊界層干擾區內分離泡特性，得到以下結論：

1) 數值模擬準確捕捉到了分離泡的三維復雜結構。分離泡流向長度明顯大于法向高度和展向寬度，同時沿展向呈現中間高兩邊低的扁平型單峰構型。

2) 數值模擬準確捕捉到了分離泡的低頻膨脹/收縮運動。與物面壓力脈動類似，分離泡低頻振蕩的時間尺度為10/～100/，其展向三維結構對峰值頻率的影響相對較小，同時分離泡兩側略滯后于中間區域。

3) 條件統計結果表明，分離泡膨脹和收縮對分離微團幾何特征沒有實質影響。研究發現，分離微團流向長度明顯大于法向高度，其高度/長度比值的概率峰值出現在0.1附近，同時分離微團面積和法向高度近似符合二次方分布。

4) POD結果表明，低階模態峰值頻率出現在低頻區，與分離泡低頻膨脹/收縮運動密切相關，而高階模態則以高頻小尺度脈動為主，總體貢獻相對較小。采用前10個低階模態可以準確重構出分離泡的低頻振蕩過程。

致 謝

感謝國家超級計算廣州中心、中國空氣動力研究與發展中心計算中心提供計算機時。