段志偉， 肖志祥

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

?

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

段志偉， 肖志祥*

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

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

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

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

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

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

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

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

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

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

對(duì)于大尺寸的粗糙元，Iyer和Mahesh[16]對(duì)平板上的單個(gè)半球形粗糙元進(jìn)行了研究。他們對(duì)這種孤立式的粗糙元轉(zhuǎn)捩機(jī)理進(jìn)行了定性描述：由于粗糙元的存在，邊界層內(nèi)會(huì)產(chǎn)生三維的流動(dòng)分離，進(jìn)而會(huì)形成剪切層和渦結(jié)構(gòu)；在粗糙元下游，不穩(wěn)定的反轉(zhuǎn)渦結(jié)構(gòu)會(huì)不斷沖擊剪切層從而形成發(fā)卡渦。Muppidi和 Mahesh[17]對(duì)來(lái)流馬赫數(shù)2.9的平板上的分布式粗糙元進(jìn)行了直接數(shù)值模擬(DNS)，他們認(rèn)為尾跡渦上洗作用會(huì)促使流向渦結(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)方法研究了孤立的長(zhǎng)方體粗糙元導(dǎo)致的超聲速平板邊界層轉(zhuǎn)捩問(wèn)題，來(lái)流馬赫數(shù)為2.5。他們認(rèn)為粗糙元誘導(dǎo)轉(zhuǎn)捩是整個(gè)三維剪切層不穩(wěn)定性的結(jié)果而不是傳統(tǒng)的K-H不穩(wěn)定性。他們發(fā)現(xiàn)粗糙元改變了基本流的分布，從而導(dǎo)致流動(dòng)穩(wěn)定性的極大改變，而其中最重要的因素就是尾跡中的反轉(zhuǎn)渦對(duì)將物面附近的低能流體抬升。

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

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

1　數(shù)值方法

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

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

具體來(lái)說(shuō)，本文采用了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)捩的問(wèn)題，計(jì)算結(jié)果和試驗(yàn)吻合較好[21-22]。

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

2　模型驗(yàn)證

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

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

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

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

圖2　試驗(yàn)和計(jì)算的壓力脈動(dòng)頻譜Fig.2　Frequency spectrum of pressure fluctuation in test and calculation

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

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

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

圖3　計(jì)算模型Fig.3　Calculation model

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

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

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

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

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

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

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

在對(duì)稱(chēng)面位置(z/d=0)處，在0 mm

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

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

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

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

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

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

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

圖8　粗糙元附近的無(wú)量綱熱流分布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

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

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

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

3.2.1斜坡角度

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

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

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

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

在斜坡形粗糙元對(duì)稱(chēng)面內(nèi)，不同流向位置的湍動(dòng)能最大值如圖11所示，其可以表征尾跡擾動(dòng)的發(fā)展過(guò)程。3個(gè)粗糙元的擾動(dòng)發(fā)展過(guò)程類(lèi)似，經(jīng)過(guò)線(xiàn)性和非線(xiàn)性發(fā)展階段之后，擾動(dòng)進(jìn)入飽和狀態(tài)，標(biāo)志著流動(dòng)發(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　不同斜坡角度的斜坡形粗糙元對(duì)稱(chēng)面的擾動(dòng)沿流向發(fā)展Fig.11　Disturbance streamwise evolution on symmetry plane of ramp roughness elements with different ramp angles

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

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

3.2.2斜坡間距

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

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

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

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

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

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

4　結(jié)　論

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

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

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

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

致　謝

感謝清華大學(xué)信息科學(xué)與技術(shù)國(guó)家實(shí)驗(yàn)室提供計(jì)算資源。感謝中航工業(yè)產(chǎn)學(xué)研工程項(xiàng)目“航空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)氣道邊界層強(qiáng)制轉(zhuǎn)捩試驗(yàn)[J]. 實(shí)驗(yàn)流體力學(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]周玲， 閻超， 孔維萱. 高超聲速飛行器前體邊界層強(qiáng)制轉(zhuǎn)捩數(shù)值模擬[J]. 航空學(xué)報(bào)， 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.

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

Tel.: 010-62795411

E-mail: dzw15@tsinghua.edu.cn

肖志祥男，博士， 副研究員， 博士生導(dǎo)師。主要研究方向： RANS-LES混合方法、 高精度湍流預(yù)測(cè)、 流動(dòng)轉(zhuǎn)捩和計(jì)算氣動(dòng)聲學(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;

時(shí)間： 2016-04-2915:37

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

國(guó)家自然科學(xué)基金 (11372159)； 中國(guó)博士后科學(xué)基金 (2015M571029); 國(guó)家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目 (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é)報(bào), 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