高速风洞摆动叶片式阵风发生器非定常流场数值模拟与试验验证

2022-12-09张颖刘南

实验流体力学 2022年5期

张颖，刘南

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

0 引言

现代飞机，如军用和民用运输机、高空长航时无人机、远程轰炸机等，为了追求较高的气动效率和较低的结构重量，通常采用大展弦比机翼，在结构材料上大量使用复合材料，具有较大的柔性，这使得飞机在跨声速阶段的阵风问题凸显。阵风载荷，特别是垂直离散阵风载荷，有可能成为飞行载荷中最严重的载荷。因此，准确预测阵风载荷并采取一定措施进行减缓控制通常是飞机设计阶段的一项重要工作[1-5]。

阵风风洞试验是进行飞行器阵风响应特性和阵风载荷减缓控制技术研究的重要手段之一，而阵风发生器是开展阵风试验的重要设施。在高速阵风发生器的设计和研制方面，具有代表性的英国ARA 和美国NASA 均在其大型连续式跨声速风洞上配备了阵风发生器，如英国ARA TWT 风洞的尾缘吹气式阵风发生器[4]以及美国TDT 风洞的摆动叶片式阵风发生器[5]。尾缘吹气式阵风发生器无运动部件，无惯性载荷，控制精度高，可以产生频率较宽和形式多样的阵风，但由于其尾缘射流通过大规模电磁阀进行控制，叶片内部结构复杂，控制系统庞大，代价高，实现难度较大。摆动叶片式阵风发生器机械结构和控制形式相对简单，代价较低，可以通过叶片的机械摆动在叶片下游产生正弦周期性变化的阵风速度场。利用该阵风发生器，TDT 风洞针对多种型号飞机（如B-52E 飞机、C-5A 飞机、传感器飞机以及超声速民机等）开展了大量的阵风响应及减缓控制研究工作，为飞机设计过程中阵风响应和载荷减缓控制律设计奠定了良好的基础[5-10]。

国内在阵风试验技术研究方面起步相对较晚，且大多数研究局限于低速领域[11-14]，大型高速风洞没有配备阵风发生器，尚不具备阵风试验能力，而高速巡航阶段的阵风响应对飞行安全影响极大，因此有必要开展高速风洞阵风试验能力建设。

阵风发生器作为阵风试验系统的重要组成部分，其设计难点主要在于：1）阵风发生器形式选择；2）阵风发生器气动设计。影响阵风发生器气动性能的参数众多，如叶片截面形状、弦长、展长、数量、间距等，需要进行多参数优化迭代以使其性能达到最优，同时还要兼顾叶片的刚、强度要求。本文在综合考虑国外高速风洞2 种常用形式阵风发生器优缺点的基础上，采用自研ENSMB 流场计算软件对摆动叶片式阵风发生器设计方法进行了探索性研究。依托中国航空工业空气动力研究院0.6 m 连续式跨声速风洞（以下简称FL-61 风洞）建立了一套摆动叶片式阵风发生器，并进行了阵风流场校测试验，分析了该阵风发生器下游阵风速度场形成机理及分布特性，重点开展了叶片摆动频率和最大摆动幅值等参数对叶片下游阵风速度幅值影响规律研究。

1 摆动叶片式阵风发生器

摆动叶片式阵风发生器采用单叶片等直机翼形式，叶片垂直安装于风洞喷管入口的上壁面（图1）。叶片截面形状为NACA0012 翼型，弦长0.2 m，展长0.2 m（图2）。通过电机驱动曲柄摇杆机构带动阵风发生器叶片绕1/4 弦线做周期性正弦运动（图3），进而在叶片下游产生正弦形式的周期性阵风速度场。叶片能实现的摆动频率f 和最大摆动幅值Amax范围分别为f=0～20 Hz，Amax=−20°～20°。数值计算模型包括风洞收缩段、喷管段、试验段以及阵风发生器叶片，如图4 所示。

图1 阵风发生器在风洞中的安装示意图Fig. 1 Gust generator mounted in the tunnel

图2 阵风发生器结构示意图Fig. 2 Structure of gust generator

图3 阵风发生器驱动机构示意图Fig. 3 Driven mechanism of gust generator

图4 计算模型Fig. 4 Computational model

2 数值模拟方法

2.1 计算网格

计算采用结构化N–S 网格进行，网格单元总数约407 万，节点总数约412 万。叶片周围边界层第一层网格高度为叶片弦长的十万分之一。为了更好地捕捉叶片后方阵风场，对叶片后方重点关注区域的网格沿流向进行了局部加密，如图5 所示。图中：x 轴沿气流方向，向后为正；y 轴沿叶片展长方向，向上为正；z 轴符合右手坐标系定则；坐标原点为喷管入口端的风洞中心点。

