自然层流机翼的翼套试验及数值方法

陈艺夫，王一雯，邓一菊，王波，白俊强，卢磊

1. 西北工业大学 航空学院，西安 710072

2. 西北工业大学 无人系统技术研究院，西安 710072

3. 航空工业第一飞机设计研究院，西安 710089

4. 中国科学院 工程热物理研究所，北京 100190

近年来，全球气候变暖愈发严重，联合国气候变化大会要求全球积极采取措施应对气候变化。进而导致各大航空公司对飞机的经济性和环保性要求逐步提升。在众多的减阻技术之中，层流减阻技术具有巨大的应用前景[1-2]。研究表明，平板边界层内层流区的摩擦阻力比湍流区的要低90%，而摩擦阻力占商用运输机总阻力约50%[3]，因此，延长飞机表面层流区将带来很大的减阻收益[4]。

层流技术的核心是推迟转捩、延长层流区，从而实现减阻增升的目标。对于小后掠角在中等雷诺数运营环境中的构型，主要是通过自然层流(Natural Laminar Flow，NLF)技术，借助有利压力分布抑制边界层内不稳定波的增长，达到推迟转捩的目的。围绕自然层流技术，欧美进行了大量研究。早期美国NACA(National Aeronautics and Space Administration)设计了包含NACA1、NACA2-5、NACA6、NACA7系列的层流翼型[5]。

将数值方法与风洞试验相结合是NLF机翼设计的一种广泛使用的研究方法。目前基于线性稳定性理论的eN方法已经在工程设计中得到了广泛应用。基于eN方法，Campbell和Lynde[6]设计了雷诺数较高的跨声速NLF机翼风洞模型。Cella等[7]成功使用eN方法设计了UW-5006跨声速层流机翼，并通过风洞测试验证了设计结果。Hue等[8]进行了风洞试验，并结合eN方法进行了CFD研究，重点是层流机翼的转捩预测。Shi等[9-10]针对三维构型耦合雷诺平均纳维埃斯托克斯方程(Reynolds Average Navier Stokes, RANS)求解器建立了基于线性稳定性理论的有效转捩预测工具，转捩预测结果和风洞试验数据匹配良好。

国内的风洞试验研究方面，耿子海等[11]对二维NACA0006翼型进行了混合层流控制风洞试验研究，王菲[12-13]、邓双国[14]等分别基于NACA64A-204翼型后掠翼模型和双圆弧对称翼型后掠翼模型进行了层流控制试验研究，对层流流动的机理有了进一步理解。张彦军[15]和Xu[16]等完成了超临界自然层流机翼的设计并用风洞试验研究了层流机翼在不同马赫数和雷诺数下的转捩特性，试验最大雷诺数可达1 000万。李强[17]、陈苏宇[18]等针对高超声速边界层转捩问题进行了风洞测量试验。

虽然风洞试验的数据精度较高，但是受风洞试验段大小的限制，其雷诺数通常低于1 000万，这远小于运输机的典型雷诺数。尽管利用加压和低温技术可以使风洞雷诺数增加到约2 000万，但与此同时湍流水平也变得比常规大气飞行环境中的湍流水平更高。例如，NASA基于通用研究模型(Common Research Model, CRM)设计了超临界NLF机翼[19-20]，并在雷诺数约为2 000万的国家跨声速风洞(National Transonic Facility, NTF)中进行了测试验证。但是由于湍流度相对较高，T-S波临界N因子大约为6，这比常规飞行条件下使用的N因子小得多。因此，这些数据对于验证具有高雷诺数的层流机翼设计的CFD结果不是直接有用的，这促使研究人员进行飞行测试以验证和改进CFD工具。

