D形柱对摆动水翼水动力性能影响分析

舒蒋鹏, 石小涛, 陈小龙, 周柳青, 王思莹, 胡晓, 张立胜*

(1.三峡大学 湖北省鱼类过坝技术国际科技合作基地, 湖北 宜昌 443002;2.三峡大学 水利与环境学院, 湖北 宜昌 443002;3.武汉理工大学 新材料力学理论与应用湖北省重点实验室, 湖北 武汉 430063)

0 引言

自然界中的鱼类经过数十亿年的自然选择,已经具备了高机动性、低功耗的游动能力[1],其游动能力是仿鱼类水下机器人所无法比拟的。为了探究鱼类高效的游动性能,许多学者对此展开了大量研究。前期研究主要通过理论和实验的方法对仿生尾鳍开展研究。Lighthill[2]研究了尾鳍的摆动,并提出了大摆幅细长体理论。Freymuth[3]通过烟线法研究摆动水翼在风洞中的流场变化。Triantafyllou等[4]研究发现,当摆动水翼在斯特劳哈尔数St为0.25～0.35时,其推进效率最大。Anderson等[5]和Read等[6]通过实验研究发现,摆动水翼在特定频率下和特定攻角下的推进效率最高。Lauder等[7]用数字粒子图像测速(digital particle image velocimetry, DPIV)技术测量了摆动水翼的尾部流场,分析了摆动水翼运动参数对尾涡结构的影响。王肇等[8-9]通过水洞试验研究了摆动水翼的水动力性能。Schouveiler等[10]通过实验研究Sr和攻角对摆动水翼推进性能的影响。Iverson等[11]在循环水道内对摆动水翼的力进行直接测量,并结合粒子图像测速(particle image velocimetry, PIV)技术对水翼尾流的涡结构进行评估。近些年,计算流体力学(computational fluid dynamics,CFD)方法的发展弥补了理论和实验的不足,可以对流场结构变化、游动机制以及推进性能等进一步深入研究。张晓庆等[12]采用数值方法模拟了二维刚性和柔性水翼并与实验值进行对比,结果较好。杨亮[13]采用雷诺平均的NS方程(reynolds averaged navier strokes, RANS)求解器和动网格方法研究了运动参数对摆动水翼水动力性能的影响。文敏华等[14]基于动网格技术分析了摆动水翼推力中黏性力和压差力的占比情况。刘焕兴[15]采用RANS求解器研究了Sr和最大攻角对摆动水翼推进性能的影响。

上述研究均将摆动水翼放置在均匀流场中,而鱼类真实的生存环境常常存在各种漩涡,其流场较为复杂,并且摆动水翼作为鱼类尾鳍的一种简化模型,研究流场中涡对摆动水翼的推进性能影响至关重要。另外有研究指出,鱼类可以从漩涡中吸收能量,提高其推进性能。为此本文中基于CFD,利用Fluent的动网格技术和用户自定义函数(user define functions,UDF),通过在摆动水翼前方放置一个D形柱模拟复杂流场,分析摆动水翼的升沉振幅和摆动频率对摆动水翼推进性能的影响,并与均匀流场中的摆动水翼对比,所得结果可为仿鱼类水下机器人的研制提供理论参考。

1 数值计算模型

1.1 数值计算模型

本文定义2种运动模式:模式1为在均匀流场下的摆动水翼俯仰-升沉耦合运动;模式2为D形柱尾流场(复杂流场)下的摆动水翼俯仰-升沉耦合运动。摆动水翼前缘与D形柱圆心距离S=4d,d表示D形柱直径,d=c/2(c为弦长)。模式2示意图如图1所示。二维摆动水翼采用NACA0012水翼,枢轴点为c/3处,其满足的运动规律[16]可以表示为

图1 模式2示意图Fig.1 Schematic of mode 2

y(t)=h0sin(2πft),

(1)

θ(t)=θ0sin(2πft+φ),

(2)

式中:y(t)为升沉速度;θ(t)为俯仰速度;h0为升沉运动振幅;θ0为俯仰运动振幅;φ表示耦合运动相位差;f为水翼运动频率;t为时间。有研究表明摆动水翼在φ=-90°时推进性能最优,故本研究φ值恒为-90°。D形柱和摆动水翼的周期性运动图如图2所示。

图2 D形柱和摆动水翼的周期性运动图Fig.2 Periodic motion diagram of D-cylinder and flapping hydrofoil

1.2 水动力系数

