0 引言

长期以来，由于风洞尺寸和试验条件的限制，风洞试验雷诺数远低于飞行雷诺数，对于大型飞机研制工作而言，试验雷诺数不足的问题更加突出。雷诺数的不同将对边界层流动产生显著影响，同时改变激波位置，造成不同雷诺数下飞机气动特性的明显差异^[1-3]。

与常规翼型相比，超临界翼型的上翼面平坦、下翼面后部内凹，这也就决定了超临界翼型典型的压力分布对雷诺数的影响较为敏感。雷诺数的差异对流动产生影响，导致风洞试验结果与飞行试验结果之间产生较大的差别^[4]。因此，对风洞试验数据必须进行雷诺数影响的修正。然而，与采用薄翼的战斗机仅需对风洞试验数据的阻力系数和最大升力系数进行修正不同，雷诺数变化对中等厚度和大厚度超临界机翼气动特性的影响十分复杂，影响阻力大小和最大升力系数的同时，也影响升力曲线斜率和俯仰力矩随迎角的变化规律^[5]，造成使用变雷诺数试验方法把试验数据外插到飞行值非常困难。

为了进行雷诺数的相关研究各航空大国都相继建立起高雷诺数风洞，如德国宇航研究院的KKK风洞、美国的NTF^[6]风洞，以及欧洲的ETW^[7]风洞。这些风洞试验取得了非常丰富的结果，为雷诺数效应研究奠定了基础，但是这些风洞无论是造价还是使用费用都非常昂贵。因此人们开始通过数值模拟方法来研究雷诺数效应。国外，Karl Pettersson^[8]等人开展了运输机雷诺数影响数值研究，Noah M. Favaregh^[9]等人通过数值试验和CFD方法开展了雷诺数效应研究，Mazuriah Said^[10]等人通过数值模拟方法研究了雷诺数对三角翼流场和气动特性的影响。国内张培红^[11]等采用CFD方法研究了雷诺数对飞机气动性能的影响，张辉^[12]等人开展了运输机雷诺数修正的研究，李忠武^[13]等开展了低速增升装置雷诺数计算研究，张耀兵^[14]等也开展了小展弦比模型的雷诺数影响的数值研究。通过这些研究进一步加深了人们对雷诺数效应的认识。然而这些研究的雷诺数范围不够大，不能反映真实飞行雷诺数状态，其次研究结果都是单方面的数值方法或试验方法，未能将两者有效结合、综合考虑，为这两种方法提供有效的综合结论。

本文综合采用数值模拟和风洞试验对比分析方法，研究了雷诺数效应问题。采用风洞试验结果结合CFD数值模拟将低雷诺数风洞试验数据修正至飞行雷诺数，并与高雷诺数风洞试验结果进行了对比，为大型飞机的研制提供有用的参考。

1 风洞试验描述

试验采用典型上单翼布局运输机模型，该模型在德国的ETW风洞进行了高速高雷诺数测力试验，部分试验车次采用如图 1所示的支撑方式。

ETW风洞是以氮气为试验介质的连续式回流跨声速风洞，通过低温和增压达到高的试验雷诺数。风洞的试验参数范围如下：

试验段尺寸：2.4 m×2.0 m

马赫数：0.15~1.3

总压：(1.25~4.5)×10⁵ Pa

总温：110~313 K

全模最大雷诺数(基于c=0.2 m)：50×10⁶

本次试验模型主要采用自然转捩的方式，只有低雷诺数试验项目另外开展了固定转捩试验。所有项目在机身头部都贴有转捩带。除了机翼和机头以外的其他模型部件都是自然转捩，这是为了避免不同试验条件下边界层过度转捩或转捩不足^[15]现象发生。机翼和机头的固定转捩带的尺寸规格、粘贴位置及使用范围如表 1所示。图 2给出了机翼、机头模型上贴有粗糙带的照片。

模型共进行了三个雷诺数的试验，分别为Re=3.3×10⁶(固定转捩和自由转捩)、25×10⁶和35×10⁶，参考弦长为c=0.2 m。翼身组合体构型的两个试验雷诺数Re=25×10⁶、35×10⁶是通过采用相同的气流总压和不同的气流总温实现的。对相同马赫数在不同雷诺数时保持了q/E相同，可以认为模型的气动弹性变形相当，因此在相同马赫数这两个雷诺数之间气动特性变化是纯粹的雷诺数效应。计算模型是绝对刚体，不存在变形的问题，三个雷诺数之间气动特性变化由雷诺数效应造成。本次试验雷诺数变化范围大，从常规的试验雷诺数到飞行雷诺数，为雷诺数影响分析奠定了基础。

