基于k-ω湍流模型的翼尖涡演化过程数值模拟

钱 宇,蒋 皓

(中国民用航空飞行学院飞行技术学院,广汉 618307)

当机翼产生正升力时,由于气体流过上下翼面产生不同的压强,从而导致气流的纵向偏移,并在翼尖脱落时形成翼尖涡[1],翼尖涡向后发展继而形成尾流迹。由于大气环境的特性,翼尖涡所产生的尾流会持续很长一段时间,其扰动的不稳定空气具有较大的动能,为后机飞行安全造成很大的事故隐患[2]。在中国经济高速发展的现在,有限的空域资源逐渐不能满足民航运输的需求,而尾流将直接影响到终端区的容量,影响飞机的起降架次,造成航班的延误。在旋翼的运动过程中,翼尖涡的存在将造成机翼的颤振,从而产生严重的噪声污染[3]。为解决上述问题,对于翼尖涡演化过程的研究显得尤为重要。

目前,对于翼尖涡的研究方法主要有风洞实验、激光雷达以及计算流体力学等一系列手段[4]。其中计算流体力学较其他方法成本低、周期短、计算精度高,且能较好地捕捉到翼尖涡的变化过程。对于翼尖涡演化过程的仿真,不同的模型方法运用于不同的计算环境,由此中外学者进行了大量的研究。杨祥生等[5]采用双方程SSTk-ω模型对风力机尾流流场进行了数值分析,以期证明该模型在翼尖涡预测方面具有较好的应用效果。Anusonti-Inthra等[6]利用雷诺平均和基于粒子的涡量输运方法,用来预测低速流动条件下孤立机翼的性能和尾流参数。Uhl等[7]采用DLR TAU流体求解器与雷诺应力模型结合,研究了不同翼尖形状对翼尖涡的演化过程的影响,并在此基础上预测了涡核轴向速度的变化规律。中国对于翼尖涡流迹的研究现大多数还处于试验或验证阶段[8],在仿真计算方面也多为研究单一翼型。

研究基于SSTk-ω湍流模型,对整机模型进行模拟仿真与数值计算,充分考虑机身、水平安定面、垂直安定面对于机翼近场翼尖涡的影响,得到机翼上两个共转融合涡的形成过程及其在近翼流场转化为尾流迹的连续演化过程,从气动力性能方面为近场及远场翼尖涡流迹的研究提供相应的理论基础,也能为民航制定相应的规章制度提供一定的技术参考。

1 模型与方法

1.1 计算模型

数值计算前需对计算模型进行三维建模,为保证结果的准确性,在采用SolidWorks对计算模型进行整体建模,所得计算模型如图1(a)所示,建模参数及运行实况如表1所示。考虑到计算流体动力学( computational fluid dynamics，CFD)方法对计算资源的要求较大,且主要研究近翼流场的翼尖涡,为提高计算效率,缩短计算时间,因此简化计算模型,即采用半模,保留机身、机翼、水平安定面及垂直安定面,如图1(b)所示。

图1 A320飞机等比例缩放与计算模型Fig.1 Proportional scaling model and computational model of A320 airplane

表1 模型尺寸大小及其运行实况Table 1 Model size and operating conditions

1.2 网格的划分

网格的划分将计算模型离散为若干网格,通过对每一个网格进行计算及分析,整合后可得计算模型的整体气动力性能。利用ICEM-CFD对计算模型进行结构化网格划分。结构化网格较非结构化网格而言更为精细,对于区域边界的拟合能更好地实现,且在网格生成的过程中速度快、质量好、表面光滑、更加贴近实际模型,能很好地适用于流体和表面应力的计算。结构化网格的网格质量较好,较好的网格质量能充分地提高数值计算过程中各项参数的收敛性,使得计算所得的结果更为准确、合理。划分后的计算网格如图2所示。

图2 计算流场及计算模型局部放大图Fig.2 Calculation flow field and local enlarged drawing of calculation model

在划分网格前,为保证模拟真实着陆过程,屏蔽其他因素的干扰。因此设置流场入口距离模型4c(c为翼展长度),流场出口离模型20c,上壁面离模型5c,下表面离模型1c。为了使计算结果更为精确,则机翼处需生成更加精细的网格。考虑到单纯增加网格节点数或者减小网格增长率会大大增加网格数量,为数值计算带来相应的难度。为了解决这个问题,采用两次O-Block方法对机翼附面层网格进行相应的加密处理。该方法不仅能对局部进行加密,而且可以较好地避免复杂形状的Block顶点处发生的网格扭曲,从而在机翼表面得到所需的边界加密层。O-Block方法划分的机翼网格如图3所示,参数如表2所示。

