不同B-V频率下的飞机尾涡数值模拟研究

潘卫军,姜沿强,张钰沁

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

0 引言

飞机在飞行中会形成与大气环境和自身机型特征相关的、强度逐渐减弱的尾涡流场,这会对后机的飞行安全产生威胁,制约机场的起降效率。因此,研究不同气象条件下的尾流的演化规律对保障飞机安全运行以及提升机场容量具有重要意义[1]。

尾涡的发展可以分成近区阶段、扩展近区阶段、中、远区阶段和衰减区域5个阶段[2],近区阶段主要是指在上下翼面压强差的作用下机翼、尾翼处涡旋的卷起形成。随后在翼尖涡的主导下,不断发展形成1对稳定的反向涡对,并在气象因素的作用下不断衰减[3]。目前尾涡演化的研究主要依靠激光雷达的探测、数值模拟和流体实验等手段进行[4]。大量的数值模拟研究填补了尾涡发展规律的诸多空白[5-6],并促进了尾涡快速预测模型的搭建与完善[7]。尾涡演化的模拟是分析尾流遭遇安全和尾流间隔的重要基础。Gerz等[8]在对商用飞机尾涡的研究中,给出了尾流遭遇模型和尾流安全评价模型。Ruith等[9]通过直接数值模拟研究证明了尾涡耗散过程是一个自身流场稳定性决定的轴对称模型。Ghalandari等[10]通过对典型机翼的研究发现对于整机尺度而言,机翼颤振行为引起的运动较小。环境因素是尾涡演化的重要影响因素,诸多学者对不同气象条件下的尾涡演化进行了研究[11-12]。Misaka等[13]通过大涡模拟研究了飞机尾涡从卷起到衰减的演变过程,得到了在不同环境湍流的情况下主翼脱落的涡面详细卷起的过程以及涡流场的特征参数之间的相关性。周金鑫等[14]基于欧拉多相流模型,研究得出降雨会加速尾涡强度衰减,降低尾流遭遇的危险性。潘卫军等[15-16]通过对机翼后缘网格加密提高网格质量,使用添加旋转修正的SST-RC模型模拟全机尾涡场结构等方法对侧风影响下的飞机尾涡演变规律进行了深入的探究。在现代科学中,浮力(buoyancy or brunt-väisälä,B-V)频率多应用于层状自然介质的动力学分析中,其在温跃层研究之中具有良好的表现[17],也称浮力振荡频率,是流体对垂直位移(例如对流引起的位移)的稳定性的量度。近地大气情况复杂,分层现象也较为多样。对于近地飞行而言,大气分层、风速风向以及湍流耗散率的强弱都对尾涡演变有影响[18]。Greene[19]最先将B-V频率、湍流和雷诺数组合在同一个尾涡耗散预测模型中,并且预测模型与飞行测试的结果良好吻合。

部分尾涡耗散预测模型考虑了B-V频率对飞机尾流耗散的影响。但由于大气分层情况复杂,难以提出一个适用于所有耗散环境的模型。本文中通过数值模拟对其进行研究,使用UDF(user-defined functions)将B-V频率引入飞机尾流的环境场中。以空客A333为研究对象,使用ANSYS选择SSTk-ω湍流模型进行飞机尾流数值模拟。由于温度分层的存在,模拟也考虑了初始温度对尾涡耗的影响散。以得到不同初始温度和B-V频率组合下的运动与衰减特征,探究大气分层对尾涡演化的影响规律。

1 数值模拟方法

1.1 控制方程

本文中研究集中在尾涡耗散阶段,为节约计算资源采用雷诺平均(reynolds-averaged navier-stokes,RANS)方法进行瞬态计算。RANS方程如式(1)和式(2)所示。

(1)

(2)

