李國(guó)俊， 白俊強(qiáng)， 劉 南， 徐家寬， 喬 磊

(西北工業(yè)大學(xué) 航空學(xué)院， 西安 710072)

李國(guó)俊， 白俊強(qiáng)， 劉 南， 徐家寬， 喬 磊

(西北工業(yè)大學(xué) 航空學(xué)院， 西安 710072)

轉(zhuǎn)捩；時(shí)域；顫振；全湍；跨音速

機(jī)翼表面的邊界層轉(zhuǎn)捩現(xiàn)象在現(xiàn)代飛行器的飛行過(guò)程中較為常見(jiàn)，邊界層轉(zhuǎn)捩對(duì)摩擦阻力大小、流動(dòng)的分離位置和跨音速激波位置等有很大的影響，使得氣動(dòng)力的非線(xiàn)性特性和黏性效應(yīng)更為復(fù)雜[1]，進(jìn)而影響機(jī)翼的顫振特性。然而，目前基于RANS(Reynold Averaged Navier-Stockes)方程對(duì)顫振問(wèn)題進(jìn)行研究時(shí)大多數(shù)采用全湍假設(shè)，忽視了邊界層轉(zhuǎn)捩對(duì)顫振特性的影響。如果在飛行器的顫振設(shè)計(jì)校核階段不對(duì)轉(zhuǎn)捩加以考慮，這可能導(dǎo)致設(shè)計(jì)出來(lái)的飛行器的顫振特性過(guò)于保守，難以充分發(fā)揮飛行器的氣動(dòng)性能；或者設(shè)計(jì)出來(lái)的飛行器存在顫振安全隱患，導(dǎo)致結(jié)構(gòu)破壞，嚴(yán)重時(shí)會(huì)造成機(jī)毀人亡。因此，開(kāi)展考慮邊界層轉(zhuǎn)捩影響的顫振問(wèn)題研究具有重要的實(shí)際意義。

1 流場(chǎng)求解

1.1 非定常氣動(dòng)力求解

本文采用課題組自研的CFD代碼-TeAM求解非定常氣動(dòng)力，其控制方程是三維可壓縮非定常積分形式的N-S方程，其直角坐標(biāo)系的守恒形式積分方程為

(1)

式中:Q=[ρρuρvρwe]T為守恒向量;ρ，(u,v,w)和e分別為密度、直角坐標(biāo)系下的速度分量和單位質(zhì)量氣體的總能量;F,G,H為三個(gè)方向的無(wú)黏矢通量;Fv,Gv,Hv為三個(gè)方向的黏性矢通量。其中無(wú)黏項(xiàng)采用Roe格式進(jìn)行離散，黏性項(xiàng)采用二階中心差分進(jìn)行離散，時(shí)間推進(jìn)采用雙時(shí)間步法進(jìn)行迭代求解，為封閉方程引入k-ω SST湍流模型。為了進(jìn)一步提高計(jì)算效率，采用多重網(wǎng)格加速收斂技術(shù)、并行計(jì)算技術(shù)。

轉(zhuǎn)捩的觸發(fā)與轉(zhuǎn)捩區(qū)發(fā)展的預(yù)測(cè)主要在間歇因子輸運(yùn)方程中完成。該方程由Menter等在2004年提出，Langtry等在2006年進(jìn)行了一些改進(jìn)，解決了轉(zhuǎn)捩區(qū)過(guò)短、駐點(diǎn)間歇因子生成項(xiàng)過(guò)大等問(wèn)題。Langtry等在2006年給出的間歇因子輸運(yùn)方程為

(2)

式中:γ為間歇因子;Pγ為生成項(xiàng);Eγ為耗散項(xiàng)。若“當(dāng)?shù)剞D(zhuǎn)捩”的判據(jù)滿(mǎn)足，Pγ啟動(dòng)即γ開(kāi)始增長(zhǎng)，Eγ保證在層流邊界層中γ趨近于0，這也為預(yù)測(cè)再層流化現(xiàn)象提供了條件。

(3)

1.3 轉(zhuǎn)捩預(yù)測(cè)精度驗(yàn)證

圖1 NLR7301翼型網(wǎng)格示意圖Fig．1 Mesh for NLR7301

圖2 升力系數(shù)和阻力系數(shù)隨馬赫數(shù) 變化曲線(xiàn)Fig．2 Lift，drag coefficients as a function of Mach number

