基于MSFLE与Liutex的涡合并实验研究

2020-08-08郭延昂董祥瑞蔡小舒

空气动力学学报 2020年3期

郭延昂, 董祥瑞, 蔡小舒,*, 周骛

(1. 上海理工大学颗粒与两相流测量研究所, 上海 200093; 2. 上海市动力工程多相流动与传热重点实验室, 上海 200093)

0 引言

涡是湍流边界层拟序结构的核心，影响条带、猝发、上升流、下扫流等其它拟序结构的发展和演化[1]，著名空气动力学家Küchemann曾说，“旋涡是流体运动的肌腱”[2]，这句话深刻概括了涡的重要性。自从1952年Theodorsen[3]提出马蹄涡(又称发卡涡)模型以来，无数科研工作者致力于研究涡的相关机理。精确可靠的实验测量是诸多机理研究的前提与保障。目前，国内外已发展了多种实验测量和流场可视化方法，如氢气泡法、烟流法、染料法等。Kline等人[4]采用氢气泡法发现了湍流边界层近壁流场中的快慢斑和长带条结构，利用染料发现了猝发现象。Brodkey等人[5]使用示踪粒子和立体摄影发现了上升流和呈手指状的下扫流，并认为下扫流是产生雷诺应力的主要来源。Lian[6]利用氢气泡法对不同单个流向涡进行了研究，结果显示流向涡的旋转运动会导致下扫流和上升流的产生，流向涡结构沿流向拉伸后会从周围吸收能量，导致旋转速度增大。虽然这些可视化方法是基于定性研究的，但因其经济且易于实现，所以成为了观察大尺度湍流结构很好的实验方法[7]，为人们认识和理解湍流边界层中拟序结构的宏观特性提供了很大的帮助。

随着光学技术及计算机技术的快速发展，粒子图像测速(particle image velocimetry, PIV)和粒子追踪测速(particle tracking velocimetry, PTV)等得到了广泛应用。Adrain等人[8-10]利用PIV测量到了边界层外区存在着沿流向排列且具有相同速度的发卡涡，它们组成了大尺度发卡涡包结构。近年来，三维(3D)测量技术得到了迅速发展。3D PIV或PTV测量的基本原理是结合多个2D PIV或PTV，采用多个相机进行同步测量，再利用算法重构得到3D速度场。如Kinzel等人[11]利用两个具有不同空间分辨率的PTV开展3D测量，低空间分辨率的PTV测量速度场及大尺度流动特征，高空间分辨率的PTV测量与大尺度流动特征相关的涡量和应力等小尺度物理量。Hutchins等人[12]使用两个成一定角度布置的相机开展体视PIV研究，测量与流向成45°和135°平面中的涡结构，其中在135°平面中测量到了一对反向旋转的涡结构，两涡中间发生上喷现象，这与发卡涡的涡腿特征一致。将PIV与高速相机和高能连续激光或高频脉冲激光结合，形成了时间分辨PIV(time-resolved PIV, TRPIV)，可以测量流场随时间的连续变化过程[13-14]。层析PIV技术可以得到瞬时三维三分量(3D3C)速度场，已经被应用于测量湍流边界层和深入分析拟序结构机理[15-17]。如Schroder等人[15]使用六个CMOS相机，采用层析TRPIV测量了与负雷诺应力事件(Q2和Q4)相关的拟序结构及其在湍流边界层对数律层中的时间演化过程，研究了高雷诺数湍流边界层中发夹涡的初始发展机理。Katz等人[18]详细介绍了数字全息3D流场测量的原理和应用。上述测量方法在定量测量3D流场以及研究湍流机理方面发挥了巨大作用，不过测量系统非常复杂，且存在测量分辨率和测量区域大小的矛盾，难以测量涡随时空的动态演变过程。对此，蔡小舒等人[19]提出了运动单帧长曝光(moving single-frame and long-exposure, MSFLE)图像法，测量系统简单方便，可以得到直观的低雷诺数湍流边界层中涡随时空的演变过程。

