基于改进混合长度尺度的机翼延迟分离涡模拟
作者： 丘德新 王良军 张武
摘 要：对于中小迎角分离问题，延迟分离涡模拟（Delayed Detached Eddy Simulation，DDES）的方法在计算初期，流动分离区需要一定时间的发展，通常会延迟RANS模式到LES模式的切换，从而低估流动分离。本文采用改进RANS-LES混合长度尺度lhyb_new的方法，对网格尺度和涡探测函数VTM进行修正，基于计算流体力学开源平台OpenFOAM，将修正后的混合网格尺度、涡探测函数VTMnew以及二维剪切层识别函数FKH（VTMnew）植入到S-A DDES模型中，研究了大迎角NACA0012翼型以及12°迎角下的NACA0012半展长矩形机翼绕流问题，并对比了采用原S-A DDES和改进RANS-LES混合长度尺度后的S-A DDES湍流模型计算的机翼背风区涡结构和翼尖涡。结果表明，采用本文方法可以在机翼背风区观察到剪切层的快速失稳，并且能够更加精细地模拟机翼背风区的大尺度非定常涡系结构以及翼尖涡的形成细节，有效地加快DDES方法的RANS模式切换到LES模式。此外，本文还计算了该矩形机翼在18°迎角的分离流问题，进一步验证了所提方法的有效适应性。
关键词：延迟分离涡模拟;非定常流动;RANS; LES;翼尖涡;OpenFOAM; 空气动力学
中图分类号：TJ760;V211.41
文献标识码：A
文章编号：1673-5048（2022）04-0070-07
DOI：10.12132/ISSN.1673-5048.2021.0239
0 引 言
飞机在大迎角飞行时，往往伴随复杂的分离流问题，机翼表面会出现明显的涡破裂现象[1]。飞机前飞过程中，机翼翼尖会产生持续且较强的涡，在机翼后缘远场处能观察到较为明显的翼尖涡旋。近些年，关于翼尖涡的风洞实验研究相对较多，Dghim等[2]采用实验的方法研究了NACA0012机翼的翼尖涡与远场格栅生成湍流的相互作用。García-Ortiz等[3]在固定迎角α=9°的条件下，研究了展向连续射流对NACA0012机翼尾涡的影响。Lee等[4]研究了近地面处机翼翼尖涡的流动特性。Qiu等[5]采用PIV技术，研究了NACA0015矩形机翼在尾迹六倍弦长范围内所产生翼尖涡的演化。除了实验方法，计算流体力学（Computational Fluid Dynamics，CFD）方法也是重要的研究手段之一，其不仅可以降低经济成本，还可以缩短研究周期。Lombard等[6]介绍了基于大涡模拟方法的翼尖涡形成和演化的数值研究。Pereira等[7]分别采用雷诺平均Navier-Stokes（Reynolds Averaged Navier-Stoke，RANS）方程和尺度解析模拟（Scale-Resolving Simulation，SRS）模型对翼尖涡流动进行了数值模拟。
目前工程上主要的研究方法是通过求解RANS方程进行数值模拟，但RANS只能对湍流的平均矢量场和标量场进行预测[8]，而且RANS计算的涡粘性较大，不能预测湍流的脉动场，对研究细小涡结构和大迎角分离流等问题仍具有局限性。直接数值模拟（Direct Numerical Simulation，DNS）能获得几乎所有尺度的流动信息，但其要求极高的网格分辨率，计算量巨大。而大涡模拟（Large Eddy Simulation，LES）方法与RANS和DNS不同，该方法可以将流动分为大尺度流动和小尺度脉动。大尺度区流动受外部因素的影响较为强烈，需进行直接求解，而对小尺度区的湍流脉动进行模化。LES采用的是空间滤波技术，过滤掉小于截止尺度的涡，能够有效地处理圆柱绕流、大迎角分离流等复杂湍流问题。然而，在边界层内层涡尺度较小，难以模化，其大部分湍流能量主要集中在小尺度涡中[9]，因此，计算量仍然十分巨大。考虑到计算资源限制及节省计算成本，同时为了精确地模拟复杂流动问题，采用RANS-LES混合方法是一种较好的选择。该方法结合了RANS在边界层内层计算效率较高以及LES在大分离区精度较高的优点，是目前计算大迎角分离流、自由射流等问题应用较为广泛的方法。
分离涡模拟（Detached Eddy Simulation，DES）类方法是RANS-LES混合方法中比较有代表性的一种方法。1997年，Spalart等首次提出基于Spalart-Allmaras（S-A）湍流模型的DES方法[10]。但DES存在一定的缺陷[11]，若网格加密不当，在极端情况下会出现网格诱导分离（Grid Induced Separation，GIS）等现象[12]。随着研究的不断深入，当分离区较窄或边界层较厚时，由于LES对边界层的入侵，可能会引发模化应力衰减（Modeled-Stress Depletion，MSD）等现象。为了消除此类现象，Spalart等对原来的DES进行了改进，提出延迟分离涡模拟（Delayed Detached Eddy Simulation，DDES）的方法[13]，在一定程度上防止了LES入侵到边界层。航空兵器 2022年第29卷第4期
丘德新，等：基于改进混合长度尺度的机翼延迟分离涡模拟
大量的研究表明[14-16]，DES在复杂分离流等问题中的应用较为成功，相比于LES，其计算成本较小。虽然原始DDES在一定程度上改善了“灰区”问题，在计算中小迎角分离流、强剪切流等问题时，壁面边界层以及邻近区域的网格通常为各项异性网格，其表现为法向网格尺寸远小于流向和展向的网格尺寸。由于DDES的网格尺度采用的是网格边长的最大值，有可能造成剪切层区域模化的涡粘性过大，进而导致RANS到LES转换的显著延迟。因此，为了改善DDES从RANS过渡到LES的特性，本文基于开源平台OpenFOAM，对DDES的混合长度尺度进行改进，研究NACA0012半展长矩形机翼在产生分离流时的湍流涡结构以及翼尖涡的生成发展。
1 数值方法
1.1 DDES方法
基于S-A湍流模型的DES方法是在S-A湍流模型的基础上，对破坏项进行了修改[10]，将壁面距离dw替换为广义长度尺度lDES：
lDES=mindw， CDESΔ， Δ=max（Δx， Δy， Δz）（1）
修改后的长度尺度lDES使得S-A湍流模型只在壁面附近起作用，而在其他区域转变为与网格尺度Δ相关的类Smagorinsky亚格子应力模型。
基于S-A湍流模型的DES方程为
图4为同一时刻不同位置的涡量（z-vorticity）等值面，根据图4可以发现翼型背风面为大分离区。由于翼型表面附近具有更密集的网格，可以解析出更精细的涡结构。背风区的涡旋运动较为复杂且具有随机性。在背风区下游可以清晰观察到大涡结构以及正、负涡量交替脱落现象。
2.2 半矩形翼计算
为研究小迎角分离以及翼尖涡的形成与演化过程，计算模型基于JAXA的NACA0012矩形机翼风洞实验模型[17]，机翼弦长为0.4 m，半翼展长为1 m，其中翼尖为不带翼梢小翼的钝体翼尖，雷诺数为1.8×106，马赫数为0.175。根据实验研究表明，机翼迎角为12°时，在机翼尾缘附近能够观察到流动分离，因此，计算也采用12°迎角。
为了验证NACA0012矩形机翼的网格无关性，计算了三套不同疏密程度的非结构化网格，如表2所示，细网格的网格量约为粗网格的4.5倍。
由表2可知，RANS方法计算的CL最大相对误差约为2.44%，基于lhyb_new的DDES方法计算的CL最大相对误差约为0.72%。对于细网格，两种方法计算的CL相对误差约为0.24%。可以认为计算结果达到了网格收敛，在之后的计算中主要采用中等网格和细网格。
为了提高网格类型复杂度，计算采用混合网格，在机翼表面附近采用六面体以及三棱柱网格，其他区域采用四面体以及金字塔形网格。图5（a）为对称面和机翼表面网格，图5（b）为机翼翼端处表面网格。为了更好地捕捉在机翼翼尖和后缘附近涡结构，对机翼所在区域尤其在机翼翼尖周围进行局部加密，网格数大约为210万，即中等密度网格。为了保证在机翼表面湍流附面层计算的准确性，第一层网格点与机翼表面之间的距离满足y+小于1。
首先，采用RANS对NACA0012机翼进行计算，以获得机翼稳态流动特性。采用S-A湍流方程模型;对流项采用高斯迎风格式;扩散项采用高斯线性格式;矩阵求解器采用预条件共轭和稳定化预条件双共轭求解器;速度-压力耦合采用SIMPLE算法。
为减小计算时间和保证计算精度，将RANS的稳态结果作为DDES的初始条件，分别采用原S-A DDES方法和改进RANS-LES混合长度后的S-A DDES方法进行模拟。对流项采用混合格式，在涡主导的流动分离区域采用低耗散线性格式，而在远场无粘无旋区采用迎风格式，保持足够耗散以提高计算稳定性。扩散项采用高斯线性格式;时间项采用Crank-Nicolson方法进行离散;速度-压力耦合采用PIMPLE算法。为确保计算的稳定性，采用自动调整时间步方法并保证最大库朗数小于1。
图6为12°迎角下RANS、原S-A DDES、基于lhyb_new的DDES方法机翼表面及翼尖处的Q涡量等值面云图（Q=30 000 s-2），采用流向的速度分量进行着色。其中，图6（a）RANS方法是在计算达到收敛时候的结果，而原DDES和基于lhyb_new的DDES方法的Q涡量等值面是在计算初期相同时刻的结果。比较图6（b）和图6（c）可以发现，分离初期由于原DDES方法在剪切层区域模化的涡粘性较强，抑制了RANS到LES的转换，以至于在机翼表面的涡结构与RANS方法计算的涡结构类似，仅在机翼尾缘处有少量细小涡结构。而基于lhyb_new的DDES方法实现了RANS到LES的快速转换，在机翼背风区可以清楚地观察到剪切层的快速失稳，旋涡已不再是一个整体而是分散的细小涡结构。对比图6在翼尖处的Q涡量等值面发现，三种方法的计算结果相似，都清晰地显示了涡在翼尖处卷曲并旋转拖出较长的尾涡。由于在翼尖区域的网格较为粗糙，基于lhyb_new的DDES方法的计算结果并未展现出更多的细节。
因此，为了进一步验证改进混合长度尺度的有效性，对计算网格进行加密，尤其在机翼背风区及翼尖处，减小了网格尺寸并扩大了加密区域。加密后的网格数大约为540万。图7为网格加密后的机翼翼根处截面压力系数计算结果对比。结果表明，基于lhyb_new的DDES方法和RANS方法的计算结果与实验结果较为吻合。图8（a）为基于lhyb_new的DDES方法计算得到的在翼尖周围的流线图。其中，红线表示主涡再附着线，蓝线表示次涡分离线。与图8（b）实验油流图对比可以看出，计算结果与实验结果较为吻合。
图9给出了细网格基于lhyb_new的DDES方法的机翼表面及翼尖处的Q涡量等值面云图（Q=30 000 s-2）。从图中可以看出，在翼尖处进行网格加密后得到的涡结构展现出更多的细节。根据已有研究可知，翼尖周围的流场与双涡结构有关。机翼上表面的涡称为主涡旋，翼尖和下表面的涡称为次级涡以及次级尾迹涡。机翼前缘翼尖下表面处的涡管比上表面的涡管强，下表面涡管逐渐在尾缘处分成数根涡管并逐渐上翻演变，在尾缘处与上表面涡管“拧”成螺旋状并在后方形成一条大涡管。翼尖涡结构脱离机翼表面并向下游远场处传播，在尾迹区某一位置形成闭环结构。