2 数值模拟描述 2.1 数值计算概述

数值计算求解雷诺平均N-S方程，采用全湍流计算，不模拟流动转捩，与试验流动状态有所区别，在低雷诺数下区别明显，下文将详细说明。

2.2 控制方程

三维积分形式的雷诺平均N-S方程^[16]可写为：

$ \frac{\partial }{{\partial t}}\iiint {\mathit{\boldsymbol{Q}}{\text{d}}V} + \iint {\mathit{\boldsymbol{f}} \cdot \mathit{\boldsymbol{n}}{\text{d}}S} = 0 $

(1)

式中，V为控制体体积，S为控制体表面面积，Q为守恒变量，f为通过表面S的无黏通量和黏性通量之和，n为控制体表面S的外法向单位矢量。

以有限体积法构造空间半离散格式，无黏通量项采用二阶Roe迎风通量差分格式离散，黏性通量项采用中心差分格式离散，隐式时间推进，采用多重网格技术加速收敛^[17]。

2.3 湍流模型

采用当前通用的SST湍流模型进行湍流模拟。该模型将κ-ε和κ-ω模型进行调和，在固体壁面附近采用κ-ω模型，在自由流和边界层流外边界采用标准的κ-ε模型^[18]。SST模型吸收了κ-ε和κ-ω优点，其特点是：(1)加权模型系数；(2)限制快速应变流中的湍流黏性系数增长；(3)该模型对湍流黏性函数进行了修改，以加强对分离流的预测。SST湍流模型适用于中等分离的流动，尤其是在逆压梯度较大的流场模拟中表现出色。其湍动能方程和湍流频率方程的具体形式为：

$ \frac{{\partial \left( {\rho \kappa } \right)}}{{\partial t}} + \nabla \cdot \left( {\rho \mathit{\boldsymbol{u}}\kappa } \right) = \nabla \cdot \left[ {\left( {\mu + \frac{{{\mu _t}}}{{{\sigma _k}}}} \right)\nabla \kappa } \right] + {P_\kappa } - \beta '\rho \kappa \omega $

(2)

$ \begin{gathered} \frac{{\partial \left( {\rho \omega } \right)}}{{\partial t}}\nabla \cdot \left( {\rho \mathit{\boldsymbol{u}}\omega } \right) = \nabla \cdot \left[ {\left( {\mu + \frac{{{\mu _t}}}{{{\sigma _\omega }}}} \right)\nabla \omega } \right] + \hfill \\ 2\rho \left( {1 - {F_1}} \right)\frac{1}{{{\sigma _{\omega 2}}\omega }}\nabla \kappa \nabla \omega + \alpha \frac{\omega }{\kappa }{p_\kappa } - \beta \rho {\omega ^2} \hfill \\ \end{gathered} $

(3)

$ {P_\kappa } = - \rho \overline {{{u'}_i}{{u'}_j}} \nabla \cdot \mathit{\boldsymbol{u}} $

(4)

模式中参数Φ的值是两组参数值的“混合”，由下式确定：

$ \mathit{\Phi } = {F_1}{\mathit{\Phi }_1} + \left( {1 - {F_1}} \right){\mathit{\Phi }_2} $

(5)

第一组参数用于κ-ω模型：

β′=0.09, α₁=5/9, β₁=0.075，

σ_k1=2, σ_ω1=2, κ=0.41。

第二组参数用于κ-ε模型：

β′=0.09, α₂=0.44, β₂=0.0828,

σ_k2=1, σ_ω2=1/0.856, κ=0.41。

湍流黏性系数${\mu _t} = \frac{{{\alpha _1}\rho \kappa }}{{\max \left( {{\alpha _1}, S{F_2}} \right)}}$，

函数F₁、F₂定义为

$ {F_1} = \tanh \left( {\mathit{arg}_1^4} \right), {F_2} = \tanh \left( {\mathit{arg}_2^2} \right)，$

其中，

