張彥軍，趙軻，張同鑫，陳利麗

(第一飛機(jī)設(shè)計研究院 總體氣動設(shè)計研究所，西安 710089)

0 引 言

在航空航天領(lǐng)域，失速問題對飛行器設(shè)計和使用造成惡劣影響，甚至導(dǎo)致機(jī)毀人亡。為了滿足飛行安全和改善起降性能，失速問題一直是氣動設(shè)計師必須面對和解決的基礎(chǔ)科學(xué)問題。探索飛行器的失速問題隱含的空氣動力學(xué)機(jī)理是明確失速機(jī)制及建立控制措施的必由之路，但是復(fù)雜的部件干擾和雷諾數(shù)變化增加了研究失速機(jī)理的難度。幸運的是，雖然難以給出飛行器失速問題精確的解析關(guān)系式，但組成這種強非線性現(xiàn)象的核心機(jī)制仍然是翼型的失速問題，因此研究翼型的失速問題是建立飛行器失速問題理論體系的基石。

從20世紀(jì)上半葉開始，學(xué)者們采用試驗方法對翼型失速的表面氣動力和邊界層變化進(jìn)行了諸多分析，E.N.Jacobs[1]第一次推測：升力線的拐點意味著翼型表面出現(xiàn)了分離現(xiàn)象，即提出失速問題；隨后，E.N.Jacobs等[2]第一次提出雷諾數(shù)與翼型分離之間存在關(guān)系，重點關(guān)注層流轉(zhuǎn)捩現(xiàn)象對翼型前緣和后緣分離的影響；Lissaman[3]闡述低速翼型與雷諾數(shù)之間的影響作用，并分析雷諾數(shù)影響下分離泡的延伸與破裂是導(dǎo)致翼型失速的主要原因；V.O.Michael等[4]基于PIV技術(shù)研究低雷諾數(shù)翼型分離泡的相關(guān)特性；Hu H等[5]使用PIV檢測GA(W)-1翼型流場和表面壓力分布，分析層流分離泡與壓力梯度的相互影響；H.Yamato等[6]使用動態(tài)控制裝置影響前緣分離泡破碎，從而抑制翼型分離。國內(nèi)， Wang N等[7]通過試驗研究雷諾數(shù)對RAE2822翼型氣動特性的影響；孫智偉等[8]在某NF-3和NF-6風(fēng)洞中進(jìn)行超臨界翼型的設(shè)計和實驗研究，給出超臨界翼型在高、低速時的分離演化特性。

在數(shù)值模擬方面，A.C.Aranake等[9]使用SA湍流模型耦合γ-Reθt轉(zhuǎn)捩模型對S827和S809翼型分別進(jìn)行失速特性數(shù)值模擬，相比于全湍流計算結(jié)果，無論是升力特性還是阻力特性均與試驗結(jié)果更加接近，且考慮邊界層轉(zhuǎn)捩時計算所得氣動力系數(shù)能夠反應(yīng)失速狀態(tài)下的輕失速和深失速特性；Ke J等[10]采用耦合轉(zhuǎn)捩模型的RANS方法和混合方法對A-Airfoil及NACA0012翼型靜態(tài)失速特性進(jìn)行模擬，分析前緣層流分離泡及后緣分離發(fā)展，其結(jié)果與試驗結(jié)果吻合良好。國內(nèi)，馮濤等[11]針對典型氣動負(fù)載渦輪葉片在不同雷諾數(shù)下邊界層特性變化進(jìn)行敏感性研究；靳允立等[12]采用不同湍流模型研究翼型失速及風(fēng)速變化引起的雷諾數(shù)改變對翼型和風(fēng)力機(jī)氣動性能的影響；王強等[13]采用基于SST模型的γ-Reθt轉(zhuǎn)捩模型研究S809翼型動態(tài)失速過程中分離泡的作用；袁尚科[14]通過理論分析、數(shù)值計算及外場實驗相結(jié)合的方法，分別對風(fēng)力機(jī)常用翼型及葉片的靜態(tài)失速、動態(tài)失速及三維旋轉(zhuǎn)效應(yīng)等進(jìn)行研究。