早在20世纪60年代，美国实施了X-21层流飞行试验，验证了飞行条件下实现层流的可能性[21]。而后德国也相继开展了类似的飞行试验，在LFC-205和VFW-614飞机上实现了自然层流[22]。20世纪80年代以来，随着制造条件及试验技术的进步，NASA进行了大量的层流机翼的飞行试验研究，如F-111/TACT跨声速飞行试验[23]、F-14A飞机变后掠飞行试验[24]以及波音757飞行平台的层流翼套试验[25]等，研究了后掠角、马赫数以及发动机噪声水平对机翼表面层流的影响。此外，Belisle[26]和Roberts[27]等通过在湾流-Ⅲ飞机的机翼上设计翼套装置进行离散粗糙元的被动层流控制，并进行了数值模拟分析。基于几十年的风洞试验及飞行试验的技术积累，波音公司逐步将层流技术应用到波音737的翼梢小翼、波音787-9的尾翼及短舱等部件的气动设计中[28]。欧洲也对层流机翼进行了大量的研究。欧盟进行了“清洁天空”计划[29]，在此计划下，空客公司在A340-300的外翼上加装了自然层流翼套进行了层流验证飞行，并进一步计划将混合层流尾翼安装在远程客机上进行试飞试验。总之，国外对自然层流技术的研究已非常深入，并进行了广泛的飞行试验，部分成果已进行工程转化。

结合以上的总结和讨论，自然层流风洞试验和飞行试验都能为转捩预测数值模拟方法提供验证性数据。两种试验方法在来流湍流度和雷诺数上存在一定差异，这对转捩预测方法的可靠性和鲁棒性提出了较高的要求，同时需要探明风洞试验对飞行试验的前期指导意义。对于eN方法，如果在风洞试验条件标定的N值下能够达到较高的转捩预测精度，那么需要探究在飞行条件下是否仍能获得较好的预测效果。先前的飞行试验并未涵盖所有飞行场景。实际上，大多数以前的自然层流控制飞行试验都是在巡航条件下进行的，这不能完全反映飞行条件下自然层流控制的转捩现象。因此，需要开展具有高雷诺数和不同攻角的飞行测试试验，可为相同条件下的自然层流设计提供参考和指导，也能更有效地验证转捩预测方法的鲁棒性。最后，在公开的文献和国内的研究中，鲜有同时利用风洞试验数据和飞行试验数据对数值转捩预测方法的验证研究。

因此，基于某型公务机飞行试验平台，对设计的自然层流翼套进行层流特性研究。以层流机翼翼套飞行试验为最终的目标，需要首先采用风洞试验手段对翼套的层流特性进行研究；其次，基于试验结果验证本文采用的eN转捩预测方法的预测精度，并结合风洞试验和数值模拟两种手段分析气动攻角、压力分布形态对T-S波失稳主导的转捩的影响机制；最后，对飞行条件下层流翼套转捩特性进行研究并与飞行试验数据进行对比，探明转捩预测方法具有较好的可靠性和鲁棒性。

1 eN转捩预测方法及其求解

1.1 eN转捩预测方法

eN法是当前工程界中较为广泛应用的一种预测转捩的半经验方法，又叫N因子法。eN法基于稳定性理论，分析小扰动在边界层内的演化，进而判断出转捩的位置。

采用eN方法预测转捩。首先利用边界层方程求解得到边界层信息(速度型、温度、密度、压力等)，然后将这些信息作为稳定性分析的输入，进行边界层内小扰动的稳定性求解，当扰动随着弦向慢慢增大直至达到临界阈值失稳即发生转捩。基于线性稳定性理论(Linear Stability Theory, LST)进行稳定性分析，其基本思想是将流场变量假设为一个定常的时均流动和一个非定常扰动，即

(1)

(2)

(3)

式中：A、B、C由时均流场参数及α、β决定。其求解采用当地化反Rayleigh商方法，详细参考Cebeci[30]的论文。通过稳定性方程可以得到扰动的时间增长率ωi，进而得到空间角度的扰动增长率αi[31-32]，使用式(4)在流向积分即可得到放大N因子，ds为弦向微元弧长：

(4)

N因子的积分求解有不同的方法，采用NTS/NCF方法进行T-S波N因子NTS、CF涡N因子NCF放大因子积分。之后，辅以转捩阈值即可得到转捩位置。由于本文研究构型后掠角较小，边界层内T-S波扰动占主导作用，可忽略CF涡扰动。因而，主要关注流向扰动放大因子及其临界阈值。转捩阈值一般由风洞试验或飞行试验标定得到，与飞行条件(如湍流度)有关。Mack[33]给出了T-S波临界N因子NTScr与自由来流湍流度Tu的关系：

