超聲速平板邊界層旁路轉(zhuǎn)捩直接數(shù)值模擬

朱海濤，單 鵬

(北京航空航天大學(xué)能源與動力工程學(xué)院，北京 100191)

超聲速平板邊界層旁路轉(zhuǎn)捩直接數(shù)值模擬

朱海濤，單 鵬*

(北京航空航天大學(xué)能源與動力工程學(xué)院，北京 100191)

采用具有七階精度的有限差分格式和基于Euler方程組的特征邊界條件，通過直接數(shù)值模擬研究了由固定相位周期性吹吸擾動誘導(dǎo)的、Ma∞=2.25、Re∞=2.5×107/m(635 000/inch)空間發(fā)展平板邊界層的旁路轉(zhuǎn)捩過程，并建立了相應(yīng)的直接數(shù)值模擬數(shù)據(jù)庫。數(shù)值結(jié)果表明:擾動相位的隨機性并不是轉(zhuǎn)捩發(fā)生的必要條件，固定相位周期性吹吸擾動即可誘發(fā)邊界層的轉(zhuǎn)捩;當(dāng)采用固定相位時，隨著吹吸頻率的增加，轉(zhuǎn)捩位置向下游移動;當(dāng)擾動頻率足夠大時，周期性吹吸無法誘導(dǎo)出邊界層的轉(zhuǎn)捩。隨后，在所得數(shù)據(jù)庫中提取了速度及熱力學(xué)變量的平均值，與已有的數(shù)據(jù)庫和理論分析吻合良好;速度脈動均方根與公認(rèn)的不可壓縮平板邊界層實驗基本吻合;熱力學(xué)變量脈動均方根和雷諾應(yīng)力與已有數(shù)據(jù)庫的統(tǒng)計結(jié)果及實驗數(shù)據(jù)具有相同的分布特征。對比發(fā)現(xiàn)，不同轉(zhuǎn)捩的誘導(dǎo)因素對平均轉(zhuǎn)捩位置和熱力學(xué)變量均方根、雷諾應(yīng)力等高階統(tǒng)計量的影響較大。

平板;超聲速;直接數(shù)值模擬;湍流邊界層;雷諾應(yīng)力;轉(zhuǎn)捩

0 引 言

可壓縮壁湍流在航空航天工程中廣泛存在?？蓧嚎s壁湍流的轉(zhuǎn)捩、非定常分離、激波及其與邊界層的相互作用等復(fù)雜流動已成為超聲速、高超聲速先進飛行器及其動力裝置設(shè)計的關(guān)鍵問題［1］。與之相關(guān)的風(fēng)洞，特別是超聲速風(fēng)洞，實驗研究受到較大的限制，不同實驗結(jié)果之間也有一定的差異。隨著高性能計算機及高精度數(shù)值格式的發(fā)展，直接數(shù)值模擬已成為探索復(fù)雜流動機理，發(fā)現(xiàn)新的流動現(xiàn)象及流動規(guī)律的重要工具。直接數(shù)值模擬不僅可以給出各個流動時刻清晰的流動圖畫，而且可以計算包括脈動運動在內(nèi)的所有湍流瞬時運動量在三維空間中的時空演化過程。因此，可以通過建立可靠的直接數(shù)值模擬數(shù)據(jù)庫來很好地彌補上述實驗研究的不足。

轉(zhuǎn)捩主要可分為兩種類型［2］:自然轉(zhuǎn)捩和旁路轉(zhuǎn)捩。前者為小擾動誘發(fā)的轉(zhuǎn)捩，經(jīng)過感受性階段、不穩(wěn)定波的線性增長階段和非線性增長階段，層流破碎轉(zhuǎn)捩到湍流。后者為較大的局部干擾導(dǎo)致的轉(zhuǎn)捩，邊界層內(nèi)的不穩(wěn)定波直接進入非線性增長階段，層流轉(zhuǎn)捩為湍流。目前，對第二類轉(zhuǎn)捩知之甚少，主要依靠經(jīng)驗處理，且風(fēng)洞特別是超聲速風(fēng)洞實驗研究比較有限［3］。

