风驱动下平板水膜表面波动特性实验研究

2022-04-27胡玖琦王敬鑫朱春玲朱程香刘森云

南京航空航天大学学报 2022年2期

胡玖琦，王敬鑫，朱春玲，朱程香，赵宁，刘森云

（1.南京航空航天大学航空学院，南京 210016；2.南京航空航天大学机械结构力学及控制国家重点实验室，南京 210016；3.中国空气动力研究与发展中心结冰与防除冰重点实验室，绵阳 621000）

当飞机穿过含有过冷水滴的云层时，机体表面容易发生结冰现象，造成飞机性能显著恶化，严重威胁飞行安全［1⁃2］。在霜冰条件下，撞击水滴在迎风面瞬间冻结。然而，在明冰条件或防冰系统工作状态下时，持续撞击的过冷水滴在迎风表面相变并释放潜热，在表面留存未冻结的液态水，并在气流驱动下向后发生回流，影响结冰过程［3］。气流在水膜表面沿流动方向施加压力，当气液界面剪切力克服液体表面黏性时，水膜流动开始形成表面波［4］。水膜自由表面在气流扰动下产生不稳定性，引起连续水膜的表面波动，同时表面波动反作用于气流边界层的流动，进而影响了结冰或防冰表面的热量分布。因此，飞机表面液态水膜的存在不仅增加了积冰的危害，同时增加了冰形预测和防冰系统设计的复杂性。

近几年，随着实验技术进步，研究液态水膜流动的实验逐渐增多。对于连续水膜流动的厚度分布，Muzik 等基于实验和数值方法研究了翼型表面的水膜流动特性及其破碎现象［5］。冷梦尧等利用色散共焦位移计方法测量了风驱动下平板水膜流动，并得到同一位置不同时刻的厚度变化［6⁃7］。对于水膜表面波动特性，Cherdantsev 等基于激光诱导荧光技术（Laser induced fluorescence，LIF）研究了在水平矩形管内高速气流剪切下液膜的流动特性，提出了在气流剪切作用下液膜流动形成波纹三维结构的不同机制［8］。Chang 等利用平面激光诱导荧光方法（Planar laser induced fluorescence，PLIF）分析了水膜厚度和表面波特性［9］。尽管利用各种实验手段对液膜问题进行了广泛的研究，但气流驱动下平板水膜流动是典型的三维问题，波形及分布特征是研究水膜波动特性的重要参数，因此需要一种非接触式测量方法实现水膜流动过程中全视域内三维信息测量和较高的时间分辨率。

数字图像投影（Digital image projection，DIP）技术可以实现气流驱动下水膜流动厚度和表面波动的非侵入式测量。DIP 技术基于立体视觉系统的结构光三角测量原理［10］，用数字图像投影仪将特定的网格图像投影到待测区域，当被测区域发生变化时，网格图像也会随之变动，同时利用高速相机对变形网格图像进行捕捉，通过图像关联算法［11］将实验图像与参考图像对应点的位移矢量转化为三维形状表面的高度数据。Zhang 等利用DIP技术实现了翼型和平板表面水膜/溪流厚度分布测量［12⁃13］。Liu 等利用DIP 技术研究了带有初始结冰粗糙度的翼型对气流驱动下水膜流动的影响［14］。综上所述，大部分研究侧重于水膜厚度测量，对于气流驱动下平板水膜流动表面波动特性研究较少。以往的测量手段针对于单点位置的信号获取，并不能完全展现水膜流动在不同条件下波形变化的特征，因此，需要对全视域范围内水膜的三维信号特征进行分析研究。

本文针对气流驱动下平板水膜表面波形成和发展过程进行了测量，通过搭建特定的平板水膜流动实验台产生连续的流动水膜，使用DIP 技术记录一段时间内的水膜流动图像序列，基于MATLAB图像后处理实现表面波连续变化的三维复原，分析表面波在流动过程中波动特性随风速和水膜雷诺数的关系，从而揭示水膜流动过程中波形特征变化的物理机制，提升冰形预测和防冰系统的精确性。

1 实验设置

