不同侧风类型影响下的飞机尾涡数值模拟研究

2023-03-15潘卫军罗昊天罗玉明王靖开姜沿强

科学技术与工程 2023年4期

潘卫军，罗昊天，罗玉明，王靖开，姜沿强

(中国民航飞行学院空中交通管理学院，广汉 618307)

航空器在进入前机尾流区易出现俯仰、滚转等影响安全的现象。尾流作为影响民航安全的重要因素，尽管对着陆过程有着不良影响，但相关研究表明，侧风可以迅速将尾涡吹离跑道反而能减少了飞机进近阶段所需的着陆间隔距离，提高了机场的空域容量。侧风对尾涡的扰动十分复杂，包括输送以及对涡强度和衰减率的明显干扰，由于进近阶段各种阻力和湍流的影响，以及温度变化的影响，复杂和非线性的风切变梯度侧风在大气边界层的最低高度很常见，因此对侧风场中尾流的演变与探测进行更详细的研究是必不可少的。中外学者借助计算流体力学(computational fluid dynamics, CFD)方法对飞机尾流做了较为深入的研究。Holzäpfel等[1]在维也纳机场进行了在飞机尾涡在不同风速、热分层等条件下的尾涡耗散研究评估当前欧洲尾流间隔重新分类标准(European wake turbulence categories and separation minima on approach and departure, RECAT-EU)以缩减飞机间隔；Zholtovski等[2]通过可压缩雷诺平均(Reynolds-averaged Navier-Stokes, RANS)解析了矩形 NACA0012几何机翼涡流周围的近场及其边界层问题；Landa等[3]应用剪切应力(shear stress transport, SST)涡流黏度模型和(Speziale-Sarkar-Gatski/Launder-Reece-Rodi, SSG/LRR-ω)雷诺应力模型的数值模拟结果验证湍流模型捕捉涡流发展的能力；Crouch等[4]采用二阶矩RANS CFD方法对尾涡初始耗散阶段和演化阶段进行了模拟并通过尾涡强度评估机场终端区飞机之间安全间隔；Knopp等[5]提出了一个在显著的正压力梯度和高雷诺数下的湍流边界层流动实验，使用SSTk-ω(k为湍动能，ω为耗散率)、SSG/LRR-ω和JHh-v2雷诺应力模型进行RANS模拟；文献[6-8]将SST双方程k-ω湍流模型和双方程间歇性-过渡起始动量厚度雷诺数(γ-Reθ)过渡模型在CFD求解器中的传输方程一般框架下进行整合进行了关于多级低压涡轮机中的非稳态湍流和过渡效应的模拟；Imiela等[9]采用SSTk-ω湍流模型和C型网格的组合实现了对空气动力系数的预测，并提出了一种通过CFD和计算气动声学(computational aeoro acoustics，CAA)生成360°翼型极曲线和气动声学特性的方法；Maruyama等[10]引入了使用机器学习方法对封闭系数进行统计推断的框架，以提高CFD在飞行包线限制范围内可靠湍流建模的预测能力；Sedlacek等[11]通过数值模拟对三角翼和双三角翼尾涡主导流的网格灵敏度和建模误差依赖性分析了不同湍流闭合的影响；钱宇等[12]通过建立飞机着陆状态仿真模型，采用结构化网格对模型进行网格划分，利用转捩修正的SSTk-ω湍流模型，通过求解不可压缩的N-S(Narier-Stokes)方程对生成的网格进行数值计算，得到了着陆状态下机翼及近翼流场翼尖涡的连续演化过程；谷润平等[13]为研究尾流特性，降低飞机运行风险，基于数值模拟的研究情况，采用大涡模拟的方法，借助ANSYS软件对尾流进行仿真模拟；张宇轩等[14]以NACA0012机翼为对象，采用大涡模拟方法，研究了3组不同的马赫数和雷诺数下翼尖涡中主涡和次级涡的演化特性以及其对机翼气动力的影响;张钧铎等[15]应用升力面模型和自适应网格大涡模拟方法，模拟了国产ARJ21飞机尾涡在3种侧风条件下的演化与衰减过程，并对尾涡流场结构特点进行了分析；周金鑫等[16]引入了多相流模型，利用数值模拟方法研究了不同降雨强度条件下，尾涡演化过程中的特征量演变规律；Xu等[17]通过引入地面吹吸区增强机场跑道上的飞机尾涡衰减；Li等[18]通过数值模拟研究均匀侧风对流和线性垂直切变侧风对流对尾涡强度的影响；Zhou等[19]采用基于解的动态自适应网格方法计算尾涡演化并采用基于自适应网格的大涡模拟(large eddy simulation, LES)对环境湍流强度不同的3种情况进行了研究；潘卫军等[20]使用添加旋转修正的SST-RC模型对A330-200飞机进行全机数值模拟研究侧风下的翼尖涡耗散。在飞机尾涡探测研究中，Köpp等[21]在塔布斯机场进行的现场试验，验证了多普勒相干雷达(coherent Doppler lidar,CDL)实现了尾流从产生的时刻到尾涡衰减阶段在长时间内的精确测量；Smalikho等[22]开发了一种通过1.5 μm脉冲相干多普勒激光雷达“流线”进行测量的策略，并提出了一种根据激光雷达数据估计飞机尾涡环量与涡核位置的方法；王筱晔等[23]提出了一种基于CDL频谱宽度和径向风速的尾涡快速识别方法,并基于此方法分析了典型机型的尾涡演化过程；Wu等[24]在近地效应(near-ground effect，NGE)下使用脉冲相干多普勒激光雷达(pulsed coherent Doppler lidar, PCDL)评估尾涡特性。为了实时可视化尾涡，开发了尾涡可视化演示器(vortex visualization demonstrator, V2D)。结合径向速度分布和快速傅里叶变换(fast Fourier transform, FFT)谱表征来识别尾涡的涡核位置。同时基于速度包线和Burnham-Hallock模型修正用于反演NGE下的尾涡环量[24]。因此研究侧风下的尾涡演化规律，对CDL激光雷达实地尾涡探测与数值模拟对比验证有重要指导意义。