$ \begin{gathered} \mathit{ar}{\mathit{g}_1} = \min \left( {\max \left( {\frac{{\sqrt \kappa }}{{\beta '\omega y}}, \frac{{500v}}{{{y^2}\omega }}} \right), \frac{{4\rho \kappa }}{{C{D_{\kappa \omega }}{\sigma _{\omega 2}}{y^2}}}} \right), \hfill \\ \mathit{ar}{\mathit{g}_1} = \max \left( {\frac{{\sqrt \kappa }}{{\beta '\omega y}}, \frac{{500v}}{{{y^2}\omega }}} \right), \hfill \\ C{D_{\kappa \omega }} = \max \left( {2\rho \frac{1}{{{\sigma _{\omega 2}}\omega }}\nabla k\nabla w, 1.0 \times {{10}^{ - 10}}} \right), \hfill \\ \end{gathered} $

y是点到固壁的距离，ν是运动黏性系数。

2.4 计算模型与计算网格

本文的数值模拟采用了多块对接结构网格技术，计算区域的远场边界向前、向上、向下和展向方向取为80倍平均气动弦长，向后取140倍。壁面的第一层网格达到了0.5×10^-6m，增长比率为1.2左右，网格在各个剪切层附近均进行了适当的加密，以保证附面层内和剪切层的数值模拟精度，半模网格规模2235万。图 3给出了物面网格的局部图。

3 结果分析

根据风洞试验状态，确定数值计算状态为：Ma=0.76。分别模拟了Re=3.3×10⁶、25×10⁶和35×10⁶三个雷诺数状态。对试验构型Re=3.3×10⁶的固定转捩状态、Re=25×10⁶和35×10⁶的自然转捩状态进行深入研究，主要是因为这些状态下能够保证机翼前缘或前缘附近是湍流边界层，有利于和全湍流计算结果进行对比。

图 4为Ma=0.76时不同雷诺数下翼身组合体构型风洞试验结果与计算结果的对比。可以看出，计算结果和试验结果趋势一致，特别是在高雷诺数下吻合良好：雷诺数从3.3×10⁶增大到25×10⁶时，零升迎角减小，升力线斜率增大，升力线拐点相当，最大升力系数增大；雷诺数从25×10⁶增大到35×10⁶时，最大升力系数略有增加，升力线斜率略有增大；随雷诺数增大，阻力有所减小，最大升阻比增大；随雷诺数增大，低头力矩增大，力矩曲线拐点对应的升力系数增大。

从差量上来看，雷诺数从3.3×10⁶增大到25×10⁶，最小阻力系数风洞试验结果减小约0.0043，计算结果为0.00424，试验结果低头力矩增大约0.026，计算结果增大0.016；雷诺数从25×10⁶增大到35×10⁶，最小阻力系数风洞试验结果减小约为0.0005，计算结果为0.00056，试验结果低头力矩增大约0.002，计算结果增大0.0021；雷诺数从3.3×10⁶增大到35×10⁶，试验结果最大升阻比增大约3.8，计算结果增大3.83。

计算结果和试验结果之间的差异在于：虽然升力线斜率一致，升力线拐点相当，但是相同迎角下计算结果升力系数略大，随着雷诺数增大，相同迎角下升力系数差量逐渐减小；力矩曲线斜率和力矩曲线拐点对应的升力系数基本一致，但在大升力系数下，试验的力矩曲线斜率变大，计算的力矩曲线斜率略有减小，两者规律出现差异，可能与翼面的分离有关，有待进一步研究验证。

从以上分析来看，对于湍流边界层的机翼，试验所得的机翼的气动力系数随雷诺数的变化趋势与计算的趋势是一致的，量值也较为接近，计算结果在这种情况下具有很高的参考价值。

图 5给出了不同雷诺数下机翼剖面压力分布比较。在相同迎角下，雷诺数增加使得机翼翼面压力分布形态有所变化。雷诺数增加，机翼上表面负压峰值增加，中后部负压略微提高。激波的弦向位置稍微靠后，激波强度略为增加。由于外翼的当地雷诺数较内翼小，所以外翼变化比内翼较明显。机翼下表面压力值提高。后缘附近压力值提高，有助于减弱后缘附近的附面层分离，从而减小了压差阻力。

在相同升力系数下，雷诺数增加，机翼上表面负压峰值减小，中后部负压略微提高。激波的弦向位置略微前移，强度减小。机翼下表面压力值提高。后缘附近压力值提高。但总体而言，雷诺数变化使相同升力系数下翼面压力分布形态变化很小。由此可说明，小雷诺数风洞测压结果与实际飞行中的压力分布在相同升力系数时基本一致。

4 巡航使用范围阻力和升阻比雷诺数修正研究

从第3节可知，在前缘固定转捩或者湍流边界层的前提下，计算可准确模拟升力、阻力、俯仰力矩随雷诺数变化的规律性，变化量值也与风洞试验较为接近，可利用风洞试验结果结合CFD手段将低雷诺数风洞试验数据修正至飞行雷诺数。具体方法如下：

1) 进行风洞试验雷诺数下的飞机气动特性数值模拟；

2) 进行飞行雷诺数下的飞机气动特性数值模拟；