对于摆动水翼运动,定义水翼斯特劳哈尔数St=2fh0/U,式中U为来流速度[17]。水翼的瞬时推力系数Ct、瞬时升力系数Cy和瞬时力矩系数Cm分别为

(3)

式中:Fx(t)、Fy(t)、Mz(t)分别是摆动水翼的瞬时推力、瞬时升力、瞬时力矩;ρ为水的密度。

相应地,水翼一个运动周期内的平均推力系数CT、平均升力系数CL、平均力矩系数CM分别为

(4)

式中T为整个运动周期。

水翼的推进效率η为[18]

(5)

式中CP为一个周期内的水翼平均输入功率系数,其表达式为

(6)

(7)

2 数值计算方法

2.1 控制方程

本文采用二维不可压缩黏性纳维-斯托克斯(N-S)方程作为流体的控制方程,可以表示为

·u=0,

(8)

(9)

式中:u是流体运动速度;p是流体压力;ν是流体动力学黏性系数;湍流模型选择SSTk-ω两方程模型。

2.2 计算域网格划分及参数设定

计算域网格划分及参数设定如图3所示。计算域的尺寸为10c×8c(长度×宽度)。计算域左侧为速度入口,右侧为压力出口,上下为对称边界,摆动水翼边界设为无滑移壁面。为了提高计算效率,采用结构化网格跟非结构化网格混合的方法进行网格划分。计算域分为左、右域,交界处设为Interface,又将左、右2个流体区域分为内、外域,内、外域交界处设为Interior,内域分别为半径为c和半径为2c的圆形区域,采用非结构化四边形网格划分,整个计算域除去内部计算域的区域为外域,采用结构化三角形网格划分,并采用弹簧光顺和网格重构法对畸变较大的网格进行重新划分以确保网格的合格性。

图3 计算域网格划分及参数设定Fig.3 Computing domain grid division and parameter setting

采用SIMPLE算法对连续方程中的压力和速度进行耦合,为了提高计算精度,离散方式中选择时间项采用一阶隐式,动量项采用二阶迎风格式,压力项选择二阶格式。

2.3 网格无关性验证

为了验证网格的无关性,以模式1为验证算例进行网格划分,划分3种不同数量的网格,分别为8.95×104、12.10×104、16.10×104。时间步长选择Δt=T/200,运动参数选择,h0/c=1.0,St=0.25,θ0=25°,f=0.5 Hz,φ=-90°,不同网格数量时的CT值和相对误差见表1。从表1中可以看出,网格B的相对误差最小,因此为了节省计算时间,本文后续计算采用网格B。

表1 网格无关性验证Tab.1 Mesh agnostic validation

2.4 模型有效性验证

为了验证数值方法的有效性,数值计算了斯特劳哈尔数St分别是0.15、0.20、0.25、0.30、0.35、0.40,最大攻角amax=15°,相位差φ=-90°,h0/c=1.0,雷诺数Re=40 000下的CT值,并将计算值与文献[6]实验结果进行比较,结果如图4所示。由图4可见,平均推力系数CT的计算值与实验值吻合良好,具有相同的变化趋势,说明该数值方法是有效的。

图4 平均推力系数CT计算值与文献[6]结果的对比Fig.4 Comparison between the calculated values in this paper and those in Reference[6]

3 计算结果与分析

本文主要研究摆动频率f和升沉振幅h0对2种模式下摆动水翼的推进性能影响以及流场变化。

3.1 摆动频率对摆动水翼推进性能的影响

设置运动参数:f分别取0.5c、1.0c、1.5c、2.0c、2.5cHz,h0=0.5c,θ0=30°,U=2.0c。

2种模式下的摆动水翼的平均推力系数CT和增长率δ如图5所示。由图5可以看出,2种模式下摆动水翼的CT均随着f的增加而增大,并且模式2在任何f下始终最大,说明D形柱产生的尾流场可以提升摆动水翼的推进性能。此外,从图5可以看出,随着f的增加,δ呈先减小后增大再减小的趋势,f=0.5cHz时,δ=37.2%,f=1.0cHz时,δ=2.13%,研究发现,与增大f所带来的对摆动水翼推力改善相比,D形柱引起的尾流场在摆动水翼推力改善中所占的百分比先减小后增大,然后再减小。

图5 摆动水翼的平均推力系数CT及增长率δFig.5 Average thrust coefficient CT and growth rate δ of flapping hydrofoil

