基于转捩模型的二维翼型跨音速颤振边界预测

李国俊， 白俊强， 刘 南， 徐家宽， 乔 磊

(西北工业大学 航空学院， 西安 710072)

李国俊， 白俊强， 刘 南， 徐家宽， 乔 磊

(西北工业大学 航空学院， 西安 710072)

转捩；时域；颤振；全湍；跨音速

机翼表面的边界层转捩现象在现代飞行器的飞行过程中较为常见，边界层转捩对摩擦阻力大小、流动的分离位置和跨音速激波位置等有很大的影响，使得气动力的非线性特性和黏性效应更为复杂[1]，进而影响机翼的颤振特性。然而，目前基于RANS(Reynold Averaged Navier-Stockes)方程对颤振问题进行研究时大多数采用全湍假设，忽视了边界层转捩对颤振特性的影响。如果在飞行器的颤振设计校核阶段不对转捩加以考虑，这可能导致设计出来的飞行器的颤振特性过于保守，难以充分发挥飞行器的气动性能；或者设计出来的飞行器存在颤振安全隐患，导致结构破坏，严重时会造成机毁人亡。因此，开展考虑边界层转捩影响的颤振问题研究具有重要的实际意义。

1 流场求解

1.1 非定常气动力求解

本文采用课题组自研的CFD代码-TeAM求解非定常气动力，其控制方程是三维可压缩非定常积分形式的N-S方程，其直角坐标系的守恒形式积分方程为

(1)

式中:Q=[ρρuρvρwe]T为守恒向量;ρ，(u,v,w)和e分别为密度、直角坐标系下的速度分量和单位质量气体的总能量;F,G,H为三个方向的无黏矢通量;Fv,Gv,Hv为三个方向的黏性矢通量。其中无黏项采用Roe格式进行离散，黏性项采用二阶中心差分进行离散，时间推进采用双时间步法进行迭代求解，为封闭方程引入k-ω SST湍流模型。为了进一步提高计算效率，采用多重网格加速收敛技术、并行计算技术。

转捩的触发与转捩区发展的预测主要在间歇因子输运方程中完成。该方程由Menter等在2004年提出，Langtry等在2006年进行了一些改进，解决了转捩区过短、驻点间歇因子生成项过大等问题。Langtry等在2006年给出的间歇因子输运方程为

(2)

式中:γ为间歇因子;Pγ为生成项;Eγ为耗散项。若“当地转捩”的判据满足，Pγ启动即γ开始增长，Eγ保证在层流边界层中γ趋近于0，这也为预测再层流化现象提供了条件。

(3)

1.3 转捩预测精度验证

图1 NLR7301翼型网格示意图Fig．1 Mesh for NLR7301

图2 升力系数和阻力系数随马赫数 变化曲线Fig．2 Lift，drag coefficients as a function of Mach number

图3 翼型表面转捩位置随马赫数变化曲线Fig．3 Mach number effect on transition onset location

图4 NACA64A010 翼型网格示意图Fig．4 Mesh for NACA64A010

图5 考虑转捩影响的NACA64A010翼型升力系数随攻角变化曲线Fig．5 Lift coefficient versus AOA of NACA64A010 airfoil incorporating transition modeling

图6 考虑转捩影响的NACA64A010翼型力矩系数随攻角变化曲线Fig．6 Moment coefficient versus AOA of NACA64A010 airfoil incorporating transition modeling

2 结构运动方程

具有浮沉和俯仰两自由度的二维翼型结构运动方程为

(4)

(5)

式中:m为机翼质量;Sα为二维翼型对刚心的质量静矩;Dh为浮沉阻尼;Kh为翼型关于刚心的沉浮刚度;Iα为二维翼型对刚心的质量惯性矩;Dα俯仰阻尼;Kα为翼型关于刚心的俯仰刚度;L为升力;My为俯仰力矩。

针对上述二维两自由度翼型结构运动方程，经无量纲化可得

(6)

无量纲颤振速度定义为

(7)

式中:U∞为自由来流速度；μ=m/πρ∞b2为质量比;ρ∞为自由来流密度。

(8)

