多机型组合下的尾流遭遇仿真研究

2023-12-08潘卫军张钰沁姜沿强王靖开罗昊天

科学技术与工程 2023年31期

潘卫军, 张钰沁, 姜沿强, 王靖开, 罗昊天

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

飞机在飞行过程中会形成扰乱原本空气流场的一对反向涡旋。随着前机尾涡的不断演化与下沉,机翼后缘会产生较大的诱导力,后机若在此刻进入前机的尾涡流场,机翼会被不断施加诱导力,导致升力急剧变化,从而使机体发生过度的滚转、俯仰,甚至失控。

国外率先对尾流遭遇风险进行研究,Speijker等[1]提出使用诱导滚转角速度的大小来计算后机遭遇尾涡的影响。Campos等[2]将最大滚转角速度作为评价尾涡遭遇严重程度的指标,提出了一种关于前机尾流对后机滚转稳定性影响的理论。Visscher等[3]提出了一种用于地面附近尾涡行为预测的快速时间模型,使用滚转力矩系数对尾流遭遇的安全性进行评估。Compos等[4]使用滚转力矩计算模型计算后机所受影响,并使用最大滚转角速度来评价尾流遭遇的严重程度。Holzäpfel等[5]使用两种不同的大涡模拟(large eddy simulation, LES)代码对各种湍流和稳定分层大气环境中的飞机尾涡演化进行了大涡模拟。中国对于尾流遭遇的研究有:谷润平等[6]提出了一种新型的尾流遭遇危险区的评估模型,同样使用滚转角速度来评价尾涡遭遇风险。潘卫军等[7-8]通过在最后进近航迹与侧风影响下后机对尾流的承受能力,建立了配对进近条件下尾涡间隔优化模型。张钧铎等[9]使用数值模拟的方法,采用自适应网格大涡模拟技术数值研究ARJ21客机尾涡在侧风条件下的近地演化过程,并分析在不同侧风条件下尾涡的演化与衰减特性。林孟达等[10]采用自适应网格技术,大幅减少了所需网格量,提高了飞机尾涡演变模拟中的计算效率。潘卫军等[11]通过模拟全机尾涡场结构,对侧风影响下的飞机尾涡四涡系结构进行了进一步的探究。潘卫军等[12]通过A330-300机翼后缘网格加密的方法,提高了网格质量,采用大涡模拟研究总结出了进近阶段不同侧风条件对尾涡耗散的影响规律。潘卫军等[13]对ARJ21飞机近地阶段的安全性进行研究;潘卫军等[14]应用ANSYS FLUENT UDF(user defined function)编译环境侧风不同的七种情况进行尾涡耗散机理的数值模拟。

目前尾涡的演化多以尾涡耗散模型进行计算,将尾涡的耗散阶段分为扩散阶段和快速衰减阶段,使用分段函数进行表示。相比于这种传统的尾涡耗散模型,现通过使用计算流体力学的方法进行尾涡数值模拟,从而得到全流场内的时空尾流信息。并将其与后机空气动力学响应模型相结合,进行全面的尾涡演化机制与进行多机型下的尾流遭遇风险研究。

1 尾涡演化仿真

1.1 控制方程

尾流数值模拟采用雷诺平均法RANS(Reynolds average Navier-Stockes),雷诺平均数值模拟将雷诺平均方程作为控制方程,进行数值模拟。其连续性方程与动量方程分别为

(1)

(2)

式中:ρ为流体密度;t为时间;ui表示计算域中xi方向的速度;uj表示计算域中xj方向的速度;p为流体压力;μ为流体黏性系数;δij为应力张量分量;τij为亚格子Reynolds应力。

湍流模型选取SSTk-ω模型,能够准确及时预测分离的特性,适用于旋转流动的情况。湍流动能k和比耗散率ω可从以下SSTk-ω方程中获得,即

Gk-Yk+Sk+Gb

(3)

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

(4)

式中:Gk表示湍流动能;Gω代表ω的产生;Γk和Γω分别代表k和ω的有效扩散率;Yk和Yω代表k和ω由于湍流的耗散;Sk与Sω为用户自定义项,分别为湍动能项与湍流耗散原项;Gb和Gωb为浮力项的解释。

