基于直角网格法的单个和阵列布置下柔性水翼绕流数值模拟*

辛建建 陈振雷† 石凡 石伏龙

1) (宁波大学海运学院, 宁波 315211)

2) (长沙理工大学水利工程学院, 长沙 410114)

研究柔性水翼在不可压缩流体中的水动力特性, 对于船舵和减摇鳍等海洋结构物的设计和性能优化具有重要意义.本文将自主开发的径向基函数虚拟网格法求解器扩展到模拟绕单个或多个柔性水翼的不可压缩流动问题.数值模型基于虚拟网格有限差分法考虑浸入边界对流场的影响, 引入紧支径向基函数(compact supported radial basis function, CSRBF)以物面 Lagrangian 质点追踪复杂的柔性动边界.基于该方法, 首先模拟了均匀流中主动拍动的柔性水翼, 升阻力系数良好的网格收敛性结果验证了本文方法的精度和可靠性.并研究了柔性水翼在不同振荡频率下的水动力特性, 阐述了柔性水翼的推力生成机制.然后模拟了绕阵列布置柔性水翼的流动现象, 研究了不同间距和不同振荡频率下水翼表面的升阻力系数变化规律和尾涡特性, 观察到紧密布置的柔性水翼在高频振荡下推力系数存在显著的放大效应, 同时推力为零时的临界频率提前.

1 引 言

鱼类游动是学者和工程师开发先进高效推进设备的重要灵感来源.例如, 仿生机器鱼以类似游鱼摆动尾鳍的方式获得推力, 相比传统螺旋桨推进器, 具有机动性强、噪声低及隐蔽性好等优点, 在海洋勘测和水下侦查等民用和军用等领域具有重要的应用价值.另外, 许多鱼类趋向于成群游动,通过从邻近同伴产生的涡流中获得能量, 以提高游动效率和节省体力.因此, 研究游鱼或鱼群在水中的游动特性, 阐明鱼尾鳍横向摆动的推力产生机制和水动力相互作用机制, 有助于船舵、减摇鳍和仿生推进器等海洋结构物的设计和性能优化, 并为小型自主潜航器或无人机群的最优排布提供技术指导.然而, 游鱼或鱼群游动是一个复杂的水动力现象, 涉及周期性的旋涡脱落, 涡流与鱼体之间的能量交换以及游鱼之间的水动力相互干涉作用.

随着计算机硬件和数值算法的迅速发展, 计算流体力学 (computational fluid dynamics, CFD)技术为鱼类等柔性翼型结构物的水动力特性研究提供了强有力的手段.然而, 数值模拟此流固耦合现象对于CFD研究是极具挑战性的, 由于涉及复杂几何形状、界面大变形等难点.典型的求解技术是贴体网格方法[1], 界面附近网格随物体一起运动以精确地描述运动界面, 适合于模拟高雷诺数流动.然而, 贴体网格生成对经验要求高, 拟合柔性多体运动边界涉及复杂的网格生成和重构过程, 也会极大增加计算量和引入额外数值误差.另一方面, 浸入边界法[2]在近些年受到极大的关注.浸入边界法在固定背景(一般为直角)网格上求解Navier-Stokes方程, 通过在流体动量方程中附加力源项[3]或重构浸入边界附近插值模板[4]以考虑浸入边界对流场的影响.其网格生成简单, 无需网格重构, 适于模拟复杂几何形状、多体动边界流动问题.例如, Huang 等[5]采用惩罚浸入边界法,Tullio 和Pascazio[6]采用直接力浸入边界法研究柔性丝带的流致振荡特性, 吴晓笛等[7]采用基于浸入边界的格子玻尔兹曼求解器模拟了双圆柱两自由度涡激振动现象.

针对鱼类游动的水动力机理研究, 许多学者基于直角网格方法对绕单个或多个鱼形水翼的流动问题进行了数值模拟, 并取得卓有成效的进展.Yeo等[8]基于混合直角网格方法研究了单个鱼体的自主推进性能和鱼体大变形时的操纵性, 他们展示了仿生鲤鱼在体长范围内执行70°急转弯的能力.Tian等[9]基于贴体动网格技术的空间时间有限元法模拟了一对母子鱼在均匀流中的振荡特性,研究了鱼体间相对距离、雷诺数、相位差和相速度对力系数、能效和涡结构的影响.Bergmann等[10]基于虚拟区域法模拟了单个和并列布置两鱼体的游动特性.Khalid 等[11,12]采用 Mittal等[13]的多维虚拟网格法研究了串列布置下非同步振荡的两游鱼的水动力性能, 他们发现不同振荡频率的两游鱼产生非线性相互作用, 生成了新的水动力分量, 而且, 在下游鱼体的导边处尾涡会发生分裂, 导致上游鱼体的阻力降低.国内方面, 王亮和吴锤结[14]研究了三鱼群游动以揭示“槽道效应”在鱼群游动中的节能机制, 王文全[15]研究了鱼体自主游动的尾涡模式和鱼体运动特性.