在直接數(shù)值模擬過程中，為了促進轉(zhuǎn)捩的發(fā)生，通常需要添加一些人工擾動。常采用的擾動形式有:在壁面上增加隨機吹吸擾動或具有隨機相位的周期性吹吸擾動、在入口處增加周期性尾跡或格柵湍流形式的擾動。這兩種擾動經(jīng)常用來誘導(dǎo)旁路轉(zhuǎn)捩的發(fā)生。對于自然轉(zhuǎn)捩，通常根據(jù)線性穩(wěn)定性理論或線性Navier-Stokes(N-S)方程組來確定所需的小擾動，例如在入口處增加二維或三維Tollmien-Schlichting(TS)波形式的擾動等。在外流中，通常采用壁面上的吹吸擾動誘導(dǎo)旁路轉(zhuǎn)捩的發(fā)生。例如，Rai［4］等采用高精度有限差分格式實現(xiàn)了Ma∞=2.25的平板邊界層流動的直接數(shù)值模擬，發(fā)現(xiàn)速度脈動均方根的分布與不可壓縮平板邊界層湍流實驗結(jié)果很吻合，從而驗證了Morkovin假設(shè)的正確性;Rai采用的擾動即為在平板壁面上增加具有隨機相位的周期性吹吸擾動。采用相同的擾動，Gatski等［5］和Pirozzoli等［6］分別對超聲速平板邊界層湍流的高階統(tǒng)計量和湍動能的生產(chǎn)和耗散過程進行了細(xì)致的研究。在入口處增加尾跡或格柵湍流形式的擾動在內(nèi)流中比較常用。例如Rodi和Wissink等［7-8］采用直接數(shù)值模擬的方法研究了在周期性尾跡作用下低壓渦輪葉柵中分離流、分離泡、轉(zhuǎn)捩等復(fù)雜的流動現(xiàn)象，并研究了這些流動現(xiàn)象對渦輪葉片表面熱傳導(dǎo)效應(yīng)的影響，數(shù)值結(jié)果和實驗結(jié)果吻合良好。Zaki等［9］采用基于交錯網(wǎng)格的有限體積法，通過求解不可壓縮N-S方程組，對自由湍流作用下和周期性尾跡作用下V103葉柵流場的非定常分離及旁路轉(zhuǎn)捩過程進行了直接數(shù)值模擬，數(shù)值結(jié)果反應(yīng)的平均流場、分離及轉(zhuǎn)捩等特征與實驗結(jié)果是一致的。Egorov等［10］采用根據(jù)線性N-S方程組得到小擾動波，研究了超聲速平板邊界的感受性問題。國內(nèi)陳健［11］等完成了周期性吹吸擾動作用下光滑圓管流動轉(zhuǎn)捩過程的直接數(shù)值模擬。傅德熏等［1，12-14］采用具有隨機相位的周期性吹吸擾動完成馬赫數(shù)Ma∞=2.25和Ma∞=6的平板邊界層湍流的直接數(shù)值模擬，并研究了Ma∞=2.25平板湍流邊界層中的聲效應(yīng)，以及高馬赫數(shù)下湍動能的生成和耗散機制。李新亮等［15］采用隨機吹吸擾動對Ma∞=6的鈍錐繞流進行了直接數(shù)值模擬，輔以線性穩(wěn)定性理論分析了鈍錐表面局部區(qū)域轉(zhuǎn)捩位置相對靠后的原因。

盡管已經(jīng)得到了若干直接數(shù)值模擬的數(shù)據(jù)庫，但是仍然需要新的、采用不同數(shù)值方法、邊界條件和轉(zhuǎn)捩誘導(dǎo)因素的直接數(shù)值模擬結(jié)果對已有的數(shù)據(jù)庫進行驗證和檢驗，以期對所建立的DNS數(shù)據(jù)庫和轉(zhuǎn)捩的發(fā)生有較為全面的認(rèn)識。本文采用具有七階精度的有限差分格式和基于Euler方程組的特征邊界條件［16］，以固定相位周期性吹吸為轉(zhuǎn)捩的誘導(dǎo)因素，對Ma∞=2.25、Re∞=2.5×107/m的空間發(fā)展平板邊界層湍流進行了直接數(shù)值模擬，并研究了吹吸擾動的相位和頻率對轉(zhuǎn)捩發(fā)生的影響。

