粗糙元誘導(dǎo)的高超聲速邊界層轉(zhuǎn)捩

段志偉， 肖志祥

清華大學(xué) 航天航空學(xué)院， 北京　100084

?

粗糙元誘導(dǎo)的高超聲速邊界層轉(zhuǎn)捩

段志偉， 肖志祥*

清華大學(xué) 航天航空學(xué)院， 北京100084

基于有限體積方法，直接數(shù)值模擬了高超聲速邊界層內(nèi)不同形狀粗糙元導(dǎo)致的強制轉(zhuǎn)捩現(xiàn)象；為了能夠深入探究強制轉(zhuǎn)捩機理，解析小尺度運動，同時又能夠較好地捕捉激波，采用高階色散最小耗散可調(diào)(MDCD)格式對Navier-Stokes方程組對流項進(jìn)行重構(gòu)。計算結(jié)果表明，數(shù)值結(jié)果與對應(yīng)的實驗值吻合較好；該方法能解析小尺度的流動結(jié)構(gòu)以及規(guī)則結(jié)構(gòu)的破碎與失穩(wěn)過程，可揭示粗糙元引起的強制轉(zhuǎn)捩機理，即此類強制轉(zhuǎn)捩主要由粗糙元頂部的三維剪切層失穩(wěn)導(dǎo)致。對多種粗糙元的轉(zhuǎn)捩效果進(jìn)行了定量研究，影響因素包括粗糙元形狀、幾何參數(shù)等。

高超聲速； 粗糙元； 邊界層轉(zhuǎn)捩； 參數(shù)化影響； 直接數(shù)值模擬

在高超聲速飛行器繞流中，邊界層轉(zhuǎn)捩是非常重要的流動現(xiàn)象。它主要影響飛行器的摩擦阻力、氣動加熱、進(jìn)氣道質(zhì)量捕獲、發(fā)動機內(nèi)燃料摻混和燃燒等[1]。粗糙元可以改變邊界層基本流動剖面，甚至產(chǎn)生較大擾動，對邊界層流動的轉(zhuǎn)捩產(chǎn)生較大影響；由于粗糙元改變了流動狀態(tài)，其自身也面臨比較嚴(yán)峻的氣動熱環(huán)境等。粗糙元誘導(dǎo)轉(zhuǎn)捩的影響因素非常多，對其轉(zhuǎn)捩機理與具體效果的研究尚未十分明確，是目前國際上的研究熱點之一。因此本文針對粗糙元形狀、幾何參數(shù)等，對其誘導(dǎo)轉(zhuǎn)捩的效果進(jìn)行了定量研究。

大多數(shù)情況下，粗糙元會促使轉(zhuǎn)捩提前，并且隨著粗糙元高度的增加，其促進(jìn)轉(zhuǎn)捩效果會更明顯。早在20世紀(jì)60年代， Van Driest等[2]對粗糙元在馬赫數(shù)Ma=2.71的超聲速邊界流動中進(jìn)行了觸發(fā)轉(zhuǎn)捩研究，通過大量的試驗數(shù)據(jù)擬合得到轉(zhuǎn)捩位置與單位雷諾數(shù)的關(guān)系曲線，他們發(fā)現(xiàn)粗糙元存在著“臨界高度”與“有效高度”。雖然該試驗是在常規(guī)風(fēng)洞內(nèi)進(jìn)行，風(fēng)洞噪聲比較大，所獲得的轉(zhuǎn)捩曲線與實際情況稍有差別，但是轉(zhuǎn)捩趨勢相同，該試驗對后續(xù)的粗糙元研究有極大的啟發(fā)意義。在近期的風(fēng)洞試驗中，趙慧勇等[3]也發(fā)現(xiàn)了粗糙元的“臨界高度”和“有效高度”，有效高度大約為邊界層厚度的0.7～0.8倍。

粗糙元高度是影響其轉(zhuǎn)捩的重要因素之一。當(dāng)然，馬赫數(shù)、雷諾數(shù)、來流湍流度、粗糙元形狀等對轉(zhuǎn)捩效果也有較大影響。試驗研究表明轉(zhuǎn)捩雷諾數(shù)隨馬赫數(shù)增加而增加[4]，這可能與高馬赫數(shù)對剪切層的穩(wěn)定作用有關(guān)。即便在粗糙元足夠高時，來流湍流度對轉(zhuǎn)捩效果仍然有很大影響。Casper等[5]的靜音風(fēng)洞試驗表明，在靜音條件下的轉(zhuǎn)捩距離比噪聲條件下增長了一倍左右。

粗糙元的形狀不僅影響轉(zhuǎn)捩效果，同時對流場中的擾動產(chǎn)生、氣動力熱環(huán)境等也有著較大影響。理想的粗糙元應(yīng)該能夠很快觸發(fā)轉(zhuǎn)捩，并且對流場的擾動要盡量小，同時自身不產(chǎn)生過大的阻力和熱流，尤其在粗糙元上的熱流不至于過大，否則易燒毀。因此，關(guān)于粗糙元形狀對轉(zhuǎn)捩的影響研究由來已久。