綜上所述，國內(nèi)外諸多學(xué)者將翼型失速的研究對象重點集中于失速點的精確預(yù)測，對于失速時邊界層的變化，層流轉(zhuǎn)捩現(xiàn)象的影響作用和雷諾數(shù)對邊界層演化的完整影響過程并未有過多涉足。

本文使用Menter SST湍流模型和擾動放大因子輸運模型對轉(zhuǎn)捩翼型邊界層進(jìn)行數(shù)值模擬，分析雷諾數(shù)變化時邊界層的發(fā)展，總結(jié)雷諾數(shù)變化對翼型失速的影響規(guī)律。

1 數(shù)值模擬方法

流場模擬采用Navier-Stokes方程組進(jìn)行計算，其守恒積分形式為

(1)

式中：Q為守恒向量；S為控制體V所有表面S的集合；F為通過控制體的S面上的通量；n為S面的單位外法向矢量。無粘項、粘性項的空間離散方法分別為二階迎風(fēng)Roe格式及二階中心差分方法。時間離散則采用近似因子分解方法。

Menterk-ωSST湍流模型[15]是一個兩方程k-ω渦粘性模型，在工程界廣泛應(yīng)用。SST模型是一個耦合模型，在近壁面保持Wilcoxk-ω模型[16]的特性，而在湍流邊界層外及自由來流中保持Launder-Sharmak-ε模型[17]的特性。通過混合函數(shù)保留兩個模型的優(yōu)點并克服模型各自的不足。其方程形式為

(2)

(3)

渦粘性由下式計算得到

(4)

公式中模型參數(shù)(α，β，σ，F(xiàn)1，F(xiàn)2)的數(shù)值及處理方法可以參照文獻(xiàn)[15]。

在模擬層流-湍流轉(zhuǎn)捩邊界層時，采用Coder等[18]提出的擾動放大因子輸運模型。該輸運方程模型形式如式(5)所示。

(5)

式中的Fcrit，F(xiàn)growth，dn/dRδ2均與邊界層形狀因子相關(guān)，關(guān)系式的具體表示可參考文獻(xiàn)[18]。因此模型的核心點為使用層流邊界層的當(dāng)?shù)貕毫μ荻葏?shù)HL計算非當(dāng)?shù)貐?shù)形狀因子H12，后者常用于描述速度型的特征。當(dāng)?shù)貕毫μ荻葏?shù)定義為

(6)

式中：d為壁面法向距離；S為當(dāng)?shù)丶羟袘?yīng)變率；Ue為邊界層邊緣速度。

由Falkner-Skan二維邊界層相似性解可以得到不同壓力梯度下HL最大值與形狀因子的函數(shù)關(guān)系，即

(7)

從而實現(xiàn)在流場中形狀因子H12的求解，邊界層中max(H12)為精確值。

將擾動放大因子輸運方程與SST模型耦合，其湍動能輸運方程為

(8)

(9)

ft2=[1-exp(2(n-Ncrit))]exp(ct4χ2)

(10)

式中：模型參數(shù)F1,ct4,ct5,χ,Ncrit等的設(shè)置與文獻(xiàn)[19]一致。

2 翼型轉(zhuǎn)捩邊界層發(fā)展及失速影響分析

2.1 翼型和網(wǎng)格

使用NLF(1)-0416自然層流翼型繞流[20]研究雷諾數(shù)的影響。該翼型有較長的層流區(qū)且具備高升力特性，為典型后緣失速翼型，一般應(yīng)用于通用航空飛機(jī)翼型和風(fēng)力機(jī)葉型等領(lǐng)域。對其進(jìn)行轉(zhuǎn)捩數(shù)值模擬，計算網(wǎng)格分布如圖 1所示，計算域半徑為100倍弦長，翼型表面布置481個網(wǎng)格單元，法向布置181個網(wǎng)格單元，網(wǎng)格量約為18萬。壁面法向第一層網(wǎng)格高度約為1.0×10-6m，使y+小于1。風(fēng)洞試驗雷諾數(shù)工況分別為2×106,4×106,9×106。風(fēng)洞自由來流馬赫數(shù)為0.1，湍流度設(shè)置為0.03%。

圖1 翼型繞流數(shù)值模擬網(wǎng)格分布

