基于CFD的仿生水动力学数值计算方法及其验证

冯亿坤，苏玉民，徐小军，刘焕兴，王兆立

（1.国防科技大学智能科学学院，长沙 410073；2.哈尔滨工程大学水下机器人技术重点实验室，哈尔滨 150001；3.北京特种机械研究所，北京 100143；4.北京控制工程研究所，北京 100094）

0 引 言

海洋覆盖了地球表面的71%，蕴藏着丰富的自然资源，但海洋的平均深度达到三千多米，最深处超过一万米，这是人类所无法直接到达的。因此，各种各样的水下运载器应运而生，并且得到了快速的发展。桨-舵推进系统作为最成熟的水下推进与操纵装置，被广泛地应用于水下运载器的研制，然而其存在体积大、噪音高、推进效率低以及操纵性能差的缺点。受水中鱼类游动的启发，研制具备鱼类优异游动性能的仿生水下机器鱼在近三十年得到了广泛的研究，从麻省理工学院在1994年研制成功世界上第一条真正意义上的仿生机器鱼Robo-Tuna[1-2]开始，国内外科研工作者研制了各类仿生机器鱼，进行仿生机器鱼研究的首要任务是对鱼类游动机理的研究。由于鱼类游动性能的复杂性，采用活鱼进行试验无法精确地得到鱼体的水动力性能和流场信息，因而难以对鱼类的游动机理进行详细的参数化研究。随着计算流体流体力学（CFD）方法的发展，研究者对鱼类的游动开展了大量的数值模拟研究。Yang 等[3]对仿生金枪鱼尾鳍的水动力性能进行了数值计算，对尾鳍的刚性和柔性进行了对比；张曦等[4]数值计算了仿生金枪鱼尾鳍、仿海豚尾鳍和仿生鲸鱼尾鳍的水动力性能，并进行了对比分析；Li和Su[5]对仿生胸鳍的推进性能进行了数值计算，并对胸鳍运动过程中的三维涡结构与胸鳍表面压力分布进行了详细分析；张纪华等[6]采用实验的方法，对柔性胸鳍的水动力性能进行了水池试验；Zhu 等[7]采用边界元法（BEM）对仿生鱼在均匀来流中的水动力性能进行了数值计算；Yang 和Su[8]对仿生金枪鱼在均匀来流中水动力性能进行了数值计算，并分析了流场涡结构。考虑到自然界中的鱼类是在自身推进下的自主游动，Borazjani 和Sotiropoulos[9]对仿生鱼的单自由度（1-DoF）自主游动进行了数值计算。由于鱼类在进行周期性游动时，身体躯干及尾鳍的运动是周期对称的，因此可以看作是尾鳍运动平面内的自主游动；Xia等[10]、Li等[11]和Feng等[12]对仿生金枪鱼在水平面内的三自由度（3-DoF）自主游动进行了数值计算，并且分析了仿生鱼产生的三维尾涡结构以及尾鳍表面的前缘涡（LEV）在推力产生中的作用；刘焕兴等[13]对仿生机器鱼的摆动-滑行游动模式进行了数值计算，结果表明滑行比大于等于0.6时，摆动-滑行游动模式比连续摆尾游动模式消耗更少的能量。

考虑到鱼类游动的复杂性，通常为多个方向的平移/旋转运动耦合而成的刚性运动或者三维空间中的柔性运动，因而难以针对CFD 方法在求解鱼类游动性能的有效性开展实验验证。本文首先建立适用于求解仿生水动力学问题的CFD 网格划分策略，其次基于CFD 方法建立仿生鱼水动力学中的关键运动部件和游动问题的数值计算方法，最后针对仿生尾鳍、仿生胸鳍和仿生鱼在均匀来流中的水动力性能以及仿生鱼的3-DoF 自主游动性能展开数值计算，并与文献中的实验与数值计算结果进行对比，验证计算方法的有效性。本文系统地给出了基于CFD 二次开发的仿生水动力学数值计算方法以及数值算例和验证结果，对仿生水动力学在工程中的应用具有重要的意义。

1 CFD数值计算方法

1.1 控制方程

自然界中鱼类游动的介质为不可压缩粘性流体，其运动满足时均连续性方程和雷诺平均纳维-斯克托斯(RANS)方程：

式中，ui和uj为流体速度分量的时均值（i，j=1，2，3），p为流体压力时均值，μ为流体动黏性系数，ui'uj'为速度脉动分量乘积的时均值，Si为广义源项。