图3 O-Block划分机翼网格Fig.3 Generating wing mesh of O-Block

表2 网格划分参数Table 2 Meshing parameters

最后一定要保证附面层网格的高度,这将直接影响计算的收敛性,附面层网格高度y的计算公式[9]为

(1)

式(1)中:y+取1；μ为流体的动力黏度；Uτ为流体的估计速度；ρ为流体密度。

通过式(1)的计算,可以得到附面层网格高度为1.979×10-3mm,因此取y=1.5×10-3mm，满足计算要求。最终生成的结构化网格数总量为7.057×106。

上述网格数量是在满足附面层网格要求的基础上,通过网格增长率计算得到的网格数量。为验证网格的无关性,适量增加节点数目、减小网格增长率,以此增大网格数量,计算增大后的758×104网格,并与705×104网格所得计算结果做比较。结果显示758×104网格与705×104网格对于近翼流场涡量的相对误差小于5%,如图4所示,说明网格数量对于结果的影响已不敏感,故可采用705×104网格进行数值计算。

图4 不同网格数量近翼流场涡量对比图Fig.4 Vorticity diagram of near wing flow field with different mesh numbers

1.3 控制方程

控制方程采用常黏度条件下不可压缩流体的N-S(Navier-Stokes)方程,该矢量形式[10]为

(2)

1.4 湍流模型

对于湍流的计算方法,目前主要有雷诺时均(RANS)、尺度解析以及直接数值模拟(DNS)。选取雷诺时均,该方法满足计算要求且具有较为高效的计算速度。基于雷诺时均方法下常用的湍流模型有如下几种:Spalart-Allmaras模型、SSTk-ω模型、标准k-ω模型、Reynolds Stress等。

在数值计算的过程中,不同的模型具有不同的适应环境。Spalart-Allmaras模型广泛用于航空领域对空气动力学的计算,对于有壁面边界空气动力学流动有较好的计算结果,但是该模型为单一运输方程模型,有较大能量耗散,可以运用于粗网格的计算,所得计算精度不高。标准k-ε模型计算量适中,具有较高的计算精度,但在模拟旋流和扰流时有缺陷,整个计算过程中伴随较大的能量耗散,且在计算压力梯度较大的复杂流场时易失准,导致前端涡旋的过度生成至整片机翼上,从而影响计算结果的准确性。Reynolds Stress多用于强涡流,通过直接使用运输方程来求解雷诺应力,避免了其他模型的黏性假设。RANS模型最复杂,收敛性较差,占用过多的CPU、运算时间及运算内存,对于较为复杂的三维强流动较为适用。

对于SSTk-ε模型,该模型使用混合函数将标准的k-ε和标准的k-ε相结合,在近翼面区域计算时保留了标准k-ε湍流模型的特性,即适合于存在逆压梯度情况的边界层流动；在远场区域则按照标准k-ε湍流模型进行计算。由此该模型在计算流动分离和逆压方面有较好的表现,适用于高逆压梯度及剪切流的计算[11]。其数学表达式为

Sω+Dω

(4)

式中：Γ为有效扩散率；G、Y、S、D分别为生成项、扩散项、用户自定义项以及交叉扩散项；ρ为流体密度；k为紊流动能；ω为耗散率，其表达式为

(9)

σk=σω=2

(10)

式中：ε为紊流耗散项。将N-S方程与SSTk-ω湍流模型相结合,保证数值计算过程的整体封闭,即可对生成的结构化网格进行数值计算。

2 仿真计算与分析

2.1 方法可行性验证

采用Fluent对结构化网格进行数值求解,选用基于涡度的转捩修正的SSTk-ω湍流模型；流场采用理想气体,黏度选择sutherland；边界条件采用远场压力入口,设置标准大气压为101.325 kPa,速度为130 n mile/h,迎角为5.5°,边界出口选用压力出口。求解方法采用基于耦合的隐式求解,又由于一阶迎风夸大了扩散的影响,容易偏离真正的场分布,因此以一阶迎风迭代400步作为初始计算值,再转为二阶迎风迭代500步,总迭代步长为900步。最后对网格标准初始化后开始计算。

在数值计算的过程中,为了证明数值方法的可行性,运用上述方法计算来流速度124 n mile/h,雷诺数4.6×106,迎角10°工况下NACA0012翼型的静压力系数CP,并绘制翼展方向1/3处曲线图,与Chow等[12]实验进行对比验证,如图5所示。

图5 翼型静压力系数曲线图Fig.5 Diagram of static pressure coefficient on airfoil

