| Issue |
JNWPU
Volume 43, Number 1, February 2025
|
|
|---|---|---|
| Page(s) | 76 - 83 | |
| DOI | https://doi.org/10.1051/jnwpu/20254310076 | |
| Published online | 18 April 2025 | |
Numerical simulation of manta ray's self-propulsion
蝠鲼自主游动数值模拟研究
1
Air Defence and Antimissile School, Air Force Engineering University, Xi'an 710038, China
2
School of Marine Science and Technology, Northwestern Polytechnical University, Xi'an 710072, China
Received:
10
January
2024
A manta ray is a typical representation of its pectoral fin propulsion mode, which is stable and conducive to wide area cruise, thus being suitable for bionic targets of a bionic vehicle. This paper established a computational model of the two-degree-of-freedom self-propulsion of the manta ray. Its motion equation and fluid dynamic equation were coupled to numerically simulate its self-propulsion process from stationary-state start-up and acceleration to steady-state cruise. The time history changes of the manta ray's swimming speed, hydrodynamic force, pressure distribution and three-dimensional flow field structure were analyzed. The simulation results show that, during the self-propulsion process of the manta ray, its speed, acceleration and displacement in its forward direction are completely determined by the net thrust generated during the flexible deformation of the pectoral fins and the net resistance encountered during its forward swimming. In the acceleration stage, the net thrust is superior and mainly generated near the tips of the pectoral fins. When the balance between the thrust and the resistance is reached, the manta ray is in its steady cruise stage. The flexible deformation of the spanwise and chordwise superposition of the pectoral fins may produce complex three-dimensional vortex structures. The numerical simulation method proposed in this paper and the study of the manta ray's self-propulsion process lay the foundation for further revealing its swimming mechanisms.
摘要
蝠鲼是胸鳍推进模式的典型代表, 稳定性好, 利于广域巡游, 非常适合作为仿生航行器的仿生目标。建立蝠鲼两自由度自主游动的计算模型, 耦合求解蝠鲼运动方程与流体动力方程, 数值模拟了蝠鲼从静止状态启动加速, 最终到达稳态巡游的自主游动过程, 分析了蝠鲼游动速度、水动力、压力分布以及三维流场结构的时间历程变化。计算结果表明: 在蝠鲼自主游动过程中, 前游方向的速度、加速度、位移等完全由胸鳍柔性变形时产生的净推力以及前游时受到的净阻力共同决定, 在加速阶段净推力占优, 且推力主要由靠近鳍尖部分产生, 推阻力平衡时蝠鲼到达稳态巡游阶段, 蝠鲼胸鳍展向及弦向叠加的柔性变形会产生复杂的三维涡结构。文中提出的数值仿真方法以及对蝠鲼自主游动过程的研究为后续揭示蝠鲼游动机理奠定了基础。
Key words: manta ray / self-propulsion / numerical simulation / wake structure
关键字 : 蝠鲼 / 自主游动 / 数值模拟 / 涡结构
© 2025 Journal of Northwestern Polytechnical University. All rights reserved.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
鱼类经过上亿年演化得到低能耗与高效率的推进方式, 这是与自然和谐共生的结果, 是传统推进器无法企及的。蝠鲼以中间鳍-对鳍(MPF)作为稳态巡游的方式[1], 机动性好, 效率高于鳗科、鲹科等水中生物[2], 具有很高的仿生价值。
Liu等[3]采用浸没边界法(IBM)计算了一组运动参数下蝠鲼的流体动力性能和尾涡结构, 发现净推力与胸鳍弯曲角度相关, 且推力主要是由胸鳍远端部分产生。杨少波[4]根据牛鼻鲼外形相似、质量及刚度等效的原则建立三维有限元模型, 采用鳍条模仿肌肉的内力作用, 利用多鳍条组合实现胸鳍拍动, 结果表明鳍面上自前向后的波动对牛鼻鲼有推进作用。张栋[5]提出一种改进的IB-SGKS方法研究了牛鼻鲼主动推进的水动力特性。Menzer等[6]研究了蝠鲼的水动力性能和前缘涡的变化, 发现胸鳍俯仰角对推力和效率都有显著影响。Sharma等[7]用圆盘模型模拟鳐鱼运动, 改变展弦比和波长类比鳐鱼和牛鼻鲼的运动, 研究尾涡的变化。
以上研究均是非自主游动, 采用均匀来流进行仿真, 不符合真实生物的运动, 应该让蝠鲼依靠其所受到的力找到自己的“路”。Wu[8]强调非自主游动与自主游动的受力与流场是不同的, Lin等[9]的研究结果也明确体现了这个结论。
鳐科的胸鳍推进模式一直是鱼类研究中的薄弱点, 起步较晚, 蝠鲼在鳐科中属于拍动模式, 或mobuliform模式, 相关的自主游动研究非常少, 且多考虑前进方向的性能。伍志军[10]对晶吻鳐模型的直线前游进行数值仿真研究。Fish等[11]计算了蝠鲼在不同运动频率下的前游速度与自主游动效率, 随着频率增大, 前游速度增大, 稳态速度无量纲化后均为1.5倍体长。Bianchi等[12]采用嵌套网格, 通过OpenFoam计算牛鼻鲼由静止加速至稳态速度的单自由度自主游动过程。
本文建立与真实蝠鲼运动较为相似的蝠鲼胸鳍变形方程, 提出自主游动的数值计算方法, 对流体动力学方程与蝠鲼动力学方程的耦合过程进行求解, 对蝠鲼2-DOF自主游动过程中的速度特性、水动力特性、表面压力分布和流场涡结构展开研究, 以揭示蝠鲼自主游动的机理。
1 计算模型与数值计算方法
本文使用的蝠鲼物理模型根据真实蝠鲼外形[13]绘制而成, 对吻鳍、眼睛、尾部进行了简化, 最终的物理模型如图 1所示, 展长和体长分别为LS和LB, 按照成年蝠鲼比例约束物理模型, LS/LB=0.78, 与真实生物外形具有较高的的一致性。
 |