为了探究气流驱动下水膜流动表面波动特性，设计并搭建了开放式水膜流动实验台，如图1 所示。高速气流由轴流风机提供，经过蜂窝入口段、稳定段和收缩段加速后进入实验段，同时为了消除风机震动对实验产生的误差，在风道末端与风机之间采用软连接。实验段截面尺寸为长300 mm、宽70 mm，侧面使用亚克力玻璃板，上下表面均有开槽口，顶部安装可拆卸透明观测窗，底部安装铝板作为实验平板。实验基板为铝板，尺寸为长500 mm、宽250 mm，上表面经过抛光处理，表面粗糙度为0.2 μm，下表面与精密升降台固定。在DIP系统标定过程中，调节升降台实现实验基板进行法向移动，从而获取高度位移转换系数［12］。在实验过程中，微调升降台用于保持实验基板与风道实验段底部齐平，以便对实验段进行密封处理。

图1 实验装置Fig.1 Test setup

观测所使用的水源为掺混白色乳胶漆的蒸馏水（体积比例为100∶1），用于增强液体表面漫反射，通过精密齿轮泵驱动，水源经蓄水池从供水箱底部注入实验段。供水口为凹槽，长度70 mm、宽度1 mm，下游侧面边缘倒有圆角，并利用密封圈增加气密性。水膜在气流驱动下沿实验基板向前流动，在下游回水段收集并排出。本研究中用于测量水膜表面波动特性的实验参数如表1 所示。

表1 水膜波动特性测量的实验参数Table 1 Experimental parameters for water film fluctuation characteristic measurement

本文使用DIP 系统定量测量了观测区域内水膜厚度分布。如图2 所示，数字投影仪用于在实验基板表面投影网格图案，当水膜流经实验基板时，投影网格图案出现变形。CCD 高速相机安装在实验段顶部，用于获取网格图像序列。投影仪和相机对齐，垂直于水膜流动方向，以抑制水膜表面强烈的镜面反射。所有实验均在常压室温条件下进行，并待气流和水膜运动稳定一段时间后再进行记录。在本研究中，为实现更精确记录水膜流动瞬态特性，高速相机的拍摄帧率设置为1 000 帧/s。

图2 DIP 系统示意图Fig.2 DIP system

同时为了测定本实验所搭建的DIP 系统的测量精度，对1 mm 标准厚度板进行了测量，测量结果如图3 所示。定义绝对偏差值计算式为

式中：hi为单点的测量计算值，hs为标准件实际厚度值。由图3 表明，被测面积内中间区域误差较小，绝对偏差值在6%以内，边缘区域误差较大，整体均值统计测量误差为9.5%。这是因为投影单元在中心区域投射出的网格点焦距适中，图像清晰，边缘区域聚焦较差导致图像模糊，引起测量误差。因此，在水膜实验和数据处理过程中选取网格清晰的中心区域。

图3 标准板厚度测量误差Fig.3 Error of standard plate thickness measurement

2 结果与分析

2.1 水膜波动的瞬态特性

针对气流驱动下的水膜流动问题，当气流和液膜之间存在速度差时，在气液界面产生扭转应力，发生位移和动量交换，形成界面波动现象，即Kelvin⁃Helmholtz 不稳定性。基于气⁃液两相分层流动的流域划分方式［15⁃16］，Bae 等通过粒子图像测速（Particle image velocimetry，PIV）技术对矩形管道中气⁃液两相流动进行观测［17］，实验结果表明不同风速下，液膜表面呈现4 种不同的波动形态。气流驱动下水膜流动区域主要分为：二维波（2D⁃wave）、三维波（3D⁃wave）、滚动波（Roll⁃wave）和卷吸（Entrainment）。在中等水膜雷诺数条件下，气流速度对水膜流动区域内表面波的典型影响，如图4 所示。定义水膜雷诺数公式为

图4 不同风速下流域波形Fig.4 Shapes of waves at different air speeds

式中：Qm为水膜质量流量，L为水膜润湿区域宽度，μw为水膜的动力黏度系数。

低风速下，气⁃液剪切力克服流体黏性，形成表面波，进入发展段后表现为具有完整周期性的二维波，并在几个周期内保持其形状统一。当风速增大，表面波在横向上发生形状明显改变，发展成为三维波，这是因为气液界面剪切应力增大，表面波界面曲率增大，不稳定性增强，表面波被拉伸变形展现了多样的表面波形态。随着气流速度的持续增大，水膜自由表面不稳定性极大加剧，表面波最终发生破碎，形成大量小股的滚动波，这是由于气液界面剪切应力持续增大，下层液体交换获得更大动量，直至形成漩涡，其瞬态信号内相位速度高，分布随机性强。当气液界面剪切力增加至表面张力无法维持液体保留在液膜内时，表面波顶部的液体会被气流卷吸脱离水膜表面。