采用通用CFD 软件包FLUENT 对流场进行求解，基于有限体积法对式（1）和式（2）进行离散，湍流模型采用SSTk-ω两方程模型，选择基于压力的求解器，采用标准壁面函数法对近壁面流动进行处理，采用SIMPLE 算法对压力和速度进行耦合，离散方式中的梯度项选择Least Squares Cell Based，压力项选择Second Order，动量项选择Second Order Upwind，时间项采用1st-Order Implict，计算物体表面采用无滑移壁面条件，计算的边界条件针对不同的游动问题进行不同的设置，这在下文将会详细地给出。

1.2 动网格方法

考虑到鱼类的游动涉及复杂的刚性或柔性变形运动，在数值计算过程中仿生鱼鳍或者仿生鱼表面的网格节点随着计算时间的递进不断地发生着变化，这将改变计算域的几何拓扑结构，因此需要改变计算域中的网格拓扑结构以满足1.1 节的CFD 计算要求。在本文的研究中采用基于弹性光顺的网格变形以及基于尺度和歪斜率标准的局部网格重构的动网格方法。

在基于弹性光顺的网格变形方法中，将任意两个相邻的网格节点之间的连接理想化为相互连接的弹簧。在给定的边节点(计算物体网格节点)上的位移将产生一个力，这个力与沿连接到该节点的所有弹簧的位移成比例。利用胡克定律，作用于网格节点i上的力Fi可以写为

式中，ni为连接于网格节点i的相邻节点数量，kij为网格节点i与相邻节点j之间的弹簧常数，Δxi和Δxj为网格节点i与相邻节点j的位移。kij表示为

式中，kfa为弹簧常数系数。

在平衡时，连接到节点i上的所有弹簧的合力必须为零，产生一个迭代方程为

式中，q为迭代次数。

式（5）的收敛条件由收敛容差ε和迭代次数q控制，当满足ε或m时，式（5）在每个时间步长的迭代解停止。收敛容差ε为

由于移动边界（计算物体表面）的位移已知，因此在所有内部节点上使用雅可比扫描来求解式（5）。当收敛条件满足时，可得到内部所有节点的位置：

式中，n+1和n分别表示下一个时间步和当前时间步的位置。

当计算物体表面的网格运动幅度较大时，例如胸鳍的大角度旋转拍动运动以及身体躯干和尾鳍的大幅度柔性变形运动，基于弹性光顺的网格变形方法不再足以保证网格质量，甚至产生负体积。为解决这一问题，基于歪斜率和尺度标准对质量差的网格单元进行聚类和局部重构。尺寸标准是根据最小长度尺度和最大长度尺度指定的。长度尺度小于最小长度尺度或者大于最大长度尺度以及歪斜率大于设定的歪斜率的网格单元格被标记进行重构。整个过程如下：

（1）由歪斜率和尺度标准对网格单元进行标记；

（2）删除被标记的网格单元；

（3）在被删除的网格单元的空间区域生成扭曲度优于先前网格的新网格单元；

（4）在重构后的网格单元进行物理量的差值计算，保证与之前的物理量相同。

2 仿鱼类游动数值计算方法

2.1 数值计算网格划分策略

正如1.2 节所述，针对仿鱼类运动体的水动力数值模拟之前，需要划分适用于动网格方法的数值计算网格，因此应针对鱼类游动的特点，建立合理的数值计算网格划分策略，使得生成的网格能够满足动网格方法的要求，在保证计算精度的同时减少网格数量。以三维矩形水翼为对象，给出适用于仿生鱼水动力学数值计算的网格划分策略，如图1 所示。整个计算域为一个长方体，水翼上下、左右对称布置于计算域内，与计算域下游边界的距离大于到上游边界的距离，比值通常为3～5。将整个计算域划分为内外两个域，其中内部计算域包含着整个水翼，通常也为长方体，外部计算域为整个计算域减去内部计算域剩余的几何空间。内部计算域的尺寸应保证水翼在运动过程中不会与其边界发生接触。在确定了计算域的几何拓扑结构之后，接下来进行网格的划分。

图1 水翼的数值计算网格Fig.1 Numerical mesh of the hydrofoil

首先对内部计算域进行网格的划分，在水翼表面划分三角形网格，通过控制网格的稀疏程度对曲率小的区域的几何特征进行有效的捕捉，在内部计算域的六个表面划分均匀的四边形网格，以水翼和内部计算域边界的面网格为基础，采用Delaunay三角剖分法对内部计算域进行四面体非结构网格划分，其中内部计算域边界处的网格为金字塔网格。