2种模式下的摆动水翼的平均输入功率系数CP和推进效率η如图6所示。从图6可以看出,摆动水翼的CP随着f的增加而增加,而η随着f的增加而减小。当f<1.5cHz时,2种模式下的CP值没有明显区别;当f>1.5cHz后,两者差值逐渐变大,说明D形柱的存在能够增加摆动水翼的功率消耗。此外,2种模式下摆动水翼的η随着f的增加而减小,在相同f下,模式2的η始终比模式1的大;在f=0.5cHz时,模式2的η比模式1的η大29.5%;在f=2.5cHz时,模式2的η比模式1的η大4.9%,说明D形柱的存在可以提升摆动水翼的推进效率,但增加效果是逐渐降低。

图6 摆动水翼的平均输入功率系数CP和推进效率ηFig.6 Average input power coefficient CP and propulsive efficiency η of flapping hydrofoil

3.2 升沉振幅对摆动水翼推进性能的影响

设置运动参数:h0分别取0.5c、1.0c、1.5c、2.0c、2.5c,h0=d,f=1.0 Hz,θ0=30°,U=2.0c。

2种模式下的摆动水翼的平均推力系数CT和增长率δ如图7所示。从图7中可以看出,2种模式下摆动水翼的CT的变化特征,与3.1节中的变化特征一致,不再进行叙述。此外,从图7可以看出,随着h0的增加,δ呈先增大后减小再增大的趋势,在h0/c=0.5时,有最小增量δ=2.1%,在h0/c=2.5时,有最大增量δ=8.2%,表明与增大h0对摆动水翼推力改善相比,D形柱引起的尾流场在摆动水翼推力改善的增长率是先增大后减小再增大,与3.1节中的结果正好相反。

图7 摆动水翼的平均推力系数CT及增长率δFig.7 Average thrust coefficient CTof flapping hydrofoil

2种模式下的摆动水翼的平均输入功率系数CP和推进效率η如图8所示。从图8可以看出,摆动水翼的CP均随着h0的增加而增加,与3.1节的变化规律一致,但值得注意的是,在h0>1.5c时,模式1的CP值均在模式2的上方,说明D形柱的存在能够减少摆动水翼的功率消耗。从图8还可以看出,2种模式下摆动水翼的η随着h0的增加而减小,在相同f下,模式2的η始终比模式1的大;在h0/c=0.5时,模式2的η比模式1的η大0.9%;在h0/c=2.0时,模式2的η比模式1的η大26.1%;在h0/c=2.5时,模式2的η比模式1的η大15.4%,表明D形柱的存在可以提升摆动水翼的推进效率,总体呈先增大后减小的趋势,在h0/c=2.0时增加最明显。

图8 摆动水翼的平均输入功率系数CP和推进效率ηFig.8 Average input power coefficient CP and propulsive efficiency η of flapping hydrofoil

3.3 摆动水翼的流场变化

分析对比2种模式下从D形柱脱落于流场中的涡对摆动水翼流场结构的影响。考虑参数取值分别为f=1.0cHz,h0=1.5c,θ0=30°,U=2.0c。

2种模式下,摆动水翼在一个周期下的流场压力分布云图如图9所示。由图9可见,当t/T=0时,此时摆动水翼的摆角最大,摆动水翼的上、下表面均覆盖高压区和低压区,值得注意的是,在模式1和模式2中摆动水翼摆角为最大时(t/T=0),摆动水翼下表面靠近前缘的区域存在一个明显凸起的强度较大的低压区,并且模式2的区域比模式1的更大。随着摆动水翼运动过程中摆角的减小,这个凸起的低压区向摆动水翼的尾缘移动,并在摆动水翼摆角为0时(t/T=1/4)从摆动水翼下表面分离。由于摆动水翼在一个运动周期的对称性,因此在t/T=1/2时,摆动水翼的上表面靠近前缘的区域也存在一个明显凸起的低压区,此时摆动水翼上表面为低压区,下表面为高压区,与摆动水翼在t/T=0时的情况相反。模式1和模式2中摆动水翼上、下两侧的压力分布区域一致,不同之处在于模式2中压力的幅值更大。摆动水翼两侧压力差的变化决定摆动水翼的推进性能,又因为2种模式下摆动水翼的摆角相同,所以模式1和模式2的摆动水翼的平均推力系数逐渐增大。

(a) 模式1

(b) 模式2图9 流场压力分布云图Fig.9 Flow field pressure distribution cloud

2种模式下,摆动水翼在一个周期下的摆动水翼尾流场的涡结构云图如图10所示。设ωz为流场涡量在z方向的分量,其计算公式[19]为