目前中外学者已对线性分布的侧风和尾涡之间的相互作用有了一定的研究，但非线性分布或均匀的侧风究竟是如何影响涡流迹的并没有得到较好的模拟与验证，因此对于侧风下的激光雷达尾涡探测与数值模拟结果的对比验证非常重要。为此，建立了空客A330-300的尾涡模型，使用ICEM绘制计算域网格，在ANSYS Fluent中采用SSTk-ω湍流模型进行RANS计算，通过Tecplot和MATLAB对数值模拟后的流场结构进行后处理。在上述研究的基础上，还研究了均匀或非线性变化的侧风中的尾涡的运动。在成都双流机场进行了CDL尾涡探测实验，并结合数值模拟结果与实地探测结果进行对比。更加全面地验证了不同侧风影响下的尾涡演化规律。

1 控制方程

数值模拟使用ANSYS Fluent的双精度求解器进行，控制方程中的非定常项使用二阶隐式格式离散化；对流项和扩散项分别使用二阶迎风格式和二阶中心差分格式进行离散。SSTk-ω模型[25]用于求解不可压缩流动。

在雷诺平均中，瞬时(精确)Navier-Stokes方程中的解变量被分解为平均值(整体平均或时间平均)和波动分量。对于速度分量ui有

(1)

同样，对于压力和其他标量，可表示为

(2)

(3)

湍动能k和比耗散率ω可由式(4)、式(5)获得。

Gk-Yk+Sk+Gb

(4)

Gω-Yω+Sω+Gωb

(5)

式中：Gk为由于平均速度梯度产生的湍动能；Gω为ω的梯度；Гk和Гω分别为k和ω的有效扩散率；Yk和Yω分别为k和ω由于湍流的耗散；Sk和Sω为用户定义源项；Gb和Gωb为浮力项。

SSTk-ω模型包括BSLk-ω模型的所有增强，此外在湍流黏度的定义中还考虑了湍流切应力的传递。这些特点使SSTk-ω模型对于更广泛的流量类别(如逆压梯度流、翼型、跨音速冲击波)比标准k-ω模型和SSTk-ω模型更准确和可靠。

前面描述的BSL模型结合了Wilcox模型和k-ω模型的优点，但仍无法正确预测从光滑表面开始的边界层分离和数量。其主要原因是这两个模型都没有考虑到湍流剪切应力的传输。这导致了涡流黏度的过度预测，适当的传输行为可以通过对涡流黏度表述的限制器得到。

(6)

式(6)中：μt为涡流黏度；S′为应变率；常数α1= 0.31；系数α*抑制湍流黏度导致低雷诺数修正；F2可表示为

(7)

式(7)中：μ为流体黏性系数；y为到下一个表面的距离。

2 前期处理工作

2.1 尾涡模型

当飞机垂直方向受力平衡时，飞机尾涡垂直动量等于飞机所受重力，尾涡初始环量可表示为[26]

(8)

式(8)中：Γ0为飞机尾涡初始环量；MLW为飞机最大着陆重量；g为当地加速度；B为翼展；V为飞机速度。