除实验测量外，测量结果中的涡提取手段也是湍流研究的重要议题。Epps[20]对涡识别方法进行了全面总结，其中基于欧拉体系的方法有λ2[21]、Q[22]、λci[10]、Ω[23]等，基于拉格朗日体系的方法有有限时间Lyapunov指数法[24]、MZ[25]等。上述涡识别方法大都需要人为设定阈值，不同的阈值可能会得到不同的涡结构，而且不能提供涡的准确定义。近几年，刘超群团队基于涡的数学定义提出了一种新的涡识别方法Liutex[26-30]，相比于其它涡识别方法，Liutex向量可以准确表征流体旋转强度及当地旋转轴方向。

综上所述，本文结合MSFLE图像测量法和Liutex涡识别方法对湍流边界层涡结构的时空演变过程开展了实验研究，得到了直观的涡结构特征及其演化过程，并对涡合并现象进行了量化分析。

1 实验方法及装置

1.1 MSFLE图像测量方法

MSFLE图像测量方法是从单帧长曝光(single-frame and long-exposure, SFLE)图像法[31]发展而来。在SFLE测量方法中，流场中布撒跟随性良好的示踪粒子，激光器连续照明流场，通过设置相机较长的合适的曝光时间，可以将粒子的运动轨迹清晰地记录在单帧图像中。

图像中运动轨迹的长度代表了曝光时间内示踪粒子的运动距离，进而可以得到流体运动速度，通过前后两帧图像可得到流体运动方向。示踪粒子的运动速度可以由下式得到：

(1)

式中：v为示踪粒子的运动速度，s表示图像中示踪粒子运动轨迹长度，m为相机镜头的放大倍率，Δt是相机曝光时间。

相比于PIV和PTV采用两次曝光和两帧图像来获得粒子运动的始末位置，SFLE法采用一次长曝光可以在单帧图像中非常直观地显示粒子运动轨迹。而且SFLE法对实验条件要求较低，仅需较小功率的连续激光器和工业相机，测量装置简单，实验成本低。

通常的流场测量方法是通过欧拉视角研究流场，如PIV等。在这些方法中，测量系统固定在某个位置，只能得到在该固定位置的流场信息，无法测量湍流边界层中涡结构随时空演变过程。另外图像法测量中相机视场和图像分辨率是一对矛盾，采用高分辨率的相机镜头测量涡精细结构时，因图像视场较小，无法捕捉快速运动的涡结构。针对上述问题，本文在SFLE的基础上，采用MSFLE图像测量方法，通过拉格朗日视角研究湍流边界层中涡结构特征及其演化过程。

MSFLE图像测量方法具有运动和长曝光两大特点。一方面，实验中相机测量系统匀速运动以跟踪捕捉到与其速度相同或相近的涡结构；另一方面，采用较长曝光时间将示踪粒子运动轨迹记录在单帧图像中，得到的连续测量图像可以清晰反映涡结构及周围流场的发展变化过程。

MSFLE是一种拉格朗日型测量方法，可得到涡结构随时空运动演变过程的直观图像，但得到的流体速度是相对于相机运动速度的相对速度。在该方法中高于相机移动速度的流体在图像中向下游方向运动，低于相机移动速度的流体向上游方向运动，与相机基本相同移动速度的流体结构则在图像中位置基本不变。这样可以清晰地显示出高低速流体之间的剪切作用。如图1所示，连续拍摄的图像清晰地显示出相对相机运动速度的高/低速流动间的剪切现象及涡是如何“搓”出来的。图1(a)为低速流体(黄色)与高速流体(红色)发生剪切的现象，剪切导致顺时针旋转涡的产生如图1(b)所示。图1直观地显示了涡与上下高/低速流体的联系，即涡是由高/低速流体的剪切作用产生的，这正如著名流体力学家陆士嘉提出的观点，“流体的本质就是涡，因为流体经不住搓，一搓就搓出了涡”，“搓”就是高/低速流体之间的剪切。