圖3 翼型表面轉(zhuǎn)捩位置隨馬赫數(shù)變化曲線(xiàn)Fig．3 Mach number effect on transition onset location

圖4 NACA64A010 翼型網(wǎng)格示意圖Fig．4 Mesh for NACA64A010

圖5 考慮轉(zhuǎn)捩影響的NACA64A010翼型升力系數(shù)隨攻角變化曲線(xiàn)Fig．5 Lift coefficient versus AOA of NACA64A010 airfoil incorporating transition modeling

圖6 考慮轉(zhuǎn)捩影響的NACA64A010翼型力矩系數(shù)隨攻角變化曲線(xiàn)Fig．6 Moment coefficient versus AOA of NACA64A010 airfoil incorporating transition modeling

2 結(jié)構(gòu)運(yùn)動(dòng)方程

具有浮沉和俯仰兩自由度的二維翼型結(jié)構(gòu)運(yùn)動(dòng)方程為

(4)

(5)

式中:m為機(jī)翼質(zhì)量;Sα為二維翼型對(duì)剛心的質(zhì)量靜矩;Dh為浮沉阻尼;Kh為翼型關(guān)于剛心的沉浮剛度;Iα為二維翼型對(duì)剛心的質(zhì)量慣性矩;Dα俯仰阻尼;Kα為翼型關(guān)于剛心的俯仰剛度;L為升力;My為俯仰力矩。

針對(duì)上述二維兩自由度翼型結(jié)構(gòu)運(yùn)動(dòng)方程，經(jīng)無(wú)量綱化可得

(6)

無(wú)量綱顫振速度定義為

(7)

式中:U∞為自由來(lái)流速度；μ=m/πρ∞b2為質(zhì)量比;ρ∞為自由來(lái)流密度。

(8)

本文采用基于預(yù)估-校正技術(shù)的四階隱式Adams線(xiàn)性多步法[15]對(duì)方程(8)進(jìn)行時(shí)域推進(jìn)求解

(9)

該方法既保證了方程的求解效率，又具有較好的魯棒性。

3 全湍流動(dòng)下的跨音速顫振邊界驗(yàn)證

跨音速顫振凹坑[16-17]與空氣的可壓縮性和激波相位滯后效應(yīng)密切相關(guān)，因此跨音速顫振邊界預(yù)測(cè)的關(guān)鍵在于激波捕捉的準(zhǔn)確性。本文首先對(duì)Isogai案例A模型在全湍流動(dòng)下的跨音速顫振邊界進(jìn)行預(yù)測(cè)，以驗(yàn)證本文顫振計(jì)算方法的可靠性。其中全湍計(jì)算采用k-ω SST湍流模型，結(jié)構(gòu)參數(shù)采用文獻(xiàn)[11]中的參數(shù)，雷諾數(shù)按照文獻(xiàn)[12]中對(duì)所有馬赫數(shù)給定為6×106。計(jì)算網(wǎng)格為上述非定常氣動(dòng)力驗(yàn)證所用網(wǎng)格，如圖4所示。

圖7、圖8分別為采用k-ω SST湍流模型計(jì)算得到的全湍條件下的顫振速度邊界和顫振頻率比邊界，與文獻(xiàn)[18]的結(jié)果吻合很好。數(shù)值結(jié)果表明，跨音速凹坑最低點(diǎn)的Ma≈0.83，隨后顫振速度急劇增大；在Ma=0.87時(shí)顫振速度再次減小，在Ma=0.9時(shí)出現(xiàn)第二個(gè)凹坑，隨后顫振速度繼續(xù)增大。

跨音速凹坑最低點(diǎn)后顫振速度的急劇增大與翼型表面分離流動(dòng)的出現(xiàn)和擴(kuò)大密切相關(guān)。圖9展示了四個(gè)不同馬赫數(shù)下的表面摩阻分布對(duì)比，圖10展示了顫振邊界和準(zhǔn)定常升力線(xiàn)斜率隨馬赫數(shù)變化曲線(xiàn)。從圖中結(jié)果可以看出，當(dāng)0.83≤Ma≤0.87時(shí)，隨著馬赫數(shù)增大，翼型表面的分離區(qū)擴(kuò)大，使得升力線(xiàn)斜率急劇減小，導(dǎo)致顫振邊界急劇增大。

圖7 顫振速度邊界(全湍)Fig.7 Flutter speed index boundary(fully turbulent)

圖8 顫振頻率比邊界(全湍)Fig.8 Frequency ratio boundary(fully turbulent)