图1 蝠鲼数值计算模型 |
全局坐标系为O-XYZ, 随体坐标系为o-xyz, 纵轴ox沿蝠鲼中心线指向尾部, 横轴oy沿右侧胸鳍展向方向, 垂直轴oz指向胸鳍上挑方向, 本文主要讨论xoz平面内前游与上浮下潜2个自由度的自主游动, 因此设置蝠鲼随体坐标系与全局坐标系的Y轴方向相同。
在随体坐标系下, 初始时刻胸鳍上任意一点的坐标为(xf, yf, zf), 胸鳍变形至任意时刻下坐标为(x(xf, yf, t), y(xf, yf, t), z(xf, yf, t)), 胸鳍对称运动时鳍尖最大幅值无量纲化表示为A, 转动过程中胸鳍上不同展向位置对应的最大转动角θmax不同, 从鳍尖到胸鳍根部递减, 另外采用系数k控制身体的拉伸量保持在一个较小的范围, 蝠鲼胸鳍运动控制方程为
式中: ω=2πf, f是胸鳍拍动的频率, 波数W是在胸鳍弦向上完整行波的个数, 表示为W=LB/λ, λ是行波的波长。根据以上运动方程, 当蝠鲼做上下幅值对称运动时, 一个周期内蝠鲼胸鳍运动如图 2所示, 可以看出与真实蝠鲼的游动姿态有较高的相似性。
 |