一般來說，對于高超聲速邊界層轉(zhuǎn)捩，三維粗糙元效果要優(yōu)于二維粗糙元[6]。對于不同形狀的三維粗糙元，其有著類似的流動拓?fù)浣Y(jié)構(gòu)：在粗糙元上游會形成幾對反轉(zhuǎn)的渦，在正后方形成不穩(wěn)定尾跡渦和側(cè)邊的馬蹄渦(或次渦)結(jié)構(gòu)[7-13]。如Danehy等[11]采用一氧化氮平面激光可視技術(shù)顯示了詳細(xì)的尾跡渦結(jié)構(gòu)隨著粗糙元鈍度增加，其上游分離渦與表面之間的距離也增大，轉(zhuǎn)捩效果更加明顯。當(dāng)然綜合考慮轉(zhuǎn)捩效果和其對流場的影響后，斜坡形單元可能是比較好的選擇[8]。

風(fēng)洞試驗可以研究粗糙單元轉(zhuǎn)捩，但是由于受到顯示技術(shù)的限制，其給出的流動數(shù)據(jù)較少。隨著計算機技術(shù)的發(fā)展，大規(guī)模的數(shù)值計算已成為可能。采用數(shù)值計算方法可以給出風(fēng)洞試驗所不能觀測到的更豐富的流動結(jié)構(gòu)，近年來很多學(xué)者對粗糙元導(dǎo)致的高超聲速邊界層轉(zhuǎn)捩進(jìn)行了數(shù)值模擬。Bernardini等[14]的數(shù)值計算表明，當(dāng)雷諾數(shù)足夠高時，在所有馬赫數(shù)工況下(0，0.25，1.1，2，3，4)粗糙元尾跡區(qū)都會出現(xiàn)不穩(wěn)定的發(fā)卡渦結(jié)構(gòu)，而這正是充分發(fā)展湍流的顯著特征之一。他們的計算結(jié)果表明，壓縮性對流動結(jié)構(gòu)的影響并不明顯。雖然壓縮性對尾跡區(qū)的流動結(jié)構(gòu)影響不大，但是在高超聲速工況下，當(dāng)粗糙元的高度超過聲速線時，粗糙元上游會形成明顯的脫體激波，其前方的分離泡也會導(dǎo)致分離激波的出現(xiàn)。Bartkowicz等[15]的計算結(jié)果清晰地展示了粗糙元前的激波與渦結(jié)構(gòu)；粗糙元前的分離渦隨主流向下游發(fā)展形成了極強的馬蹄渦和尾跡渦結(jié)構(gòu)。這些數(shù)值模擬結(jié)果與風(fēng)洞試驗[10]所展示的流場結(jié)構(gòu)幾乎完全相同。

對于粗糙元促進(jìn)高超聲速邊界層轉(zhuǎn)捩的研究已經(jīng)持續(xù)超過了50年，而且多種形式的粗糙元也已經(jīng)被成功應(yīng)用到高超聲速飛行器上，但是對于其促進(jìn)轉(zhuǎn)捩的機理卻一直沒有統(tǒng)一完整的解釋。 Schneider[4]在他的綜述性文章中將粗糙元的作用總結(jié)為：① 粗糙元會在尾跡區(qū)產(chǎn)生流向渦結(jié)構(gòu)和不穩(wěn)定的剪切層。當(dāng)粗糙元大于“有效高度”后，轉(zhuǎn)捩會在下游立即發(fā)生；當(dāng)其尺寸較小時，轉(zhuǎn)捩距離可能會在更下游位置發(fā)生，尾跡失穩(wěn)可能是由剪切層失穩(wěn)或者尾跡渦失穩(wěn)主導(dǎo)的。② 當(dāng)粗糙元尺寸較小時(基于粗糙元高度h的雷諾數(shù)Reh<10)，粗糙元后的流向渦可能通過一些失穩(wěn)機制發(fā)展，例如橫流、G?rtler渦或者瞬態(tài)增長等。③ 粗糙元也可能通過感受性機理，與聲波或者其他的來流擾動相互作用進(jìn)而發(fā)展出失穩(wěn)波。

對于大尺寸的粗糙元，Iyer和Mahesh[16]對平板上的單個半球形粗糙元進(jìn)行了研究。他們對這種孤立式的粗糙元轉(zhuǎn)捩機理進(jìn)行了定性描述：由于粗糙元的存在，邊界層內(nèi)會產(chǎn)生三維的流動分離，進(jìn)而會形成剪切層和渦結(jié)構(gòu)；在粗糙元下游，不穩(wěn)定的反轉(zhuǎn)渦結(jié)構(gòu)會不斷沖擊剪切層從而形成發(fā)卡渦。Muppidi和 Mahesh[17]對來流馬赫數(shù)2.9的平板上的分布式粗糙元進(jìn)行了直接數(shù)值模擬(DNS)，他們認(rèn)為尾跡渦上洗作用會促使流向渦結(jié)構(gòu)與剪切層不斷相互作用，從而導(dǎo)致邊界層轉(zhuǎn)捩。

