馬祎蕾, 余 平, 劉 璟, 梁偉棟, 王建林

(中國運載火箭技術研究院空間物理重點實驗室, 北京 100076)

引 言

高速飛行器邊界層流場發(fā)生轉捩時, 飛行器氣動特性及表面熱環(huán)境均會發(fā)生顯著改變, 正確預測轉捩位置, 對飛行器設計與優(yōu)化具有十分重要的工程意義[1]. 轉捩誘因繁多, 觸發(fā)機制復雜, 受到轉捩Reynolds數(shù)、 攻角、 鈍度、 邊界層邊緣Mach數(shù)等多重因素影響[1-4], 其中, 外形鈍度是影響邊界層流動穩(wěn)定性的一個重要因素.

以往關于高速邊界層轉捩鈍度影響的研究很多, 但多為針對鈍度半徑、 前緣曲率的研究, 文獻[4]指出, 對于軸對稱外形, 適度的端頭鈍度會顯著地推遲轉捩, 對于面對稱的升力體外形, 適度的前緣鈍度亦會顯著地推遲轉捩; Stetson等[5]通過地面試驗發(fā)現(xiàn)在Mach數(shù)5.5的激波風洞中, 鈍錐端頭鈍化可有效推遲轉捩, 然而, 鈍度一旦超過某一臨界值, 轉捩將提前, 出現(xiàn)“轉捩反轉”, 并在1984年發(fā)現(xiàn)這種轉捩延遲效應與鈍度帶來的熵層有關, 存在“熵吞”現(xiàn)象[6]; Malik等[7]采用PNS方法研究發(fā)現(xiàn)鈍度使邊界層穩(wěn)定. Rosenboom等[8]發(fā)現(xiàn)當鈍度增加時第2模態(tài)臨界Reynolds數(shù)單調(diào)增加, 但第1模態(tài)臨界Reynolds數(shù)在大鈍度和小鈍度情況下單調(diào)性不同; Kara等[9]采用eN法, 將N=10作為轉捩判據(jù), 發(fā)現(xiàn)隨著鈍度增大, 轉捩Reynolds數(shù)變大; 萬兵兵等[10], 歐吉輝等[11]研究了鈍度對激波后熵層基本流特性及熵層不穩(wěn)定性的影響.

實際上, 除了鈍度半徑會影響轉捩外, 鈍度形狀也是構成鈍度的一部分, 據(jù)了解, 改變端頭形狀亦會對下游的邊界層轉捩產(chǎn)生影響, 選擇合適的端頭形狀能夠起到推遲轉捩的作用, 然而, 這方面的資料相對很少. Lin等[12]研究發(fā)現(xiàn)前緣橢球比增大時, 邊界層內(nèi)被激發(fā)的T-S波幅值減小, 反之幅值增大; Lin等[13]研究發(fā)現(xiàn), 前緣頂端曲率可通過前緣Reynolds數(shù)影響邊界層流動穩(wěn)定性; 沈露予等[14]通過DNS方法研究了自由來流湍流下不同橢前緣形狀的平板邊界層感受性問題.

目前, 對高速邊界層常用的轉捩預測方法主要有3類: 轉捩準則、 轉捩模式法與eN法, 其中eN法[15-17]的理論基礎即線性穩(wěn)定性理論(linear stability theory, LST), 并最大限度地利用了邊界層擾動演化的理論預測[18-19]. 盡管其忽略了邊界層的感受性、 非線性作用等因素, 仍被Cebeci等波音工程師認為是最有效的轉捩預測方法.

考慮到采用升力體外形的高速飛行器, 多包含鈍化的端頭和翼前緣、 大后掠角的翼面形狀及較為扁平的迎風面構型等, 三維幾何特征明顯, 因此, 本文選擇平板鈍三角翼外形, 圍繞典型高速風洞實驗狀態(tài)開展流場數(shù)值計算和線性穩(wěn)定性分析, 并采用eN法獲取表面N值分布進行轉捩預測, 研究在LST理論背景下, 不同橢前緣形狀對邊界層流動穩(wěn)定性特征的影響, 并探究前緣鈍度形狀對邊界層轉捩的影響.

1 計算方法

1.1 數(shù)值計算方法

N-S方程數(shù)值求解中, 采用有限體積法, 對流項采用3階迎風Roe-FDS格式離散, 黏性項采用中心差分格式離散, 時間推進為AF-LU隱式方法.