1.2 前机尾涡速度模型

采用Hallock-Burnhan模型描述前机产生的尾流的切向速度[15],其公式为

(5)

式(5):r为空间任意一质点到单涡涡心的半径;r0为初始涡核半径;V(r)为尾涡在该点处的切向速度;Γ0为尾涡初始环量。

尾涡初始环量Γ0的计算公式为

(6)

(7)

式中:m为飞机质量;g为重力加速度;ρ为空气密度;b为飞机翼展;v为前机飞行速度;b0为初始涡核间距。

1.3 网格初始化

使用ICEM CFD(the integrated computer engineering and manufacturing code for computational fluid dynamics)构建尺寸大小为(-150～150)m×(+50～-350)m的结构化网格,网格的(0,0)点为飞机重心的坐标点,向上为正方向,向下为负方向。

1.4 参数设置

将尾涡流场设置为不可压理想气体,一般大气湍流强度,计算方法选择瞬态,求解方法选取稳定性较好的couple算法。使用用户自定义函数编译前机尾涡场模型参数,设置机型数据参数进行计算。选取中国民航使用率较高的飞机B777和A330作为前机尾涡数值模拟的对象,其ICAO RECAT分类皆属于B类机,以便更好地进行同机型之间的比较[16]。相关参数如表1所示。

表1 前机机型参数表Table 1 Parameter table of previous airccrafts

将前机机型参数、环境参数等导入Fluent,模拟出各机型以时间为计量单位的尾涡演化及消散过程,图1所示为A330-200尾涡演化至30 s时的切向速度场V(r)分布图。

2 后机遭遇响应

后机遭遇前机尾流时,会受到不同程度的上洗或下洗的作用力使其机翼受力发生变化,相应地,飞机的高度、俯仰、滚转等姿态也会改变。后机纵向跟随前机飞行时,主要分为纵向进入前机单涡区域与纵向进入前机双涡区域两种方式。如图2 A所示,进入前机单涡时,左右机翼受到的上下洗气流方向相反,这将使得后机主要发生滚转变化。如图2 B所示,进入前机双涡时,后机主要受到下洗气流,这将引起后机明显的掉高度与俯仰变化。下文将主要对后机跟随前机飞行时的遭遇尾流受力模型,进行具体的受力变化、滚转变化以及危险区分析。而当后机横穿前机尾涡时,主要发生较大幅度的俯仰变化,引发机体大幅度颠簸,甚至使机体失控。

图2 后机跟随前机飞行时的尾流遭遇情况Fig.2 Wake encounter when the rear aircraft follows the leading aircraft

2.1 飞机的升力变化量

飞机的升力变化量计算公式为

(8)

(9)

(10)

2.2 机翼诱导滚转力矩

后机机翼诱导滚转力矩的计算方式基于传统力矩计算方法,采用力与距离的乘积来表示单一受力的力矩,其计算公式为

Z=DF

(11)

式中:Z表示滚转力矩;D表示力臂;F表示尾涡的诱导滚转力;Δβ表示迎角变化量。

采用条带法对后机机翼所受到的滚转力矩进近计算,条带状法是一个一个分布式的空气动力学模型,适用于尾流遭遇的不均匀速度场气动计算。同时它实现了最大迎角的特殊限制,以防止局部迎角超过最大的迎角限制,其计算公式为

(12)

(13)

(14)

式中:dΓ(y)表示局部滚转力;ΔZ表示局部滚转力矩;Δβ表示迎角变化量。

通过积分可得尾涡对飞机产生的诱导力矩:

(15)

2.3 滚转力矩系数

基于其对于飞机遭遇的良好判别性和无量纲性,滚转力矩系数可用于划分后机遭遇危险区的安全指标,其计算公式为

(16)

式(16)中:RMC为滚转力矩系数;V为后机的飞行速度;S为后机机翼面积;B为后机翼展。

2.4 后机机型组合

为了研究的实用性考量,将选取中国数量较多、投入使用时间较长的机型B737、A320以及国内自主设计制造的支线客机ARJ21作为响应后机。由于机型布局与参数的不同,各机型所能承受的最大滚转力矩也不同。表2所示为各类机型具体的参数。

