任意运动变形边界流场测量PIV算法的研究：薄膜涡激振动

曹洪才，施圣贤，刘应征

（上海交通大学 动力机械及工程教育部重点实验室，上海 200240）

0 引 言

随着人们对仿生机械、微型飞行器等研究的日益深入，以前只关注机械结构设计的单一研究思路被摈弃，研究人员开始结合流体力学、结构力学等多学科进行优化设计，深入研究边界变形与流体流固耦合的物理机制及其对气动性能的影响［1］。此外，流固耦合现象在自然界和工程实践中很常见，如鸟类振翅、旗帜飘动、心脏瓣膜开合过程、机翼颤振等［2-4］。在针对这类问题的数值建模和计算分析中，由于涉及流体和固体控制方程，两种非线性控制方程存在强烈耦合，使得这类问题的理论分析和数值计算面临很大的挑战，所建立的数值模型及其计算结果也需要高精度的实验数据进行验证［5］。因此，精确测得强烈非线性流固耦合过程中的流场信息，尤其是运动变形边界处的局部流动特征，对于流固耦合机制的理论分析和数值计算都尤为重要。激光粒子图像测速技术（PIV）由于具有全场、非接触测量等特点，被广泛应用于流场实验研究中。然而，将传统的PIV 测量技术用于存在任意运动变形边界问题的流固耦合流场测试时，互相关计算窗口中可能包含流固边界，会导致不正确甚至错误的测量结果［6］。此外，传统PIV 算法在用于流固耦合流场测试时，无法获得流固耦合边界的运动特征。因此，将PIV 应用于流固耦合这种复杂流场测试时，需要充分考虑这些特点，发展能智能处理任意运动变形边界条件下的PIV 流场测试技术［7］。对于强流固耦合问题，由于存在边界时变、变形大、边界运动复杂等特点，主要难点在于大变形运动边界的智能识别、精确数字化，以及边界附近局部流场信息的精确定量分析。针对这一问题，采用基于Radon变换的滑动窗口边界识别方法以精确获取任意运动变形边界信息，并应用贴体自适应互相关计算窗口测得瞬态流场信息，大大提高了流固耦合问题的PIV 测试能力。最后，将提出的新型PIV 算法应用于柔性薄膜涡激共振流场的PIV 实验测试，获得了薄膜在各种大变形状态下的高分辨率瞬态流场数据。

1 任意运动变形边界流场测量PIV算法

1.1 算法基本原理

传统PIV 算法一般只适合于模型固定的流场测量，而流固耦合相关问题具有模型边界时变、变形量大等特点，对传统PIV 流场测量提出了挑战。采用基于边界信息的贴体自适应PIV 互相关计算窗口可以获得高精度近壁流场信息，但这一方法的前提是需要从“高噪声”粒子图像中提取运动变形边界的精确位置信息，即任意运动变形边界的识别问题。对于运动边界的识别而言，这里所说的“噪声”就是示踪粒子本身。Radon变换是一种检测直线的图像处理算法，对灰度图像I（x，y）进行Radon变换，即在灰度图像平面内沿不同的直线做线积分，可以理解为灰度图像在ρ-θ空间的投影。

直接使用Radon变换难以处理信噪比低的边界信息，因而需要对算法进行必要的修正和改良。提出了基于Radon变换的滑动窗口边界识别方法（如图1所示）：（1）根据初态条件确定粒子图像的初始窗口，进行Radon变换，根据灰度值积分侦测初始窗口中的直线特征；（2）基于前一个窗口识别的直线特征，预测下一个检测窗口的中心位置，并基于模板匹配的互相关方法对预测的窗口中心位置进行修正；（3）基于先期位置信息，建立窗口中边界位置的概率模型，进行图像灰度值的加权，以降低图像中示踪粒子的干扰；（4）经过预测—修正—加权的窗口图像，图像信噪比得到极大提升。对灰度值进行Radon变换，检测低信噪比图像中的直线特征。如此循环，可以识别出边界分割开的若干小直线线段，进行曲线拟合，可以得到粒子图像中运动变形边界的精确位置。

图1 基于Radon变换边界识别方法示意图Fig.1 Interface tracking based on Radon transformation

为验证边界识别算法的精确性，根据已知流场信息和边界运动信息数字合成粒子图像对。在数字合成粒子图像中，边界位置已预先知道，可和采用Radon变换算法识别出来的边界信息进行误差对比，验证边界识别的正确性及精度。图2所示为数字合成粒子图像，红色圆圈为边界识别程序识别出来的边界信息。

图2 数字合成粒子图像及边界识别示意图Fig.2 Sketch of synthetic PIV image and interface tracking

为考察粒子浓度大小也即“噪声”对边界识别精度的影响，生成了不同粒子浓度的数字合成图像来模拟不同信噪比环境下边界识别的适用性。如图3所示为32pixel×32pixel窗口中粒子数量分别为5 个和20个的数字合成图像示意图。

图3 不同信噪比的数字合成粒子图像Fig.3 Synthetic PIV images with different SNR

图4 粒子浓度与信噪比、识别精确性的关系曲线Fig.4 SNR and tracking accuracy under different particle density