(5)

式中：自由来流湍流度定义为湍流脉动速度均方根U′和自由来流时均速度U∞之比，其表达式为

(6)

1.2 RANS方程与eN方法的耦合求解

采用RANS方程耦合eN方法的迭代求解进行转捩预测，图1给出了求解的迭代流程。首先，RANS求解器得到的压力分布数据以及计算状态输入到边界层方程，边界层方程求解得到边界层信息，作为线性稳定性分析的输入，稳定性分析可以模拟边界层内扰动发展情况，进而通过建立的eN方法得到转捩位置，以间歇因子方程的形式对层流进行模化，并在RANS求解器中实现固定转捩计算。eN方法和RANS求解器进行松耦合迭代直至转捩位置收敛即可实现对转捩的数值求解。

图1 RANS方程与eN方法的耦合求解过程

2 自然层流翼套的风洞试验方案及测量技术

2.1 风洞试验模型

选取翼套翼型剖面，拉伸成等直机翼并使后掠角与翼套前缘一致(均为5°)，形成图2所示的翼身组合体试验构型。试验机翼弦长0.3 m，展长0.9 m。在机翼段距对称面450、600、750 mm处布置了3排测压孔，每排测压孔设置25个测压点。为避免上游测压孔产生的尾流对下游测压孔产生影响，测压孔与流动方向呈15°排列。图2 中黄色区域是可更换的机翼前缘壁板，通过更换不同类型的壁板可进行自然层流或混合层流机翼风洞试验。试验仅对自然层流机翼进行研究，试验观测区域在机翼展向内、外侧的测压剖面之间。

图2 风洞试验模型平面

风洞试验在中国航空工业空气动力研究院沈阳院区的FL-3风洞进行(图3)。FL-3风洞是直流暂冲下吹式亚、跨、超三声速风洞。试验段截面尺寸为1.5 m×1.6 m，试验马赫数范围是0.3～2.25。试验中，使用开槽壁试验段，试验段上、下壁板开有非线性槽缝，左、右侧壁为实壁。

图3 FL-3风洞

试验采用半模构型，需要在机身与洞壁间添加附面层隔板，支撑方式为半模转窗。试验构型在风洞中的安装情况见图4。

图4 风洞试验构型及安装

2.2 红外热像转捩探测及其验证

利用红外热像技术进行层流机翼边界层转捩探测。由于湍流边界层的对流热系数远大于层流边界层，当模型表面与气流之间存在温度差时，湍流区的模型表面温度将更快地接近气流温度，进而使得模型表面的层流区与湍流区之间产生温度差异。红外热像技术利用这一特性，通过红外热图获取模型表面温度分布，进而根据温度梯度判定转捩位置。

具体测量方法为：在风洞模型表面观测区域喷涂亚光黑漆，提高表面红外发射率，并在试验吹风前利用卤素灯对模型进行烘烤加热，提高模型表面初始温度。由于吹风时气流温度较低，与模型间存在温差，在吹风过程中使用红外相机获取待测机翼表面红外热图，最终根据表面温度分布来判断转捩位置。

为了对模型转捩位置的测量精度进行验证，进行了固定转捩的测试试验。测试方法为：在机翼弦向10%位置布置粗糙带(Roughness Ribs)，在粗糙带局部留有一定间隙用于显示自然层流的状态，方便直观对比现象。在试验条件下对不同攻角下的模型进行吹风试验并获取红外热图。固定转捩使用的粗糙带参数如下：高度0.1 mm，粗糙元直径1.2 mm，相邻单元圆心距2.5 mm。根据边界层理论以及Braslow[34]的研究结论，诱发边界层在当地转捩所需的粗糙元高度对应的临界雷诺数Rek应达到600～650。根据试验条件参数估算，Rek的范围是700～800，满足诱发当地转捩的条件，因此本文方法进行转捩探测结果的验证是可靠的。