3) 将数值模拟得到的飞行雷诺数与风洞试验雷诺数下气动特性差量叠加到风洞试验数据，即可得到飞机飞行雷诺数下的气动特性数据。本文采用该方法对雷诺数3.3×10⁶试验结果进行修正，试图得到雷诺数25×10⁶和35×10⁶下的相应数据，并且与相同雷诺数下的试验数据进行对比。图 6为修正数据阻力极曲线与试验结果对比，在巡航升力系数范围(C_L=0.45~0.55)，阻力系数C_D修正结果与试验结果的相差不超过0.0004。图 7为修正数据升阻比曲线与试验结果对比，升阻比K修正结果与试验结果的最大差量约为0.2。

图 8给出了力矩系数的修正结果对比，可见经修正后的结果精度很差，误差达到0.01以上，不能满足飞机设计对力矩系数的精度要求。

经过分析，升力系数差量也会对力矩产生影响，因此对于力矩系数，在考虑力矩数据差量的基础上，需要进一步考虑升力系数差量对力矩的影响。本文采用下式进行修正。

$ \begin{gathered} {C_{m{\text{flight}}}} = {C_{m\exp }} + \left( {{C_{m{\text{cal - flight}}}} - {C_{m{\text{cal - exp}}}}} \right) + \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\left( {{C_{L{\text{cal - flight}}}} - {C_{L{\text{cal - exp}}}}} \right)d/L \hfill \\ \end{gathered} $

(6)

d为焦点到力矩参考点的距离，亦即静稳定性裕度$\partial {C_m}/\partial {C_L}$。为了对比分析，分别选取了两个静稳定性裕度系数，试验状态和飞行计算状态。

图 9给出了力矩改进修正的结果，其中Exp._Re3.3M+ (Cal._Re25M-Cal._Re3.3M)_1为采用试验状态纵向静稳定性裕度的结果，Exp._Re3.3M+(Cal._Re25M-Cal._Re3.3M)_2为采用计算状态纵向静稳定性裕度的结果，可见修正精度大幅提高，误差降为0.001。对于雷诺数25×10⁶，采用计算纵向静稳定性裕度修正结果精度在线性段较高，在非线性段与试验纵向静稳定性裕度修正结果更接近。对于雷诺数35×10⁶，试验纵向静稳定性裕度修正结果更接近真实试验值。

5 结论

本文通过试验和数值模拟对大型飞机的雷诺数影响规律进行相关研究，讨论了结合数值模拟手段对风洞试验数据进行雷诺数修正的方法，得到以下结论：

1) 对于湍流边界层的机翼, 试验所得的机翼的气动力系数随雷诺数变化的趋势与计算的趋势是一致的，量值也较为接近；

2) 雷诺数增大对大型飞机气动特性影响显著，从本文的试验和数值模拟结果来看，相同迎角下升力系数和升力线斜率增大，阻力系数减小，俯仰力矩系数增大；

3) 小雷诺数风洞测压结果与实际飞行在相同升力系数时基本一致；

4) 采用风洞试验结果结合CFD手段将低雷诺数风洞试验数据修正至飞行雷诺数，修正结果与相应雷诺数下的试验结果相比，阻力系数相差不超过0.0004，升阻比最大差量约为0.2，力矩系数最大差量约为0.001。该修正方法为未来风洞试验数据应用于飞机设计、计算提供了一种可基本满足要求的雷诺数修正方法。