高超聲速圓錐邊界層轉(zhuǎn)捩反轉(zhuǎn)數(shù)值研究

雷娟棉，曹家偉

（北京理工大學(xué) 宇航學(xué)院，北京 100081）

高超聲速邊界層從層流到湍流的轉(zhuǎn)捩問題是空氣動力學(xué)領(lǐng)域的難點(diǎn)和熱點(diǎn)問題之一，邊界層轉(zhuǎn)捩后飛行器表面的摩擦阻力系數(shù)、熱流會急劇上升，對整個飛行器的氣動力、氣動熱特性產(chǎn)生重要影響，進(jìn)而影響著飛行器的飛行效率和飛行安全[1]. 因此，關(guān)于邊界層轉(zhuǎn)捩問題的研究一直是空氣動力學(xué)領(lǐng)域的研究熱點(diǎn). 在相關(guān)的風(fēng)洞實(shí)驗(yàn)、數(shù)值模擬計(jì)算、理論分析中，已探明主要原因的轉(zhuǎn)捩現(xiàn)象和轉(zhuǎn)捩規(guī)律被認(rèn)為是一般性結(jié)論，而與一般性結(jié)論中的轉(zhuǎn)捩趨勢相反的現(xiàn)象被稱為轉(zhuǎn)捩反轉(zhuǎn)[2]. 圓錐體外形作為高超聲速飛行器的典型常用外形，被廣泛應(yīng)用于邊界層轉(zhuǎn)捩問題研究. 目前，已有部分實(shí)驗(yàn)研究結(jié)果表明，在某些情況下圓錐表面邊界層呈現(xiàn)轉(zhuǎn)捩反轉(zhuǎn)的現(xiàn)象，但關(guān)于邊界層轉(zhuǎn)捩反轉(zhuǎn)的認(rèn)識非常有限，邊界層轉(zhuǎn)捩反轉(zhuǎn)現(xiàn)象仍然是目前轉(zhuǎn)捩研究中尚未探明原因的一個重要問題[3].

對細(xì)長圓錐邊界層轉(zhuǎn)捩行為最早、最全面的研究之一是STETSON[4]在1967 年進(jìn)行的風(fēng)洞實(shí)驗(yàn)研究，實(shí)驗(yàn)研究結(jié)果表明當(dāng)圓錐攻角從0°開始增大時，其背風(fēng)面邊界層轉(zhuǎn)捩位置向前移動，而迎風(fēng)面邊界層轉(zhuǎn)捩位置向后移動，在有攻角條件下圓錐表面的轉(zhuǎn)捩陣面呈現(xiàn)為背風(fēng)面轉(zhuǎn)捩位置靠前、迎風(fēng)面轉(zhuǎn)捩位置靠后的周向不對稱分布，該現(xiàn)象被認(rèn)為是圓錐邊界層轉(zhuǎn)捩典型的一般性規(guī)律. 在某些情況下，亦有可能出現(xiàn)圓錐迎風(fēng)面與背風(fēng)面轉(zhuǎn)捩位置同時前移的現(xiàn)象，但背風(fēng)面轉(zhuǎn)捩前移的幅度要大于迎風(fēng)面，使轉(zhuǎn)捩陣面沿周向呈現(xiàn)同樣的非對稱特性分布規(guī)律.MUIR 等[5]在1972 年的風(fēng)洞實(shí)驗(yàn)中發(fā)現(xiàn)當(dāng)圓錐頭部鈍度增大到特定尺度時，在高雷諾數(shù)條件下圓錐迎風(fēng)面轉(zhuǎn)捩位置隨攻角增大而向前運(yùn)動，而背風(fēng)面轉(zhuǎn)捩位置先向后運(yùn)動再向前運(yùn)動，整體呈現(xiàn)出隨攻角增大迎風(fēng)面轉(zhuǎn)捩先于背風(fēng)面轉(zhuǎn)捩的邊界層轉(zhuǎn)捩反轉(zhuǎn)現(xiàn)象. 陳久芬等[6]在對7°尖錐的紅外測量轉(zhuǎn)捩實(shí)驗(yàn)中發(fā)現(xiàn)當(dāng)攻角α從8°增大到10°時，迎風(fēng)面轉(zhuǎn)捩位置從推遲變?yōu)樘崆?，在部分攻角條件下發(fā)生了邊界層轉(zhuǎn)捩“逆轉(zhuǎn)”的現(xiàn)象.

早在20 世紀(jì)60 年代，在風(fēng)洞實(shí)驗(yàn)中便發(fā)現(xiàn)了圓錐的鈍度變化引起轉(zhuǎn)捩位置反轉(zhuǎn)現(xiàn)象，即在小范圍內(nèi)增大錐體頭部鈍度，轉(zhuǎn)捩位置向下游移動；然而當(dāng)錐體頭部鈍度增大到某一臨界值時，轉(zhuǎn)捩位置隨鈍度增加由推遲轉(zhuǎn)變?yōu)樘崆? 之后，PARADES 等[7]通過在圓錐頭部增加粗糙帶的實(shí)驗(yàn)方法發(fā)現(xiàn)圓錐鈍錐邊界層轉(zhuǎn)捩對頭部粗糙度十分敏感，認(rèn)為轉(zhuǎn)捩提前的原因是隨著頭部鈍度增大圓錐表面對于粗糙度的敏感性增強(qiáng)，從而促進(jìn)了轉(zhuǎn)捩提前. 陳蘇宇等[8]在FD-14A 激波風(fēng)洞中通過實(shí)驗(yàn)在不同馬赫數(shù)下研究了圓錐邊界層轉(zhuǎn)捩鈍度反轉(zhuǎn)的現(xiàn)象，復(fù)現(xiàn)了“Λ”型和“N”型兩類反轉(zhuǎn)，并且通過功率譜密度分析，從熵層影響邊界層結(jié)構(gòu)的角度分析了圓錐鈍度反轉(zhuǎn)的原因. 陳堅(jiān)強(qiáng)等[9]從鈍前緣的感受性以及熵層效應(yīng)兩個角度初步分析了鈍度反轉(zhuǎn)產(chǎn)生的可能原因.

