三种典型翼型边界层稳定性对比分析

王 菲，额日其太，王 强，苏沛然

(北京航空航天大学能源与动力工程学院，北京 100191)

0 引言

摩擦阻力可占亚音速运输机总阻力的50%，而层流的阻力比相同雷诺数湍流的阻力要小90%［1］，因此减小表面摩擦阻力可以有效地降低飞行器的总阻力［2］。为了计算机翼的摩擦阻力从而进一步更合理地预测机翼飞行特性，转捩位置的预测成为了重要前提工作。对于高空飞行外界扰动较小的特点，线性稳定性方法的eN方法仍然是目前工程应用中最常用的转捩预测方法。L.J.Runyan和 D.George-Falvy通过修正 Mack的稳定性方法计算了翼型在不同工况下的转捩位置并与实验对照，得到相应的最大增长率的值［3］，为之后的相关研究提供了丰富的数据。朱万琳等人在二维不可压流动的稳定性数值分析基础上进行了转捩的判定［4］。为了便于工程上翼型的选取及优化，针对多个典型翼型边界层稳定性的对比分析具有重要的参考意义。首先本文采用萘升华实验方法验证对称翼型SD8020数值模拟的准确性，并拟定稳定性计算的幅值放大因子的值，然后针对常规翼型(NACA0012)，自然层流翼型(NACA64-204)和超临界翼型(RAE2822)这三种典型翼型采用线性理论方法计算在不同弦长雷诺数情况下转捩位置的变化，研究流向压力梯度对转捩点雷诺数、转捩点弦长百分数和其相应的扰动频率的影响，并分析影响这三种翼型稳定性增长的原因。

1 数值模拟方法

1.1 边界层平均流求解

采用无粘流和边界层方法求解基本流，利用F-S法变换基本流方程［5］，取，这里 ue表示自由流速度。相对于流函数ψ(x，y)的无量纲流函数为:

由此可得连续方程和动量方程为:

边界条件为 η =0，f=f'=0;η =ηe，f'=1。其中导数表示对η求导，m表示无量纲梯度:

1.2 线性稳定性方程

扰动的发展模式可以分成时间模式和空间模式，时间模式的扰动在空间某点仅随时间规律变化，空间模式的扰动传播则沿空间某一方向变化，而后者是可以通过实验结果验证的［6］，所以本文采用更接近实验条件的空间模式方法进行扰动增长计算。将流函数定义为:

这里φ表示流函数的型函数。α和β分别表示流向和展向的扰动波数，ω表示扰动的圆频率。对应于空间模式 α ，β 为复数，α = αr+iαi，β = βr+iβi;ω 为实数。这里定义无量纲扰动相速度c:

利用线性稳定性方程(O-S方程)计算扰动波在流场中的传播，由于在二维基本流条件下三维扰动的失稳点大于二维扰动的失稳点，因此引入的二维扰动也是合理的［7］，即假定β=0。这里二维流无量纲化的O-S方程可以表示为一个四阶常微分方程［1］:

边界条件为 y=0，φ =0，φ'=0;y=δ，φ =0，φ'=0。R 表示弦长雷诺数，这里c表示弦长。

转捩预测采用工程上最为常用的eN方法预测翼型上表面的转捩位置。扰动增长率N为:

其中－αi为给定频率下的空间放大率。

2 实验验证

本实验在北京航空航天大学流体力学研究所低速风洞中完成。实验利用萘的易升华特性，将萘溶于丙酮溶液并喷涂于机翼模型表面。由于转捩点前后气流对壁面的剪切力变化很大，使萘在湍流区内升华速度加快，因此模型表面湍流区的白色萘涂层先消失，从而在层流和湍流之间形成了明显的边界。选择实验状态为Rec=7.8 ×105，Ma=0.087，攻角为 －2°以保持较长的层流区，该翼型转捩位置在78%弦长左右(图1)。