Tullio等[18]綜合應(yīng)用DNS、三維穩(wěn)定性方程(Three Dimensional Parabolized Stability Equation, PSE-3D)和二維全局穩(wěn)定性(BiGlobal Stability)方法研究了孤立的長方體粗糙元導(dǎo)致的超聲速平板邊界層轉(zhuǎn)捩問題，來流馬赫數(shù)為2.5。他們認(rèn)為粗糙元誘導(dǎo)轉(zhuǎn)捩是整個三維剪切層不穩(wěn)定性的結(jié)果而不是傳統(tǒng)的K-H不穩(wěn)定性。他們發(fā)現(xiàn)粗糙元改變了基本流的分布，從而導(dǎo)致流動穩(wěn)定性的極大改變，而其中最重要的因素就是尾跡中的反轉(zhuǎn)渦對將物面附近的低能流體抬升。

雖然對強制轉(zhuǎn)捩機理存在著多種理論解釋，但多數(shù)人認(rèn)為，粗糙元對基本流的改變以及其自身提供的擾動在轉(zhuǎn)捩中占據(jù)主導(dǎo)作用。在眾多影響強制轉(zhuǎn)捩效果的因素中，粗糙元的高度是最重要的參數(shù)之一，對其研究也較為透徹；粗糙元的幾何外形對轉(zhuǎn)捩效果也有較大的影響，目前研究相對較少。

評價粗糙元誘導(dǎo)的轉(zhuǎn)捩，必須綜合考慮其轉(zhuǎn)捩效果及其帶來的氣動熱、氣動力增量。因此本文著重針對粗糙元的幾何外形對轉(zhuǎn)捩的影響進(jìn)行系統(tǒng)研究，給出各參數(shù)的定量影響，為工程應(yīng)用提供指導(dǎo)。具體來說，是利用實驗室開發(fā)的三維Navier-Stokes方程求解程序UNITs (Unsteady Navier-Stokes solver)，采用高精度的數(shù)值方法，對粗糙元誘導(dǎo)的高超聲速邊界層轉(zhuǎn)捩進(jìn)行了直接數(shù)值模擬研究。

1　數(shù)值方法

本文采用基于有限體積方法的DNS求解Navier-Stokes方程。對流項采用任玉新等[19-20]發(fā)展的色散最小耗散可調(diào)(MDCD)格式進(jìn)行重構(gòu)，該格式可以對色散和耗散特性分別進(jìn)行控制，在保證色散最小的前提下，優(yōu)化耗散特性。

MDCD格式的基本思路為，采用P個模板點進(jìn)行數(shù)值重構(gòu)時，最多可得到P階精度結(jié)果。如果降低精度為P-2階，則可得到2個自由參數(shù)。這2個自由參數(shù)分別對應(yīng)了格式的色散特性和耗散特性，即可對該格式的色散和耗散分別進(jìn)行優(yōu)化控制。MDCD為線性重構(gòu)格式，為了能夠計算高超聲速流動，將其與加權(quán)緊致非線性(WENO)格式進(jìn)行耦合，在流動間斷處采用WENO格式，而在光滑流場采用MDCD格式，即構(gòu)成了MDCD-WENO格式。具體過程可見文獻(xiàn)[19-20]。

具體來說，本文采用了4階MDCD-WENO格式。其最小耗散可等于零，即與6階中心格式相同；最大耗散則要大于傳統(tǒng)WENO格式。在耗散特性不變的基礎(chǔ)上，其色散特性同樣進(jìn)行了優(yōu)化，達(dá)到了該格式的理論最小值。

MDCD-WENO格式重構(gòu)已經(jīng)成功應(yīng)用于直接數(shù)值模擬粗糙元觸發(fā)邊界層轉(zhuǎn)捩的問題，計算結(jié)果和試驗吻合較好[21-22]。

另外，采用隱式上下分解對稱高斯賽德爾偽時間推進(jìn) (Lower-Upper Symmetric-Gauss-Seidel pseudo Time Sub-iteration， LU-SGS-TST)方法進(jìn)行時間推進(jìn)。全局統(tǒng)一時間步長(無量綱時間t=1×10-5)以獲得非定常流動特征，通過消息傳遞界面 (Message Passing Interface， MPI)實現(xiàn)并行計算。

2　模型驗證

計算模型是一個高h(yuǎn)=10.2 mm，直徑D=5.97 mm的圓柱。該試驗由Wheaton和Schneider[23]在普渡大學(xué)的靜音風(fēng)洞中進(jìn)行，風(fēng)洞直徑Dt=240 mm。當(dāng)流場中沒有粗糙元時，該處的邊界層厚度約為8.3 mm。數(shù)值計算中，自由來流馬赫數(shù)Ma=6，雷諾數(shù)Re=6.657×107m-1，來流溫度T∞=53 K， 壁面溫度Tw=300 K。三維計算網(wǎng)格總數(shù)約為1億，粗糙元尾跡區(qū)的網(wǎng)格在流向和展向幾乎各向同性，尺度為0.4 mm，網(wǎng)格距壁面的最小距離為20 μm。計算域為粗糙元上游15D處至下游66D處。

