为了研究鱼类优异游泳性能背后的物理机制，国内外学者进行了大量活体观测实验，通过高速相机记录下鱼类在运动过程中鳍和身体的变形情况，最终总结出大量的运动方程和仿生外形模型[1-6]．通过活体实验来精确测量鱼类瞬时运动及变形难度较大，往往不能达到预期效果．目前最常见的研究方法是使用降阶物理模型，如俯仰或升沉扑翼代替鱼鳍运动和外形，并在扑翼前缘植入一个周期性的输入信号来控制扑翼运动[7-10]．文献[7]采用浸入边界法研究了雷诺数(Re)和翼型厚度对自推进翼型的水动力特性及尾流结构的影响．文献[8]研究了非正弦输入信号对俯仰翼尾流结构的影响，使用可调参数K来控制输入信号波形，随着K的增大，易使摆动翼产生的尾迹结构失稳．然而，简单的控制信号并不能准确代表生物运动机制，鱼类在推进过程中鳍和身体的运动会产生复杂波形，在鳍和身体上出现不同频率和振幅的多波．从生物控制系统来看，输出控制信号由基本输出信号和由各种感官输入产生的调制信号共同决定．从外部影响因素来看，无论是鱼还是水下航行器所处的海洋环境都瞬息万变，水下复杂的流场环境会在一定程度上影响运动控制．本研究采用数值模拟的方法，将一个高频同相正弦扰动引入到升沉扑翼的基础正弦输入信号中，对升沉扑翼在不同基频和扰频组合下的推力系数、推进效率及尾涡结构展开研究，以期为水下航行器设计提供参考．1 数值计算方法1.1　问题描述数值计算模型采用NACA 0012翼型，如图 1所示，刚性翼型放置在速度为U的均匀流场中，x方向指向来流方向，y方向垂直于来流方向，图中：A为振幅；c为弦长；h(t)为运动方程．10.13245/j.hust.220421.F001图 1数值计算模型输入信号由基础信号和扰动信号两部分构成，其中基础信号的升沉振幅A为0.25c(弦长)，频率范围为0.50~1.75 Hz；扰动信号的升沉振幅为0.1A，频率范围为0.5~12.0 Hz．输入信号的波形图见图2所示，定义式为h(t)=Asin(2πfbt)+0.1Asin(2πfpt+φ)，(1)式中：t为时间；fb和fp分别为基频频率和扰动频率；φ为基础信号与扰动信号的相位差，设为0．图 2为输入信号波形图， 图中：h为无量纲升沉幅值；T为周期时长．10.13245/j.hust.220421.F002图 2输入信号波形图本研究引入无量纲参数来描述扑翼运动，即h，Re，无量纲频率k．为了衡量升沉扑翼推进性能，引入符号CT和CL表示瞬时推力系数、瞬时升力系数．定义式分别为：h=A/c；Re=U/υ；k=2πfc/U；CTM=1T∫tt+TCT(t)dt;CLM=1T∫tt+TCL(t)dt,式中CTM和CLM分别为平均推力系数、平均升力系数．从文献[11]得知，瞬时输入功率系数CP(t)=CL(t)h˙(t)．(2)因此，时均输入功率系数(CPM)计算公式为CPM=1T∫tt+TCL(t)h˙(t)Udt．升沉扑翼的推进效率定义为η=CTM/CPM．(3)1.2　计算方法采用 Fluent软件进行计算，采用有限体积法离散纳维-斯托克斯方程，采用二阶迎风法进行动量离散，并基于格林-高斯节点计算导数，采用 PISO (压力的隐式算子分割)算法对连续方程中的压力和速度进行耦合，湍流模型采用SST k-ω模型，通过用户定义函数(UDF)控制输入信号．计算域由内域和外域两部分组成，计算域总尺寸为20c×12c的长方形，内域尺寸为10c×3c．内域采用非结构化网格通过动态扩散网格法减少计算资源，避免网格明显变化，外域采用结构化网格．计算区域边界条件设置为：入口边界为速度入口；出口边界为压力出口；扑翼表面及上下边界均采用无滑移边界条件．2 收敛性与数值验证2.1　收敛性验证为验证网格的收敛性，选取大、中、小三种网格尺寸，网格1~3包含网格数量分别为9×104，2×105和3×105．以升沉扑翼在h=0.25，fb=0.5 Hz，fp=0 Hz情况下的CT曲线进行网格收敛性验证，网格1~3的收敛性验证如图 3(a)所示．由图 3(a)可知：三种网格尺度下的CT变化无明显差异，表明数值计算是网格收敛的．为保证计算精度及节约计算资源，后续结果分析采用网格2的计算结果．在此基础上还进行了时间步验证．以升沉扑翼在h=0.25，fb=0.5 Hz，fp=0 Hz情况下的CTM曲线进行时间步长(ΔT)验证如图3(b)所示，将时间步从每周期1 000步增加到1 500步不会改变推力系数．因此，在后续研究中时间步采用1 000．10.13245/j.hust.220421.F003图 3网格收敛性验证2.2　数值计算方法验证为了验证计算结果，将结果与文献[12-13]的结果进行了比较．结果表明：当kh＜1.5时，层流模拟结果与湍流模拟结果相符合且与实验数据相近．从图 4可以看出：本文方法的计算结果与实验结果符合较好，且优于文献[12]计算结果．10.13245/j.hust.220421.F004图4不同kh下的平均推力系数3 数值计算结果与分析3.1　推进性能fb和fp存在较大差异，不能仅根据其中任何一个频率来确定计算周期，因此本研究采用选取基础信号和扰动信号运动周期的最小公倍数作为本研究的计算周期．例如，当基频为0.5 Hz，扰频为10 Hz时，将计算周期设置为2 s．总计算时间选为计算周期的10倍，并在最后两个计算周期取平均作为最终结果．图5展示了升沉扑翼在不同fb下平均推力系数CTM随扰动频率fp的变化情况．当扰频为0 Hz时，表示输入信号仅由基础信号决定，作为对照组，便于探究扰动信号对扑翼推力系数的影响．由图5可以看出：当fb＜1.25 Hz时，扰动信号对推力系数有显著影响，尤其在最低基频0.5 Hz处，当引入一个12 Hz的扰频时，它的推力系数变化为对照组(fp=0 Hz)的8.435倍．然而，当fb＞1.25 Hz时，尽管CT总体上呈上升趋势，但在某些扰频上存在一些显著波动．例如，当基频为1.5 Hz、扰频为0.5 Hz以及基频为1.75 Hz、扰频为4 Hz时，CTM均出现明显下降．10.13245/j.hust.220421.F005图5不同基频下的平均推力系数和扰动频率升沉扑翼在不同fb下η随fp的变化情况如图6所示．从图中可见η随着fp的增加总体呈现下降趋势，最值得注意的是：在较低基频(fb=0.5，1.0 Hz)处，η随着fp的增加显著下降．从而在图上呈现出：当扰频fp＜4 Hz时，低基频扑翼推进效率远远大于高基频扑翼推进效率；当fp≥4 Hz时，高基频扑翼推进效率开始反超低基频扑翼推进效率．本研究还发现：当fb=1.5 Hz，fp=1 Hz时，扑翼推进效率出现骤降．10.13245/j.hust.220421.F006图6不同基频下的推进效率和扰动频率结合图5和6发现：在推力系数巨幅提升的同时，推进效率却在随之减小．解释这一现象须要从效率的定义(见式(3))出发，效率反映了推力系数与功率系数的相对关系，它并不是由推力系数单一决定．为了探索效率下降的物理机理，图7给出CPM在不同fb下随fp的变化情况．从图7可见：随着扰频逐渐增大，输入功率系数也逐渐增大．由式(2)可知功率系数与升力系数及运动速度有关，随着扰动频率的增加，运动速度必然急速增加，从而导致输入功率系数增加．观察图5和6的纵坐标数量级可以看出输入功率系数增加幅度远远大于推力系数的增大幅度，这必然导致扑翼推进效率随着扰频的增加而逐渐减小．10.13245/j.hust.220421.F007图7不同基频下的平均输入功率系数和扰动频率由式(1)和(2)可知：CP与h.(t)紧密相关，且h.(t)由基础信号和扰动信号决定．当fp＜4 Hz时，扰频对输入功率系数影响甚微，主要由基频决定输入功率系数大小，随着基频增加，输入功率系数显著增加；当fp≥4 Hz时，h.(t)中两部分分量(基础信号和扰动信号)趋向同阶次，因此各基频下扑翼输入功率差异减小．由图5看出：当fp＜4 Hz时，推力系数相差不大；当fp≥4 Hz时，基频较高扑翼推力系数远超基频较低扑翼的推力系数．这也很好地解释了为何扑翼推进效率在fp=4 Hz处出现分界．扑翼推进效率在fb=1.5 Hz，fp=1 Hz处出现骤降的原因是扑翼在此运动参数下推力系数出现下降，然而输入功率系数相较于相邻两种工况并没出现明显变化，从而导致了其推进效率突然大幅下降．3.2　流场结构为了深入探究扰动信号对扑翼推进性能的影响，本研究对5组基频扑翼在不同扰频下的流场进行了分析．分析结果显示5组扑翼流场变化大致分为三类：当fb=0.50 Hz时，随着扰频增加，扑翼尾涡结构由无明显涡脱类型转变为2S模式(即每个周期有两个单独的漩涡生成)；当fb=0.75 Hz时，随着扰频增加尾涡结构由2P模式(即每个周期有两对漩涡生成且每对漩涡由两个方向相反的漩涡组成)演变为2S模式；在其他4种基频下，扑翼尾涡的结构均由反卡门涡街演变为杂乱模式，并伴随有前缘涡的脱落．因为流场结构相似，所以本研究仅列举当基频为0.50，0.75，1.75 Hz时，扑翼处于不同扰频在一个基础信号周期内5个不同时刻(0，T/4，T/2，3T/4，T)的尾涡结构，初始运动方向均沿y轴的正向．图8展示了升沉扑翼当fb=0.5 Hz时，尾涡结构随扰频的变化情况．当fp＜4 Hz时，尾流场中并未发现明显的旋涡，这很好地解释了为何当fp=0.0，0.5，1.0，2.0 Hz时，扑翼推力系数值较低且没有太大变化．当fp≥4 Hz时，从图8中很清晰地看见越来越多的反向旋涡开始从扑翼尾端脱落，形成反卡门涡街，且旋涡密集程度随着扰频增大而增大．同时观察到随着扰频增大，旋涡在下游保持更远的距离后才消散，这大大提升了扑翼推进性能．由图8还可以看出：当扑翼由中心位置向上运动时，涡对向y负向偏移；当扑翼到达上极位往中心位置运动时，涡对向y正向偏移，涡街偏转角度随扑翼运动改变而发生变化．10.13245/j.hust.220421.F008图8fb=0.50 Hz处的扑翼尾涡结构图图9展示了当fb=0.75 Hz时，升沉扑翼尾涡结构随扰频的变化情况．和图8对比，可以看出：当没有扰动信号及扰频较低时，扑翼涡量场存在明显差异．结合图8和9中fp=0 Hz处的涡量场来描述差异，从图9中可以看出：在0和T/4时刻，第一、二个负涡开始从扑翼尾端脱落，在T/2和3T/4时刻发生了第一、二个正涡的脱落，即是说在每个半周期中会出现一对旋转方向相同的涡，且尾迹中两个涡与中心线完全对称，无偏转，与2P模式相似．随着fp的增大，尾涡结构同fb=0.5 Hz时类似，呈现为2S模式，且扰频越大涡对越密集，涡街强度也越高，从而提高了推力系数．同时可以看出：涡对偏移方式和图8相同，当扑翼由中心位置向上运动时，涡对向y负向偏移；当扑翼到达上极位往中心位置运动时，涡对向y正向偏移．10.13245/j.hust.220421.F009图 9fb=0.75 Hz处的扑翼尾涡结构图图10展示了当fb=1.75 Hz时，升沉扑翼尾涡结构随扰频的变化情况．当扰频较低时，不仅在扑翼尾端产生了明显的反卡门涡街，并且在扑翼前缘上下两侧产生了前缘涡，同时可以观察到没有偏转涡街的出现．随着扰频的增大，前缘涡向尾流扩散打乱了反卡门涡街，使得尾涡结构由2S模式转换为混乱模式．在基频较高的情况下，扰频增加促使扑翼推力系数增大的来源有两个：一是扰频增大使得涡对越密集，涡街强度也越高；另一个重要因素是前缘涡的发展过程，前缘涡环绕前缘从一侧向另一侧运动，在扑翼前端产生一个低压区，进而产生一个吸力区，使得扑翼推力进一步提高．10.13245/j.hust.220421.F010图 10fb=1.75 Hz处的扑翼尾涡结构图纵观图8~10，基频的不同造成了扑翼初始尾涡结构的不同，当fb=0.50 Hz时，初始尾涡结构并没表现出明显模式；当fb=0.75 Hz时，初始尾涡结构展现为2P涡；当fb＞0.75 Hz时，初始尾涡结构转变为了反卡门涡街．这解释了为何当fp=0 Hz时，随着基频的增加，推力系数逐渐增大．在引入扰动信号后，当基频较低时(fb=0.50，0.75 Hz)，尾涡结构由初始形态(无明显模式或2P模式)转换为2S模式，随着扰频的逐渐加大，涡对变得更加密集，涡能量也逐渐增加，这给扑翼的推进性能带来了大幅增益；当基频较大时(fb=1.00，1.25，1.50，1.75 Hz)，涡街能在下游保持更远的距离才开始消散，前缘涡向下游扩散使得尾涡结构由2S涡演变为杂乱模式，也正因为前缘涡的发展进一步提升了扑翼推进性能．4 结论a. 随着扰动频率的增加，升沉扑翼的平均推力系数大幅提升，尤其当fb=0.50 Hz时，引入12 Hz扰频使得平均推力系数提升7.43倍．b. 随着扰动频率的增加，扑翼在相同时间内位移量大幅提升，导致输入功率系数增加，从而引起升沉扑翼推进效率降低．c. 基础频率决定了扑翼尾涡结构的初始形态，当fb=0.50 Hz时，尾涡结构初始形态为非明显模式；当fb=0.75 Hz时，尾涡结构初始形态为2P模式；当fb＞0.75 Hz时，尾涡结构初始形态为2S模式．d. 引入扰动频率，促使fb=0.50，0.75 Hz的低基频扑翼尾涡结构转换为2S模式，大幅提升了扑翼推进性能．当扑翼处于fb=1.00，1.25，1.50，1.75 Hz等高基频时，扰动信号的引入导致了前缘涡产生和脱落，导致扑翼尾涡结构演变为杂乱模式，并在扑翼前端形成吸力区，进一步提升扑翼推进性能．