1 數(shù)值方法及問題描述

考慮如下三維非定常N-S方程組:

補充完全氣體狀態(tài)方程、溫度與總能量、熱傳導(dǎo)系數(shù)和Prandtl數(shù)之間的關(guān)系式及粘性系數(shù)滿足的Sutherland公式，

上述公式中用到的常數(shù)值如下:

至此，完成方程組(1)的封閉過程。在下文的計算過程中，采用了在一個方向上進行加密的非均勻網(wǎng)格。此時，需要在方程組(1)中的對應(yīng)方向上無粘和粘性空間導(dǎo)數(shù)項乘以該方向上的網(wǎng)格變換系數(shù)。此外，粘性應(yīng)力張量和熱通量中該方向上的空間導(dǎo)數(shù)項也要做類似的變換，此處不再詳述。本文計算過程中參考長度L0=0.0254 m(1 inch)，基于L0的雷諾數(shù)Re∞= 635 000;自由來流的馬赫數(shù)Ma∞=2.25，溫度T∞= 169.44 K。上述方程組根據(jù)自由來流的密度ρ"、速度V∞、溫度T∞及上述參考長度L0進行無量綱化，并以方程組(1)進行直接數(shù)值模擬。

對式(1)中的無粘項，首先采用Steger-Warming方法進行流通矢量分裂，然后對正負(fù)流通矢量均采用守恒型的七階迎風(fēng)差分格式計算其空間導(dǎo)數(shù)。粘性項采用具有八階精度的中心差分格式離散。最后通過具有TVD性質(zhì)的三步三階Runge-Kutta公式進行時間積分。

本文的直接數(shù)值模擬分兩步進行，首先進行二維層流計算，然后以所得的穩(wěn)態(tài)二維層流解為初場，通過在壁面增加固定相位的周期性吹吸擾動進行三維邊界層轉(zhuǎn)捩的直接數(shù)值模擬。二維和三維計算采用的數(shù)值方法是相同的。

二維計算區(qū)域如圖1所示，平板長度取為Lx=0.281 m，入口取在平板前緣前方0.0762 m的地方，上邊界與平板的距離為Ly=0.127 m。入口邊界條件為自由來流邊界;平板壁面為等溫壁面條件，Tw=322 K;從入口到平板前緣采用對稱邊界條件(圖中虛線);平板上方邊界采用基于Euler方程組的特征邊界條件;出口邊界當(dāng)流向馬赫數(shù)小于1時采用特征邊界條件，反之，由插值獲取出口邊界上的流場信息;二維計算的初場取為自由來流。在進行一定時間的流場推進后即可獲得穩(wěn)態(tài)的二維層流解。

圖1 二維計算區(qū)域Fig.1 Two-dimensional calculation area

在三維計算中，平板長度取為Lx=0.14224 m，法向高度取為 Ly=0.00508 m，展向長度為 Lz= 0.00889 m。壁面條件與二維計算一致;展向采用周期性邊界條件;壁面上方邊界采用如下基于Euler方程組的特征邊界條件其中下標(biāo)b、e、"分別表示邊界處、在邊界處進行外插和自由來流對應(yīng)物理量的值;n為邊界處指向計算域內(nèi)部的單位法向量;c為當(dāng)?shù)芈曀?γ=1.4。出口邊界條件與二維計算一致:若流向馬赫數(shù)小于1.0，則由式(3)所示的特征邊界條件確定，其中，pb為沿法向最后一個流向馬赫數(shù)小于1.0的網(wǎng)格點處的壓力值;若流向馬赫數(shù)大于等于1.0則直接進行外插確定出口邊界上的流場信息。