目前, 学者对单鱼体、串列或并列布置的双鱼体的流体特性进行了大量研究.然而, 尚未知串列或并列布置时形成的力模式和尾涡结构是否适用于鱼群游动, 当并列和串列布置同时存在时是否有新的力模式和尾涡结构出现.因此, 研究流体中以阵列布置的鱼群游动能更好地再现鱼群游动的水动力特性.鉴于此, 本文拟数值研究单个或阵列布置下鱼形水翼在流场中的水动力特性, 以揭示鱼体游动的推力产生机制和鱼群游动的水动力优势.数值模型基于自主开发的径向基函数虚拟网格法[16,17],在固定直角交错网格上求解不可压缩Navier-Stokes方程, 引入紧支径向基函数到虚拟网格法以处理任意复杂的柔性边界, 一个面积分数表示方法通过提高浸入界面的局部质量守恒性以抑制动边界的压力振荡问题.采用该方法模拟了不同激励频率下以鱼形游动的水翼在均匀流中的水动力特性,并进行了网格收敛性的验证分析.然后, 研究了9个以阵列布置的鱼形水翼在不同间距和不同激励频率下的力系数变化规律和尾涡形态.

2 数学模型和数值方法

2.1 控制方程和数值离散

非定常黏性不可压缩流体的控制方程如下:

其中 u 是速度向量, 具有直角网格分量 (u,v) , ρ 是密度, t 是时间, p 是压力, ν 是流体的运动黏性系数, V和a是固体的速度和加速度; n是固体表面的外法向向量; Γ 表示固体表面.

流体控制方程(1)和方程(2)采用有限差分法在交错直角网格上离散, 采用虚拟网格法施加浸入动边界的无滑移边界条件(3)式.在本文求解器中,以分步法解耦速度和压力, 以二阶Total Variation Diminishing Runge-Kutta(TVD-RK2)格式进行时间积分, 采用时间显式格式处理对流项, 黏性项以 Crank–Nicolson 格式半隐式处理.对于空间离散格式, 采用TVD-MUSCL(Monotonic Upstream-centered Scheme for Conservation Laws)[18]格式离散对流项, 黏性项以中心差分格式离散.采用不完全Cholesksy预处理的双共轭梯度法 (Incomplete Cholesksy Biconjugate gradient stabilized method, BiCGSTAB)[19]求解压力泊松方程离散形成的大型稀疏线性方程组.

2.2 虚拟网格法

为了考虑静止或运动物体对流场的影响, 本文采用虚拟网格法.该方法的基本理念是, 通过在浸入物体内部布置有限数量的虚拟网格以表示静止或运动边界的存在, 以适当的插值模块及技术(如双线性插值方法), 施加浸入动边界条件.

图1 虚拟网格法的二维插值格式 (a)无虚拟网格; (b) 一个虚拟网格; (c) 两个虚拟网格Fig.1.Two-dimensional interpolation stencil of the ghost cell method: (a) No ghost cells; (b) one ghost cell; (c) two ghost cells.

虚拟网格的插值重构过程可细分为三步: 标记镜像点、插值重构和满足镜像条件.首先, 对于每个虚拟网格(GC), 沿着物体表面的外法向从虚拟网格单元中心延伸到流体域以确定唯一的一个镜像点(MP), 如图1所示.然后是采用双线性插值格式计算每个镜像点的流体变量, 为了避免内插模板点是虚拟网格自身的情况, 镜像点的计算需要分如图1中的三种情况进行讨论.最后是根据镜像点变量通过满足镜像条件得到虚拟网格变量值, 更多细节参见文献[17].需要注意的是由于交错网格布置, 两个方向速度分量或压力的离散单元位置并不重合, 因此, 需要对每个方向的速度分量或压力分别应用以上的重构过程.在整个计算域(液体和固体)求解Navier-Stokes方程, 边界条件便得到隐式满足.