之后，采用分块结构网格划分的方法对外部计算域进行网格划分。首先，使用等比数列函数将内部计算域六个面上均匀的四边形网格拉伸到外部域的边界生成体网格，控制合理的增长率以保证网格对流场细节的捕捉。在宽度和高度方向上网格迅速拉伸；在纵向上，计算物体上游区域网格的拉伸速度较快，而在下游的尾流区域，网格的增长率保持在2%以下，以保持较高的尾涡分辨率。最后，对于外部计算域的其他区域，以周围的表面网格为条件直接切分生成六面体网格。整个计算域的网格划分情况如图1所示，在本文后续的数值计算中均采用这种网格划分策略。

2.2 表面网格更新方法

自然界中的鱼类在游动过程中，身体上的各类鳍及身体躯干随着时间的推进不断地发生着变化，因此在数值计算过程中它们的几何形状或者空间位置同样需要发生不断的变化，具体表现为表面网格节点坐标的变化。FLUENT 的动网格功能提供了强大的二次开发能力，研究者可以基于FLUENT的用户自定义函数（UDF）功能，采用C 语言编写二次开发程序，对数值计算模型的刚性或者柔性运动进行控制。

针对动网格的UDF 程序的开发，FLUENT 提供了多种宏命令，本文的数值计算中，借助DE⁃FINE_CG_MOTION 宏命令实现各种复杂的刚性体运动，例如拍动翼、尾鳍以及胸鳍的刚性运动；借助DEFINE_GRID_MOTION 宏命令实现各种柔性变形运动，例如身体躯干以及各类鳍的柔性运动。采用DEFINE_CG_MOTION 对采用刚性运动的胸鳍或尾鳍等鱼鳍进行控制时，可以通过控制鱼鳍的三个线速度和角速度实现对鱼鳍表面网格节点的三维坐标进行整体更新，具体过程如图2 所示。采用DEFINE_GRID_MOTION 对仿生鱼的柔性运动进行控制时，首先对面网格单元的所有节点进行遍历，在遍历过程中，首先读取当前时间步下仿生鱼表面网格节点的三维坐标，再通过给定的运动模型计算得到下一时刻仿生鱼表面网格节点的坐标，在时间步更新之后对仿生鱼表面网格节点的三维坐标进行更新，具体过程如图3所示。

图2 刚性运动求解过程Fig.2 Solving process of rigid motion

最后，由1.2 节的动网格方法对计算域网格进行变形与重构，由1.1 节的数值方法对流场进行求解，得到鱼鳍或身体躯干受到的水动力。涉及自主游动时，在得到水动力之后需要运行鱼体-流体耦合求解程序，对仿生鱼的游动速度进行求解。

2.3 鱼体-流体耦合求解方法

鱼类依靠身体上的各类鱼鳍能够进行各种各样的游动，然而依靠身体和尾鳍的游动通常可以看作是平面内的游动，这时尾鳍和身体躯干的运动通常为平面内的运动。在平面内的运动包含两个平移运动和一个旋转运动，考虑到鱼类的向前游动性能是水下航行器所追求的，因此相关学者对平面的游动进行进一步的简化，只保留前进方向上的运动并进行研究。针对鱼类的自主游动数值模拟，通常不考虑鱼体柔性变形过程中的内力，仅考虑鱼体在外界水动力作用下的运动[9-13]。对于平面内的运动，由质心运动定理：

式中，m为仿生鱼的质量，UX和UY分别为XOY平面下的两个线速度，aX和aY分别为XOY平面下的两个线加速度，FX和FY分别为XOY平面下的两个水动力，ωz、αz、Iz和Mz分别为质心坐标系下的角速度、角加速度、转动惯量和水动力矩。FX、FY和Mz分别通过对仿生鱼表面的粘性力和压力进行积分得到：

式中，p为局部面元dS的压力矢量，nj为j方向单位矢量，τij为黏性应力张量，r为相对于质心的位置矢量。

将式（9）计算得到的水动力和力矩代入式（8），即可对仿生鱼的运动进行求解。针对这一问题的求解通常分为两种方法：一种是求解式（8）得到仿生鱼的线加速度和角加速度，再对加速度进行积分得到游动速度；另一种是由式（8）进行直接求解得到仿生鱼的游动速度。这两种方法在数值求解精度上均能够满足数值计算的需要，在下文的数值计算中将对这两种方法计算得到的游动速度进行对比。两种方法的计算流程如下：