表2 后机机型参数表Table 2 Parameter table of following aircrafts

3 结果分析

将数值模拟实验的结果与后机响应相结合,计算机型遭遇尾流时的可接受安全水平。下面均以双涡涡核中心位置为坐标原点对尾涡演化及危险区进行分析。

3.1 尾涡耗散

使用环量Γ表示尾涡的强度,反映尾涡消散速率。提取以涡心为圆心,半径为5～15 m区域的11个圆面环量面积分的平均值进行环量计算,再使用初始环量进行量纲统一的方法表示环量变化,依次提取0～100 s的涡心下沉位置的变化图与环量数据来反应尾涡的耗散,Γ′5～15表示半径为5～15m区域的11个圆面环量面积分的平均值,如图3所示。

采用数值模拟的方式更有利于分析尾涡耗散过程中强度的变化。从图3可以看出,通过数值模拟得出的尾涡的消散速率符合尾涡耗散的两个阶段,即0～20 s为尾涡的初始耗散阶段,尾涡依靠自身衰减,耗散速度较慢,而20 s后为尾涡的快速耗散阶段,耗散速度明显加快。虽然属于同类机型,但B777的初始环量明显大于A330。

图4为A330与B777两种机型所产生的尾涡涡心下沉速率示意图,可以看出B777的尾涡下沉速率明显高于A330。虽然A330与B777均为B类机型,但其尾涡涡心下沉的速率有明显的差异,故对于同类机型同样需要分别研究其尾流演化。

图3 环量变化图Fig.3 Circulation change chart

图4 涡心下沉速率示意图Fig.4 Schematic diagram of vortex sinking rate

3.2 升力变化量

根据现行的国际民航组织的尾流间隔标准所规定的跟随各类机型的升力变化量最大值,可以确立不同机型组合下的模拟安全距离。故分别提取不同机型组合下尾涡演化至0、30、60、90 s时所对应的涡核附近处的升力变化量,结果如图5所示。图5是以A330为前机,不同后机在涡核水平线上的升力变化图,两涡涡心初始生成点的横坐标为±25.45 m,而升力变化量所处位置在最大值所在位置在±40 m左右。由于没有添加侧风的影响,升力变化量曲线的趋势呈现明显的对称性。以右涡为例,在涡核外侧,升力不断增大,变化量先增大至最大值,后逐渐减小并趋向于0。在双涡内侧,升力变化量急剧减小。

对比0、30、60、90 s的曲线可知,随着尾涡的演化,升力变化量的最大值取值逐渐减小。特别地,0 s的升力变化量曲线的最小值在并非在两涡中心处而是两涡中心两侧取得,而30、60、90 s的曲线升力变化量在两涡中心处取到极值,且随着不断演化变化量的大小逐渐减小。这种差异是由于尾涡演化初期,涡核半径较小,左右涡影响程度有限导致两涡中心处的下洗气流弱于其两侧。而随着尾涡的演化,涡量逐渐扩散,两涡中心的下洗气流会逐渐强于其两侧。

图5 不同后机跟随A330的升力变化图Fig.5 Lift change diagram of different following aircraft with A330

对于相同前机,不同后机所受的升力变化量大小不同,对于D类机型,如A320,B737等,其最大升力变化量在1.5×105～2×105,而F类机型,如ARJ21,其最大升力变化量在1×105～1.5×105,同样的,处于两涡之间的区域,选取的D类机型的升力变化量的大小也大于E类机型。

图5(a)与图5(d)展示了对于不同前机,同类后机的机型组合的升力变化量曲线。B777-A320组合的后机升力变化量最大值大于A330-A320的机型组合。由图3可知,B777的初始环量大于A330,且随着尾涡的演化,B777的环量始终大于A330。由此可知,前机尾涡强度的越大,后机所受的升力变化量的最大值也会越大。

3.3 危险区划分