圖9 不同馬赫數(shù)下的表面摩阻分布對(duì)比Fig.9 Comparison of skin friction distribution between different Mach number

圖10 顫振邊界和準(zhǔn)定常升力線(xiàn)斜率隨馬赫數(shù)變化曲線(xiàn)Fig.10 Flutter boundary and quasi-steady lift curve slope versus Mach number

4 考慮轉(zhuǎn)捩影響的跨音速顫振邊界預(yù)測(cè)

為了探究轉(zhuǎn)捩引起跨音速顫振凹坑加深的原因，對(duì)跨音速凹坑最低點(diǎn)附近翼型的定常壓力分布和表面摩擦阻力分布進(jìn)行對(duì)比分析，如圖13～圖15所示。結(jié)果表明自由轉(zhuǎn)捩下激波位置較全湍靠后，激波強(qiáng)度大于全湍，且自由轉(zhuǎn)捩較全湍流動(dòng)在翼型表面提前發(fā)生分離。激波增強(qiáng)會(huì)降低翼型的穩(wěn)定性，使得翼型的顫振速度降低。結(jié)合以上分析結(jié)果可以得出如下結(jié)論：自由轉(zhuǎn)捩較全湍使得翼型在跨音速凹坑附近具有更低的穩(wěn)定性，從而具有更低的顫振邊界，導(dǎo)致凹坑程度加深。

圖11 顫振速度邊界Fig.11 Flutter speed index boundary

圖12 顫振頻率比邊界Fig.12 Frequency ratio boundary

圖13 全湍和自由轉(zhuǎn)捩壓力及表面摩阻分布對(duì)比( Ma = 0． 82)Fig.13 Comparison of pressure and skin friction distribution between fully turbulent and free transition( Ma = 0． 82)

圖14 全湍和自由轉(zhuǎn)捩壓力及表面摩阻分布對(duì)比( Ma = 0． 83)Fig.14 Comparison of pressure and skin friction distribution between fully turbulent and free transition( Ma = 0． 83)

圖15 全湍和自由轉(zhuǎn)捩壓力及表面摩阻分布對(duì)比( Ma = 0． 84)Fig.12 Comparison of pressure and skin friction distribution between fully turbulent and free transition( Ma = 0． 84)

Bendiksen針對(duì)跨音速顫振提出了“跨音速穩(wěn)定性法則[19]”(Transonic Stabilization Laws)，從氣動(dòng)力做功的角度研究了氣動(dòng)力的幅值及相位和顫振穩(wěn)定性之間的關(guān)系，此處的相位指的是機(jī)翼俯仰力矩相對(duì)于俯仰位移的滯后相位角。本文采用類(lèi)似的分析方法對(duì)Isogai案例A模型在全湍和自由轉(zhuǎn)捩條件下的跨音速穩(wěn)定性進(jìn)行分析，其中k為減縮頻率。

圖16、圖17分別為俯仰運(yùn)動(dòng)中力矩系數(shù)的幅值和相位隨馬赫數(shù)變化曲線(xiàn)。計(jì)算結(jié)果表明，轉(zhuǎn)捩較全湍使得翼型在跨音速凹坑附近具有更大的力矩系數(shù)幅值和更大的超前相位，這意味著此時(shí)氣動(dòng)力對(duì)翼型的做功增大，使得系統(tǒng)的穩(wěn)定性降低，從而導(dǎo)致跨音速凹坑程度加深。值得注意的是，當(dāng)Ma=0.86時(shí)，轉(zhuǎn)捩的相位角表明系統(tǒng)處于不穩(wěn)定狀態(tài)，氣動(dòng)力對(duì)系統(tǒng)做正功，使得系統(tǒng)的穩(wěn)定性降低；而全湍的相位角表明系統(tǒng)處于穩(wěn)定狀態(tài)，氣動(dòng)力對(duì)系統(tǒng)做負(fù)功，使得系統(tǒng)的穩(wěn)定性增大。因此當(dāng)Ma=0.86時(shí)，轉(zhuǎn)捩較全湍具有更小的顫振速度，這也使得轉(zhuǎn)捩的跨音速凹坑范圍較全湍有所擴(kuò)大。

圖16 俯仰運(yùn)動(dòng)中力矩系數(shù)的幅值隨馬赫數(shù)變化曲線(xiàn)Fig.16 Amplitude of moment coefficient versus Mach numbers for pitching motion