方法1：

（1）初始时刻，仿生鱼的速度和加速度均为零；

（2）第n个时间步，将水动力和力矩代入式（8），计算得到仿生鱼的线加速度和角加速度；

（3）采用二阶龙格库塔积分获得第n个时间步迭代结束时仿生鱼的线速度和角速度，

（4）时间步更新，进入下一次求解。

方法2：

（1）初始时刻，仿生鱼的速度和加速度均为零；

（2）第n个时间步，将水动力和力矩代入式（8），采用三阶向后差分方法，直接计算获得第n个时间步迭代结束之后仿生鱼的线速度和角速度，

（3）时间步更新，进入下一次求解。

在得到仿生鱼的速度信息之后，通过无滑移壁面边界条件将仿生鱼表面的速度施加于流场。无滑移壁面边界条件为

式中，U为仿生鱼表面流体的速度，Us为仿生鱼表面的速度。

式中：UT=(UX,UY)为仿生鱼的平移速度；UR=ωz×r为仿生鱼的旋转角速度；UDEF为仿生鱼表面给定的柔性变形速度，由给定的动力学模型得到。随后进行下一时间步的数值计算，通过时间步的不断递进，实现仿生鱼在自身推进下的自主游动。

3 数值计算案例

3.1 单个鳍在均匀来流中的水动力计算

鱼类身体上有着多个种类的鳍，每种鳍都有着特定的功能，例如尾鳍主要进行快速推进，胸鳍主要进行灵活的操纵。因此，为了便于计算以及机理分析，学者们以单个鳍为研究对象，对水动力性能及机理进行数值研究，以求了解鱼类的游动机理。尾鳍的高效推进性能和胸鳍的灵活操纵性能是研究者所密切关心的，因而得到较为广泛的研究，并形成了若干经典的计算模型。

（1）仿尾鳍运动模型

首先以尾鳍为例，学者们通常以作升沉-俯仰运动的拍动翼模拟尾鳍的运动，并且已经成为尾鳍的水动力数值研究的经典验证算例。Dong[14]以三维矩形翼为研究对象，针对水翼在不同运动参数下的推进性能开展了详细的水池实验，水翼在均匀来流Uin中升沉-俯仰运动，水翼的展长Sfin和弦长Cfin分别为0.6 m 和0.1 m，水翼数值计算模型如图4 所示。关于水翼的详细实验研究见文献[14]，本文不做过多描述，只对水翼的运动模型和相关运动参数进行介绍。水翼的升沉运动hfin(t)和俯仰运动θfin(t)为

图4 水翼数值计算模型Fig.4 Computational model of the hydrofoil

式中，h0和θ0分别为水翼的升沉运动幅值和俯仰运动幅值，f为运动频率，ϕ为俯仰运动与升沉运动的相位差。水翼在运动时的瞬时有效攻角α(t)定义为

水翼在作升沉-俯仰耦合运动时，影响水翼推进性能的一个关键无因次参数为斯特劳哈尔数St，定义为

水翼在给定的运动参数下进行运动时的推力系数CT定义为

式中，Tfin为水翼一个运动周内推力的均值。

在使用CFD 方法对流体中的物体进行水动力数值计算时，用于数值计算的网格的合适与否直接影响着数值计算结果的精度。采用前文给出的数值计算方法和网格划分策略，其中数值计算网格与图1 给出的一致，计算域的出口设置为压力出口，计算域的其他面均为速度入口。由2.1 节给出的网格划分策略，对内部计算域和外部计算域划分三种不同数量的网格，按照网格由密集到稀疏的顺序分别命名为mesh 1、mesh 2 和mesh 3。采用这三种网格对水翼在运动参数为h0/Cfin=0.8、St=0.5、ϕ=90°、Uin=0.4 m/s和θ0=20°下的推进性能进行数值模拟，由计算得到的水翼的推力系数CT对网格数量的无关性进行分析，如表1所示。

表1 不同网格计算值Tab.1 Computation by different meshes

基于Stern 等[15]提出基于Richard 外推（Richard extrapolation,RE）法给出网格数量的收敛性验证的方法和步骤。在进行网格无关性验证时，首先需要给出网格细化率rg。一般来说，rg的定义是粗细网格的单元格大小之比。Baek 等[16]给出了采用非结构网格进行CFD 数值模拟时网格细化率的计算方法，这个方法已经被广泛地应用于船舶螺旋桨、水下航行器以及水面船等多种涉及CFD网格数量的无关性验证中[17]。rg的计算公式为