目前對于邊界層轉(zhuǎn)捩反轉(zhuǎn)現(xiàn)象的認(rèn)識主要是基于有限的風(fēng)洞實(shí)驗(yàn)結(jié)果，但由于轉(zhuǎn)捩試驗(yàn)難度大、影響因素多且費(fèi)用高，很難開展大量系統(tǒng)性的轉(zhuǎn)捩反轉(zhuǎn)研究. 而在各種轉(zhuǎn)捩預(yù)測方法中，基于雷諾平均NS (RANS)方程的轉(zhuǎn)捩預(yù)測模型方法具有計(jì)算穩(wěn)定、求解效率高等優(yōu)點(diǎn)，已在高超聲速邊界層轉(zhuǎn)捩預(yù)測中得到了廣泛應(yīng)用[10]. 楊云軍等[11]利 用基于k-ω湍流模型耦合的轉(zhuǎn)捩模型在高超聲速小攻角條件下對圓錐的典型非對稱轉(zhuǎn)捩進(jìn)行了預(yù)測，分析了非稱轉(zhuǎn)捩效應(yīng)下的圓錐氣動力特性. 蔡林峰等[12]通過基于SSTk-ω湍流模型耦合的γ轉(zhuǎn)捩模型，在添加可壓縮性修正后，具備預(yù)測高超聲速流動轉(zhuǎn)捩的能力. 王衛(wèi)星等[13]在針對Ma=4.5 的錐型曲面進(jìn)氣道進(jìn)行數(shù)值模擬時發(fā)現(xiàn)，隨著錐尖鈍度增加，邊界層轉(zhuǎn)捩呈“推遲-提前-推遲”的轉(zhuǎn)捩反轉(zhuǎn)現(xiàn)象.

因此，本文通過求解基于SSTk-ω湍流模型耦合γ轉(zhuǎn)捩模型的RANS 方程，并進(jìn)行可壓縮性修正，對高超聲速圓錐邊界層轉(zhuǎn)捩的攻角反轉(zhuǎn)和鈍度反轉(zhuǎn)問題進(jìn)行數(shù)值模擬研究. 旨在通過數(shù)值模擬方法研究高超聲速圓錐邊界層轉(zhuǎn)捩位置隨攻角、頭部鈍度變化出現(xiàn)的邊界層轉(zhuǎn)捩反轉(zhuǎn)現(xiàn)象，通過不同條件下轉(zhuǎn)捩位置的變化規(guī)律及對流場的影響，對比分析轉(zhuǎn)捩反轉(zhuǎn)與正常一般性轉(zhuǎn)捩規(guī)律的流動特點(diǎn)，為高超聲速飛行器邊界層轉(zhuǎn)捩、特別是轉(zhuǎn)捩反轉(zhuǎn)研究提供支撐.

1 數(shù)值模擬方法

1.1 基本數(shù)值方法

本文數(shù)值研究采用有限體積法求解定常可壓縮RANS 方程組，時間推進(jìn)采用隱式格式，無黏通量采用Roe-FDS 格式，黏性通量采用中心差分格式進(jìn)行離散.

1.2 γ 轉(zhuǎn)捩模型

本文數(shù)值模擬計(jì)算采用的湍流模型為SSTk-ω湍流模型，利用MENTER 等[14]根據(jù)實(shí)驗(yàn)數(shù)據(jù)重新得到的Reθc經(jīng)驗(yàn)關(guān)系式，加入間歇因子輸運(yùn)方程，封閉成一方程γ轉(zhuǎn)捩模型. SSTk-ω[15]模型建立了湍流流體中湍動能k和比耗散率ω的輸運(yùn)方程：

式中Pk和Dk分別為湍動能的生成項(xiàng)和耗散項(xiàng).

SSTk-ω模型不能預(yù)測層流流動的原因是Pk和Dk無法判別邊界層流動狀態(tài). 在SSTk-ω模型的基礎(chǔ)上建立間歇因子γ輸運(yùn)方程并作用于Pk和Dk，從而實(shí)現(xiàn)對流動中從層流到湍流轉(zhuǎn)捩過程的數(shù)值模擬.γ輸運(yùn)方程：

利用MENTER 提出的Reθc與當(dāng)?shù)貕毫μ荻圈甩群彤?dāng)?shù)赝牧鞫萒uL的經(jīng)驗(yàn)關(guān)系式，從而避免了建立Reθc的輸運(yùn)方程：

利用文獻(xiàn)[16]給出的基于來流馬赫數(shù)的可壓縮性修正方法，針對高超聲速可壓縮流轉(zhuǎn)捩問題進(jìn)行可壓縮修正. 基于當(dāng)?shù)伛R赫數(shù)，對臨界動量厚度雷諾數(shù)進(jìn)行修正，具體形式如下：

1.3 計(jì)算外形與邊界條件