計算的流場相干結(jié)構(gòu)如圖1所示，用Q的等值面表示，并以無量綱速度u/U∞著色。本文的計算結(jié)果與Bartkowicz等[15,24]的結(jié)果相近。

圖1　流場中的相干結(jié)構(gòu)(以速度著色的Q等值面表示)Fig.1　Coherent structures showed by Q contour surface (colored by velocity)

風(fēng)洞試驗中在粗糙元上游1.5D處的風(fēng)洞壁面布置了脈動壓力傳感器，觀測到粗糙元的存在會造成主頻率約為21 kHz的壓力脈動，而在數(shù)值計算中獲得的主頻為19.23 kHz，如圖 2所示。Bartkowicz等[15,24]對同樣算例采用2.7億網(wǎng)格進(jìn)行了數(shù)值模擬，其計算的頻率約為18 kHz。圖中縱坐標(biāo)為歸一化壓力脈動F。

圖2　試驗和計算的壓力脈動頻譜Fig.2　Frequency spectrum of pressure fluctuation in test and calculation

3　研究結(jié)果及討論

3.1不同形狀粗糙元誘導(dǎo)轉(zhuǎn)捩的效果

計算模型為3種不同形狀的粗糙元，基準(zhǔn)流動為平板邊界層流動。粗糙元的形狀如圖 3所示，對于斜坡形粗糙元，其底邊長d1=3.2 mm，高度h與當(dāng)?shù)氐倪吔鐚雍穸然鞠嗤?，?.8 mm，斜坡角θ=100；鉆石形粗糙元的對角線長度d2=3.2 mm，圓柱形粗糙元的直徑d3=3.2 mm，其高度與斜坡高度相同，同為0.8 mm。粗糙元距平板前緣Lu=60 mm。在數(shù)值計算中，將粗糙元下游長度Ld設(shè)為100 mm，平板寬度Wp設(shè)為34.1 mm，如圖3所示。

圖3　計算模型Fig.3　Calculation model

Tirtey等[10]在馮·卡門研究所(Von Karman Institute, VKI)的高超聲速風(fēng)洞中進(jìn)行了該試驗。根據(jù)試驗數(shù)據(jù)，來流馬赫數(shù)為6，雷諾數(shù)為2.6×107m-1，來流溫度約為70 K，物面設(shè)為等溫壁面，溫度為300 K。

本文研究的斜坡形、鉆石形和圓柱形粗糙元的高度和寬度都相同，但是轉(zhuǎn)捩效果卻有較大的差別。圖 4給出了距離平板0.4 mm(0.5h)處的溫度T分布云圖，在粗糙元尾跡區(qū)，可以看到極強的剪切層和流向渦結(jié)構(gòu)，隨著流動向下游發(fā)展，流向渦逐漸扭曲形成規(guī)則的渦結(jié)構(gòu)，在更下游接近出口處，流動結(jié)構(gòu)變得更加復(fù)雜無序，流動轉(zhuǎn)捩為湍流。

對于斜坡形粗糙元，其馬蹄渦很弱，尾跡集中在較窄的范圍內(nèi)；對于鉆石形和圓柱形粗糙元，除了在其正后方的流動區(qū)域內(nèi)形成尾跡渦，在展向較遠(yuǎn)距離處，會形成極強的馬蹄渦。鉆石形粗糙元的馬蹄渦比較穩(wěn)定，在計算域內(nèi)沒有失穩(wěn)；而圓柱粗糙元的馬蹄渦更加不穩(wěn)定，在其下游20 mm處就形成了規(guī)則的渦結(jié)構(gòu)，在下游60 mm處與尾跡渦糾纏在一起，流動轉(zhuǎn)捩為湍流。

圖4　距平板0.4 mm處的溫度分布云圖Fig.4　Temperature distribution contour on y=0.4 mm plane

圖5給出了粗糙元附近的壓力脈動頻譜。壓力樣本點分別位于粗糙元上下游分離泡的起始位置附近。對于斜坡形單元，其上游壓力樣本點位置為(-4.5,0,0) mm，下游位置為(2.3,0,0) mm；對于鉆石形和圓柱形粗糙元，上下游壓力樣本點位置分別為(-2.4,0,0) mm和(3.4,0,0) mm。斜坡形粗糙元的不穩(wěn)定分離泡脈動主頻為55.27 kHz，同時還存在著明顯的諧頻(110.54, 165.82，221.09 kHz等)，而鉆石形和圓柱形粗糙元的脈動主頻則分別為30.19 kHz和42.58 kHz，沒有明顯的諧頻出現(xiàn)。 圓柱形粗糙元的脈動能量最高，鉆石形粗糙元次之，斜坡形粗糙元的脈動能量最低。

圖5　粗糙元上游和下游位置處的壓力脈動頻譜Fig.5　Frequency spectrum of pressure fluctuationbefore and after roughness elements

圖6　不同展向位置截面的最大湍動能分布Fig.6　Maximum TKE distribution on different streamwise cross-section planes

在對稱面位置(z/d=0)處，在0 mm