水膜在气液界面的质量分布对传热传质、空气边界层流动有显著影响，同时气流边界层的扰动加剧改变自由液面波动。为了更加清晰地表述气流驱动下水膜流动表面波动演化过程，基于高速相机，可以清晰地捕捉水膜流动过程中表面波的传质过程，如图5 所示。可以发现，水膜表面波以孤立行波形式沿气流方向发展，并呈现弯曲弓形。

图5 表面波发展序列图及厚度分布图(Ref=48.34, Ua=33.7 m/s)Fig.5 Sequential images of dynamic surface wave development and thickness distribution (Ref=48.34, Ua=33.7 m/s)

针对流域内不同位置处表面波发展过程分析，结果如图6 所示，定义波速公式为

图6 不同截面处的波速变化(Ref=48.34, Ua=33.7 m/s)Fig.6 Variations of wave velocity at different sections(Ref=48.34, Ua=33.7 m/s)

式中：f为高速相机帧率，xt-Δt为t时刻前一帧波峰极值点所在横坐标（沿流动方向下），xt+Δt为t时刻后一帧波峰极值点所在横坐标。

随着区域厚度的增大，波速明显提高，这是因为气液界面剪切力在不同位置处对水膜表面波动影响差异较大，同时在液膜表面张力作用下，导致水膜厚度分布不均匀，质量在Y=55 mm 所在的轴线区域集中，因而惯性力增大，同时增大迎风面积，气流压力进而在该区域集中，导致表面波速从高质量向低质量方向衰减，进而波形向弓形发展。

2.2 风速对平板水膜表面波动特性影响

相关文献针对平板水膜流动的测量主要针对不同工况下气液界面不稳定性的临界条件展开研究，Chang 等通过PLIF 技术观测了较大风速区间内水膜波动特性，发现水膜平均厚度随风速增大而降低，频率随风速增大而增大［9］，但是缺少对界面波动其他特征参数清晰的描述，如振幅、波速和波峰间距等。因此，本文基于水膜厚度的时间信号和空间信号，通过控制单一变量法，设定水膜雷诺数Ref为48.34 保持不变，得到不同来流风速下平板水膜的波动特性，如图7 所示。选定在供水口附近第一次产生完整波形的位置为观测点。基于水膜自由液面的厚度测量，可以提取关于表面波特征的信息，如图8 所示。结果表明，在供水量不变的条件下，随着风速的增加，气流将带走更多水量，水膜表面波峰幅值高度显著降低，水膜整体变薄，最后趋于稳定。当气流来流速度处于33.7～39.5 m/s 条件下时，水膜厚度较薄，静水面接近底部黏性边界层，表面波表现为多股滚动波形式在湿表面上快速流动，如图4（d，e）所示。同时观察发现，随着减小风速变化区间，频率出现极值点，在中等的来流速度下（＜28 m/s），水膜表面波频率随着风速的增加而增加，呈现由2D 波向3D 波转化趋势，这是因为随着来流风速增大，界面剪切力增强，表面波波长缩短，同时表面波速提升。当来流速度高于28 m/s 时，水膜表面波频率变化不再单调递增，水膜表面出现滚动波，破坏了波动的周期性和均匀性，表面波频率出现下降。当来流速度达到45.5 m/s 时，水膜厚度极为稀薄，此时表面波为高频滚动波随机分布和快速移动，并且向卷吸形式发展。

图7 不同风速下表面波厚度时域变化Fig.7 Wave thickness time traces at different air speeds

图8 风速对表面波波峰和频率影响Fig.8 Effect of air speeds on crest and frequency of surface wave