本文數(shù)值模擬計(jì)算的外形為圓錐體，圖1 為計(jì)算外形示意圖，其中Rn為圓錐頭部鈍度半徑，θ為圓錐的半錐角，L為圓錐長度，坐標(biāo)系的原點(diǎn)取圓錐的頭部頂點(diǎn). 模型頭部鈍度Rn、半錐角θ、錐長L隨下文中不同研究條件而有所變化. 定義周向角φ為縱向?qū)ΨQ面與周向位置順時針的夾角，背風(fēng)面子午線φ=0°. 計(jì)算域選取圓錐前向距離為0.2 倍錐長,法向距離為2 倍底部半徑的流場范圍. 采用結(jié)構(gòu)化網(wǎng)格對圓錐流場進(jìn)行離散，由于圓錐為對稱外形，為了減小計(jì)算量取半模進(jìn)行計(jì)算. 圖2 為計(jì)算網(wǎng)格示意圖. 壁面采用無滑移等溫壁面條件，來流條件為自由來流遠(yuǎn)場邊界條件，底部為超聲速外推邊界條件.

圖1 圓錐幾何模型示意圖Fig. 1 Schematic diagram of cone geometry model

圖2 圓錐計(jì)算域及網(wǎng)格示意圖Fig. 2 Schematic diagram of cone calculation domain and grid

1.4 數(shù)值方法驗(yàn)證

邊界層轉(zhuǎn)捩前后往往伴隨著表面熱流以及摩阻的突躍變化，因此圓錐表面的熱流或者摩阻分布可以作為轉(zhuǎn)捩位置的一種表征方法. 本節(jié)基于文獻(xiàn)[17]中的高超聲速圓錐邊界層轉(zhuǎn)捩風(fēng)洞實(shí)驗(yàn)的圓錐表面熱流分布數(shù)據(jù)，對本文采用的γ轉(zhuǎn)捩模型對圓錐邊界層轉(zhuǎn)捩預(yù)測能力進(jìn)行驗(yàn)證. 在網(wǎng)格無關(guān)性驗(yàn)證中，分別對網(wǎng)格量為144 萬的A 網(wǎng)格、網(wǎng)格量為230 萬的B 網(wǎng)格、網(wǎng)格量為370 萬的C 網(wǎng)格進(jìn)行了無關(guān)性驗(yàn)證，其中A 網(wǎng)格對于試驗(yàn)熱流值的捕捉相差較大，B 網(wǎng)格與C 網(wǎng)格對試驗(yàn)熱流的吻合度較好. 為了在滿足計(jì)算精度的同時提高計(jì)算效率，下文數(shù)值驗(yàn)證采用200（流向）×120（法向）×100（周向），總量約為230 萬的B 網(wǎng)格進(jìn)行計(jì)算. 計(jì)算外形和計(jì)算條件與文獻(xiàn)[17]中風(fēng)洞試驗(yàn)的條件下相同，圓錐頭部鈍度Rn=2.5 mm，半錐角θ=7°，圓錐長度L=1.2 m，自由流湍流度Tu=0.6%，湍流黏性比μt/μ=10，計(jì)算條件如表1 所示.

表1 算例驗(yàn)證計(jì)算條件Tab. 1 Calculation conditions for verification of calculation examples

本文采用基于γ轉(zhuǎn)捩模型的數(shù)值模擬方法，通過計(jì)算得到了圓錐子午線上的熱流分布曲線，如圖3所示. 為了對比分析，圖3 中也給出了分別采用Langtry和Menter 發(fā)展的γ-Reθt轉(zhuǎn)捩模型、層流模型、湍流模型計(jì)算得到的結(jié)果. 由圖3 可見，采用γ轉(zhuǎn)捩模型計(jì)算得到的圓錐轉(zhuǎn)捩起始位置和熱流峰值的預(yù)測結(jié)果與實(shí)驗(yàn)測得的轉(zhuǎn)捩位置和熱流峰值基本一致，但預(yù)測得到的圓錐下游表面熱流有所偏高.γ-Reθt轉(zhuǎn)捩模型預(yù)測得到的峰值熱流值與實(shí)驗(yàn)值基本吻合，對圓錐轉(zhuǎn)捩起始位置的預(yù)測過于提前，這可能是由于在對經(jīng)驗(yàn)參數(shù)標(biāo)定時未對高速流進(jìn)行修正. 以上結(jié)果表明，基于雷諾平均的N-S 方程，采用SSTk-ω湍流模型耦合γ轉(zhuǎn)捩模型的數(shù)值方法能夠比較準(zhǔn)確的預(yù)測圓錐邊界層轉(zhuǎn)捩位置，適合進(jìn)行對圓錐轉(zhuǎn)捩位置的數(shù)值研究.

圖3 圓錐子午線熱流分布計(jì)算結(jié)果對比Fig. 3 Comparison of calculation results of heat flow distribution on conical meridian

2 轉(zhuǎn)捩反轉(zhuǎn)問題數(shù)值研究與分析

2.1 隨攻角變化的轉(zhuǎn)捩反轉(zhuǎn)

本節(jié)采用基于耦合γ轉(zhuǎn)捩模型的RANS 求解方法，對全長L=452mm，半錐角θ=8°，頭部鈍度半徑分別為Rn=0.063 5 mm、20.32 mm 的尖錐和鈍錐，在馬赫數(shù)Ma=6，攻角α=0°、2°、4°、6°、8°的條件下的流場分別進(jìn)行模擬計(jì)算，給出圓錐邊界層轉(zhuǎn)捩位置隨攻角增大的變化規(guī)律，研究隨攻角增大圓錐邊界層轉(zhuǎn)捩的反轉(zhuǎn)現(xiàn)象. 計(jì)算條件參照文獻(xiàn)[5]，自由流湍流度Tu=0.1%，湍流黏性比μt/μ=10，具體計(jì)算參數(shù)如表2 所示.