2.2 數(shù)值模擬結(jié)果

為了與試驗狀態(tài)保持一致，便于對比驗證和分析，本文選擇了與試驗雷諾數(shù)相同的三個雷諾數(shù)進(jìn)行數(shù)值模擬研究，鑒于雷諾數(shù)對失速特性的影響一直是一個復(fù)雜的問題，對于更大雷諾數(shù)，我們將在后續(xù)工作中開展研究。使用擾動放大因子輸運模型對三個雷諾數(shù)下翼型繞流進(jìn)行模擬，得到升力系數(shù)-阻力系數(shù)極曲線并與全湍數(shù)據(jù)和試驗數(shù)據(jù)進(jìn)行對比，結(jié)果如圖2所示。

(a) Re=2×106

(b) Re=4×106

(c) Re=9×106

從圖2可以看出：相比全湍模型的數(shù)值模擬結(jié)果，轉(zhuǎn)捩模型的預(yù)測結(jié)果與試驗數(shù)據(jù)更加貼合，與真實流動更加吻合。

雷諾數(shù)為4×106，0°迎角工況的翼型表面壓力系數(shù)分布如圖3所示。

圖3 Re=4.0×106，0°迎角翼型表面的壓力系數(shù)分布

從圖3可以看出：在該工況下，翼型繞流轉(zhuǎn)捩類型為層流分離泡轉(zhuǎn)捩，且轉(zhuǎn)捩模式對該分離再附著過程描述較好。模型對轉(zhuǎn)捩過程的準(zhǔn)確預(yù)測為分析邊界層發(fā)展奠定基礎(chǔ)。

2.3 雷諾數(shù)影響分析

三個雷諾數(shù)下翼型升力系數(shù)隨迎角變化曲線如圖4所示，可以看出：雷諾數(shù)對翼型的氣動特性有較大的影響，當(dāng)雷諾數(shù)增大時，翼型的失速迎角變大，同時最大升力系數(shù)增大。產(chǎn)生此現(xiàn)象的主要原因是雷諾數(shù)變化對翼型邊界層轉(zhuǎn)捩及發(fā)展有重要影響，因此通過對比不同雷諾數(shù)下的轉(zhuǎn)捩過程及邊界層發(fā)展情況，對雷諾數(shù)的影響規(guī)律進(jìn)行分析。

圖4 不同雷諾數(shù)下翼型升力系數(shù)隨迎角變化曲線

不同雷諾數(shù)下翼型上下表面的轉(zhuǎn)捩位置隨升力系數(shù)變化曲線如圖5所示。

圖5 不同雷諾數(shù)翼型轉(zhuǎn)捩位置隨升力系數(shù)的變化

從圖5可以看出：隨著迎角增大，升力系數(shù)增大，翼型上表面的轉(zhuǎn)捩位置逐漸前移，而下表面變化趨勢相反。相同升力系數(shù)下，當(dāng)雷諾數(shù)增大時，上、下表面的轉(zhuǎn)捩位置均逐漸前移。

對0°迎角工況下翼型邊界層參數(shù)隨雷諾數(shù)的變化進(jìn)行分析。0°迎角時翼型的表面壓力分布圖如圖6所示，可以看出：雷諾數(shù)變化對壓力系數(shù)的整體影響較小，其主要影響體現(xiàn)在壓力系數(shù)突變位置的變化，即出現(xiàn)層流分離泡、發(fā)生轉(zhuǎn)捩的位置。

圖6 0°迎角時翼型表面壓力系數(shù)分布對比

不同雷諾數(shù)在0°迎角時翼型的表面摩擦力系數(shù)分布如圖7所示，可以看出：在雷諾數(shù)為2×106和4×106的工況下，翼型上下表面的摩擦力系數(shù)均出現(xiàn)負(fù)值，表面附面層中出現(xiàn)層流分離泡，并伴隨轉(zhuǎn)捩的發(fā)生；而雷諾數(shù)為9×106時，翼型繞流發(fā)生層流到湍流的轉(zhuǎn)捩過程，但是并沒有層流分離泡出現(xiàn)。初步的結(jié)論與傳統(tǒng)的研究結(jié)論一致，低雷諾數(shù)環(huán)境雖然轉(zhuǎn)捩位置比較靠后，但是層流非常容易被逆壓梯度誘導(dǎo)產(chǎn)生層流分離泡，觸發(fā)轉(zhuǎn)捩。該現(xiàn)象在高空長航時無人機(jī)表面比較常見，因為數(shù)萬米高空中空氣相對稀薄，雷諾數(shù)較小，容易產(chǎn)生層流分離泡。