在x>40 mm之后，擾動增長速度逐漸下降，之后TKEmax達(dá)到峰值，意味著湍動能達(dá)到飽和狀態(tài)。在此區(qū)間范圍，擾動經(jīng)歷了非線性失穩(wěn)和“breakdown”過程。通常擾動在“breakdown”之后流動迅速轉(zhuǎn)捩為湍流，因此本文以擾動達(dá)到峰值位置為轉(zhuǎn)捩位置，下文將對此進(jìn)行比較。

斜坡形粗糙元在x≈80 mm處擾動達(dá)到峰值，即發(fā)生轉(zhuǎn)捩，而鉆石形和圓柱形的轉(zhuǎn)捩位置為x≈70 mm。

在z/d=0.5和1.0位置處，圓柱形粗糙元的TKEmax均能到達(dá)峰值，說明在這些位置其尾跡都發(fā)生了轉(zhuǎn)捩；而斜坡形和鉆石形粗糙元在這些位置處沒有發(fā)生轉(zhuǎn)捩。這和圖 4中的瞬時溫度分布情況一致。

圖7給出了平板上的無量綱熱流斯坦頓數(shù)St分布，可以看到斜坡形的尾跡最窄，而圓柱形粗糙元的尾跡最寬，與之前的分析一致。

圖7　平板的無量綱熱流分布Fig.7　Normalized heat flux distribution on wall surface

粗糙元附近的無量綱熱流分布如圖8所示，圓柱粗糙元的頭激波強度最大，波后溫度最高(見圖4)，因此壁面熱流最大；而斜坡形粗糙元的頭激波最弱，其相應(yīng)的壁面熱流也最低；鉆石形粗糙元的熱流居中。另外，由于圓柱形和鉆石形粗糙元的鈍度較大，其上游會形成馬蹄渦，在壁面上則表現(xiàn)為粗糙元周圍的高熱流條帶；當(dāng)前的斜坡形粗糙元則沒有形成明顯的馬蹄渦結(jié)構(gòu)。

表1給出了不同粗糙元引起的最大熱流Hmax及其位置分布。對于斜坡形粗糙元，其最大熱流在粗糙元的頂點位置(即粗糙元最高處的中點位置)，其熱流最大值約為該處層流熱流Hlam的16.7倍；鉆石形粗糙元的最大熱流在迎風(fēng)頂點處，約為層流的133.7倍；圓柱形粗糙元的最大熱流位置與鉆石形粗糙元相同，其約為層流的173.8倍。

圖8　粗糙元附近的無量綱熱流分布Fig.8　Normalized heat flux distribution near roughness elements

表1　粗糙元引起的最大熱流及位置

Table 1Maximum heat flux caused by roughness element and its location

TypeHmaxHmax/HlamLocationRamp4.10616.691VertexDiamond32.888133.691WindwardvertexCylinder42.753173.793Windwardvertex

計算結(jié)果表明，粗糙元尾跡區(qū)內(nèi)的三維剪切層失穩(wěn)是流動轉(zhuǎn)捩的主要原因。對于不同形狀的粗糙元，斜坡形的湍流尾跡最窄，圓柱形的湍流尾跡最寬；鉆石形和圓柱形粗糙元下游的轉(zhuǎn)捩位置更靠前，斜坡形的轉(zhuǎn)捩位置稍微靠后。

另一方面，鉆石形和圓柱形粗糙元引起的熱流峰值遠(yuǎn)高于斜坡形粗糙元的熱流峰值；鉆石形粗糙元的熱流峰值稍小于圓柱形。綜合轉(zhuǎn)捩效果與引起的熱流峰值來看，斜坡形粗糙元可以作為較好的轉(zhuǎn)捩觸發(fā)裝置。

3.2斜坡形粗糙元對誘導(dǎo)轉(zhuǎn)捩的參數(shù)化

3.2.1斜坡角度

斜坡形粗糙元的寬度和高度與前文的模型相同，而斜坡角度則分別取為10°、20°和30°。保持粗糙元的安裝位置和高度不變，相當(dāng)于粗糙元長度變小，分別為4.54、2.20和1.39 mm。

3種不同斜坡角度的斜坡形粗糙元尾跡渦結(jié)構(gòu)如圖9所示，以無量綱Q=40 000的等值面表示，并以無量綱流向速度著色。不同角度的粗糙元尾跡渦結(jié)構(gòu)相似，其渦管均失穩(wěn)形成“發(fā)卡渦”結(jié)構(gòu)，標(biāo)志著流動開始轉(zhuǎn)捩為湍流。在20° 斜坡的尾跡區(qū)，相干結(jié)構(gòu)的尺度要明顯小于10° 斜坡，并且尾跡的寬度更大，渦管開始變形的位置更靠近上游。對于30° 斜坡而言，其尾跡區(qū)渦結(jié)構(gòu)與20° 斜坡幾乎相同，渦管失穩(wěn)的位置更靠前。

