小攻角对后掠机翼边界层稳定性及转捩的影响

靖振荣，孙朋朋，黄章峰*，2

(1.天津大学 力学系，天津 300072;2.中国空气动力研究与发展中心 空气动力学国家重点实验室，绵阳 621000)

机翼作为航空飞行器上产生升力的重要部件之一，其设计的水平很大程度上决定了飞机性能的优劣.在飞机飞行阻力中，摩擦阻力是一个重要的组成部分.对于亚声速飞机，表面摩擦阻力占总阻力的50%左右［1］.使飞机在飞行时机翼表面尽可能长地保持层流状态是降低摩擦阻力的最有效的办法之一.因此，深入研究转捩机理、预测转捩位置、防止或推迟转捩对机翼的优化设计和气动减阻有着十分重要的意义.

飞机一般采用后掠机翼，其目的是提高飞行的临界马赫数、克服因接近音速飞行而急剧增大的空气阻力、提高飞行稳定性.在后掠机翼头部边界层内，由于压力梯度和后掠角的共同作用，使得与势流方向相垂直的方向有速度分量，称为横流［2］.横流失稳是引起后掠机翼边界层流动转捩的主要因素［3］.攻角是指机翼翼弦与自由来流方向的夹角，是机翼产生升力的主要因素之一，一般不为零.因此，有攻角的后掠翼成为人们的研究热点.

在国外，Boltz等［4］在低湍流风洞中对后掠机翼边界层进行了研究，发现随着攻角的减小，同一位置出现转捩所需的雷诺数逐渐增大，转捩提前发生在机翼的迎风面.Haynes［5］采用非线性抛物化方程(NPSE)研究了横流涡的非线性特性和扰动幅值饱和特性，发现在小攻角下，横流扰动波的流向增长范围随着攻角的增加几乎呈线性增长.Dagenhart和 Saric［6］实 验 发 现 － 4°攻 角 的NLF(2)-0415后掠机翼会产生较合适的压力梯度，对Tollmien-Schlichting(T-S)波有稳定作用，对横流涡的增长却有促进作用，采用eN方法预测的转捩位置的 N 值范围为6.4～6.8.Reibert和Saric［7］对后掠机翼边界层横流失稳及转捩的实验研究进行了综述.Bushnell等［8］总结了线性理论预测转捩位置的成果，认为在粗糙度和背景湍流度较小时，eN方法可以应用到马赫数、攻角、压力梯度、后掠变化等情况，而且当来流湍流度小于0.05%时转捩对应的 N 值为 9 ～11.Arnal等［9］利用eN方法研究了有无抽吸时的后掠翼的转捩，发现当横流转捩占主导时，N值为6～10，并且抽吸可以减小N值.

在国内，孙朋朋和黄章峰［10］研究了后掠角对横流不稳定性的影响，发现在机翼头部靠近前缘位置，后掠角为45°时横流强度最强，驻波(ω=0°)在后掠角为30°～45°时不稳定区域的范围最大.李峰等［11］采用数值模拟方法研究了低速大攻角下机翼的绕流特性，发现在某些马赫数和攻角下NACA0012翼型的湍流解具有周期性.吴鋆、等［12-13］采用实验方法研究了攻角对翼型绕流结构的影响，发现分离点和分离剪切层形成旋涡的位置随迎角的增大而向上游移动，同时翼型上表面流动分离后形成的回流区尺寸随着翼型迎角的增加而增大.袁湘江等［14］采用直接数值模拟(DNS)和线性稳定性(LST)分析方法，对高超声速小攻角钝锥边界层的失稳机制和转捩特点进行了研究，发现攻角为2°的高超声速钝锥边界层存在多支不稳定模态.对于机翼边界层稳定性问题，王斌等［15］采用风洞试验研究了攻角对后掠机翼转捩位置的影响，发现机翼上表面转捩位置随着攻角的增大而前移，当达到某一攻角时，转捩快速接近前缘.

因此攻角对横流的产生和演化都有很大的影响，而且在后掠机翼头部的迎风面与背风面边界层内的压力梯度有所不同，攻角对迎风面和背风面的横流稳定性影响也将不同.本文以适航的后掠机翼为研究对象，在给定后掠角及雷诺数的情况下，结合数值模拟和线性稳定性理论，研究不同攻角对后掠机翼边界层在机翼迎风面和背风面的流动稳定性的影响，并采用转捩预测的eN方法进行转捩预测，为机翼的优化减阻设计提供借鉴和参考.