尾涡初始流场采用Burnham-Hallock模型[27]模拟，切向速度Vθ(r)由初始环量Γ0、涡核半径rc和到涡心距离r决定，其表达式为

(9)

Gerz等[28]对初始涡核半径rc0通过初始涡核间距b0定义为

(10)

rc0=0.052b0

(11)

时间尺度t0描述了飞机或飞机模型产生下沉的尾涡对向下沉降一个初始涡核间距的时间，可表示为

(12)

式(12)中：s为翼展载荷系数，取值π/4。

无量纲时间t*可表示为

(13)

式(13)中：t为尾涡耗散时间。

(14)

2.2 网格和初始化

选择A330-300作为研究对象，其机型参数如表1所示。

计算后得到A330-300尾涡参数如表2所示。

表 1 A330-300机型性能和尺寸参数

表 2 A330-300尾涡参数

为提高尾涡数值模拟精度，迎风耗散格式采用二阶迎风，湍流模型采用SSTk-ω模型。选取固定位置涡核周围处涡量，进行的网格无关性验证如图1所示，全局网格尺寸0.7 m较1 m网格尺寸特征位置处涡量缩减13.8%；在0.5 m网格尺寸时涡量值较0.7 m网格尺寸时特征位置处涡量缩减3.84%，并在0.3 m网格尺寸时趋于平缓，较于0.5 m的网格尺寸涡量减少1.69%，远小于5%，符合网格无关性验证条件，在此基础上继续增加网格数量对结果影响较小。

图1 网格无关性收敛曲线

同时采用基于有限体积法的结构化网格完成计算域网格划分，相比于非结构化网格，结构化网格更为精细，可以更好地提高参数收敛性与计算精度，设置网格尺寸为0.3 m。以涡核间距中心点为坐标原点;以垂直于地面方向左侧为y轴正方向；以速度入口来流方向为x轴正方向。在Fluent UDF中解释H-B模型，计算域的初始涡核位置和生成的初始尾涡速度场分别如图2和图3所示。

（一）疫病防控将保险所用资金用于疫病防控，除疫苗注射和圈舍定期消毒投资以外，当出现重大动物疫病需要对病死牛羊进行无害化处理时，在条件允许的情况下可以进行补贴。

图2 计算域和初始涡位置

图3 初始切向速度剖面图

2.3 边界条件和计算参数

风切变效应的非线性函数可采用对数律分布或指数律分布来进行描述，但指数率计算的风速值与实测值偏差较小，且用指数律分布计算风速轮廓线比较简便。风切变的指数律分布可表示为[29]

(15)

根据民航局CAAC-25文件，对陆上飞机应制定在干跑道上对起飞和着陆演示是安全的90°侧风分量，该分量必须至少为20节或0.2倍基准失速速度(0.2VSRO)，取大者，但不必超过25节。因此最大侧风风速取25节。

侧风切变的垂直梯度可以优先降低涡旋的下降速度，导致它们倾斜，增加分离，有时甚至向上上升，环境侧风高度的垂直二阶导数S定义为[31-37]

(16)

(17)

式(17)中：S*为无量纲化的S；v0为尾涡初始下降率。

图4为尾涡对侧风切变效应示意图。

绿色曲线为侧风切变梯度为负或正时的侧风垂直剖面；左涡由红色圆圈表示；右涡由蓝色圆圈表示；θ为尾涡对的倾角

为研究尾涡在不同侧风条件下的演化规律，分别取1、3、7、10 m/s的4种均匀侧风、线性侧风和非线性侧风，速度随高度变化的具体函数关系如表3所示。流体域环境变量如表4所示。

表3 不同类型侧风的函数关系

表4 流体域环境变量

3 仿真结果

3.1 数值模拟结果

尾涡运动示意图如图5所示。

图5 侧风下的尾涡运动示意图

二维Q准则(Qcriterion)可表示为[38]

(18)

式(18)中：u、v分别为x、y的速度方向。

图6 均匀侧风1 m/s

图7 均匀侧风3 m/s

图8 均匀侧风7 m/s

图9 均匀侧风10 m/s

图10 线性侧风(S*=0)

图11 非线性侧风(S*<0)

图12 非线性侧风(S*>0)

对以下变量进行无量纲化处理。

(19)

(20)

(21)

(22)

由图13可以看出，左右涡横向偏移量受侧风风场大小影响，侧风越大偏移量越大。在S*>0或S*<0的垂直非线性侧风条件下，横向偏移量更大。

图13 不同侧风类型下的尾涡无量纲横向偏移随无量纲时间变化示意图

