可压缩湍流边界层近壁区马蹄涡的演化

史万里，葛 宁，陈 林，唐登斌

(1.南京航空航天大学能源与动力学院，江苏南京 210016;2.南京航空航天大学航空宇航学院，江苏南京 210016)

0 引言

在流体力学的研究领域中，湍流边界层问题是一个至今尚未解决的重大难题，而拟序结构问题则对于湍流边界层的进一步深入研究和探索至关重要［1］，尤其在湍流边界层近壁区域，流动的许多宏观特性，如湍流的生成和维持、湍流能量的输运和耗散、阻力减小、热量传递等都与之密切相关［2］。大量实验研究与数值模拟结果表明，流向涡结构是拟序结构的一种重要典型特征［3］，马蹄涡(或发卡涡)更受到了广泛的重视［4－5］，其生成与演化规律的研究，能够为湍流边界层内的流动现象提供很多合理的解释，并有助于深入理解和解决工程中遇到的流动问题。

在湍流边界层中，马蹄涡与发卡涡在结构上稍有不同，通常在动量厚度雷诺数大于800时称为发卡涡，低于此值时涡头部分通常较宽，因此称为马蹄涡［6－7］。以前人们认为这种马蹄涡(或发卡涡)结构可用对称的一对反向旋转的流向涡对来表示［8］。最近的研究发现［2，6－7，9］，湍流近壁区的拟序结构不仅存在对称的流向涡对，而且绝大多数流向涡对更是以非对称形式存在，这种结构有时甚至是以一个流向涡的形式出现，这种现象引起人们的很大兴趣。

在不同的数值方法中，大涡模拟(LES)方法相比于雷诺平均方法(RANS)，能够显示出流场的细节和瞬态过程［10］，而对计算资源的需求大大小于直接数值模拟(DNS)方法，所以本文是采用LES方法研究三维马蹄涡的结构。

对于采用LES方法的湍流模拟，如果直接从一般层流开始转捩到湍流，那么计算湍流场则需要附加很大的计算量。因此本文采用进口给定三维扰动波方法［11－12］，作为大涡模拟的湍流进口边界的给定方法，模拟了可压缩三维槽道湍流流动。根据计算结果，与DNS解的数据进行了比较，并对近壁湍流的马蹄涡结构的演化进行了深入的研究和分析。

1 控制方程

对预处理下［13－14］的可压缩流 Navier-Stokes方程，采用Faver提出的密度加权过滤方法，得到曲线坐标系下，无量纲通量形式的LES控制方程为:

式中，上标“－”表示系统平均，“～”表示细网格过滤，转换矩阵 M= ∂Q/∂q，原始变量为为守恒变量。Re为特征雷诺数表示亚格子粘性应力为亚格子热流通量，雅克比矩阵转换矩阵

需要注意的是，Γ=diag{1 1 1 1 β(Mr)}为预处理矩阵［13－14］，预处理参数 β 的选取:若 Mr＜1.0，β(Mr)=Mr2;若Mr≥1.0，β(Mr)=1.0。其中为特征马赫数。这里采用预处理技术，主要是为了改善边界层近壁区低马赫数流动情况下的计算格式的稳定性以及加快计算结果的收敛性。

2 数值方法

将对控制方程使用有限体积法进行离散，粘性项采用二阶中心格式，对流项采用六阶精度对称WENO格式，时间推进采用三步Runge-Kutta法进行精确求解。其亚格子模型、计算域和进口边界条件分述如下:

2.1 亚格子模型

建立合理的亚格子模型是大涡模拟的关键，Smagorinsky模型是广泛使用的一种亚格子模型，表达式如下［15］:

相应的亚格子热流通量为:

但是该模型存在着某些缺陷，如系数C实际最佳值取决于运动流体的性质与状态而不是常数，且模型未体现流体压缩性的影响等，因此本文采用Germano提出的动态亚格子模型，具体如下［15］:

动态模型是动态确定亚格子涡粘系数C，对湍流场做粗细两次过滤，并假设粗网格上最小脉动产生的应力等于粗细网格分别过滤产生的亚格子应力之差，在此基础上用最小二乘法推导出亚格子模型中的系数C如下:

其中:

式中，上标“^”表示粗网格过滤，“＜ ＞”表示在统计方向上平均量。

2.2 计算域

计算域见图1(a)，计算中选取进口湍流边界层的位移厚度δ*作为特征长度，计算域无量纲尺度取为50×10×8，分别对应流向x、法向y和展向z的长度。三个方向上计算网格点数分别取为192×96×64，在流向和展向上计算网格平均分布，而在法向上的近壁区域内进行加密。按照进口边界条件确定的网格尺度为△x+=13.0，△z+=6.25，第一层网格的y+=0.23，粘性底层区域(即 y+=10)内有22个网格点，x-y面局部网格如图1(b)。为了方便并行计算将网格平均分为4份，每份在一个CPU上计算，编号依次从0到3，如图1(a)，利用MPI来传递各区之间信息交换，并行计算。

流场初始化在边界层内采用Spalding湍流平均速度剖面。在计算域的出口给定反压，两边的展向位置给定周期性边界条件，在上边界采用对称边界条件，平板壁面处采用绝热无滑移边界条件。因此图1中所给出的仅是法向y上的半个槽道区域。

图1 计算域及网格示意图Fig.1 Schematic of the parallel computational domain decomposition

2.3 进口边界条件

在计算域进口采用了进口给定三维扰动波方法，这是基于目前湍流研究者普遍认同的湍流是拟序结构的组合这一假设，在进口产生包括时间和空间变化的拟序结构，即扰动波。这种方法相比于其它常用的湍流进口给定方法，是简单可行的，其进口发展段更短，更适用于工程上使用。

进口边界条件，主要是确定进口处的瞬时速度。进口处瞬时速度通常由三部分组成:瞬时值=平均值+扰动值+随机噪音。平均值是采用Spalding湍流平均速度剖面，随机噪音的最大幅值不超过平均自由流速度的1%。

对于扰动值的确定，本文直接给定三维波。通常湍流边界层按动态特征可分为内层和外层，这种方法是认为内层中存在低速条带，用一道三维波表示，而把外层中的扰动波用三维基波和亚谐波组合的涡结构［11－12］来表示。其公式［11］及参数的选择如下:

式中，展向速度w按照自由散度条件由u和v解出，基本参数参照表1。表中的ωj是脉动频率，βj是展向波数，φj是相位。

表1 湍流边界层进口湍流脉动参数Table 1 Parameters of inflow turbulent fluctuations for the turbulent boundary-layer

3 结果与讨论

湍流边界层的结构是三维的，边界层外缘具有弱的自由剪切层的性质，内外层之间有相互作用，内层又受到壁面的抑制，因此湍流边界层内的流动结构非常复杂。为了检验计算结果的可靠性，这里首先把计算结果与Spalart的DNS解［16］等进行了比较，然后对LES计算得到的马蹄涡结构进行分析。文中算例的马赫数Mr=0.6，基于动量厚度和自由流速度的雷诺数为670(对应于基于位移厚度和自由流速度的雷诺数为1000)。

3.1 时间平均速度和表面摩擦系数

图2(a)是LES计算得到的平均速度剖面图，同时给出了Spalart的DNS的结果［16］。由图可以看出，本文的LES计算结果与相同雷诺数下的DNS精确解符合的很好，说明了所采用的进口剖面在经过很短的发展段之后，已发展成稳定的湍流。

图2(b)是湍流脉动速度，即湍流强度分布图，其纵坐标是脉动速度的绝对值(已用自由流速度无量纲化)。该图也将LES计算值与Spalart的湍流DNS解进行了比较，发现各个方向上的湍流强度的幅值与DNS解基本一致，仅在流向脉动速度的幅值稍有偏离，其原因可能是DNS解为不可压缩流计算的结果，而本文的LES解是预处理下可压缩流计算的结果。