图2 蝠鲼在一个周期内的运动形态 |
数值计算中, 本文采用三维不可压缩黏性Navier-Stokes方程作为流体运动的控制方程, 可以表示为
式中: u是流体速度; p是流体压力; ν是流体运动学黏性系数。
本文基于有限体积法(FVM)求解控制方程, 采用SIMPLEC算法对压力和速度进行耦合, 压力项为二阶格式, 动量项为二阶迎风格式, 时间项采用一阶隐式格式。
本文对蝠鲼自主游动的研究只包括整体平动运动, 因此其速度可分解为整体平动 uc=(ux, uz)以及自身柔性变形运动u′, 表示为
蝠鲼的平动由外力决定, 根据动量定理得到蝠鲼动力学方程
式中：m是蝠鲼质量；F为求解流体动力方程得到的水动力。将蝠鲼表面的压力和切应力积分得到x, z方向的瞬时作用力
式中：dA为蝠鲼表面微元; p为dA上的压力; ni为dA上法向向量的第i个分量; τ1i, τ3i为黏性应力张量。根据Fx(t)在周期变化中正负的变化, 将其分为净推力FT(t)和净阻力FD(t), 分别表示为
根据本文坐标系设置, 蝠鲼前游方向为全局坐标系X轴的负方向, 因此定义蝠鲼前游速度uf=-ux, 前游方向的纵向力Ff=-Fx。蝠鲼由静止加速前游, 最终到达稳定巡游状态, 此时uf的时均值恒定, 即为巡游速度Us。由于Us不是一个独立值, 所以选择特征速度U0=1LB/s将水动力参数无量纲化, 选择特征加速度a0=1LB/s2将游动加速度无量纲化
由图 2可见蝠鲼两侧胸鳍柔性变形较大, 因此采用混合网格, 网格细节如图 3所示, 将靠近蝠鲼的内域划分为四面体非结构网格, 动网格变形时使用弹簧光顺法和局部网格重构法, 其余部分划分为六面体结构网格, 不参与网格变形, 捕捉尾涡更为精确, 这样的网格设置使计算更稳定并提高计算效率。
 |
图3 内部网格细节 |
2 计算方法验证
本文为了验证网格的无关性, 选取3种网格尺寸, 网格数量分别为2.58×106, 3.82×106, 5.22×106。图 4为f=1 Hz, A=0.3LB, W=0.4时前游速度的变化曲线, 可以看到网格2与网格3的结果无明显差异, 综合计算精度与计算效率采用网格2进行数值仿真。
 |
图4 网格无关性验证 |
本文选择3个时间步长分别为T/200, T/500, T/800, 且每个计算步中迭代50步, 前游速度曲线变化如图 5所示, T/500与T/800步长下的速度曲线较为接近, 满足时间步长无关性, 故选择Δt=T/500作为本文计算的时间步长。
 |
图5 时间步长无关性验证 |
为验证数值计算方法的准确性, 与文献[14]中仿生鱼的两自由度自主游动结果进行对比, 本文计算的前游速度U和侧向速度V如图 6所示, 最终巡游速度为0.43LS, 侧向速度围绕零值上下波动, 与文献[14]结果0.4LS基本吻合, 故本文的数值计算方法正确有效。
 |
图6 数值计算方法验证 |
3 数值计算结果与分析
本文主要讨论蝠鲼自主游动过程中运动特性、水动力特性、表面压力以及流场结构随时间的变化, 选择参数组合幅值A=0.3、频率f=1、波数W=0.4进行分析, 该参数组合确定蝠鲼胸鳍变形。
3.1 自主游动过程中的运动特性
图 7给出了蝠鲼自主游动全过程中加速度、速度、位移随时间的变化曲线, 本文关注点为蝠鲼前游运动, 故只需研究在x-z平面内的运动。
 |
