小河道扑动水翼装置推水流动特性研究

华尔天 陈万前 汤守伟 谢荣盛 郭晓梅 徐高欢

(1.浙江工业大学机械工程学院， 杭州 310023； 2.先进水利装备浙江省工程研究中心， 杭州 310018；3.浙江水利水电学院机械与汽车工程学院， 杭州 310018)

0 引言

平原河网地区地势平缓、水体流速较小，水动力不足引发的水体污染对居民生活质量与身体健康造成了巨大的影响[1]。调水引流作为常用的治理措施之一，能够提升河网水动力，有效改善水环境质量[2-4],但对分流能力小的支杈的水动力状况帮助较小[5]。而且降低转速、增大叶轮直径虽然可以改善泵站在超低扬程下运行效率低、稳定性差[6-8]等问题，却依旧无法适应平原小河道扬程几乎为零的情况。文献[9-10]模拟鱼尾摆动形式研制的仿尾鳍式微型压电泵，具有结构简单、流量大、扬程低等优点，验证了鱼尾摆动形式对水体的推动作用。受此启发，本文提出利用扑动水翼进行水体推动的方法，以提升平原小河道的水动力条件。

扑动水翼是由鱼类游动方式简化而来的，具有阻力小、噪声低、效率高、机动性好等特点。自文献[11]提出扑动水翼可以作为一种新式螺旋桨以来，扑动水翼的推进性能成为了研究热点[12-15]，扑动水翼的应用范围也被不断扩展，比如利用串联水翼从非恒定场中获得推进力的波浪滑翔机[16-17]，以及利用水翼从波浪中获取能量的海洋能量采集装置等[18-19]。文献[20]对比了扑动、摆动以及纯俯仰3种运动形式，发现3种方式在合适的参数下均能在尾流中形成双列反卡门涡街；文献[21]则对不同运动参数下的扑动水翼进行了研究，分析了升沉以及俯仰之间的相位差、攻角的偏置以及运动方程中高次项对扑动水翼推进性能的影响；文献[22]研究了斯特劳哈尔数和最大攻角的变化对推力和流体力学效率的影响，确定了在一定的参数范围内能达成高效率与高推力的共存，并且试验中测得的效率峰值超过了70%。文献[23]对常规正弦运动的4种组合方式进行了模拟分析，通过分析各形式下的尾流结构以及升阻系数，提出俯仰的相位比升沉超前90°时推进效果最佳。文献[24]利用改变扑动水翼的运动形式的方法增加了扑动水翼的平均迎角，有效提升了扑动水翼在较大斯特劳哈尔数St时的推力以及效率。上述研究均表明在合适的运动参数下，扑动水翼能够从尾流的双列反卡门涡街获得强大的推进力，有着优异的水动力性能。

以上对扑动水翼的研究主要集中于扑动水翼的推进性能，关于扑动水翼对水体的推动效果的研究较少。为了研究扑动水翼的推水性能，本文利用CFD方法建立扑动水翼的二维计算模型，对不同频率以及不同来流速度情况下扑动水翼的推水性能进行模拟研究，并搭建扑动水翼试验装置对模拟结果进行验证。

1 运动模型

1.1 运动描述

水翼的研究分为柔性水翼以及刚性水翼两种，其中柔性水翼自由度较多，运动复杂，且其优势在较高扑动频率区域；而刚性水翼的运动更为简单、便于控制、成本低廉，有着更为实际的工程意义[25]，因此本文选择刚性水翼为研究对象。扑动水翼的运动由升沉与俯仰两种运动结合而成，根据文献[23]选择常规正弦运动，如图1所示。图中，Amax表示扑动水翼的升沉幅值；θmax表示扑动水翼的俯仰幅值；T表示运动周期。

图1 扑动水翼的扑动形式示意图Fig.1 Schematic of flapping form of hydrofoil

扑动水翼的基本运动方程为

(1)

式中y(t)——扑动水翼升沉位移，m

θ(t)——扑动水翼俯仰位移，rad

ω——扑动角频率，rad/s