表2 攻角反轉(zhuǎn)數(shù)值模擬計(jì)算條件Tab. 2 Numerical simulation calculation conditions of attack angle reversal

圖4 給出了兩種不同鈍度圓錐在迎風(fēng)面和背風(fēng)面子午線上的轉(zhuǎn)捩位置隨攻角變化曲線，圖中也給出了文獻(xiàn)[5]中的風(fēng)洞實(shí)驗(yàn)結(jié)果. 圖4（a）中左側(cè)曲線上向上的箭頭表示圓錐迎風(fēng)面轉(zhuǎn)捩已推遲至計(jì)算外形之外，說明攻角α≥4°時圓錐迎風(fēng)面邊界層表現(xiàn)為層流狀態(tài).

圖4 圓錐迎風(fēng)面和背風(fēng)面對稱面子午線上轉(zhuǎn)捩位置隨攻角變化曲線（Ma∞=6）Fig. 4 The curve of the transition position on the meridian in the symmetry plane of the windward and leeward side of the cone with the angle of attack（Ma∞=6）

由圖4（a）可見，對于尖錐，當(dāng)α≤θ/2 時，隨著攻角增大迎風(fēng)面轉(zhuǎn)捩位置向后移動，背風(fēng)面轉(zhuǎn)捩位置向前移動且前移速隨著攻角增大而減??；當(dāng)α＞θ/2 時，迎風(fēng)面轉(zhuǎn)捩推遲至模型以外，背風(fēng)面轉(zhuǎn)捩位置由向前運(yùn)動轉(zhuǎn)變?yōu)橄蚝筮\(yùn)動. 在攻角不為0 時，頭部鈍度半徑為Rn=0.0635 mm 的圓錐邊界層呈現(xiàn)一般的非對稱轉(zhuǎn)捩特征. 由圖4（b）可見，對于鈍錐，當(dāng)α≤θ/4 時，背風(fēng)面轉(zhuǎn)捩位置后移，當(dāng)α＞θ/4 時，背風(fēng)面轉(zhuǎn)捩位置由后移轉(zhuǎn)變?yōu)橄蚯耙苿?；而迎風(fēng)面轉(zhuǎn)捩位置則隨攻角增大呈單調(diào)前移的趨勢；在攻角不為零時，頭部鈍度半徑為Rn=20.32 mm 的圓錐，迎風(fēng)面的轉(zhuǎn)捩位置都比背風(fēng)面的轉(zhuǎn)捩位置靠前，呈現(xiàn)轉(zhuǎn)捩反轉(zhuǎn)的現(xiàn)象，攻角α=2°時轉(zhuǎn)捩反轉(zhuǎn)現(xiàn)象最顯著. 鈍錐外形在背風(fēng)面的轉(zhuǎn)捩預(yù)測位置變化趨勢與實(shí)驗(yàn)趨勢一致，但轉(zhuǎn)捩預(yù)測位置比實(shí)驗(yàn)值有所提前. 推測可能原因?yàn)棣棉D(zhuǎn)捩模型中的橫流準(zhǔn)則對于背風(fēng)面橫流速度較為敏感，使得轉(zhuǎn)捩提前，但依然較好地預(yù)測出實(shí)驗(yàn)所提到的攻角轉(zhuǎn)捩反轉(zhuǎn)現(xiàn)象.

圖5 和圖6 分別給出了馬赫數(shù)Ma∞=6，攻角α=2°時，圓錐頭部鈍度分別為Rn=0.063 5 mm 的尖圓錐和Rn=20.32 mm 的鈍圓錐表面間歇因子分布圖. 圖5 和圖6 中可見，有攻角情況下尖圓錐呈一般的非對稱轉(zhuǎn)捩特征，而鈍圓錐呈現(xiàn)轉(zhuǎn)捩反轉(zhuǎn)現(xiàn)象，這與上面分析給出的結(jié)論一致. 圖7 給出了尖錐與鈍錐的周向摩阻分布曲線，從圖7（a）中可知，在有攻角的條件下，尖錐表面不同子午線上的轉(zhuǎn)捩位置不同且變化范圍很大，背風(fēng)面中心線（?=0°）上的轉(zhuǎn)捩位置最靠前，迎風(fēng)面中心線（?=180°）上的轉(zhuǎn)捩位置最靠后，沿周向從背風(fēng)面到迎風(fēng)面轉(zhuǎn)捩位置逐漸后移；在周向角?從0°到90°的背風(fēng)區(qū)轉(zhuǎn)捩位置后移較慢，在周向角?從90°到180°的迎風(fēng)區(qū)轉(zhuǎn)捩位置很快后移；尖錐頭部附近的摩阻系數(shù)比較大，且頭部迎風(fēng)面的摩阻明顯比背風(fēng)面摩阻大；邊界層轉(zhuǎn)捩后表面摩阻系數(shù)顯著增大，迎風(fēng)面中心線和背風(fēng)面中心線上轉(zhuǎn)捩前后摩阻系數(shù)增大一倍以上. 從圖7（b）可看出，沿周向鈍錐表面轉(zhuǎn)捩位置不同，迎風(fēng)面中心線（?=180°）上的轉(zhuǎn)捩位置最靠前，背風(fēng)面中心線（?=0°）上的轉(zhuǎn)捩位置最靠后，沿周向從迎風(fēng)面到背風(fēng)面轉(zhuǎn)捩位置逐漸后移，出現(xiàn)了與圖7（a）中轉(zhuǎn)捩位置沿周向變化方向相反的規(guī)律，即出現(xiàn)了轉(zhuǎn)捩反轉(zhuǎn)，但鈍錐表面轉(zhuǎn)捩位置變化范圍較尖錐的??；鈍頭圓錐頭部的摩阻系數(shù)明顯比尖錐頭部的摩阻系數(shù)小，且周向分布比較均勻；與尖錐相同，邊界層轉(zhuǎn)捩后表面摩阻系數(shù)顯著增大，但轉(zhuǎn)捩前后摩阻系數(shù)的相對變化量更大，轉(zhuǎn)捩后的摩阻系數(shù)約為轉(zhuǎn)捩前摩阻系數(shù)的4～5 倍. 摩阻系數(shù)的突躍對應(yīng)著物體壁面熱流的急劇增大，因此可以判斷，在處于攻角反轉(zhuǎn)對應(yīng)的鈍錐來流條件下，迎風(fēng)面的表面熱防護(hù)將會是其熱設(shè)計(jì)的重點(diǎn).

