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

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

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

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

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

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

(2)

式中:γ為間歇因子;Pγ為生成項(xiàng);Eγ為耗散項(xiàng)。若“當(dāng)?shù)剞D(zhuǎn)捩”的判據(jù)滿足，Pγ啟動(dòng)即γ開始增長(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ù) 變化曲線Fig．2 Lift，drag coefficients as a function of Mach number

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

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

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

圖6 考慮轉(zhuǎn)捩影響的NACA64A010翼型力矩系數(shù)隨攻角變化曲線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線性多步法[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)定常升力線斜率隨馬赫數(shù)變化曲線。從圖中結(jié)果可以看出，當(dāng)0.83≤Ma≤0.87時(shí)，隨著馬赫數(shù)增大，翼型表面的分離區(qū)擴(kuò)大，使得升力線斜率急劇減小，導(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)定常升力線斜率隨馬赫數(shù)變化曲線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ì)于俯仰位移的滯后相位角。本文采用類似的分析方法對(duì)Isogai案例A模型在全湍和自由轉(zhuǎn)捩條件下的跨音速穩(wěn)定性進(jìn)行分析，其中k為減縮頻率。

圖16、圖17分別為俯仰運(yùn)動(dòng)中力矩系數(shù)的幅值和相位隨馬赫數(shù)變化曲線。計(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ù)變化曲線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)類型，不再予以展示。從圖17中可以得知，自由轉(zhuǎn)捩與全湍在“凍結(jié)區(qū)域”的相位角均接近-180°，這意味著此時(shí)氣動(dòng)力基本不做功，系統(tǒng)處于臨界穩(wěn)定狀態(tài)，翼型的顫振速度隨馬赫數(shù)變化很小。

圖17 俯仰運(yùn)動(dòng)中力矩系數(shù)的相位隨馬赫數(shù)變化曲線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ì)顫振特性幾乎沒有影響，隨馬赫數(shù)的增大，顫振速度的變化很小。

轉(zhuǎn)捩現(xiàn)象對(duì)激波強(qiáng)度和翼型表面分離區(qū)大小均有影響，使得翼型的氣動(dòng)力非線性較全湍增強(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] 錢煒祺, 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] 葉正寅, 張偉偉, 史愛明, 等. 流固耦合力學(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