图2(c)是时间平均和展向平均的表面摩擦系数Cf。可以看到从计算域入口仅经过很短的入口过渡距离，大约在流向Rex为139000时，计算的Cf就符合了湍流的经验值。在图2中，湍流摩擦系数Cf并不是完全符合湍流的经验值，而是存在波动，这是由于湍流的非定常性引起的。

3.2 湍流边界层的马蹄涡结构

本节对计算得到的湍流边界层内马蹄涡结构的演化进行分析，其中给出的时间是以进口标准周期T为单位。因为涡核处存在强的涡量及最小压力值，所以采用拉普拉斯压力［17］P，kk=(ωiωi/2 － SijSij)的正的等值面是很容易识别拟序结构的，以下图中涡结构显示均为P，kk=0.6 的等值面。

3.2.1 马蹄涡的形成和演化

图3是t=6.62T时的速度云图和流场涡结构，图中深蓝色区域为低速区域。图中可以清楚地看到流场存在很多的流向涡结构，这些流向涡结构与低速条带是相对应的，流向涡在x-z平面内是左右摆动、交替出现的，因而高低速条纹是左右摆动而不单纯是流向的。这和文献［9］和文献［18］中的结论一致。

图4中显示了对称流向涡对形成马蹄涡结构的过程，这个结构在后面的图8中的标号为E1。正如图3中提到的，流向涡结构左右摆动，在图4中可以看到，当两条流向涡相互靠近时，两涡腿的头部之间会形成桥状的展向涡结构，进而合并融合而成的，这种对称产生的马蹄涡结构形成过程与文献［5］的DNS解一致。

图2 LES计算结果与DNS解和经验公式的对比Fig.2 Comparison of result by LES，DNS and experiential formula

图3 流场涡结构与低速条带俯视图Fig.3 Top view of vortex structures and low-speed streaks in flow field

图4 对称流向涡形成的马蹄涡结构Fig.4 Generation of the horseshoe-shaped structures by symmetric streamwise vortex

两条涡腿相互靠近会形成带两条长涡腿的马蹄涡结构，文中的计算结果也清楚地展示了这一点，但是如果两条涡腿距离较远，则马蹄涡结构的涡头不会封闭，出现中间开口的现象，这样最终形成由单侧流向涡腿形成的马蹄涡结构。图5显示了F1结构由沿流向倾斜的流向涡发展为马蹄涡的过程，F1马蹄涡的涡腿是在右侧。图5(a)中，F1为流向涡，其图中显示的流向涡的涡头部分在自诱导的作用下向上抬头，同时向一侧的展向倾斜。图5(b)、图5(c)中可以看到这个向展向倾向的涡头沿展向拉伸，因此涡头部分的涡管变细，沿展向的长度变长。在图5(c)、图5(d)中可以看到，涡头的展向部分在拉伸的同时，在扰动的作用下，形成小的马蹄涡。在图5(e)～图5(h)可以看到这个单侧涡腿的马蹄涡较短涡腿的形成过程，它正如Robinson所述的带有单侧涡腿的马蹄涡形成过程那样［6］(如图6(a))，马蹄涡的颈部向壁面移动，逐渐下降并拉伸成一个加长的涡腿。在图5(e)、图5(f)，可以看到在马蹄涡的短涡腿侧出现了细碎的涡管a1，这是在颈部的移动、下降和拉伸时，这个被加长的涡腿的涡量较弱，在图中不能完全显示出来，因此看上去这个涡腿是断开的。图6(b)是Delo等人［7］在实验中得到了类似马蹄涡和流向涡的结果。图中标号为作者对流场中可能的涡结构的注释，其中V6是指流向涡，V4是指可能的单腿马蹄涡，这个马蹄涡的涡腿在左侧。