图5 计算网格Fig. 5 Computational grid

2.2 计算方法

本文采用中国航空工业空气动力研究院自研的ENSMB 流场计算软件[15-16]求解三维可压缩非定常N–S 方程，湍流模型选用k–ωSST 模型，空间方向采用二阶迎风格式进行离散，时间方向采用双时间方法进行离散，时间步长的选取与叶片摆动频率相关，一个周期T 内进行1000 次外迭代，计算总时长为5 个周期。

图6 为来流马赫数Ma=0.5，f=4 Hz，Amax=6°工况下，采用不同内迭代步计算得到的阵风速度特性对比。图中t 为计算时间，v 为阵风速度，I 为内迭代步，vgust为阵风速度幅值。阵风速度幅值vgust的定义如下：

图6 不同内迭代步下阵风速度特性对比Fig. 6 Comparison of gust velocity characteristics with different iteration steps

式中，vmax和vmin分别为一个周期内阵风速度的峰值和谷值。从图中可以看出，在时间步长保持一致的条件下，不同内迭代步计算得到的阵风速度幅值几乎无差异。考虑计算效率问题，本文的内迭代步设置为10 步。

计算时的边界条件设置如下：风洞入口设置为总温总压边界条件；风洞出口设置为压力出口边界条件；风洞壁面和阵风发生器叶片均设置为黏性无滑移壁面。阵风发生器叶片绕1/4 弦线位置按正弦规律做周期性运动，其方程为：

式中，A（t）为叶片的瞬态摆角。

叶片的运动通过动网格方式实现，为了提高网格的变形效率和质量，采用Radial Basis Function（RBF）和Linear Transfinite Interpolation （TFI）混合方式进行。

3 结果及分析

本文重点开展了叶片摆动频率和最大摆动幅值对叶片下游阵风速度幅值的影响规律研究，并与试验结果进行了对比分析。计算状态设置为Ma=0.5，Amax=2°～18°，f=2～20 Hz。

为便于分析叶片后方阵风场，在叶片后方选取124 个监测点（沿x 轴布置4 列，沿y 轴布置31 行），具体位置如图7 所示。监测点的x 轴坐标分别为x=2.75 、3.55 、4.25 、4.75 m，监测点的y 轴坐标位于–0.30～0.30 m 之间，间隔0.02 m。监测点编号记为Pij，下标i 为监测点所处的列号，下标j 为监测点所处的行号。为与试验进行对比，本文仅对图7 中第二列（x=3.55 m）监测点位置处的阵风速度特性进行分析。

图7 监测点示意图Fig. 7 Monitor point

3.1 阵风速度特性分析

对Ma=0.5，f=4 Hz，Amax=6°典型工况进行分析。图8 给出了该工况一个周期内垂直方向的阵风速度云图。图9 给出了该工况叶片下游阵风速度特性曲线。图9（a）为叶片下游某一监测点（x=3.55 m，y=0 m，z=0 m）处阵风速度的时间历程曲线。图9（b）为阵风速度幅值沿叶片展向分布曲线，其中虚线为叶片所在位置。从图9 中可以看出，叶片的周期性正弦摆动使得在其下游产生正弦周期性变化阵风速度场，阵风速度场频率与叶片摆动频率一致。结合图8 可以看出，叶片尾涡是叶片下游阵风速度产生的原因之一，阵风速度幅值沿叶片展长方向存在较大波动。图9（c）为不同流向位置阵风速度幅值对比，从图中可以看出，随着x 的增大，阵风速度幅值呈减小趋势。

图8 一个周期内阵风速度云图Fig. 8 Contours of computed amplitude of vertical gust velocity

图9 阵风速度特性曲线Fig. 9 Characteristics of gust field

3.2 参数影响分析

图10 给出了Ma=0.5 时，叶片下游某一监测点（x=3.55 m，y=0 m，z=0 m）处阵风速度幅值与叶片摆动频率和最大摆动幅值的关系曲线。从图10（a）中可以看出，同一叶片摆动频率（f=4 Hz）下，阵风速度幅值随叶片最大摆动幅值的增大呈先增大然后基本不变的趋势，当Amax=10°时（即叶片失速迎角附近），阵风速度幅值达到最大，10°之后阵风速度幅值基本保持不变，这可能是叶片最大摆动幅值增大时叶片失速所致。从图10（b）和（c）中可以看出，同一叶片最大摆动幅值（Amax=6°）下，叶片摆动频率改变导致叶片下游阵风速度频率改变，而阵风速度幅值随叶片摆动频率变化很微弱，呈略微减小趋势，摆动频率f 从2 Hz 增大至20 Hz，阵风速度幅值仅减小3%左右。