φ——升沉与俯仰运动的相位差，rad

对式(1)求导，得扑动水翼的速度方程为

(2)

1.2 运动参数与力学模型

在扑动水翼的研究中，瞬时推力系数Ct以及瞬时升力系数Cy是衡量扑动水翼水动力性能的关键参数，计算公式分别为

(3)

式中Ft(t)——水平方向瞬时推力，N

Fy(t)——竖直方向瞬时升力，N

ρ——流体密度，kg/m3

s——翼型的展长，m

平均推力系数及平均升力系数定义为

(4)

式中n——计算总周期数

k——取值周期数，取2

除此以外，为表征扑动水翼的推水性能，还需要计算扑动水翼装置的流量、扬程以及效率。

流量公式定义为

(5)

b——流道宽度，取0.8 m

平均扬程可由出入口的压力差换算获得，计算公式为

(6)

g——重力加速度，m/s2

扑动水翼运动的平均输入功率计算公式为

(7)

M(t)——翼型绕转轴的瞬时扭矩，N·m

(8)

水力推进效率ηpr则可以表示为

(9)

式中CP——瞬时功率系数

(10)

由此，可得装置泵效率ηpu计算公式为

(11)

2 数值方法

2.1 控制方程与湍流模型

本文利用商业CFD软件Fluent进行数值计算，考虑黏性二维不可压缩流动的雷诺平均纳维-斯托克斯方程(RANS)和时均连续性方程，可知其运动控制方程为

(12)

(13)

xi、xj——i、j方向空间坐标，m

p——流体压强，Pa

υ——层流运动粘性系数，Pa·s

υt——湍流粘性系数，Pa·s

本文采用Realizablek-ε湍流模型求解N-S方程的模型，相应方程参见文献[26]。

2.2 计算域与网格生成

计算域尺寸以及远场边界条件对计算结果的影响不可忽略，为使翼型后续流场能够充分发展，常取翼型后流域长为20c，故本文中设置x方向总长为8 m，并设置x正方向长度为6 m，翼型旋转中心置于原点位置处。为了模拟河流推水时河道两岸对流场的影响，设置流道宽度为0.8 m。计算域及初始时刻无俯仰角的网格划分如图2所示。

图2 计算域及网格划分Fig.2 Computational domain and grid partition

由于扑动水翼需要进行起伏以及俯仰结合的运动，还需要结合动网格技术实现不同时刻的网格重构。为了提高计算效率，利用非结构化三角形网格划分以翼型起始位置为中心的长1 m、宽0.8 m的长方形变形域，利用结构化网格划分其余计算域。为了更好地捕获翼型以及壁面的边界流场，在其周围划分边界层，第1层网格厚度为0.000 1 m(y+<2.71)。此外，为保持翼型周围网格的质量，将翼型边界层及其向外数层网格单独分为一个固定域，该区域随翼型运动，不参与网格重构。

为验证网格数量的无关性，使用不同网格数量进行数值模拟计算。结果表明，在网格数量达到8.3万后，其对计算结果的影响可以忽略，因此为减少计算的时间，选择网格数为82 994。

2.3 数值计算

本文基于ANSYS Fluent软件，结合Realizablek-ε湍流模型以及动网格技术进行扑动水翼推水性能的数值研究。设置近翼型侧为压力入口，远翼型侧为压力出口，两侧壁面与翼型表面为无滑移壁面，3个计算域的交界面设为Interface，近壁计算采用增强壁面处理。

动网格模型采用弹性光顺法以及局部重构法。对于弹性光顺法，为不影响较远区域的网格，设定弹性系数为0.8，为不影响边界处网格节点分布，设定边界点松弛因子为0.000 6。对于局部重构法，设置最大网格扭曲率为0.7以及合适的最大最小网格长度标准。最后利用用户自定义指令(UDF)进行运动域运动形式的设定，得到不同时刻的动网格如图3所示。求解器设置中，选用Coupled算法耦合压力场与速度场，采用二阶迎风格式来离散时间。而步长设置应小于最小网格尺度与流场速度的比值，需要根据不同工况的计算情况进行调整。

