基于改进混合长度尺度的机翼延迟分离涡模拟

2022-10-11丘德新王良军

航空兵器 2022年4期

丘德新，王良军，张武, 3*

(1. 上海大学力学与工程科学学院，上海 200072； 2. 上海大学信息化工作办公室，上海 200444；3. 上海大学上海市应用数学和力学研究所，上海 200072)

0 引言

飞机在大迎角飞行时，往往伴随复杂的分离流问题，机翼表面会出现明显的涡破裂现象。飞机前飞过程中，机翼翼尖会产生持续且较强的涡，在机翼后缘远场处能观察到较为明显的翼尖涡旋。近些年，关于翼尖涡的风洞实验研究相对较多，Dghim等采用实验的方法研究了NACA0012机翼的翼尖涡与远场格栅生成湍流的相互作用。García-Ortiz等在固定迎角= 9°的条件下，研究了展向连续射流对NACA0012机翼尾涡的影响。Lee等研究了近地面处机翼翼尖涡的流动特性。Qiu等采用PIV技术，研究了NACA0015矩形机翼在尾迹六倍弦长范围内所产生翼尖涡的演化。除了实验方法，计算流体力学(Computational Fluid Dynamics，CFD)方法也是重要的研究手段之一，其不仅可以降低经济成本，还可以缩短研究周期。Lombard等介绍了基于大涡模拟方法的翼尖涡形成和演化的数值研究。Pereira等分别采用雷诺平均Navier-Stokes(Reynolds Averaged Navier-Stoke，RANS)方程和尺度解析模拟(Scale-Resolving Simulation，SRS)模型对翼尖涡流动进行了数值模拟。

目前工程上主要的研究方法是通过求解RANS方程进行数值模拟，但RANS只能对湍流的平均矢量场和标量场进行预测，而且RANS计算的涡粘性较大，不能预测湍流的脉动场，对研究细小涡结构和大迎角分离流等问题仍具有局限性。直接数值模拟(Direct Numerical Simulation，DNS)能获得几乎所有尺度的流动信息，但其要求极高的网格分辨率，计算量巨大。而大涡模拟(Large Eddy Simulation，LES)方法与RANS和DNS不同，该方法可以将流动分为大尺度流动和小尺度脉动。大尺度区流动受外部因素的影响较为强烈，需进行直接求解，而对小尺度区的湍流脉动进行模化。LES采用的是空间滤波技术，过滤掉小于截止尺度的涡，能够有效地处理圆柱绕流、大迎角分离流等复杂湍流问题。然而，在边界层内层涡尺度较小，难以模化，其大部分湍流能量主要集中在小尺度涡中，因此，计算量仍然十分巨大。考虑到计算资源限制及节省计算成本，同时为了精确地模拟复杂流动问题，采用RANS-LES混合方法是一种较好的选择。该方法结合了RANS在边界层内层计算效率较高以及LES在大分离区精度较高的优点，是目前计算大迎角分离流、自由射流等问题应用较为广泛的方法。

分离涡模拟(Detached Eddy Simulation，DES)类方法是RANS-LES混合方法中比较有代表性的一种方法。1997年，Spalart等首次提出基于Spalart-Allmaras(S-A)湍流模型的DES方法。但DES存在一定的缺陷，若网格加密不当，在极端情况下会出现网格诱导分离(Grid Induced Separation，GIS)等现象。随着研究的不断深入，当分离区较窄或边界层较厚时，由于LES对边界层的入侵，可能会引发模化应力衰减(Modeled-Stress Depletion，MSD)等现象。为了消除此类现象，Spalart等对原来的DES进行了改进，提出延迟分离涡模拟(Delayed Detached Eddy Simulation，DDES)的方法，在一定程度上防止了LES入侵到边界层。