本文采用基于预估-校正技术的四阶隐式Adams线性多步法[15]对方程(8)进行时域推进求解

(9)

该方法既保证了方程的求解效率，又具有较好的鲁棒性。

3 全湍流动下的跨音速颤振边界验证

跨音速颤振凹坑[16-17]与空气的可压缩性和激波相位滞后效应密切相关，因此跨音速颤振边界预测的关键在于激波捕捉的准确性。本文首先对Isogai案例A模型在全湍流动下的跨音速颤振边界进行预测，以验证本文颤振计算方法的可靠性。其中全湍计算采用k-ω SST湍流模型，结构参数采用文献[11]中的参数，雷诺数按照文献[12]中对所有马赫数给定为6×106。计算网格为上述非定常气动力验证所用网格，如图4所示。

图7、图8分别为采用k-ω SST湍流模型计算得到的全湍条件下的颤振速度边界和颤振频率比边界，与文献[18]的结果吻合很好。数值结果表明，跨音速凹坑最低点的Ma≈0.83，随后颤振速度急剧增大；在Ma=0.87时颤振速度再次减小，在Ma=0.9时出现第二个凹坑，随后颤振速度继续增大。

跨音速凹坑最低点后颤振速度的急剧增大与翼型表面分离流动的出现和扩大密切相关。图9展示了四个不同马赫数下的表面摩阻分布对比，图10展示了颤振边界和准定常升力线斜率随马赫数变化曲线。从图中结果可以看出，当0.83≤Ma≤0.87时，随着马赫数增大，翼型表面的分离区扩大，使得升力线斜率急剧减小，导致颤振边界急剧增大。

图7 颤振速度边界(全湍)Fig.7 Flutter speed index boundary(fully turbulent)

图8 颤振频率比边界(全湍)Fig.8 Frequency ratio boundary(fully turbulent)

图9 不同马赫数下的表面摩阻分布对比Fig.9 Comparison of skin friction distribution between different Mach number

图10 颤振边界和准定常升力线斜率随马赫数变化曲线Fig.10 Flutter boundary and quasi-steady lift curve slope versus Mach number

4 考虑转捩影响的跨音速颤振边界预测

为了探究转捩引起跨音速颤振凹坑加深的原因，对跨音速凹坑最低点附近翼型的定常压力分布和表面摩擦阻力分布进行对比分析，如图13～图15所示。结果表明自由转捩下激波位置较全湍靠后，激波强度大于全湍，且自由转捩较全湍流动在翼型表面提前发生分离。激波增强会降低翼型的稳定性，使得翼型的颤振速度降低。结合以上分析结果可以得出如下结论：自由转捩较全湍使得翼型在跨音速凹坑附近具有更低的稳定性，从而具有更低的颤振边界，导致凹坑程度加深。

图11 颤振速度边界Fig.11 Flutter speed index boundary

图12 颤振频率比边界Fig.12 Frequency ratio boundary

图13 全湍和自由转捩压力及表面摩阻分布对比( Ma = 0． 82)Fig.13 Comparison of pressure and skin friction distribution between fully turbulent and free transition( Ma = 0． 82)

图14 全湍和自由转捩压力及表面摩阻分布对比( Ma = 0． 83)Fig.14 Comparison of pressure and skin friction distribution between fully turbulent and free transition( Ma = 0． 83)

图15 全湍和自由转捩压力及表面摩阻分布对比( Ma = 0． 84)Fig.12 Comparison of pressure and skin friction distribution between fully turbulent and free transition( Ma = 0． 84)

Bendiksen针对跨音速颤振提出了“跨音速稳定性法则[19]”(Transonic Stabilization Laws)，从气动力做功的角度研究了气动力的幅值及相位和颤振稳定性之间的关系，此处的相位指的是机翼俯仰力矩相对于俯仰位移的滞后相位角。本文采用类似的分析方法对Isogai案例A模型在全湍和自由转捩条件下的跨音速稳定性进行分析，其中k为减缩频率。