图7 蝠鲼自主游动过程 |
图 7a)展示了蝠鲼前进方向的运动特性。在t=0时刻, 蝠鲼与流体均为静止状态, 随后蝠鲼开始运动并作用于流体, 在相互作用下流体作用在蝠鲼上的力并不平衡, 使蝠鲼获得一定的平动加速度及速度。当0 < t/T≤9时, 蝠鲼处在加速阶段, 加速度ax振荡变化, 但平均加速度仍大于0, 与之对应的前游速度uf由零逐渐阶梯式增大, 当t/T>9时, 流场逐渐趋于稳定, 作用在蝠鲼上的力也达到平衡, 加速度围绕零值稳定波动, 使uf以一定的幅值稳定波动, 即进入稳态巡游阶段。其稳定后的周期平均前游速度即为巡游速度Us=0.998LB/s, 根据雷诺数的定义Re=UsLB/ν, 蝠鲼在该参数组合游动时的雷诺数为3.4×106。与此同时, 蝠鲼启动后保持前游状态, 前游位移持续增大。
图 7b)为蝠鲼自主游动过程中z方向的运动变化。z向加速度az和速度uz在启动后迅速转为稳定波动变化, 周期平均速度和平均加速度为零。从位移曲线(t-sz)可以看出, 蝠鲼胸鳍上挑时, 身体稍微下沉, 向下拍动时身体稍微上浮, 并一直保持小幅度上下波动, 这与蝠鲼实际游动完全一致, 并且随着时间的推进, 整体位移略微增大, 蝠鲼逐渐上浮。这与蝠鲼上下不对称的外形有关, 下表面较上表面更平坦。其次与蝠鲼和水之间的质量比有关, 若质量比等于1, 即蝠鲼密度与水密度相同, 游动时会轻微上浮; 若质量比小于1, 即蝠鲼密度小于水密度, 游动时上浮位移更大; 若质量比大于1, 即蝠鲼密度大于水密度, 游动时上浮位移更小, 但一定存在一个临界质量比使蝠鲼保持在固定深度前游, 当大于该临界值时会逐渐下沉。在实际应用方面, 自主游动计算得到垂向位移的变化，对于仿生航行器的深度控制有一定帮助, 调节浮力即为调节质量比, 能够帮助判定上浮下沉的深度。
3.2 自主游动过程中的水动力特性
图 8展示了蝠鲼自主游动过程中水动力系数随时间变化的曲线。蝠鲼前进方向的速度、加速度、位移完全由胸鳍拍动过程中产生的推力及蝠鲼前进过程中受到的流体阻力共同决定。在蝠鲼初始运动阶段, 蝠鲼所受摩擦力等阻力CD较小, x方向合力系数Cf大于零则推力占优, 游速迅速增大, 蝠鲼加速向前游动, 随着游动速度uf的增大, 摩擦阻力随之逐渐增大, Cf减小到负值时, 阻力大于推力使游速减小, 当游速减小到最小值时, 推力大过阻力, 蝠鲼又开始加速。如此循环往复, 后期游动收敛到稳态后, Cf均值围绕零值上下波动, 推力CT与阻力CD达到动态平衡, 虽然波动幅值不同, 但平均值相等, 从而使Cf呈周期性变化, 且合力均值为零, 使蝠鲼在稳定游动时的周期平均加速度ax为零, 蝠鲼以稳定波动的速度巡游。
 |
图8 自主游动过程水动力系数的时间历程 |
对于非自主游动情况, 在均匀来流下, 计算得到的x向作用力均值通常不为零, 这与稳态游动时需要满足的受力平衡条件不符, 因此非自主游动计算得到的受力是不准确的。
对于z向作用力, 即升力系数CL, 启动后, 随着蝠鲼在z向的上下波动, CL逐渐增大并很快稳定下来, 且均值为零, 其巡游阶段的波动幅值约为Cf的10倍。蝠鲼左右两侧外形相同, 且运动相同, 两侧胸鳍拍动产生的沿y方向上的合力始终为零, 故y向作用力系数为零。
在蝠鲼游动收敛到稳态以后, 由于自身周期性运动而使流体作用力也周期性交替变化。图 9给出了蝠鲼稳态游动时，x向、z向作用力系数随时间的变化曲线。从图上可以看出, Cf波动周期为蝠鲼运动周期的2倍, CL波动周期与蝠鲼运动周期相同。当胸鳍从幅值最小处上挑时, 升力单调递增, 在幅值最大处到达峰值, 下拍时升力减小, 最大推力出现在上挑下拍的平衡位置之前。
 |
图9 稳态巡游时受力的时间历程 |
3.3 表面压力分布
为了揭示蝠鲼自主游动的推进机理, 对蝠鲼在启动加速和巡游阶段的身体表面压力分布云图进行分析。
图 10为启动加速阶段的蝠鲼表面压力分布云图。胸鳍从最低点上挑时, 胸鳍远端上表面为明显的高压区, 下表面则为明显的低压区。由于胸鳍运动产生柔性弯曲, 使得高压区与低压区之间形成压力差并产生一个与前游方向相同的分力以及向下的分力, 分别表现为向前的正推力和向下的升力。随后胸鳍运动至最高点时, 虽然上表面为高压区, 下表面为低压区, 但此时胸鳍变形呈迎水斜面, 压差产生一个与前游方向相同的分力以及向上的分力, 分别表现为向前的正推力和向上的升力。所以蝠鲼在启动阶段鳍尖部分是主要产生推力的部位。
 |