1 数值计算方法

采用有限体积法，求解三维可压缩无量纲Navier-Stokes方程:

式中:t为时间;x和y为流向和法向坐标;U为守恒型通量，由密度、3个方向的速度(流向、法向和展向)和温度组成;E和F为非线性项通量，其中包含压力项;Ev和 Fv为黏性项通量［16］.由于机翼在展向无限长，因此机翼和流动在不同的展向位置均相同，其在展向的导数恒为零.

无量纲黏性系数μ由Sutherland关系式确定:

式中:T 为温度;C=110.4/T∞，T∞为来流温度.

数值计算中，时间导数采用LU-SGS离散;对流项采用Steger-Warming格式进行通量分裂，采用3阶WENO格式进行重构.网格外边界采用远场边界条件;壁面采用无滑移绝热边界条件;展向采用周期边界条件.

得到基本流场后，将Navier-Stokes方程写成扰动形式，并进行线化就得到了经典的Orr-Sommerfeld方程.采用2点4阶的Malik差分法，将Orr-Sommerfeld方程离散后得到特征关系式:

式中:α和β分别为扰动在流向和展向波数;ω为扰动频率;f为与基本流有关的量.对于空间模式;β、ω为实数，α为复数，其虚部 －αi为扰动的增长率.给定β、ω及f，采用Mueller法，就可以求出特征值α和特征函数φ.通过变动β、ω及f使αi=0，即扰动的增长率为零，从而得到中性曲线.对结果进行稳定性分析，进而采用转捩预测eN方法进行预测［17］.

2 计算模型和条件

采用NACA0012模型，弦长1 m，后掠角取25°.计算模型坐标示意图如图1所示.其中x、y、z为直角坐标，η、ξ为贴体坐标;Λ 为后掠角;为贴体坐标下的速度分量为势流方向速度;为横流速度.选取高空10 km处气体参数作为参考量，温度为223.3K，来流马赫数为 0.7，单位雷诺数为 5.945 ×106/m.攻角 α 选取 0°、1°、2°、3°、4°、5°.若无注明，文中的量均为无量纲量，其中速度用来流速度U∞无量纲化，长度采用1mm无量纲化，相应的无量纲雷诺数为5945.

图1 计算模型坐标示意图Fig.1 Sketch map of the model coordinate

计算域及网格如图2所示.其中，x方向为机翼弦长方向，y方向为机翼法向方向.机翼流向网格最小间距小于0.1mm，法向上一个边界层内布有200个点以上，法向网格最小间距小于0.01mm.

对韩步璋等［18］实验工况的基本流进行数值模拟，得到的压力系数与实验结果的对比如图3所示.可以看出，计算结果与实验结果符合得很好，说明本研究计算的结果可靠.

图2 计算网格示意图Fig.2 Sketch map of the computational domain

图3 压力系数沿弦长方向分布Fig.3 Distribution of pressure coefficient along chord

图4给出了在横流坐标系下，21.6万、40万、60万3种网格模型在不同位置处的横流速度剖面，其中，为边界层外沿处的势流速度，方向与边界层外沿势流方向垂直.可以看出，3种网格模型计算得到的结果基本一致，说明模型网格已经够密.

图4 横流速度剖面Fig.4 Profile of the crossflow velo city

3 数值结果分析

3.1 基本流分析

图5给出了不同攻角情况下机翼的压力系数曲线.从曲线中可以看出，随着攻角的增大，机翼两侧面的压力系数之差逐渐增大，说明产生的升力逐渐增加.当攻角α为4°和 5°时，在 x/c=0.1～0.3位置，机翼背风面压力系数出现突降现象.这是由于当攻角α≥4°时，在机翼的背风面一定范围内会有激波出现［18］.

图5 不同攻角下的压力系数Fig.5 Pressure coefficients with different attack angles

选取数值模拟计算的流场截面，建立新的坐标系:ξ的大小定义为在翼型剖面上距离翼型前缘的壁面弧长距离，方向沿翼型剖面的弧线方向;η的大小定义为该点到壁面的距离，方向与壁面法向方向相同.重新划分网格并加密处理，得到可用于稳定性分析的基本流场.