(a) Re=2×106

(b) Re=4×106

(c) Re=9×106

翼型上、下表面層流分離泡形態(tài)如圖8所示，可以看出：當(dāng)雷諾數(shù)增大時，翼型上下表面的層流分離泡均逐漸減小，在雷諾數(shù)為9×106時，層流分離泡已經(jīng)消失。

(a1) 上表面分離泡

(a2) 下表面分離泡

(b1) 上表面分離泡

(b2) 下表面分離泡

(c1) 上表面分離泡

(c2) 下表面分離泡

對比不同雷諾數(shù)下翼型上表面邊界層的邊界層動量損失厚度δ2、形狀因子H12、渦量雷諾數(shù)Reν等參數(shù)分布，如圖9～圖11所示。

圖9 0°迎角上表面邊界層動量損失厚度分布

圖10 0°迎角翼型上表面邊界層形狀因子分布

圖11 翼型上表面邊界層渦量雷諾數(shù)分布

從圖9可以看出：當(dāng)雷諾數(shù)減小時，動量損失厚度增大，這是由于雷諾數(shù)減小使流動粘性作用相對增強，有利于邊界層的發(fā)展。

邊界層形狀因子的定義為

(11)

式中：δ1為邊界層位移厚度。

從圖10可以看出：當(dāng)層流邊界層通過轉(zhuǎn)捩點變?yōu)橥牧鬟吔鐚樱吔鐚雍穸葧蝗辉龊?。湍流邊界層由于流速分布更趨均勻化而使邊界層位移厚度減小，同時由于阻力增加而使邊界層動量損失厚度增大，于是邊界層形狀因子在轉(zhuǎn)捩后將減小。翼型上表面在三個雷諾數(shù)下的轉(zhuǎn)捩位置均在40%弦長附近，因此上表面邊界層形狀因子在40%弦長處出現(xiàn)峰值，而且由于雷諾數(shù)為2×106和4×106時翼型出現(xiàn)層流分離泡，所以這兩個雷諾數(shù)對應(yīng)的形狀因子峰值更高；而雷諾數(shù)為9×106時邊界層未出現(xiàn)分離泡，且轉(zhuǎn)捩位置前移，40%弦長處流動已基本發(fā)展為湍流，對應(yīng)的形狀因子值較低。在轉(zhuǎn)捩發(fā)生后，形狀因子的發(fā)展情況與前文的分析對應(yīng)：三組形狀因子均迅速下降，雷諾數(shù)越大，最終邊界層形狀因子的值越小。

從圖11可以看出：相同弦長站位處，雷諾數(shù)越大，邊界層中的渦量雷諾數(shù)越高。因為渦量雷諾數(shù)與常用轉(zhuǎn)捩判定參數(shù)動量厚度雷諾數(shù)幾乎是正比關(guān)系，因此可以作為邊界層發(fā)展的標(biāo)志性參數(shù)進(jìn)行分析。邊界層轉(zhuǎn)捩觸發(fā)位置隨雷諾數(shù)升高而向前緣移動，在轉(zhuǎn)捩后，渦量雷諾數(shù)迅速增大，使邊界層保持在湍流狀態(tài)。

0°迎角時翼型上表面不同站位流向速度型分布對比如圖12所示。

(a) 0.2c

(b) 0.3c

(c) 0.4c

(d) 0.5c

(e) 0.6c

(f) 0.7c

(g) 0.8c

(h) 0.9c

從圖12可以看出：在翼型前緣，轉(zhuǎn)捩前較高雷諾數(shù)對應(yīng)邊界層的速度型更加均勻、飽滿，邊界層更薄。在40%弦長處，由于翼型邊界層在較低雷諾數(shù)時出現(xiàn)層流分離泡，速度型分布在壁面附近呈凹曲線，速度型出現(xiàn)拐點。轉(zhuǎn)捩后，三個雷諾數(shù)對應(yīng)的速度型更加相近，高雷諾數(shù)的速度型分布更加飽滿。

