0 引　言

现阶段陆地资源日益稀少，人类开发海洋资源的步伐也在加快，以应对日益增长的需求^[1]。水下机器人是人类进行海洋探索的重要工具，拥有广泛的应用潜力。传统的螺旋桨驱动水下机器人存在噪声大、稳定性差、推进效率低等缺点^{[2 − 3]}。故有必要开发改进的水下机器人。水下生物经过长期的进化，具有独特的游动能力，所以各种仿生机器鱼的高效推进就成了研究者们关注的焦点。

根据形态特征和游动方式的不同，鱼类在海洋中的运动可分为体/尾鳍（BCF）推进模式和中/对鳍（MPF）推进模式^[4]。海洋中80%左右的鱼类采用BCF推进方式，具有游动速度快、推进效率高的特点^[5]。

随着计算机技术与计算流体力学（CFD）的迭代发展，数值模拟方法已成为流体力学研究领域的核心工具。张天宇等^[6]采用流体力学方法，对三维仿生鲀鱼模型的2种尾鳍运动模式（S型和C型）开展水动力性能研究。Zhong等^[7]以金枪鱼为研究对象，采用有限体积法对均匀流场中2种斯特劳哈尔数（St）数下的4种柔性尾鳍推进性能进行数值模拟，结果表明中等柔性尾鳍的平均推力系数较高。Zhou等^[8]利用数值模拟分析了展向柔性对尾鳍游泳性能影响，发现了展向柔性可显著增强推力并提升推进效率。Su等^[9]针对金枪鱼仿生机器人的流体动力学性能展开研究，系统探讨了柔性尾鳍运动参数对机器人游泳性能的影响。Mannam等^[10]结合实验和数值模拟方法，研究了扑翼在不同参数下的推力生成规律与效率特征，并对刚性和柔性扑翼的水动力性能进行了对比研究。耿文豹等^[11]以金枪鱼、海豚和旗鱼为仿生对象，构建了3种机器鱼模型，借助动网格技术设定边界条件，通过数值方法分析了模型在流场中的速度分布与受力特性。张栋等^[12]对近年来生物游动与飞行过程中的涡结构模拟及力学特性规律进行了分析，提出了真实雷诺数下生物运动精细数值模拟及涡结构对生物体推进机制的研究展望。Zhu等^[13]基于边界元法，分析了金枪鱼在BCF推进模式身体涡和尾鳍涡的耦合作用。结果发现两者的相互干扰可有效提升尾鳍的推进效率。

综上所述，研究人员通过理论分析、实验测试、数值模拟等方法对鱼类游泳进行了广泛的研究，取得了丰富的研究成果。然而，上述研究缺乏对尾鳍柔性特性的全面研究。本研究以鲣鱼为仿生对象，通过模拟鲣鱼的运动模式，基于双向流固耦合（FSI）数值模拟方法，系统研究了仿生鲣鱼尾鳍在不同结构特征、材料性能和运动参数下的水动力性能。研究结果为仿生推进系统的设计提供了新的思路，并为其在推进器巡航运动场景中的应用提供了实践指导。

1 模型和方法 1.1 尾鳍几何模型

本文参考生物鲣鱼尾鳍标本如图1(a)所示，建立仿生鲣鱼尾鳍模型。尾鳍的三维模型如图1(b)所示，弦长C=0.09 m，翼展E=0.15 m，迎水面面积S=0.0086 m²。

1.2 尾鳍运动模型

在BCF推进模式下，尾鳍运动主要由尾柄末端摆动和尾柄驱动的横向振荡组成^[14]。计算模型的坐标系如图2所示。其中来流速度方向与x轴正方向一致。尾鳍端部沿z轴方向做横向往复运动，同时绕通过尾鳍端部且平行于y轴的旋转轴线摆动^[14]。

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

$ z(t)={z}_{0}\sin (2\text{π} ft) ，$

(1)

$ \theta \left(t\right)={\theta }_{0}\sin \left(2\text{π} ft-{\phi }_{0}\right) 。$

(2)

式中：$ f $为尾鳍运动频率；$ {z}_{0} $为横移运动幅值；$ {\theta }_{0} $为摆动运动的角幅度；$ {\phi }_{0} $为横移运动和摆动运动之间的相位差。根据文献[15]，尾鳍的横移与摆动的相位差为90°时，尾鳍摆动对推进作用的促进最大，因此成为尾鳍仿真建模设计中所需考虑的关键问题之一。

为实现仿生尾鳍推进性能的定量评估，基于仿真模型构建的XYZ坐标系，分别定义无量纲推力系数$ {C}_{x}(t) $、无量纲侧向力系数$ {C}_{z}(t) $和力矩系数$ {C}_{m}(t) $如下：