在斜坡形粗糙元1/2高度(y=0.4 mm)截面內(nèi)的溫度分布能更清晰地顯示尾跡寬度及其失穩(wěn)過程，如圖10所示。粗糙元尾跡區(qū)內(nèi)兩側(cè)出現(xiàn)穩(wěn)定的低溫條帶結(jié)構(gòu)，隨著流動向下游發(fā)展，條帶結(jié)構(gòu)逐漸失穩(wěn)形成規(guī)則扭曲，隨后破碎成更小的流動結(jié)構(gòu)，進(jìn)入無序狀態(tài)，流動轉(zhuǎn)捩為湍流。對于10° 斜坡，條帶結(jié)構(gòu)扭曲的位置約為x=31.8 mm，20° 和30° 斜坡分別約為21.1 mm和16.4 mm。到達(dá)計算域出口(x=100 mm)位置，10° 斜坡的尾跡區(qū)寬度約為8.1 mm，而其他兩個斜坡的尾跡寬度幾乎相同，約為10.9 mm。

圖9　不同斜坡角度的斜坡形粗糙元尾跡渦結(jié)構(gòu)Fig.9　Wake vortex structures of ramp roughnesselements with different ramp angles

在斜坡形粗糙元對稱面內(nèi)，不同流向位置的湍動能最大值如圖11所示，其可以表征尾跡擾動的發(fā)展過程。3個粗糙元的擾動發(fā)展過程類似，經(jīng)過線性和非線性發(fā)展階段之后，擾動進(jìn)入飽和狀態(tài)，標(biāo)志著流動發(fā)生轉(zhuǎn)捩。

圖10　不同斜坡角度的斜坡形粗糙元y=0.4 mm截面溫度分布Fig.10　Temperature distribution on y=0.4 mm cross-section of ramp roughness elements with different ramp angles

圖11　不同斜坡角度的斜坡形粗糙元對稱面的擾動沿流向發(fā)展Fig.11　Disturbance streamwise evolution on symmetry plane of ramp roughness elements with different ramp angles

在流向位置x=63.0 mm之前，10° 斜坡的擾動幅值最小，20° 斜坡居中，30° 斜坡擾動最大，可見隨著斜坡角度增大，其逆壓梯度增大，形成的擾動也隨之增大。10° 斜坡的擾動線性增長區(qū)約為22.0 mm

在流向位置x=79.5 mm之后，3個斜坡造成的擾動幅值不再增長，幾乎趨于統(tǒng)一水平，約為0.013。通過判斷湍動能最大值的峰值位置，可以看出3個斜坡尾跡對稱面處的轉(zhuǎn)捩位置分別為x=79.0、54.0和52.5 mm。

3.2.2斜坡間距

斜坡形粗糙元的寬度和高度與3.1節(jié)的模型相同，斜坡角度為10°。選取兩個粗糙元模型，粗糙元對稱面之間的距離稱為間距，記作Δ，取Δ=1.0d1、1.5d1和2.0d1這3種間距情況進(jìn)行計算，如圖12所示。展向設(shè)為周期邊界條件。

圖12　斜坡形粗糙元的間距示意圖Fig.12　Schematic of interval space between ramproughness elements

圖13　不同間距的斜坡形單元對稱面的擾動發(fā)展Fig.13　Disturbance evolution on symmetry plane of ramp roughness elements with different interval space

圖13給出了尾跡中的擾動發(fā)展過程。在粗糙元正后方平面內(nèi)(z=0 mm)，分布式粗糙元的尾跡內(nèi)擾動幅值更大，隨著間距越小擾動幅值增加，但整體而言擾動幅值相差不大。對于Δ/d1=1.0的工況，其尾跡擾動在x=68.5 mm即到達(dá)峰值，轉(zhuǎn)捩位置更靠前；而對于Δ/d1=1.5和2.0的工況，峰值位置與單個粗糙元情況幾乎相同，約為75.5 mm。分布式粗糙單元尾跡擾動的線性增長率約為0.074，小于單個粗糙元的增長率0.099。

在z=0.8 mm平面內(nèi)，分布式粗糙元尾跡擾動幅值大于單個粗糙元情況，Δ/d1=1.0工況的擾動幅值與其他3種情況差距明顯，在擾動線性發(fā)展階段其幅值約高出其他3種工況1～2個數(shù)量級，并且其擾動飽和位置更靠前，約為x=55.5 mm；Δ/d1=1.5和2.0工況的擾動幅值相近，轉(zhuǎn)捩位置約為x=61.5 mm和66.0 mm；單個粗糙元的轉(zhuǎn)捩位置最為靠后，約為x=77.0 mm。

z=1.6 mm截面內(nèi)的擾動發(fā)展情況與z=0.8 mm 截面類似，Δ/d1=1.0工況擾動的幅值和飽和位置均領(lǐng)先于其他工況，轉(zhuǎn)捩位置為x=68.0 mm，Δ/d1=1.5和 2.0工況的轉(zhuǎn)捩位置為74.5 mm和82.5 mm，單個粗糙元的轉(zhuǎn)捩位置約為x=91.5 mm。

4　結(jié)　論