大量的研究表明，DES在复杂分离流等问题中的应用较为成功，相比于LES，其计算成本较小。虽然原始DDES在一定程度上改善了“灰区”问题，在计算中小迎角分离流、强剪切流等问题时，壁面边界层以及邻近区域的网格通常为各项异性网格，其表现为法向网格尺寸远小于流向和展向的网格尺寸。由于DDES的网格尺度采用的是网格边长的最大值，有可能造成剪切层区域模化的涡粘性过大，进而导致RANS到LES转换的显著延迟。因此，为了改善DDES从RANS过渡到LES的特性，本文基于开源平台OpenFOAM，对DDES的混合长度尺度进行改进，研究NACA0012半展长矩形机翼在产生分离流时的湍流涡结构以及翼尖涡的生成发展。

1 数值方法

1.1 DDES方法

基于S-A湍流模型的DES方法是在S-A湍流模型的基础上，对破坏项进行了修改，将壁面距离替换为广义长度尺度：

=min(,),=max(Δ, Δ, Δ)

(1)

修改后的长度尺度使得S-A湍流模型只在壁面附近起作用，而在其他区域转变为与网格尺度相关的类Smagorinsky亚格子应力模型。

基于S-A湍流模型的DES方程为

(2)

为了防止LES提前进入边界层，Spalart等在原DES长度尺度中引入边界层识别函数，得到混合长度尺度：

=-max(0,-)

(3)

(4)

在计算分离流问题时，RANS到LES的转换需要一定的时间和空间的发展，尤其在分离初期，由于二维剪切层区域模化的涡粘性较大，非定常的特性较弱，延迟了RANS到LES的转换。为了改善分离初期RANS到LES的转换特性，降低特定区域的涡粘性，Shur等对RANS-LES混合方法中的LES长度尺度进行改进，首先对网格尺度进行修正，得到新的网格尺度：

(5)

式中：=/||×，为涡量。

(6)

(7)

(8)

1.2 RANS-LES混合长度尺度修正

(1) 对于边界层附近的“扁平”网格，采用上述网格尺度：

(9)

(2) 对于四面体、金字塔形、三棱柱、正六面体等各项同性网格，则将网格尺度定义为网格单元体积的立方根, 即

(10)

图1为网格尺度类型判定流程图。计算网格为任意网格类型的非结构化网格，为了识别出边界层附近的各向异性网格，首先对每个控制体单元的面数进行判定。若控制体单元面数为6，则根据网格长宽比判定是否为各向异性网格。网格长宽比大于阈值，则采用网格尺度，否则均采用网格尺度。

图1 网格判断及混合网格尺度流程图

1.2.2 涡探测函数修正

为了避免()在远场无粘区及湍流边缘区开启，Shur等对进行了修正：

(11)

考虑到边界层附着流的数值安全性，避免()在边界层区域中开启，引入式(4)延迟函数中的，对式(11)进行修正：

tanh[()]

(12)

本文对OpenFOAM中现有的S-A DDES湍流模型进行修改，植入了和能够分辨二维剪切层的函数()，并实现了混合网格尺度方法。

2 数值实验与分析

在128核的高性能计算服务器下进行计算，总物理内存为256 GB，处理器型号为AMD EPYC 7702。

2.1 NACA0012翼型验证

为验证改进混合长度尺度的有效性，采用基于修正长度尺度的S-A DDES方法对NACA0012翼型进行计算。数值实验马赫数为0.15，来流迎角为18°，弦长及展长均为1 m。基于弦长的雷诺数为6×10，计算网格如图2所示，采用非结构混合网格，并在NACA0012局部区域进行网格加密，以更好地捕捉流动细节，网格单元总数约为200万。

图2 NACA0012翼型网格

图3给出了18°迎角的NACA0012翼型在某时刻和()的分布情况。根据分布云图可知，在背风区，对准二维剪切层的判断起到了作用，且有效避免了()在边界层区域开启，保证了边界层附着流的数值安全性。由于需要避免()在远场无粘区及湍流边缘开启，在这些区域引入的判断函数max(1.0,)值较大，因此，值也较大，在湍流边缘可以发现明显的边界。根据()分布云图可知，在流动分离初期，剪切层区域的值趋于阈值=0.1，在流动分离的下游大部分涡破碎区域及远场区域的值趋于阈值=1。