2.3 CSRBF界面追踪方法

浸入边界法的网格生成过程简单, 但由于浸入边界本身并未包含在计算网格中, 精确表示背景网格中浸入边界的存在是个难点.本文采用基于CSRBF的隐式等值面方法重构任意界面, 以光滑的基函数拟合离散的物面控制点构造任意复杂物体的等值曲面, 根据等值面的特性识别场点属性.本文CSRBF界面描述方法简单、以有限的物面控制点获得较高的精度, 适于描述不规则、大变形甚至凹形界面.

针对任意浸入界面, 基于有限数量的界面位置信息的离散控制点, 引入径向基函数构造如下场函数:

其中M为物面控制点的数量, βj为权重系数,P(xi)为线性多项式.在文献 [16]中, 我们采用MQ (Multi-quadrics)核函数拟合复杂刚性界面.然而, 当涉及复杂薄壁边界时, 由于矩阵病态, 拟合过程可能失败.为了解决该问题, 本研究采用紧支 Wendland’s C2 函数[20], 得到的矩阵为正定稀疏矩阵, 矩阵的病态性得到缓解, 收敛性得到提高.Wendland’s C2 函数表达如下:

3 数值结果和讨论

本文径向基函数虚拟网格法的计算程序采用Fortran 90语言编写, 数值算法的流程图见图2.以圆柱绕流、翼型绕流和圆柱振荡等算例验证了本文方法的精度、可靠性和模拟刚性动边界的能力[16,17].

3.1 鱼形水翼在均匀流中的振荡运动

这里模拟了均匀流中柔性水翼的横向振荡运动以验证本文方法模拟柔性动边界流动的能力.水翼剖面为NACA0012机翼, 采用径向基函数拟合141个样品点进行追踪.为了模拟水中游鱼的运动,柔性水翼的横向振荡采用一个波浪函数描述:

图2 静止或运动边界流动求解算法的流程图Fig.2.Flowchart of numerical algorithm for stationary or moving boundary flows.

y(x,t)是水翼表面质点在y方向的位移, ω 是振荡角频率, 波数 k =2π .振荡幅值 A (x) 表示为A(x)=A1+A2x+A3x2, 其中系数分别为 A1=0.04 , A2=−0.16和 A3=0.32 .因此, 水翼运动在导边和随边的最大幅值分别为 0.04 和 0.2 m.雷诺数 (Re)和Strouhal数(St)是两个重要的流动参数, Strouhal数定义为 St=2Amaxf/U∞, 其中 Amax是随边的最大幅值, f 是振荡幅值 (f = ω/2π).根据 Sui等[21]的工作, 雷诺数设置为 Re = 400, 选择了五个振荡角频率 (ω = 0.15, 0.30, 0.45, 0.60 和 0.75)以研究振荡频率对力系数的影响, 分别对应的Strouhal数 为 St= 0.191, 0.382, 0.573, 0.764 和0.955.计算区域设置为 [–3b, 9b] × [–3b, 3b], 翼弦长度 b = 1 m.入口边界为均匀流, 均匀流速度 U∞=0.05 m/s, 在出口采用出流边界条件, 上下边界采用辐射边界条件.对于压力, 在出口指定压力为零,其余边界采用Neumann边界条件, 选择固定的时间步长 ∆ t = 0.002 s.

首先对激励频率 ω = 0.45时鱼体振荡运动进行了网格收敛性研究, 采用非均匀网格对局部区域 [–0.2b, 1.3b] × [–0.3b, 0.3b]进行网格细化, 细化区域分别采用 150 × 80, 300 × 150, 600 × 300,对应在x-方向的最小网格间距 ∆ x 分别为0.01b,0.005b和0.0025b, 对应在y-方向的最小网格间距∆y分别为 0.0075b, 0.004b 和 0.002b.表1 比较了本文三套网格下计算得到的力系数与Sui等[21]的参考结果.从表1知, 随着网格细化, 平均阻力系数和均方根升力系数逐渐趋近Sui等[21]的Lattice Boltzmann方法结果, 而且, 两套细网格上的结果基本重合, 说明本文方法获得收敛解, 其中中间网格(300 × 150)上计算的平均阻力系数相比参考结果略大2.22%, 均方根升力系数略小4.17%.

