自然层流飞行测试翼套的仿真和试验

王浩，钟敏，华俊，钟海，杨体浩，王猛，雷国东

1. 中国航空研究院 技术研究二部，北京 100012

2. 中国飞行试验研究院 中航工业飞行仿真航空科技重点实验室，西安 710089

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

4. 航空工业空气动力研究院 高速高雷诺数航空科技重点实验室，沈阳 110034

绿色发展已经成为航空业发展的趋势和要求，国际航空运输协会(The International Air Transport Association, IATA)在2020年初发布了《2050飞机技术路线图》，概述和评估了未来可能的绿色航空技术，层流控制技术被认为是最有前景的技术之一，并有望达到较高的技术成熟度[1]。跟据对大量计算和飞行试验的结果分析，摩擦阻力可占民机总阻力的50%，层流流动的摩擦阻力远小于湍流流动，因而采用层流控制技术扩大层流区是减阻的一个重要途径[2]。在层流边界层控制的研究中，逐步形成了自然层流控制、全层流控制和混合层流控制等3种技术[3-4]，其中自然层流技术是一种被动控制技术，不需要引入额外的能量，适用于机翼后掠角较小的飞行器中。

针对层流控制技术的发展，国内外已经进行了大量的数值仿真研究，目前常用的转捩预测方法有eN方法[5-7]，经验关系式法[8]、低雷诺数湍流模型[9]和γ-Reθ t转捩模型[10-12]，其中eN方法基于线性稳定性理论，着重于从物理上描述边界层中小扰动TS(Tollmien-Schlichting)波沿流向的线性放大过程，并根据经验选定判定转捩发生的放大因子临界值，可以预测低湍流度下沿流向的自然转捩和分离流转捩，在工程实践中的应用较为广泛。

在自然和混合层流机翼方面国外开展了大量的风洞和飞行试验，并进行了一定的商业化尝试。美国国家航空航天局和波音公司于2018年基于CRM(Common Research Model)标模分别在美国NTF(National Transonic Facility)风洞和欧洲ETW(European Transonic Wind Tunnel)风洞开展中等后掠角机翼的自然层流试验，研究了不同雷诺数下TS波和CF(Cross Flow)波对转捩的影响[13-14]。ATR(Avions de Transport Regional)公司和法国航空航天中心合作开展了未来支线飞机低阻层流机翼的试验测试，在S1MA风洞中进行了雷诺数接近于真实飞行状态的风洞试验[15]。波音公司曾在B757飞机上开展了混合层流减阻研究[16]，NASA开展的SARGE(Subsonic Aircraft Roughness Glove Experiment)项目，在G-II飞机上进行了亚声速飞机离散粗糙元翼套试验研究[17-18]。欧盟在SAAB2000飞机的机翼上进行了混合层流翼套控制飞行试验[19]，并在A320垂尾上进行了混合层流控制的飞行验证[2]。空客在BLADE(Breakthrough Laminar Aircraft Demonstrator in Europe)项目中，采用A340作为载机，将外翼段替换为后掠角较小的自然层流机翼，评估大型客机上引入层流机翼技术的可行性。Honda Jet轻型公务机采用自然层流机翼，在2003年进行了首飞，达到了预期的设计目标[20-21]。国内张彦军[22]、艾梦琪[23]等开展了一系列风洞试验研究，但国内之前尚未开展层流机翼的飞行试验，对相关数据的对比分析研究也较为罕见，特别是缺少对仿真、风洞试验和飞行试验三者间的对比分析。

近期，国内层流翼套飞行测试团队进行了自然和混合层流机翼翼套的设计、数值仿真、风洞试验、飞机改装和飞行试验的相关探索并陆续发表有关研究内容[24-26]，数值仿真可以为试验研究提供预评估结果。本文主要针对该项目中的自然层流机翼翼套的数值仿真、风洞试验和飞行试验的结果开展对比分析，并且通过数值分析结果与试验结果、特别是中国首次自主开展的层流飞行实测数据的对比和数据相关性分析，验证本项目的数值方法，为今后开展自然层流机翼的仿真、设计和试验研究及工程化开发提供参考。

1 飞行试验机、风洞试验模型和数值仿真方法

1.1 飞行试验机和风洞试验模型

飞行试验使用了某小型民用飞机作为原型机，为直线型下单翼，后掠式垂直尾翼带大背鳍，发动机安装形式为尾吊，便于翼套的安装和转捩的测量。在原型机的单侧机翼上安装自然层流翼套，翼套的设计基本保持两侧机翼载荷的平衡。加装翼套的层流机翼翼套飞行试验机半模如图1所示，翼套的控制面位于展向3.36 m和展向4.36 m，翼套展向长度为1 m，展向3.36 m位置处的弦长为2.37 m，展向4.46 m位置处的弦长为2.12 m，相较原始机翼前伸0.25 m；过渡段的两个控制剖面分别位于2.91 m和4.69 m，内过渡段长度为0.45 m，外过渡段长度为0.33 m，翼套与过渡段的长度之和为1.78 m，约占半展长的21%。翼套后掠角与机翼外翼段后掠角一致均为5度，试验中翼套的当地雷诺数约为12×106～18×106，转捩主要受TS波扰动造成，横流的影响较小。