入口邊界和初場的設(shè)置如下所述:獲得二維穩(wěn)態(tài)結(jié)果后，取二維平板Lx=0.1778m處的數(shù)值解(從平板前緣算起，此處已遠(yuǎn)離平板前緣的激波)沿展向復(fù)制作為三維計算的入口邊界條件，展向速度為0，然后把入口邊界處的物理量沿流向復(fù)制，作為三維直接數(shù)值模擬的初場。入口邊界條件在三維瞬態(tài)計算過程中保持不變。

計算網(wǎng)格在流向和展向為均勻網(wǎng)格，分別取Nx=1681和Nz=320個網(wǎng)格點。法向網(wǎng)格沿壁面按雙曲正切函數(shù)進行加密，網(wǎng)格數(shù)為Ny=71。本文三維計算的網(wǎng)格數(shù)及無量綱化的空間步長和相關(guān)文獻的對比如表1所示。相較于其他文獻，本文在流向和展向上的空間步長較大，但由于本文采用高精度的計算格式，在一定程度上可以降低對空間步長的要求［2］。

為了加快層流到湍流的轉(zhuǎn)捩，計算過程中在平板流向Lx=0.0127～0.0254 m的區(qū)域添加了如式(4)所示的固定相位周期性吹吸擾動:vbs為平板壁面法向的吹吸速度;xa、xb為擾動區(qū)域沿流向上的起始和終止坐標(biāo);φl、φm分別為擾動的空間和時間相位，β為擾動的基頻。在本文的計算中，A=0.04，xa=0.5，xb=1.0，lmax=10，mmax=5。φl、φm取固定值，它們和β的取值在后文詳細(xì)給出。

表1 本文及相關(guān)文獻的網(wǎng)格參數(shù)Table 1 Non-dimensional mesh parameters of present and references

2 直接數(shù)值模擬結(jié)果

2.1 擾動相位和頻率的影響

對于周期性的吹吸擾動，通常時間和空間相位φl、φm均?。?，1］內(nèi)的隨機數(shù)。在本文的直接數(shù)值模擬過程中，采用固定相位的吹吸擾動，即誘導(dǎo)轉(zhuǎn)捩發(fā)生的擾動不再具有隨機性。為了較為快速地檢驗不同的時間和空間相位組合能否導(dǎo)致轉(zhuǎn)捩的發(fā)生，本小節(jié)的計算中采用了較為稀疏的網(wǎng)格:流向上取Nx= 1121，展向上取Nz=160，法向上仍取Ny=71，吹吸頻率β=2.5。圖2給出了i)φl=0，φm=0;ii)φl=1/4，φm=0;iii)φl=0，φm=1/2三種相位組合情形下，平板表面摩擦系數(shù)的時均分布圖。在ii)和iii)情形下，時間和空間擾動中幅值最大的擾動波的相位分別相差π/2、π?？梢钥闯鋈N相位組合均導(dǎo)致了轉(zhuǎn)捩的發(fā)生。圖3給出了z=0.0055截面上轉(zhuǎn)捩過程中和充分發(fā)展湍流場中展向渦量的分布云圖。可以看出經(jīng)過邊界層平均流和周期性擾動的相互作用，誘導(dǎo)出了不穩(wěn)定的剪切層;在x=2.2附近，該剪切層失穩(wěn)導(dǎo)致了層流的破碎;破碎的邊界層不斷向下游發(fā)展，逐漸形成充分發(fā)展的湍流場，如圖3所示。因此，對于周期性吹吸擾動，相位的隨機性并不是轉(zhuǎn)捩發(fā)生的必要條件，湍流場所具有的隨機性是其內(nèi)在的本質(zhì)特征，與擾動的隨機性并無直接的聯(lián)系。

圖2 擾動相位對平板表面時均摩擦系數(shù)的影響Fig.2 Time averaged surface friction with different time and space disturbance phases

圖3 不同時刻z=0.0055截面上展向渦量的分布云圖Fig.3 Distribution of spanwise vorticity on the plane of (x，y)at z=0.0055，t=4.2 and t=8.0