图6(a)和图6(b)分别给出了不同攻角工况下，在ξ=100位置处迎风面和背风面的流向速度剖面.由图可知，随着攻角的增大，边界层附近的流向速度在迎风面逐渐减小，在背风面逐渐增大.在η＞1处流向速度都接近于常值.在迎风面边界层中流向速度小于背风面边界层中的流向速度.

图6 在ξ=100处的流向速度 剖面Fig.6 Profile of flow velocity atξ=100

图7 在ξ=100处横流速度剖面Fig.7 Profile of crossflow velocity at ξ=100

图7(a)和图7(b)对比了迎风面和背风面在横流坐标系下，不同攻角工况时ξ=100位置处的横流速度剖面.可以看出，横流速度在其拐点处达到峰值，且随着攻角的增大，峰值的大小在迎风面逐渐增大，在背风面逐渐减小.攻角不为零时，迎风面的横流速度大于背风面的横流速度.在迎风面，横流速度的峰值均在η=0.2mm位置左右，而在背风面，横流速度的峰值所在位置随着攻角的增加而增加.

3.2 中性曲线

图8(a)和图8(b)给出了背风面和迎风面在不同攻角情况下，展向波数β=3的扰动波在(ξ，ω)的中性曲线(ω为扰动频率).可以看出，在迎风面，随着攻角的增大，展向波数β=3的扰动波的中性曲线的流向范围逐渐增大.频率较大的扰动波对攻角的变化不太敏感，而频率较小的扰动波受到攻角的影响很大，而且对于攻角小于4°的小攻角，中性曲线在流向的最大位置随着攻角的增加而线性增加，与 Haynes［5］的结果一致.在背风面，中性曲线的流向范围随攻角的增大而减小，而且中性曲线在流向的最大位置随着攻角的增加而线性减小.随着攻角的增加，在背风面出现了高频的不稳定波.迎风面和背风面在(ξ，ω)平面上的中性曲线与ξ坐标轴都有交点，而且不同的攻角时交点的位置不同，这说明不稳定区域存在频率ω=0的驻波，而且攻角对驻波的稳定性影响很大.

图8 中性稳定曲线(β=3)Fig.8 Neutral stability curves(β =3)

图9(a)和图9(b)给出了迎风面和背风面在不同攻角下驻波(ω=0)在(ξ，β)平面上的中性曲线.可以看出，展向波数β较大的驻波的流向增长区域随攻角的增大变化较小，而展向波数较小的波(1＜β＜4)的流向增长区域随攻角的增大变化较大.随着攻角的增大，迎风面驻波的流向增长区域逐渐增大，背风面驻波的流向增长区域逐渐减小，且在迎风面的增长范围大于背风面.说明相对于背风面，迎风面更不稳定.类似(ξ，ω)平面上的中性曲线，小攻角时(ξ，β)平面上迎风面/背风面的中性曲线在流向的最大位置随着攻角的增加而线性增加/减小.

图9 中性稳定曲线(ω=0)Fig.9 Neutral stability curves(ω =0)

3.3 N值曲线及转捩预测

图10给出了迎风面和背风面在不同攻角工况下采用转捩预测eN方法计算的N值曲线.可以看出，迎风面的N值相比背风面增长较快，且在迎风面的N值要远大于背风面的N值.对比可知，当机翼攻角不为零时，扰动波在迎风面的N值最大，说明扰动波在迎风面传播时被放大的倍数最大.如图10所示，当机翼攻角为3°时，扰动波的N值可达到 10，即扰动波可被放大 e10倍(约2.2万倍)，表明迎风面更容易发生转捩.

攻角增大时，在迎风面，来流速度方向与迎风面壁面切线方向的夹角增大，气流在壁面附近的偏转角也增大，因此气流沿着流向的速度减小，更多的流体流向展向方向，导致迎风面的横流速度增强，而在背风面反之.从图6可以看出，随着攻角的增加，迎风面的流向速度逐渐减小，而背风面的流向速度逐渐增大，从图7可以看出随着攻角的增加，迎风面的横流速度明显增大，而背风面的横流速度明显减小.横流不稳定性是由于拐点产生的，其不稳定性的强度与横流速度的最大值及其位置息息相关.横流速度的最大值越大，横流不稳定性越强.图8和图9中迎风面的中性曲线范围所包含的面积要大于背风面，而且图10中迎风面的N值要明显大于背风面的N值，说明相对于背风面，迎风面更不稳定.