图1 自然层流机翼翼套飞行试验机

层流机翼翼套采用复合材料加工，通过红外热成像的方式进行边界层转捩的测量，试验过程中通过表面加热的方式获得机翼上表面质量较高的红外热图[25]。飞行测试翼套如图2所示。

图2 自然层流机翼翼套测试区域

在展向3.623 m位置处从前缘至60%当地弦长布置19个测压点，记为S3测压剖面；在展向4.303 m处从前缘至60%当地弦长布置18个测压点，记为S2测压剖面；考虑到测压孔的扰动会引起转捩，在S3测压剖面从16.2%当地弦长到60%当地弦长沿斜向与S3剖面成18°角布置了9个测压点，记为S1测压剖面。由于更换混合层流控制的吸气前缘所需，层流翼套在18%当地弦长处，存在一条接缝，通过手工处理使其表面粗糙度和波纹度尽量满足自然层流的要求。根据试验机的大气数据系统测量所得，自然层流机翼翼套飞行试验的部分典型工况如表1所示。

表1 自然层流机翼翼套飞行试验典型工况

在飞行试验前，截取翼套中间部位翼型开展了风洞试验研究，风洞试验模型机翼翼型与自然层流翼套S2和S3测压剖面翼型的对比如图3所示(为便于对比，厚度方向有放大)，为适应机翼扭转角，内侧的S3测压剖面的当地攻角略大于外侧的S2测压剖面。设计加工了简单后掠机翼的风洞试验模型(图4)。

图3 自然层流机翼翼套翼型

图4 自然层流机翼翼套风洞试验模型

风动试验模型机翼展长为900 mm，弦长为300 mm，前缘后掠角为5°。为与飞行试验和仿真结果更好的进行对比，在展向450、600、750 mm 处分别沿斜向与机身轴线成18°角布置测压孔。该模型在FL-3风洞开展了自然层流风洞试验，在内外测压剖面之间的区域采用红外热图的方式进行转捩测量，风洞试验工况如表2所示。

表2 自然层流机翼翼套风洞试验工况

1.2 数值仿真方法

数值仿真流程采用的计算流体力学(CFD)求解器为中国航空研究院的雷诺平均Navier-Stokes (RANS)求解器AVICFD-Y，采用的湍流模型为SST。

进行转捩分析时采用的是基于线性稳定性理论的eN方法，结合CFD求解器得到的机翼截面几何信息、压力信息、当地后掠角、马赫数(Ma)、雷诺数，基于锥形流假设，求解层流边界层方程，获得边界层的基本解[27]，进而求解四阶可压缩Orr-Sommerfeld方程进行稳定性分析，求解得到扰动放大因子(N)，当N达到经验所得的转捩阈值(Ncr)时可以认为发生转捩，从而获得转捩位置，计算中风洞试验的Ncr定为6.5，而飞行试验的Ncr定为10；将求解得到的转捩信息通过强制给定转捩位置的方式代回到CFD求解器中，重新进行流场计算获得压力信息，并重新进行转捩分析获得新的转捩位置，直到转捩位置的变化满足收敛条件为止[6]。计算流程如图5所示。

图5 计算流程图

数值计算采用结构网格，飞行试验机网格节点数为1 000万，边界层内第一层网格高度为0.000 001 m，增长率为1.15，机翼沿弦向一周的网格节点数为201，飞行试验机的计算模型网格如图6所示；风洞试验模型网格节点数为400万，边界层内第一层网格高度为0.000 002 m，增长率为1.2，机翼沿弦向一周的网格节点数为265，风洞试验模型的计算模型网格如图7所示。计算状态以翼套的设计状态马赫数0.6，攻角0°为主，兼顾其他试验和飞行状态。

图6 飞行试验机计算模型网格

图7 风洞试验模型计算模型网格

2 数值仿真和风洞试验结果对比

2.1 压力系数分布

机翼各截面数值仿真与风洞试验的压力系数(Cp)分布对比如图8和图9所示，其中Pt为总压，风洞测压结果与数值仿真的一致性较好。

图8 Pt=0.1 MPa、Ma=0.6时机翼截面压力系数分布

图9 Pt=0.2 MPa、Ma=0.6时机翼截面压力分布

2.2 转捩位置