圖5 尖錐表面間歇因子分布圖（Ma∞=6，α=2°）Fig. 5 Distribution of intermittency on the surface of the sharp cone（Ma∞=6，α=2°）

圖6 鈍錐表面間歇因子分布圖（Ma∞=6，α=2°）Fig. 6 Distribution of intermittency on the surface of the blunt cone（Ma∞=6，α=2°）

圖7 不同周向角對應(yīng)的圓錐表面子午線上流向摩阻系數(shù)分布圖（Ma∞=6，α=2°）Fig. 7 Distribution diagram of friction coefficient of flow direction on the meridian of cone surface corresponding to different circumferential angles（Ma∞=6，α=2°）

為了進(jìn)一步分析隨攻角變化圓錐邊界層出現(xiàn)的轉(zhuǎn)捩反轉(zhuǎn)現(xiàn)象，取反轉(zhuǎn)現(xiàn)象最明顯的攻角α=2°時的尖錐和鈍錐流場進(jìn)行分析. 沿圓錐軸向分別取x/L=0.111、x/L=0.221、x/L=0.332、x/L=0.553、x/L=0.885 的5 個軸向橫截面位置. 根據(jù)圖4 和圖7 中的轉(zhuǎn)捩預(yù)測結(jié)果，可給出不同軸向位置處尖錐和鈍錐迎風(fēng)面和背風(fēng)面邊界層內(nèi)的流態(tài)，見表3. 圖8 和圖9 分別給出了軸向5 個不同橫截面內(nèi)圓錐表面的流向摩阻系數(shù)Cf和名義邊界層厚度δ沿周向的分布曲線.

圖8 圓錐不同橫截面周向摩阻系數(shù)分布（Ma∞=6，α=2°）Fig. 8 Circumferential friction coefficient distribution of different cross-sections of cones（Ma∞=6，α=2°）

圖9 圓錐不同橫截面周向名義邊界層厚度（Ma∞=6，α=2°）Fig. 9 Nominal circumferential boundary layer thickness of different cone cross-sections（Ma∞=6，α=2°）

表3 圓錐不同軸向位置邊界層狀態(tài)（Ma∞=6，α=2°）Tab. 3 The state of boundary layer at different axial positions of cone（Ma∞=6，α=2°）

從圖8 可看到，鈍錐和尖錐在各橫截面內(nèi)表面摩阻系數(shù)沿周向的分布規(guī)律不同，從摩阻系數(shù)的變化規(guī)律可看出鈍錐表面的轉(zhuǎn)捩呈現(xiàn)與尖錐不同的反轉(zhuǎn)現(xiàn)象；結(jié)合圖7 可以看到，發(fā)生轉(zhuǎn)捩反轉(zhuǎn)的鈍錐迎風(fēng)面摩阻系數(shù)明顯大于背風(fēng)面的摩阻系數(shù)，但總體來說鈍錐表面摩阻系數(shù)小于尖錐表面的摩阻系數(shù).