图1 萘升华实验Fig.1 Naphthalene sublimation experiment

对应上述实验条件，采用eN方法对SD8020翼型进行转捩预测。图2为翼型上表面扰动增长曲线，计算20条不同频率的扰动增长曲线，选取常用的扰动放大因子N=9。当某一频率增长率最先达到N=9时认为计算得到转捩点，用虚线与频率增长率曲线的交点对应的弦长位置来表示。计算得到转捩位置在79%左右，与实验结果吻合，因此以下方法均采用e9方法进行转捩预测。

图2 翼型SD8020在弦长雷诺数Rec=7.8×105下不同频率N因子曲线及转捩位置Fig.2 N factor curves and transition locations of SD8020 airfoils for Rec=7.8 ×105

3 计算结果及分析

利用上述eN方法，选用三种典型翼型NACA0012，NACA64-204，RAE2822(图3)，计算来流马赫数 Ma=0.3，攻角为0的条件下这三种翼型压力分布及转捩位置。图4为三种翼型上表面压力系数Cp分布曲线，其中NACA64-204和RAE2822分别对应自然层流翼型和超临界翼型，这两种翼型压力分布特点是都具有较长的顺压梯度区，在前缘附近剧烈的顺压梯度之后，有较长一段弱顺压区，最后达到最低压力点。NACA64-204和RAE2822最低压力点分别在弦长的40%和55%位置处。

图5(a)(b)(c)为三种翼型上表面在相同弦长雷诺数(Rec=9×106)下扰动增长曲线，采用e9方法进行转捩预测。三个翼型中NACA0012的转捩点雷诺数最小，也即在相同弦长雷诺数下最早发生转捩。由于在二维流动中T-S波对转捩的影响起主导作用，而顺压梯度又对T-S波有稳定作用。对应于图4压力系数分布，NACA0012仅有10%弦长的顺压区，远小于后两者，因此其扰动也越不稳定。三种翼型的高频扰动增长起始位置总是比低频更靠近前缘，而受顺压梯度影响，高频增长速度也更容易减慢甚至衰减。相对的低频扰动增长点起始位置更靠后，而其受顺压梯度影响的区域也减少。若低频扰动增长率在进入逆压区时并未开始衰减，则逆压将会导致扰动迅速增长。这种现象与中性曲线的分布规律是一致的。正是由于逆压梯度使速度剖面产生拐点，有拐点的速度剖面的稳定性界限比无拐点的剖面低得多(拐点准则)［7］。

图6(a)、图5(b)和图6(b)表示 NACA64-204翼型分别在弦长雷诺数为6×106、9×106和12×106不同频率扰动增长率。由这三幅图可以看出随弦长雷诺数增长，在层流区内由于顺压区与逆压区的长度比变大，其转捩点雷诺数也变大。当弦长雷诺数增大到12×106时，在顺压区扰动增长率已经达到了N=9，因此可以判定在最小压力点之前转捩已经发生。NACA64-204翼型扰动增长的特点是扰动通过翼型前缘后部大范围弱顺压梯度后，逆压梯度对其影响非常明显，即使很小的逆压梯度也会使低频扰动迅速增长。

弦长雷诺数的增长使高频扰动增长率变大，因此使最不稳定频率(对应转捩位置的扰动频率)也相应提高(图7)，而不同翼型相同弦长雷诺数下对应最不稳定频率也存在很大差别。

图5 三种翼型在弦长雷诺数Rec=9×106下不同频率N因子曲线及转捩位置Fig.5 N factor curves and transition locations of three airfoils for Rec=9×106

如图8所示，同一翼型，随弦长雷诺数增长，中性曲线位置向雷诺数更大频率更高的方向偏移，中性曲线的不稳定区域也变大。相比高弦长雷诺数，低弦长雷诺数下相同频率的扰动更早进入中性曲线的不稳定区域，即扰动增长的更早，且沿频率增长方向扰动不稳定区范围更小。因此随弦长雷诺数的增长，同一翼型最不稳定扰动波频率也变大。