风洞试验经过处理可以得到规整的红外热图，利用梯度算法计算热图的灰度梯度场，并计算弦向和展向的合成梯度绝对值，利用梯度最大值判断每个展向剖面的转捩位置，得到一组转捩点站位数据。利用统计方法得到这些转捩点站位在整个弦向坐标(0～1.0)内的概率密度分布曲线，根据最大概率密度指示转捩位置。根据数值仿真分析的结果，在中间翼段沿展向转捩点的弦向位置变化不大，以y=600 mm的截面为例，对比不同状态时数值仿真和风洞试验机翼转捩位置(xtr/c)随攻角的变化曲线，如图10所示。各种工况下，在小攻角和大攻角时，转捩位置随攻角变化较小，数值仿真得到的转捩位置较风洞试验结果靠前，更接近翼套设计转捩位置。该仿真方法在风洞状态小雷诺数时会得到层流分离泡，从而转捩预测结果较风洞试验靠前。随着攻角的增大，从某个攻角开始转捩位置会迅速前移，数值仿真得到的转捩位置开始迅速前移的攻角较风洞试验结果更大。对比不同风洞总压即不同雷诺数时的结果，可以发现当总压即雷诺数增大后，在较小攻角时风洞试验测得的转捩位置就开始较为显著的随攻角增大而前移，在0°攻角时已前于飞行雷诺数的设计转捩位置。风洞试验和数值仿真两者转捩位置变化的趋势存在一定差异，原因可能为转捩位置对风洞品质、表面质量等因素较为敏感，导致风洞试验中的转捩位置较仿真结果在更小的攻角时即发生前移，增压后风洞品质进一步的变化促使转捩位置前移的攻角进一步减小，目前仿真分析方法尚无法对暂冲式增压风洞紊流度进行准确的模拟。后续对红外热图的分析也可以得到类似结论，有必要进一步进行研究确认。

图10 转捩位置随攻角的变化曲线

对比Ma=0.6、攻角=0°时不同总压风洞试验得到的红外热图与数值仿真得到的摩阻云图，如图11所示，不同总压下仿真得到的转捩位置基本一致，且沿展向基本无变化。对于风洞试验得到的红外热图，机翼表面存在较多湍流楔，主要为模型表面测压孔和灰尘等引起，在总压为0.1 MPa 时可以观察到明确的转捩线，获得可靠的转捩位置；在增压至0.2 MPa时机翼表面湍流楔迅速增多，推断此时试验测得的转捩位置较早前移可能与风洞增压后品质改变有关，有待进一步研究确认。

图11 风洞试验红外热图与数值仿真摩阻系数云图对比

3 数值仿真和飞行试验结果对比

3.1 压力系数分布

飞行试验中飞行器的气动力需要通过换算得出，因此通过对比压力系数分布可以更好的保证数值仿真与飞行试验状态的一致性，但飞行试验中压力系数分布的获取需要确定更为准确的参考静压。试验过程中在机身最大截面处侧面舷窗开孔测量的机身静压和机头处空速管测量的机头静压，作为对翼套压力数据进行无量纲化的参考压力，同时根据飞行试验机的大气数据系统测得的气压高度结果可以求得相应的理论静压，也可以作为参考压力。几种不同参考压力的对比如表3所示,Ma=0.6、攻角=0°时机身静压与理论静压较为接近，但Ma=0.46，攻角=2°时机身静压减小较为明显；机头静压在攻角变化时相对稳定，显著大于机身静压和理论静压。

表3 不同方法得到的参考压力

以外侧的S2测压剖面为例，对比不同静压值作为参考压力时仿真和飞行试验的压力系数分布，如图12和图13所示，其中H表示飞行高度。在Ma=0.6，攻角=0°的设计状态下，机身静压和理论静压作为参考压力，数值仿真和飞行试验的压力系数分布均可以取得较好的一致性，采用机头静压作为参考压力时，仿真和飞行试验的压力系数分布差异较大。当Ma减小到0.46，攻角增大为2°时，相较于机头静压和机身静压，理论静压作为参考压力时飞行试验和数值仿真的压力系数分布的一致性更好。仿真与飞行试验的攻角差值约为0.3°。

图12 H=7 000 m、Ma=0.60、攻角=0°时数值仿真与飞行试验S2截面压力系数分布

3.2 转捩位置