圖4給出了情形i)下，不同吹吸頻率對平板表面平均摩擦系數(shù)的影響?？梢园l(fā)現(xiàn)，當(dāng)β=3.5時，摩擦系數(shù)不斷地減小，與層流對應(yīng)的摩擦系數(shù)值很相近，表明邊界層內(nèi)沒有轉(zhuǎn)捩的發(fā)生。β=2.5、1.5時，摩擦系數(shù)在平均轉(zhuǎn)捩位置處突然變大，隨后保持在一個顯著地大于層流摩擦系數(shù)的水平。顯然β=1.5對應(yīng)的轉(zhuǎn)捩位置比較靠前。因此，當(dāng)吹吸頻率增大時，轉(zhuǎn)捩位置向平板的下游移動，但當(dāng)頻率繼續(xù)增大時，固定相位的周期性吹吸擾動無法誘導(dǎo)出轉(zhuǎn)捩的發(fā)生，這與高慧等人［13］的結(jié)論是一致的。

圖4 不同吹吸頻率下平板表面平均摩擦系數(shù)的分布Fig.4 Averaged friction with different frequencies of blow and suction disturbances

2.2 統(tǒng)計結(jié)果分析

下面給出固定相位周期性吹吸擾動作用下，在湍流充分發(fā)展區(qū)所得的平均及高階統(tǒng)計量與相關(guān)參考文獻、理論分析和實驗數(shù)據(jù)的對比結(jié)果。在以下的計算過程中，取φl=0，φm=0，β=2.5。

由相關(guān)文獻的計算可知，在流向坐標(biāo)0.12192 m處，邊界層已成為充分發(fā)展的湍流邊界層。圖5給出了y+=5.1323，流向坐標(biāo)為0.12192 m處密度、速度和溫度脈動的二階自相關(guān)量的統(tǒng)計平均值沿展向的分布?？梢钥闯?，各個二階自相關(guān)量的統(tǒng)計平均值均趨向于0。因此，本文計算域展向長度是足夠的。

圖5 二階自相關(guān)量統(tǒng)計平均值沿展向的分布Fig.5 Correlation tensor distribution along spanwise

圖6(a)為t=8.0時，流場z=0.0055處x-y截面上密度分布云圖;上方為層流邊界層失穩(wěn)及破碎時的密度分布;下方為充分發(fā)展的湍流區(qū)中密度的分布云圖。可見，盡管所添加的擾動不再具有隨機性，但邊界層中仍然誘導(dǎo)出了非線性增長的不穩(wěn)定波，使層流逐漸失穩(wěn)、破碎，最終發(fā)生轉(zhuǎn)捩，發(fā)展為充分發(fā)展的湍流邊界層。在充分發(fā)展的湍流區(qū)中可以發(fā)現(xiàn)低密度流體的上拋和高密度流體的下卷，并由此形成銳利的湍流邊界。圖6(b)和圖6(c)分別是在流向坐標(biāo)為0.12192m處熱力學(xué)變量的統(tǒng)計平均值沿法向的分布和經(jīng)過Van Driest變換［2］后的平均速度剖面。圖6 (b)中的參考值為Gatski等［5］人的直接數(shù)值模擬結(jié)果?？梢园l(fā)現(xiàn)，二者基本吻合。圖6(c)中的兩個參考值分別為符合不可壓縮平板邊界層壁面律的= y+和不可壓縮平板邊界層對數(shù)律的=1/0.41 ln(y+)+5.1，其中為經(jīng)過Van Driest變換后的速度?？梢钥闯?，本文直接數(shù)值模擬結(jié)果與理論分析符合甚好。由上可知:盡管本文采用固定相位周期性吹吸擾動，誘導(dǎo)轉(zhuǎn)捩發(fā)生的擾動不再具有隨機性，但所得湍流場與具有隨機性的擾動誘導(dǎo)的湍流場具有相同的平均特征。

圖6 邊界層轉(zhuǎn)捩及湍流充分發(fā)展區(qū)熱力學(xué)變量和速度平均值的統(tǒng)計分析Fig.6 Density contours and the statistics of thermodynamical variables and velocity

圖7 摩擦系數(shù)的分布Fig.7 Averaged surface friction distribution