鉴于300 × 150的中间网格上已计算得到相对合理的收敛解, 将300 × 150的中间网格用于其他工况的计算.图3将五个频率下的力系数与Sui等[21]的Lattice Boltzmann方法的计算结果进行了比较.平均阻力系数和均方根升力系数与参考结果合理地符合.对于角频率 ω ≤0.45 , 本文力系数与参考结果较好地符合, 当 ω >0.45 , 可观察到力系数与参考结果的微小偏差.以角频率ω=0.75为例, 本文平均阻力系数略大于参考结果7.1%, 均方根升力系数略小于参考结果3.3%.从图3可观察到, 随着振荡频率的增加, 平均阻力系数迅速降低.临界频率为 ω =0.6 , 此时平均阻力系数从正值变为负值, 即水翼在流体中运动所遭受的阻力变为推力, 这是由于高频振荡产生的推力即为水中鱼类向前运动的动力来源.当角频率从0.6增加到0.75时, 推力系数从大约为0显著增加到0.21.另外, 当角频率从0.15增加到0.75时, 均方根升力系数呈现近乎线性的增加趋势(从0.21 到 3.1).

表1 本文力系数结果与 Sui等[21]的参考结果比较Table 1.Comparison of the force coefficients between the present method and Sui’s method[21].

图3 五 个 频 率 下 的 力 系 数 比 较 (a) 平 均 阻 力 系 数 ;(b) 均方根升力系数Fig.3.Comparison of force coefficients under five oscillation frequencies: (a) Average drag coefficient; (b) root mean square lift coefficient.

图4展示了三个角频率下 (ω =0.15 , 0.45,0.75)升阻力系数随无因次时间的变化趋势.初始阶段经历一些不稳定的波动后, 三个频率下的升力和阻力系数都呈现出周期性的变化趋势.升力系数的振荡周期是阻力系数的两倍, 该特性与圆柱绕流的卡门涡街现象类似.随着振荡频率从 ω =0.15 增加到 ω =0.75 , 平均阻力系数的平衡位置从0.36降低到–0.3左右, 但是升力系数的振荡幅值从0.4显著增加到4.3左右, 与图3的平均阻力系数和均方根升力系数的变化趋势一致.对于 ω =0.15 ,阻力和升力系数呈现小幅值和规则的正弦曲线变化.当频率增加到 ω =0.45 , 阻力系数曲线的波峰变得略微平坦, 波谷变得较为尖锐, 同时, 升力系数曲线的正弦曲线特性发生轻微偏斜.当频率进一步增加到 ω =0.75 , 阻力系数观察到更加平坦的波峰和更加尖锐的波谷, 类似于脉冲形状, 而且, 升力系数曲线观察到更明显的波峰波谷偏斜现象.

图4 三个频率下升阻力系数随无因次时间的变化趋势(a) 阻力系数; (b) 升力系数Fig.4.Variation of drag and lift coefficients with the dimensionless time under three oscillation frequencies: (a)Drag coefficient; (b) lift coefficient.

进一步分析了柔性水翼振荡的水动力特性.图5—图7分别展示了在三个振荡频率下相位角θ= π/2和 3 π/2 时的局部瞬时涡等值线, 根据(5)式,相位角 θ = π/2 和 3 π/2 分别对应随边位移处于正的和负的最大位移时.由于柔性水翼周期性的摆动,沿着水翼表面的剪切流变得不稳定, 一对涡从水翼表面交替脱落并沿着流动方向流向尾流远方.当水翼尾部在前半个周期摆动到正的最大位移时, 一个顺时针旋转的涡从水翼上表面生成.当水翼尾部在后半个周期摆动到负的最大位移时, 一个逆时针旋转的涡从水翼下表面脱落.这对从水翼尾部交替脱落的反向旋转的涡列就是卡门涡街, 因而产生周期性的升力系数变化.在这个涡列的尾流形成一个能量缺损区, 由此产生一个正向阻力.随着振荡频率的 增 加, 观 察 到 ω =0.45 时尾涡 变 换为 逆 von-Karman涡街形态, 此时, 在水翼尾流区产生一个能量增强区, 阻力降低, 当振荡频率增加到ω=0.6时, 阻力甚至降低到负值, 即为推力, 如图3所示.在更高的振荡频率如 ω =0.75 时, 水翼快速的拍动运动导致逆卡门涡街过渡为偏斜尾涡, 如图5所示.另外, 随着振荡频率增加, 附着在水翼尾部的涡列随着水翼的加速运动以更快的速度脱落, 如图7 所示.在 ω =0.15 , 0.45, 0.75 时分别可观察到两对、四对和七对涡列, 当水翼以高频振荡时(ω=0.75), 迅速脱落的涡列在尾流区紧致排列, 而且,涡之间发生相互排挤, 导致涡形状的变形, 在尾流远方, 涡列甚至发生偏斜.