以H=7 000 m、Ma=0.60、攻角=0°，H=7 000 m、Ma=0.46、攻角=2°和H=6 000 m、Ma=0.60、攻角=0°的状态为例，飞行试验机翼上表面红外热图(左)和数值仿真摩阻系数云图(右)的对比如图14所示。两图中黑点及虚线从前缘向后分别标记0.05、0.16、0.20、0.25、0.30、0.35、0.40和0.45倍当地弦长的位置，红外热图中蓝色实线标识出飞行试验测量得到的转捩位置。可以发现，在Ma=0.6，攻角=0°的设计状态，飞行中S2、S3截面前缘处的测压孔诱导出两个湍流楔，在湍流楔的影响区域之间存在大范围的层流区，此时仿真结果和试验结果一致性很好，转捩均发生在约0.45倍弦长位置附近；在Ma=0.46，攻角=2°的非设计状态，转捩发生在接缝之前，此时仿真结果和飞行试验结果也较为接近。对于图14(c) 所示的设计马赫数和攻角，当飞行高度降低为6 000 m时，雷诺数增大到17.8×106，此时层流机翼翼套仍然可以保持0.45倍弦长左右的层流区，同时也进一步说明暂冲式风洞在增压时的流场品质对转捩位置测量具有一定影响。

图14 飞行试验红外热图与数值仿真摩阻系数云图对比

选取翼套展向y=3.932 m位置的中间截面记为T1截面，如图14(b)右图所示，分析其扰动放大因子的变化过程，不同状态时的结果如图15所示。对应不同高度即不同雷诺数，Ma=0.6，攻角=0°的设计状态下，受较大的顺压梯度的影响，TS波扰动放大因子(NTS)增长均较慢，直到0.45倍弦长之后，较大的逆压梯度才导致其迅速增长引发转捩；此时在上翼面顺压区NTS均较小，多处于2以下。在马赫数0.46，攻角=2°时(图15(b))，在前缘附近NTS即迅速增长导致在接缝前发生转捩，稳定性分析的结果与飞行测量也有很好的一致性。结合图12～图15结果进行分析，攻角=0°时顺压较为显著，此时扰动的放大速度较慢，故图14(a)和图14(c)中翼套上表面沿展向均在45%弦长附近发生转捩；攻角=2°时顺压梯度较小，此时扰动的放大速度较快，转捩对压力分布较为敏感，因此图14(b)中数值仿真得到的翼套上表面转捩位置会沿展向变化。

图15 T1截面NTS变化过程

3.3 截面阻力

选取部分自然层流区较长的飞行试验状态，在T1截面后方布置尾流耙测量截面阻力，分别进行自然转捩和6%弦长处强制转捩的飞行试验，研究自然转捩技术的减阻效果，涉及的飞行试验状态如表4所示。

表4 截面阻力系数测量飞行试验工况

根据数值仿真结果进行分析时，考虑到翼套仅对上翼面进行了特殊处理，进行红外测量，因此仅在翼套上翼面进行自然转捩分析，其余部分均按照全湍进行处理。采用尾迹动量损失积分法计算阻力，阻力系数计算公式为

(1)

式中：c为截面弦长；Cf为尾迹截面单位长度上的无量纲动量损失，具体可以参考文献[28]。为与飞行试验保持一致，在截面后方0.12倍弦长的位置截取数据求解截面阻力，数值仿真结果与飞行试验结果的对比如表5所示。从表中可以看出，同样状态时，数值仿真的自然转捩和强制转捩的截面阻力系数差值在0.001 2～0.001 6之间，飞行试验结果在0.001 3～0.001 6之间，二者结果一致性很好，差别约为0～0.000 2。

表5 飞行试验与数值仿真截面阻力系数减小量

4 结 论

1) 建立了自然层流机翼数值仿真的CFD-流动稳定性分析流程，针对自然层流机翼翼套飞行试验机和缩比风洞试验模型开展了设计验证和典型试验状态仿真，多数状态下数值仿真和风洞试验的压力系数分布具有较好的一致性，和飞行试验压力分布的对比取得了更好的一致性。

2) 研究中发现，风洞总压为0.1 MPa时，雷诺数为300万量级，数值仿真得到的转捩位置随攻角的变化规律与风洞试验结果一致，转捩位置相对风洞试验靠前，更接近翼套的设计转捩位置；风洞增压为0.2 MPa时，雷诺数增大为600万量级，风洞来流紊流度也随之增大，机翼表面出现大量的湍流楔，影响了风洞试验对转捩位置的测量。

3) 对于飞行测试的压力分布，对比了不同参考静压源的结果：相较于机头处静压和机身最大截面处静压，飞行试验机大气数据系统测量的气压高度对应的理论静压作为参考压力时得到的压力系数分布与数值仿真结果更为一致，作为对比数据更为合理，此时数值仿真得到的压力系数分布和飞行试验结果具有良好的一致性。

4) 在翼套设计马赫数0.6，攻角0°附近，数值仿真得到的转捩位置和飞行试验结果高度一致；在雷诺数达到1 800万量级时机翼上翼面也可以实现45%左右的层流区；在较低的马赫数较大攻角时，仿真预测的转捩发生在接缝之前，也与飞行测试重合良好。

5) 同样状态下数值仿真计算得到的采用层流技术带来的截面阻力减小量与飞行测量结果基本一致，两者差别在0～0.000 2之间。