圖18為定常流動(dòng)0度攻角下自由轉(zhuǎn)捩和全湍的激波位置對(duì)比圖，當(dāng)Ma≥0.91時(shí)，激波到達(dá)翼型后緣，流場(chǎng)進(jìn)入“凍結(jié)區(qū)域[20]”(Freeze region)，全湍下的空間馬赫數(shù)云圖如圖19～圖22所示。由于此時(shí)自由轉(zhuǎn)捩與全湍的翼型流場(chǎng)類(lèi)型，不再予以展示。從圖17中可以得知，自由轉(zhuǎn)捩與全湍在“凍結(jié)區(qū)域”的相位角均接近-180°，這意味著此時(shí)氣動(dòng)力基本不做功，系統(tǒng)處于臨界穩(wěn)定狀態(tài)，翼型的顫振速度隨馬赫數(shù)變化很小。

圖17 俯仰運(yùn)動(dòng)中力矩系數(shù)的相位隨馬赫數(shù)變化曲線(xiàn)Fig.17 Phase angle of moment coefficient versus Mach numbers for pitching motion

圖18 激波位置Fig.18 The location of shock

圖19 馬赫數(shù)云圖( Ma = 0． 91)Fig.19 Mach contour( Ma = 0． 91)

圖20 馬赫數(shù)云圖( Ma = 0． 93)Fig.20 Mach contour( Ma = 0． 93)

圖21 馬赫數(shù)云圖( Ma = 0． 95)Fig.21 Mach contour( Ma = 0． 95)

圖22 馬赫數(shù)云圖( Ma = 0． 98)Fig.22 Mach contour( Ma = 0． 98)

5 結(jié) 論

(1)本文采用全湍計(jì)算得到的顫振邊界和文獻(xiàn)的結(jié)果吻合很好，驗(yàn)證了本文建立的方法的可靠性。

(2)轉(zhuǎn)捩現(xiàn)象在跨音速范圍內(nèi)使得翼型表面激波強(qiáng)度較全湍增強(qiáng)，氣動(dòng)力對(duì)翼型做功增加，降低了翼型的顫振穩(wěn)定性；跨音速凹坑最低點(diǎn)的位置基本保持不變，該點(diǎn)的顫振速度減小了41.6%，且凹坑范圍擴(kuò)大。

(3)當(dāng)激波到達(dá)翼型后緣，進(jìn)入凍結(jié)區(qū)域，氣動(dòng)力對(duì)翼型幾乎不做功，系統(tǒng)處于臨界穩(wěn)定狀態(tài)，顫振速度隨馬赫的增大而變化很小。

綜上所述，當(dāng)馬赫數(shù)小于跨音速凹坑最低點(diǎn)的馬赫數(shù)時(shí)，激波的強(qiáng)度在顫振的影響因素中占據(jù)主導(dǎo)地位；而在跨音速凹坑附近，激波和翼型表面的分離現(xiàn)象共同作用，影響翼型的顫振特性，其中激波增強(qiáng)會(huì)降低系統(tǒng)的穩(wěn)定性，分離區(qū)域的擴(kuò)大會(huì)增加系統(tǒng)的穩(wěn)定性；當(dāng)馬赫數(shù)接近“凍結(jié)區(qū)域”，激波靠近翼型后緣，此時(shí)激波和分離對(duì)顫振特性的影響逐漸減小，翼型的相位角逐漸接近臨界狀態(tài)，系統(tǒng)逐漸趨于穩(wěn)定，翼型的顫振速度增大；進(jìn)入凍結(jié)區(qū)域后，激波到達(dá)翼型后緣，而翼型表面幾乎不存在分離區(qū)，此時(shí)激波和分離對(duì)顫振特性幾乎沒(méi)有影響，隨馬赫數(shù)的增大，顫振速度的變化很小。

轉(zhuǎn)捩現(xiàn)象對(duì)激波強(qiáng)度和翼型表面分離區(qū)大小均有影響，使得翼型的氣動(dòng)力非線(xiàn)性較全湍增強(qiáng)，顫振特性與全湍條件下的有明顯區(qū)別。因此針對(duì)機(jī)翼表面存在轉(zhuǎn)捩現(xiàn)象的飛機(jī)在進(jìn)行顫振設(shè)計(jì)時(shí)，需要對(duì)轉(zhuǎn)捩現(xiàn)象加以考慮，以便得到精確的顫振邊界，確保飛機(jī)的飛行安全。