$ {C}_{x}(t)=\frac{2{F}_{x}(t)}{\rho S{V}^{2}}，{C}_{z}(t)=\frac{2{F}_{z}(t)}{\rho S{V}^{2}}，$

(3)

$ {C}_{m}(t)=\frac{2{M}_{y}(t)}{\rho S{V}^{2}} 。$

(4)

式中：$ {F}_{x}(t) $为推进力；$ {F}_{z}(t) $为侧向力，为运动中绕y轴的瞬时力矩；$ \rho $为流体密度；V为来流速度，S为尾鳍的迎水面面积。从能量转换的角度分析；输入功率P_in用于克服侧向力$ {F}_{z}(t) $及力矩$ {M}_{{y}}(t) $；而输出功率P_out则定义为推力$ {F}_{x}(t) $以速度V向前推动尾鳍的能量。

$ P_\text{out }=\frac{1}{T} \int_{0}^{T} F_{\mathrm{e}}(t) V \mathrm{d} t，$

(5)

$ p_\text{in}=\frac{1}{T} \int_{0}^{T}\left[F_{0}(t) \dot{z}(t)+M_{y}(t) \dot{\theta}^{\prime}(t)\right] \mathrm{d} t。$

(6)

式中：z′(t)和θ′(t)分别为尾鳍沿z轴的横向运动速度和绕y轴摆动的速度。

平均推力系数$ {C}_{T} $、平均侧向力系数$ {C}_{E} $、平均力矩系数$ {C}_{M} $和平均输入功率系数$ {C}_{P} $定义如下：

$ {C}_{T}=\frac{1}{T}\int _{0}^{T}{C}_{x}(t){\mathrm{d}}t，$

(7)

$ {C}_{E}=\frac{1}{T} \int_{0}^{T} C_{Z}(t) \mathrm{d} t，$

(8)

$ C_{M}=\frac{1}{T} \int_{0}^{T} C_{m}(t) \mathrm{d} t，$

(9)

$ {C}_{p}=\frac{{P}_{{\mathrm{in}}}}{0.5\rho S{V}^{3}}。$

(10)

尾鳍运动的推进效率$ \eta $可以表示为平均推力系数与$ {C}_{T} $平均输入功率系数$ {C}_{P} $的比值。

$ \eta =\frac{{P}_{\mathrm{{out}}}}{{P}_{\mathrm{{in}}}}=\frac{{C}_{T}}{{C}_{P}} 。$

(11)

本文中尾鳍的斯特劳哈尔数(St)定义如下：

$ St=\frac{2{z}_{0}f}{V} 。$

(12)

1.3 数值方法

采用的雷诺平均纳维-斯托克斯（RANS）模型、连续性方程和动量方程。

连续性方程为：

$ \frac{\partial {u}_{i}}{\partial {x}_{i}}=0 ，$

(13)

动量方程为:

$ \rho \frac{\partial {u}_{i}}{\partial t}+\rho {u}_{j}\frac{\partial {u}_{i}}{\partial {x}_{j}}=-\frac{\partial p}{\partial {x}_{i}}+\frac{\partial }{\partial {x}_{j}}\left(\mu \frac{\partial {u}_{i}}{\partial {x}_{j}}-\rho \overline{u_{i}^{{}^{\prime}}u_{j}^{{}^{\prime}}}\right)。$

(14)

式中：$ {u}_{i} $为平均速度分量($ {u}_{1}=u $，$ {u}_{2}=v $，$ {u}_{3}=w $)；$ \rho $为流体质量密度；$ p $为平均压力，$ \mu $为流体分子粘性系数；$ -\rho \overline{{u}^{\prime}{}_{i}{u}^{\prime}{}_{j}} $为雷诺应力项。

湍流模型采用RNG k-ε湍流模型^[16]，该湍流模型在处理分离流动方面表现出色。

k方程：

$ {\begin{split}&\rho \frac{\partial k}{\partial t}+\rho {u}_{j}\frac{\partial k}{\partial {x}_{j}}=\frac{\partial }{\partial {x}_{j}}\left[\left(\mu +\frac{{\mu }_{i}}{{\sigma }_{k}}\right)\frac{\partial k}{\partial {x}_{j}} \right]+\\ &{\mu }_{t}\frac{\partial \overline{{u}_{i}}}{\partial {x}_{j}}\left(\frac{\partial \overline{{u}_{i}}}{\partial {x}_{j}}+\frac{\partial \overline{{u}_{j}}}{\partial {x}_{i}}\right)-\rho \varepsilon，\end{split}} $