如图4所示为粒子浓度与数字合成图像信噪比的关系曲线，以及边界信息识别误差大小。随着32pixel×32pixel像素窗口中粒子浓度从5个逐渐增加到20个，数字合成图像的信噪比随之降低。在图像中加入白噪声后，图像的信噪比进一步降低，如红色曲线所示。对边界信息识别的误差绝对值进行了统计，可见识别误差在0.45个像素以内。

图5 传统矩形和贴体自适应PIV 互相关计算窗口示意图Fig.5 Sketch of traditional and body-fitted window for cross-correlation

传统PIV 算法采用规则正交的矩形互相关计算窗口，计算具有变形边界的流固耦合流场时，会出现对边界附近不同侧粒子进行互相关计算的问题，引起当地流速的误判。如图5所示，为提高具有运动边界流场的测量精度，在识别出流固耦合运动边界信息之后，基于运动变形边界信息生成贴体自适应互相关计算窗口，然后采用标准互相关算法计算当地速度矢量。采用此算法，可精确获取流固耦合边界附近流场信息。

1.2 基于数字合成粒子图像的算法验证

为验证基于边界信息的贴体自适应窗口改进PIV 算法的正确性，首先用计算流体手段计算了一个具有弯曲边界的流场，其速度分布是已知的，称为标准速度场。然后，采用数字合成PIV 粒子图像对的方法，数字合成了两帧粒子图像。通过互相关运算计算出两帧粒子图像的速度场，与标准速度场进行对比，可验证改进PIV 算法的正确性。

图6 改进PIV 算法的数字合成图像验证Fig.6 Velocity accuracy for improved PIV algorithm by synthetic PIV images

如图6（a）所示为利用改进PIV 算法计算出的速度云图和矢量图，而图6（b）是通过数值计算获得的标准流场速度云图和矢量图。图中的小方格是生成的贴体自适应互相关计算窗口，通过互相关运算在每个窗口生成一个速度矢量，如图中黑色小箭头所示。对比两图对应窗口的速度矢量，可以发现改进PIV算法计算的流场流态与标准流场基本一致。图中颜色代表该互相关计算窗口计算出的速度值，在运动边界附近，贴体自适应窗口经过互相关运算的速度值相对误差在5%以内。总体来看，改进PIV 算法计算出的速度云图与标准的速度场非常接近，计算结果基本正确，尤其在大变形边界附近没有出现传统PIV 算法的明显问题。

2 实验应用

2.1 水槽及实验模型

为展示所发展的任意运动变形边界条件下流场的PIV 实验测量效果，选择布置在钝体尾迹区域中的柔性薄膜涡激共振现象作为实验对象。流场测量实验在低速闭式循环水洞［8］中进行，由进口水箱、整流段、收缩段、实验段、循环水箱、磁力泵、变频仪以及管路系统等组成。进口水箱处设置有孔板，以削弱水泵出水的振荡，达到初步整流的效果。水槽收缩段长450mm，收缩型线采用三次方曲线。实验段横截面为正方形，截面尺寸为100mm×100mm。为了减小水泵振动对流场的影响，采用变频仪控制的磁力泵驱动，可以通过改变变频仪的频率来调节流速。水槽的具体设计及结构详见文献［9］。

实验时，在封闭水槽实验段的中央放置一个刚性钝体，钝体采用直径D=8mm 的方柱。钝体下游一倍D处放置一端固定的柔性薄膜，柔性薄膜的上下端和另一端自由。柔性压电薄膜的尺寸为：长度L=41mm，宽度b=10mm，厚度h=0.1mm。自由来流速度u0=1.2m／s，雷诺数Re=u0D／μ=9930。柔性薄膜与钝体尾流之间的流固耦合作用示意图如图7所示。流场上游经过方柱绕流后，在下游所自然形成的旋涡脱落造成柔性薄膜变形运动，驱使柔性薄膜振动。由于柔性薄膜的结构振动使流场发生改变，流场变化和薄膜变形运动在此过程中形成强耦合作用。

图7 钝体尾流柔性薄膜涡激振动示意图Fig.7 Sketch of vortex induced vibration of flexible film

PIV 实验时，示踪粒子采用密度为1.03g／mm3的空心玻璃微珠，直径为20～30μm。照明采用半导体连续式激光器，测量区域片光厚度小于1mm。为提高空间分辨率，实验中采用了1600万像素的CCD相机。通常，相机的时间分辨率和空间分辨率是矛盾的。空间分辨率高的相机，其时间分辨率一般较低。采用了高空间分辨率的CCD相机，其采样频率较低为1f／s，采集窗口大小为4904pixel×3280pixel。测得薄膜各种运动变形状态时的粒子图像后，采用Radon变换对粒子图像中的薄膜形态进行图像处理精确获取薄膜形态曲线的数字化信息，然后生成贴体自适应图像互相关分析窗口，采用北京立方天地公司的Micro Vec粒子图像分析软件计算得到速度矢量。计算中采用图像偏置及迭代算法［10］。互相关计算的判读窗口大小为64pixel×64pixel，步长为32pixel。互相关运算的峰值采用高斯拟合方法确定［11］，计算精度可达到±0.1pixel。

