尾鳍摆动推进水动力与流动涡结构特征分析

苏博越，翟树成，白亚强，张 军

（中国船舶科学研究中心，江苏 无锡 214082）

0 引 言

大多数鱼类通过身体和（或）尾鳍的摆动，推动水流，使鱼体向前运动，人们将这种运动仿生定义为BCF 模式（body-caudal fin）。这种模式通常具有游动速度快、推进效率高的特点，从而成为仿生推进研究的热点。1994年，MIT研究制造出世界上第一条机器鱼，仿金枪鱼机器鱼Robotuna[1]，此后各国高校与科研院所纷纷开展仿生机器鱼研制和仿生推进技术研究。英国ESSEX 大学的胡豁生团队[2-6]研制出了G 系列机器鱼，能以90°/0.2 s 的最大速度实现大角度、高频率的运动。2009 年，哈尔滨工程大学[7]以蓝鳍金枪鱼为蓝本研制出仿生水下机器鱼“仿生-Ⅰ号”，依靠月牙形尾鳍摆动获取推力，尾鳍摆动频率为1.33 Hz，游动速度最高可达1.2 m/s。随着计算机科学、计算方法和流体仿真软件的发展，数值模拟和仿真方法已成为研究鱼类复杂游动现象的一个有效手段。日本Chiba大学的刘浩[8-10]在20世纪90 年代中后期引入人工伪压缩因子，采用有限体积法对Navier-Stokes 方程进行离散化，实现计算流体动力学（CFD）的解算、流场显示及动力分析等后处理操作，并将CFD 方法成功应用于幼鱼游动、飞蛾悬停等生物游动、飞行的复杂流动机理研究。Wolfgang[11]等结合活鱼的试验数据，使用三维面元法计算了金枪鱼的游动，并分析了尾流中涡旋之间的相互作用。夏全新（2005）[12]采用有限体积法和动网格技术，对鱼类的波动游动方式进行数值模拟。通过对二维零厚度板的波状摆动推进进行计算，研究了在高雷诺数鱼类匀速运动状态下的流场分布，并在尾迹中发现了“之”字形排列的涡环结构。2012 年张曦[13]对蓝鳍金枪鱼尾鳍摆动推进水动力性能进行了多角度的分析，包括在振动半圆柱扰流场中二维振荡水翼的推进性能、尾鳍形状和尾鳍柔性以及非对称运动对摆尾推进性能的影响等。2015年，白亚强等[14]基于弹性光顺模型和局部网格重构模型动网格技术，建立了在大摆幅运动下柔性长鳍扭波推进水动力及流场涡结构的数值计算方法，并通过水动力试验从流场和水动力方面验证了数值计算方法的正确性。2017 年李宁宇等[15]采用基于RANS 方程的数值模拟方法研究了金枪鱼自主游动的水动力特性，并通过模拟金枪鱼运动姿态以及巡游状态受到的阻力与游动速度的关系，求解了金枪鱼从静止到巡游状态下的自主游动过程。以上研究大多着眼于宏观水动力结果，对鱼游流场的精细化空间涡结构的计算和涡结构与水动力相互关系的研究尚不够深入。2021年张栋[16]对近年生物游动和飞行等过程中的涡结构模拟及力学特征规律进行了综述，并提出了开展真实雷诺数下生物游动飞行精细数值模拟和涡结构对生物体游动与飞行过程中受力贡献方面研究的展望。

本文主要从BCF 推进运动模式出发，研究建立仿金枪鱼尾鳍摆动推进粘性流场精细化数值计算方法以及水动力、涡结构提取与分析方法，并通过网格无关性验证及与相关文献结果的比较等验证本文方法的有效性。采用本方法数值研究不同运动参数组合下尾鳍摆动推进水动力随St数的变化规律，结合尾鳍摆动运动、推力、表面压力等，分析尾鳍摆动流动涡结构特征及其演化过程。

1 计算方法介绍

1.1 BCF运动模式

在BCF 推进模式中，尾鳍运动通常由围绕尾柄末端的摇摆和由尾柄驱动的横向振荡运动组成。计算模型的坐标系如图1所示。来流速度的方向与x轴的正方向相同。尾鳍端部沿z方向横向往复移动，同时围绕通过尾鳍端部平行于y轴的旋转轴线摇摆[17]。