(a) 模式1

(b) 模式2图10 摆动水翼尾流场的涡结构云图Fig.10 Vortex structure cloud image of the wake field of a flapping hydrofoil

(10)

式中:ux、uy分别表示流场在x轴和y轴方向的速度分量。

由图10可知,摆动水翼在2种运动模式下的尾流场均为反卡门涡街形式。对于模式1,当摆动水翼处于最大摆角时(t/T=0),尾缘处存在一个尾缘涡M-1(顺时针),并与旋向相同的附着涡LM-1相连接,上端存在一个FM-1涡,与水翼的下表面前缘涡TM-1旋向相同。当水翼摆动至平衡位置时(t/T=1/4),此时摆角为0,TM-1和M-1不断扩大,但M-1并未脱落。此后,摆动水翼反向运动到最大摆角时(t/T=1/2),M-1脱落至尾流场中,形成一个高强度的涡核,TM-1不断向尾缘移动,并与LM-2相连接,此刻,水翼上、下表面存在旋向相同的TM-2和FM-2涡。当水翼反向摆动至平衡位置时(t/T=3/4),TM-1和TM-2不断扩大,但并未完全从尾缘脱落,最后水翼运动至初始位置,TM-1发生脱落,与M-1形成旋向相反的涡对。

对于模式2,摆动水翼前方存在D形柱,D形柱不断脱离2个旋向相反的涡1和涡2,摆动水翼的摆角最大时(t/T=0),D形柱产生的顺时针涡1靠近水翼前缘,此时水翼尾缘处存在一个尾缘涡M-1(顺时针),上端存在一个FM-1涡,与水翼的下表面前缘涡TM-1旋向相同。当水翼摆动至平衡位置时(t/T=1/4),TM-1不断扩大,此时D形柱产生的顺时针涡1与TM-1相互作用,在D形柱产生的顺时针涡1的挤压下,相同时刻下,TM-1涡明显大于模式1的,M-1在尾缘的带动下也逐渐拉长扩大,但并未完全从尾缘脱落。此后,摆动水翼反向运动到最大摆角时(t/T=1/2),D形柱产生的逆时针2涡靠近水翼前缘,TM-1不断扩大,摆动水翼上表面存在逆时针涡TM-2,水翼下表面存在一个FM-2涡,此时M-1完全脱落于尾流场中,形成一个强度高度集中的涡核,与在相同时刻下的模式1中的M-1相比耗散得更快。当水翼反向摆动至平衡位置时(t/T=3/4),此时摆角为0,TM-2不断扩大,TM-1在尾缘的带动下也逐渐拉长扩大,但并未完全从尾缘脱落,此时D形柱产生的逆时针涡2与TM-2相互作用,在D形柱产生的顺时针涡2的挤压下,相同时刻时,TM-2涡明显大于模式1的,最后水翼运动至初始位置,TM-1发生脱落,与M-1形成旋向相反的涡对。D形柱产生的涡1和涡2对摆动水翼的前缘涡TM-1和TM-2产生挤压作用,使得摆动水翼的尾涡区域增大,脱落的涡强度更大,涡耗散得更快,涡对之间诱导作用更强,产生更大的推力。

4 结论

本文基于CFD方法并结合Fluent的动网格技术,通过改变摆动水翼的摆动频率和升沉振幅计算摆动水翼在均匀流场下进行俯仰-升沉运动(模式1)和摆动水翼在D形柱尾流场中进行俯仰-升沉运动(模式2)中的推进性能,并结合流场的压力分布云图和涡结构变化分析了D形柱对摆动水翼的水动力性能影响机制,得到了以下结论:

① 2种模式下,摆动水翼的平均推力系数和平均输入功率系数均随着摆动频率和升沉振幅的增加而不断增加,而推进效率则相反。

② 2种模式下,随着摆动频率的增加,D形柱尾流场中的摆动水翼的平均推力系数始终大于均匀流场中摆动水翼的平均推力系数,并且增长率是先减小后增大再减小。

③ 2种模式下,随着升沉振幅的增加,D形柱尾流场中的摆动水翼的平均推力系数始终大于均匀流场中摆动水翼的平均推力系数,并且增长率是先增大后减小再增大。

④ 2种模式下,摆动水翼尾流场均为反卡门涡街形式,但D形柱尾流场中,D形柱脱落的涡对水翼前缘涡发生挤压作用,使摆动水翼的尾涡区域增大,脱落的涡耗散得更快,涡对之间诱导作用更强,产生更大的推力。