图16、图17分别为俯仰运动中力矩系数的幅值和相位随马赫数变化曲线。计算结果表明，转捩较全湍使得翼型在跨音速凹坑附近具有更大的力矩系数幅值和更大的超前相位，这意味着此时气动力对翼型的做功增大，使得系统的稳定性降低，从而导致跨音速凹坑程度加深。值得注意的是，当Ma=0.86时，转捩的相位角表明系统处于不稳定状态，气动力对系统做正功，使得系统的稳定性降低；而全湍的相位角表明系统处于稳定状态，气动力对系统做负功，使得系统的稳定性增大。因此当Ma=0.86时，转捩较全湍具有更小的颤振速度，这也使得转捩的跨音速凹坑范围较全湍有所扩大。

图16 俯仰运动中力矩系数的幅值随马赫数变化曲线Fig.16 Amplitude of moment coefficient versus Mach numbers for pitching motion

图18为定常流动0度攻角下自由转捩和全湍的激波位置对比图，当Ma≥0.91时，激波到达翼型后缘，流场进入“冻结区域[20]”(Freeze region)，全湍下的空间马赫数云图如图19～图22所示。由于此时自由转捩与全湍的翼型流场类型，不再予以展示。从图17中可以得知，自由转捩与全湍在“冻结区域”的相位角均接近-180°，这意味着此时气动力基本不做功，系统处于临界稳定状态，翼型的颤振速度随马赫数变化很小。

图17 俯仰运动中力矩系数的相位随马赫数变化曲线Fig.17 Phase angle of moment coefficient versus Mach numbers for pitching motion

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

图19 马赫数云图( Ma = 0． 91)Fig.19 Mach contour( Ma = 0． 91)

图20 马赫数云图( Ma = 0． 93)Fig.20 Mach contour( Ma = 0． 93)

图21 马赫数云图( Ma = 0． 95)Fig.21 Mach contour( Ma = 0． 95)

图22 马赫数云图( Ma = 0． 98)Fig.22 Mach contour( Ma = 0． 98)

5 结 论

(1)本文采用全湍计算得到的颤振边界和文献的结果吻合很好，验证了本文建立的方法的可靠性。

(2)转捩现象在跨音速范围内使得翼型表面激波强度较全湍增强，气动力对翼型做功增加，降低了翼型的颤振稳定性；跨音速凹坑最低点的位置基本保持不变，该点的颤振速度减小了41.6%，且凹坑范围扩大。

(3)当激波到达翼型后缘，进入冻结区域，气动力对翼型几乎不做功，系统处于临界稳定状态，颤振速度随马赫的增大而变化很小。

综上所述，当马赫数小于跨音速凹坑最低点的马赫数时，激波的强度在颤振的影响因素中占据主导地位；而在跨音速凹坑附近，激波和翼型表面的分离现象共同作用，影响翼型的颤振特性，其中激波增强会降低系统的稳定性，分离区域的扩大会增加系统的稳定性；当马赫数接近“冻结区域”，激波靠近翼型后缘，此时激波和分离对颤振特性的影响逐渐减小，翼型的相位角逐渐接近临界状态，系统逐渐趋于稳定，翼型的颤振速度增大；进入冻结区域后，激波到达翼型后缘，而翼型表面几乎不存在分离区，此时激波和分离对颤振特性几乎没有影响，随马赫数的增大，颤振速度的变化很小。

转捩现象对激波强度和翼型表面分离区大小均有影响，使得翼型的气动力非线性较全湍增强，颤振特性与全湍条件下的有明显区别。因此针对机翼表面存在转捩现象的飞机在进行颤振设计时，需要对转捩现象加以考虑，以便得到精确的颤振边界，确保飞机的飞行安全。

[ 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. 考虑转捩影响的翼型动态失速数值模拟[J]. 空气动力学学报， 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] 乔磊. 考虑转捩判定的分离流动数值模拟研究[D]. 西安： 西北工业大学，2013.

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

[15] 叶正寅, 张伟伟, 史爱明, 等. 流固耦合力学基础及其应用[M]. 哈尔滨: 哈尔滨工业大学出版社, 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

国家“973”计划(2014CB744804)

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

李国俊 男，硕士生，1992年生

白俊强 男，博士，教授，1971年生

V215.34

A

10.13465/j.cnki.jvs.2017.22.032