图1 尾鳍运动坐标系Fig.1 Coordinate system of caudal fin motions

尾鳍横向运动的位移z(t)和尾鳍运动坐标系中尾鳍轴线的摇摆角度θ(t)满足正弦变化规律，即

式中：A为横移运动的幅值；θ0为摇摆运动的角幅度；f为尾鳍摆动频率，在尾鳍摆动过程中，横移运动和摇摆运动的频率相等；ϕ0为横移运动和摇摆运动的相位差。

在尾鳍运动过程中，影响摆动尾鳍水动力性能的无量纲参数——斯特罗哈尔数（St）表示为

式中：A0表示尾鳍运动过程中脱落尾涡的宽度，近似表达为横移运动幅值的两倍，即A0=2A；V0为来流速度。

通过数值计算可以得到尾鳍在不同运动状态下的水动力系数。水动力系数表示为

式中：Fx、Fz分别为尾鳍所受到的x方向、z方向的分力，Fx方向与推力Ft方向相反，将推力系数记为Kt，则Kt=-Cx；S为尾鳍弦长与展长的乘积；C0为尾鳍弦长。

1.2 数值计算

本文采用商用软件Fluent V17.0进行数值模拟，应用动网格技术来处理尾鳍运动边界随时间运动变形的问题。本文尾鳍边界的运动由UDF程序控制。数值计算中网格动态变化过程采用弹簧光顺模型（spring-based smoothing）和局部重构模型（local remeshing）。

1.2.1 控制方程

本文采用不可压缩粘性流体的雷诺平均方程（RANS）作为流体控制方程，方程如下：连续性方程为

动量方程为

式中，ui为平均速度分量（u1=u,u2=v,u3=w），ρ为流体质量密度，p为平均压力，μ为流体分子粘性系数为雷诺应力项。

湍流模型采用SSTk-ω两方程模型[18]。

本文对流项和空间项均采用二阶精度格式进行离散。

1.2.2 弹簧光顺模型

弹簧光顺模型基于弹簧的平滑，将网格中任意两个连接的网格节点之间的边看作是连接两点的弹簧，所有网格节点的边形成弹簧网络。边界运动前，各个节点的初始间距为网格的平衡状态，当某一边界节点发生位移时，系统将产生一个与该节点连接的所有弹簧位移成正比的力，会引起整个弹簧网络的变化，通过求解弹簧网络重新达到平衡时节点的新位置即可相应地得到新的网格节点位置。

在尾鳍运动过程中，由于运动范围较大，运动速度快，尾鳍表面的网格很容易产生严重的变形。在数值模拟过程中，通过对弹簧刚度系数的不断调整，可以把尾鳍周围的网格变形传递到离尾鳍较远的网格，从而避免尾鳍周围网格的较大变形，保证周围网格的质量。

1.2.3 局部网格重构模型

弹簧光顺模型通常用于网格区域中的动态网格生成。但如果运动边界的位移远大于网格尺寸，或者边界运动方向造成网格的扭曲，即使采用弹簧光顺模型仍可能导致网格质量下降，甚至出现负体积网格，导致网格失真过大引起计算不收敛。采用局部网格重构技术可以很好地控制网格质量，防止因位移变化过大而产生负体积网格。

1.3 计算参数选取与边界条件设置

本文计算网格如图2 所示，在尾鳍附近建立加密区域（黄色区域），其尺寸为5C0×5C0×7.5C0，外围流场尺寸为25C0×12.5C0×12.5C0。图中蓝色表面为速度入口，红色表面为压力出口，加密区域内尾鳍壁面设置为无滑移边界。

图2 网格划分Fig.2 Mesh distributions

尾鳍摆动推进流场为非定常流动，选取时间步长为0.005 s，每个时间步最大迭代次数为50，总时间步数为4000。

1.4 网格无关性验证

在网格无关性验证的数值计算中，尾鳍型值选自文献[13]，其中尾鳍弦长C0=129 mm，展长D=241 mm，尾鳍横移幅值A=100 mm，摇摆角幅度θ0=25°，横移运动和摇摆运动的相位差ϕ0=-π/2，运动频率f=0.5 Hz，来流速度V0=0.31 m/s。