选取shear-stress transport(SST)k-ω湍流模型[20]进行方程封闭,相比于baseline(BSL)k-ω该模型考虑了湍流剪切应力的传递,有着更高的计算精度。适用于旋转流动等情况。使模型封闭RANS方程。湍流动能k和比耗散率ω的求解如下:

Gk-Yk+Sk+Gb

(3)

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

(4)

式(3)和式(4)中:Gk为由于平均速度梯度而产生的湍流动能;Gω为ω的产生;σk和σω分别为k和ω的湍流prandtl数;ut为湍流粘度;Yk和Yw分别为k和ω由于湍流引起的耗散;Dw为交叉扩散项;Sk与Sω分别为湍动能项与湍流耗散原项;Gb和Gωb为浮力项的解释。

1.2 尾涡模型

飞机尾流是指飞机在飞行过程中,由于飞机自身升力的存在,机翼上下表面的压力差通过对机翼后方气流的扰动形成的1对反向旋转的漩涡。飞机尾涡参数一般由初始尾涡环量Γ0,初始涡核半径r0以及初始涡间距b0,3个基本参数描述。由kutta-joukowski圆柱绕流定理将环量与升力结合起来,如式(5),得到Γ0的计算公式:

(5)

式(5)中:m为飞机质量;g为重力加速度ρ为空气密度;V为飞行速度。Gerz等[8]研究发现尾涡的初始间距与飞机翼展b长度之间的比例关系与机翼载荷系数有关,并且该比例接近于π/4。初始涡核半径与初始涡核间距的比值在0.01～0.05之间[9],本文中取0.035,如式(6)所示。

b0=πb/4,r0=0.035b0

(6)

Holzäpfezl研究发现使用hallock-burnham(HB)涡模型进行尾涡耗散快速预测的结果更加符合雷达探测的实际值[7]。如式(7),借助HB模型可以良好表示单个尾涡环量和切向速度的分布,将左右双涡流场叠加可得飞机尾涡双涡分布模型:

(7)

式(7)中:vθ表示尾涡的切向速度;Γl、Γr为左涡和右涡的初始环量,其大小等于Γ0;rl、rr为流场中一点到左涡和右涡心的距离。根据Γ0与b0可以进一步定义尾涡演化的特征速度v0与特征时间t0。

v0=Γ0/2πb0,t0=b0/v0

(8)

式(8)中,v0是根据biot-savart定律所得尾涡在相互诱导下的初始下降速度;t0为尾涡以v0下降距离b0所需时间。

本文中以A333作为尾涡演化的对象,其翼展为60.3 m,航空器质量m为最大着陆质量(MLW) 187 000 kg,飞行速度70 m/s。尾涡初始位置处的初始温度分别设置为263、283、303 K,其对应的相关尾涡参数如表1。

表1 不同初始温度下的尾涡参数

使用飞机翼展,初始环量和特征时间对物理量进行归一化处理,得到无量纲化时间t′=t/t0,无量纲化各向位移x′=x/b,h′=h/b,无量纲化环量Γ′=Γ/Γ0。x为翼展方向上一点距离原点的距离,h为距初始涡核位置的高度差。

(9)

1.3 环境场

当尾流在大气层中下降时,会受到粘性力与浮力的影响。表示流体分层的参数是B-V频率N,表示流体粒子从平衡位置位移的假设振荡频率,即垂直浮力振荡固有频率[19]。对于理想气体,B-V频率可以表示为

(10)

式(10)中:Ti表示大气温度;h为飞行高度;A为尾流椭积;θ为位温,表示干气块干绝热压缩或膨胀至1 000 hpa所具有的温度;Cp为定压比热。

N*=Nt0,283

(11)

如式(11),使用T=283 K时的特征时间t0对B-V频率无量纲化,选取无量纲B-V频率N*为0,0.5和1.0的3种情况进行研究。

2 尾涡的演化与验证

2.1 尾涡演化设置