计算结果显示,仿真值与实验值的平均误差小于0.05,所得结果置信度高,误差在接受范围之类,因此该数值计算方法可行。

2.2 机翼表面翼尖涡的演化

图6所示为机翼侧表面不同截面处涡量分布。通过该图可知翼尖涡在机翼上的演化大体可以分为4个阶段。

图6 机翼不同截面处涡量分布Fig.6 Vorticity distribution at different sections of airfoil

第1阶段,次涡旋于机翼前缘形成后,涡量沿前缘到后缘有不断增加的趋势,当涡旋位于x/b=0.20(b为翼弦长度)时,次涡旋开始从机翼表面脱落且向上翼面不断移动,涡量表现出减小的趋势。第2阶段,主涡旋在x/b=0.20时形成,且有一定的增加趋势,但涡旋的当量值总体来说小于次涡旋。第3阶段,当次涡旋持续后移,在位于x/b=0.6时次涡旋开始和主涡旋相互纠缠融合,逐渐形成同一个涡旋,即第1共转融合涡。第4阶段为第2共转融合涡的形成,在位于x/b=0.90时机翼表面次涡旋和主涡再一次形成,并于x/b=0.95时再一次纠缠融合,形成第2共转融合涡。

图7所示为机翼表面涡核静压力系数变化曲线,从图7可以看出,次涡旋与次级尾迹涡的涡核静压力系数有一个当量先减小后增加的过程,原因是因为机翼前缘附近有较大的逆压梯度差,在位于x/b=0.05和x/b=0.10时存在一对旋向相反的涡旋。而当位于x/b=0.15时,正向的涡旋消失,从而导致负向涡旋增加,至此次涡旋的涡核静压力系数开始不断上升。

图7 机翼表面涡核静压力系数变化曲线Fig.7 Variation curve of vortex core static pressure coefficient on wing surface

第1共转融合涡在形成后其涡核的静压力系数的数值变化较小,且涡核的静压力系数要小于第2共转融合涡,主要是因为第1共转融合涡形成的时间及距离较第2共转融合涡更长,因此有一定的能量耗散,且在形成第2共转融合涡时,主涡所提供的能量要大于第1次融合时所提供的能量。

综上,由主涡和次级涡两次纠缠融合形成的两个共转融合涡旋与次级尾迹涡将共同影响近翼流场涡旋的形成。

2.3 近场翼尖涡的演化过程

通过上述涡旋的相互作用产生的负向涡旋流动,最终在x/c=0.10处形成一个完整的翼尖涡,如图8所示。在模拟的过程中也可以发现,来流对近场涡旋移动有较大的影响,无风情况下,涡旋应有下降的趋势,而5.5°的迎风来流使得翼尖涡有一定程度的向上偏移,因此翼尖涡在近翼流场的演化也很大程度上受环境因素的影响。

图8 近翼流场翼尖涡的涡量分布Fig.8 Vorticity distribution of wingtip vortex at near wing flow field

由图4、图9可知,翼尖涡在近场的形成及传递会有一定的能量转换。当涡旋从翼面脱离后,在近翼流场形成涡旋的过程中,近场涡核的静压力系数呈二次曲线减少。在x/c=0.40的位置时,涡核的静压力系数才逐渐趋于平稳值,直到在x/L=0.5(L为机身长度)的位置时,其数值变化较小呈现稳定状态。而翼尖涡也随着涡核静压力系数的减小,其体积在不断膨胀增加,继而逐步生成一个稳定的翼尖涡。

图9 不同截面处静压力系数分布云图Fig.9 Static pressure coefficient distribution at different sections

综上,在机翼表面形成的融合涡旋受次涡旋、主涡旋共同作用。且融合后涡旋直至离开机翼表面无明显的能量损耗。在负向涡流脱离机翼后,翼尖涡在近翼流场的形成与运动过程受环境因素影响较大。翼尖涡的生成过程中,伴随着明显的能量转换。

3 结论

通过仿真计算,获得了机翼表面及近翼流场翼尖涡的演化过程,得到了如下结论。

(1)O-Block方法能得到较好的近壁面网格,有利于近翼流场的计算；SSTk-ω湍流模型将标准k-ε模型和标准k-ω模型相结合,能较好地适用于涡旋演化过程的数值计算。

(2)在近翼流场的数值计算过程中,翼尖涡演化至x/L=0.50时才有较为稳定的涡核静压力,该处的计算参数则能成为远涡流场的初始计算值。在着陆阶段,主要的影响因素是大气的质量与空气的流动。对于远场尾迹的计算则需要考虑地效等因素,这也是今后研究远场尾涡演化过程的重点。