1) 粗糙元尾跡區(qū)內(nèi)的三維剪切層失穩(wěn)是流動轉(zhuǎn)捩的主要原因。對于不同形狀的粗糙元，斜坡形的湍流尾跡最窄，圓柱形的湍流尾跡最寬；鉆石形和圓柱形粗糙元下游的轉(zhuǎn)捩位置更靠前，斜坡形的轉(zhuǎn)捩位置稍微靠后。

2) 鉆石形和圓柱形粗糙元引起的熱流峰值遠(yuǎn)高于斜坡形粗糙元的熱流峰值；鉆石形粗糙元的熱流峰值稍小于圓柱形。綜合轉(zhuǎn)捩效果與引起的熱流峰值來看，斜坡形粗糙元可以作為較好的轉(zhuǎn)捩觸發(fā)裝置。

3) 隨著斜坡粗糙元角度的增大，其鈍度隨之增大，上游逆壓梯度增強，尾跡區(qū)寬度增大，尾跡擴張角增大，但是20° 斜坡與30° 斜坡幾乎沒有差距；尾跡區(qū)擾動幅值增大，線性增長率減小，轉(zhuǎn)捩位置提前，流動轉(zhuǎn)捩位置相同

4) 隨著粗糙元之間的間距減小，其造成的擾動幅值增大；擾動線性增長率幾乎相同，但小于單個粗糙元工況；尾跡區(qū)轉(zhuǎn)捩距離提前。綜合比較，粗糙元的安裝間距Δ/d1=1.0時，觸發(fā)轉(zhuǎn)捩效果最優(yōu)。

致　謝

感謝清華大學(xué)信息科學(xué)與技術(shù)國家實驗室提供計算資源。感謝中航工業(yè)產(chǎn)學(xué)研工程項目“航空CFD共用軟件體系中的若干關(guān)鍵技術(shù)研究”的資助。

[1]WHITEHEAD A, JR. NASP aerodynamics: AIAA-1989-5013[R]. Reston: AIAA, 1989.

[2]VAN DRIEST E R, BLUMER C B. Boundary-layer at supersonic speeds-three dimensional roughness effects (spheres):SID 61-275[R]. Wichita: North American Aviation, Inc., 1961.

[3]趙慧勇， 周瑜， 倪鴻禮, 等. 高超聲速進(jìn)氣道邊界層強制轉(zhuǎn)捩試驗[J]. 實驗流體力學(xué)，2012, 26(1): 1-6.

ZHAO H Y, ZHOU Y, NI H L, et al. Test of forced boundary-layer transition on hypersonic inlet [J]. Journal of Experiments in Fluid Mechanics, 2012, 26(1): 1-6 (in Chinese).

[4]SCHNEIDER S P. Effects of roughness on hypersonic boundary-layer transition[J]. Journal of Spacecraft and Rockets, 2008, 45(2)： 193-209.

[5]CASPER K M, WHEATON B M, JOHNSON H B, et al. Effect of freestream noise on roughness-induced transition at Mach 6: AIAA-2008-4291[R]. Reston: AIAA, 2008.

[6]VAN DRIEST E R, MCCAULEY W. The effect of controlled three-dimensional roughness on boundary layer transition at supersonic speeds[J]. Journal of Aerospace Science, 1960, 27(4): 261-271

[7]HICKS R M, HARPER W R, JR. A comparison of spherical and triangular boundary-layer tips on a flat plate at supersonic speeds: NASA, TM-X-2146[R]. Washington, D.C.: NASA, 1970.

[8]BERRY S A, AUSLENDER A H. Hypersonic boundary-layer trip development for hyper-X[J]. Journal of Spacecraft and Rockets, 2001, 38(6): 853-864.

[9]WHITEHEAD A H, JR. Flowfield and drag characteristics of several boundary-layer tripping elements in hypersonic flow: NASATND-5454[R]. Washington, D.C.: NASA, 1969.

[10]TIRTEY S C, CHAZOT O, WALPOT L. Characterization of hypersonic roughness-induced boundary-layer transition [J]. Experiments in Fluids, 2011, 50(2): 407-418.

[11]DANEHY P M, IVEY C B, INMAN J A, et al. High-speed PLIF imaging of hypersonic transition over discrete cylindrical roughness: AIAA-2010-0703[R]. Reston: AIAA, 2010.

[12]DANEHY P M, BATHEL B F, IVEY C B, et al. NO PLIF study of hypersonic transition over a discrete hemispherical roughness element: AIAA-2009-0394[R]. Reston: AIAA, 2009.

[13]周玲， 閻超， 孔維萱. 高超聲速飛行器前體邊界層強制轉(zhuǎn)捩數(shù)值模擬[J]. 航空學(xué)報， 2014， 35(6): 1487-1495.