图5 单侧流向涡形成的较小马蹄涡对结构Fig.5 Generation of the horseshoe-shaped vortex pair by the one-side streamwise vortex

图6 湍流边界层中马蹄涡Fig.6 The horseshoe-shaped vortex in the turbulent boundary layer

由于文中进口边界条件所产生的是反向旋转的流向涡对，因此计算结果中这种单侧涡腿形成的马蹄涡结构总是成对出现的，如在后面的图8中标识的F1和F2及F3和F4结构。图7(a)是 t=5.94T时刻，F1和F2马蹄涡的俯视图和侧视图，俯视图中的背景为速度云图，图中深蓝色区为低速区。在俯视图中可以看到F1和F2结构都骑跨在低速条带上，在侧视图中可以看到马蹄涡头部的抬起。

图7(b)是t=5.98T时刻的典型马蹄涡结构涡腿处速度矢量场。可以看到，F1结构中环绕马蹄涡结构两条流向涡腿的流体是反向旋转的，因此在两条流向涡腿之间的流体存在强烈的上升，而在两条涡腿外侧，流体下扫。同时还注意到，在两条涡腿之间还存在一些碎的流向涡结构，而流体在流向涡两侧同样也形成上升和下扫运动。

3.2.2 亚谐波影响下的马蹄涡结构

图7 单腿马蹄涡细节图Fig.7 Particular view of one-side horseshoe vortices

图8显示的是流场中存在的拟序结构，为更清楚起见，把计算结果沿展向扩展了一个周期。我们注意到，图8中计算域的前面部分(即图中点划线之前)为进口发展段，这里的拟序结构只是研究在这之后(即点划线之后)的结果。图中显示了在湍流场中存在很多马蹄涡，且这些马蹄涡结构的形状并不是规则一致的，也不存在严格的对称性。文献［2］、文献［9］指出，这是由于湍流的非定常性影响，即湍流脉动会扭曲涡管的发展，从而很难在湍流边界层中得到完整一致的马蹄涡。

图8(a)、图8(b)分别显示了两个不同时间的流场拟序结构图，图中标示的大写英文字母序列表示在亚谐波影响下出现的马蹄涡结构。如图8中所示，可以发现近壁区存在几种对称流向涡对形成的马蹄涡的复杂结构(如E1到E3所标识的拟序结构，这种结构的涡头部分的形成如图4中所示)，但是基本的流动特征还是表示为亚谐波失稳引起的对称产生的马蹄涡交错出现。即沿流向x方向看，前一个对称产生的马蹄涡结构中心位置与后一个对称产生的马蹄涡结构中心位置沿流向并不是在一条直线上的，前排结构与后排结构的展向位置呈叉排交错分布。另外，在图8中还可以发现，在两组对称产生的马蹄涡结构之间通常存在着由单侧涡腿产生的较小马蹄涡结构(如F1到F4所标识的拟序结构，这种结构的涡头部分的形成如图5中所示)。此外，还可以看到在两组结构之间沿展向不均匀分布着流向涡结构，这些流向涡结构总是左右摆动，而不是单纯的沿流向向前发展。

图8 湍流场中存在的涡结构Fig.8 Vortex structures at different time instants

由上面的论述可知，流向涡结构在向下游发展过程中，既能单独演化，又有相互作用。单独演化时往往会出现单侧涡腿产生的较小马蹄涡，相互作用则可发现两个流向涡结构相互靠近合并成较大马蹄涡结构。同时这些马蹄涡形状结构并非完全对称的。

4 结论

本文采用动态亚格子模型和预处理的大涡模拟方法对槽道可压缩流湍流进行了精确的数值模拟，得到了满意的结果。通过对结果的分析比较，可以得出以下结论:

(1)文中LES计算得到的湍流平均速度剖面、湍流强度及平均表面摩擦系数，与精确的DNS解和经验公式符合得很好，验证了计算结果的可靠性;