构建分辨率为0.5 m的结构化网格,以初始尾涡左右2涡涡核连线中点为网格的原点。通过UDF将HB尾涡模型,在ANSYS Fluent中编译添加到搭建的环境场之中,得到A333的初始切向速度分布,如图1所示。

图1 初始切向速度分布

求解器选择基于压力基(分离式)的求解器。算法选择PISO(pressure-implicit with splitting of operators),该算法将SIMPLE(semi-implicit method for pressure linked equations)所需的重复计算转化为求解压力修正方程,使修正后的速度满足连续性方程和动量方程,从而大大减少收敛所需迭代次数,适用于瞬态计算[21]。梯度离散方案选择基于单元体的最小二乘法插值。为保障计算精度,选择二阶隐式瞬态公式,压力项、能量项以及其他项的离散格式均选择二阶迎风格式。

飞机尾涡场流体速度远小于1.0马赫数,压缩性效应可以忽略不计,所以气体选择为不可压缩流体。无量纲大气湍流耗散率ε′=(εb0)1/3/v0设置为0.23(ε为大气湍流耗散率)[18]。

2.2 数值模拟与雷达探测对比

提取激光雷达探测到的尾涡信息,与T=283 K,N*=0时的数值模拟计算结果相比较,进一步验证数值模拟方法的可靠性。

图2为雷达探测得到的单涡涡核附近的速度分布曲线和数值模拟得到的双涡速度分布曲线。对于单涡而言,雷达探测结果和数值模拟结果都体现出涡核附近的速度最大,远离涡核的速度逐渐减小的特征。对于无侧风影响的数值模拟,左右涡的变化基本对称,所以在后续结果分析中环量变化只给出了单涡的结果。受雷达扫描周期和角度的限制,对于同一尾涡,通常只能得到3～4个数据。由于雷达探测受到环境场的影响,所以尾涡径向速度分布与模拟的数值结果不完全相同。

图2 尾涡速度分布

3 计算结果分析

3.1 涡结构与涡强度的变化

图3为不同温度场下保存步长为20 s的涡量切面图。由图3(a)(b)(c)可得,不同温度条件下,相同B-V频率下的大气层结稳定性一致,尾涡结构变化趋势基本相同。从图3(b)(d)(e)可得,同一温度,随着B-V频率的增大尾涡下沉速率明显减小。在尾涡演化至t′=1.5附近时,尾涡涡形发生较大变化,N*=1时,尾涡主涡上方出现了反向的小型次级涡,使尾涡涡心间距出现了小幅的减小。这是由于不稳定的大气层结,可以加速尾涡的耗散,促进涡形的改变。同时B-V频率的增大会加强尾涡在下沉时所受的惯性振荡,抑制尾涡下沉,从而加快尾涡耗散,减小尾涡垂向风险。

图4中无量纲量均以T=283 K时对应的初始量进行归一化处理。如图4所示,演化后期各温度下的环量值趋近一致。不同的大气分层对初期环量下降的作用较小,差值在0.5%以内。但由于不同温度下的初始环量不同,难以判断环量初始下降速率的不同的变化是由温度还是由初始尾涡强度引起的。为控制变量,选择尾涡初始环量为443.41 m2/s进行各初始温度下的补充计算。

图3 尾涡结构变化

图4 不同温度场下的环量变化

由图5可得,各初始温度下的尾涡环量变化趋势差别不大,初始温度对尾涡强度减弱影响较小。尾涡耗散前期,具有一定的结构稳定性,大气分层的影响较小,环量的下降主要受两涡互诱导效应[18]。随着涡量的扩散以及尾涡强度的下降,大气分层对耗散的影响逐渐加强,B-V频率的增大会抑制尾涡强度的减小。t′=3时,初始温度283 K,无量纲浮力频率N*=0、0.5、1条件下的尾涡环量分别减小了71.97%、73.87%、75.58%。

图5 相同初始尾涡强度下的环量变化

3.2 尾涡位移的变化