1.2 实验装置

实验测量系统如图2所示。实验中片光位于矩形管道展向中间位置，平行于来流并垂直矩形管道入射，以尽可能避免测量区受到管道两侧壁面边界层的影响。实验测量了矩形管道底部湍流边界层流向-法向二维平面内的涡结构演变过程。图2中x轴为流向方向，y轴为法向方向，z轴为展向方向。实验系统由恒位水箱、回水箱、水泵、流量计、流量调节阀、实验管道、相机镜头、激光器和电动导轨等构成。水泵将回水箱中的水输送至恒位水箱以保证管道入口处压力恒定，从而确保实验段来流的稳定，水流过实验管道之后经流量调节阀和流量计回到水箱，形成循环。

图2 实验测量系统示意图Fig.2 Schematic diagram of the experimental measurement system

有机玻璃实验管道长2200 mm，管道横截面为80 mm×80 mm的正方形。长1500 mm、宽78 mm、厚4 mm有机玻璃平板水平放置在实验管道底部，故流体流通横截面为80 mm×76mm的矩形，平板前缘距管道入口处700 mm，正对来流方向的平板前缘设计为半椭圆形，长短轴比例为4∶1，以减少平板前缘对流动的干扰。将直径为4 mm的拌线布置于距平板前缘50 mm处以加速边界层的转捩。示踪粒子选用平均直径5 μm、比重1.04和折射率为1.584的聚酰胺树脂颗粒(polyamide resin particle, PSP)。图像测量系统由分辨率1280 pixel×1024 pixel、像元大小4.8 μm的XIMEA相机，放大倍率为0.14倍的远心镜头和输出1 mm厚片光、450 nm波长、4 W功率的连续半导体激光器组成。图像视场范围为44 mm×35 mm，空间分辨率为34 μm。相机曝光时间为100 ms，帧率9.995 fps，帧间间隔50 μs，基本实现连续测量。

图像测量系统安装在沿流向平行布置的导轨上，伺服电机控制图像测量系统匀速运动以捕捉矩形管道底部边界层中与此运动速度相近的涡结构。当流量为1.6 m3/h时，采用SFLE图像法由公式(1)计算得到来流速度约为U∞=100 mm/s。令Uc表示图像测量系统运动速度，U*=Uc/U∞。通过大量实验发现，当U*∈[0.8,0.9]时测量到的边界层内涡结构现象最为丰富，这表明涡大都以这个速度运动，故图像测量系统速度在此区间内选取。令x表示距平板前缘的流向位置，实验测量范围从x=400 mm到x=1300 mm，两个位置对应的边界层厚度分别为10 mm和20 mm，动量厚度θ分别为0.97 mm和1.94 mm，基于动量损失厚度的雷诺数分别为Reθ=θU∞/ν=97和194，其中ν为流体运动黏性系数。

需要指出的是本研究中得到的是3D涡结构在流向-法向2D截面的涡特性，反映的是涡在该截面的时空演变过程。该2D时空演变过程与3D涡结构时空演变过程密切相关。

2 基于骨架提取的图像处理方法

从图像获取流场速度信息的图像处理方法主要包括粒子轨迹识别和速度计算两部分[32]。粒子轨迹识别图像处理流程如图3所示，由去噪锐化、自适应阈值分割(OTSU)、去除小颗粒和骨架提取四部分组成。自适应阈值分割算法是一种多阈值自动分割方法，通过分析像素周围局部范围的灰度特性来决定该像素点的阈值，不需要人为设定阈值，这样能够使阈值分割后的图片保留更多信息。对图像进行骨架提取后粒子轨迹是单像素宽度的轨迹，即像素精度的迹线。可以更准确计算轨迹长度，且对于弯曲的轨迹也能区分各像素点处速度方向的细微区别，同时保持图像的形态学特征。

