由近壁条带失稳诱导的发卡涡包的形成机制

李 健，董 刚，张 辉，陈正寿，张建雷

（1.浙江海洋大学 船舶与机电工程学院，浙江 舟山316022；2.南京理工大学 瞬态物理国家重点实验室，南京210094；3.浙江省近海海洋工程技术重点实验室，浙江 舟山316022）

0 引 言

发卡涡结构是Theodorsen[1]于1952年首次提出的概念模型，它是指以“Λ”形状或“Ω”形状呈现出的旋涡结构。随着实验测试技术和计算机性能的高速发展，转捩和湍流边界层中存在发卡涡结构这一现象逐渐得到了证实[2-7]，尤其是直接数值模拟的开展，提供了三维瞬时流场的丰富数据并给出了发卡涡存在的直接证据[7]。该结构不仅可以被用来解释湍流雷诺应力和近壁条带的形成机理，而且对壁面摩擦阻力的形成和流动传热传质过程有重要的影响。

发卡涡的形成可来源于多种途径，如转捩过程中三维扰动波的发展[8-9]、发卡涡之间的碰撞[10]以及低速条带的失稳[11-15]。对于转捩后期出现的发卡涡结构，吴介之等[8]指出发卡涡的形成起源于层流边界层中的展向涡层，当扰动波指数型增长到一定阈值时，展向涡层通过Kelvin-Helmholtz不稳定卷曲形成展向涡管。在非线性扰动背景下，展向涡管发生变形，由于远离壁面的涡管沿下游的运动速度较快，展向涡管沿流向拉伸最终形成发卡涡。Adrain等[10]通过数值模拟发现发卡涡之间的相互作用可以诱导出新的发卡涡。由于发卡涡在沿下游的发展过程中其尺度逐渐增大，这会导致多个发卡涡之间沿展向相互碰撞，涡腿与涡腿、涡脖与涡腿、涡脖与涡脖之间的碰撞都可合并形成尺度更大的发卡涡。发卡涡也可由低速条带的失稳诱导产生。Asai等[11]通过平板实验发现肿胀模式（varicose mode）低速条带失稳可诱导产生发卡涡，该实验结果随后被Brandt[12]用数值模拟证实。Brandt和De Lange[13]在无背景噪音的环境下模拟了高、低速条带沿流向相互碰撞的演化过程，并在条带对称碰撞的算例中发现了发卡涡结构。在该算例中，低速条带先是发展成肿胀模式低速条带，随后失稳形成发卡涡结构。除肿胀模式低速条带失稳可诱导出发卡涡结构外，Konishi和Asai[14]通过平板实验发现亚谐型弯曲模式（subharmonic sinuous mode）低速条带失稳同样可以诱导产生发卡涡。本文作者曾对亚谐型弯曲模式低速条带进行了直接数值模拟[15]，验证了Konishi和Asai的实验结果[14]并总结了发卡涡形成的三个阶段：（1）低速条带之间沿展向发生碰撞导致低速条带之间的展向截面中同时出现高、低速流体；（2）高、低速流体之间沿流向剪切形成法向涡量和负展向涡量；（3）法向涡量和负展向涡量与位于其下部上游的流向涡量完成搭接形成发卡涡结构。

在湍流边界层中，发卡涡可以以单个涡结构的形式出现，但大多以发卡涡包结构（hairpin vortex packet，一种由多个发卡涡组成的发卡涡链）的形式出现[3，16]。Zhou等[16]模拟了由单个发卡涡演化成发卡涡包的完整过程，发现在初始发卡涡（PHV）的发展过程中，PHV腿部相互靠近的部位与其上方的高速流体之间相互作用可产生弓形涡头，弓形涡头形成后与PHV的腿部搭接形成二次发卡涡（SHV）。此外，他们还发现初始发卡涡的下游也可诱导出新的发卡涡，称为下游发卡涡（DHV），该涡的形成与初始发卡涡头部下游侧的突起有关。Adrain[3]认为壁湍流实际是由不同尺度的涡包结构构成且每个涡包通常包含5～10个发卡涡，其尺度和迁移速度均随壁面距离增大。成璐和姜楠[17]运用高时间分辨率粒子图像测速技术，识别和提取了湍流边界层中发卡涡包结构的空间特征，发现在湍流边界层中不同法向位置处均存在发卡涡包结构，并指出近壁和外区的发卡涡包结构不是孤立存在的，而是和外区发卡涡包及其所夹带的高、低速条带流体构成了紧密联系的湍流边界层整体。本文作者在之前亚谐型弯曲模式低速条带失稳的研究中观察到了由PHV到SHV的演化过程[15]，但由于所给的初始扰动振幅较小，低速条带的碰撞强度较弱，SHV的形成总是发生在PHV的耗散后，故不能明显观察到发卡涡包结构。