圖7給出了本文平板表面平均摩擦系數(shù)的分布與層流和湍流邊界層理論預(yù)測值及相關(guān)參考文獻計算結(jié)果的對比?？梢园l(fā)現(xiàn)，在層流區(qū)本文的計算結(jié)果與層流邊界層理論預(yù)測值符合甚好，但在湍流區(qū)本文的計算結(jié)果大于湍流邊界層理論的預(yù)測值。這些差異與他人的直接數(shù)值模擬結(jié)果是一樣的。與相關(guān)文獻結(jié)果相比，本文的計算結(jié)果更接近理論預(yù)測值，原因可能與近壁區(qū)網(wǎng)格的分辨率相關(guān)。

摩擦系數(shù)曲線從層流值開始突然顯著地增大，然后慢慢減小，逐漸靠近湍流邊界層理論的預(yù)測值。轉(zhuǎn)捩過程即發(fā)生在摩擦系數(shù)突然增大的區(qū)域?？梢钥闯?，三個直接數(shù)值模擬結(jié)果的平均轉(zhuǎn)捩位置有明顯的差異，本文的轉(zhuǎn)捩位置最為靠后。Rai［4］和Gao［13］的結(jié)果均是基于具有隨機相位周期性擾動的，由此可知:不同的轉(zhuǎn)捩誘導(dǎo)方式對轉(zhuǎn)捩平均位置具有一定的影響。

圖8是在流向坐標(biāo)為0.12192 m處湍流充分發(fā)展區(qū)速度和熱力學(xué)變量脈動均方根(Root Mean Square，RMS)沿法向的分布曲線。圖8(a)中的參考數(shù)據(jù)為不可壓縮平板邊界層的實驗結(jié)果［17］?？梢园l(fā)現(xiàn)，盡管自由流的Ma∞=2.25，但在平板邊界層內(nèi)，速度脈動均方根沿法向的分布與不可壓縮平板邊界層的分布特征十分接近。說明當(dāng)來流馬赫數(shù)為2.25時，流場的可壓縮性對速度脈動的影響也是比較小的。此時，可壓縮效應(yīng)對湍流脈動動力學(xué)特征的直接影響較弱，其影響主要反映在對流場平均流動特征的改變。因此，通過可以把平均密度變化影響考慮進去的Van Driest變換，流場的平均速度分布與不可壓縮平板邊界層的壁面律和對數(shù)律速度剖面是相吻合的，Morkovin假設(shè)也是成立的。圖8(b)中的參考值為Gatski等［5］人直接數(shù)值模擬結(jié)果。可以發(fā)現(xiàn)，盡管本文熱力學(xué)變量的統(tǒng)計平均值與Gatski等的計算結(jié)果符合甚好，但這些變量脈動均方根值明顯小于文獻的計算結(jié)果，但二者的分布趨勢是一致的。由此可以推測:與轉(zhuǎn)捩的平均位置相似，不同的轉(zhuǎn)捩誘導(dǎo)方式可以影響到熱力學(xué)變量脈動均方根的大小。

圖8 速度及熱力學(xué)變量脈動均方根的分布Fig.8 Distributions of the RMS of the fluctuations of velocity and thermodynamics variables

圖9是雷諾正應(yīng)力和剪切應(yīng)力的統(tǒng)計結(jié)果和相關(guān)的實驗數(shù)據(jù)和參考文獻計算結(jié)果的對比。圖9(a)的參考值為文獻［17］不可壓縮平板邊界層湍流DNS的結(jié)果。可以發(fā)現(xiàn)，流向雷諾正應(yīng)力明顯大于不可壓縮流直接數(shù)值模擬的結(jié)果，法向和展向分量的差異相對較小。圖9(b)的參考值為文獻［17］中雷諾剪切應(yīng)

圖9 雷諾應(yīng)力的分布Fig.9 Distributions of Reynolds stress in normal direction