图5 ω = 0.15 时在两个时刻的局部瞬时涡等值线 (a) θ=/2; (b) θ =3π/2Fig.5.Local instantaneous vortex contours for ω = 0.15 at wo moments: (a) θ = π/2 ; (b) θ =3π/2 .

图6 ω = 0.45 时在两个时刻的局部瞬时涡等值线 (a) θ=π/2; (b) θ =3π/2Fig.6.Local instantaneous vortex contours for ω = 0.45 at two moments: (a) θ = π/2 ; (b) θ =3π/2 .

图7 ω = 0.75 时在两个时刻的局部瞬时涡等值线 (a) θ=π/2; (b) θ =3π/2Fig.7.Local instantaneous vortex contours for ω = 0.75 at two moments: (a) θ = π/2 ; (b) θ =3π/2 .

3.2 阵列布置鱼形水翼在均匀流中的振荡运动

本节模拟了9个以阵列布置水翼在均匀流中的振荡运动以研究鱼群游动中振荡水翼表面的力模式和尾涡形态.为了简化计算模型, 研究聚焦于中心位置的水翼, 该水翼以(6)式指定的模式进行振荡运动, 其余水翼保持静止.几何模型如图8所示, 计算区域设置为 [–4b, 12b] × [–4b, 4b].水翼之间的横向间距和纵向间距相同, 研究了三个间距下鱼群的力系数特性, 间距 G分别为 G = 0.25,0.5 和 1 m.根据 3.1 节的网格收敛性结果, 采用均匀网格对水翼附近 [–2.5b, 2.5b] × [–1.5b, 1.5b]的区域进行细化, 在x和y方向的最小网格间距分别为 0.005 和 0.004 m, 对应于 3.1 节的收敛网格.边界条件和计算参数的设置类似于3.1节, 本节研究了五个振荡角频率 (ω = 0.15, 0.30, 0.45, 0.60和0.75)下中间位置水翼的力模式和尾涡形态.

图8 阵列布置水翼游动的几何模型Fig.8.Geometrical model of the undulating hydrofoil in array arrangement.

首先对中间位置水翼表面的升阻力系数进行了定量分析.图9比较了三个位置下群游中间位置水翼和单水翼游动中振荡频率和水翼表面力系数之间的关系.在水翼间距 G = 1 m 时, 中间位置水翼的平均阻力系数和升力均方根系数与3.1节单水翼绕流算例近似, 由于距离远, 邻近水翼对中间位置水翼的影响较为微弱, 中间水翼保持相对独立的游动模式.当间距缩小为 G = 0.5 m 时, 水翼之间存在显著的水动力相互作用, 力系数对振荡频率的敏感性增加, 当振荡频率从ω = 0.15增加到ω =0.75 时, 平均阻力系数从 0.32 降低到–0.38, 升力均方根系数从0.41增加到4.11.当缩小到最窄间距G = 0.25 m时, 此时水动力相互作用最为显著,在振荡频率从 ω =0.15 增加到 ω =0.75 时, 平均阻力系数从0.32降低到–0.67, 升力均方根系数从0.75增加到10.10.值得注意的是, 相比单水翼绕流, 水翼表面阻力变为推力的临界频率从大约0.60减小到0.45左右, 说明水翼群游时, 水翼可以以更少的能量消耗获得同样推力.在低振荡频率ω =0.15时, 三个间距下的水动力相互作用力系数变化不大, 因为此频率下脱落的尾涡结构相对稳定独立, 受到周围水翼的影响小.在振荡频率 ω =0.75 时, G = 0.25 m 时群游中间水翼相比单水翼绕流可获得2.5倍的平均推力增加, 升力均方根系数可获得3倍增加, 说明增大激励频率可放大水翼之间的水动力相互作用, 也说明群游水翼在相同的运动形式下能获得更高的游动效率, 这反映出鱼类群游的水动力学优势.

图9 三个间距下群游中间位置水翼和单水翼在五个频率下的力系数比较 (a)平均阻力系数; (b)均方根升力系数Fig.9.Comparison of force coefficients between the central hydrofoil in array arrangement for three gap distances and the single hydrofoil under five frequencies: (a) Average drag coefficient; (b) root mean square lift coefficient.