從圖8（a）可見，沿軸向尖錐的流向摩阻系數(shù)的周向分布規(guī)律變化較大. 在靠近頭部x/L=0.111 的橫截面，尖錐周向流態(tài)為層流，此時的迎風(fēng)面摩阻高于背風(fēng)面摩阻，從迎風(fēng)面的中心線沿周向到背風(fēng)面摩阻系數(shù)一直在減??；在x/L=0.221 的橫截面內(nèi)，沿周向尖錐表面摩阻系數(shù)從迎風(fēng)面到背風(fēng)面摩阻系數(shù)逐漸減小，分布規(guī)律與x/L=0.111 橫截面內(nèi)摩阻系數(shù)分布規(guī)律基本一致，但摩阻系數(shù)值明顯減小，此時尖錐背風(fēng)面將要發(fā)生轉(zhuǎn)捩，摩阻水平尚未發(fā)生突躍；當(dāng)軸向從x/L=0.221 到x/L=0.332 時，可以看到尖錐背風(fēng)面流向摩阻發(fā)生了突躍，這說明此時背風(fēng)面已從層流發(fā)展為湍流，摩阻系數(shù)顯著提高，同時可以看到在此橫截面內(nèi)迎風(fēng)面的摩阻系數(shù)明顯比背風(fēng)面的小，沿周向從迎風(fēng)面到背風(fēng)面摩阻系數(shù)在靠近兩側(cè)最突出的位置有陡增的現(xiàn)象；在x/L=0.553 的橫截面內(nèi)，尖錐表面摩阻系數(shù)沿周向的分布與x/L=0.332 時類似，沿周向摩阻陡增處為周向轉(zhuǎn)捩位置；沿流向，尖錐轉(zhuǎn)捩位置逐漸向迎風(fēng)面移動，在x/L=0.885 的橫截面內(nèi)轉(zhuǎn)捩位置已移到迎風(fēng)面中心線附近，除了尖錐迎風(fēng)面中心線兩邊很小的區(qū)域外其他都已轉(zhuǎn)捩為湍流流動，因此迎風(fēng)面摩阻小于其它周向位置；在x/L=0.332之后的橫截面內(nèi)，尖錐背風(fēng)面的均為湍流流動，摩阻系數(shù)值比較接近. 由圖8（b）可知，鈍錐表面流向摩阻分布呈現(xiàn)了迎風(fēng)面始終高于背風(fēng)面，且沿流向摩阻逐漸增大的特點(diǎn)；當(dāng)軸向位置從x/L=0.221 到x/L=0.332時，迎風(fēng)面摩阻陡增，背風(fēng)面摩阻明顯減小，此時迎風(fēng)面邊界層發(fā)生轉(zhuǎn)捩，而背風(fēng)面處于層流區(qū)；當(dāng)軸向位置從x/L=0.332 到x/L=0.553 時，鈍錐背風(fēng)面摩阻陡增，此時鈍錐背風(fēng)面發(fā)生轉(zhuǎn)捩；鈍錐表面在x/L=0.885與x/L=0.553 截面的摩阻沿軸向分布規(guī)律基本一致，迎風(fēng)面摩阻系數(shù)比背風(fēng)面的大，該兩個截面內(nèi)邊界均是已轉(zhuǎn)捩后的湍流流動.

從圖9 可以看到，在有攻角的條件下尖錐與鈍錐的邊界層厚度沿流向一直在增大，只是在不同的周向位置邊界層厚度增大的幅度不同，導(dǎo)致在各橫截面內(nèi)呈現(xiàn)沿周向邊界層厚度不相等的分布. 對于尖錐，其在周向各橫截面內(nèi)沿周向邊界層厚度呈現(xiàn)基本相似的不對稱分布規(guī)律，即背風(fēng)面邊界層厚度大于迎風(fēng)面厚度；沿軸向從x/L=0.111 到x/L=0.221 時，圖9（a）中尖錐背風(fēng)面中心線附近邊界層厚度陡增，此時對應(yīng)背風(fēng)面開始發(fā)生邊界層分離；而當(dāng)沿軸向從x/L=0.553 到x/L=0.885 時，尖錐迎風(fēng)面中心兩側(cè)出現(xiàn)邊界層厚度陡增的現(xiàn)象，這與前面分析的邊界層轉(zhuǎn)捩規(guī)律相同. 對于鈍錐，在軸向x/L=0.111 到x/L=0.332的橫截面內(nèi)迎風(fēng)面邊界層厚度大于于背風(fēng)面邊界層厚度；在軸向x/L=0.553 到x/L=0.885 時，邊界層厚度沿流向明顯增厚，且背風(fēng)面邊界層厚度大于迎風(fēng)面邊界層厚度，圓錐在這兩個截面沿周向均已發(fā)展為湍流，邊界層厚度呈現(xiàn)的是有攻角條件下的湍流邊界層特征. 沿流向鈍錐最大邊界層厚度從迎風(fēng)到背風(fēng)的這種變化規(guī)律，與其邊界層轉(zhuǎn)捩反轉(zhuǎn)緊密相關(guān).

2.2 隨鈍度變化的轉(zhuǎn)捩反轉(zhuǎn)

本節(jié)基于上文的數(shù)值模擬方法，對半錐角θ=7°，長度L=0.8 m 的圓錐通過改變頭部鈍度半徑，研究頭部鈍度對圓錐邊界層轉(zhuǎn)捩的影響，頭部鈍度Rn變化范圍為0.01 mm～20 mm，來流靜溫T∞=57.81 K，來流靜壓P∞=696.7 Pa，壁面采用等溫壁面邊界條件，Tw=300 K.采用與文獻(xiàn)[8]在FD-14A 中進(jìn)行高超聲速鈍頭體邊界層轉(zhuǎn)捩試驗(yàn)相同的條件，在0°攻角下進(jìn)行數(shù)值模擬計(jì)算，取自由流湍流度Tu=0.2%，湍流黏性比μt/μ=10，具體的數(shù)值模擬計(jì)算條件見表4. 表4 中也給出了由數(shù)值模擬得到的圓錐表面邊界層轉(zhuǎn)捩位置距頭部頂點(diǎn)的距離xt. 表中的ReR和Re xt分別為以頭部鈍度、轉(zhuǎn)捩位置至坐標(biāo)原點(diǎn)距離為特征長度的雷諾數(shù).在0°攻角下，對具有不同頭部鈍度的圓錐的轉(zhuǎn)捩繞流場進(jìn)行數(shù)值模擬，研究頭部鈍度對圓錐邊界層轉(zhuǎn)捩的影響.