力直接數(shù)值模擬結(jié)果和公認(rèn)的不可壓縮平板邊界層湍流實驗數(shù)據(jù)［17］。圖9(c)的參考值分別為公認(rèn)的不可壓縮平板邊界層湍流實驗數(shù)據(jù)［17］和文獻［4］直接數(shù)值模擬的計算結(jié)果。圖9(b)、圖9(c)中的實驗數(shù)據(jù)均沒有考慮密度的影響。可以看出，本文的直接數(shù)值模擬的計算結(jié)果稍小于兩個參考文獻DNS的結(jié)果，但三個DNS結(jié)果的分布趨勢是一致的，說明直接數(shù)值模擬采用的數(shù)值方法可以影響到雷諾應(yīng)力的計算。與不可壓縮平板邊界層湍流實驗值相比，考慮密度影響的雷諾剪切應(yīng)力的峰值大于實驗值，反之，則小于實驗值。相對而言，考慮密度變化后的剪切應(yīng)力分布與實驗吻合較好。

3 結(jié) 論

采用高精度有限差分格式和基于Euler方程組的特征邊界條件，通過直接數(shù)值模擬研究了周期性吹吸擾動的相位和頻率對超聲速平板邊界層旁路轉(zhuǎn)捩的影響，并建立了由固定相位周期性吹吸擾動誘導(dǎo)的、Ma∞=2.25、Re∞=2.5×107/m超聲速平板邊界層湍流的直接數(shù)值模擬數(shù)據(jù)庫。數(shù)值計算結(jié)果表明:

(1)對于周期性吹吸擾動，相位的隨機性不是導(dǎo)致旁路轉(zhuǎn)捩發(fā)生的必要條件;采用固定相位，盡管擾動不再具有隨機性，但邊界層內(nèi)旁路轉(zhuǎn)捩仍然會發(fā)生。

(2)采用固定相位周期性吹吸擾動，隨著吹吸頻率的增大，轉(zhuǎn)捩的平均位置向下游移動;但頻率超過某一臨界值時，邊界層不再發(fā)生轉(zhuǎn)捩。

(3)無論采用隨機相位或固定相位，湍流的平均特征是一致的;二階統(tǒng)計量具有相同的變化趨勢。相對于平均參數(shù)而言，平均轉(zhuǎn)捩位置和二階統(tǒng)計量對轉(zhuǎn)捩的誘導(dǎo)因素、數(shù)值方法和邊界條件較為敏感。

［1］ Li Xinliang，F(xiàn)u Dexun，Ma Yanwen.Direct numerical simulation of compressible turbulence flows［J］.Acta Mechanica Sinica，2010，26(6):795-806.

［2］ Fu Dexun，Ma Yanwen，Li Xinliang，et al.Direct numerical simulation of compressible turbulent flow［M］.Beijing:Science Press，2010.(in Chinese)

傅德熏，馬延文，李新亮，等.可壓縮湍流直接數(shù)值模擬［M］.北京:科學(xué)出版社，2010.

［3］ Egorov I V，F(xiàn)edorov A V，Soudakov V G.Direct numerical simulation of unstable disturbances in supersonic boundary layer［R］.AIAA 2004-588.

［4］ Rai M M，Gatski T B，Erlebacher G.Direct numerical simulation of spatially evolving compressible turbulent boundary layers［R］.AIAA 95-0583.

［5］ Gatski T B，Erlebacher G.Numerical simulation of a spatially evolving supersonic turbulent boundary layer［R］.NASA TM 2002-211934.

［6］ Pirozzoli S，Grasso F.Direct numerical simulation and analysis of a spatially evolving supersonic turbulent boundary layer at M=2.25［J］.Physics of Fluids，2004，16(3):530-545.

［7］ Rodi W.DNS and LES of some engineering flows［J］.Fluids Dynamics Research，2006，38(2-3):145-173.

［8］ Wissink J G.DNS of separation，low Reynolds number flow in a turbine cascade with incoming wakes［J］.Engineering turbulence modeling and experiments，2002，(5):731-740.

［9］ Zaki T A，Wissink J G，Durbin P A，et al.Direct computations ofboundary layers distorted by migrating wakes in a linear compressor cascade［J］.Flow Turbulence Combust，2009，83:307-322.

［10］Egorov I V，F(xiàn)edorov A V，Nechaev A V.Receptivity of supersonic boundary layer on a blunt plate to acoustic disturbances［R］.AIAA 2004-249.