式中，Nfine和Ncoarse分别表示相邻尺度的网格系统中细密和粗糙的网格的总数目，d表示求解问题的空间维度。rg在本节的验证中约为1.2。

以si(i=1,2,3)分别表示mesh 1、mesh 2 和mesh 3 计算得到的结果，则不同网格系统计算结果的差值为

网格系统的收敛率Rg定义为

Stern等[15]依据Rg的值将网格系统的收敛情况分为以下四种：

（a）0＜Rg＜1时，单调收敛；

（b）Rg＜0且 ||

Rg＜1时，振荡收敛；

（c）Rg＞1时，单调发散；

（d）Rg＜0且 ||

Rg＞1时，振荡发散。

三组网格的Rg约为0.104，表明本文所采用的网格划分策略得到的网格在求解仿鱼类运动体的水动力性能时具有较好的收敛性。在本文后续的研究中所采用的数值计算网格均采用这种方法进行网格的独立验证。

水翼的运动参数为h0/Cfin=1.0，St以0.05为间隔从0.25变化至0.4，ϕ=90°，Uin=0.4 m/s，α(t)的最大值等于30°，获得每个St下水翼的CT。图5显示了本文计算结果与Dong[14]的实验结果和数值计算结果的对比，从图5 中可以看出本文的计算结果与实验结果吻合较好，并且优于Dong[14]采用的前缘涡分离（LVS）的BEM 所计算得到的结果。由于BEM 不考虑流体的粘性，无法对水翼受到流体作用的粘性阻力进行精确的计算，因而得到的推力系数较实验结果和本文的计算结果偏高。从图中可以看出，本文基于CFD技术建立的仿生水动力学求解方法能够对仿生拍动翼的水动力性能进行精确的求解。本文采用Hunt 等[18]提出的Q准则对流场中的三维涡结构进行可视化处理，通过Q等值面对流场中的三维涡结构进行显示，Q的定义为Q=(‖W‖2-‖S‖2)/2，式中W和S分别表示速度梯度的反对称张量和对称张量，‖ ‖为矩阵的欧几里德范数。根据Hunt等[18]的定义，当Q大于零时，Q等值面为旋涡结构，本文中对流场的三维涡结构的显示均采用Q准则。图6 给出了水翼在运动过程中某一瞬间流场的三维涡结构。从图中可以看出水翼上表面靠近前缘的区域产生了一个强度较大的前缘涡，水翼上表面被前缘涡附着的区域被强度较大的低压所覆盖，如图7 所示；在水翼的两个梢部有明显的梢涡产生，两个梢涡向上卷曲，与前缘涡之间有着明显的分界；在水翼的尾缘处存在明显的尾缘涡，尾缘涡与水翼的尾缘相连并脱落于尾流场中。由三维涡结构特征可以看出，本文采用非结构网格与结构网格混合的分块网格划分方法，对包含计算物体的流场区域的加密能够有效地捕捉水翼表面产生并脱落于流场中的涡结构特征。

图6 水翼产生的三维涡结构Fig.6 Three-dimensional vortex structure generated by hydrofoil

图7 水翼表面的压力分布云图Fig.7 Pressure distribution on hydrofoil surface

（2）仿胸鳍运动模型

接下来以胸鳍为例，鱼类的胸鳍有着多个旋转运动组合而成的复杂运动，与尾鳍的运动模式有着明显的区别，基于前人对鱼类胸鳍的实验研究中给出的几何模型作为数值计算模型，对本文基于CFD技术建立的仿生水动力学求解方法在求解胸鳍水动力性能上的有效性进行验证。Kato和Furushima[19]针对黑鲈鱼的实验研究建立了真实黑鲈鱼胸鳍的几何模型，如图8所示。与Suzuki等[20]对黑鲈鱼胸鳍的水动力实验结果进行对比，在Suzuki 等[20]的实验中，胸鳍的展长Spec为0.1 m，弦长Cpec为0.08 m，厚度为0.7 mm，展弦比（AR=，Apec=0.005 834 m2为胸鳍的投影面积）约为1.7。靠近鱼体的一侧为鳍根，远离鱼体的一侧为鳍梢，迎流的一侧为前缘，去流的一侧为尾缘，鳍根与转动中心的距离为0.055 m。

图8 胸鳍数值计算模型Fig.8 Computational model of the pectoral