图5和表1分别给出了4°攻角下固定转捩试验的实测图和不同攻角下的固定转捩测试结果，根据实测数据可知机翼表面在粗糙带当地发生了转捩，而在未布置粗糙带的区域层流正常发展,达到20%c左右，c为弦长。从表1的结果可以看出，不同攻角下的固定转捩测试结果与粗糙带的位置偏差都在1%以内，证明试验的转捩测试方法偏差较小，转捩探测方法是可信的。

图5 4°攻角固定转捩试验实测效果

表1 不同攻角下的固定转捩试验转捩位置测试结果

3 自然层流翼套风洞试验及数值模拟

3.1 风洞试验结果

风洞试验在马赫数Ma=0.5、0.6、0.7，雷诺数Re=3.37×106的条件下进行，对不同攻角下机翼的转捩位置进行了测量，测量结果见图6。图中给出了攻角从-2°到4°变化时，机翼表面的层流和湍流的分布情况，其中颜色较亮的区域是层流，较暗的区域是湍流，转捩的观测效果非常明显。试验中监测的层流区呈楔形及尖峰形分布，部分区域提前发生了转捩，其原因可能是模型加工过程中表面部分区域不够光滑或在试验过程中翼面上粘附了微小的污染物，但该现象并不影响转捩位置的判断，转捩位置在层流段最长的区域进行选取。图中试验观测区域内两个菱形的点是设置的位置标志点，其中左、右标记分别在55%c和50%c的位置。

图6 机翼表面转捩实测结果

从试验结果可以看出，-2°、0°、2°这3个攻角下层流区范围随攻角的增大基本不变，转捩位置都在50%c左右。当攻角增大到4°时，转捩提前至20%c，层流区大大缩短。

3.2 数值模拟与风洞试验结果的对比分析

使用1.2节中介绍的RANS耦合eN的转捩预测方法对风洞试验构型进行转捩计算。图7给出了试验构型的RANS计算网格。由于是等直段机翼，压力分布沿展向变化不大，所以在翼展方向仅布置73个网格点，而翼面和前缘曲率及压力分布变化明显，因此在机翼上下翼面各布置101个点，在前缘布置33个点，以满足CFD计算精度要求。附面层网格采用O型拓扑，物面第一层网格高度1×10-6m，网格法向增长率1.15，最终生成网格量700万。CFD计算的条件与风洞试验状态保持一致。

图7 风洞试验构型网格分布

风洞试验与数值计算的机翼压力系数Cp分布对比见图8，考虑外翼段受三维效应影响较大，仅选取内翼段450 mm测压剖面结果进行对比。从压力分布对比中可以看出：仅在-2°攻角状态下，数值结果在机翼前缘顺压区域对压力分布捕捉与试验结果略有差别，但整体趋势一致；在机翼中段及后段，各攻角的结果对比都符合较好；整体来看，数值计算的压力分布与试验结果吻合良好，证明本文的求解器精度较高。

图8 数值计算与风洞试验压力系数分布对比

使用eN方法求解转捩位置，在稳定性分析阶段，由于试验构型机翼后掠角只有5°，故只考虑T-S波的增长，计算结果见图9，xtr/c表示弦向转捩位置。由于FL-3风洞为暂吹式风洞，湍流度较高为0.4%，因此转捩阈值较低。根据Mack公式计算得到的T-S波临界N因子为

图9 不同攻角下T-S波增长包络线

(7)

稳定性分析结果显示，NTScr=5对各攻角下转捩位置预测的精度较高。eN方法与试验结果的转捩位置对比见表2。可以看出，eN方法得出的转捩位置与试验结果趋势一致且误差较小，在0°攻角下误差最大，两者结果相差在3.3%c以内，符合工程应用的精度要求。

表2 eN方法与风洞试验转捩位置对比