(15)

ε方程：

$ {\begin{split}&\rho \frac{\partial \varepsilon }{\partial t}+\rho {u}_{j}\frac{\partial \varepsilon }{\partial {x}_{j}}=\frac{\partial }{\partial {x}_{j}}\left(\mu +\frac{{\mu }_{t}}{{\sigma }_{k}}\right)+\\ &{C}_{1}\frac{\varepsilon }{k}{\mu }_{i}\frac{\partial {u}_{i}}{\partial {x}_{j}}\left(\frac{\partial {u}_{i}}{\partial {x}_{j}}+\frac{\partial {u}_{j}}{\partial {x}_{i}}\right)-{C}_{2}\rho \frac{{\varepsilon }^{2}}{k}。\end{split}} $

(16)

固体区域求解控制方程如下^[17]：

$ {\rho }_{s}\frac{{{\mathrm{d}}}^{2}{d}_{s}}{{\mathrm{d}}{t}^{2}}=\nabla \cdot {\sigma }_{s}+{\rho }_{s}g ，$

(17)

$ \nabla \times {\Gamma }_{s}\times \nabla =0。$

(18)

式中：$ {d}_{s} $为位移矢量；$ {\sigma }_{s} $为柯西应力张量；$ g $为重力加速度矢量；$ {\rho }_{s} $为固体密度。根据广义胡克定律，应变张量$ {\Gamma }_{s} $：

$ {\Gamma }_{s}=\frac{1+{\nu }_{s}}{E}{\sigma }_{s}-\frac{{\nu }_{s}}{E}ITr({\sigma }_{s})。$

(19)

式中：$ E $为杨氏模量；$ {v}_{s} $为泊松比；$ Tr({\sigma }_{s}) $为应力张量的迹。

双向流固耦合相互作用数据流程图如图3所示。在每个时间步内Fluent和Transient Structural分别对流场和结构场进行迭代求解，并通过预先设定的流固耦合交界面传递数据，待流固耦合交界面传递数据收敛后，进行下一个时间步计算直至最终计算完成，该方法在保证求解精度的同时极大提高了求解速度和收敛性。

以上过程需要基于应力连续性和速度连续性的2个数学边界条件。固体力学计算的应力连续边界条件^[13]：

$ -pn+\tau \cdot n={\sigma }_{s}\cdot n。$

(20)

方程的左侧是FSI表面上流体域传递的压力和剪切的有效分量，右侧是FSI面上固态柯西应力张量的有效分量。

流体计算的速度连续边界条件^[18]：

$ u=\frac{{\mathrm{d}}{d}_{s}}{{\mathrm{d}}t} 。$

(21)

1.4 计算域与网格生成

数值模拟中，研究对象表面网格因刚性和柔性变形导致计算域网格拓扑结构变化，需采用动态格子法实时更新网格。利用弹簧光顺法控制变形区域网格整体变形，对于局部网格质量下降或边界扭曲的区域，采用局部网格重构方法。在瞬态模块中，通过添加运动关节模拟鲣鱼尾鳍的摆动和横向运动，依据其实际运动规律设定摆动角度范围及振幅频率。

如图4(a)所示，计算区域为11C×8C×5C矩形，左边界为流体入口（固定流入速度），右边界为压力出口（压力设为0）。原点位于尾鳍前缘中心，前后边界距尾鳍前、后缘比为3∶8，留足尾涡显示空间；上、下边界距前缘中心2.5C，避免地面效应。尾鳍表面为无滑移条件，周边壁面为滑移条件。图4(b)为流体域网格，尾鳍表面细化（稀疏15 mm，精细0.5 mm）；图4(c)显示瞬态结构网格，采用1 mm四面体划分实体，结合独立网格算法（最小0.1 mm）提升计算效率。

2 数值方法的验证

采用数值仿真计算之前，首先需要验证网格无关性，本次计算的模型和参数采用标准算例zhong等^[7]：$ {z}_{0}= 0.1\;\text{m} $，$ \theta =25{^{\circ}} $，$ {\phi }_{0}=90{^{\circ}} $，$ V=0.33\;{\mathrm{m/s}} $，$ f=0.5\;{\mathrm{Hz}} $。

2.1 独立性验证

采用文献中粗、中、细（40万、80万、160万）3种不同计算网格，系统分析了流体域网格数量对计算结果的影响，并对网格无关性进行了验证。图5所示为一个周期内3种网格密度下瞬时推力系数、瞬时侧向力系数和力矩系数随时间变化曲线。结果表明，中等网格、精细网格和粗网格的计算结果与上述文献中的数据呈现出基本相同的变化规律，中等网格、精细网格误差值在7%以内；粗网格的计算结果误差在11%以内。