由图14可知，在任何环境侧风条件下，左涡下沉量较右涡大，左右涡下沉趋势不变。对比分析左右涡纵向运动趋势发现尾涡的高度变化受不同侧风影响较大，侧风风速的改变将对尾涡的纵向偏移产生截然不同的运动轨迹和趋势。这是因为左涡直接和侧风接触，相互作用并融合，导致左涡的空间高度较右涡在垂直方向上有着更强的波动。在施加的持续侧风风场扰动下尾涡的速度场将在短时间内发生剧烈改变，尾涡对在空间上的分布也趋于复杂。

图14 不同侧风类型下的尾涡的无量纲下降高度随无量纲时间变化示意图

由图15可以看到，侧风风场对尾涡存在明显的输运作用，且随着风速的增加输运作用更为强烈。在均匀侧风和线性侧风条件下尾涡涡流迹在等距离X下沉更多，而非线性侧风条件下则在等距离Y下横向侧移更多。

图15 不同侧风类型下的涡流迹示意图

3.2 实地探测数据对比

实地雷达探测在成都双流国际机场获得的A330-300着陆时的各参数反演值与假设的A330-300尾流的数值模拟的结果接近。由于尾涡迅速减弱，可能只在有限的时间内被探测到。

图16为根据H-B模型计算出尾涡初始的切向速度分布与实地探测结果对比。可以看出，由于左涡(顺时针)、右涡(逆时针)旋转，使得尾涡的切向速度带有方向性，负值代表尾涡切向速度方向为y轴负方向，对比雷达尾涡切向速度变化曲线可以看出，随径向距离增加切向速度呈现出先增后减的变化趋势。尾涡模型在切向速度分布上的差异与探测结果基本一致。

图16 H-B尾涡模型切向速度分布与CDL雷达实测分布

现场测量通过雷达回波判断尾涡是否存在并根据回波数据给出尾涡位置、强度等信息。基于CDL观测到的非线性侧风条件下的A330-300型客机尾涡演化的径向风速图如图17、图18所示。可以看出，对应时刻的尾涡强度较大，正负风速绝对值显著大于周围环境风场。据此可以判断，尾涡已经产生，经过数据处理后，可以判断为尾涡初始位置。

图17(e)、图18(e)分别为为同时段尾涡演化的频谱宽度组图，其呈现的尾涡基本演化特征与径向风速图和数值模拟径向速度云图[图17(c)、图18(c)]一致。因初始涡核间距较大，如图17(d)、图18(d)所示，可明显区分左右尾涡及涡核位置，频谱宽度较大。

图17 径向风速(2018年9月9日 08：37,机型：A330-300,S*<0)

图18 径向风速(2018年9月19日 15：19,机型：A330-300,S*>0)

从图17(c)可以看出，在受左侧非线性侧风(S*<0)条件下左右涡核的位置距离地面约37 m和42 m，距离CDL分别约517 m和529 m，涡核间距约为11 m。随后左右尾涡均从高空下降并向外扩散，其中上风涡低于下风涡，符合数值模拟结果。

从图18(c)可以看出，在受右侧非线性侧风(S*>0)条件下左右涡核的初始位置距离地面分别约20 m和21 m，距离CDL分别约254 m和288 m，涡核间距约为34 m。随后左右尾涡均从高空下降并向外扩散，其中上风涡高于下风涡，与数值模拟结果左侧非线性侧风结果一致。

从图19数值模拟结果中可以看出，尾涡对耗散开始之前，两种不同非线性环境侧风下的环量都会增加。左涡和右涡的环量在开始时是一致的。在S*<0的情况下，左涡(上风)环量在t*=0.1不断增加，在t*=1.1之前，右涡(下风)环量大于左涡环量，随后发生改变；CDL激光雷达探测反演结果显示右涡环量在t*=0.8时刻之前右涡环量一直大于左涡环量，与数值模拟结果基本一致。相反，在S*>0的情况下，左涡(上风)环量大于右涡(下风)环量。

图19 CDL探测和数值模拟结果的尾涡无量纲环量与无量纲时间的关系

通过数值模拟可以看到由于诱导速度与尾涡强度成正比，上面讨论的环量差异导致两个尾涡以不同的速率下降。尾涡对的倾角θ定义如图4所示。从图20可以明显看出，当S*为负时倾斜率逐渐增加，在t*=1.0时刻后几乎呈线性增长；同时，当S*为正时，倾斜率也逐渐增加；这种情况下，尾涡对倾斜率在t*=0.8内增长幅度相对较小，然后迅速增加并最终也呈线性增长趋势。尽管在两个非线性侧风模拟中|S*|被设置为相同，但倾斜率的变化却是不同的。这可能有两个原因；首先，有其他环境因素会影响下降过程。其次，侧风不会对称地改变尾涡对的倾斜。