滚转力矩系数通常用来衡量飞机尾涡遭遇的危险程度,根据Lang等[17]的实验结论,0.05～0.07为仅使用副翼控制飞机的滚转的最大值,若超过此安全阈值,飞机将失控。根据尾涡附近速度场分布情况,可计算得出该区域内的滚转力矩系数分布。图6所示机型A330-B737组合下尾涡耗散初期流场的滚转力矩系数分布图。滚转力矩系数为正表示后机滚转方向为顺时针方向,滚转力矩系数取值为负表示后机滚转方向为逆时针方向。

将图6滚转力矩系数进行绝对值处理,根据B737的滚转力矩系数最大取值范围,选择绝对值为0.047的滚转力矩系数的线性范围作为机型B737危险区的边界,得到如图7所示的A330-B737危险区示意图。其中外部深蓝色及以内为滚转力矩系数大于0.047的空间,表示当该飞机在其内运行时,滚转力矩系数将超过0.047,飞机的姿态会受到一定的影响难以改出,此时后机会处于不安全的运行状态。内部黄色区域表示滚转力矩系数大于0.07的空间区域,表示当该飞机在其内运行时,滚转力矩系数将超过阈值,后机一旦进入该区域,将会面临失控的巨大危险,基本无法改出。随着尾流的下沉运动,危险区的范围在空间上也会呈现出下沉的趋势。

图6 A330-B737滚转力矩系数分布图Fig.6 Distribution of A330-B737 rolling-moment coefficient

图7 A330-B737危险区模拟图Fig.7 A330-B737 hazard zone simulation diagram

由于不同的机型具有不同的参数,滚转力矩系数的限制也会不同;不同前后机组合下的危险区域也会有所差异。图8为机型组合A330-ARJ21的尾涡危险区三维可视切片图,提取间隔为10 s,直至危险区消失。其中,以初始双涡为起点,涡中心点为坐标原点,H表示涡的下沉高度,x表示双涡距离中心的距离,t表示尾涡演化至该位置所需的时间,可以看出在尾涡演化初期,尾涡强度高,尾流危险区呈现出两个主危险区与四个副危险区,四个副危险区分别位于两个主危险区两侧。两涡危险区之间的区域虽然满足滚转力矩系数阈值的限制,但其区域空间较小且不可容纳航空器在其间飞行,所以仍需将其视作危险区进行划分。根据模拟仿真将得到多种机型组合的危险区进行比较。

图8 A330-ARJ21尾涡危险区三维可视图Fig.8 A330-ARJ21 tail vortex hazard zone three-dimensional diagram

根据不同机型组合危险区的持续时间,可计算出不同机型前后机尾流纵向安全间隔,结果如表3所示。

通过计算结果得出,通过数值模拟尾涡危险区的纵向距离符合RECAT-CN和RECAT-EU的尾流安全间隔标准。选择B型机为前机,后机选择D型机时,其尾流安全间隔与RECAT-CN相近。后机ARJ21时,模拟仿真的尾流遭遇危险区安全间隔远远小于RECAT-CN的标准,这与其可承受的滚转力矩系数较大、下单翼结构以及气动分布特性有关。

比较图9、图10所示A330与B777危险区俯视图,可以观察到尾涡耗散至0～20 s内,危险区范围较大,同时包含危险主区与副区,危险区的宽度为60～100 m,随着尾涡强度的逐渐减小,危险副区逐渐消失。此后危险区的宽度保持为60 m左右,小幅度减小。

表3 尾流纵向安全间隔表Table 3 Wake longitudinal safety seperation table

图9 A330危险区分布图Fig.9 A330 hazard zone distribution

图10 B777危险区分布图Fig.10 B777 Hazard zone distribution

对于同类型后机,比较B777-ARJ21与A330-ARJ21的前后机组合,在0～20 s内,两种机型组合的危险区宽度减小速率不同,前者在尾涡演化20 s左右时危险副区消失,而后者危险副区在10 s左右消失。其后宽度减少至60 m左右,但A330-ARJ21的机型组合的危险区宽度始终大于B777-ARJ21。最终的危险区消失时间A330-ARJ21也略早于B777-ARJ21。

对于不同后机而言,ARJ21的危险区在80 s左右消失,明显早于A320与B737。A320与B737的危险区在100 s左右消失,但由于这两种后机气动特性的差距不大,所对应的危险区结束时间也几乎一致。