表4 不同鈍度圓錐轉(zhuǎn)捩計(jì)算條件及轉(zhuǎn)捩數(shù)值預(yù)測結(jié)果Tab. 4 Calculation conditions for transition of cones with different bluntness and numerical prediction results of transition

圖10 給出了不同頭部鈍度和馬赫數(shù)下的轉(zhuǎn)捩位置數(shù)值模擬結(jié)果和實(shí)驗(yàn)結(jié)果的對比，為了方便將影響轉(zhuǎn)捩位置的變量進(jìn)行耦合研究，橫坐標(biāo)為Ma2sinθ以10 為底的對數(shù)，縱坐標(biāo)為轉(zhuǎn)捩雷諾數(shù)Rext以10 為底的對數(shù). 需要說明的是，因文獻(xiàn)[8]給出的頭部鈍度和半錐角數(shù)據(jù)為范圍量，本文所涉及的計(jì)算條件均在實(shí)驗(yàn)雷諾數(shù)Re和馬赫數(shù)Ma范圍以內(nèi)，頭部鈍度和半錐角數(shù)值也在對應(yīng)范圍以內(nèi)，因此圖10 中給出數(shù)值計(jì)算結(jié)果時的橫坐標(biāo)與實(shí)驗(yàn)結(jié)果的橫坐標(biāo)并不完全一致.

圖10 不同馬赫數(shù)下頭部鈍度雷諾數(shù)與轉(zhuǎn)捩雷諾數(shù)的關(guān)系（α=0°）Fig. 10 The relationship between the head dullness Reynolds number and the transition Reynolds number under different Mach numbers

由圖10（a）可知，在來流馬赫數(shù)Ma∞=6，攻角α=0°時，圓錐轉(zhuǎn)捩位置隨著頭部鈍度的增大呈現(xiàn)了“推遲-提前-推遲”的非線性變化，即前文所述由鈍度變化引起的“N”型轉(zhuǎn)捩反轉(zhuǎn)趨勢，該趨勢和實(shí)驗(yàn)轉(zhuǎn)捩結(jié)果十分相似. 圖10（b）中實(shí)驗(yàn)所給的轉(zhuǎn)捩位置分布呈“Λ”型反轉(zhuǎn)，即隨著頭部鈍度的增大，轉(zhuǎn)捩位置變化趨勢呈“提前-推遲”的變化. 而在來流馬赫數(shù)Ma∞=9 時，對應(yīng)工況的數(shù)值模擬預(yù)測結(jié)果仍然呈現(xiàn)“N”型反轉(zhuǎn)的分布，并且兩次反轉(zhuǎn)的鈍度雷諾數(shù)更加接近，但計(jì)算結(jié)果中未能復(fù)現(xiàn)“Λ”型反轉(zhuǎn)的轉(zhuǎn)捩位置分布，下面針對Ma∞=6 條件下的轉(zhuǎn)捩鈍度反轉(zhuǎn)工況作進(jìn)一步分析.

圖11 給出了馬赫數(shù)Ma∞=6 時對應(yīng)的不同頭部鈍度圓錐子午線摩阻系數(shù)Cf分布圖. 從圖11 中摩阻系數(shù)分布可看出轉(zhuǎn)捩位置變化規(guī)律，在頭部鈍度Rn從0.01 mm 增大到0.05 mm，圓錐表面邊界層轉(zhuǎn)捩位置明顯向下游移動，摩阻系數(shù)峰值和湍流區(qū)摩阻水平基本維持不變，層流區(qū)的摩阻系數(shù)最小值隨著轉(zhuǎn)捩位置的后移而減小；在頭部鈍度Rn從0.05 mm 增大到2.5 mm 時，圓錐轉(zhuǎn)捩位置轉(zhuǎn)變?yōu)橄蛏嫌我苿?，摩阻系?shù)峰值和湍流區(qū)摩阻水平依舊基本不變，層流區(qū)的摩阻系數(shù)最小值不斷增大，并且隨著鈍度的增大，摩阻系數(shù)最小值的增長幅度減??；在頭部鈍度Rn從2.5 mm 增大到10 mm 時，轉(zhuǎn)捩位置第二次向下游移動，本次移動過程中的摩阻變化特征與第一次轉(zhuǎn)捩推遲相似；Rn=20 mm 時，轉(zhuǎn)捩位置繼續(xù)向下游移動至模型外，并且移動速度顯著增快，使得圓錐表面的邊界層維持在層流狀態(tài). 可見，隨著圓錐頭部鈍度的變化，邊界層轉(zhuǎn)捩出現(xiàn)“N”型反轉(zhuǎn)的過程中圓錐表面的摩阻系數(shù)的大小及陡增位置也發(fā)生了相應(yīng)的變化.

圖11 Ma∞=6，不同頭部鈍度圓錐子午線摩阻分布圖（α=0°）Fig. 11 Ma∞=6, the distribution diagram of the cone meridian friction with different head bluntness（α=0°）

為了進(jìn)一步從流場分析頭部鈍度對圓錐邊界層轉(zhuǎn)捩的影響規(guī)律，圖12 給出了頭部鈍度Rn=0.01～10 mm的圓錐周圍速度分布云圖，圖中還給出了表示熵層外緣的流線，右側(cè)圖為圓錐頭部頂端附近局部放大圖. 根據(jù)文獻(xiàn)中的熵層外緣確定方法[18]，定義激波后熵變ΔS，其表達(dá)式為

