張瑞民,于金玲
(中國航天空氣動力技術研究院第二研究所,北京100074)
目前,在飛行器研制過程中,關于升力的計算精度很高,能夠達到實際要求,但阻力的計算精度相對較低。尤其對大展弦比機翼和超臨界機翼來說,摩擦阻力占整個阻力的主要部分,而影響摩擦阻力特性的關鍵因素是邊界層轉捩和分離。因此,準確預測邊界層的轉捩點位置對于提高飛行器的氣動特性預測精度具有重要意義。
由于對邊界層轉捩的物理機理缺乏認識,能有效預測邊界層轉捩的理論還不成熟。目前,工程上常用的轉捩預測方法主要有[1-2]:(1)eN方法[3-4];(2)經驗關系式;(3)低雷諾數模型方法。eN方法以線性穩定性理論為基礎,同時與風洞試驗相結合,能夠預測低雷諾數下二維翼型的自然轉捩和分離流轉捩,卻難以模擬bypass轉捩,且難以推廣到三維流動;經驗關系式則基于大量的試驗數據,借助適當的關系式可以預測自然轉捩和bypass轉捩等,但因該方法引入邊界層動量厚度,致使其與CFD方法不相容[5];低雷諾數湍流模型是從N-S方程推導出來的流場中脈動量的描述方程,具有模擬脈動量發展、捕捉層流到湍流流態轉變的能力,不僅能夠反映轉捩現象,而且僅涉及當地變量,與現代CFD技術完美相容。
本文以k-ωSST兩方程湍流模型為研究對象,引入間歇函數對其所帶的Wilcox轉捩模式進行了修正,進而對傳統湍流翼型的氣動特性進行了研究。
流場計算采用基于壓力和速度耦合的SIMPLE算法來求解定常可壓縮的N-S方程,即:

式中,ρA為氣流的密度;U為氣流速度矢量;Γφ為擴散系數;Sφ為源項;φ為通用項;邊界條件選用壓力遠場和無滑移壁面邊界條件;方程的離散采用有限體積法;為了保證計算求解穩定和收斂,所有求解方程的對流項采用二階迎風格式,擴散項采用中心差分格式。
湍流模型選用k-ωSST兩方程模型,湍流運動方程如下:

式中,ρ為空氣密度;Gk,Gω為k和ω的生成項;Yk,Yω為由湍流引起的k和ω的耗散項;Γk,Γω為k和ω的有效擴散系數,即:

式中,μ為層流黏性系數;μt為湍流黏性系數,表達式為:

式中,S為應變率幅值;α1為常數;F2為混合函數。在高雷諾數下,間歇函數α*=1,此時流動為湍流,如果考慮低雷諾數轉捩的影響,Wilcox的間歇函數表達式為:

式中,α*0=0.024;Ret=ρk/(μω);Rk=6。
由于基于k-ωSST兩方程湍流模型的原始Wilcox轉捩模式本身對擾動過于敏感,使得邊界層轉捩位置比實際情況明顯靠前,致使邊界層流動很快從轉捩區域進入到全湍流區域,進而影響了翼型的流場特性和氣動性能。因此,本文從轉捩流動的物理特征出發,針對湍流黏性系數中的間歇函數進行調整,使其既能模擬全湍流流動,又能模擬邊界層轉捩的流動特性。由于間歇函數僅僅是湍流雷諾數Ret的函數,因此只對 Ret進行修正,修正后的Ret為[6]:

本文分別在未考慮轉捩和修正轉捩模式的情況下,研究了NACA0012翼型的流場特性和氣動性能。翼型周圍網格以及局部放大圖如圖1和圖2所示,壁面第一層網格滿足y+≤1.0,計算域四周邊界距翼型表面均為弦長的20倍。

圖1 NACA0012翼型的網格劃分

圖2 翼型網格的局部放大圖
計算條件為:Ma=0.3;Re=3 × 106;ω =k/μt,湍流參數取值為 k=(0.001×U∞)2,μt=0.009μ。
圖3和圖4分別給出了在原始的Wilcox轉捩和修正的Wilcox轉捩模式下翼型的升力系數和阻力系數隨迎角的變化關系,并與文獻[4]中的試驗值進行了比較。由圖可知,與原始的Wilcox轉捩模式相比,基于修正后的Wilcox轉捩模式的阻力預測結果與試驗值相比更接近,但仍有較大的差距。
圖5給出了翼型在不同迎角下基于原始的Wilcox轉捩模式和修正的Wilcox轉捩模式的邊界層轉捩位置,并與文獻[7]中的試驗結果進行了比較。其中,xtr/c為翼型上表面的轉捩位置與弦長之比,α為氣流迎角。從圖中可以看出,基于修正的轉捩模式,對阻力的預測精度有了一定程度的提高。

圖3 升力系數隨迎角的變化關系

圖4 阻力系數隨迎角的變化關系

圖5 不同迎角下翼型上、下表面的轉捩位置
圖6和圖7分別給出了在0°和4°迎角下轉捩模式修正前后翼型上表面壓力分布的計算結果。通過比較可以看出,在不同的Wilcox轉捩模式下,翼型表面的壓力分布基本一致。圖8和圖9分別給出了在0°和4°迎角下轉捩模式修正前后翼型上表面摩擦阻力分布的計算結果。與0°迎角相比,迎角為4°時的翼型上表面的順壓梯度區減小,逆壓梯度區增加,導致邊界層流動提前轉捩,這與試驗結果相符;此外通過比較可以看出,當迎角為0°和4°時,基于本文的轉捩模式修正,翼型表面的轉捩位置預測結果與試驗值更加接近。

圖6 翼型上表面壓力系數分布(α=0°)

圖7 翼型上表面壓力系數分布(α=4°)

圖8 翼型上表面摩擦系數分布(α=0°)

圖9 翼型上表面摩擦系數分布(α=4°)
邊界層轉捩是決定翼型氣動特性的重要因素。本文從轉捩流動的物理特征出發,引入間歇函數對SST湍流模型的Wilcox轉捩模式進行了修正,并選取了典型的有壓力梯度的NACA0012翼型進行了計算和驗證。研究表明,改進后的模型對轉捩位置具有更好的預測能力;在采用修正邊界層轉捩模型的情況下,翼型的阻力預測精度有了一定程度的提高,但與試驗結果相比,仍存在差距,有待進一步改進。參考文獻:
[1]徐星仲.轉捩流動的數值方法研究[D].北京:中國科學院工程熱物理研究所,1996.
[2]是勛剛.湍流[M].天津:天津大學出版社,1994:136-138.
[3]Smith A M O,Gamberoni N.Transition,pressure gradient and stability theory[R].USA:Douglas Aircraft Company,1956.
[4]Van Ingen JL.A suggested semi-empirical method for the calculation of the boundary layer transition region[R].Delft,Holland:University of Technical Aerospace Engineering,Report VTH-74,1956.
[5]Langtry R B,Menter F R.Transition modeling for general CFD applications in aeronautics[R].AIAA-2005-0522,2005.
[6]錢煒祺,詹浩.一種基于湍流模式的轉捩預測方法[J].空氣動力學學報,2006,24(4):502-507.
[7]Johansen J,Sorensen JN.Prediction of laminar/turbulent transition in airfoil flow[R].AIAA-98-0702,1998.