不考慮對底部流場的模擬, 物面采用無滑移邊界條件, 出口設為超聲速外插邊界條件. 計算域為半模, 經(jīng)網(wǎng)格收斂性驗證, 確定網(wǎng)格量為221×241×351, 如圖1所示, 法向網(wǎng)格布置點數(shù)較多, 原因在于后續(xù)穩(wěn)定性分析需要較為精細的邊界層流場.

圖1 計算域及網(wǎng)格Fig. 1 Computational domain and meshes

1.2 穩(wěn)定性計算方法

根據(jù)線性穩(wěn)定性理論, 二維局部平行流假設下, 分析O-S方程可將小擾動q表示為行進波形式

(1)

(2)

式中, (x0,z0)為頻率為ω的擾動波開始失穩(wěn)處(αi,βi由正轉為負), 其位置在中性曲線上; 增長率-αi， -βi取當?shù)仡l率為ω的最不穩(wěn)定擾動波(增長率最大)的增長率, 擾動傳播路徑方向即該擾動波的群速度方向, 由于群速度方向與勢流方向非常接近, 本文采用勢流方向作為擾動波增長率積分方向.

預設轉捩判據(jù)為Ntr(通常由經(jīng)驗給出), 即對于不同頻率ω的擾動波, 若在(xtr,ztr)處幅值放大因子達到Ntr, 預示轉捩發(fā)生. 不同頻率擾動波擾動幅值達Ntr的位置不同, 一般取最上游位置作為轉捩發(fā)生位置, 對應N值稱為包絡N值. 轉捩位置曲線應滿足

(3)

1.3 算例驗證

文獻[20]針對HIFiRE-1飛行工況進行了數(shù)值模擬, 鈍錐半錐角7°, 頭部半徑2.5 mm, 計算工況見表1.

為考核三維流場計算的可靠性, 特別是邊界層流場和橫流效應, 本文采用目前的方法也進行了計算, 計算所得流向x=850 mm, 背風面θ=135°位置點的基本流速度剖面與文獻[20]結果比較如圖2所示, 其中橘色速度剖面為計算結果, 與算例基本相符, 驗證了邊界層流場計算的正確性.

表1 算例工況

(a) Streamwise velocity

(b) Crossflow velocity

(c) Wall-normal velocity圖2 算例對比Fig. 2 Comparison of numerical examples

2 幾何模型與計算工況

幾何模型為一平板鈍三角翼外形[21], 長度取400 mm, 鈍度半徑3 mm, 后掠角75°. 定義坐標系原點O為鈍三角板頂點,x,z軸指向如圖3(a)所示,y軸符合右手坐標系.

定義橢圓前緣截面的長短軸比為形狀因子m, 則m=a/b(如圖3 (c)所示), 保持端頭鈍度半徑為3 mm, 即保持b恒為3 mm, 分別改變前緣形狀因子為0.5, 1和2, 以m050,m100和m200表示. 圖3(b)為m200外形迎風面示意圖, 圖3(c)為(b)中m-n截面上3種前緣形狀的橢圓截面分布.

根據(jù)靜風洞實驗結果, 選取典型風洞實驗狀態(tài)(如表2所示)為計算工況.

(a) Model and coordinate system

(b) Windward side of m200

(c) m-n圖3 平板鈍三角翼外形及坐標系示意圖Fig. 3 Model and coordinate system of flat-plate blunt delta wing

表2 計算工況

3 橢前緣形狀影響分析

3.1 基本流分析

圖4(a), (b)分別為m100前緣鈍三角翼不同流向位置處Mach數(shù)分布圖, 其迎、 背風邊界層均從前緣向中心逐漸增厚, 并在對稱中心線附近出現(xiàn)流向渦結構, 有攻角狀態(tài)下, 迎風面流場特征相對簡單, 流向渦尺度明顯小于背風面, 故本文僅對迎風面進行分析.

(a) Windward side

(b) Leeward side圖4 不同流向位置Mach云分布Fig. 4 Mach number contours at different streamwise locations

定義最大橫流速度為沿法向分布的橫流速度最大值, 以來流速度U∞= 860.838 m/s對橫流速度無量綱化, 則迎風面最大無量綱橫流速度分布如圖5所示. 可知正攻角工況下, 鈍三角翼流向前部, 靠近前緣位置處最大橫流速度較大; 沿流向發(fā)展, 前緣處最大橫流速度逐漸減小; 同一流向位置截面上, 展向中間區(qū)域(去除前緣和中心區(qū)域)最大橫流速度較大.