本文首先设计了网格稀疏程度不同的三套网格G1、G2、G3，其网格数分别为120 万、286 万和541万，其加密区域网格数分别为90万、137万和185万。通过对比不同网格数下的计算结果（见表1、图3）可知，采用网格G1的计算结果与其它两个结果相差较大，而采用网格G3的水动力系数均值计算结果与网格数量相关性较小。其中网格G3 的Cxm计算结果与网格G2 的结果偏差小于4%，而对于侧向力系数均值Czm，由于尾鳍运动的对称性，在数值计算达到稳定后，一个周期内的Czm理论值为0，从计算结果来看，可以认为其已达到无关性要求。图3 为三套网格的水动力系数达到稳定后一个周期内水动力系数随时间的变化曲线。因此，后续将选取网格G3的水动力计算结果进行比较分析。

表1 不同数量网格计算结果Tab.1 Results of different meshes

图3 不同网格数下水动力系数随时间变化Fig.3 Hydrodynamic coefficients varying with time for different meshes

1.5 涡结构提取分析方法

本文采用Q准则判据提取尾鳍摆动推进过程中流场的涡特征结构[19]。图4为采用网格G3计算的流场涡结构Q=0.01 s-1结果。可以看到，尾流中涡环的轮廓清晰光顺，较小尺度涡的特征也有保留。虽然继续增加网格会使涡结构更加精细，但同时对计算资源的要求更大，因此，本文选用G3网格进行尾鳍摆动推进流场数值计算，并对涡结构进行提取分析。

图4 尾鳍摆动推进涡结构（网格G3）Fig.4 Vortex structure of caudal fin(Mesh G3)

1.6 水动力性能计算结果验证

2012 年张曦[13]采用RANS 方法计算尾鳍摆动推进水动力，同时开展了尾鳍摆动推进水动力试验，水动力数值计算结果得到了试验结果的验证。从图5中可以看到，本文Cx数值计算结果与文献[13]计算结果趋势一致，量值基本吻合，峰值位置稍有差异，可见本文方法水动力的数值计算结果是合理有效的。

图5 计算结果与文献结果比较Fig.5 Comparison of calculation results with the results in literature

2 水动力随St数的变化

采用前文的方法对黄鳍金枪鱼尾鳍摆动推进水动力与流场进行了计算分析，黄鳍金枪鱼尾鳍弦长C0为0.4 m，展长D为0.66 m（如图1）。

为了研究St数对水动力系数的影响，本文通过系列变化来流速度、摆动频率、横移幅度，来实现St数的改变，具体计算工况参数见表2。如图6 所示，在均匀来流下单尾鳍摆动推进水动力与流动涡结构由参数St主控，在所研究的参数范围内，随着St增大，水动力系数均值Cxm逐渐减小，也即平均推力系数Kt逐渐增大。

表2 计算工况参数Tab.2 Parameters for calculation

图6 Cxm随St变化曲线Fig.6 Calculated results of Cxm with different St

本文将尾鳍尾柄位置横移为0且下一时刻向z轴正向运动的时刻作为一个周期起始，记为t=0。

图7为计算工况3下（St=0.2）一个周期内尾鳍水动力系数Cx、Cz、Cm瞬时变化曲线。可以发现，在t=(1/8)T、(5/8)T时刻附近，尾鳍的瞬时推力最大；在t=(3/8)T、(7/8)T时刻附近，尾鳍瞬时推力最小。由于尾鳍运动是关于x轴镜像对称，因而推力系数Cx变化周期为尾鳍运动周期的一半；而尾鳍运动关于z轴不存在对称性，所以侧向力系数Cz和力矩系数Cm变化周期与尾鳍运动周期相同。

图7 一个周期内尾鳍瞬时水动力变化曲线（工况3）Fig.7 Hydrodynamic coeffcients of caudal fin within one period

3 尾涡结构特征