图3 VTMnew和FKH(VTMnew)分布云图

NACA0012翼型在18°迎角时已经为失速状态，表1给出了RANS以及基于的DDES方法与实验结果的升力系数对比，DDES方法所得的为计算稳定后的平均值，实验结果来自McCroskey所提供的数据。根据对比结果可知，由于流场具有明显的非定常特性，RANS对的预测有一定的延迟，而基于的DDES计算结果虽然略小于实验值，但误差相对较小。

表1 升力系数CL对比

图4为同一时刻不同位置的涡量(-vorticity)等值面，根据图4可以发现翼型背风面为大分离区。由于翼型表面附近具有更密集的网格，可以解析出更精细的涡结构。背风区的涡旋运动较为复杂且具有随机性。在背风区下游可以清晰观察到大涡结构以及正、负涡量交替脱落现象。

图4 瞬时涡量云图

2.2 半矩形翼计算

为研究小迎角分离以及翼尖涡的形成与演化过程，计算模型基于JAXA的NACA0012矩形机翼风洞实验模型，机翼弦长为0.4 m，半翼展长为1 m，其中翼尖为不带翼梢小翼的钝体翼尖，雷诺数为1.8×10，马赫数为0.175。根据实验研究表明，机翼迎角为12°时，在机翼尾缘附近能够观察到流动分离，因此，计算也采用12°迎角。

为了验证NACA0012矩形机翼的网格无关性，计算了三套不同疏密程度的非结构化网格，如表2所示，细网格的网格量约为粗网格的4.5倍。

表2 NACA0012机翼网格无关性验证

由表2可知，RANS方法计算的最大相对误差约为2.44%，基于的DDES方法计算的最大相对误差约为0.72%。对于细网格，两种方法计算的相对误差约为0.24%。可以认为计算结果达到了网格收敛，在之后的计算中主要采用中等网格和细网格。

为了提高网格类型复杂度，计算采用混合网格，在机翼表面附近采用六面体以及三棱柱网格，其他区域采用四面体以及金字塔形网格。图5(a)为对称面和机翼表面网格，图5(b)为机翼翼端处表面网格。为了更好地捕捉在机翼翼尖和后缘附近涡结构，对机翼所在区域尤其在机翼翼尖周围进行局部加密，网格数大约为210万，即中等密度网格。为了保证在机翼表面湍流附面层计算的准确性，第一层网格点与机翼表面之间的距离满足+小于1。

首先，采用RANS对NACA0012机翼进行计算，以获得机翼稳态流动特性。采用S-A湍流方程模型；对流项采用高斯迎风格式；扩散项采用高斯线性格式；矩阵求解器采用预条件共轭和稳定化预条件双共轭求解器；速度-压力耦合采用SIMPLE算法。

图5 NACA0012机翼网格

为减小计算时间和保证计算精度，将RANS的稳态结果作为DDES的初始条件，分别采用原S-A DDES方法和改进RANS-LES混合长度后的S-A DDES方法进行模拟。对流项采用混合格式，在涡主导的流动分离区域采用低耗散线性格式，而在远场无粘无旋区采用迎风格式，保持足够耗散以提高计算稳定性。扩散项采用高斯线性格式；时间项采用Crank-Nicolson方法进行离散；速度-压力耦合采用PIMPLE算法。为确保计算的稳定性，采用自动调整时间步方法并保证最大库朗数小于1。