图6反映了涡核运动情况,在重力、浮力和空气黏性力的作用下尾涡呈现出不同程度的下沉。由图6(a)可得,B-V频率对尾涡运动的影响较为明显。图6(b)表示尾涡横向位移随时间的变化,各演化场景下的横向位移的范围为-0.24～0.26个翼展。可以看出,B-V频率对尾涡横向运动影响较大,而初始温度对尾涡横向移动的作用较小。N*=0时两涡向内侧运动,N*=0.5和N*=1两涡向外侧运动。

图6(c)则体现出初始温度与B-V频率共同作用于尾涡的下沉运动。初始温度越高所对应的尾涡下沉速度越大,但影响程度较小在6.2%以内。受B-V频率的影响涡核下沉速率明显减小,t′=3时N*=0所对应的下沉距离为2.7个翼展,为N*=1时下沉距离的2.71倍。

3.3 温度场的变化

由于N*=1,N*=0 2种情况下温度分层存在明显差异,本文中选取初始温度为288 K,不同B-V频率进行温度变化的研究。提取两涡涡心连线中点左右各6.6倍翼展的温度变化,结果如图7所示。涡旋形成初期会扰乱原本的温度分层,尾涡内部温度低于环境温度。随着尾涡自身的向内向下的涡旋以及尾涡的下沉,会导致N*=1条件下涡核间的温度显著大于左右两涡,整体尾涡温度大于环境温度;N*=0条件下涡核间的温度显著小于左右两涡,整体尾涡温度小于环境温度。由PROCTOR的三阶段耗散模型[22]可知,随着尾涡长时间的发展,与外界不断地热交换,尾涡内部温度会逐渐趋近环境温度,直至尾涡消散。

结合图6和图7可得,尾涡的横向位移由于无重力影响,主要受到大气分层中浮力效应影响。N*=1时,外部温度小于尾涡内部,涡核间距不断减小;N*=0时,外部温度小于尾涡内部,两涡向外侧移动,间距不断增大。

图6 尾涡位移变化

图7 尾涡横向温度场变化

3.4 安全性分析

近地面气象条件复杂,B-V频率的值较为多变,这是建立尾流间隔预测系统时需要重点考虑的。由前面的环量分析可得,虽然较大的B-V频率会加快尾涡强度的减小,提高了后机遭遇前机尾涡时的安全性,但是减弱的程度十分有限。图8为受B-V频率的影响尾涡下沉速率不同所引起的不同前后机尾流安全间隔的示意图。淡蓝色与淡绿色的区域表示N*=0和N*=1大气条件下,随尾涡下沉前机尾涡危险区向下位移的情况。可得,N*=1条件下小于N*=0条件下的危险区下沉速率,相应的尾流间隔更大。这是由于较大的B-V频率会导致尾涡下沉大幅度的减慢,促使前机尾流在后机的飞行航路上停留更久,造成后机需要与前机保持更大的安全距离以避免尾流遭遇对后机飞行安全的影响。

图8 不同B-V频率下的尾流间隔

4 结论

1) 不同初始温度主要会引起尾涡初始强度的不同,但对尾涡特征变化的影响较小。B-V频率对尾涡下沉运动的影响较大。随B-V频率的增大,尾涡下沉速率明显减小,导致前后机的尾流安全间隔增大。

2) 大气分层主要作用于t′=1.5之后,会对尾涡的结构变化产生较大影响。无量纲B-V频率N*=1时,诱使尾涡上方出现小尺度的反向涡旋,促进尾涡形变溃散,加快尾涡强度衰减。

3) 尾涡发展前期,由于涡旋作用两涡间的温度会出现一定时间的明显高于或低于尾涡温度。随着尾涡的下沉运动,尾涡温度与外界温度差值逐渐增大,这会导致两涡涡核的运动呈现出不同方向。B-V频率越大,尾涡演化后期的涡核间距越小。