取激波線上滿足ΔS=0.01ΔS0的點(diǎn)作為熵層外緣起點(diǎn)，ΔS0為頭部駐點(diǎn)熵增值. 利用激波后沿流線等熵的性質(zhì)，取激波后沿熵層外緣的流線來表示熵層外緣. 從圖12 可以看到，在圓錐頭部鈍度很小時，基于數(shù)值模擬計(jì)算得到的流場數(shù)據(jù)根據(jù)以上經(jīng)驗(yàn)公式求得的熵層很小，因此可以近似認(rèn)為在尖錐情況下不存在熵層，這一結(jié)果與文獻(xiàn)[1]中的結(jié)論一致；隨著圓錐頭部鈍度增大，圓錐頭部激波后的熵層高度也隨之基本成比例增加. 沿軸向在距圓錐頭部頂點(diǎn)距離xn=Rn處取橫截面.

圖12 Ma∞=6，不同頭部鈍度圓錐速度分布云圖以及熵層外緣示意圖（α=0°）Fig. 12 Ma∞=6, cone velocity distribution cloud diagram and schematic diagram of the outer edge of the entropy layer（α=0°）

圖13 給出了該橫截面內(nèi)不同鈍度圓錐表面的邊界層厚度周向分布圖，結(jié)合圖12 可以看到當(dāng)頭部鈍度相對很小時，圓錐的名義邊界層厚度以及熵層高度都很小，此時圓錐頭部邊界層還未發(fā)展起來. 當(dāng)圓錐頭部鈍度增大到一定值時，其頭部的邊界層厚度也逐漸增厚，同時頭部激波后的熵層高度也隨之增大，隨著圓錐頭部激波后熵層區(qū)域增大，熵層擾動越增強(qiáng)[9]，當(dāng)熵層高頻擾動進(jìn)入邊界層時，會使得邊界層失穩(wěn)，導(dǎo)致圓錐邊界層轉(zhuǎn)捩提前，這可能是圓錐鈍度轉(zhuǎn)捩反轉(zhuǎn)的內(nèi)在原因.

圖13 不同頭部鈍度錐前緣名義邊界層厚度周向分布圖Fig. 13 Circumferential distribution of the nominal boundary layer thickness of the leading edge of cones with different head bluntness

3 結(jié) 論

本文基于雷諾平均的N-S 方程，采用γ轉(zhuǎn)捩模型，在高超聲速條件下對圓錐體的繞流場進(jìn)行了數(shù)值模擬，研究了圓錐邊界層轉(zhuǎn)捩位置隨攻角、鈍度的變化規(guī)律，分析了圓錐邊界層分別隨攻角、頭部鈍度變化出現(xiàn)的轉(zhuǎn)捩反轉(zhuǎn)現(xiàn)象及流動機(jī)理，主要結(jié)論如下：

①隨攻角增大，頭部鈍度較小的尖錐邊界層轉(zhuǎn)捩呈現(xiàn)一般的迎風(fēng)面轉(zhuǎn)捩延遲、背風(fēng)面轉(zhuǎn)捩提前的非對稱轉(zhuǎn)捩規(guī)律；當(dāng)頭部鈍度增大到特定值時，圓錐邊界層轉(zhuǎn)捩呈現(xiàn)迎風(fēng)面轉(zhuǎn)捩提前、背風(fēng)面轉(zhuǎn)捩推遲的轉(zhuǎn)捩反轉(zhuǎn)現(xiàn)象.

②頭部鈍度較大的鈍錐隨攻角增大出現(xiàn)邊界層轉(zhuǎn)捩反轉(zhuǎn)的情況下，迎風(fēng)面摩阻沿流向一直處于周向最高水平，并在轉(zhuǎn)捩過后較于層流區(qū)摩阻系數(shù)約提高了4～5 倍. 相較于一般常規(guī)的非對稱轉(zhuǎn)捩情況，攻角反轉(zhuǎn)情況下的迎風(fēng)面轉(zhuǎn)捩前后摩阻變化更大，因此在對應(yīng)條件下的熱防護(hù)工作重點(diǎn)在于錐體迎風(fēng)面.

③高超聲速零攻角下隨頭部鈍度增大圓錐邊界層轉(zhuǎn)捩位置會發(fā)生 “N”型反轉(zhuǎn)現(xiàn)象，即圓錐表面邊界層轉(zhuǎn)捩位置呈現(xiàn)“推遲-提前-推遲”的變化規(guī)律.隨著來流馬赫數(shù)增大轉(zhuǎn)捩位置呈現(xiàn)“N”型的兩次反轉(zhuǎn)鈍度雷諾數(shù)更加接近.

④圓錐頭部鈍度增大會引起激波后熵層高度和邊界層厚度的增大，圓錐邊界層轉(zhuǎn)捩隨頭部鈍度增大出現(xiàn)的反轉(zhuǎn)現(xiàn)象的可能原因是當(dāng)頭部鈍度增大到一定值時引起頭部激波后熵層高度顯著增大造成的.

本文基于轉(zhuǎn)捩模型通過數(shù)值模擬方法對圓錐邊界層轉(zhuǎn)捩的隨攻角和頭部鈍度變化出現(xiàn)的轉(zhuǎn)捩進(jìn)行了初步研究，下一步需要結(jié)合高精度數(shù)值模擬方法對轉(zhuǎn)捩反轉(zhuǎn)的產(chǎn)生機(jī)理作進(jìn)一步分析，研究熵層高頻擾動進(jìn)入邊界層時對邊界層轉(zhuǎn)捩的具體影響.