張彥軍,趙軻,張同鑫,陳利麗
(第一飛機設計研究院 總體氣動設計研究所,西安 710089)
在航空航天領域,失速問題對飛行器設計和使用造成惡劣影響,甚至導致機毀人亡。為了滿足飛行安全和改善起降性能,失速問題一直是氣動設計師必須面對和解決的基礎科學問題。探索飛行器的失速問題隱含的空氣動力學機理是明確失速機制及建立控制措施的必由之路,但是復雜的部件干擾和雷諾數變化增加了研究失速機理的難度。幸運的是,雖然難以給出飛行器失速問題精確的解析關系式,但組成這種強非線性現象的核心機制仍然是翼型的失速問題,因此研究翼型的失速問題是建立飛行器失速問題理論體系的基石。
從20世紀上半葉開始,學者們采用試驗方法對翼型失速的表面氣動力和邊界層變化進行了諸多分析,E.N.Jacobs[1]第一次推測:升力線的拐點意味著翼型表面出現了分離現象,即提出失速問題;隨后,E.N.Jacobs等[2]第一次提出雷諾數與翼型分離之間存在關系,重點關注層流轉捩現象對翼型前緣和后緣分離的影響;Lissaman[3]闡述低速翼型與雷諾數之間的影響作用,并分析雷諾數影響下分離泡的延伸與破裂是導致翼型失速的主要原因;V.O.Michael等[4]基于PIV技術研究低雷諾數翼型分離泡的相關特性;Hu H等[5]使用PIV檢測GA(W)-1翼型流場和表面壓力分布,分析層流分離泡與壓力梯度的相互影響;H.Yamato等[6]使用動態(tài)控制裝置影響前緣分離泡破碎,從而抑制翼型分離。國內, Wang N等[7]通過試驗研究雷諾數對RAE2822翼型氣動特性的影響;孫智偉等[8]在某NF-3和NF-6風洞中進行超臨界翼型的設計和實驗研究,給出超臨界翼型在高、低速時的分離演化特性。
在數值模擬方面,A.C.Aranake等[9]使用SA湍流模型耦合γ-Reθt轉捩模型對S827和S809翼型分別進行失速特性數值模擬,相比于全湍流計算結果,無論是升力特性還是阻力特性均與試驗結果更加接近,且考慮邊界層轉捩時計算所得氣動力系數能夠反應失速狀態(tài)下的輕失速和深失速特性;Ke J等[10]采用耦合轉捩模型的RANS方法和混合方法對A-Airfoil及NACA0012翼型靜態(tài)失速特性進行模擬,分析前緣層流分離泡及后緣分離發(fā)展,其結果與試驗結果吻合良好。國內,馮濤等[11]針對典型氣動負載渦輪葉片在不同雷諾數下邊界層特性變化進行敏感性研究;靳允立等[12]采用不同湍流模型研究翼型失速及風速變化引起的雷諾數改變對翼型和風力機氣動性能的影響;王強等[13]采用基于SST模型的γ-Reθt轉捩模型研究S809翼型動態(tài)失速過程中分離泡的作用;袁尚科[14]通過理論分析、數值計算及外場實驗相結合的方法,分別對風力機常用翼型及葉片的靜態(tài)失速、動態(tài)失速及三維旋轉效應等進行研究。
綜上所述,國內外諸多學者將翼型失速的研究對象重點集中于失速點的精確預測,對于失速時邊界層的變化,層流轉捩現象的影響作用和雷諾數對邊界層演化的完整影響過程并未有過多涉足。
本文使用Menter SST湍流模型和擾動放大因子輸運模型對轉捩翼型邊界層進行數值模擬,分析雷諾數變化時邊界層的發(fā)展,總結雷諾數變化對翼型失速的影響規(guī)律。
流場模擬采用Navier-Stokes方程組進行計算,其守恒積分形式為
(1)
式中:Q為守恒向量;S為控制體V所有表面S的集合;F為通過控制體的S面上的通量;n為S面的單位外法向矢量。無粘項、粘性項的空間離散方法分別為二階迎風Roe格式及二階中心差分方法。時間離散則采用近似因子分解方法。
Menterk-ωSST湍流模型[15]是一個兩方程k-ω渦粘性模型,在工程界廣泛應用。SST模型是一個耦合模型,在近壁面保持Wilcoxk-ω模型[16]的特性,而在湍流邊界層外及自由來流中保持Launder-Sharmak-ε模型[17]的特性。通過混合函數保留兩個模型的優(yōu)點并克服模型各自的不足。其方程形式為
(2)
(3)
渦粘性由下式計算得到
(4)
公式中模型參數(α,β,σ,F1,F2)的數值及處理方法可以參照文獻[15]。
在模擬層流-湍流轉捩邊界層時,采用Coder等[18]提出的擾動放大因子輸運模型。該輸運方程模型形式如式(5)所示。
(5)
式中的Fcrit,Fgrowth,dn/dRδ2均與邊界層形狀因子相關,關系式的具體表示可參考文獻[18]。因此模型的核心點為使用層流邊界層的當地壓力梯度參數HL計算非當地參數形狀因子H12,后者常用于描述速度型的特征。當地壓力梯度參數定義為
(6)
式中:d為壁面法向距離;S為當地剪切應變率;Ue為邊界層邊緣速度。
由Falkner-Skan二維邊界層相似性解可以得到不同壓力梯度下HL最大值與形狀因子的函數關系,即
(7)
從而實現在流場中形狀因子H12的求解,邊界層中max(H12)為精確值。
將擾動放大因子輸運方程與SST模型耦合,其湍動能輸運方程為
(8)
(9)
ft2=[1-exp(2(n-Ncrit))]exp(ct4χ2)
(10)
式中:模型參數F1,ct4,ct5,χ,Ncrit等的設置與文獻[19]一致。
使用NLF(1)-0416自然層流翼型繞流[20]研究雷諾數的影響。該翼型有較長的層流區(qū)且具備高升力特性,為典型后緣失速翼型,一般應用于通用航空飛機翼型和風力機葉型等領域。對其進行轉捩數值模擬,計算網格分布如圖 1所示,計算域半徑為100倍弦長,翼型表面布置481個網格單元,法向布置181個網格單元,網格量約為18萬。壁面法向第一層網格高度約為1.0×10-6m,使y+小于1。風洞試驗雷諾數工況分別為2×106,4×106,9×106。風洞自由來流馬赫數為0.1,湍流度設置為0.03%。
圖1 翼型繞流數值模擬網格分布
為了與試驗狀態(tài)保持一致,便于對比驗證和分析,本文選擇了與試驗雷諾數相同的三個雷諾數進行數值模擬研究,鑒于雷諾數對失速特性的影響一直是一個復雜的問題,對于更大雷諾數,我們將在后續(xù)工作中開展研究。使用擾動放大因子輸運模型對三個雷諾數下翼型繞流進行模擬,得到升力系數-阻力系數極曲線并與全湍數據和試驗數據進行對比,結果如圖2所示。
(a) Re=2×106
(b) Re=4×106
(c) Re=9×106
從圖2可以看出:相比全湍模型的數值模擬結果,轉捩模型的預測結果與試驗數據更加貼合,與真實流動更加吻合。
雷諾數為4×106,0°迎角工況的翼型表面壓力系數分布如圖3所示。
圖3 Re=4.0×106,0°迎角翼型表面的壓力系數分布
從圖3可以看出:在該工況下,翼型繞流轉捩類型為層流分離泡轉捩,且轉捩模式對該分離再附著過程描述較好。模型對轉捩過程的準確預測為分析邊界層發(fā)展奠定基礎。
三個雷諾數下翼型升力系數隨迎角變化曲線如圖4所示,可以看出:雷諾數對翼型的氣動特性有較大的影響,當雷諾數增大時,翼型的失速迎角變大,同時最大升力系數增大。產生此現象的主要原因是雷諾數變化對翼型邊界層轉捩及發(fā)展有重要影響,因此通過對比不同雷諾數下的轉捩過程及邊界層發(fā)展情況,對雷諾數的影響規(guī)律進行分析。
圖4 不同雷諾數下翼型升力系數隨迎角變化曲線
不同雷諾數下翼型上下表面的轉捩位置隨升力系數變化曲線如圖5所示。
圖5 不同雷諾數翼型轉捩位置隨升力系數的變化
從圖5可以看出:隨著迎角增大,升力系數增大,翼型上表面的轉捩位置逐漸前移,而下表面變化趨勢相反。相同升力系數下,當雷諾數增大時,上、下表面的轉捩位置均逐漸前移。
對0°迎角工況下翼型邊界層參數隨雷諾數的變化進行分析。0°迎角時翼型的表面壓力分布圖如圖6所示,可以看出:雷諾數變化對壓力系數的整體影響較小,其主要影響體現在壓力系數突變位置的變化,即出現層流分離泡、發(fā)生轉捩的位置。
圖6 0°迎角時翼型表面壓力系數分布對比
不同雷諾數在0°迎角時翼型的表面摩擦力系數分布如圖7所示,可以看出:在雷諾數為2×106和4×106的工況下,翼型上下表面的摩擦力系數均出現負值,表面附面層中出現層流分離泡,并伴隨轉捩的發(fā)生;而雷諾數為9×106時,翼型繞流發(fā)生層流到湍流的轉捩過程,但是并沒有層流分離泡出現。初步的結論與傳統(tǒng)的研究結論一致,低雷諾數環(huán)境雖然轉捩位置比較靠后,但是層流非常容易被逆壓梯度誘導產生層流分離泡,觸發(fā)轉捩。該現象在高空長航時無人機表面比較常見,因為數萬米高空中空氣相對稀薄,雷諾數較小,容易產生層流分離泡。
(a) Re=2×106
(b) Re=4×106
(c) Re=9×106
翼型上、下表面層流分離泡形態(tài)如圖8所示,可以看出:當雷諾數增大時,翼型上下表面的層流分離泡均逐漸減小,在雷諾數為9×106時,層流分離泡已經消失。
(a1) 上表面分離泡
(a2) 下表面分離泡
(b1) 上表面分離泡
(b2) 下表面分離泡
(c1) 上表面分離泡
(c2) 下表面分離泡
對比不同雷諾數下翼型上表面邊界層的邊界層動量損失厚度δ2、形狀因子H12、渦量雷諾數Reν等參數分布,如圖9~圖11所示。
圖9 0°迎角上表面邊界層動量損失厚度分布
圖10 0°迎角翼型上表面邊界層形狀因子分布
圖11 翼型上表面邊界層渦量雷諾數分布
從圖9可以看出:當雷諾數減小時,動量損失厚度增大,這是由于雷諾數減小使流動粘性作用相對增強,有利于邊界層的發(fā)展。
邊界層形狀因子的定義為
(11)
式中:δ1為邊界層位移厚度。
從圖10可以看出:當層流邊界層通過轉捩點變?yōu)橥牧鬟吔鐚樱吔鐚雍穸葧蝗辉龊?。湍流邊界層由于流速分布更趨均勻化而使邊界層位移厚度減小,同時由于阻力增加而使邊界層動量損失厚度增大,于是邊界層形狀因子在轉捩后將減小。翼型上表面在三個雷諾數下的轉捩位置均在40%弦長附近,因此上表面邊界層形狀因子在40%弦長處出現峰值,而且由于雷諾數為2×106和4×106時翼型出現層流分離泡,所以這兩個雷諾數對應的形狀因子峰值更高;而雷諾數為9×106時邊界層未出現分離泡,且轉捩位置前移,40%弦長處流動已基本發(fā)展為湍流,對應的形狀因子值較低。在轉捩發(fā)生后,形狀因子的發(fā)展情況與前文的分析對應:三組形狀因子均迅速下降,雷諾數越大,最終邊界層形狀因子的值越小。
從圖11可以看出:相同弦長站位處,雷諾數越大,邊界層中的渦量雷諾數越高。因為渦量雷諾數與常用轉捩判定參數動量厚度雷諾數幾乎是正比關系,因此可以作為邊界層發(fā)展的標志性參數進行分析。邊界層轉捩觸發(fā)位置隨雷諾數升高而向前緣移動,在轉捩后,渦量雷諾數迅速增大,使邊界層保持在湍流狀態(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可以看出:在翼型前緣,轉捩前較高雷諾數對應邊界層的速度型更加均勻、飽滿,邊界層更薄。在40%弦長處,由于翼型邊界層在較低雷諾數時出現層流分離泡,速度型分布在壁面附近呈凹曲線,速度型出現拐點。轉捩后,三個雷諾數對應的速度型更加相近,高雷諾數的速度型分布更加飽滿。
為了分析雷諾數對失速迎角附近邊界層發(fā)展的影響,取12°迎角工況下翼型邊界層進行觀察。12°迎角時翼型的表面摩擦力系數分布如圖13所示,隨雷諾數增高翼型轉捩點前移,轉捩處層流分離泡逐漸縮小消失,僅雷諾數為2×106時存在比較明顯的層流分離泡。
圖13 12°迎角翼型表面摩擦力系數分布
12°迎角時翼型上表面不同站位流向速度型分布對比如圖14所示。由于雷諾數為2×106時在約7%弦長處產生層流分離泡,其對應的速度分布出現明顯的拐折。雷諾數為9×106時的速度型更早趨于均勻,因其轉捩點更加靠近前緣。轉捩后,隨雷諾數的增高,流速分布變得更加均勻而壁面附近流速梯度增大。
翼型上表面邊界層內各站位處的湍動能分布如圖15所示。從翼型表面摩擦力系數分布可知,三個雷諾數下翼型上表面均已在10%弦長前完成轉捩,因此下圖中的結果均為轉捩后邊界層湍動能的分布。
(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可以看出:轉捩后相同站位低雷諾數對應的湍動能峰值更高。這是由于翼型前緣發(fā)生轉捩,在較低雷諾數時伴隨有層流分離泡形成,這一過程使邊界層流動的脈動特性得到加強,其對應的湍動能、渦粘性增大。但是層流分離泡會消耗邊界層的動能,使流動自持能力即抗分離能力減弱,從而產生分離。當雷諾數增高,分離泡減小甚至消失,流動能量消耗減小,動能更充沛,壓力分布可以維持較長距離的梯度,后緣分離滯后發(fā)生,從而增大翼型升力系數及失速迎角。
最后對比相同迎角不同雷諾數下翼型流場的后緣分離情況,結果如圖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可以看出:在相同迎角下,隨著雷諾數增大,流動穩(wěn)定性增強,抵抗逆壓梯度干擾的能力提高,流動分離趨勢減弱。后緣分離逐漸減小,這與本節(jié)中對邊界層發(fā)展及其參數變化趨勢的觀察結果相一致。
(1) 對于層流-湍流轉捩翼型邊界層,雷諾數增大,使流動粘性作用相對減弱,會增大邊界層當地渦量雷諾數,從而使翼型轉捩位置前移,減小層流分離泡尺寸。
(2) 分離泡減小使邊界層能量損耗減小,同時轉捩位置前移使湍流區(qū)增長,這些使得流動動能更加充沛,流動穩(wěn)定性得到加強,抵抗逆壓梯度干擾的能力增強,流動分離趨勢減弱。
(3) 當雷諾數增大時,翼型后緣分離將延遲,失速迎角及最大升力系數提高。