圖5 迎風面最大橫流速度分布Fig. 5 Distribution of maximum cross flow velocity on windward surface

由于中心流向渦結構區(qū)域流動復雜, 流場變化劇烈, 因此, 采用線性穩(wěn)定性理論分析時應避開該區(qū)域. 故基于m100鈍三角翼外形, 特別是中心流向渦尺度特征和邊界層橫流流動特征, 選擇橫流效應相對較強, 位于展向中間區(qū)域的point(300,-3,40)作為典型點, 來分析不同前緣形狀對鈍三角翼不同模態(tài)流動穩(wěn)定性特征影響.

圖6為不同前緣形狀模型流向x=300 mm的迎風面邊界層厚度分布, 邊界層厚度由總焓確定, 為由壁面開始過峰值后等于來流總焓的1.005倍位置處[10]. 分析可知前緣形狀僅對前緣區(qū)域邊界層厚度產(chǎn)生影響, 翼身區(qū)域的邊界層厚度線性分布區(qū)域基本相同.

圖6 x=300 mm邊界層厚度Fig. 6 Boundary layer thickness at x=300 mm

不同前緣形狀三角翼外形的迎風面最大橫流速度分布云圖見圖7(a)～(c), 可知前緣形狀主要影響鈍三角翼流向前部, 靠近前緣的具有強橫流效應區(qū)域, 且前緣形狀越尖, 最大橫流速度越大, 具有較強橫流效應區(qū)域越廣; 前緣形狀對三角翼翼身最大橫流速度分布影響較小.

(a) m050

(b) m100

(c) m200圖7 迎風面最大橫流速度Fig. 7 Distributions of maximum cross flow velocity

3.2 穩(wěn)定性分析

由圖5, 6可知, point(300,-3,40)位于展向40 mm處, 不在中心流向渦影響區(qū)域內(nèi), 故可作為后續(xù)線性穩(wěn)定性特征分析點.

具有不同橢前緣截面形狀外形過point點的邊界層外緣勢流方向如圖8(a)所示, 勢流方向基本重合. 圖8(b)～(d)分別為m050,m100和m200 這3種外形沿該勢流方向不同頻率擾動波增長率分布曲線.

(a) Streamline passing point 2

(b) m050

(c) m100

(d) m200圖8 勢流方向與增長率沿勢流分布Fig. 8 Streamline and streamwise distribution of grouth rate

比較可知, 不同前緣形狀主要影響失穩(wěn)起始區(qū)域(前緣附近區(qū)域)的穩(wěn)定性特征分布, 形狀越尖, 前緣附近第1模態(tài)增長率越大, 而勢流方向下游遠離前緣區(qū)域的第1模態(tài)增長率分布基本相同; 第2模態(tài)增長率則隨著前緣變尖而逐漸增大. 前緣形狀對不同模態(tài)開始出現(xiàn)擾動增長的流向位置產(chǎn)生不同影響, 前緣越尖, 第1模態(tài)的流向(x向)中性點位置越靠后, 第2模態(tài)的中性點位置則略微前移.

選取f=0 kHz(橫流駐波)， 30 kHz(第1模態(tài))和200 kHz(第2模態(tài))3條擾動波, 其增長率沿圖8(a)中的勢流方向的分布如圖9所示. 分析可知, 同一頻率擾動波沿勢流方向發(fā)展, 前緣形狀不同, 增長率在前緣的分布完全不同.

(a) f=0 kHz

(b) f=30 kHz

(c) f=200 kHz圖9 增長率分布Fig. 9 Distributions of growth rate

圖9(a)為橫流駐波(f=0 kHz), 前緣較尖時(m200), 橫流效應強, 前緣處橫流駐波增長率大, 導致增長率沿勢流方向出現(xiàn)局部峰值點, 隨后緩慢降低, 而前緣較鈍時(m050), 從前緣過渡到翼身區(qū)域, 橫流駐波增長率平緩增長再下降, 未出現(xiàn)增長率局部峰值點; 圖9(b)f=30 kHz擾動波在不同前緣形狀的增長率分布趨勢與f=0 kHz類似, 但增長率幅值差異有所降低; 以圖9(c)f=200 kHz 表征第2模態(tài)增長率分布, 則前緣變尖, 增長率變大, 同時中性點位置略微前移.