利用发卡涡包模型可以很好地解释近壁条带的形成机制，即近壁条带是发卡涡包迁移后遗留下的“痕迹”[3]，而近壁条带自身的失稳能否反之诱导出发卡涡包结构目前还未得到证实。基于此，本文以仅包含低速条带而没有涡结构的小尺寸区域为初始条件，通过给定的扰动方程激发条带失稳，模拟了三维流场的演化过程并揭示了近壁小尺度发卡涡包的形成机理。

1 数值方法与初始流场的构造

以槽道流形成的不可压缩湍流边界层为基本研究对象，控制方程如下：

式中：u为速度矢量，p为压强，ν为运动粘度。利用傅里叶-切比雪夫谱方法（Fourier-Chebyshev spectral method）[18]离散上述方程。由于湍流脉动在流向（x）和展向（z）具有统计平均特性，故采用傅里叶-伽辽金方法变换并使用周期性边界条件；而法向（y）则采用切比雪夫-τ方法变换并使用无滑移壁面条件。方程（1）中的时间项采用三阶精度的半隐式向后差分格式处理；方程（1）左边第二项，即非线性项，采用3/2规则以消除混淆误差；方程（1）右边的压力项和线性项采用切比雪夫-τ方法和影响矩阵法联立求解，以消除流场出现的残余散度。与其它基于亚网格模型的大涡模拟方法（Large Eddy Simulation）[19-20]以及基于各种湍流模型的平均N-S方程的方法（Reynolds Average Navier-Stokes Equations）[21]模拟湍流不同，采用基于谱方法的DNS没有引入任何模型假设，而是在谱空间中直接对控制方程进行求解，因而具有很高的计算精度。本文作者曾多次使用上述方法进行壁湍流减阻[22]以及湍流稳定性方面的研究[15，23-24]，其有效性和可靠性已得到验证。

湍流边界层中条带的数学形式采用文献[25]提供的壁面单位下的初始条带分布：

为使两根低速条带之间沿展向发生碰撞，采用亚谐型弯曲模式的初始扰动形式[15]激发条带失稳，其表达式如下：

2 结果与讨论

2.1 “X”形、“Λ”形和“Ω”形涡的形成

图1 低速条带和涡结构初期的演化过程Fig.1 Evolutions of the low-speed streaks and the vortex structures during the early stage

图2 “Ω”形涡结构的涡量分布（对应图1 （c））Fig.2 Distribution of the vorticity on the‘Ω’-like vortex structure(corresponding to Fig.1(c))

图3 z+=100截面中低速流体（u′+＜0）和负展向涡量（＜0）的统计平均值历史曲线（u′+统计范围：y+=100～165，统计范围：y+=0～165）Fig.3 Time history of the statistical average value of the low-speed fluid(u′+＜0)and the negative spanwise vorticity(＜0)in y-z plane at z+=100(The statistical region of u′+：y+=100～165,the statistical region of：y+=0～165)

与文献[15]中Am=0.02条件下拟序结构初期的演化相比，Am=0.1条件下初始发卡涡的形成过程在形态上与Am=0.02条件下的一致，由亚谐型弯曲模式低速条带诱导的三维涡结构都是从“X”形涡演化成“Λ”形涡和“Ω”形涡，且涡结构在形态上的改变都与低速条带的弯曲、抬升以及条带之间的碰撞、反弹密切相关。然而，不同初始扰动振幅条件下的涡量强度不同。初始扰动振幅的增大会加剧条带的振荡幅度，进而使流向涡的强度增强。进一步地，流向涡的增强可加剧低速条带之间的碰撞，进而使环形剪切层的剪切强度增强。因此，初始扰动振幅的增大可导致流向涡量、法向涡量和展向涡量同时增强。

2.2 发卡涡包的形成

图4 t+=109时刻低速条带和涡结构的等值面分布（等值面同图1 ）Fig.4 Distribution of the iso-surfaces of the low-speed streaks and the vortex structures at t+=109(The iso-surfaces are the same as described in Fig.1)