图10 蝠鲼在启动加速阶段的压力云图 |
图 11为巡游阶段的蝠鲼压力分布云图。相较于加速阶段的云图, 巡游时蝠鲼身体表面的低压区面积更大, 头部和胸鳍前缘依旧是阻力来源的高压区, 形成逆压梯度, 胸鳍后缘与尾部的高压区形成的正压梯度有效克服了头部的逆压梯度, 使蝠鲼能够保持稳态运动。
 |
图11 蝠鲼在巡游阶段的压力云图 |
3.4 流场结构演化
为了进一步理解蝠鲼自主游动机理, 需要对前游产生的三维流场结构进行分析。本文采用Hunt等[15]提出的Q准则可以计算流场中每个位置的Q值，并使用等值线图绘制Q值的分布，可以直观地展示涡旋的位置、形状和强度。本文采用Q=0.12的等值面图。
图 12为蝠鲼启动加速阶段的三维流场涡结构。在1.0T时胸鳍在幅值最低点, 产生前缘涡LEV和鳍尖涡T1, 同时在胸鳍后缘及身体后端生成后缘涡TEV。随着胸鳍上挑鳍尖涡向下脱落, 同时生成第一个涡环V1, 在缓慢前游时身体后端形成的涡逐渐脱落, 但不构成一个环路或涡环, 称之为C形涡Vc, 这2个涡流向相反。胸鳍上挑至最高点时产生第2个鳍尖涡T2, 下拍时T2向上脱落, 身体后端形成3个Vc, 下拍到最低点时形成第2个涡环V2。在启动阶段蝠鲼前游速度较小, 尾流场中的涡来不及相互分离挤压在一起, 该阶段主要以前缘涡、后缘涡、鳍尖涡和尾端C形涡的生成与脱落为主。
 |
图12 蝠鲼在启动加速阶段的流场涡结构 |
图 13为蝠鲼巡游阶段的三维流场涡结构。从图 13中可以看出与启动加速阶段的涡结构有很大不同, 在这种情况下, 蝠鲼处于巡游阶段时, 尾流场中的涡结构会转变为单排的涡环结构, 每个涡环沿着流体运动方向被拉长, 使涡环变得更细, 这是由于蝠鲼在此时具有较大的前游速度, 在胸鳍上产生的涡环通过其自身产生的射流作用与蝠鲼持续前进相互作用, 迅速分离开来, 相邻的2个涡环不再相互挤压和交错。从侧视图观察, 所有的涡环形成了一个“正弦式”曲线, 这种排列方式减弱了相邻涡环之间的干扰, 更有利于周围流体向后运动, 因此蝠鲼能够以较大的速度向前游动。另外从俯视图能够看到最早生成的涡环在蝠鲼快速游动下消散较快, 后缘涡脱落更快, 最晚生成的涡环更靠近于鳍尖。
 |