3.3 轉捩預測

目前工程上對于T-S波轉捩與橫流轉捩同時存在的情況, 將取兩個不同的轉捩N值, 即NTS,tr與NCF,tr, 再分別積分NTS與NCF, 二者之一達臨界值時即預示轉捩發(fā)生[22]. 本文將NTS,tr與NCF,tr取為同一個值, 并統(tǒng)一歸至第1模態(tài)進行轉捩預測.

如圖10所示, (a)為包括橫流模態(tài)的第1模態(tài)包絡N值分布, (b)為第2模態(tài)包絡N值分布. 顯然, 第1模態(tài)N值遠高于第2模態(tài), 即迎風面轉捩為第1模態(tài)導致. 經(jīng)與靜風洞溫敏漆轉捩測量結果對比, 轉捩位置處Ntr約為3. 上文中點point位于轉捩區(qū)域, 作為參考, 該位置處N值約為4, 失穩(wěn)波頻為30 kHz, 為第1模態(tài)失穩(wěn).

(a) First modes(f=0～76 kHz)

(b) Second modes(f=77～328 kHz)圖10 不同模態(tài)包絡N值分布Fig. 10 N-factor contours for different mode disturbances

圖11(a)～(c)分別為不同橢前緣外形, 頻率f=0, 30, 200 kHz擾動波N值沿勢流方向分布. 由圖可知, 隨著前緣變尖,f=0 kHz擾動波出現(xiàn)一定量值后同一N值的流向位置逐漸前移, 頻率f=30 kHz 的擾動波出現(xiàn)一定量值后同一N值的流向位置逐漸后移, 如表3所示為該擾動波的N值達到轉捩N值(Ntr=3)的流向坐標, 而頻率f=200 kHz的擾動波N值隨著前緣變尖而增大.

(a) f=0 kHz

(c) f=200 kHz圖11 N值分布Fig. 11 Distributions of N-factor

表3 f=30 kHz擾動波Ntr=3流向位置

結合圖9可知, 不同橢前緣截面形狀外形, 前緣處增長率差異較大, 同時沿勢流方向不穩(wěn)定區(qū)域長度(積分區(qū)間)也不同, 前緣越尖增長率越大, 但相應擾動中性點越靠后, 積分區(qū)間越短, 而N值是增長率沿勢流方向積分的結果, 故f=0 kHz時增長率起主導作用, 前緣變尖加速失穩(wěn),f=30 kHz時增長率差異小, 積分區(qū)間起主導作用, 前緣變尖積分區(qū)間變短, 可推遲轉捩. 即前緣形狀變尖, 將抑制第1模態(tài)失穩(wěn), 促進橫流模態(tài)和第2模態(tài)失穩(wěn).

不同橢前緣截面形狀外形的轉捩陣面(Ntr=3)如圖12所示. 3種外形均為第1模態(tài)轉捩, 沿展向呈“M型”分布, 隨著前緣形狀變尖, 流向轉捩位置逐漸后移, 展向轉捩位置幾乎無變化.

圖12 轉捩陣面Ntr=3Fig. 12 Transition locations at Ntr=3

4 結論

本文研究了在典型風洞實驗狀態(tài)(Re∞= 3×106,Ma∞= 6,α= 10°)下, 不同橢前緣截面形狀對迎風面平板鈍三角翼高速邊界層流動穩(wěn)定性及轉捩的影響, 所得結論如下:

(1)前緣形狀將影響前緣附近區(qū)域邊界層流場特征, 對遠離前緣區(qū)域影響較小, 且前緣形狀越尖, 橫流速度越大.

(2)前緣形狀對不同模態(tài)影響不同: 前緣形狀變尖, 一定頻率的第1模態(tài)和橫流模態(tài)擾動波的中性點位置后移, 擾動增長率變大; 一定頻率的第2模態(tài)擾動波增長率變大, 中性點位置基本不變.

(3)由靜風洞實驗結果標定, 本計算工況下轉捩N值Ntr=3, 轉捩預測結果為第1模態(tài)導致轉捩; 隨著前緣形狀變尖, 受第1模態(tài)中性點位置后移影響, 流向轉捩位置逐漸后移, 展向失穩(wěn)位置幾乎不變.