综合压力系数分布特点与eN方法计算结果对风洞试验构型的层流特性进行分析。当攻角为-2°时，上翼面压力分布呈现较大的顺压力梯度，在45%c左右出现大的逆压梯度进行压力恢复。从T-S波扰动的增长来看，在大顺压力梯度区，T-S波的增长被有效抑制，当逆压出现时，T-S波迅速放大，直至发生转捩。0°攻角的压力分布与T-S波增长特征与-2°攻角基本一致。随着攻角的增大，上翼面顺压力梯度逐步减小，2°攻角下，机翼前缘出现了弱逆压梯度，T-S波出现了较快增长但并未达到转捩阈值，之后在20%c后压力分布转为顺压梯度，T-S波增长被抑制，而在45%c后的较大逆压梯度又使T-S波迅速放大，达到转捩阈值发生转捩。当攻角到达4°，上翼面压力分布形态已经转为明显的逆压梯度，T-S波也从前缘就迅速增长，并在20%c附近很快达到转捩阈值。

通过试验结果与数值计算结果的对比分析，可以看出，本文的数值方法在压力分布、转捩位置预测方面的模拟精度较高，可以反映风洞试验构型真实的转捩位置变化趋势。在多个攻角状态下，按照T-S波临界N因子为5所确定的转捩位置均与试验结果吻合较好，证明了利用湍流度和Mack公式确定风洞试验条件下T-S波临界N因子的可行性和本文采用的转捩预测方法的鲁棒性。

4 自然层流翼套飞行试验

4.1 基于某公务机飞行试验平台的层流翼套构型

为了对飞行试验条件下的自然层流进行研究，基于某公务机飞行平台，在其右侧机翼上设计改装了自然层流翼套，设计结果见图10。图中红色区域是设计的自然层流翼套，翼套区域占当地弦长的70%，展向宽度1 m，且展向翼型保持一致，两边的绿色区域是与主翼面的过渡区域，翼套的前缘后掠角与机翼相同，均为5°。相对于原始机翼，翼套在机翼前缘突出并包裹住原始机翼。图11给出了翼套与当地机翼的翼型对比。翼套外形与主翼面在后梁处相接。相比主翼原始翼型，翼套翼型的前缘半径较小，最大厚度后移，上表面最大弯度位置在50%弦长处，具有典型的自然层流翼型设计特点。

图10 某型公务机上自然层流翼套设计

图11 翼套翼型和当地原始翼型的对比

4.2 试验设备和数据测量

飞行试验需要同时获取翼套测试部分上表面的压力分布和转捩位置。压力分布是通过测压孔测量的。如图12所示，在测试段的上表面，沿流向分布有两排测压孔，以及倾斜的一排测压孔。分别标记为Sec1、Sec2和Sec3。

图12 试验段测压孔布置

转捩信息通过红外成像技术检测。如图13所示，将红外热像仪放置在机身上。在飞行测试期间，翼套通过安装在内部的加热丝加热，以增加测试部分和背景环境之间的温差。由于不同的热交换系数，层流和湍流表现出不同的温度分布。在图14的(红外)IR图像中，亮黄色区域表示层流，其他区域表示湍流。

图13 红外IR相机安装位置

图14 IR相机转捩成像图

流向测压孔的存在会触发从前缘附近的层流到湍流的过早转捩，从而导致表面上出现湍流楔形。因此，两个湍流楔形之间的区域用于转捩检测。在翼套的上表面上沿弦向方向设置8条等距的线。这些标记线可以帮助确定层流范围的长度。

5 飞行条件下层流翼套的转捩特性

5.1 飞行条件及计算模型的简化

由于对转捩的计算主要集中在层流翼套区域，出于节省计算资源的考虑，对翼套的计算构型做了部分简化。文献[26]表明，平、垂尾由于在机翼的尾流区域且距机翼较远，对机翼周围的流场影响几乎可以忽略，而吊舱距离机翼较近，对机翼表面的流场会产生较大的影响，因此将飞机的平、垂尾去除，仅保留机身、机翼、吊舱及其挂架进行层流翼套的转捩预测分析。计算使用的网格见图15，其中机翼环向布置201个点、展向布置113个点，在翼套部分进行了加密，全机网格量1 100万。

图15 层流翼套飞行构型计算网格

选定2个飞行试验工况，包含马赫数、飞行高度H、基于弦长的雷诺数Rec和攻角的不同变化。主要通过调节攻角来改变测试段的压力系数分布。飞行试验选定的测试点见表3。