2.2 实验结果及分析

图8所示为Re=9930下钝体尾涡柔性薄膜涡激振动不同时刻的PIV 粒子图像，图中蓝色所示是程序所识别的边界位置信息。这些PIV 粒子图可以使我们对柔性薄膜涡激振动这种强烈的流固耦合现象有直观的认识。从图中可以看出，方柱脱落的旋涡，由于柔性薄膜的阻隔作用，上下剪切层无法交换能量，旋涡贴附薄膜向下游输运。这样，上下剪切层能量存在差异，柔性薄膜两侧受到压力差而发生弯曲变形，在不同时刻呈现不同的振动状态。

图8 柔性薄膜连续振动的不同时刻粒子图Fig.8 Particle images of voetex induced vibration in different frames

图9 柔性薄膜涡激振动不同时刻的速度云图Fig.9 Contour plot of velocity field and streamline

图9为柔性薄膜涡激振动过程中4个不同时刻的流场速度云图（Re=9930）。图9（a）钝体尾流中出现两个大尺度旋涡，薄膜上方为顺时针旋涡，下方为逆时针旋涡。薄膜上方的旋涡处于发展阶段，距离方柱的位置相对较远。图9（b）中，上方的旋涡向下游输运过程，在其影响下，柔性薄膜发生向下的弯曲变形。同时，在上游已经生成新的旋涡。图9（c）中，钝体上、下缘脱落的旋涡贴附柔性薄膜输运，在柔性薄膜前半段出现两个尺度相同的反向旋涡。在下游，柔性薄膜后半段被下方旋涡迅速向上扬起。图9（d）中，在柔性薄膜上方存在两个旋涡，旋涡向下游输运过程中迫使柔性薄膜发生弯曲变形。而柔性薄膜下方只有一个逆时针旋涡，其能量刚好可以平衡上方的顺时针旋涡。这与文献［12］的计算流体动力学数值模拟结论是一致的。

3 结 论

提出了一种用于测量任意运动变形边界流场的PIV 计算方法，可以智能识别强烈流固耦合问题中的运动变形边界，继而生成基于边界信息的贴体自适应互相关计算窗口，获得任意运动变形边界附近的流场精确信息。采用这一算法，对已知运动变形边界流场信息的数字合成图像进行系统的PIV 分析，发现所获得的变形边界附近的流场信息与已知流场信息高度一致，证明了所提出的PIV 算法的合理性。最后，利用PIV 测量了低速闭式循环水洞中钝体尾涡柔性薄膜涡激振动现象，采用改进PIV 算法对4 种典型薄膜形态下的流场进行了分析。实验测试结果表明，钝体尾流周期脱落的旋涡与尾流中的柔性结构形成强烈的非线性耦合作用，流态特征与薄膜形态变化的关联特征与相关文献结果一致。这说明所提出的任意运动边界流场测量PIV 算法可用于强流固耦合问题的流场实验测量，为流固耦合物理机制理论分析和数值计算方法的验证提供了有力的实验方法。

参考文献：interaction［J］.Annual Review of Fluid Mechanics，2001，33（1）：445-490.

［3］ SARPKAYA T.A critical review of the intrinsic nature of vortex-induced vibrations［J］.Journal of Fluids and Structures，2004，19（4）：389-447.

［4］ GABBAI R D，BENAROYA H.An overview of modelingand experiments of vortex-induced vibration of circular cylinders［J］.Journal of Sound and Vibration，2005，282：575-616.

［5］ 王文全，张立翔，闫妍，等.方柱绕流诱发的弹性薄板流固耦合特性研究［J］.工程力学，2011，28（3）：17-22.

［6］ ADRIAN R J.Particle image techniques for experimental fluid mechanics［J］.Annual Review of Fluid Mechanics，1991，23：261-304.

［7］ JERRY Westerweel.Particle image velocimetry for complex and turbulent flows［J］.Annual Review of Fluid Mechanics，2013，45：409-436.

［8］ 陈建民.钝体尾流压电薄膜涡激共振流动特性的实验研究［D］.［硕士学位论文］.上海：上海交通大学，2013.

［9］ SHI Sheng-xian，LIU Ying-zheng.Flapping dynamics of a low aspect-ratio energy-harvesting membrane immersed in a square cylinder wake［J］.Experimental Thermal and Fluid Science，2013，46：151-161.

［10］WESTERWEEL J，DABIRI D，GHARIB M.The effect of a discrete window offset on the accuracy of crosscorrelation analysis of digital PIV recordings［J］.Experiments in Fluids，1997，23（1）：20-28.

［11］YASUHIKO S，NISHIOS，OKUNO T.A highly accurate iterative PIV technique using a gradient method［J］.Measurement Science and Technology，2000，11（12）：1666-1673.

［12］LIEW K M.A computational approach for predicting the hydroelasticity of flexible structures based on the pressure Poisson equation［J］.International Journal for Numerical Methods in Engineering，2007，72：1560-1583.