图3 1个周期内4个主要时刻的动网格更新Fig.3 Dynamic mesh updating at four main points in a cycle

2.4 方法验证

为验证数值计算方法的有效性，选取NACA0012翼型，在最大攻角相同的情况下，对St不同的多种工况进行计算，并与文献[21]的试验结果进行对比。基本参数设置与文献[21]相同，U=0.4 m/s，φ=π/2，c=Amax=0.1 m，θmax=π/12，Re=4×104，其中U表示计算域入口的流速。

将扑动水翼的平均推力系数随St变化的数值计算结果与文献[21]试验数据进行对比，结果如图4所示。

图4 数值模拟结果与文献[21]试验数据的对比Fig.4 Comparison of numerical simulation results with experimental data in reference[21]

由图4可知，平均推力系数的计算值与试验值相符合，具有较好的一致性，本文使用的数值计算方法是正确有效的。

3 数值结果与分析

计算模型的各参数设置如下：弦长c为0.3 m，扑动水翼的旋转中心与前缘距离l为0.2c，起伏幅值Amax为0.5c，俯仰幅值θmax为π/6，相位角φ为-π/2，计算时展长s为0.3 m。

3.1 扑动频率对扑动水翼装置水动力性能的影响

为了研究扑动频率对扑动水翼装置水动力性能的影响，本文设置扑动水翼摆动频率从0.1 Hz逐渐增加到1 Hz，每0.1 Hz变化一次；此后每0.5 Hz增加一次，最高频率为5 Hz。从其中选择较为典型的4组进行对比分析。

图5为不同扑动频率下扑动水翼的瞬时推力系数随时间变化的曲线。为了便于对比，将各自周期的无量纲时间的相对值作为横坐标。由图5可以看出，当扑动水翼从平衡位置开始向上扑动，瞬时推力系数先增大再减小，扑动至最大扑动幅值时瞬时推力系数到达谷点，随后扑动水翼开始反向扑动，瞬时推力系数再次开始增加，过平衡位置后迅速达到另一波峰。而且在一个运动周期内，瞬时推力系数总是大于零，即在扑动的全过程扑动水翼总是对水流起到推动作用。

图5 不同频率下扑动水翼的瞬时推力系数变化曲线Fig.5 Variation of instantaneous thrust coefficient of flapping hydrofoil at different frequencies

图6为不同扑动频率下扑动水翼的瞬时升力系数随时间变化的曲线。由图6可以看出，瞬时升力系数曲线基本沿轴线对称分布，所有工况下的平均升力系数均接近零。对比图5以及图6可以发现，在瞬时推力系数的极大值处(扑动平衡位置)瞬时升力系数取极值，而在瞬时推力系数的极小值处(即最大扑动位置处)瞬时升力系数为零。

图6 不同频率下扑动水翼的瞬时升力系数变化曲线Fig.6 Variation of instantaneous lift coefficient of flapping hydrofoil at different frequencies

为分析扑动水翼尾流的结构形态对其性能的影响，选择涡量在z方向分量Ωz为特征量，并由此绘制各工况下的涡量等值线图，z向涡量分量的计算公式为

(14)

式中ux、uy——x、y方向的速度，m/s

图7为不同频率下各时刻的涡量等值线图。可以看出，扑动水翼在有限通道内进行扑动推水时，尾流中的反卡门涡街可能会发生偏移，被压缩侧受到的来自壁面的作用更为强烈，消散更为迅速，导致双列反卡门涡街两侧的旋涡强度差异较大。综合观察图5～7可以发现，瞬时推力系数与瞬时升力系数的大小与反卡门涡街的旋涡强度密切相关。以f=1 Hz为例，该工况下的尾部射流偏向斜上方，导致反卡门涡街下方的旋涡强度明显强于上方旋涡。当翼型下扑时，下方强旋涡对翼型的作用远大于翼型上扑时上方弱旋涡对翼型的作用，导致翼型下扑时对水流的推动力明显强于上扑时的推动力，故图5中f=1 Hz时瞬时推力系数曲线中每周期的第2个波峰远高于第1个波峰。此外，尾涡作用也是导致瞬时推力系数极大值偏离平衡位置的原因之一。

