張瑞民，于金玲

(中國(guó)航天空氣動(dòng)力技術(shù)研究院第二研究所，北京100074)

0 引言

目前，在飛行器研制過(guò)程中，關(guān)于升力的計(jì)算精度很高，能夠達(dá)到實(shí)際要求，但阻力的計(jì)算精度相對(duì)較低。尤其對(duì)大展弦比機(jī)翼和超臨界機(jī)翼來(lái)說(shuō)，摩擦阻力占整個(gè)阻力的主要部分，而影響摩擦阻力特性的關(guān)鍵因素是邊界層轉(zhuǎn)捩和分離。因此，準(zhǔn)確預(yù)測(cè)邊界層的轉(zhuǎn)捩點(diǎn)位置對(duì)于提高飛行器的氣動(dòng)特性預(yù)測(cè)精度具有重要意義。

由于對(duì)邊界層轉(zhuǎn)捩的物理機(jī)理缺乏認(rèn)識(shí)，能有效預(yù)測(cè)邊界層轉(zhuǎn)捩的理論還不成熟。目前，工程上常用的轉(zhuǎn)捩預(yù)測(cè)方法主要有［1-2］:(1)eN方法［3-4］;(2)經(jīng)驗(yàn)關(guān)系式;(3)低雷諾數(shù)模型方法。eN方法以線性穩(wěn)定性理論為基礎(chǔ)，同時(shí)與風(fēng)洞試驗(yàn)相結(jié)合，能夠預(yù)測(cè)低雷諾數(shù)下二維翼型的自然轉(zhuǎn)捩和分離流轉(zhuǎn)捩，卻難以模擬bypass轉(zhuǎn)捩，且難以推廣到三維流動(dòng);經(jīng)驗(yàn)關(guān)系式則基于大量的試驗(yàn)數(shù)據(jù)，借助適當(dāng)?shù)年P(guān)系式可以預(yù)測(cè)自然轉(zhuǎn)捩和bypass轉(zhuǎn)捩等，但因該方法引入邊界層動(dòng)量厚度，致使其與CFD方法不相容［5］;低雷諾數(shù)湍流模型是從N-S方程推導(dǎo)出來(lái)的流場(chǎng)中脈動(dòng)量的描述方程，具有模擬脈動(dòng)量發(fā)展、捕捉層流到湍流流態(tài)轉(zhuǎn)變的能力，不僅能夠反映轉(zhuǎn)捩現(xiàn)象，而且僅涉及當(dāng)?shù)刈兞?，與現(xiàn)代CFD技術(shù)完美相容。

本文以k-ωSST兩方程湍流模型為研究對(duì)象，引入間歇函數(shù)對(duì)其所帶的Wilcox轉(zhuǎn)捩模式進(jìn)行了修正，進(jìn)而對(duì)傳統(tǒng)湍流翼型的氣動(dòng)特性進(jìn)行了研究。

1 數(shù)值方法

1.1 控制方程

流場(chǎng)計(jì)算采用基于壓力和速度耦合的SIMPLE算法來(lái)求解定?？蓧嚎s的N-S方程，即:

式中，ρA為氣流的密度;U為氣流速度矢量;Γφ為擴(kuò)散系數(shù);Sφ為源項(xiàng);φ為通用項(xiàng);邊界條件選用壓力遠(yuǎn)場(chǎng)和無(wú)滑移壁面邊界條件;方程的離散采用有限體積法;為了保證計(jì)算求解穩(wěn)定和收斂，所有求解方程的對(duì)流項(xiàng)采用二階迎風(fēng)格式，擴(kuò)散項(xiàng)采用中心差分格式。

1.2 湍流模型

湍流模型選用k-ωSST兩方程模型，湍流運(yùn)動(dòng)方程如下:

式中，ρ為空氣密度;Gk，Gω為k和ω的生成項(xiàng);Yk，Yω為由湍流引起的k和ω的耗散項(xiàng);Γk，Γω為k和ω的有效擴(kuò)散系數(shù)，即:

式中，μ為層流黏性系數(shù);μt為湍流黏性系數(shù)，表達(dá)式為:

式中，S為應(yīng)變率幅值;α1為常數(shù);F2為混合函數(shù)。在高雷諾數(shù)下，間歇函數(shù)α*=1，此時(shí)流動(dòng)為湍流，如果考慮低雷諾數(shù)轉(zhuǎn)捩的影響，Wilcox的間歇函數(shù)表達(dá)式為:

式中，α*0=0.024;Ret=ρk/(μω);Rk=6。

1.3 轉(zhuǎn)捩模式修正

由于基于k-ωSST兩方程湍流模型的原始Wilcox轉(zhuǎn)捩模式本身對(duì)擾動(dòng)過(guò)于敏感，使得邊界層轉(zhuǎn)捩位置比實(shí)際情況明顯靠前，致使邊界層流動(dòng)很快從轉(zhuǎn)捩區(qū)域進(jìn)入到全湍流區(qū)域，進(jìn)而影響了翼型的流場(chǎng)特性和氣動(dòng)性能。因此，本文從轉(zhuǎn)捩流動(dòng)的物理特征出發(fā)，針對(duì)湍流黏性系數(shù)中的間歇函數(shù)進(jìn)行調(diào)整，使其既能模擬全湍流流動(dòng)，又能模擬邊界層轉(zhuǎn)捩的流動(dòng)特性。由于間歇函數(shù)僅僅是湍流雷諾數(shù)Ret的函數(shù)，因此只對(duì) Ret進(jìn)行修正，修正后的Ret為［6］:

2 計(jì)算結(jié)果與分析

本文分別在未考慮轉(zhuǎn)捩和修正轉(zhuǎn)捩模式的情況下，研究了NACA0012翼型的流場(chǎng)特性和氣動(dòng)性能。翼型周?chē)W(wǎng)格以及局部放大圖如圖1和圖2所示，壁面第一層網(wǎng)格滿足y+≤1.0，計(jì)算域四周邊界距翼型表面均為弦長(zhǎng)的20倍。

圖1 NACA0012翼型的網(wǎng)格劃分

圖2 翼型網(wǎng)格的局部放大圖

計(jì)算條件為:Ma=0.3;Re=3 × 106;ω =k/μt，湍流參數(shù)取值為 k=(0.001×U∞)2，μt=0.009μ。

圖3和圖4分別給出了在原始的Wilcox轉(zhuǎn)捩和修正的Wilcox轉(zhuǎn)捩模式下翼型的升力系數(shù)和阻力系數(shù)隨迎角的變化關(guān)系，并與文獻(xiàn)［4］中的試驗(yàn)值進(jìn)行了比較。由圖可知，與原始的Wilcox轉(zhuǎn)捩模式相比，基于修正后的Wilcox轉(zhuǎn)捩模式的阻力預(yù)測(cè)結(jié)果與試驗(yàn)值相比更接近，但仍有較大的差距。

圖5給出了翼型在不同迎角下基于原始的Wilcox轉(zhuǎn)捩模式和修正的Wilcox轉(zhuǎn)捩模式的邊界層轉(zhuǎn)捩位置，并與文獻(xiàn)［7］中的試驗(yàn)結(jié)果進(jìn)行了比較。其中，xtr/c為翼型上表面的轉(zhuǎn)捩位置與弦長(zhǎng)之比，α為氣流迎角。從圖中可以看出，基于修正的轉(zhuǎn)捩模式，對(duì)阻力的預(yù)測(cè)精度有了一定程度的提高。

圖3 升力系數(shù)隨迎角的變化關(guān)系

圖4 阻力系數(shù)隨迎角的變化關(guān)系

圖5 不同迎角下翼型上、下表面的轉(zhuǎn)捩位置

圖6和圖7分別給出了在0°和4°迎角下轉(zhuǎn)捩模式修正前后翼型上表面壓力分布的計(jì)算結(jié)果。通過(guò)比較可以看出，在不同的Wilcox轉(zhuǎn)捩模式下，翼型表面的壓力分布基本一致。圖8和圖9分別給出了在0°和4°迎角下轉(zhuǎn)捩模式修正前后翼型上表面摩擦阻力分布的計(jì)算結(jié)果。與0°迎角相比，迎角為4°時(shí)的翼型上表面的順壓梯度區(qū)減小，逆壓梯度區(qū)增加，導(dǎo)致邊界層流動(dòng)提前轉(zhuǎn)捩，這與試驗(yàn)結(jié)果相符;此外通過(guò)比較可以看出，當(dāng)迎角為0°和4°時(shí)，基于本文的轉(zhuǎn)捩模式修正，翼型表面的轉(zhuǎn)捩位置預(yù)測(cè)結(jié)果與試驗(yàn)值更加接近。

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

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

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

圖9 翼型上表面摩擦系數(shù)分布(α=4°)

3 結(jié)束語(yǔ)

邊界層轉(zhuǎn)捩是決定翼型氣動(dòng)特性的重要因素。本文從轉(zhuǎn)捩流動(dòng)的物理特征出發(fā)，引入間歇函數(shù)對(duì)SST湍流模型的Wilcox轉(zhuǎn)捩模式進(jìn)行了修正，并選取了典型的有壓力梯度的NACA0012翼型進(jìn)行了計(jì)算和驗(yàn)證。研究表明，改進(jìn)后的模型對(duì)轉(zhuǎn)捩位置具有更好的預(yù)測(cè)能力;在采用修正邊界層轉(zhuǎn)捩模型的情況下，翼型的阻力預(yù)測(cè)精度有了一定程度的提高，但與試驗(yàn)結(jié)果相比，仍存在差距，有待進(jìn)一步改進(jìn)。參考文獻(xiàn):

［1］徐星仲．轉(zhuǎn)捩流動(dòng)的數(shù)值方法研究［D］．北京:中國(guó)科學(xué)院工程熱物理研究所，1996．

［2］是勛剛．湍流［M］．天津:天津大學(xué)出版社，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］錢(qián)煒祺，詹浩．一種基于湍流模式的轉(zhuǎn)捩預(yù)測(cè)方法［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2006，24(4):502-507．

［7］Johansen J，Sorensen JN．Prediction of laminar/turbulent transition in airfoil flow［R］．AIAA-98-0702，1998．