胸鳍的运动由前后的拍动运动和纵倾运动组成，胸鳍的拍动运动和纵倾运动如图9 所示。胸鳍的拍动运动和纵倾运动方程[20]为

图9 胸鳍的旋转运动Fig.9 Rotational motion of the pectoral fin

式中，θRC为拍动角的周期平均值，θRA为拍动角幅值，ωpec=2πf为胸鳍运动的圆频率，θFC为纵倾角的周期平均值，θFA为纵倾角幅值，ΔθF为纵倾运动与拍动运动的相位差。

胸鳍在均匀来流Uin=0.201 m/s 中进行运动，运动参数设置为：θRC=30°，θRA=30°，θFC=30°，θFA=30°，ΔθF=90°，f=2.0 Hz。通过数值计算获得胸鳍的纵向力系数CXpec和侧向力系数CYpec，与Suzuki 等[20]的实验结果进行对比，验证本文的数值计算方法在求解胸鳍水动力性能上的有效性，纵向力系数CXpec和侧向力系数CYpec定义为

式中，FXpec为胸鳍的纵向力，FYpec为胸鳍的侧向力。

采用前文给出的网格划分策略生成胸鳍的数值计算网格，如图10 所示，整个计算域的网格数量为3 335 564，采用与水翼相同的数值计算方法。

图11给出了本文的数值计算方法得到的胸鳍纵向力系数CXpec和侧向力系数CYpec与Suzuki 等[20]的实验结果的对比。从图中可以看出，实验值的变化趋势发生改变时会出现小幅的波动，并且计算结果与实验值的幅值存在一定的差异，但是两者的总体变化趋势是一致的，并且吻合得较好。由以上分析可以得出，本文基于CFD 技术建立的仿生水动力学求解方法能够对具有多旋转运动的胸鳍的水动力性能进行有效的求解。

图11 胸鳍的水动力系数对比Fig.11 Comparison of hydrodynamic coefficients of pectoral fin

图12 和图13 分别给出了胸鳍运动过程t/T=0.875 时流场的三维涡结构和表面压力分布云图。从图12 可以看出在胸鳍的前缘附着一个强度较大的前缘涡，这导致了一个强度较大的低压区域（图13）；在胸鳍的梢部、根部以及尾缘处均产生了强度较大的涡，并且这些涡与胸鳍连接在一起并脱落于尾流场中，在尾流场中存在着先前运动周期所产生的涡环以及梢涡。尽管胸鳍的运动为空间内复杂的旋转运动，会给网格带来高度的扭曲变形，但是本文的网格划分策略同样能够对流场中涡结构进行较好的捕捉。胸鳍上表面的前缘区域附着一个强度较大的前缘涡，前缘涡附着的区域被强度较大的低压所覆盖。

图12 胸鳍产生的三维涡结构Fig.12 Three-dimensional vortex structure generated by pectoral fin

图13 胸鳍表面的压力分布云图Fig.13 Pressure distribution on pectoral fin surface

3.2 仿生鱼在均匀来流中的水动力计算

鱼类在游动过程中身体上的各类鳍必定与躯体产生相互干扰，采用单个鳍进行鱼类游动机理的研究无法计入这种干扰的影响。因此以整鱼为数值计算模型，研究各类鳍在躯体干扰下的水动力性能对进一步掌握鱼类游动机理是十分有必要的。整鱼的运动涉及躯体的柔性波动运动，因此难以开展精确的实验研究，目前针对这类问题的数值计算方法验证通常以3.1 节中的水翼或者胸鳍作为替代，这不能充分地证明计算方法是合理有效的。为有效开展仿生鱼在均匀来流中的水动力性能数值验证，本节以Zhu 等[7]给出的RoboTuna 仿生机器鱼几何模型为验证算例。国内外多个科研团队以Ro⁃boTuna为数值计算模型，采用BEM 或者CFD 方法对其游动过程中的水动力性能进行了数值计算及试验验证，通过与这些数值计算结果的对比，能够有效地验证本文所采用的数值计算方法在求解整鱼水动力性能上的有效性。

仿生鱼数值计算模型如图14 所示，图中红色虚线表示鱼体的中线，O-XYZ为全局坐标系，o-xyz为固定在鱼体上的随体坐标系，x轴与鱼体的中线重合由头部指向尾鳍，坐标原点位于距离鱼体前端30%身体躯干（不包含尾鳍）长度处，y轴指向鱼体右侧，z轴满足右手准则。仿生鱼几何外形上下左右对称，身体躯干的长度L为1 m，外形为纺锤形，垂直于x轴的截面为椭圆形，椭圆的长短轴之比为3:2，长轴与z轴平行。因此，由身体躯干在中纵剖面上的轮廓线可以确定身体躯干的几何外形，轮廓线的曲线方程为