為了分析雷諾數(shù)對失速迎角附近邊界層發(fā)展的影響，取12°迎角工況下翼型邊界層進(jìn)行觀察。12°迎角時翼型的表面摩擦力系數(shù)分布如圖13所示，隨雷諾數(shù)增高翼型轉(zhuǎn)捩點前移，轉(zhuǎn)捩處層流分離泡逐漸縮小消失，僅雷諾數(shù)為2×106時存在比較明顯的層流分離泡。

圖13 12°迎角翼型表面摩擦力系數(shù)分布

12°迎角時翼型上表面不同站位流向速度型分布對比如圖14所示。由于雷諾數(shù)為2×106時在約7%弦長處產(chǎn)生層流分離泡，其對應(yīng)的速度分布出現(xiàn)明顯的拐折。雷諾數(shù)為9×106時的速度型更早趨于均勻，因其轉(zhuǎn)捩點更加靠近前緣。轉(zhuǎn)捩后，隨雷諾數(shù)的增高，流速分布變得更加均勻而壁面附近流速梯度增大。

翼型上表面邊界層內(nèi)各站位處的湍動能分布如圖15所示。從翼型表面摩擦力系數(shù)分布可知，三個雷諾數(shù)下翼型上表面均已在10%弦長前完成轉(zhuǎn)捩，因此下圖中的結(jié)果均為轉(zhuǎn)捩后邊界層湍動能的分布。

(a) 0.005c(b) 0.010c(c) 0.030c

(d) 0.060c(e) 0.100c(f) 0.200c

(g) 0.400c(h) 0.600c(i) 0.900c

圖14 12°迎角翼型上表面不同站位速度型分布

Fig.14 The velocity distribution of different stations on the upper surface of airfoil at 12° angle of attack

(a) 0.2c

(b) 0.4c

(c) 0.6c

(d) 0.9c

從圖15可以看出：轉(zhuǎn)捩后相同站位低雷諾數(shù)對應(yīng)的湍動能峰值更高。這是由于翼型前緣發(fā)生轉(zhuǎn)捩，在較低雷諾數(shù)時伴隨有層流分離泡形成，這一過程使邊界層流動的脈動特性得到加強，其對應(yīng)的湍動能、渦粘性增大。但是層流分離泡會消耗邊界層的動能，使流動自持能力即抗分離能力減弱，從而產(chǎn)生分離。當(dāng)雷諾數(shù)增高，分離泡減小甚至消失，流動能量消耗減小，動能更充沛，壓力分布可以維持較長距離的梯度，后緣分離滯后發(fā)生，從而增大翼型升力系數(shù)及失速迎角。

最后對比相同迎角不同雷諾數(shù)下翼型流場的后緣分離情況，結(jié)果如圖16～圖17所示。

(a) Re=2×106

(b) Re=4×106

(c) Re=9×106

(a) Re=2×106

(b) Re=4×106

(c) Re=9×106

從圖16～圖17可以看出：在相同迎角下，隨著雷諾數(shù)增大，流動穩(wěn)定性增強，抵抗逆壓梯度干擾的能力提高，流動分離趨勢減弱。后緣分離逐漸減小，這與本節(jié)中對邊界層發(fā)展及其參數(shù)變化趨勢的觀察結(jié)果相一致。

3 結(jié) 論

(1) 對于層流-湍流轉(zhuǎn)捩翼型邊界層，雷諾數(shù)增大，使流動粘性作用相對減弱，會增大邊界層當(dāng)?shù)販u量雷諾數(shù)，從而使翼型轉(zhuǎn)捩位置前移，減小層流分離泡尺寸。

(2) 分離泡減小使邊界層能量損耗減小，同時轉(zhuǎn)捩位置前移使湍流區(qū)增長，這些使得流動動能更加充沛，流動穩(wěn)定性得到加強，抵抗逆壓梯度干擾的能力增強，流動分離趨勢減弱。

(3) 當(dāng)雷諾數(shù)增大時，翼型后緣分離將延遲，失速迎角及最大升力系數(shù)提高。