3 数值模拟结果分析 3.1 摆角幅度对水动力参数的影响

图6为尾鳍摆角幅度对水动力参数的影响规律。变厚度仿生尾鳍材料（杨氏模量为9×10⁷ Pa，泊松比为0.45）设置运动频率为1 Hz、横向位移为90 mm，摆角幅度15°、25°、30°、35°。

由图6(a)和图6(b)可知，瞬时推力系数曲线的波动频率是瞬时侧向力系数曲线的2倍。当尾鳍摆角幅度增加，其瞬时推力系数曲线和瞬时侧向力系数曲线的波动幅值均增加。由图6(c)和图6(d)可知，随着摆角幅度的逐渐增大，平均推力系数$ {C}_{T} $和平均侧向力系数$ {C}_{E} $均呈增大趋势变化。上述结论显示，尾鳍的摆角幅度的增加，可以有效提高推进力，但当摆角幅度达到30°以上时瞬时推力的峰值与谷值差异更显著，曲线震荡更剧烈，导致稳定性下降。

3.2 结构参数对水动力参数的影响

为研究结构特征（外形与材料）对仿生鲣鱼尾鳍推进性能的影响，设计等厚度尾鳍作为外形对照组（与变厚度结构对比），选用仿生材料、聚乙烯（PE）和结构钢（材料属性见表1）。对ES（等厚结构钢）、EP（等厚聚乙烯）、EF（等厚仿生材料）、VS（变厚结构钢）、VP（变厚聚乙烯）、VF（变厚仿生材料）6种尾鳍在不同运动频率下的水动力性能进行仿真分析。设定数值模拟中各参数为：$ {z}_{0}=c $，$ \theta = 25{^{\circ}} $，$ {\phi }_{0}=-90{^{\circ}} $，$ V=0.5\;{\mathrm{m/s}} $，$ f=0.5、0.8、1、1.25\;{\mathrm{Hz}} $。

如图7和图8所示，系统分析了4种不同结构特征的仿生尾鳍在不同运动频率下的水动力性能，主要集中研究了瞬时推力系数、平均推力系数及推进效率的变化趋势。

研究结果表明，变厚度尾鳍（VS/VP/VF）的推进性能在所有频率范围内均优于等厚度尾鳍（ES/EP/EF），验证了外形仿生设计在优化水动力特性方面的重要性，符合鲣鱼尾鳍的真实特性。

值得注意的是，3种材料尾鳍的性能表现出显著的频率依赖性。在低频（如0.5 Hz）时，刚性尾鳍和PE尾鳍凭借其高效的动力传递能力，其瞬时推力、平均推力和效率均高于仿生材料尾鳍；在0.8 Hz和1 Hz的运动频率下，此时仿生尾鳍的st数分别为0.288和0.36，均在真实的鲣鱼的st数范围（0.2～0.4）内，尾鳍的推进性能表现为VF>VP>VS和EF>EP>ES。此外，研究还发现，即使是等厚度设计的柔性在高频环境下也比刚性和PE等厚度尾鳍更具性能优势，进一步表明柔性仿生材料在高频工况下的普适性价值。

通过对各尾鳍推进效率分析发现，VS和VP的高效区间为0.8～1 Hz；0.8 Hz时效率分别为0.4856和0.4945，峰值1 Hz（0.5067、0.5123）较初始效率提升约12.07%和15.44%，高频段（1.25 Hz）因流固振动失谐导致效率锐减25.24%和25.61%；VF得益于动态变形能力在保留低频表现（0.5 Hz效率0.4144）的同时，进一步突破刚性材料限制，实现1 Hz更高峰值效率（0.5282，较VS高4.07%）并向高频延展，尽管1.25 Hz效率仍下降24.29%，但其在1.25 Hz仍比VS高3.97%，证明柔性材料的频带拓展潜力。等厚度尾鳍因几何简化降低了低频涡旋捕获效率，0.5 Hz时推进效率仅为变厚度尾鳍的1/2，却在0.8～1.25 Hz中高频区间展示出超常稳定性，揭示了简约构型在高扰动流场中的鲁棒性特征。

综上所述，虽然低频时刚性尾鳍和PE尾鳍具备一定优势，但随着频率升高（≥0.8 Hz），柔性尾鳍（VF和EF）逐渐展现出更优的推进性能。特别是VF表现出不同运动频率段下的水动力适应性较好，说明仿生设计结构特征可以有效优化仿生鲣鱼尾鳍的推进性能，提高推进效率。