图3 轨迹识别处理流程图Fig.3 Process flow chart of trajectory identification

获得像素精度的迹线后，通过连通区域分割函数识别出每条轨迹上的所有像素点，传统算法采用迹线起点与终点之间的距离作为长度，在处理弯曲轨迹时会产生较大误差。本文采用求线积分的长度计算算法[32]，先求出每两个像素点之间的距离，然后迭加得到总的轨迹长度，两种算法的对比如图4所示。得到轨迹长度后由公式(1)计算流场速度，速度方向由连续的两帧图像得到。

图4 不同算法求轨迹长度的对比[32]Fig.4 Comparison of trajectory length between two algorithms[32]

3 Liutex涡识别方法

通过骨架提取算法计算得到速度场后，要进行涡结构的识别。由于涡量在近壁强剪切区域值很大，而在涡旋转区域较小，但不能够表征流体旋转，本文采用Liutex理论[26-30]进行涡的识别和表征。Liutex理论可以系统地描述当地流体的转动强度大小及旋转轴方向，其基本思想为：如果速度梯度张量V存在一对复特征值λcr±λci和一个实特征值λr，那么当地流体旋转轴方向即为V的实特征向量的方向，即Liutex的矢量方向。Liutex的大小与方向分别以R和r表示：

V·r=λrr

(2)

R的显示计算公式如下：

(3)

因此，Liutex向量R被定义为R=Rr。

综上，本文所采用的涡识别与表征方法是通过骨架提取图像处理方法得到速度场，进一步使用三角剖分法进行插值，插值大小为图像像素数，之后插入第三维数据，令z=1，w=0，再计算速度梯度张量，最终通过公式(3)计算得到Liutex以实现涡强度的表征。

4 MSFLE测量结果与分析

4.1 湍流边界层实验验证

为验证矩形管道底部边界层为湍流边界层，利用SFLE方法分别测量了x=400 mm与x=1300 mm处边界层沿法向平均速度分布。测量时片光位于展向中心处，垂直矩形管道平行来流方向布置，相机垂直片光进行测量。对测得的瞬时速度场进行统计分析，再根据公式y+=yuτ/v，u+=u/uτ将平均速度u和法向高度y用摩擦速度uτ无量纲化。摩擦速度uτ通过壁面切应力τw获得。τw根据黏性底层的速度梯度获得。

(4)

式中ρ为流体密度，μ为动力黏性系数。最终计算得到x=400 mm处uτ=6 mm/s，x=1300 mm处uτ=5.48 mm/s。两处边界层法向平均速度分布结果见图5，该无量纲速度分布与湍流边界层内层[33]速度分布公式吻合，且两个流向位置的平均速度分布呈现自相似，表明在矩形管道底部的实验测量区域内为湍流边界层。

图5 x=400 mm与1300 mm处边界层法向平均速度分布Fig.5 Dimensionless velocity profile of experimental results at x=400 mm and x=1300 mm

4.2 涡合并结果与分析

采用MSFLE测量，获得了在矩形管道底部湍流边界层内涡结构合并的整个过程。图6给出了其中获得的两对涡合并过程中的部分时刻的图像，以及计算处理结果。其中左图为原始图像，中图为图像处理后得到的u*云图及流线图(u*=u/U∞，u为相对流向速度)，右图为计算得到的R云图及流线图。图中横坐标和纵坐标均用δ0无量纲化，δ0为实验初始位置(即x=400 mm)处的边界层厚度10 mm，图像下边缘对应矩形管道底层上表面，白色轨迹线是示踪粒子在曝光时间下的运动轨迹。流体从图像左侧流向图像右侧，由于测量系统以速度Uc匀速运动，当在连续测量的图像中示踪粒子轨迹向右侧运动表示该处流体流动速度大于相机运动速度Uc，如u*云图中上部区域流线所示；当流体速度小于Uc时，示踪粒子轨迹向左侧移动，如u*云图中下部区域流线所示。定义相对于Uc运动快慢的流体分别为高速流体和低速流体。当示踪粒子轨迹在图像中仅为亮点时，则表明该粒子为稳定匀速运动状态，速度大小为Uc。