图10 eN方法的N值曲线Fig.10 N factor in eN method

图11为在不同位置处，扰动波N值随攻角的变化曲线，其中攻角为正值表示为迎风面，攻角为负值表示为背风面.可以看出，在同一位置处迎风面的N值随着攻角的增大而增大，背风面的N值随攻角的增大而减小.这说明在攻角不为零的情况下，扰动波的增长在迎风面得到了加强，在背风面受到抑制.攻角的大小与扰动波受到增强和抑制的程度有关，攻角越大，扰动波受到抑制和得到加强的程度越大.

图11 N值随攻角的变化曲线Fig.11 Curves of N factor changing with attack angle

研究表明，当扰动波的无量纲幅值达到来流的20%就会发生转捩［19］.因此当飞行器在空中飞行时，机翼边界层流动中扰动的初始幅值A0的大小对流动能否转捩有着极其重要的影响.图12给出了不同攻角下ξ位置处发生转捩时的初始位置处(ξ=25)的扰动幅值A0.可以看出，随着攻角的增大，迎风面引起转捩所需的初始扰动幅值逐渐减小，说明转捩越容易发生;而在背风面，随攻角的增大引起转捩所需的初始扰动幅值逐渐增大，说明转捩越难以发生.

图12 初始扰动幅值A0随攻角的变化曲线Fig.12 Curves of amplitude A0 changing with attack angle

在工程应用中，判断转捩发生的预设值需要由实验测得并需要进一步假定.若采用eN方法计算的扰动的N值达到预设值时，则认为转捩发生了.根据图13给出了当计算的扰动波的N值达到预设值发生转捩的位置.可以看出，随着攻角的增大迎风面转捩位置向机翼前缘移动，背风面则与之相反.

图13 转捩位置ξ随攻角的变化曲线Fig.13 Curves of transition position ξ changing with attack angle

本文的来流参数选取为高空10 km处，来流马赫数为0.7，预计来流中扰动速度约为0.1m/s(约为来流速度的0.05%)，当扰动波的无量纲幅值达到来流的20%就会发生转捩［19］，因此预计转捩发生的N值为6左右.从图13可以看出，此时转捩发生的位置在0.1～0.2个弦长之间.实际发生的转捩位置与来流扰动的幅值息息相关.

4 结论

在采用数值模拟方法得到无限展长后掠机翼基本流场的基础上，通过求解O-S方程，研究了攻角对后掠机翼边界层稳定性的影响并进行了转捩预测，得到结论如下:

1)随着攻角的增大，流向速度在迎风面逐渐减小，而横流速度线性增大，在背风面反之.

2)攻角对中性驻波的展向波数影响较小，但是驻波的中性曲线在流向的最大位置对攻角很敏感，在迎风面/背风面随着攻角的增加呈线性增加/减小.

3)扰动波的增长在迎风面得到加强，在背风面受到抑制.

4)迎风面转捩最容易发生，随着攻角的增大，转捩位置向机翼前缘移动.

5)当扰动速度为来流速度的0.05%时，转捩发生的 N值在6左右，转捩发生的位置在0.1 ～0.2 个弦长之间.

References)

［1］ Joslin R D.Overview of laminar flow control［M］.Virginia:National Aeronautics and Space Administration，Langley Research Center，1998.

［2］周恒，赵耕夫.流动稳定性［M］.北京:国防工业出版社，2004:157-158.Zhou H，Zhao G F.Hydrodynamic stability［M］.Beijing:National Defense Industry Press，2004:157-158(in Chinese).

［3］吴永健.横流不稳定性实验研究［D］.南京:南京航空航天大学，2002.Wu Y J.Experimental study on crossflow instabilities in the boundary-layer of swept wing［D］.Nanjing:Nanjing University of Aeronautics and Astronautics，2002(in Chinese).

［4］ Boltz FW，Kenyon G C，Allen C Q.Effects of sweep angle on the boundary-layer stability characteristics of an untapered wing at low speeds，Technical Note:D-338［R］.Moffett Field:National Aeronautics and Space Administration，1960.

［5］ Haynes T S.Nonlinear stability and saturation of crossflow vortices in swept-wing boundary layers［D］.Tempe:Arizona State University，1996.