3.3 结构变形和游动参数对流场结构影响分析

1）结构变形对流场结构影响分析

为研究结构变形对流场的影响，在其他参数不变的情况下，仅改变材料参数以分析不同材料尾鳍对流场及推力的作用机制。由图9和图10可知，尾鳍两侧高压与低压之间的压力差呈现出仿生材料>聚乙烯材料>结构钢材料的顺序。结构变形的程度显著影响水动力性能：仿生材料因其卓越的柔性变形性能，增强了流场压力梯度，推力最大；聚乙烯材料次之，其适中变形能力提升了部分推力输出；而结构钢因变形受限，流场扰动较弱，推力最低。

2）摆频对流场结构影响分析

如图10所示，当摆频增大时，尾鳍两侧压差带的形成速度显著加快，压力及压力差值同步增大，推力随之增加。同时，全流场范围内的压力分布呈现结构性重构：高频运动促使尾鳍前缘高压区与后缘低压区的梯度差强化，伴随涡流脱落频率与流场扰动强度的协同增强，表明频率变化不仅影响压力参数，更通过改变流体与尾鳍的动态耦合模式，对涡系结构、流场能量耗散及推力生成机制产生系统性影响。

3.4 涡量图分析

在尾鳍摆动推进过程中，采用Q标准判据（Q=0.0002 s⁻¹）提取流场涡特征结构，图11为计算的流场涡结构结果。分析尾鳍周期摆动的涡旋结构发现：(N+1/4)T时刻，尾鳍处于上极限位置，后方涡量结构较弱，新涡量尚未形成，先前涡量逐渐脱落至尾迹，涡量生成速率低，尾鳍对流体作用力最小。(N+1/2)T时刻，尾鳍返回初始位置，后方生成强烈涡量对，涡量生成速率与尾鳍运动速度均达峰值，显著的剪切作用促使强涡对结构形成并开始脱落，流体扰动最大，对应推力峰值。(N+3/4)T时刻，尾鳍到达另一侧极限位置，涡对结构远离尾鳍进入尾迹区，新涡量生成减弱，涡量场稀疏，推力接近0或负值，系统处于反向换向阶段。(N+1)T时刻：尾鳍换向后加速运动，涡量生成逐渐增强，涡量场趋于清晰，尾鳍回到中间位置并重新推动水流形成新涡量对，推力开始上升。

如图12所示，在0.5 Hz摆动频率下，涡量分布稀疏，涡结构间距大且形态规整，呈现类周期性排列，涡街有序，涡量扩散缓慢且均匀，远离尾鳍区域衰减平缓。当频率增至0.8 Hz，涡量场显著变化，单位时间内涡量生成量增加，结构间距缩小，分布更紧凑；尾鳍边缘等关键区域的涡量强度与空间梯度急剧提升，表现出局部集中特征，反映流体扰动能量的显著聚集。频率进一步升至1 Hz时，涡量场展现高度动态特性，涡量生成效率大幅提高，空间分布极其密集，结构间距趋近于0，形成紧实涡团簇；涡形态非周期性且波动明显，涡量在尾鳍附近高度集中，扩散受抑，远离区域衰减迅速，难以形成稳定的扩散模式。

4 结　语

在这项研究中，本文采用双向流固耦合技术模拟了鲣鱼尾鳍的运动，对仿生鲣鱼尾鳍的水动力特性进行了数值分析。结果发现：

1）尾鳍摆角幅度在25°时推力相比15°提升约20%，当摆角幅度超过30°时推力增幅有限，侧向力波动加剧，稳定性下降。摆动频率在0.8～1 Hz范围内柔性尾鳍游动效率最优。

2）变厚度尾鳍全频率范围内推进性能均优于等厚度尾鳍。低频（0.5 Hz和0.8 Hz）时，结构钢材料尾鳍推力最大，高频（1.0 Hz和1.25 Hz）时，仿生材料推力最大。

3）结构变形与摆频显著影响流场与推力。仿生材料因优异柔性变形特性，增强压力梯度，推力最大，优于聚乙烯及结构钢；后者因变形受限推力最低。摆频提升加速压差带形成，增强涡流脱落与流场扰动。

4）涡量分析显示，摆动周期内，(N+1/2)T涡对形成推力达峰，(N+3/4)T涡场稀疏推力近0。频率从0.5 Hz增至1 Hz，涡量分布从稀疏规整变为密集动态，1 Hz时形成紧实涡团簇，反映频率对流场结构及推力输出的系统性影响。