图6为流向位置接近实验终止位置、U*=0.9时测量到的位于边界层内的涡合并过程中的4帧典型实验图像及处理结果。图6的u*云图和R云图显示的涡结构的法向高度y/δ0在1～1.5之间，这里采用的δ0=10 mm，实验终止位置处的边界层厚度为20 mm。该系列图像捕捉到了三个顺时针旋转的A涡、B涡和C涡。A涡和B涡在流动过程中发生合并形成尺寸更大的D涡，而C涡在此期间位置及形状基本保持不变，表明运动速度与Uc基本一致。

令t*=tU∞/δ0,t为测量时刻。图6(a)为t*=80时，两个尺寸基本相同的A涡和B涡开始发生合并，从u*云图及流线图可知此时A涡和B涡涡心之间的距离约为0.65。t*=81时，如图6(b)所示，两涡涡心之间的距离减小到0.55。从测量可知B涡运动速度逐渐变慢，A涡运动速度基本不变，逐渐追上B涡，两涡水平靠近。到图6(c)t*=83时，A涡和B涡开始接触且互相绕对方沿顺时针方向发生了旋转，此时两涡不仅具有流向方向的速度差，也具有法向方向的速度差。当t*=85时，如图6(d)所示，新的顺时针旋转的D涡形成。

(a) t*=80

图7 A涡和B涡合并过程中RInt随时间的变化Fig.7 Temporal distributions of RInt during the merging process of A vortex and B vortex

图8为U*=0.8时测量到的另一组涡合并过程中的4帧典型的实验图像及处理结果。t*=34时，由图8(a)R云图得到三个涡的ΔS分别为0.075、0.08和0.12，E涡和F涡开始发生合并，但F涡没有同G涡合并。分析认为由于F涡和E涡的尺寸更加接近，而F涡尺寸小于G涡，且G涡诱使下方流体向上运动到了F涡和G涡之间，阻碍了F涡和G涡的接触。t*=34时由图8(a)u*云图可知E涡和F涡两涡心之间的距离为0.9，t*=35时涡心距离减小到了0.7，由连续两张u*云图可知E涡运动速度变快，F涡速度基本不变，两个涡在流向方向具有速度差，因此E涡逐渐平移靠近F涡。到t*=36时，由图8(c)可知E涡和F涡发生接触互相绕对方沿顺时针方向发生旋转，表明两涡此时也具有了法向方向的速度差。当t*=38时，如图8(d)所示，新的顺时针旋转的G涡生成，其ΔS为0.18，即新生成的G涡尺寸约为合并初始时E涡和F涡尺寸之和。

E涡和F涡在合并过程中RInt随时间的变化如图9所示。合并过程可以分为三个阶段。第一个阶段为从t*=34到t*=35，由于流向速度差，在两个涡水平靠近过程中E涡不断吸收F涡能量，从而E涡RInt逐渐增大，F涡RInt则逐渐减小，此阶段称为平移靠近阶段。从t*=35到t*=37为第二阶段，两涡相互接触并在水平和法向方向速度差共同作用下绕对方沿顺时针旋转，在旋转过程中两个涡的能量发生扩散，能量高的E涡会将能量传输给F涡，因此E涡RInt逐渐降低，F涡RInt逐渐上升，此阶段称为相互旋转阶段。从t*=37到t*=38为第三阶段，RInt值为0.42的H涡生成，初始合并时E涡和F涡的RInt值分别为0.22和0.18，可见H涡强度基本为初始合并时两涡强度之和，此阶段为最终合并阶段。