［11］ Chen Jian，Cui Guixiang，Xu Chunxiao，et al.Direct numerical simulation of spatial evolutional transition in pipe flow［J］.Chinese Journal of Theoretical and Applied Mechanics，2003，35(1):6-13.(in Chinese)

陳健，崔桂香，許春曉，等.空間發(fā)展圓管轉(zhuǎn)捩的直接數(shù)值模擬［J］.力學(xué)學(xué)報，2003，35(1):6-13.

［12］Li Xinliang，F(xiàn)u Dexun，Ma Yanwen，et al.Acoustic calculation for supersonic turbulent boundary layer flow［J］.Chinese Physics Letters，2009，26(9):094701.

［13］Gao Hui，F(xiàn)u Dexun，Ma Yanwen，et al.Direct numerical simulation of supersonic turbulent boundary layer flow［J］.Chinese Physics Letters，2005，22(7):1709-1702.

［14］Li Xinliang，F(xiàn)u Dexun，Ma Yanwen.Direct numerical simulation of spatially evolving supersonic turbulent boundary layer at Ma=6［J］.Chinese Physics Letters，2006，23(6):1519-1522.

［15］Li Xinliang，F(xiàn)u Dexun，Ma Yanwen.Direct numerical simulation of hypersonic boundary layer transition over a blunt cone with a small angle of attack［J］.Physics of Fluids，2010，22(2):1-18.

［16］Hirsch C.Numerical computation of internal and external flows，Vol 2:computational methods for inviscid and viscous flows［M］.Chichester:John Wiley＆Son，1990.

［17］Spalart P R.Direct numerical simulation of a turbulent boundary layer up to Reθ=1410［J］.Journal of Fluid Mechanics，1988，187: 61-98.

Direct numerical simulation of by-pass transition of the supersonic flat plate turbulent boundary layer

Zhu Haitao，Shan Peng*
(School of Jet Propulsion，Beihang University，Beijing 100191，China)

Using the characteristic boundary conditions and high-order finite difference schemes，direct numerical simulation(DNS)was performed to investigate the effect of phase and frequency of the periodic blow and suction disturbance which was used to trigger the by-pass transition of the flat plate boundary layer with free stream Ma∞=2.25，Re∞=2.5×107/m(635 000/inch).The numerical results show that the randomness of the disturbance phase is not the necessary condition to trigger the transition，and with the deterministic disturbance phase，the average location of transition moves backward when the disturbance frequency increases.However，if the frequency is bigger than a critical value，the transition in the boundary layer will not happen.A DNS database was established under the condition of deterministic disturbance phase.The Van Direst velocity agrees well with the wall and logarithm law of incompressible flat plate boundary layer，and the turbulence intensities coincide with the incompressible flat plate experimental datum.The averaged thermodynamic variables are very close to that of reference results.These verify the present DNS results.Then，the statistics of thermodynamic variables fluctuation and Reynolds stress are extracted from the DNS database.Compared with the relative references'DNS results，the disturbance，boundary conditions and numerical methods have observably effect on the second order statistics and the averaged transition position.

flat;supersonic;DNS;turbulent boundary layer;Reynolds stress;transition

0258-1825(2015)03-0345-08

O357.4;V211.3

A

10.7638/kqdlxxb-2013.0083

2013-08-07;

2013-11-07

朱海濤(1982-)，男，山西省襄汾縣人，博士研究生，研究方向:計算流體力學(xué)，大規(guī)模科學(xué)和工程計算.E-mail:zht83@126.com

單鵬*(1955-)，男，北京人，教授，研究方向:航空發(fā)動機，葉輪機械和計算空氣動力學(xué).E-mail:pshan@buaa.edu.cn

朱海濤，單鵬.超聲速平板邊界層旁路轉(zhuǎn)捩直接數(shù)值模擬［J］.空氣動力學(xué)學(xué)報，2015，33(3):345-352.

10.7638/kqdlxxb-2013.0083 Zhu H T，Shan P.Direct numerical simulation of by-pass transition of the supersonic flat plate turbulent boundary layer［J］.Acta Aerodynamica Sinica，2015，33(3):345-352.