［6］ Dagenhart JR，Saric W S.Crossflow stability and transition experiments in swept-wing flow［M］.Virginia:National Aeronautics and Space Administration，Langley Research Center，1999:1，7.

［7］ Reibert M S，Saric W S.Review of swept-wing transition，AIAA-1997-1816［R］.Reston:AIAA，1997.

［8］ Bushnell D M，Malik M R，HarveyW D.Transition prediction in external flows via linear stability theory［C］∥Symposium Transsonicum Ⅲ.Berlin Heidelberg:Springer-Verlag，1988:225.

［9］ Arnal D，Gasparian G，Salinas H.Recent advances in theoretical methods for laminar-turbulent transition prediction［J］.AIAA，1998:98-0223.

［10］孙朋朋，黄章峰.后掠角对后掠机翼边界层稳定性及转捩的影响［J］.北京航空航天大学学报，2015，41(7):1313-1321.Sun P P，Huang Z F.Effect of the sweep angle on the stability and transition in a swept-wing boundary layer［J］.Journal of Beijing University of Aeronautics and Astronautics，2015，41(7):1313-1321(in Chinese).

［11］李锋，汪翼云，崔尔杰.翼型大攻角绕流的数值模拟［J］.航空学报，1992，13(1):17-22.Li F，Wang Y Y，Cui E J.The numerical simulation of compressible flow around an airfoil at high angle of attack［J］.Acta Aeronautica et Astronautica Sinica，1992，13(1):17-22(in Chinese).

［12］吴鋆，王晋军，李天.NACA0012翼型低雷诺数绕流的实验研究［J］.实验流体力学，2013，27(6):32-38.Wu J，Wang J J，Li T.Experimental investigation on low Reynolds number behavior of NACA0012 airfoil［J］.Journal of Experiments in Fluid Mechanics，2013，27(6):32-38(in Chinese).

［13］吴鋆，李天，王晋军.低Reynolds数NACA0012翼型绕流的流动特性分析［J］.实验力学，2014，29(3):265-272.Wu J，Li T，Wang J J.Characteristic analysis of flow around NACA0012 airfoil in a low-Reynold-number media［J］.Journal of Experimental Mechanics，2014，29(3):265-272(in Chinese).

［14］袁湘江，李国良，刘智勇，等.小攻角高超声速钝锥边界层失稳特性［J］.航空动力学报，2011，26(12):2805-2811.Yuan X J，LiG L，Liu ZY，etal.Study of the instability characteristic in the boundary layer of a hypersonic blunt cone at low angle of attack［J］.Journal of Aerospace Power，2011，26(12):2805-2811(in Chinese).

［15］王斌，白存儒，杨广郡，等.后掠机翼低速流动转捩位置的升华法测量［J］.实验力学，2009，24(3):197-201.Wang B，Bai C R，Yang G J，et al.Measurement of transition location change of swept wing in a low speed flow based on sublimation method［J］.Journal of Experimental Mechanics，2009，24(3):197-201(in Chinese).

［16］孙朋朋.马赫数、攻角及后掠角对后掠机翼边界层稳定性及转捩的影响［D］.天津:天津大学，2015.Sun P P.Effect of Mach number，attack angle and sweep angle on the stability and transition in a swept-wing boundary layer［D］.Tianjin:Tianjin University，2015(in Chinese).

［17］黄章峰，逯学志，于高通.机翼边界层的横流稳定性分析和转捩预测［J］.空气动力学学报，2014，32(1):14-20.Huang Z F，Lu X Z，Yu G T.Cross-flow instability analysis and transition prediction of airfoil boundary layer［J］.ACTA Aerodynamic Sinica，2014，32(1):14-20(in Chinese).

［18］韩步璋，黄奕裔，张其威，等.NACA0012翼型跨音速测压实验研究［J］.南京航空航天大学学报，1987，19(2):92-102.Han B Z，Hang Y Y，Zhang QW，et al.An experiment of pressure measurement for NACA0012 airfoil in a transonic wind tunnel［J］.Journal of Nanjing Aeronautical Institute，1987，19(2):92-102(in Chinese).

［19］ Huang Z F，Cao W，Zhou H.The mechanism of breakdown in laminar-turbulent transition of a supersonic boundary layer on a flat plate-temporal mode［J］.Science in China Series G:Mechanics and Astronomy，2005，48(5):614-625.