针对本研究使用DIP 技术，可以实现全视域范围内水膜的厚度特征捕捉，得到不同来流风速下平板水膜的波动特性，如图9 所示。由于记录过程中，水膜流动趋于稳定，所以选择不同时刻水膜流动状态，以保证轮廓的初始点均为谷值（静水深）。可以发现，水膜自由表面静水深随着风速的增加而降低，并且趋于平稳。此外，结果表明，随着风速增加，波峰间距逐渐缩短，波速逐渐提升，如图10 所示。在低来流速度下（＜20 m/s），水膜表面波动呈现规律的2D 波，波峰间距大而波速较低；随着风速增大（20～28 m/s），水膜表面波动由2D 波向3D 波转化，波峰间距缩短而波速升高，直至出现滚动波（＞33.7 m/s），呈现出随机分布、相位速度快的特点。

图9 不同风速下表面波空间分布Fig.9 Spatial distribution of wave at different air speeds

图10 风速对表面波波速和波峰间距影响Fig.10 Effect of air speeds on velocity and crest spacing of surface wave

2.3 水膜雷诺数对平板水膜表面波动特性影响

平板水膜流动表面波动特性受到来流风速的影响，同时受到供水流量Qm的影响，通过引入水膜雷诺数Ref，保持风速Ua为22.2 m/s 不变，建立平板水膜流动表面波动特性参数与Ref关系。在中等来流速度条件下，不同Ref对水膜流动区域内表面波的典型影响，如图11 所示。可以发现，在中等风速条件下不同Ref下水膜表面波动变化较小，呈现出具有周期性的2D 和3D 波状态，因此风速对水膜表面波动占据主导作用，而Ref影响较小，与Chang 的实验结论一致［9］。这是因为随着风速的增加，气液界面扰动不稳定性和气流湍流度增大［9］，进而引发界面失稳，而Ref增大抑制界面扰动，因此自由表面波动稳定性增强。

图11 不同Ref 下流域波形及厚度分布图Fig.11 Shapes of waves and thickness distribution at different Ref

为了探究Ref对自由表面波动的影响，基于DIP 技术，实现水膜厚度的时间信号处理，如图12所示。可以发现，随着Ref的增加，水膜表面波动峰值增大，而频率没有显著变化。这是因为随着供水量增大，静水面厚度升高，底层液膜黏性影响减弱，气液界面剪切力作用增强，导致波峰高度增大。此外，结果表明，随着Ref增加，水膜表面波波速逐渐提升导致视域内波峰数量减少，如图13 所示。这是因为供水流量增大，水膜沿流动方向速度分量明显加大，同时波峰高度较高，迎风面积增大，加快表面波移动速度。

图12 不同Ref 下表面波厚度时域变化Fig.12 Surface wave thickness time traces at different Ref

图13 不同Ref 下水膜波峰和波速变化Fig.13 Effect of Ref on crest and velocity of surface wave

不同Ref下水膜表面波厚度的空间信号，如图14 所示。可以发现，水膜自由表面静水深随着Ref的增加而增加，并且趋于平稳。此外，随着Ref增加，波峰抵达最高点的位置向远离供水口方向移动。这是因为供水流量增大，水膜厚度增大，底层黏性水膜受气流剪切作用减弱，波峰的发展段向后移动，同时表面波波速加快，进而出现波峰值延迟情况。

图14 不同Ref 下水膜表面波空间分布Fig.14 Spatial distribution of surface wave at different Ref

3 讨论

为进一步检验DIP 技术在应用于平板水膜流动中的精确性，引入量纲化为一的水膜厚度与经典的量纲化为一的水膜厚度经验公式进行对比。通过气⁃液剪切应力结合对液膜厚度进行量纲化为一处理，定义量纲化为一的水膜厚度h+为

式中：ρa为气流密度；Ua为来流风速；C为气⁃液界面波的相位速度，即波速；fi为界面剪切系数，本文使用平板层流界面剪切系数［19］，得到气⁃液界面剪切力公式为

由于气⁃液界面波的相位速度，即表面波波速与来流风速相比小于0.01，为简化计算，则在计算中忽略了波速。

表2 列出了相关研究的量纲化为一的水膜厚度计算模型，各模型的计算值与实验值的比较结果如图15 所示。整体来看，量纲化为一的平均水膜厚度随着Ref增大而增大，随着风速的增大而减小。对比发现，在中等风速下（22.2～28.0 m/s）实验值与计算模型符合程度最高，但低风速与高风速下有明显差异，这是因为中风速下表面波动在水膜流动过程具有良好的周期性和稳定性，而低风速下还未完全发展，高风速下出现的滚动波打破了原有的稳定状态。