图7 频率为0.1、0.5、1 Hz不同时刻的涡量等值线图Fig.7 Contour maps of vorticity at different frequencies of 0.1 Hz, 0.5 Hz and 1 Hz

图8为扑动水翼装置的平均流量以及平均扬程随着扑动频率变化的曲线。从图8可以看出，扑动水翼装置的流量与其扑动频率基本成正比，而扬程与扑动频率的平方成正比。结合图7以及图8可以发现，在模拟过程中尾流的偏移并不会对流量以及扬程产生明显的影响，即尾涡的偏移虽然增大了推力的瞬时变化，但是对平均推力影响较小。

图8 流量以及扬程随扑动频率的变化曲线Fig.8 Variation of flow rate and lift with flutter frequency

图9为扑动水翼装置的推进效率以及泵水效率随扑动频率变化的曲线。可以看出，扑动水翼的推进效率与泵水效率随着扑动频率的增加不断增加，但是增加的速率随频率增加不断降低。对比相同频率下的推水效率与泵水效率可以发现，推水效率总是高于泵水效率，而且两者之间的差值随着频率的增加不断减小，这也许是由计算段的水力损失所占的比重随着扑动频率的增加有所降低导致的。

图9 推进效率以及泵水效率随扑动频率的变化曲线Fig.9 Variation of propulsion efficiency and pump efficiency with flutter frequency

3.2 扑动水翼装置的推水特性

为了探究不同流量下扑动水翼装置的扬程以及泵水效率的变化规律，选取出入口流速为0～1 m/s，每0.1 m/s变化一次进行模拟计算研究。计算过程中固定扑动频率为1 Hz。图10(图中Q0.7表示来流速度为0.7 m/s时的流量)给出的是以最大效率点流量为流量标度时，扑动水翼装置的扬程以及效率随着流量变化的曲线。

由图10可知，扑动水翼装置的扬程随着流量的增加不断减小，但在流量较小时有着明显的波动；而扑动水翼装置的泵水效率随着流量的增加先增大再减小。随着流量的增大，扑动水翼与水流之间的相对速度不断减小，翼型对水流的推动力越来越小，故装置的扬程随着流量的增加不断减小。除此之外，在最佳效率点时翼型在扑动过程中首次对水流起到了阻碍作用。此外，该装置具有超低的扬程，当流量为0即封闭计算域出入口时获得的平均扬程也仅有0.024 m，在河流推水这种低扬程、大流量的应用场合具有优势。

图10 频率为1 Hz时，扬程与效率随流量的变化曲线Fig.10 Variation of lift and efficiency with flow rate when frequency was 1 Hz

图11为频率1 Hz时，不同流速下3种流场的涡量等值线图。可以发现：在低流速下，尾涡上下并列并且快速消散，这是因为新尾涡的生成过程会对已经形成的旧涡造成影响，甚至是直接破坏来不及流向下游的旧涡，尾涡的能量直接损失，能量使用效率低下。而随着流速的不断增加，每个周期中生成的两个旋涡逐渐分离，排列形成双列反卡门涡街，在涡街中心形成射流加速水流的流动，尾涡能量逐渐得到利用，装置效率不断增加。当流速超过最佳效率点后，尾涡中每个周期中生成的两个反方向旋涡远离的过程中，尾涡由双列反卡门涡街转变为单列反卡门涡街，相互作用减小，无法形成有效射流推动水流的流动，装置效率开始下降。

图11 频率为1 Hz、不同流速时不同时刻的涡量等值线图Fig.11 Vorticity contours at different time with frequency of 1 Hz and different flow velocities