图5 t+=109时刻流场的展向截面分布（x+=120）Fig.5 Distribution of the flow field in y-z plane at t+=109(x+=120)

图6 t+=127时刻低速条带和涡结构的等值面分布（等值面同图1 ）Fig.6 Distribution of the iso-surfaces of the low-speed streaks and the vortex structures at t+=127(The iso-surfaces are the same as described in Fig.1)

图7 t+=127时刻流场的展向截面分布（等值线同图5 （a））Fig.7 Distribution of the flow field in y-z plane at t+=127(The iso-lines are the same as described in Fig.5(a))

随着时间的发展，PHV的头部包裹着椭圆形低速区域继续沿下游运动，位于“树杈”形低速流体分叉处的QHV逐渐增强（图6）。在此过程中，“树杈”形低速流体合并处逐渐扩大并沿下游拉伸，二次发卡涡（SHV）和三次发卡涡（THV）同时在“树杈”形低速流体拉伸处形成。图7（a）-（b）分别给出了对应于图6中THV和SHV头部位置流场的展向截面分布。在THV头部对应的展向截面中（x+=130，图7（a）），两根低速条带上部合并，低速条带的合并位置与其上方的高速流体剪切形成环形剪切层，该剪切层的出现导致相应部位环形涡量的形成。因此，环形涡与位于其下部的一对互为反向旋转的流向涡相连构成了图6中的THV。在SHV头部对应的展向截面中（x+=175，图7（b）），合并的低速区域发生反弹且已完全进入到高速区，反弹的低速流体与其周围的高速流体同样发生剪切形成环形涡。该环形涡与位于其下方互为反向旋转的流向涡共同构成了图6中的SHV。注意到此时在SHV所处的展向截面中（图7（b））反弹低速流体沿流向已完全进入到高速流体中，故高、低速流体剪切形成的环形涡量相比THV的窄。

图8 t+=109时刻和t+=127时刻流场的展向截面分布（z+=100）Fig.8 Distribution of the flow field in x-y planes at t+=109 and at t+=127(z+=100)

根据以上分析，图9总结了由亚谐型弯曲模式低速条带失稳诱导的发卡涡包的形成过程，其中，图9（a）对应图1（c）中的流场分布；图9（b）对应图4中的流场分布；图9（c）对应图6中的流场分布。在拟序结构发展的初期（图9（a）），低速条带间沿展向相互碰撞形成低速反弹区域，该反弹区域与周围的高速区域沿流向剪切形成环形涡，环形涡与位于其下部上游的流向涡相互搭接构成了严格对称的PHV。在流向涡的旋转作用下，低速条带发生破碎且破碎条带上部呈“树杈”形结构（图9（b）），该“树杈”形结构在分叉处上游附近与高速区域沿流向剪切形成非对称的QHV。随着“树杈”形低速区域的合并处逐渐沿流向拉伸，∂v/∂x和∂u/∂y在该拉伸区域相互叠加并在叠加位置处形成负展向涡量，而相应部位负展向涡量与PHV的腿部搭接最终形成SHV和THV。

图9 发卡涡包的形成示意图（图中阴影区域为下部的低速条带，绿色区域为上部的低速条带，橙色区域为涡结构）Fig.9 Schematic of the formation of the hairpin vortex packets(Shaded regions represent the lower part of the low-speed streaks,green color regions represent the upper part of the low-speed streak,orange color regions represent the vortex structures)

3 结 论

本文采用直接数值模拟方法，通过增大初始扰动振幅，模拟了小尺寸槽道湍流中亚谐型弯曲模式低速条带的失稳并分析了发卡涡包的形成过程。结果表明，发卡涡包的形成与低速条带间的碰撞和破碎有着紧密的联系。一方面，低速条带间沿展向的碰撞使条带间的高速区域中出现低速流体分布，高、低速流体沿流向剪切形成的法向涡和负展向涡与位于其下部的流向涡相互搭接构成了对称分布的发卡涡结构。另一方面，由于近壁流向涡的旋转作用，低速条带发生破碎且上部的形态沿流向呈“树杈”形结构。在“树杈”形低速流体的上游分叉处附近，流向涡将低速流体向上携带，因而使高、低速流体之间的梯度增强，形成非对称的发卡涡。同时，“树杈”形低速流体的下游部位随时间沿流向拉伸，并与其上部的高速流体剪切形成对称的发卡涡结构。最终，多个对称和非对称的发卡涡结构沿流向排列，共同构成了发卡涡包结构。