图8 为工况3 中尾鳍摆动推进一个周期内的涡结构（Q=0.01 s-1）演化过程。在t=(0～2/8)T尾鳍时刻内，尾柄向纸面外（z轴正向）运动，然后在2/8T时刻反向，在t=(2/8～6/8)T尾柄持续向纸面内运动至位移最大处，在t=(6/8)T后改变方向向初始位置运动。由于尾柄平移和尾鳍转动两种运动间存在相位差，本次计算结果中尾涡的脱落时刻也相较于尾柄运动有些滞后。比较图7 和图8 可以发现，在t=(1/8)T、(5/8)T附近尾鳍水动力系数Cx、Cz达到极值，此刻尾鳍处于梢部脱落连续梢涡的阶段，为涡环形成的中期；在t=(3/8)T、(7/8)T附近尾鳍水动力系数Cx、Cz为0，此时尾鳍表面涡完成脱落，并与梢涡组合形成涡环结构，本文定义此涡环结构为初始涡环。一个周期内尾鳍往复运动产生了两个完整的涡环，每个涡环的产生过程中Cx完成一个周期的变化，而由于前后半个周期涡环z向符号发生变化，Cz仅变化了半个周期。

图8 一个周期内尾鳍涡结构演化（工况3）Fig.8 Vortex structure development of caudal fin within one period

图9～12 为不同摆动频率流场涡结构特征（Q=0.01 s-1）随St数变化结果。如图10 所示，在St=0.2时，从尾鳍脱落的涡环在流场中交错排列，靠近尾鳍的涡环尺寸较小，且相邻涡环相互连接。随着尾流演化，相邻涡环开始发生分离。通过比较可以发现在St数较小时，涡环之间的分离较慢，涡结构呈“之”字型；在St数较大时（St≥0.3），相邻涡环之间的分离较快，尾流涡结构变为双列楔形涡；在St数达到0.5 后，由于涡环间距缩短，双列涡的相邻涡环也逐渐开始相互连接，且在尾鳍梢部后方出现上下两条小尺度的“之”字形带状涡结构，其涡强度也随St数增加而增大。将t=0、(4/8)T时刻，上半个周期产生的初始涡环与x轴正向夹角定义为初始涡环角，初始涡环角随St数变化如图13 所示。结合图6、图13 可以发现，在参数范围内，随着St数增加初始涡环角逐渐增大，尾鳍运动产生的推力也逐渐增大。

图9 摆动频率0.5 Hz涡结构（St=0.2）Fig.9 Vortex structure result of 0.5 Hz oscillating frequency(St=0.2)

图10 摆动频率0.75 Hz涡结构（St=0.3）Fig.10 Vortex structure result of 0.75 Hz oscillating frequency(St=0.3)

图11 摆动频率1Hz涡结构（St=0.4）Fig.11 Vortex structure result of 1 Hz oscillating frequency(St=0.4)

图12 摆动频率1.25 Hz涡结构（St=0.5）Fig.12 Vortex structure result of 1.25Hz oscillating frequency(St=0.5)

图13 初始涡环角随St数变化Fig.13 Initial vortex ring angle changing with St

4 结 论

本文采用动网格技术，建立了均匀来流下尾鳍摆动推进流场与水动力的数值计算方法，并通过网格无关性计算和其它文献结果比较验证了计算方法。采用Q准则判据提取尾鳍摆动推进流场的涡结构，并通过比较水动力与涡结构的周期性变化，说明了尾鳍摆动推进水动力与涡结构的关系。同时，本文还分析了涡结构特征与St数的关系。研究表明：

（1）本文数值计算方法可模拟尾鳍摆动推进流场，获取尾鳍摆动推进水动力和涡结构的演化。尾鳍摆动推进尾流精细涡结构，其整体上呈交错分布的环状结构，并随尾流演化逐渐分离。

（2）在均匀来流下，尾鳍摆动推进水动力由参数St主控，在所研究的参数范围内，随着St增大，平均推力系数Kt逐渐增加。

（3）从一个周期涡结构的演化来看，在梢涡持续脱落时尾鳍水动力系数Cx、Cz达到极值，在尾鳍表面涡完全脱落，尾流中出现完整涡环结构时，尾鳍水动力系数Cx、Cz为0。

（4）随着St数增加，初始涡环结构与x轴正向夹角逐渐增大，在St=0.3左右时，尾鳍尾流中涡环结构逐渐从单列演变为双列。