ZHOU L, YAN C, KONG W X. Numerical simulation of forced boundary layer transition on hypersonic vehicle forebody[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(6): 1487-1495 (in Chinese).

[14]BERNARDINI M, PIROZZOLI S, ORLANDI P. Compressibility effects on roughness-induced boundary layer transition[J]. International Journal of Heat and Fluids Flow, 2012, 35(35): 45-51.

[15]BARTKOWICZ M D, SUBBAREDDY P K, CANDLER G V. Numerical simulations of roughness induced instability in the purdue mach 6 wind tunnel: AIAA-2010-4723[R]. Reston: AIAA, 2010.

[16]IYER P S, MAHESH K. High-speed boundary layer transition induced by a discrete roughness element[J]. Journal of Fluid Mechanics, 2013, 729: 524-562.

[17]MUPPIDI S, MAHESH K. Direct numerical simulations of roughness-induced transition in supersonic boundary layers[J]. Journal of Fluid Mechanics, 2012, 693: 28-56.

[18]TULLIO N D, PAREDES P, SANDHAM N D, et al. Laminar-turbulent transition induced by a discrete roughness element in a supersonic boundary layer[J]. Journal of Fluid Mechanics, 2013, 735: 613-646.

[19]SUN Z S, REN Y X, LARRICQ C, et al. A class of finite difference schemes with low dispersion and controllable dissipation for DNS of compressible turbulence[J]. Journal of Computational Physics, 2011， 230(12)： 4616-4635.

[20]WANG Q J, REN Y X, SUN Z S. High resolution finite volume scheme based on minimized dispersion and controllable dissipation reconstruction[J]. Science China Physics, Mechanics & Astronomy, 2013, 56(2): 423-431.

[21]DUAN Z W, XIAO Z X, FU S. Direct numerical simulation of hypersonic transition induced by an isolated cylindrical roughness element[J]. Science China Physics, Mechanics & Astronomy, 2014, 57(12): 2330-2345.

[22]DUAN Z W, XIAO Z X, FU S. Direct numerical simulation of hypersonic transition induced by a cylindrical roughness element: AIAA-2013-3112[R]. Reston: AIAA, 2013.

[23]WHEATON B M, SCHNEIDER S P. Roughness-induced instability in a hypersonic laminar boundary layer[J].AIAA Journal, 2012, 50(6): 1245-1256.

[24]WHEATON B M, BARTKOWICZ M D, SUBBAREDDY P K, et al. Roughness-induced instabilities at Mach 6: A combined numerical and experimental study: AIAA-2011-3248[R]. Reston: AIAA, 2011.

段志偉男， 博士， 助理研究員。主要研究方向： 空氣動力學(xué)、 高精度湍流預(yù)測和邊界層轉(zhuǎn)捩。

Tel.: 010-62795411

E-mail: dzw15@tsinghua.edu.cn

肖志祥男，博士， 副研究員， 博士生導(dǎo)師。主要研究方向： RANS-LES混合方法、 高精度湍流預(yù)測、 流動轉(zhuǎn)捩和計算氣動聲學(xué)。

Tel.: 010-62797060

E-mail: xiaotigerzhx@tsinghua.edu.cn

Roughness element induced hypersonic boundary layer transition

DUAN Zhiwei, XIAO Zhixiang*

School of Aerospace Engineering, Tsinghua University, Beijing100084, China

Hypersonic boundary layer transition from laminar to turbulent induced by different isolated roughness elements is investigated using direct numerical simulation based on finite volume formulation. To explore the transition mechanism through resolving the small flow structure, and capture shock wave in hypersonic flow, the high order minimized dispersion and controllable dissipation (MDCD) scheme is used to reconstruct the convection terms of Navier-Stokes equations. The numerical results agree well with experimental data. The numerical method adopted in this article is able to resolve small flow structures and their break-up and instability procedure. And it shows that the transition is dominated by the instability of the three-dimensional shear layer on top of the roughness elements. The effect of roughness element shape and geometrical parameters on transition mechanisms, and the transition results of multiple roughness elements are quantitatively studied.

hypersonic； roughness element； boundary layer transition； parameter effects; direct numerical simulation

2016-01-18； Revised： 2016-02-17； Accepted： 2016-04-26； Published online： 2016-04-2915：37

s： National Natural Science Foundation of China (11372159); China Postdoctoral Science Foundation (2015M571029)； National Key Research and Development Project (2016YFA0401200)

. Tel.： 010-62797060E-mail: xiaotigerzhx@tsinghua.edu.cn

2016-01-18； 退修日期： 2016-02-17； 錄用日期： 2016-04-26;

時間： 2016-04-2915:37

www.cnki.net/kcms/detail/11.1929.V.20160429.1537.008.html

國家自然科學(xué)基金 (11372159)； 中國博士后科學(xué)基金 (2015M571029); 國家重點研發(fā)計劃項目 (2016YFA0401200)

.Tel.： 010-62797060E-mail: xiaotigerzhx@tsinghua.edu.cn

10.7527/S1000-6893.2016.0130

V211.3

A

1000-6893(2016)08-2454-10

引用格式： 段志偉， 肖志祥. 粗糙元誘導(dǎo)的高超聲速邊界層轉(zhuǎn)捩[J]． 航空學(xué)報, 2016, 37(8): 2454-2463. DUAN Z W, XIAO Z X. Roughness element induced hypersonic boundary layer transition[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(8): 2454-2463.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

URL： www.cnki.net/kcms/detail/11.1929.V.20160429.1537.008.html