(2)计算结果表明，湍流场内的流向涡结构与低速条带是相对应的，流向涡在x-z平面内会左右摆动、交替出现，因而快慢条纹会左右摆动、而不单纯是流向的;

(3)结果清楚地展示了湍流场中的流向涡演化成马蹄涡的过程，这些流向涡向下游的发展过程中，既可以相互作用，由对称涡腿产生马蹄涡结构，又能单独演化成单侧涡腿的马蹄涡结构。特别是，在亚谐波的作用下，湍流场中可能会出现的复杂马蹄涡结构，这些对称涡腿产生的马蹄涡相互之间是交错出现的，这些结构由于湍流扰动，会出现形状也不是严格对称的现象。

［1］RINGUETTE M J，WU M，MARTIN M P.Coherent structures in direct numerical simulation of turbulent boundary layers at Mach 3［J］.Journal of Fluid Mechanics，2008，594:59-69.

［2］连祺祥.湍流边界层拟序结构的实验研究［J］.力学进展，2006，36(3):373-388.

［3］张楠，陆利蓬，段真真，等.湍流边界层内准流向发卡涡生成的数值模拟［J］.应用数学和力学，2008，29(1):13-20.

［4］ZHOU J，ADRIAN R J，BALACHANDAR S，et al.Mechanisms for generating coherent packets of hairpin vortices in channel flow［J］.Journal of Fluid Mechanics，1999，387:353-396.

［5］陈林，唐登斌，刘小兵，等.边界层转捩过程中环状涡和尖峰结构的演化［J］.中国科学G辑，2009，39(10):1520-1526.

［6］ROBINSON S K.Coherent motions in the turbulent boundary layer［J］.Ann Rev Fluid Mech，1991，23:601-639.

［7］DELO C J，KELSO R M，SMITS A J.Three-dimensional structure of a low-Reynolds-number turbulent boundary layer［J］.Journal of Fluid Mechanics，2004，512:47-83.

［8］CANTWELL B J.Organized motions in turbulent flows［J］.Annual Review of Fluid Mechanics，1981，13:457-515.

［9］周恒，熊忠民.湍流边界层近壁区拟序结构起因的研究［J］.中国科学 A 辑，1994，24(9):941-948

［10］MOIN P.Advances in large eddy simulation methodology for complex flows［J］.International Journal of Heat and Fluid Flow，2002，23:710-720

［11］SANDHAM N D，YAO Y F，LAWAL A A.Large-eddy simulation of transonic turbulent flow over a bump［J］.Inernational Journal of Heat and Fluid Flow，2003，24:584-595.

［12］周恒，平板湍流边界层底层的不稳定波［J］.力学学报，1988，20(6):481-488.

［13］BRILEY W R，TALTLOR L K，WHITFIELD D L.Highresolution viscous flow simulations at arbitrary Mach umber［J］.Journal of Computational Physics，2003，184(1):79-105.

［14］史万里，葛宁，盛春华.计算全马赫数下粘性流动的隐式算法［J］.空气动力学学报，2009，27(5):566-571.

［15］GERMANO M，PIOMELLI U，MOIN P，et al.A dynamic subgrid-scale eddy viscosity model［J］.Physics of Fluids，1991，3:1760-1765.

［16］SPALART P R.Direct simulation of a turbulent boundary layer up to Reθ=1410［J］.Journal of Fluid Mechanics，1988，187:61-98.

［17］TYAGI M，ACHARYA S.Large eddy simulation of film cooling flow from an inclined cylindrical jet［J］.ASME，2003，125:734-742.

［18］SCHOPPA W，HUSSAIN F.Numerical study of near-wall coherent structures and their control in turbulent boundary layers［A］.6th International Conference on Numerical Methods in Fluid Dynamics，Volume 515［C］.ISBN 978-3-540-65153-6.Springer-Verlag，1998.