表3 飞行试验工况

5.2 飞行条件下层流翼套的转捩特性

飞行条件下，该构型在各工况下翼套内侧和外侧剖面的压力系数分布计算结果如图16(a)和图17(a)所示。数值模拟与飞行试验结果在翼套头部附近由于加工精度和安装问题而有所偏差，但整体吻合较好。证明本文使用的数值模拟方法能够较好的预测和模拟翼套构型的压力分布特性。

对各状态下层流翼套的边界层进行稳定性分析，得出T-S波的扰动增长包络线云图见图16(b) 和图17(b)。图中红色线代表翼套构型的前缘，黑色虚线代表Sec1和Sec2两组测压孔，黑色实线代表湍流楔边界。考虑飞行条件下大气湍流度较低，转捩阈值较风洞试验要高。文献[35]对层流翼套飞行试验大气湍流度的测量结果显示，当高度范围为5.4～6.5 km，马赫数范围为0.45～0.57时，湍流度范围保持在0.034%～0.048% 之间。由于本文飞行试验高度和马赫数略高于此范围，且考虑到气象条件和发动机噪声等影响，来流湍流度最大可达到0.0725%[35]。根据Mack公式计算得到的T-S波临界N因子为

(8)

从稳定性分析的结果来看，在工况2状态下，由于上翼面具有较强的顺压梯度，T-S波的增长得到了较好的抑制，翼套上表面可保证45%c左右的层流区。在45%c后出现逆压恢复区，T-S波迅速增长。在工况1状态下，压力分布形态呈现一定的平台状，并在头部出现了小的逆压梯度，同时顺压梯度减小，对T-S波的抑制减弱，由于较大的飞行雷诺数，T-S波得到了较快增长，在20%c前即到达了临界值9。另外，工况1状态下从N因子NTS增长云图中可以看出翼套外侧(Sec2)的增长率较内侧(Sec1)小，这是由于外侧(Sec2)压力分布具有更大的顺压梯度，N因子NTS增长更慢。从IR图像中也可看出这一趋势。

层流翼套飞行试验选择马赫数范围为0.458～0.60和雷诺数范围为1 200万～1 600万。对比图16 (b)、图17(b)和图16(c)、图17(c)，从转捩位置的预测结果来看，在较大范围的马赫数和雷诺数变化条件下，按照T-S波临界N因子为9得到的转捩位置均与飞行试验结果吻合较好。进一步证明了利用湍流度和Mack公式确定飞行试验条件下T-S波临界N因子的可行性和本文使用的转捩预测方法的鲁棒性。

图16 工况1条件下转捩预测和飞行试验结果对比

图17 工况2条件下转捩预测和飞行试验结果对比

6 结 论

1) 在Ma=0.6、Re=3.37×106的风洞试验条件下，翼套风洞试验构型转捩位置的变化主要受攻角及压力分布形态特征影响。小攻角下(-2°～2°)翼套表面维持大范围的顺压力梯度，有效地抑制了T-S波的增长，最终转捩发生在翼型中后段的压力恢复区(50%c左右)。当攻角增大到4°，头部出现明显逆压梯度，T-S波扰动迅速放大，转捩提前至20%c。

2) 在Ma=0.458、Re=12.22×106的飞行试验条件下，由于上表面前缘存在弱逆压梯度，难以有效抑制T-S波增长，最终在20%c左右触发T-S波转捩；在Ma=0.60、Re=16.12×106的飞行试验条件下，受强顺压梯度影响，转捩位置可推迟至45%c左右。

3) 使用RANS方程耦合eN的转捩预测方法对层流翼套风洞试验和飞行试验的数值模拟在压力分布、转捩位置预测方面与试验数据吻合较好，表明本文数值方法有良好的预测精度和鲁棒性。

4) 通过湍流度结合Mack公式的方法确定T-S波临界N因子，在高湍流度低雷诺数的风洞试验条件下及低湍流度高雷诺数的飞行试验条件下均可提供较为准确的预测结果。也表明风洞试验结果能够在一定程度上指导飞行试验的实施。