[ 1 ] LAPOINTE S, DUMAS G. Improved numerical simulations of self-sustained oscillations of a NACA0012 with transition modeling: AIAA-2011-3258[R]. Hawaii: AIAA, 2011.

[ 2 ] SMITH A M O, GAMBERONI N. Transition pressure gradient and stability theory: Calif.Ref.ES26388[R]. Long Beach: Douglas Aircraft Company ,1956.

[ 3 ] VAN INGEN J L. A suggested semi-empirical method for the calculation of the boundary layer transition region: VTH-74[R]. Delft: University of Delft, 1956.

[ 4 ] LANGTRY R B, MENTER F R. Aorrelation-based transition modeling for unstructured parallelized computational fluid dynamics codes[J]. AIAA Journal, 2009, 47(12): 2894-2906.

[ 5 ] MENTER F R, LANGTRY R B, LIKKI S R, et al. A correlation-based transition model using local variables, Part I-model formulation[J]. Journal of Turbomachinery，2004，128(3)： 413-422.

[ 6 ] MENTER F R, LANGTRY R B, LIKKI S R, et al. A correlation-based transition model using local variables, PartⅡ-model formulation[J]. Journal of Turbomachinery，2004，128(3)： 423-434..

[ 7 ] MAI H, HEBLER A. Aeroelasticity of a laminar wing[C]// International Forum on Aeroelasticity and Structural Dynamics. [S.l.]: IFASD, 2011.

[ 8 ] VAN ROOIJ A C L M. Flutter behaviour of a laminar supercritical airfoil-A numerical investigation into the influence of boundary layer transition[D]. Delft: Delft University of Technology, 2012.

[10] 錢(qián)煒祺, RANDOLPH C K L. 考慮轉(zhuǎn)捩影響的翼型動(dòng)態(tài)失速數(shù)值模擬[J]. 空氣動(dòng)力學(xué)學(xué)報(bào)， 2008， 26(1)： 50-55.

QIAN Weiqi, RANDOLPH C K L. Numerical simulation of airfoil dynamic stall incorporating transition modeling[J]. Acta Aerodynamica Sinica， 2008， 26(1)： 50-55.

[11] ALONSO J, JAMESON A. Fully-implicit time-marching aeroelastic solutions: AIAA-94-0056[R]. Washington, D.C.: AIAA, 1994.

[12] ZHANG Zhichao, YANG Shuichi, LIU Feng. Prediction of flutter anf LCO by an euler method on non-moving cartesian grids with boundary-layer corrections: AIAA-2005-833[R]. Reno: AIAA, 2005.

[13] 喬磊. 考慮轉(zhuǎn)捩判定的分離流動(dòng)數(shù)值模擬研究[D]. 西安： 西北工業(yè)大學(xué)，2013.

[14] NLR 7301 Airfoil. Experimental data base for computer program assessment: AGARDAR-138[R]. [S.l.]: Neuilly Sur Seine, 1979.

[15] 葉正寅, 張偉偉, 史愛(ài)明, 等. 流固耦合力學(xué)基礎(chǔ)及其應(yīng)用[M]. 哈爾濱: 哈爾濱工業(yè)大學(xué)出版社, 2010: 171-173.

[16] ISOGAI K. On the transonic-dip mechanism of flutter of a sweptback wing[J]. AIAA Journal, 1979, 17(7): 793-795.

[17] ISOGAI K. Transonic dip mechanism of flutter of a sweptback wing: Part Ⅱ[J]. AIAA Journal, 1981, 19(9): 1240-1242.

[18] TIMME S, BADCOCK K J. Searching for transonic aeroelastic instability using an aerodynamic model hierarchy: AIAA-2010-3048[R]. Orlando: AIAA, 2010.

[19] BENDIKSEN O. Transonic stabilization laws for unsteady aerodynamics and flutter: AIAA-2012-1718[R]. Honolulu: AIAA, 2012.

[20] BENDIKSEN O. Transonic single-degree-of-freedom flutter and natural mode instabilities: AIAA-2013-1593[R]. Boston: AIAA, 2013.

LI Guojun, BAI Junqiang, LIU Nan, XU Jiakuan, QIAO Lei

(School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China)

transition; time domain; flutter; fully turbulent; transonic

國(guó)家“973”計(jì)劃(2014CB744804)

2016-03-14 修改稿收到日期： 2016-09-02

李國(guó)俊 男，碩士生，1992年生

白俊強(qiáng) 男，博士，教授，1971年生

V215.34

A

10.13465/j.cnki.jvs.2017.22.032