表1统计了三种翼型在弦长雷诺数分别为6×106、9×106、12×106下的转捩点雷诺数和对应弦长百分数。从表中可以看出，工程上常采用雷诺平均方法计算全湍流流场［8］或固定转捩点雷诺数来假设转捩位置的方法是不准确的。同一翼型随着弦长雷诺数的增长，转捩点雷诺数也随之变大，但并未按线性比例增长，其对应的弦长百分数反而减小，也就是说转捩点相对提前。由于转捩位置的变化会对翼型性能计算产生很大的影响［9］，可见准确计算翼型的转捩点并分别求解层流区和湍流区是十分必要的。

图8 NACA64-204不同弦长雷诺数下的中性曲线Fig.8 The neutral curves for different Reynolds numbers for NACA64-204

4 结论

本文首先采用萘升华实验方法验证对称翼型SD8020数值模拟的准确性，并拟定eN方法的幅值放大因子的值，然后采用该数值方法分别计算了在不同的弦长雷诺数下三种典型翼型的压力分布，扰动增长率以及最不稳定扰动波频率，计算了翼型上表面转捩位置，并分析了不同翼型压力分布对扰动发展和转捩位置的影响。结果表明:(1)在相应实验条件下采用幅值放大因子N=9的转捩预测位置与实验结果吻合。(2)在三种典型翼型转捩预测结果中，相同弦长雷诺数下，NACA0012最先发生转捩，而NACA64-204和RAE2822都保持着较长的层流区。(3)压力梯度对扰动增长有很大影响，顺压梯度对二维翼型扰动有明显的稳定作用。(4)同一翼型随弦长雷诺数增长，顺压区内扰动增长更快，因此转捩点相对提前;又因为顺压区占弦长百分比固定，因此层流区内顺压区与逆压区的长度比变大。又因为中性曲线变宽并向下游方向移动，最不稳定频率也随之增大。由于转捩点雷诺数并不是工程上常用的固定值，而转捩位置的变化会对翼型性能参数计算产生很大的影响，因此准确计算翼型的转捩点，可以得到更合理的翼型气动性能参数(如摩擦阻力)，也为工程上选取和优化翼型提供了基本的参考数据。

表1 三种翼型不同弦长雷诺数下转捩点雷诺数及转捩点弦长百分数Table 1 Reynolds numbers and the percentage of chord length at transition points for three a irfoils

［1］周恒，赵耕夫.流动稳定性［M］.北京:国防工业出版社，2004:226.

［2］额日其太，沈遐龄，刘火星，等.飞机层流控制技术及其在大型运输机上的应用［A］.大型飞机关键技术高层论坛暨中国航空学会2007年学术年会［C］.中国广东深圳:2007.

［3］RUNYAN L J，GEORGE-FALVY D.Amplification factors at transition on an unswept wing in free flight and on a swept wing in wind tunnel［A］.17th Aerospace sciences meeting［C］.New Orleans，Louisiana:AIAA，1979.

［4］朱万琳.二维层流转捩计算及应用［D］.西北工业大学，2002.

［5］陈懋章编著.粘性流体动力学基础［M］.北京:高等教育出版社，2002.

［6］CRIMINALE W O，JCAKSON T L，JOSLIN R D.Theory and computation of hydrodynamic stability［M］.Cambridge University Press，2003.

［7］SCHLICHTING H.Boundary-layer theory［M］.7th ed.MCGraw-Hill，1979.

［8］DRIVER J，ZINGG D W.Numerical aerodynamic optimization incorporating laminar-turbulent transition prediction［J］.AIAA Journal.2007，45(8):1810 －1818.

［9］KRUMBEIN A.Automatic transition prediction and application to 3d wing configurations［A］.44th AIAA Aerospace sciences meeting and exhibit［C］.Reno，Nevada:2006.