图10展示了G = 0.5 m时群游中间位置水翼游动中三个频率下升力和阻力系数随无因次时间的变化趋势.群游水翼的升阻力系数变化趋势与单水翼类似, 然而相比单水翼游动, 随着激励频率的增加, 升力系数以更大的幅值和更大的偏斜呈现周期性的振荡特性, 阻力系数呈现更加明显的波谷平坦和波峰尖锐的现象.

图10 群游中间位置水翼在三个频率下升阻力系数随无因次时间的变化趋势 (a)阻力系数; (b)升力系数Fig.10.Variation of drag and lift coefficients on the central hydrofoil in array arrangement with the dimensionless time under three oscillation frequencies: (a) Drag coefficient; (b) lift coefficient.

图11—图13展示了群游水翼在间距G =0.5 m时三个频率下两个时刻的局部涡场.中间位置水翼的游动特性和力系数受两个因素影响.首先, 中间水翼的头部位于上游水翼尾部的低压区,其尾部位于下游水翼头部的高压区, 此压差增加了中间水翼的推力; 其次, 中间水翼与其并列的上下两水翼形成一个半封闭区域, 两水翼之间的流体受到限制而沿着水翼游动方向流动, 减少了能量向两侧的耗散, 提高了推进效率.当中间水翼尾部向上摆动到最大幅值时(如图11(a)), 中间水翼上表面的流体由于伯努利效应而得到加速, 导致压力降低, 下表面的流体减速, 导致压力升高, 因此, 上下表面之间压差形成向上的升力, 相比单水翼游动时增大了升力幅值, 如图4所示; 当中间水翼尾部向下摆动到最大幅值时(如图11(b)), 中间鱼体形成向下的升力.随着Sr数增加, 鱼体之间的剧烈水动力相互作用, 扭曲了周期性的对称交替涡泄, 导致升力系数的正弦振荡特性发生显著的偏离, 如图10中的升力系数曲线.而且, 由于激励频率增加, 相速度增加, 涡列以更快的速度脱落, 在下游鱼体的阻碍作用下, 涡列沿着下游鱼体上下表面发生分叉, 形成两列紧密布置的涡列向远处, 如图12和图13所示.

图11 ω = 0.15 时在两时刻的水翼群游局部涡场图(a) θ = π/2 ; (b) θ =3π/2Fig.11.Local instantaneous vortex contours of undulating hydrofoils in array arrangement for ω = 0.15 at two moments: (a) θ = π/2 ; (b) θ =3π/2 .

图12 ω = 0.45 时在两时刻的水翼群游局部涡场图 (a) θ=π/2; (b) θ =3π/2Fig.12.Local instantaneous vortex contours of undulating hydrofoils in array arrangement for ω = 0.45 at two moments: (a) θ = π/2 ; (b) θ =3π/2 .

图13 ω = 0.75 时在两时刻的水翼群游局部涡场图 (a) θ=π/2; (b) θ =3π/2Fig.13.Local instantaneous vortex contours of undulating hydrofoils in array arrangement for ω = 0.75 at two moments: (a) θ = π/2 ; (b) θ =3π/2 .

4 结 论

本文采用自主开发的径向基函数虚拟网格法研究了绕单个或阵列布置柔性水翼的不可压缩流动问题.采用虚拟网格法在浸入边界内部布置有限数量的虚拟网格以施加物面边界条件, 引入紧支径向基函数通过物面Lagrangian质点拟合物面等值面函数以追踪任意复杂的物面边界.在绕单柔性水翼流动计算中, 网格收敛性研究表明当网格细化时升阻力系数迅速趋近于参考结果, 同时, 不同激励频率下的升阻力系数与参考结果符合较好, 验证了本文方法的精度和可靠性.而且, 平均阻力系数随振荡频率的增加迅速降低, 临界频率为 ω =0.6 , 此时水翼在流体中运动所遭受的阻力变为推力, 即水翼高频横向振荡产生了自身向前运动的推力.同时, 在流场中观察到尾涡交替脱落的卡门涡街现象, 随着振荡频率增加, 涡列以更快的速度脱落,涡之间甚至发生相互排挤和变形.然后, 分析了不同间距和不同激励频率下阵列布置水翼的水动力特性, 随着间距缩小, 力系数呈现显著的放大效果,在间距 G = 0.25 m 时, 推力系数的临界频率减小为 ω =0.45 , 说明紧密布置的水翼在相同的运动形式下能获得更高的游动效率.