图6为12°迎角下RANS、原S-A DDES、基于的DDES方法机翼表面及翼尖处的涡量等值面云图(=30 000 s)，采用流向的速度分量进行着色。其中，图6(a)RANS方法是在计算达到收敛时候的结果，而原DDES和基于的DDES方法的涡量等值面是在计算初期相同时刻的结果。比较图6(b)和图6(c)可以发现，分离初期由于原DDES方法在剪切层区域模化的涡粘性较强，抑制了RANS到LES的转换，以至于在机翼表面的涡结构与RANS方法计算的涡结构类似，仅在机翼尾缘处有少量细小涡结构。而基于的DDES方法实现了RANS到LES的快速转换，在机翼背风区可以清楚地观察到剪切层的快速失稳，旋涡已不再是一个整体而是分散的细小涡结构。对比图6在翼尖处的涡量等值面发现，三种方法的计算结果相似，都清晰地显示了涡在翼尖处卷曲并旋转拖出较长的尾涡。由于在翼尖区域的网格较为粗糙，基于的DDES方法的计算结果并未展现出更多的细节。

图6 12°迎角机翼表面Q涡量等值面

因此，为了进一步验证改进混合长度尺度的有效性，对计算网格进行加密，尤其在机翼背风区及翼尖处，减小了网格尺寸并扩大了加密区域。加密后的网格数大约为540万。图7为网格加密后的机翼翼根处截面压力系数计算结果对比。结果表明，基于的DDES方法和RANS方法的计算结果与实验结果较为吻合。图8(a)为基于的DDES方法计算得到的在翼尖周围的流线图。其中，红线表示主涡再附着线，蓝线表示次涡分离线。与图8(b)实验油流图对比可以看出，计算结果与实验结果较为吻合。

图7 翼根处机翼表面压力系数

图8 翼尖周围流线图对比

图9给出了细网格基于的DDES方法的机翼表面及翼尖处的涡量等值面云图(=30 000 s)。从图中可以看出，在翼尖处进行网格加密后得到的涡结构展现出更多的细节。根据已有研究可知，翼尖周围的流场与双涡结构有关。机翼上表面的涡称为主涡旋，翼尖和下表面的涡称为次级涡以及次级尾迹涡。机翼前缘翼尖下表面处的涡管比上表面的涡管强，下表面涡管逐渐在尾缘处分成数根涡管并逐渐上翻演变，在尾缘处与上表面涡管“拧”成螺旋状并在后方形成一条大涡管。翼尖涡结构脱离机翼表面并向下游远场处传播，在尾迹区某一位置形成闭环结构。

图9 基于lhyb_new的DDES方法的Q涡量等值面

图10给出了翼尖处不同截面流向涡量。从图中能够更清楚地观察到翼尖涡的形成细节，在机翼前缘翼端表面有正涡量产生，但从机翼下表面卷起的负涡量更为明显。机翼上表面的正涡量直径沿流向有一定的增大，而下表面负涡量反向旋转至上表面并且增长更为迅速。结合图9中的涡量等值面，下表面的涡管在机翼尾缘上翻形成螺旋状的不稳定现象，此时与周围的流体产生充分的能量交换。

图10 机翼翼端不同截面流向涡

由于翼尖处上下表面涡旋的相互作用，在机翼尾缘处形成了较为明显的尾迹涡，如图11所示。在脱离翼端尾缘后尾迹涡继续向上翻卷，涡核不断增大，在距离尾缘大约两倍弦长的位置，尾迹涡与涡量层发生分离，并在远场处逐渐耗散。图12显示了尾迹流向涡随不同截面的演化过程，在尾缘处能清楚地观察到涡量层正负涡量的交替过程。在机翼下游，涡量层不断向尾迹涡进行供给并沿轴负方向移动。

图11 近翼流场尾迹流向涡

为了验证本文方法在计算半展长矩形机翼大迎角分离流问题时的有效适应性，对18°迎角的NACA0012矩形机翼进行数值模拟。图13为RANS方法及基于的DDES方法计算得到的涡量等值面云图。根据计算结果可知，机翼在18°迎角时具有更明显的非定常特性，且比12°迎角的机翼具有更明显的涡结构。对比RANS方法的涡量等值面云图可以发现，改进后的DDES方法在机翼背风区能够捕捉更精细的涡系结构，且在翼尖处能够清晰地看到涡破裂现象，也进一步证实了该方法的有效性。