图12为扑动水翼在最佳效率点平衡位置的涡量等值线图，此时扑动水翼装置的尾迹呈现为如图所示的双列反卡门涡街形态。图13为扑动水翼在最佳效率点平衡位置的压力云图，综合观察图12和图13可以发现，尾涡在向后移动的过程中逐渐耗散，尾涡能量中的一部分被转化为压力势能，压力逐渐增加。此外观察其他工况下的涡量等值线图可知，扬程在0.010～0.020 m范围内时，装置尾迹均为相似的双列反卡门涡街，且具有较好的推水性能，能适应扬程波动较大的应用场合。

图12 扑动水翼在最佳效率点平衡位置的涡量等值线图Fig.12 Vorticity contour map of flapping hydrofoil at equilibrium position of optimal efficiency point

图13 扑动水翼在最佳效率点平衡位置的压力云图Fig.13 Pressure contour of flapping hydrofoil at equilibrium position of optimal efficiency point

4 性能试验

为了验证扑动水翼装置的推水性能，设计了扑动水翼试验装置，如图14所示。伺服电机与减速箱提供装置运行所需的动力；曲柄滑块机构将旋转运动转化为直线运动，带动叶片进行升沉运动的同时，配合两组曲柄滑块机构之间的相位差实现叶片的俯仰运动，从而实现叶片的正弦扑动。

图14 扑动水翼装置示意图Fig.14 Schematic of flapping hydrofoil device1.伺服电机与减速箱 2.曲柄滑块机构 3.导轨 4.背板 5.扑动叶片 6、7.支架

图15为扑动水翼装置实物图。试验装置架设在4 m×0.5 m×0.7 m的开放流道上方。通过声学多普勒流速仪测量流道出口处各点的流速，并利用高速摄像机进行流场的捕捉。其中声学多普勒流速仪的量程为0～3 m/s，精度为0.005 m/s，采样频率为25 Hz；高速摄像机选用Phantom VEO 340L，帧速为100 f/s。

图15 扑动水翼装置试验台Fig.15 Flapping hydrofoil device test bench1.扑动水翼装置 2.水槽 3.方形流道 4.叶片 5.伺服驱动器

试验叶片展长为0.3 m，试验频率为0.1～1 Hz，每0.1 Hz变化一次。利用声学多普勒流速仪测量方形流道出口处均匀设置的6个测点上流速，单次测量不小于12个周期，多次测量取其平均值为测点流速。以6个测点的平均流速作为平面流速。将试验流速与仿真结果进行对比，如图16所示。频率为1 Hz时，试验流速为0.264 m/s，总不确定度u为0.020 1 m/s。

图16 试验与仿真的流速对比Fig.16 Comparison of flow velocity between test and simulation

由图16可知，试验与模拟的结果整体趋势相同，装置的流量随着运动频率的增加而增大。但模拟结果明显高于试验结果，造成该差异的主要原因是试验为保证叶片运动顺利进行，在叶片两侧留有较大的空隙，导致试验过程中部分水体通过两侧缝隙从压力侧转移至吸力侧，大大降低了装置对水体的推动作用。此外，叶片固定块以及泄漏也会对结果造成一定影响。

利用高速摄像机拍摄试验流场结果如图17所示。由于油墨的扩散较快，每个时刻仅有一个尾涡存留，可以发现试验流场与仿真流场基本相同，从侧面印证了仿真结果的正确性。

图17 试验与仿真流场图Fig.17 Flow field diagrams of test and simulation

5 结论

(1)扑动水翼装置的尾涡结构形态是性能影响的关键因素，当尾涡处于双列反卡门涡街时，推水性能较高。

(2)在不同频率下，扑动水翼装置的流量与扑动频率成正比，扬程与扑动频率的平方成正比。而装置的推进效率以及泵水效率均随着频率的增加而增大，但增加的速率不断降低。

(3)在固定的频率(1 Hz)下，扑动水翼装置的扬程随着频率的增加而减小，而效率随着流量的增加先增大后减小。

(4)扑动水翼装置具有超低的扬程，当频率为1 Hz时，该装置的最大平均扬程仅有0.024 m；而且扬程在0.010～0.020 m范围内，该装置尾流均呈现双列反卡门涡街，具有较好的推水性能，在河流推水这种低扬程、大流量的应用场合具有优势。