图14 仿生鱼几何模型Fig.14 Geometric model of the bionic fish

尾鳍前缘x(z)LE和尾缘x(z)TE的曲线方程分别为

仿生鱼在游动过程中，身体躯干进行侧向的柔性波动运动，尾鳍跟随身体躯干进行侧向平移运动的同时围绕着尾柄末端进行摆动运动。身体躯干的侧向柔性波动运动方程为

式中：a(x)为柔性波动运动的振幅包络线函数；k为波数，定义为k=2π/λ，λ为波长；a1和a2均为常数，取不同的值可以调整尾鳍的侧向平移幅值。

尾鳍在尾柄末端的带动下作侧向平移运动的同时绕着尾柄末端作摆动运动，尾鳍的摆动运动方程为

式中，xp为尾柄末端坐标，θm为摆角幅值，Δθ为尾鳍的摆动运动与平移运动之间的相位差。

在验证计算中，仿生鱼被置于均匀的来流Uin=0.7 m/s 中，身体躯干的运动参数为a1=0.002 36，a2=-0.114，f=2.07 Hz，λ=1.675 m，尾鳍的运动参数为θm=18.76°，Δθ=85°。整个计算域网格数量为3 693 734，如图15 所示。图16 给出了仿生鱼纵向力FT的时历曲线，并与Zhu 等[7]的BEM 计算结果以及Zhang等[21]和Chang 等[22]的CFD 计算结果进行了对比。从图16 中可以看出，本文的计算结果与文献中的计算结果吻合较好，FT的幅值小于BEM 的计算结果，这是因为BEM 没有考虑流体的粘性而导致的。此外，本文的数值计算方法得到的结果与Chang 等[22]的计算结果有着较高的吻合度，表明本文的数值计算方法可以有效地获得仿生鱼做柔性运动时的水动力。

图15 数值计算网格Fig.15 Numerical mesh

图16 仿生鱼推力曲线Fig.16 Thrust curve of the bionic fish

3.3 仿生鱼在静水中的3-DoF自主游动数值模拟

（1）鳗鱼数值计算验证

Kern 和Koumoutsakos[23]对鳗鲡游动模式的仿生鱼在水平面内的3-DoF 自主游动进行了数值计算，并且给出了仿生鱼数值计算模型的详细几何参数以及运动模型，针对仿生鱼在给定的一组运动参数下的自主游动进行了数值计算，并给出了最终的巡游速度，因此非常适用于仿生鱼自主游动的数值验证。仿生鱼体长La=1 m，仿生鱼的三维几何模型由空间上变化的椭圆横截面构成，椭圆的长短轴由仿生鱼高度和宽度方向上的轮廓线h(s)和w(s)控制，w(s)的定义为

式中，s为仿生鱼表面到头部顶点的纵向距离，wh=sb=0.04La，st=0.95La，wt=0.01La。h(s)的定义为

式中，a=0.51La，b=0.08La。

仿生鱼宽度和高度方向上的轮廓如图17所示。

图17 仿生鱼宽度和高度方向上的轮廓图Fig.17 Outline of the bionic fish in the width and height directions

仿生鱼在自主游动过程中某一瞬间身体的柔性波动如图18所示。在自主游动过程中，仿生鱼身体的侧向柔性波动方程为

图18 仿生鱼身体几何构型Fig.18 Geometry of the anguilliform fish

在验证计算中，流体的动粘滞系数μ=1.4×10-4kg/ms，流体的密度ρ和仿生鱼的密度ρfish相同且均为1 kg/m3，f=1 Hz。网格划分策略及数值计算方法与3.2节的一致，由于涉及仿生鱼的自主游动，将入口边界的流体速度和压力梯度设置为零，出口边界的速度梯度和压力梯度设置为零。仿生鱼由静止开始摆动身体进行游动，图19 给出了鳗鲡游动模式下仿生鱼的前进速度Uf的时历曲线，其中红色实线是三次向后差分方法的计算结果，黑色点虚线是二阶龙格库塔积分的计算结果，从图中可以看出两种数值算法的计算结果是一致的。本文的数值计算方法得到仿生鱼最终的巡游速度Ucf=0.381La/s，与Kern和Koumoutsakos等[23]的数值计算中得到的0.40La/s巡游速度吻合较好（相差-4.75%），表明本文的数值计算方法能够对仿生鱼的自主游动过程进行精确的模拟。图20 给出了鳗鱼在自主游动过程中某一瞬间流场的三维涡结构。从图中可以看出鳗鱼的尾涡为双排交错排列结构，漩涡在身体弯曲的一侧产生，随着身体上行进波向后的传播，逐渐由头部向尾鳍移动，在身体躯干的后半段形成一个明显的“发卡”涡，随着尾鳍运动方向的转变，脱落于为流场中，形成一个一个尾涡结构。

