Particle Image Techniques for Measuring Instantaneous Bed Shear Stress in Open Channel Flows
-
摘要: 床面切应力瞬时值的实验测量极为重要但富有挑战。为了利用粒子图像技术准确测量瞬时床面切应力，采用模拟不同图像和流动参数的合成粒子图片对粒子图像测速技术（PIV）、粒子图像示踪技术（PTV）和直线相关PIV技术（LCPIV）的测量精度进行分析，发现各技术在床面附近的测量精度均因床面影响而变差，其中，PIV技术受影响的范围距床面不小于判读窗口高度，PTV及LCPIV技术受影响的范围等于粒子图像半径；在床面影响范围外，PTV技术可实现无偏差测量，并在粒子图像密度较小且粒径适中时拥有高于PIV和LCPIV的精度。比较而言，PTV技术更适于测量瞬时床面切应力。将PTV技术用于雷诺数分别为8 007和11 462的明渠均匀紊流瞬时床面切应力试验测量，进一步验证了该技术无偏差、高精度测量的能力；发现了明渠紊流瞬时床面切应力不符合高斯分布，离散系数约为0.4，偏度系数和峰度系数随雷诺数增大而逐渐变大。研究结论为明渠紊流瞬时床面切应力的实验测量提供理论依据和实践参考。Abstract: Measurement of instantaneous bed shear stress in open channel flows is challenging but extremely important. In order to accurately measure the instantaneous bed shear stress with particle image techniques, synthesized particle images with different image and flow parameters were applied to evaluate the measurement precision of particle image velocimetry (PIV), particle tracking velocimetry (PTV), and line correlation PIV (LCPIV). The measurement precision of all techniques in the near bed region were found to be degraded towards the bed. The wall-normal height of the bed influence equals the particle image radius for PTV and LCPIV, and is not shorter than the height of the interrogation window for PIV. Outside of the bed influenced region, the measurement of PTV is unbiased and is with the highest precision when the particle image density is relatively small and the particle image diameter is medium. Therefore, PTV is more suitable for the measurement of instantaneous bed shear stress. The ability of PTV for unbiased and high precision measurement of instantaneous bed shear stress was further validated with experiments in two open channel flows with a Reynolds number of 8007 and 11462, respectively. The experiment results also indicated that the distribution of instantaneous bed shear stress in the turbulent open channel flows, with a coefficient of variation of about 0.4, is non-Gaussian. Meanwhile, both the skewness and kurtosis increase with the Reynolds number of the flow. The results provide theoretical foundation and practical reference for measuring instantaneous bed shear stress in open channel flows.
-
明渠床面切应力是单位面积河床上作用的水流剪切力，宏观上度量河床对水流施加的摩擦阻力，微观上以摩阻流速
$u_{*}=\sqrt{\tau_{{\rm{o}}} / \rho}$ （τo为床面切应力，ρ为水的密度）为载体度量近壁区紊流脉动尺度，也是描述水流与床沙、水工建筑物间相互作用的基本力学变量。因此，开展明渠床面切应力的试验测量对水力学及河流动力学研究具有重要价值。目前，试验测量床面切应力的方法主要有直接和间接测量两大类[1]。常见的直接测量方法包括测力浮板（FE）、油膜干涉法、微柱传感器（MPS）等，其中，基于微电子机械系统的FE和MPS近几年在多种壁面流动中被用于瞬时床面切应力测量[2-4]，但这些微型设备因易损性等原因不适于明渠水流。床面切应力的间接测量方法种类较多，对于明渠均匀紊流等简单流动，常用方法包括阻力平衡法、对数区时均流速对数律拟合法和外区雷诺应力线性拟合法[5-7]，但这些方法仅能实现时均床面切应力的估算。最近，传感器技术的快速发展使复杂流动瞬时床面切应力的测量成为可能，例如：徐华[8]将热膜传感器用于波浪作用下床面瞬时切应力的直接测量；Bagherimiyab和Lemmin[9]利用热膜传感器测得了明渠卵石床面的瞬时切应力；Chouila–Houri等[10]利用改进的热线传感器对无限长肋条尾流的瞬时床面切应力进行试验研究，但这些技术的测量结果受近床面热传导效应影响而存在不同程度偏差[1]。
近年来，粒子图像测量技术在床面切应力测量领域得到了应用。Kähler[11]、钟强[12]及de Silva[13-14]等先后采用显微粒子图像测速（PIV）或粒子示踪测速（PTV）方法实现了黏性底层内流速分布和时均床面切应力测量；Nguyen等[15]提出了用于床面速度梯度及时均床面切应力测量的界面PIV技术，这些成果展示了粒子图像测量技术在床面切应力测量领域的巨大潜力。但是，利用粒子图像测量技术测量瞬时床面切应力的研究相对较少。Sheng等[16]较早公开了利用数字全息显微成像和PTV技术实现紊流边界层瞬时床面切应力测量的报道；最近，Willert等[17-18]基于高放大倍率粒子成像和Nguyen等[15]提出的直线相关（line correlation）PIV技术在零压梯度边界层中实现了瞬时床面切应力测量；Wang等[19]采用高放大倍率粒子成像和大宽高比判读窗口PIV技术在水洞和水槽中实现了边界层瞬时床面切应力的测量，并分析了床面切应力的离散系数、偏度系数、迁移速度和尺度特征。尽管已有多种粒子图像测量技术被用于瞬时床面切应力的测量，但目前还缺少对这些技术适用性和测量精度的系统分析。
本文以明渠紊流瞬时床面切应力的准确测量为目标，针对已有研究中采用的PIV、PTV和直线相关PIV（以下简称LCPIV）技术，采用合成多参数粒子图像的方法系统研究各技术的测量精度和适用条件，并在此基础上选择最优技术开展明渠紊流床面切应力测量试验，初步验证方法的可靠性和适用性。
1. 床面切应力的粒子图像测量技术
1.1 测量原理
根据广义牛顿内摩擦定律，明渠均匀流中床面切应力的定义为：
$$ {\tau _{\rm{o}}} = \mu {\left. {\frac{{\partial u}}{{\partial y}}} \right|_{y = 0}} $$ (1) 式中，u为纵向瞬时流速，y为距床面的垂向高度。根据明渠均匀紊流时均流速剖面的壁面律，假定纵向瞬时流速在黏性底层（y+≤5）内同样符合线性分布规律，则床面瞬时切应力的测量可转化为如下黏性底层内瞬时流速的测量：
$$ {\tau _{\rm{o}}} \cong \mu {\left. {\frac{u}{y}} \right|_{{y^{\rm{ + }}} \le 5}} $$ (2) 在应用粒子图像测量技术时，由于实际明渠流黏性底层的厚度通常小于1 mm，拍摄床面附近粒子图像时需采用相对较大的成像倍率，流速判读算法应具有较高的空间分辨率；同时，由于每个纵向位置的流速矢量只能计算一个测点的床面切应力，为实现床面瞬时切应力的多点测量，还要求流速矢量具有尽量大的密度，即实现瞬时流场的测量。针对上述要求，以下对本文使用的PIV、PTV及LCPIV技术的主要技术特征进行介绍。
1.2 近床面PIV技术
利用自主研发的JFMeter® PIV软件[20]进行流场PIV计算，该软件采用基于多级网格迭代和图像变形技术的判读窗口互相关算法，并利用三点高斯插值计算亚像素位移。基于窗口相关算法的PIV流速测点始终位于判读窗口中心，在传统多级网格迭代技术中，随着窗口尺寸逐级减小，为确保窗口始终位于粒子图像（流动）内部，测点距床面的最小距离必须为最大一级判读窗口高度的1/2。
为利用PIV技术测量黏性底层内流速，使用宽高比较大的矩形判读窗口以提高垂直于床面方向的空间分辨率[19]。同时，为缩小测点与床面的最小距离，引入图1所示的虚拟窗口技术，该技术令测点与床面的最小距离为最小一级判读窗口高度的1/2。在对靠近床面的测点进行图像判读时，允许判读窗口延伸到图片边界之外，并在边界外的窗口内随机生成密度与实测图像一致且位移为零的虚假粒子图像，避免流速测量结果因窗口内图像分布不对称而产生显著偏差。后文分析时，判读窗口尺寸为32像素×16像素，窗口垂向重叠率为75%，测点距床面最小高度为8像素，测点垂向间距为4像素。
图 1 多重网格迭代及虚拟窗口技术示意图Fig. 1 Schematic diagram of multigrid iterative and fictitious window techniques
下载:
全尺寸图片
1.3 基于PIV的PTV技术
粒子示踪测速（PTV）本质上是对各粒子图像的形心进行跟踪，因此，PTV技术的首要环节是粒子形心的数字图像识别。目前，识别粒子图像形心的方法主要有灰度二值化法、灰度分布拟合法及粒子掩膜相关法[21]。本文采用粒子掩膜相关法，定义粒子掩膜的大小等于实际拍摄的示踪粒子图像的平均直径，图像亮度满足2维高斯分布。将粒子掩膜与粒子图片进行卷积运算，则得到的相关系数矩阵中各局部峰值点靠近粒子形心，形心坐标可通过对局部峰值进行三点高斯插值计算。
为了在黏性底层内测得较大密度的流速矢量，采用Ohmi和Li[22]提出的松弛迭代算法进行粒子形心跟踪，其核心思想为相邻粒子的运动相似。松弛迭代算法在实施过程中需要定义3个关键阈值：确定第一帧图像中待跟踪粒子的相邻粒子的邻域半径Rn；确定待跟踪粒子在第二帧图像中的候选匹配粒子的最大位移Rs；判断相邻粒子运动是否相似的容许位移偏差Rc。
邻域半径、最大位移和容许偏差的取值与粒子图像及局部流动特性有关，且大小对计算精度及效率影响较大。经典松弛迭代算法通常根据待测流场的平均速度估算上述参数，且计算域内使用相同取值。本文参考Kim和Lee[23]的研究，采用基于PIV的匹配参数确定方法，具体为：利用PIV算法根据两帧粒子图像计算2维流场，根据2维流场插值得到各示踪粒子形心处的位移估值，再令各粒子的邻域半径Rn、最大位移Rs及容许偏差Rc分别为对应位移估值的数倍（建议取值见表1）。上述方法充分考虑局部流动特征，可提高PTV计算精度和效率。
表 1 基于PIV的松弛迭代参数建议取值Table 1 Suggested values for parameters of iterative relaxation method based on PIV预估位移 邻域半径Rn 最大位移Rs 容许偏差Rc Vi >1.5Vi >1.0Vi <1.5Vi 1.4 LCPIV技术
经典PIV技术使用矩形判读窗口分割粒子图像，流速矢量的空间分辨率取决于判读窗口尺寸[24]。为了提高床面附近流速矢量的垂向空间分辨率，Nguyen等[15]提出了LCPIV技术。该技术以粒子图片中高度为1像素的判读直线代替矩形判读窗口，通过对两帧图片中同一位置的判读直线进行1维互相关运算得到水平流速分量，从而将测点垂向分辨率提高至1像素。但是，为了保证粒子图像的信噪比，判读直线的长度通常远大于常规判读窗口的宽度，使LCPIV的水平空间分辨率低于PIV技术，本文分析使用的直线宽度为64像素。
2. 基于合成粒子图片的技术对比
2.1 数据及分析方法
为了综合对比PIV、PTV及LCPIV技术测量黏性底层内瞬时流场的精度，采用合成粒子图片的方法模拟不同参数条件下近床面流场测量试验。其中，考虑试验中的主要变量和各测量技术的主要影响因素，选取粒子图像密度和直径作为图像参数进行分析；粒子图像密度以每像素的粒子数Nppp表示，取值分别为0.01、0.02、0.03、0.04，粒子图像直径dτ以像素表示，取值分别为3、5、7、9。根据文献[25]研究，PIV技术的适用范围为Nppp>0.02。为模拟黏性底层内的强剪切流动特征，合成图片对应的流场为：
$$ \left\{ \begin{gathered} u = ky, \hfill \\ v = 0 \hfill \\ \end{gathered} \right. $$ (3) 式中：u、v分别为流体在相邻两帧图片之间的纵向和垂向位移，单位为像素；y为距离床面的高度，单位为像素；k为速度梯度。本文选取速度梯度为流动参数进行分析，并根据常见试验条件分别考虑k=0.1和k=0.2两种工况。
根据上述图像和流动参数，本文共设置了32组合成算例。每种算例的图片尺寸均为512×128像素，图片底边与床面重合；利用标准粒子图像合成算法生成指定粒子图像密度和直径的第一帧图片[26]，再根据速度梯度取值得到各粒子位移，通过粒子图像平移得到第二帧图片。为进行统计分析，每组算例合成1 000对粒子图片，采用各测量技术进行计算均可测得1 000帧独立流场。根据流场实测值与理论值计算每个实测流速矢量的测量误差，分析误差均值（以下称平均偏差，表示为Eb）与方差（以下称均方根误差，表示为Erms）沿垂向分布规律，对比各技术在不同参数条件下的测量精度。
2.2 粒子图像密度对测量结果的影响
图2为PIV技术测量结果在不同粒子图像密度条件下的平均偏差及均方根误差沿垂向分布，由于其他粒子图像直径和速度梯度条件下的结果相似，本小节所有分析均以k=0.1，dτ=3测次结果为例进行说明。
图 2 不同粒子图像密度条件下PIV技术测量误差（k=0.1, dτ=3）Fig. 2 Measurement errors of PIV technique at various particle image density conditions (k=0.1, dτ=3)下载:
全尺寸图片
图2（a）中，平均偏差在远离床面区域沿垂向以20为周期正负交替，但基本不受粒子图像密度影响，各测次最大偏差值约±0.02像素。周期性偏差主要源于锁像素（pixel locking）现象，表现为测量值偏向整像素，与成像阶段粒子图像的像素化及图像判读阶段互相关函数的亚像素插值有关[27]。在距离床面约16像素（判读窗口高度）范围内，测量结果受床面影响而偏小，且距床面越近偏差越明显，最大偏差约–0.08像素。负偏差的产生与虚拟窗口技术有关，该技术在判读窗口内加入速度为零的虚拟粒子，使窗口内粒子的平均位移小于实际位移；需指出，若床面附近不使用虚拟窗口技术，PIV测量结果将显著偏大。图2（b）中，各测次均方根误差在远离床面的区域沿垂向无变化，除Nppp=0.01组因不适于PIV技术而误差稍大外，其余测次的均方根误差约0.04像素。在靠近床面的区域，均方根误差受床面影响而显著增大，最大误差0.14像素，受床面影响高度约16像素。
图3为PTV技术在不同粒子图像密度条件下的平均偏差及均方根误差。
图 3 不同粒子图像密度条件下PTV技术测量误差（k=0.1, dτ=3）Fig. 3 Measurement errors of PTV technique at various particle image density conditions (k=0.1, dτ=3)下载:
全尺寸图片
图3（a）中，除距床面2像素的测点存在约0.02像素的偏差外，其他测点的测量结果与理论值无偏差。图3（b）中，各测次均方根误差以10像素为周期沿垂向变化，且误差随粒子图像密度增加而变大，Nppp=0.01和Nppp=0.04时最大均方根误差分别约0.04和0.08像素；根据图3（b）中曲线可知，均方根误差在距床面2像素以内因床面影响而略微增大。分析PTV技术算法构成可知，其测量误差主要源于粒子中心识别算法，Lei等[28]的研究结果表明，粒子中心识别算法的误差随粒子图像之间重叠率的增加而变大，而粒子图像密度增大是提高粒子间重叠率的主要原因之一，从而合理解释了PTV均方根误差随粒子图像密度的变化规律。
图4为LCPIV技术在不同粒子图像密度条件下的平均偏差及均方根误差。根据误差分布剖面的特点可知，LCPIV测量结果受床面影响的高度与PTV技术基本一致。但是，图4（a）中，LCPIV测量结果沿垂向存在较明显的周期性偏差，且偏差幅值随粒子图像密度的变大而微弱增大，Nppp=0.01测次的最大偏差约±0.05像素。由于LCPIV技术与PIV技术均使用互相关算法，其平均偏差产生的原因也与锁像素现象有关，Wang等[19]近期利用LCPIV技术测得的边界层近壁面时均流速剖面也出现了周期性偏差。图4（b）中，各测次均方根误差同样沿垂向周期性变化，在相同垂向位置，除Nppp=0.01组次因图像密度过小而误差稍大外，其余测次的均方根误差最大值约为0.07像素。
图 4 不同粒子图像密度条件下LCPIV技术测量误差（k=0.1,dτ=3）Fig. 4 Measurement errors of LCPIV technique at various particle image density conditions (k=0.1, dτ=3)下载:
全尺寸图片
2.3 粒子图像直径对测量结果的影响
为分析粒子图像直径对各测量技术的影响，本小节选择Nppp=0.02，k=0.1测次为例进行分析，在其他粒子图像密度和速度梯度条件下，测量误差随垂向位置和粒子图像直径的变化规律相同，仅误差值略有差异，故在此未进行展示。图5为不同粒子图像直径条件下PIV技术的测量误差剖面，从图5可见各测次结果垂向分布与图2相似：在距床面约16像素范围内平均偏差和均方根误差受床面影响均显著增大；其他区域的平均偏差以20像素为周期沿垂向正负交替变化，均方根误差沿垂向基本不变。在相同垂向位置，平均偏差随粒子图像直径增大而略微减小，粒径为3像素时的最大偏差约±0.02像素，均方根误差在粒径为3像素时稍大，其他测次则维持在0.04像素以内。
图 5 不同粒子图像直径条件下的PIV技术测量误差（k=0.1, Nppp=0.02）Fig. 5 Measurement errors of PIV technique at various particle image diameter conditions (k=0.1,Nppp=0.02)下载:
全尺寸图片
图6为不同粒子图像直径条件下得出的PTV测量误差剖面。各测次在远离床面的区域均无平均偏差，均方根误差则沿垂向以10为周期变化；而在距床面约2像素范围内，平均偏差及均方根误差受床面影响均有增大，但平均偏差仅约0.02像素。在相同垂向位置，均方根误差随粒子图像直径的增大而先减小后变大，粒径为5像素时均方根误差最小，约为0.03像素，表明粒子中心识别算法的精度受粒子图像直径影响，粒径为5像素时识别精度最高，这与Kähler等[29]得到的结果一致。
图 6 不同粒子图像直径条件下的PTV技术测量误差（k=0.1, Nppp=0.02）Fig. 6 Measurement errors of PTV technique at various particle image diameter conditions (k=0.1,Nppp=0.02)下载:
全尺寸图片
图7为不同粒子图像直径条件下LCPIV技术的测量误差分布。
图 7 不同粒子图像直径条件下的LCPIV技术测量误差（k=0.1, Nppp=0.02）Fig. 7 Measurement errors of LCPIV technique at various particle image diameter conditions (k=0.1,Nppp=0.02)下载:
全尺寸图片
从图7中可看出，在远离床面影响的区域，各测次平均偏差及均方根误差沿垂向呈周期性变化，其中：平均偏差的大小随粒子图像直径的增大而减小，粒径为3像素时的最大偏差为±0.06像素；均方根误差的振幅相对较小，但同样随粒径增大而缓慢减小，粒径为3像素时的最大值约0.07像素。在靠近床面的区域，LCPIV技术受床面影响较为显著，对比不同测次的误差分布曲线可知，床面影响高度随粒径增大而逐渐扩大，最大影响高度与粒子图像的半径相当。粒径与床面影响高度的关系与LCPIV技术的原理相符：由于粒子呈刚体运动，在距床面小于粒子半径的区域内，LCPIV判读直线上更多的图像属于中心距床面更远（运动更快）的粒子，使该直线实测流速大于真实流速。
2.4 速度梯度对测量结果的影响
图8为不同速度梯度条件下PIV测量结果的平均偏差和均方根误差。根据第2.2和2.3节图像参数影响分析，PIV技术分别在dτ=3且Nppp=0.01，dτ=9且Nppp=0.04时具有最差和最好的测量精度，故图8中仅展示这两种图像参数对应的结果。在远离床面的区域，平均偏差随速度梯度无明显变化规律，基本保持在±0.02像素；均方根误差随速度梯度的增大而变大，取值范围介于0.02～0.05像素。若以Eb±2Erms度量各技术测量结果的综合精度，则在本文研究参数范围内，PIV技术在不受床面影响时的综合测量精度优于±0.12像素。在床面附近，平均偏差及均方根误差均随着速度梯度的增大而变大，最优测次的综合测量误差已达±0.25像素；同时，床面影响高度由k=0.1时的1倍判读窗口高度（16像素）进一步扩大到k=0.2时的约26像素。
图 8 不同速度梯度条件下的PIV测量误差对比Fig. 8 Measurement errors of PIV technique at different velocity gradient conditions下载:
全尺寸图片
图9展示了不同速度梯度条件下PTV测量结果的平均偏差和均方根误差。根据第2.2和2.3节图像参数影响分析，PTV技术分别在dτ=5且Nppp=0.01，dτ=9且Nppp=0.04时具有最好和最差的测量精度，故仅展示这两种图像参数对应的结果。图9中结果表明，在粒子图像参数较优的条件下，PTV测量结果的平均偏差、均方根误差和床面影响范围几乎不随速度梯度发生变化，此时，离床面2像素以外的平均偏差可以忽略，均方根误差小于0.01像素，综合测量精度可达±0.02像素。当粒子图像参数较差时，速度梯度增大会导致均方根误差和床面附近的平均偏差变大，最不利条件下的测量误差约为±0.15像素。
图 9 不同速度梯度条件下的PTV测量误差对比Fig. 9 Measurement errors of PTV technique at different velocity gradient conditions下载:
全尺寸图片
在床面附近，PTV测量结果的平均偏差和均方根误差在最不利参数条件下分别约为0.06像素和0.09像素，使PTV技术在床面附近的最大测量误差达到±0.24像素，但在粒子图像参数较优时的综合测量误差仅约±0.05像素。对比图9中不同粒子图像直径条件下的测量结果可以发现，PTV技术的床面影响高度与粒子图像直径有关，具体数值约等于粒子图像的半径。当粒子中心与床面距离小于半径时，粒子亮度受床面影响不再呈高斯分布，使基于掩膜相关法的粒子中心识别算法无法准确定位粒子中心，这是导致床面附近测量误差较大的原因。
图10为不同速度梯度条件下LCPIV测量结果的平均偏差和均方根误差。根据第2.2和2.3节粒子图像参数影响分析，由于LCPIV的平均偏差和均方根误差受图像参数影响规律不一致，此处分别选取dτ=3且Nppp=0.01，dτ=9且Nppp=0.04这两种粒子图像参数条件下的结果进行分析。在远离床面区域，平均偏差受速度梯度影响不显著，但均方根误差具有随速度梯度增大而变大的趋势，所有结果中最优和最差测次的测量精度分别约为±0.11像素和±0.30像素。在床面附近，不同速度梯度条件下测量结果受床面影响的高度约等于粒子图像半径，但速度梯度越大，相同高度处的测量精度越差，床面附近的综合测量误差大于±0.45像素。
图 10 不同速度梯度条件下的LCPIV测量误差对比Fig. 10 Measurement errors of LCPIV technique at different velocity gradient conditions下载:
全尺寸图片
综合不同粒子图像密度、直径及速度梯度条件下的测量结果可知：当不受床面影响时，PIV技术测量结果受粒子图像和流动参数的影响较小，综合测量精度优于±0.12像素；PTV技术测量结果受粒子图像参数影响较大，但测量结果无偏差，且在粒子图像密度较小和粒径适中时精度可达±0.02像素；LCPIV技术测量结果受粒子图像和流动参数的影响均较大，综合测量精度不及PIV和PTV技术。在床面附近，各测量技术产生偏差和更大的均方根误差，但不同技术受床面影响的程度和范围不同，PIV技术的受影响范围不小于判读窗口高度，PTV和LCPIV技术的受影响范围约等于粒子图像半径，这与Kähler等[29]针对PIV和PTV在近壁区测量精度的研究结论相吻合。综上，PTV技术空间分辨率高，当其测点与床面的距离大于粒子图像半径时，可以实现无偏差、高精度测量，是瞬时床面切应力的最佳测量技术。
3. 明渠紊流瞬时床面切应力测量试验
3.1 试验介绍
在长20 m、宽0.3 m的高精度明渠玻璃水槽中开展了明渠均匀紊流试验，利用自主研发的高频PIV系统开展粒子图像采集。水流示踪粒子为粒径10 μm的空心玻璃微珠，采用单脉冲能量10 mJ的脉冲激光器形成厚度1 mm的片光照亮水槽中垂面，并利用分辨率为1 920×1 080像素的高速相机搭配焦距135 mm镜头采集粒子图像，关于试验设备的详细介绍可参见文献[30]。共在两种水流条件下开展了测量，水流和图片采样参数见表2，其中：S为水槽底坡；h为试验水深；u*为根据对数律拟合方法得出的摩阻流速；Re=UR/v为雷诺数，U为断面平均流速，R为水力半径，v为水的运动黏度系数；Reτ= u*h/v为摩阻雷诺数；Rs为图像分辨率；Fs为图片采样帧频；N为独立粒子图片对数。根据表2所列水流条件和图片采样参数，C1和C2组次均为紊流流态，黏性底层对应的图像高度分别约为10和8像素。
表 2 明渠均匀紊流试验水流条件及图片采样参数Table 2 Flow condition and image sampling parameters for open channel flow experiments组次 S h/cm u*/(cm·s–1) Re Reτ Rs /(像素·mm–1) Fs/Hz N C1 0.001 2.5 1.75 8 007 510 37.10 1 200 5 000 C2 0.002 2.5 2.11 11 462 613 38.45 2 000 5 000 利用第1.3节介绍的基于PIV的PTV算法对床面附近的粒子图像进行判读，得到黏性底层内的瞬时流速矢量场，再按式（2）计算得到对应床面测点的瞬时切应力。将全部流场中的瞬时床面切应力进行统计分析，得到床面切应力的均值、离散系数、偏度系数和峰度系数。为了与PTV测量结果进行对比，利用第1.2节介绍的PIV算法对粒子图像进行判读，获得2维瞬时流场序列，并计算纵向时均流速和纵向紊动强度剖面。同时，为验证本文试验结果的准确性，使用Hoyas等[31]提供的摩阻雷诺数为550的槽道流直接数值模拟（DNS）数据与测量结果进行对比分析。
3.2 床面位置的确定
由式（2）可知，瞬时床面切应力的准确测量不仅需要高精度纵向流速分量，还需要确定床面位置以获得测点高度。床面在粒子图中的图案及其拟合确定方法如图11所示。图11（a）为一张典型粒子图片，图中底部白色亮带为示踪粒子在床面沉积后被激光照亮形成的图案，由于该亮带的宽度可达数像素，必须使用特殊处理方法才可以准确地确定床面位置。Kähler[11]及Nguyen[15]等曾分别提出从粒子图片中提取床面位置的数字图像处理方法，这些方法的精度在0.4～1.0像素之间，但对粒子图片成像条件有较高要求，难以普适应用。
图 11 床面在粒子图片中的图案及其拟合确定方法Fig. 11 Bed pattern in the particle image and the fitting method of its position下载:
全尺寸图片
本文采用迭代拟合黏性底层内时均流速剖面的方法确定床面位置。首先，根据粒子图片中的亮带得到床面位置初估值yb0，并采用对数律拟合方法确定摩阻流速初估值u*0。其次，定义初估床面以上第i个测点为起点、黏性底层顶端最后一个测点为终点，提取两点之间的实测时均流速剖面，并利用式（4）对剖面进行拟合；由于黏性底层顶端随床面位置和摩阻流速变化，需要反复迭代直至床面位置和摩阻流速的变化小于给定阈值，得到迭代结果u*i及ybi。然后，重新选择初估床面上方第i+1个测点作为剖面起点，重复上述迭代拟合，直至可提取的时均流速测点数少于3个。最后，选择所有迭代结果中摩阻流速最大时对应的床面位置作为最终床面位置。图11（b）为根据C1组次时均流速分布迭代确定床面位置的结果，通过对迭代拟合过程的分析发现本方法的精度可达±0.5像素。
$$ U{\text{ = }}\frac{{u_{\text{*}}^2}}{\nu }\left( {y + {y_{\rm{b}}}} \right) $$ (4) 3.3 床面切应力的统计特征
首先对PTV测得的纵向瞬时流速的准确性进行评估。图12为利用PTV和PIV技术测量得到的纵向时均流速和纵向紊动强度垂向分布。
图 12 实测纵向时均流速及纵向紊动强度分布Fig. 12 Measured profiles of time-averaged streamwise velocity and streamwise turbulent intensity下载:
全尺寸图片
图12中同时得出了相似摩阻雷诺数的槽道流DNS结果；为便于展示，图12（a）和12（b）中C2组次测量结果的纵坐标分别向上平移了10个单位和2个单位。图12（a）中，除靠近床面的1～2个测点外，PTV及PIV测得的时均流速剖面均与DNS结果吻合极好，验证了PTV测量结果在床面影响范围之外的无偏性。图12（b）中，除最靠近床面的少数测点外，PTV技术同样高精度的测得了纵向紊动强度的垂向分布。对比C1和C2组次的测量结果可知，随着时均速度梯度由C1组次的0.18增大到C2组次的0.27，纵向流速测量结果受床面影响的范围由1个像素增大到3个像素，由于试验拍摄的粒子图像平均直径约5像素，上述结果与根据合成图片分析得到的规律相符。由此可见，在离床面远于粒子图像半径的流层中，PTV技术无偏、高精度地测得了明渠紊流黏性底层内的纵向瞬时流速。
利用各组次位于黏性底层（y+≤5）且不受床面影响的测点的纵向瞬时流速和垂向位置计算瞬时床面切应力，并统计床面切应力的时均值τo、离散系数
$\tau _{{\rm{rms}}}^ + = {{{\tau _{{\rm{rms}}}}} \mathord{\left/ {\vphantom {{{\tau _{{\rm{rms}}}}} {{\tau _{\rm{o}}}}}} \right. } {{\tau _{\rm{o}}}}}$ 、偏度系数Sτ及峰度系数Fτ列于表3；表3中同时给出了利用对数律拟合方法[5]和阻力平衡方法τo=ρgRS得出的床面切应力时均值τo1和τo2，其中对数律卡门常数κ=0.41，积分常数B=5.29。对于时均床面切应力，表3中结果表明，对数律拟合方法具有较高精度，而阻力平衡法则由于水深及水面坡度难以准确测量而存在超过10%的误差。床面切应力离散系数的试验结果约为0.4，与前人在壁面紊流中得出的结果一致[19,32]。偏度系数和峰度系数分别大于0.7和3.4并随着摩阻雷诺数缓慢增大，表明瞬时床面切应力不符合高斯分布。Liu[4]和Gubian[33]等近期分别通过试验得出壁面紊流床面切应力的偏度系数和峰度系数随雷诺数增大并分别趋近于1.4和6.6，这与本文得到的结果吻合。上述分析进一步表明，本文提出的基于PTV的测量方法可以准确开展明渠紊流瞬时床面切应力的试验测量。表 3 实测床面切应力的统计参数Table 3 Statistics of the measured bed shear stress组次 τo/Pa τo1/Pa τo2/Pa $\tau _{{\rm{rms}}}^ + $ Sτ Fτ C1 0.31 0.32 0.25 0.41 0.71 3.46 C2 0.44 0.44 0.49 0.39 0.73 3.85 4. 结 论
利用粒子图像技术测量黏性底层内的瞬时纵向流速是间接测量瞬时床面切应力的重要方法。针对PIV、PTV和LCPIV这3种常用粒子图像测量技术的测量精度和适用条件，利用人工合成粒子图片开展影响因素分析和对比，并通过明渠紊流试验开展最优技术的初步应用和验证，得到以下主要结论：
1）各粒子图像测量技术在床面附近的测量均受床面影响，表现为：出现平均偏差和均方根误差显著增大；床面影响范围为紧贴床面的薄流层，PIV测量时流层厚度不小于判读窗口垂向尺寸，PTV及LCPIV测量时等于粒子图像半径。
2）在不受床面影响的流体区域，PIV技术测量精度受粒子图像密度、直径及速度梯度影响较小，测量精度高于±0.12像素，PTV技术的测量结果在图像密度较低且粒径适中时的精度达±0.02像素，LCPIV技术测量精度相对较差。
3）当测点远于粒子图像半径且位于黏性底层内时，利用PTV技术和合理的床面位置确定方法，可以无偏、高精度实现瞬时床面切应力的测量。
4）明渠紊流的瞬时床面切应力不符合高斯分布，在本文试验范围内，瞬时床面切应力的离散系数约0.4，偏度系数和峰度系数随雷诺数增大而逐渐变大。
-
图 1 多重网格迭代及虚拟窗口技术示意图
Fig. 1 Schematic diagram of multigrid iterative and fictitious window techniques
下载:
全尺寸图片
图 2 不同粒子图像密度条件下PIV技术测量误差（k=0.1, dτ=3）
Fig. 2 Measurement errors of PIV technique at various particle image density conditions (k=0.1, dτ=3)
下载:
全尺寸图片
图 3 不同粒子图像密度条件下PTV技术测量误差（k=0.1, dτ=3）
Fig. 3 Measurement errors of PTV technique at various particle image density conditions (k=0.1, dτ=3)
下载:
全尺寸图片
图 4 不同粒子图像密度条件下LCPIV技术测量误差（k=0.1,dτ=3）
Fig. 4 Measurement errors of LCPIV technique at various particle image density conditions (k=0.1, dτ=3)
下载:
全尺寸图片
图 5 不同粒子图像直径条件下的PIV技术测量误差（k=0.1, Nppp=0.02）
Fig. 5 Measurement errors of PIV technique at various particle image diameter conditions (k=0.1,Nppp=0.02)
下载:
全尺寸图片
图 6 不同粒子图像直径条件下的PTV技术测量误差（k=0.1, Nppp=0.02）
Fig. 6 Measurement errors of PTV technique at various particle image diameter conditions (k=0.1,Nppp=0.02)
下载:
全尺寸图片
图 7 不同粒子图像直径条件下的LCPIV技术测量误差（k=0.1, Nppp=0.02）
Fig. 7 Measurement errors of LCPIV technique at various particle image diameter conditions (k=0.1,Nppp=0.02)
下载:
全尺寸图片
图 8 不同速度梯度条件下的PIV测量误差对比
Fig. 8 Measurement errors of PIV technique at different velocity gradient conditions
下载:
全尺寸图片
图 9 不同速度梯度条件下的PTV测量误差对比
Fig. 9 Measurement errors of PTV technique at different velocity gradient conditions
下载:
全尺寸图片
图 10 不同速度梯度条件下的LCPIV测量误差对比
Fig. 10 Measurement errors of LCPIV technique at different velocity gradient conditions
下载:
全尺寸图片
图 11 床面在粒子图片中的图案及其拟合确定方法
Fig. 11 Bed pattern in the particle image and the fitting method of its position
下载:
全尺寸图片
图 12 实测纵向时均流速及纵向紊动强度分布
Fig. 12 Measured profiles of time-averaged streamwise velocity and streamwise turbulent intensity
下载:
全尺寸图片
表 1 基于PIV的松弛迭代参数建议取值
Table 1 Suggested values for parameters of iterative relaxation method based on PIV
预估位移 邻域半径Rn 最大位移Rs 容许偏差Rc Vi >1.5Vi >1.0Vi <1.5Vi 表 2 明渠均匀紊流试验水流条件及图片采样参数
Table 2 Flow condition and image sampling parameters for open channel flow experiments
组次 S h/cm u*/(cm·s–1) Re Reτ Rs /(像素·mm–1) Fs/Hz N C1 0.001 2.5 1.75 8 007 510 37.10 1 200 5 000 C2 0.002 2.5 2.11 11 462 613 38.45 2 000 5 000 表 3 实测床面切应力的统计参数
Table 3 Statistics of the measured bed shear stress
组次 τo/Pa τo1/Pa τo2/Pa $\tau _{{\rm{rms}}}^ + $ Sτ Fτ C1 0.31 0.32 0.25 0.41 0.71 3.46 C2 0.44 0.44 0.49 0.39 0.73 3.85 -
[1] Örlü R,Vinuesa R.Instantaneous wall-shear-stress measurements:Advances and application to near-wall extreme events[J].Measurement Science and Technology,2020,31(11):112001. doi: 10.1088/1361-6501/aba06f [2] Gnanamanickam E P,Nottebrock B,Grosse S,et al.Measurement of turbulent wall shear-stress using micro-pillars[J].Measurement Science and Technology,2013,24(12):124002. doi: 10.1088/0957-0233/24/12/124002 [3] Pabon R J,Ukeiley L,Sheplak M,et al.Characteristics of turbulent boundary layer large scale motions using direct fluctuating wall shear stress measurements[J].Physical Review Fluids,2018,3(11):114604. doi: 10.1103/PhysRevFluids.3.114604 [4] Liu Y,Klaas M,Schröder W.Measurements of the wall-shear stress distribution in turbulent channel flow using the micro-pillar shear stress sensor MPS3[J].Experimental Thermal and Fluid Science,2019,106:171–182. doi: 10.1016/j.expthermflusci.2019.04.02 [5] 刘春晶,李丹勋,王兴奎.明渠均匀流的摩阻流速及流速分布[J].水利学报,2005,36(8):950–955. doi: 10.13243/j.cnki.slxb.2005.08.010 Liu Chunjing,Li Danxun,Wang Xingkui.Experimental study on friction velocity and velocity profile of open channel flow[J].Journal of Hydraulic Engineering,2005,36(8):950–955 doi: 10.13243/j.cnki.slxb.2005.08.010 [6] 曹志刚,王逸伦,国振,等.考虑渗流效应的冲流带泥沙起动机理研究[J].水科学进展,2019,30(4):114–126. doi: 10.14042/j.cnki.32.1309.2019.04.012 Cao Zhigang,Wang Yilun,Guo Zhen,et al.Study on the sediment initiation considering the seepage in the swash zone[J].Advance in Water Research,2019,30(4):114–126 doi: 10.14042/j.cnki.32.1309.2019.04.012 [7] 李福建,黄胜,谢玉姝,等.顺直型展宽航槽沿程水沙运动规律试验研究[J].工程科学与技术,2018,50(4):96–103. doi: 10.15961/j.jsuese.201700963 Li Fujian,Huang Sheng,Xie Yushu,et al.Study of flow and sediment motions in a straight widening navigation channel[J].Advanced Engineering Sciences,2018,50(4):96–103 doi: 10.15961/j.jsuese.201700963 [8] 徐华.波流作用下床面切应力及挟沙能力研究[D].南京:南京水利科学研究院,2016. Xu Hua.Research on bed shear stress and sediment-carrying capacity under wave-current interaction[D].Nanjing:Nanjing Hydraulic Research Institute,2016. [9] Bagherimiyab F,Lemmin U.Large-scale coherent flow structures in rough-bed open-channel flow observed in fluctuations of three-dimensional velocity,skin friction and bed pressure[J].Journal of Hydraulic Research,2018,56(6):806–824. doi: 10.1080/00221686.2017.1416686 [10] Ghouila–Houri C,Talbi A,Viard R,et al.Unsteady flows measurements using a calorimetric wall shear stress micro-sensor[J].Experiments in Fluids,2019,60(4):1–10. doi: 10.1007/s00348-019-2714-5 [11] Kähler C J,Scholz U,Ortmanns J.Wall-shear-stress and near-wall turbulence measurements up to single pixel resolution by means of long-distance micro-PIV[J].Experiments in Fluids,2006,41(2):327–341. doi: 10.1007/s00348-006-0167-0 [12] 钟强,王兴奎,苗蔚,等.高分辨率粒子示踪测速技术在光滑明渠紊流黏性底层测量中的应用[J].水利学报,2014,45(5):513–520. doi: 10.13243/j.cnki.slxb.2014.05.002 Zhong Qiang,Wang Xingkui,Miao Wei,et al.High resolution PTV system and its application in the measurement in viscous sublayer in smooth open channel flow[J].Journal of Hydraulic Engineering,2014,45(5):513–520 doi: 10.13243/j.cnki.slxb.2014.05.002 [13] de Silva C M,Gnanamanickam E P,Atkinson C,et al.High spatial range velocity measurements in a high Reynolds number turbulent boundary layer[J].Physics of Fluids,2014,26(2):025117. doi: 10.1063/1.4866458 [14] Li Mogeng,de Silva C M,Rouhi A,et al.Recovery of wall-shear stress to equilibrium flow conditions after a rough-to-smooth step change in turbulent boundary layers[J].Journal of Fluid Mechanics,2019,872:472–491. doi: 10.1017/jfm.2019.351 [15] Nguyen C V,Nguyen T D,Wells J C,et al.Interfacial PIV to resolve flows in the vicinity of curved surfaces[J].Experiments in Fluids,2010,48(4):577–587. doi: 10.1007/s00348-010-0824-1 [16] Sheng J,Malkiel E,Katz J.Using digital holographic microscopy for simultaneous measurements of 3D near wall velocity and wall shear stress in a turbulent boundary layer[J].Experiments in Fluids,2008,45(6):1023–1035. doi: 10.1007/s00348-008-0524-2 [17] Willert C E.High-speed particle image velocimetry for the efficient measurement of turbulence statistics[J].Experiments in Fluids,2015,56(1):1–17. doi: 10.1007/s00348-014-1892-4 [18] Willert C E,Cuvier C,Foucaut J M,et al.Experimental evidence of near-wall reverse flow events in a zero pressure gradient turbulent boundary layer[J].Experimental Thermal and Fluid Science,2018,91:320–328. doi: 10.1016/j.expthermflusci.2017.10.03 [19] Wang Jianjie,Pan Chong,Wang Jinjun.Characteristics of fluctuating wall-shear stress in a turbulent boundary layer at low-to-moderate Reynolds number[J].Physical Review Fluids,2020,5(7):074605. doi: 10.1103/PhysRevFluids.5.074605 [20] 陈启刚.基于高频PIV的明渠湍流涡结构研究[D].北京:清华大学,2014. Chen Qigang.High-frequency measurement of vortices in open channel flow with particle image velocimetry[D].Beijing:Tsinghua University,2014. [21] Brady M R,Raben S G,Vlachos P P.Methods for digital particle image sizing (DPIS):Comparisons and improvements[J].Flow Measurement & Instrumentation,2009,20(6):207–219. doi: 10.1016/j.flowmeasinst.2009.08.001 [22] Ohmi K,Li Hangyu.Particle-tracking velocimetry with new algorithms[J].Measurement Science and Technology,2000,11(6):603–616. doi: 10.1088/0957-0233/11/6/303 [23] Kim H B,Lee S J.Performance improvement of two-frame particle tracking velocimetry using a hybrid adaptive scheme[J].Measurement Science & Technology,2002,13(4):573–582. doi: 10.1088/0957-0233/13/4/321 [24] Kähler C J,Scharnowski S,Cierpka C.On the resolution limit of digital particle image velocimetry[J].Experiments in Fluids,2012,52(6):1629–1639. doi: 10.1007/s00348-012-1280-x [25] Cierpka C,Scharnowski S,Kähler C J.Parallax correction for precise near-wall flow investigations using particle imaging[J].Applied Optics,2013,52(12):2923–2931. doi: 10.1364/AO.52.002923 [26] Okamoto K,Nishio S,Saga T,et al.Standard images for particle-image velocimetry[J].Measurement Science and Technology,2000,11(6):685–691. doi: 10.1088/0957-0233/11/6/311 [27] Adrian R J,Westerweel J.Particle image velocimetry[M].New York:Cambridge University Press,2011. [28] Lei Y C,Tien W H,Duncan J,et al.A vision-based hybrid particle tracking velocimetry (PTV) technique using a modified cascade correlation peak-finding method[J].Experiments in Fluids,2012,53(5):1251–1268. doi: 10.1007/s00348-012-1357-6 [29] Kähler C J,Scharnowski S,Cierpka C.On the uncertainty of digital PIV and PTV near walls[J].Experiments in Fluids,2012,52(6):1641–1656. doi: 10.1007/s00348-012-1307-3 [30] Duan Yangchong,Zhang Peng,Zhong Qiang,et al.Characteristics of wall-attached motions in open channel flows[J].Physics of Fluids,2020,32(5):055110. doi: 10.1063/5.0004323 [31] Hoyas S,Jiménez J.Scaling of the velocity fluctuations in turbulent channels up to Re=2003[J].Physics of Fluids,2006,18(1):011702. doi: 10.1063/1.2162185 [32] Örlü R,Schlatter P.On the fluctuating wall-shear stress in zero pressure-gradient turbulent boundary layer flows[J].Physics of Fluids,2011,23(2):021704. doi: 10.1063/1.3555191 [33] Gubian P A,Stoker J,Medvescek J,et al.Evolution of wall shear stress with Reynolds number in fully developed turbulent channel flow experiments[J].Physical Review Fluids,2019,4(7):074606. doi: 10.1103/PhysRevFluids.4.074606