图13 蝠鲼在巡游阶段的流场涡结构 |
4 结论
本文建立自主游动的数值计算方法, 对流体动力学方程与蝠鲼动力学方程的耦合过程进行求解, 对蝠鲼2-DOF自主游动过程中的速度特性、水动力特性、表面压力分布和流场涡结构展开研究, 得到的结论如下:
1) 采用自主游动数值计算方法可以有效地预测蝠鲼的瞬时游动速度、巡游速度、位移以及受力情况;
2) 蝠鲼自主游动过程分为加速阶段和巡游阶段, 从静止开始前游速度不断增加, 巡游时实现稳定波动, 并保持恒定的巡游速度, z向速度和加速度在启动后快速进入周期平均值均为零的稳定波动, 并且略微上浮;
3) 分离x方向受力中的推力和阻力是自主游动计算中的关键。稳态巡游时推力与阻力达到动态平衡, x方向合力为零, 完全不同于与非自主游动的情况;
4) 蝠鲼表面压力分布表明胸鳍上下表面交替出现高压区和低压区是蝠鲼向前游动的根本原因。胸鳍拍动时交替产生多个涡环, 相互诱导形成“正弦式”射流推动蝠鲼向前游动。
References
- SCHAEFER J TSUMMERS A P. Batoid wing skeletal structure: novel morphologies, mechanical implications, and phylogenetic patterns[J]. Journal of Morphology, 2005, 264(3): 298–313 [Article] [Google Scholar]
- BORAZJANI ISOTIROPOULOS F. Numerical investigation of the hydrodynamics of carangiform swimming in the transitional and inertial flow regimes[J]. Journal of Experimental Biology, 2008, 211: 1541–1558 [Article] [Google Scholar]
- LIU GREN YZHU J, et al. Thrust producing mechanisms in ray-inspired underwater vehicle propulsion[J]. Theoretical and Applied Mechanics Letters, 2015, 5(1): 54–57 [Article] [Google Scholar]
- YANG Shaobo. Hydrodynamic analysis of cownose ray's swimming and research of bionic robotic fish[D]. Changsha: National University of Defense Technology, 2010 (in Chinese) [Google Scholar]
- ZHANG Dong. Flexible deformation effect on the hydrodynamic performance of a rhinoptera javanica in different swimming behaviors[D]. Xi'an: Northwestern Polytechnical University, 2020 (in Chinese) [Google Scholar]
- MENZER AGONG YFISH F E, et al. Bio-inspired propulsion: towards understanding the role of pectoral fin kinematics in manta-like swimming[J]. Biomimetics, 2022, 7(2): 45–59 [Article] [Google Scholar]
- THEKKETHIL NSHARMA AAGRAWAL A. Three-dimensional biological hydrodynamics study on various types of batoid fishlike locomotion[J]. Physical Review Fluids, 2020, 5(2): 23101 [Article] [Google Scholar]
- WU T Y. On theoretical modeling of aquatic and aerial animal locomotion[J]. Advances in Applied Mechanics, 2001, 38: 291–353 [Google Scholar]
- LIN XWU JZHANG T. Self-directed propulsion of an unconstrained flapping swimmer at low Reynolds number: hydrodynamic behaviour and scaling laws[J]. Journal of Fluid Mechanics, 2021, 907: 1–14 [Google Scholar]
- WU Zhijun. Numerical simulation and experimental research on undulatory propulsion of raja eglanteria[D]. Harbin: Harbin Institute of Technology, 2015 (in Chinese) [Google Scholar]
- FISH F EDONG HZHU J J, et al. Kinematics and hydrodynamics of mobuliform swimming: oscillatory winged propulsion by large pelagic batoids[J]. Marine Technology Society Journal, 2017, 51(5): 35–47 [Google Scholar]
- BIANCHI GCINQUEMANI SSCHITO P, et al. A numerical model for the analysis of the locomotion of a cownose ray[J]. Journal of Fluids Engineering, 2021, 144(3): 031203 [Google Scholar]
- LI TLI GLIANG Y, et al. Fast-moving soft electronic fish[J]. Science Advances, 2017, 3(4): e1602045 [Google Scholar]
- KERN SKOUMOUTSAKOS P. Simulations of optimized anguilliform swimming[J]. Journal of Experimental Biology, 2006, 209(24): 4841 [Google Scholar]
- HUNT J C R, WRAY A A, MOIN P. Eddies, stream, and convergence zones in turbulent flows[C]//Studying Turbulence Using Numerical Simulation Databases, Stanford, CA, USA, 1988: 193–208 [Google Scholar]
All Figures
|
图1 蝠鲼数值计算模型 |
| In the text | |
|
图2 蝠鲼在一个周期内的运动形态 |
| In the text | |
|
图3 内部网格细节 |
| In the text | |
|
图4 网格无关性验证 |
| In the text | |
|
图5 时间步长无关性验证 |
| In the text | |
|
图6 数值计算方法验证 |
| In the text | |
|
图7 蝠鲼自主游动过程 |
| In the text | |
|
图8 自主游动过程水动力系数的时间历程 |
| In the text | |
|
图9 稳态巡游时受力的时间历程 |
| In the text | |
|
图10 蝠鲼在启动加速阶段的压力云图 |
| In the text | |
|
图11 蝠鲼在巡游阶段的压力云图 |
| In the text | |
|
图12 蝠鲼在启动加速阶段的流场涡结构 |
| In the text | |
|
图13 蝠鲼在巡游阶段的流场涡结构 |
| In the text | |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.