图19 前进速度Fig.19 Forward velocity

图20 鳗鱼产生的三维涡结构Fig.20 Three-dimensional vortex structures generated by the eel

（2）金枪鱼数值计算验证

接下来以哈尔滨工程大学水下机器人技术重点实验室研制的“仿生-I”号机器鱼的定向游动实验，进一步地验证本文的数值计算方法在求解仿生鱼自主游动上的有效性。“仿生-I”号机器鱼及数值计算模型如图21 所示。关于“仿生-I”号机器鱼的详细试验见文献[24，25]，关于仿生鱼的运动模型以及运动参数见文献[12]，文献[26]给出了仿生鱼数值计算模型的详细型值，包括身体躯干、第一背鳍、第二背鳍、胸鳍、臀鳍和尾鳍，可以用于数值计算模型的几何建立，以上内容本文不再给出。通过改变仿生鱼的游动频率f，得到最终的巡游速度Ucf，图22 显示了由本文3-DoF 数值计算方法得到的结果与“仿生-I”号机器鱼的试验结果[22]以及苏玉民等[27]的1-DoF 数值计算结果对比，图中虚心点为计算或试验的结果，点画线为二次拟合的结果，从图中可以看出本文的数值计算结果与试验结果的趋势一致，Ucf随f的增加均表现为非线性增长趋势，并且优于苏玉民等[27]的1-DoF 数值计算结果中Ucf随f的增加表现为线性增长趋势。另外，除了f为0.8 Hz时本文的计算结果与试验值相差较大（20.5%）以外，其他情况下的差值均在10%以下。以上分析进一步验证了本文的数值计算方法在求解鱼类自主游动方面的有效性。

图21“仿生-I”号机器鱼及数值计算模型Fig.21“FangSheng-I”robotic fish and computational model

图22“仿生-I”号机器鱼的巡游实验对比Fig.22 Comparison of the cruising experi⁃ments of“Fang Sheng-I”robotic fish

4 结 论

仿生水动力学长期都是计算流体力学中的研究热点，其关注的是水中生物高效、快速和高机动游动行为的机理，在水中推进和减阻方面有着重要的工程应用前景。本文基于CFD方法，结合动网格方法，给出了仿生水动力学数值计算方法，针对单个鱼鳍的水动力、均匀来流中仿生鱼的水动力以及静水中仿生鱼的自主游动等仿生水动力学问题开展了研究，并与文献中经典的且可对比的实验或数值结果进行了比较，验证了本文的数值计算方法在求解仿生水动力学问题上的通用性和有效性。本文的主要结论如下：

（1）基于弹性光顺和局部网格重构的动网格方法能够处理仿生尾鳍的大幅拍动运动、胸鳍的大角度旋转运动以及身体躯干的柔性波动运动，在处理各类仿生运动带来的网格大变形运动上有着较强的健壮性；

（2）采用多块网格划分策略能够有效地处理仿生水动力学的计算问题，对包含运动体的内部区域采用非结构网格划分方法，不包含运动体的外部区域采用结构网格划分方法，这种网格划分策略能够在保证计算精度的前提下，减少网格数量，提高计算效率；拍动翼、胸鳍以及整鱼的尾流场中的三维涡结构的显示结果表明，这种网格划分方法能够对流场中核心的涡结构特征进行捕捉。

（3）通过与文献中经典的拍动翼和胸鳍的水动力试验结果以及仿生鱼在均匀来流中进行约束游动时的水动力性能的对比，以及与文献中的仿生鱼在均匀来流中的约束游动或者静水中的自主游动数值计算的对比，表明基于CFD方法的二次开发能够对仿生水动力学问题进行有效的求解。

本文的数值计算方法为仿鱼类游动水动力学走上工程应用提供了有力的手段，能够服务于仿生推进器以及鱼类游动的机理研究，本领域的研究者同样能够依据本文给出的各种经典验证算例对不同的仿生水动力学数值计算方法开展验证。